Mercurial > hg > octave-nkf
annotate scripts/linear-algebra/krylov.m @ 20579:131ce8cfaa80
Relax input in functions that convert between colorspaces (bug #45456)
* scripts/image/hsv2rgb.m, scripts/image/ntsc2rgb.m, scripts/image/rgb2hsv.m,
scripts/image/rgb2ntsc.m: remove all input check and leave it up to new
private functions handled by all. This adds support for ND images, drops
check for values in the [0 1] range, allows integer images, and allows
sparse matrices. Also adjust tests and add extra ones.
* scripts/image/private/colorspace_conversion_input_check.m,
scripts/image/private/colorspace_conversion_revert.m: all of this functions
handle input check in the same way. In the same way, they need to prepare
the output in the same way based on what preparation was done during input
check (transforming image into colormap). So we create two new private
functions to avoid repeated code. All code was adapted from hsv2rgb.
| author | Carnë Draug <carandraug@octave.org> |
|---|---|
| date | Sun, 19 Jul 2015 17:41:21 +0100 |
| parents | 83792dd9bcc1 |
| children |
| rev | line source |
|---|---|
|
19898
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19232
diff
changeset
|
1 ## Copyright (C) 1993-2015 Auburn University. All rights reserved. |
| 3427 | 2 ## |
| 3 ## This file is part of Octave. | |
| 4 ## | |
| 5 ## Octave is free software; you can redistribute it and/or modify it | |
| 7016 | 6 ## under the terms of the GNU General Public License as published by |
| 7 ## the Free Software Foundation; either version 3 of the License, or (at | |
| 8 ## your option) any later version. | |
| 3427 | 9 ## |
| 7016 | 10 ## Octave is distributed in the hope that it will be useful, but |
| 11 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
| 13 ## General Public License for more details. | |
| 3427 | 14 ## |
| 15 ## You should have received a copy of the GNU General Public License | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 3427 | 18 |
| 3439 | 19 ## -*- texinfo -*- |
|
11470
eb9e0b597d61
Use common names for variables in documentation and code for a few more m-script files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
20 ## @deftypefn {Function File} {[@var{u}, @var{h}, @var{nu}] =} krylov (@var{A}, @var{V}, @var{k}, @var{eps1}, @var{pflg}) |
| 5555 | 21 ## Construct an orthogonal basis @var{u} of block Krylov subspace |
| 22 ## | |
| 23 ## @example | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
24 ## [v a*v a^2*v @dots{} a^(k+1)*v] |
| 5555 | 25 ## @end example |
| 26 ## | |
| 27 ## @noindent | |
|
20370
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
28 ## using Householder reflections to guard against loss of orthogonality. |
| 3427 | 29 ## |
|
11470
eb9e0b597d61
Use common names for variables in documentation and code for a few more m-script files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
30 ## If @var{V} is a vector, then @var{h} contains the Hessenberg matrix |
|
20370
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
31 ## such that @nospell{@tcode{a*u == u*h+rk*ek'}}, in which |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
32 ## @code{rk = a*u(:,k)-u*h(:,k)}, and @nospell{@tcode{ek'}} is the vector |
| 7248 | 33 ## @code{[0, 0, @dots{}, 1]} of length @code{k}. Otherwise, @var{h} is |
| 34 ## meaningless. | |
| 35 ## | |
|
20370
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
36 ## If @var{V} is a vector and @var{k} is greater than @code{length (A) - 1}, |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
37 ## then @var{h} contains the Hessenberg matrix such that @code{a*u == u*h}. |
| 5555 | 38 ## |
|
20370
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
39 ## The value of @var{nu} is the dimension of the span of the Krylov subspace |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
40 ## (based on @var{eps1}). |
| 5555 | 41 ## |
|
20370
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
42 ## If @var{b} is a vector and @var{k} is greater than @var{m-1}, then @var{h} |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
43 ## contains the Hessenberg decomposition of @var{A}. |
| 5555 | 44 ## |
|
20370
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
45 ## The optional parameter @var{eps1} is the threshold for zero. The default |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
46 ## value is 1e-12. |
| 5555 | 47 ## |
|
20370
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
48 ## If the optional parameter @var{pflg} is nonzero, row pivoting is used to |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
49 ## improve numerical behavior. The default value is 0. |
| 3427 | 50 ## |
|
19232
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
19047
diff
changeset
|
51 ## Reference: @nospell{A. Hodel, P. Misra}, @cite{Partial Pivoting in the |
|
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
19047
diff
changeset
|
52 ## Computation of Krylov Subspaces of Large Sparse Systems}, Proceedings of |
|
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
19047
diff
changeset
|
53 ## the 42nd IEEE Conference on Decision and Control, December 2003. |
| 3439 | 54 ## @end deftypefn |
| 55 | |
| 56 ## Author: A. Scottedward Hodel <a.s.hodel@eng.auburn.edu> | |
| 3240 | 57 |
| 4609 | 58 function [Uret, H, nu] = krylov (A, V, k, eps1, pflg); |
| 3240 | 59 |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
60 if (isa (A, "single") || isa (V, "single")) |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
61 defeps = 1e-6; |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
62 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
63 defeps = 1e-12; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
64 endif |
| 3457 | 65 |
| 66 if (nargin < 3 || nargin > 5) | |
| 6046 | 67 print_usage (); |
| 3457 | 68 elseif (nargin < 5) |
| 8506 | 69 ## Default permutation flag. |
| 70 pflg = 0; | |
| 3240 | 71 endif |
| 3457 | 72 |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
73 if (nargin < 4) |
| 8506 | 74 ## Default tolerance parameter. |
| 75 eps1 = defeps; | |
| 3240 | 76 endif |
| 3273 | 77 |
| 3457 | 78 if (isempty (eps1)) |
| 79 eps1 = defeps; | |
| 80 endif | |
| 81 | |
|
9872
72d6e0de76c7
fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents:
9843
diff
changeset
|
82 if (! issquare (A) || isempty (A)) |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
83 error ("krylov: A(%d x %d) must be a non-empty square matrix", rows (A), columns (A)); |
| 3240 | 84 endif |
|
9872
72d6e0de76c7
fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents:
9843
diff
changeset
|
85 na = rows (A); |
| 3225 | 86 |
| 7201 | 87 [m, kb] = size (V); |
| 3457 | 88 if (m != na) |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
89 error ("krylov: A(%d x %d), V(%d x %d): argument dimensions do not match", |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
90 na, na, m, kb); |
| 3240 | 91 endif |
| 3233 | 92 |
| 4030 | 93 if (! isscalar (k)) |
|
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11470
diff
changeset
|
94 error ("krylov: K must be a scalar integer"); |
| 3457 | 95 endif |
| 96 | |
| 97 Vnrm = norm (V, Inf); | |
| 3273 | 98 |
| 8506 | 99 ## check for trivial solution. |
| 3457 | 100 if (Vnrm == 0) |
| 101 Uret = []; | |
| 4609 | 102 H = []; |
| 3457 | 103 nu = 0; |
| 104 return; | |
| 3211 | 105 endif |
| 106 | |
| 8506 | 107 ## Identify trivial null space. |
| 3457 | 108 abm = max (abs ([A, V]')); |
| 109 zidx = find (abm == 0); | |
| 3240 | 110 |
| 8506 | 111 ## Set up vector of pivot points. |
| 3240 | 112 pivot_vec = 1:na; |
| 3225 | 113 |
| 3240 | 114 iter = 0; |
| 3273 | 115 alpha = []; |
| 3240 | 116 nh = 0; |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
117 while (length (alpha) < na) && (columns (V) > 0) && (iter < k) |
| 3240 | 118 iter++; |
| 3211 | 119 |
| 8506 | 120 ## Get orthogonal basis of V. |
| 3240 | 121 jj = 1; |
| 3457 | 122 while (jj <= columns (V) && length (alpha) < na) |
| 8506 | 123 ## Index of next Householder reflection. |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
124 nu = length (alpha)+1; |
| 3273 | 125 |
| 3240 | 126 short_pv = pivot_vec(nu:na); |
| 127 q = V(:,jj); | |
| 128 short_q = q(short_pv); | |
| 3211 | 129 |
| 3457 | 130 if (norm (short_q) < eps1) |
| 10549 | 131 ## Insignificant column; delete. |
| 3457 | 132 nv = columns (V); |
| 133 if (jj != nv) | |
| 134 [V(:,jj), V(:,nv)] = swap (V(:,jj), V(:,nv)); | |
|
19047
7bbe3658c5ef
maint: Use "FIXME:" coding convention in m-files.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
135 ## FIXME: H columns should be swapped too. |
|
7bbe3658c5ef
maint: Use "FIXME:" coding convention in m-files.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
136 ## Not done since Block Hessenberg structure is lost anyway. |
| 3240 | 137 endif |
| 138 V = V(:,1:(nv-1)); | |
| 10549 | 139 ## One less reflection. |
| 3457 | 140 nu--; |
| 3240 | 141 else |
| 10549 | 142 ## New householder reflection. |
| 3457 | 143 if (pflg) |
| 8506 | 144 ## Locate max magnitude element in short_q. |
| 3457 | 145 asq = abs (short_q); |
| 146 maxv = max (asq); | |
| 6024 | 147 maxidx = find (asq == maxv, 1); |
| 148 pivot_idx = short_pv(maxidx); | |
| 3273 | 149 |
| 10549 | 150 ## See if need to change the pivot list. |
| 3457 | 151 if (pivot_idx != pivot_vec(nu)) |
| 6024 | 152 swapidx = maxidx + (nu-1); |
| 3457 | 153 [pivot_vec(nu), pivot_vec(swapidx)] = ... |
| 10549 | 154 swap (pivot_vec(nu), pivot_vec(swapidx)); |
| 3240 | 155 endif |
| 156 endif | |
| 3273 | 157 |
| 10549 | 158 ## Isolate portion of vector for reflection. |
| 3240 | 159 idx = pivot_vec(nu:na); |
| 160 jdx = pivot_vec(1:nu); | |
| 161 | |
| 3457 | 162 [hv, av, z] = housh (q(idx), 1, 0); |
| 3240 | 163 alpha(nu) = av; |
| 164 U(idx,nu) = hv; | |
| 3211 | 165 |
| 8506 | 166 ## Reduce V per the reflection. |
| 3240 | 167 V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:)); |
| 168 if(iter > 1) | |
|
19047
7bbe3658c5ef
maint: Use "FIXME:" coding convention in m-files.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
169 ## FIXME: not done correctly for block case. |
| 3240 | 170 H(nu,nu-1) = V(pivot_vec(nu),jj); |
| 171 endif | |
| 172 | |
| 8506 | 173 ## Advance to next column of V. |
| 3457 | 174 jj++; |
| 3240 | 175 endif |
| 176 endwhile | |
| 3211 | 177 |
| 8506 | 178 ## Check for oversize V (due to full rank). |
| 3457 | 179 if ((columns (V) > na) && (length (alpha) == na)) |
| 8506 | 180 ## Trim to size. |
| 3273 | 181 V = V(:,1:na); |
| 3457 | 182 elseif (columns(V) > na) |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
183 krylov_V = V; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
184 krylov_na = na; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
185 krylov_length_alpha = length (alpha); |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
186 error ("krylov: this case should never happen; submit a bug report"); |
| 3273 | 187 endif |
| 188 | |
| 3457 | 189 if (columns (V) > 0) |
| 8506 | 190 ## Construct next Q and multiply. |
| 3457 | 191 Q = zeros (size (V)); |
| 192 for kk = 1:columns (Q) | |
| 3240 | 193 Q(pivot_vec(nu-columns(Q)+kk),kk) = 1; |
| 194 endfor | |
| 3273 | 195 |
| 8506 | 196 ## Apply Householder reflections. |
| 3240 | 197 for ii = nu:-1:1 |
| 198 idx = pivot_vec(ii:na); | |
| 199 hv = U(idx,ii); | |
| 200 av = alpha(ii); | |
| 201 Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:)); | |
| 202 endfor | |
| 3211 | 203 endif |
| 204 | |
| 8506 | 205 ## Multiply to get new vector. |
| 3240 | 206 V = A*Q; |
| 8506 | 207 ## Project off of previous vectors. |
| 3457 | 208 nu = length (alpha); |
| 209 for i = 1:nu | |
| 210 hv = U(:,i); | |
| 211 av = alpha(i); | |
|
20441
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20370
diff
changeset
|
212 V -= av*hv*(hv'*V); |
| 3240 | 213 H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:); |
| 7151 | 214 endfor |
| 3240 | 215 |
| 3211 | 216 endwhile |
| 217 | |
| 8506 | 218 ## Back out complete U matrix. |
| 219 ## back out U matrix. | |
| 3457 | 220 j1 = columns (U); |
| 221 for i = j1:-1:1; | |
| 3240 | 222 idx = pivot_vec(i:na); |
| 223 hv = U(idx,i); | |
| 224 av = alpha(i); | |
| 3457 | 225 U(:,i) = zeros (na, 1); |
| 3240 | 226 U(idx(1),i) = 1; |
| 227 U(idx,i:j1) = U(idx,i:j1)-av*hv*(hv'*U(idx,i:j1)); | |
| 3211 | 228 endfor |
| 3273 | 229 |
| 3457 | 230 nu = length (alpha); |
| 3240 | 231 Uret = U; |
| 3457 | 232 if (max (max (abs (Uret(zidx,:)))) > 0) |
| 233 warning ("krylov: trivial null space corrupted; set pflg = 1 or eps1 > %e", | |
| 10549 | 234 eps1); |
| 3225 | 235 endif |
| 236 | |
| 3211 | 237 endfunction |
| 9843 | 238 |
| 239 | |
| 240 function [a1, b1] = swap (a, b) | |
| 241 | |
| 242 a1 = b; | |
| 243 b1 = a; | |
| 244 | |
| 245 endfunction | |
|
17338
1c89599167a6
maint: End m-files with 1 blank line.
Rik <rik@octave.org>
parents:
17268
diff
changeset
|
246 |