Mercurial > hg > octave-lyh
annotate scripts/linear-algebra/krylov.m @ 14868:5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
* lin2mu.m, loadaudio.m, wavread.m, accumarray.m, bicubic.m, celldisp.m,
colon.m, cplxpair.m, dblquad.m, divergence.m, genvarname.m, gradient.m,
int2str.m, interp1.m, interp1q.m, interp2.m, interpn.m, loadobj.m, nthargout.m,
__isequal__.m, __splinen__.m, quadgk.m, quadl.m, quadv.m, rat.m, rot90.m,
rotdim.m, saveobj.m, subsindex.m, triplequad.m, delaunay3.m, griddata.m,
inpolygon.m, tsearchn.m, voronoi.m, get_first_help_sentence.m, which.m,
gray2ind.m, pink.m, dlmwrite.m, strread.m, textread.m, textscan.m, housh.m,
ishermitian.m, issymmetric.m, krylov.m, logm.m, null.m, rref.m,
compare_versions.m, copyfile.m, dump_prefs.m, edit.m, fileparts.m,
getappdata.m, isappdata.m, movefile.m, orderfields.m, parseparams.m,
__xzip__.m, rmappdata.m, setappdata.m, swapbytes.m, unpack.m, ver.m, fminbnd.m,
fminunc.m, fsolve.m, glpk.m, lsqnonneg.m, qp.m, sqp.m, configure_make.m,
copy_files.m, describe.m, get_description.m, get_forge_pkg.m, install.m,
installed_packages.m, is_architecture_dependent.m, load_package_dirs.m,
print_package_description.m, rebuild.m, repackage.m, save_order.m, shell.m,
allchild.m, ancestor.m, area.m, axes.m, axis.m, clabel.m, close.m, colorbar.m,
comet.m, comet3.m, contour.m, cylinder.m, ezmesh.m, ezsurf.m, findobj.m,
fplot.m, hist.m, isocolors.m, isonormals.m, isosurface.m, isprop.m, legend.m,
mesh.m, meshz.m, pareto.m, pcolor.m, peaks.m, plot3.m, plotmatrix.m, plotyy.m,
polar.m, print.m, __add_datasource__.m, __add_default_menu__.m,
__axes_limits__.m, __bar__.m, __clabel__.m, __contour__.m, __errcomm__.m,
__errplot__.m, __ezplot__.m, __file_filter__.m, __fltk_print__.m,
__ghostscript__.m, __gnuplot_print__.m, __go_draw_axes__.m,
__go_draw_figure__.m, __interp_cube__.m, __marching_cube__.m, __patch__.m,
__pie__.m, __plt__.m, __print_parse_opts__.m, __quiver__.m, __scatter__.m,
__stem__.m, __tight_eps_bbox__.m, __uigetdir_fltk__.m, __uigetfile_fltk__.m,
__uiputfile_fltk__.m, quiver.m, quiver3.m, rectangle.m, refreshdata.m,
ribbon.m, scatter.m, semilogy.m, shading.m, slice.m, subplot.m, surface.m,
surfl.m, surfnorm.m, text.m, uigetfile.m, uiputfile.m, whitebg.m, deconv.m,
mkpp.m, pchip.m, polyaffine.m, polyder.m, polygcd.m, polyout.m, polyval.m,
ppint.m, ppjumps.m, ppval.m, residue.m, roots.m, spline.m, splinefit.m,
addpref.m, getpref.m, setpref.m, ismember.m, setxor.m, arch_fit.m, arch_rnd.m,
arch_test.m, autoreg_matrix.m, diffpara.m, fftconv.m, filter2.m, hanning.m,
hurst.m, periodogram.m, triangle_sw.m, sinc.m, spectral_xdf.m, spencer.m,
stft.m, synthesis.m, unwrap.m, yulewalker.m, bicgstab.m, gmres.m, pcg.m, pcr.m,
__sprand_impl__.m, speye.m, spfun.m, sprandn.m, spstats.m, svds.m,
treelayout.m, treeplot.m, bessel.m, factor.m, legendre.m, perms.m, primes.m,
magic.m, toeplitz.m, corr.m, cov.m, mean.m, median.m, mode.m, qqplot.m,
quantile.m, ranks.m, zscore.m, logistic_regression_likelihood.m,
bartlett_test.m, chisquare_test_homogeneity.m, chisquare_test_independence.m,
kolmogorov_smirnov_test.m, run_test.m, u_test.m, wilcoxon_test.m, z_test.m,
z_test_2.m, bin2dec.m, dec2base.m, mat2str.m, strcat.m, strchr.m, strjust.m,
strtok.m, substr.m, untabify.m, assert.m, demo.m, example.m, fail.m, speed.m,
test.m, now.m: Use Octave coding conventions for cuddling parentheses in
scripts directory.
| author | Rik <octave@nomad.inbox5.com> |
|---|---|
| date | Tue, 17 Jul 2012 07:08:39 -0700 |
| parents | 72c96de7a403 |
| children | 1c21f264d26f |
| rev | line source |
|---|---|
|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
13931
diff
changeset
|
1 ## Copyright (C) 1993-2012 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 | |
| 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 |
|
13931
9de488c6c59c
doc: Spellcheck documentation before 3.6.0 release
Rik <octave@nomad.inbox5.com>
parents:
13780
diff
changeset
|
31 ## such that @nospell{@xcode{a*u == u*h+rk*ek'}}, in which @code{rk = |
|
9de488c6c59c
doc: Spellcheck documentation before 3.6.0 release
Rik <octave@nomad.inbox5.com>
parents:
13780
diff
changeset
|
32 ## a*u(:,k)-u*h(:,k)}, and @nospell{@xcode{ek'}} is the vector |
| 7248 | 33 ## @code{[0, 0, @dots{}, 1]} of length @code{k}. Otherwise, @var{h} is |
| 34 ## meaningless. | |
| 35 ## | |
|
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
|
36 ## If @var{V} is a vector and @var{k} is greater than |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
37 ## @code{length (A) - 1}, then @var{h} contains the Hessenberg matrix such |
| 7248 | 38 ## that @code{a*u == u*h}. |
| 5555 | 39 ## |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
40 ## The value of @var{nu} is the dimension of the span of the Krylov |
| 5555 | 41 ## subspace (based on @var{eps1}). |
| 42 ## | |
| 43 ## If @var{b} is a vector and @var{k} is greater than @var{m-1}, then | |
|
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
|
44 ## @var{h} contains the Hessenberg decomposition of @var{A}. |
| 5555 | 45 ## |
| 46 ## The optional parameter @var{eps1} is the threshold for zero. The | |
| 47 ## default value is 1e-12. | |
| 48 ## | |
| 49 ## If the optional parameter @var{pflg} is nonzero, row pivoting is used | |
| 50 ## to improve numerical behavior. The default value is 0. | |
| 3427 | 51 ## |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
52 ## Reference: A. Hodel, P. Misra, @cite{Partial Pivoting in the Computation of |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
53 ## Krylov Subspaces of Large Sparse Systems}, Proceedings of the 42nd IEEE |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
54 ## Conference on Decision and Control, December 2003. |
| 3439 | 55 ## @end deftypefn |
| 56 | |
| 57 ## Author: A. Scottedward Hodel <a.s.hodel@eng.auburn.edu> | |
| 3240 | 58 |
| 4609 | 59 function [Uret, H, nu] = krylov (A, V, k, eps1, pflg); |
| 3240 | 60 |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
61 if (isa (A, "single") || isa (V, "single")) |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
62 defeps = 1e-6; |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
63 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
64 defeps = 1e-12; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
65 endif |
| 3457 | 66 |
| 67 if (nargin < 3 || nargin > 5) | |
| 6046 | 68 print_usage (); |
| 3457 | 69 elseif (nargin < 5) |
| 8506 | 70 ## Default permutation flag. |
| 71 pflg = 0; | |
| 3240 | 72 endif |
| 3457 | 73 |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
74 if (nargin < 4) |
| 8506 | 75 ## Default tolerance parameter. |
| 76 eps1 = defeps; | |
| 3240 | 77 endif |
| 3273 | 78 |
| 3457 | 79 if (isempty (eps1)) |
| 80 eps1 = defeps; | |
| 81 endif | |
| 82 | |
|
9872
72d6e0de76c7
fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents:
9843
diff
changeset
|
83 if (! issquare (A) || isempty (A)) |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
84 error ("krylov: A(%d x %d) must be a non-empty square matrix", rows (A), columns (A)); |
| 3240 | 85 endif |
|
9872
72d6e0de76c7
fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents:
9843
diff
changeset
|
86 na = rows (A); |
| 3225 | 87 |
| 7201 | 88 [m, kb] = size (V); |
| 3457 | 89 if (m != na) |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
90 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
|
91 na, na, m, kb); |
| 3240 | 92 endif |
| 3233 | 93 |
| 4030 | 94 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
|
95 error ("krylov: K must be a scalar integer"); |
| 3457 | 96 endif |
| 97 | |
| 98 Vnrm = norm (V, Inf); | |
| 3273 | 99 |
| 8506 | 100 ## check for trivial solution. |
| 3457 | 101 if (Vnrm == 0) |
| 102 Uret = []; | |
| 4609 | 103 H = []; |
| 3457 | 104 nu = 0; |
| 105 return; | |
| 3211 | 106 endif |
| 107 | |
| 8506 | 108 ## Identify trivial null space. |
| 3457 | 109 abm = max (abs ([A, V]')); |
| 110 zidx = find (abm == 0); | |
| 3240 | 111 |
| 8506 | 112 ## Set up vector of pivot points. |
| 3240 | 113 pivot_vec = 1:na; |
| 3225 | 114 |
| 3240 | 115 iter = 0; |
| 3273 | 116 alpha = []; |
| 3240 | 117 nh = 0; |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
118 while (length (alpha) < na) && (columns (V) > 0) && (iter < k) |
| 3240 | 119 iter++; |
| 3211 | 120 |
| 8506 | 121 ## Get orthogonal basis of V. |
| 3240 | 122 jj = 1; |
| 3457 | 123 while (jj <= columns (V) && length (alpha) < na) |
| 8506 | 124 ## 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
|
125 nu = length (alpha)+1; |
| 3273 | 126 |
| 3240 | 127 short_pv = pivot_vec(nu:na); |
| 128 q = V(:,jj); | |
| 129 short_q = q(short_pv); | |
| 3211 | 130 |
| 3457 | 131 if (norm (short_q) < eps1) |
| 10549 | 132 ## Insignificant column; delete. |
| 3457 | 133 nv = columns (V); |
| 134 if (jj != nv) | |
| 135 [V(:,jj), V(:,nv)] = swap (V(:,jj), V(:,nv)); | |
| 10549 | 136 ## FIXME -- H columns should be swapped too. Not done |
| 137 ## since Block Hessenberg structure is lost anyway. | |
| 3240 | 138 endif |
| 139 V = V(:,1:(nv-1)); | |
| 10549 | 140 ## One less reflection. |
| 3457 | 141 nu--; |
| 3240 | 142 else |
| 10549 | 143 ## New householder reflection. |
| 3457 | 144 if (pflg) |
| 8506 | 145 ## Locate max magnitude element in short_q. |
| 3457 | 146 asq = abs (short_q); |
| 147 maxv = max (asq); | |
| 6024 | 148 maxidx = find (asq == maxv, 1); |
| 149 pivot_idx = short_pv(maxidx); | |
| 3273 | 150 |
| 10549 | 151 ## See if need to change the pivot list. |
| 3457 | 152 if (pivot_idx != pivot_vec(nu)) |
| 6024 | 153 swapidx = maxidx + (nu-1); |
| 3457 | 154 [pivot_vec(nu), pivot_vec(swapidx)] = ... |
| 10549 | 155 swap (pivot_vec(nu), pivot_vec(swapidx)); |
| 3240 | 156 endif |
| 157 endif | |
| 3273 | 158 |
| 10549 | 159 ## Isolate portion of vector for reflection. |
| 3240 | 160 idx = pivot_vec(nu:na); |
| 161 jdx = pivot_vec(1:nu); | |
| 162 | |
| 3457 | 163 [hv, av, z] = housh (q(idx), 1, 0); |
| 3240 | 164 alpha(nu) = av; |
| 165 U(idx,nu) = hv; | |
| 3211 | 166 |
| 8506 | 167 ## Reduce V per the reflection. |
| 3240 | 168 V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:)); |
| 169 if(iter > 1) | |
| 10549 | 170 ## FIXME -- not done correctly for block case. |
| 3240 | 171 H(nu,nu-1) = V(pivot_vec(nu),jj); |
| 172 endif | |
| 173 | |
| 8506 | 174 ## Advance to next column of V. |
| 3457 | 175 jj++; |
| 3240 | 176 endif |
| 177 endwhile | |
| 3211 | 178 |
| 8506 | 179 ## Check for oversize V (due to full rank). |
| 3457 | 180 if ((columns (V) > na) && (length (alpha) == na)) |
| 8506 | 181 ## Trim to size. |
| 3273 | 182 V = V(:,1:na); |
| 3457 | 183 elseif (columns(V) > na) |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
184 krylov_V = V; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
185 krylov_na = na; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
186 krylov_length_alpha = length (alpha); |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
187 error ("krylov: this case should never happen; submit a bug report"); |
| 3273 | 188 endif |
| 189 | |
| 3457 | 190 if (columns (V) > 0) |
| 8506 | 191 ## Construct next Q and multiply. |
| 3457 | 192 Q = zeros (size (V)); |
| 193 for kk = 1:columns (Q) | |
| 3240 | 194 Q(pivot_vec(nu-columns(Q)+kk),kk) = 1; |
| 195 endfor | |
| 3273 | 196 |
| 8506 | 197 ## Apply Householder reflections. |
| 3240 | 198 for ii = nu:-1:1 |
| 199 idx = pivot_vec(ii:na); | |
| 200 hv = U(idx,ii); | |
| 201 av = alpha(ii); | |
| 202 Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:)); | |
| 203 endfor | |
| 3211 | 204 endif |
| 205 | |
| 8506 | 206 ## Multiply to get new vector. |
| 3240 | 207 V = A*Q; |
| 8506 | 208 ## Project off of previous vectors. |
| 3457 | 209 nu = length (alpha); |
| 210 for i = 1:nu | |
| 211 hv = U(:,i); | |
| 212 av = alpha(i); | |
| 3240 | 213 V = V - av*hv*(hv'*V); |
| 214 H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:); | |
| 7151 | 215 endfor |
| 3240 | 216 |
| 3211 | 217 endwhile |
| 218 | |
| 8506 | 219 ## Back out complete U matrix. |
| 220 ## back out U matrix. | |
| 3457 | 221 j1 = columns (U); |
| 222 for i = j1:-1:1; | |
| 3240 | 223 idx = pivot_vec(i:na); |
| 224 hv = U(idx,i); | |
| 225 av = alpha(i); | |
| 3457 | 226 U(:,i) = zeros (na, 1); |
| 3240 | 227 U(idx(1),i) = 1; |
| 228 U(idx,i:j1) = U(idx,i:j1)-av*hv*(hv'*U(idx,i:j1)); | |
| 3211 | 229 endfor |
| 3273 | 230 |
| 3457 | 231 nu = length (alpha); |
| 3240 | 232 Uret = U; |
| 3457 | 233 if (max (max (abs (Uret(zidx,:)))) > 0) |
| 234 warning ("krylov: trivial null space corrupted; set pflg = 1 or eps1 > %e", | |
| 10549 | 235 eps1); |
| 3225 | 236 endif |
| 237 | |
| 3211 | 238 endfunction |
| 9843 | 239 |
| 240 | |
| 241 function [a1, b1] = swap (a, b) | |
| 242 | |
| 243 a1 = b; | |
| 244 b1 = a; | |
| 245 | |
| 246 endfunction |