|
3427
|
1 ## Copyright (C) 1993, 1998, 1999 Auburn University. All rights reserved. |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by the |
|
|
7 ## Free Software Foundation; either version 2, or (at your option) any |
|
|
8 ## later version. |
|
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
13 ## for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
|
17 ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. |
|
|
18 |
|
|
19 ## function [U,H,nu] = krylov(A,V,k{,eps1,pflg}); |
|
|
20 ## construct orthogonal basis U of block Krylov subspace; |
|
|
21 ## [V AV A^2*V ... A^(k+1)*V]; |
|
|
22 ## method used: householder reflections to guard against loss of |
|
|
23 ## orthogonality |
|
|
24 ## eps1: threshhold for 0 (default: 1e-12) |
|
|
25 ## pflg: flag to use row pivoting (improves numerical behavior) |
|
|
26 ## 0 [default]: no pivoting; prints a warning message if trivial |
|
|
27 ## null space is corrupted |
|
|
28 ## 1 : pivoting performed |
|
|
29 ## |
|
|
30 ## outputs: |
|
|
31 ## Uret: orthogonal basis of block krylov subspace |
|
|
32 ## H: Hessenberg matrix; if V is a vector then A U = U H |
|
|
33 ## otherwise H is meaningless |
|
|
34 ## nu: dimension of span of krylov subspace (based on eps1) |
|
|
35 ## if B is a vector and k > m-1, krylov returns H = the Hessenberg |
|
|
36 ## decompostion of A. |
|
|
37 ## |
|
|
38 ## Reference: Hodel and Misra, "Partial Pivoting in the Computation of |
|
|
39 ## Krylov Subspaces", to be submitted to Linear Algebra and its |
|
|
40 ## Applications |
|
|
41 ## written by A. Scottedward Hodel a.s.hodel@eng.auburn.edu |
|
3240
|
42 |
|
|
43 function [Uret,H,nu] = krylov(A,V,k,eps1,pflg); |
|
|
44 |
|
|
45 defeps = 1e-12; |
|
|
46 if(nargin < 3 | nargin > 5) |
|
|
47 usage("[U,nu] = krylov(A,V,k{,eps1,pflg})") |
|
|
48 elseif(nargin < 5) |
|
|
49 pflg = 0; # default permutation flag |
|
|
50 endif |
|
|
51 if(nargin < 4) |
|
|
52 eps1 = defeps; # default tolerance parameter |
|
|
53 endif |
|
|
54 if(isempty(eps1)) eps1 = defeps; endif |
|
3273
|
55 |
|
3240
|
56 na = is_square(A); |
|
|
57 if( !na ) error("A(%d x %d) must be square",rows(A),columns(A)); endif |
|
3273
|
58 |
|
|
59 [m,kb] = size(V); |
|
3240
|
60 if(m != na); |
|
|
61 error("A(%d x %d), V(%d x %d): argument dimensions do not match", ... |
|
|
62 na,na,m,kb) |
|
|
63 endif |
|
3225
|
64 |
|
3240
|
65 if( !is_scalar(k) ) |
|
|
66 error("krylov: third argument must be a scalar integer") |
|
|
67 endif |
|
3233
|
68 |
|
3240
|
69 Vnrm = norm(V,Inf); |
|
3273
|
70 |
|
3240
|
71 # check for trivial solution |
|
|
72 if(Vnrm == 0) |
|
|
73 Uret = []; nu = 0; return; |
|
3211
|
74 endif |
|
|
75 |
|
3240
|
76 # identify trivial null space |
|
|
77 abm = max(abs([A,V]')); nzidx = find(abm != 0); zidx = find(abm == 0); |
|
|
78 |
|
|
79 # set up vector of pivot points |
|
|
80 pivot_vec = 1:na; |
|
3225
|
81 |
|
3240
|
82 iter = 0; |
|
3273
|
83 alpha = []; |
|
3240
|
84 nh = 0; |
|
|
85 while (length(alpha) < na) & (columns(V) > 0) & (iter < k) |
|
|
86 iter++; |
|
3211
|
87 |
|
3240
|
88 # get orthogonal basis of V |
|
|
89 jj = 1; |
|
|
90 while(jj <= columns(V) & length(alpha) < na) |
|
|
91 nu = length(alpha)+1; # index of next Householder reflection |
|
3273
|
92 |
|
3240
|
93 short_pv = pivot_vec(nu:na); |
|
|
94 q = V(:,jj); |
|
|
95 short_q = q(short_pv); |
|
3211
|
96 |
|
3240
|
97 if(norm(short_q) < eps1) |
|
|
98 # insignificant column; delete |
|
|
99 nv = columns(V); |
|
3273
|
100 if(jj != nv) |
|
|
101 [V(:,jj),V(:,nv)] = swap(V(:,jj),V(:,nv)); |
|
3240
|
102 # FIX ME: H columns should be swapped too. Not done since |
|
|
103 # Block Hessenberg structure is lost anyway. |
|
|
104 endif |
|
|
105 V = V(:,1:(nv-1)); |
|
|
106 nu = nu - 1; # one less reflection |
|
3211
|
107 |
|
3240
|
108 else |
|
|
109 # new householder reflection |
|
|
110 if(pflg) |
|
|
111 # locate max magnitude element in short_q |
|
|
112 asq = abs(short_q); |
|
3426
|
113 maxv = max(asq); |
|
3240
|
114 maxidx = find(asq == maxv); |
|
|
115 pivot_idx = short_pv(maxidx(1)); |
|
3273
|
116 |
|
3240
|
117 # see if need to change the pivot list |
|
|
118 if(pivot_idx != pivot_vec(nu)) |
|
|
119 swapidx = maxidx(1) + (nu-1); |
|
|
120 [pivot_vec(nu),pivot_vec(swapidx)] = ... |
|
|
121 swap(pivot_vec(nu),pivot_vec(swapidx)); |
|
|
122 endif |
|
|
123 endif |
|
3273
|
124 |
|
3240
|
125 # isolate portion of vector for reflection |
|
|
126 idx = pivot_vec(nu:na); |
|
|
127 jdx = pivot_vec(1:nu); |
|
|
128 |
|
|
129 [hv,av,z] = housh(q(idx),1,0); |
|
|
130 alpha(nu) = av; |
|
|
131 U(idx,nu) = hv; |
|
3211
|
132 |
|
3240
|
133 # reduce V per the reflection |
|
|
134 V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:)); |
|
|
135 if(iter > 1) |
|
|
136 # FIX ME: not done correctly for block case |
|
|
137 H(nu,nu-1) = V(pivot_vec(nu),jj); |
|
|
138 endif |
|
|
139 |
|
|
140 # advance to next column of V |
|
|
141 jj=jj+1; |
|
|
142 endif |
|
|
143 endwhile |
|
3211
|
144 |
|
3273
|
145 # check for oversize V (due to full rank) |
|
|
146 if( ( columns(V) > na ) & ( length(alpha) == na ) ) |
|
|
147 # trim to size |
|
|
148 V = V(:,1:na); |
|
|
149 elseif( columns(V) > na ) |
|
|
150 krylov_V = V |
|
|
151 krylov_na = na |
|
|
152 krylov_length_alpha = length(alpha) |
|
|
153 error("This case should never happen; submit bug report."); |
|
|
154 endif |
|
|
155 |
|
3240
|
156 if(columns(V) > 0) |
|
|
157 # construct next Q and multiply |
|
|
158 Q = zeros(size(V)); |
|
|
159 for kk=1:columns(Q) |
|
|
160 Q(pivot_vec(nu-columns(Q)+kk),kk) = 1; |
|
|
161 endfor |
|
3273
|
162 |
|
3240
|
163 # apply Householder reflections |
|
|
164 for ii = nu:-1:1 |
|
|
165 idx = pivot_vec(ii:na); |
|
|
166 hv = U(idx,ii); |
|
|
167 av = alpha(ii); |
|
|
168 Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:)); |
|
|
169 endfor |
|
3211
|
170 endif |
|
|
171 |
|
3240
|
172 # multiply to get new vector; |
|
|
173 V = A*Q; |
|
|
174 # project off of previous vectors |
|
|
175 nu = length(alpha); |
|
|
176 for i=1:nu |
|
|
177 hv = U(:,i); av = alpha(i); |
|
|
178 V = V - av*hv*(hv'*V); |
|
|
179 H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:); |
|
|
180 end |
|
|
181 |
|
3211
|
182 endwhile |
|
|
183 |
|
3240
|
184 # Back out complete U matrix |
|
3211
|
185 # back out U matrix ; |
|
|
186 j1 = columns(U); |
|
|
187 for i=j1:-1:1; |
|
3240
|
188 idx = pivot_vec(i:na); |
|
|
189 hv = U(idx,i); |
|
|
190 av = alpha(i); |
|
|
191 U(:,i) = zeros(na,1); |
|
|
192 U(idx(1),i) = 1; |
|
|
193 U(idx,i:j1) = U(idx,i:j1)-av*hv*(hv'*U(idx,i:j1)); |
|
3211
|
194 endfor |
|
3273
|
195 |
|
3240
|
196 nu = length(alpha); |
|
|
197 Uret = U; |
|
|
198 if( max(max( abs(Uret(zidx,:)) )) > 0) |
|
3225
|
199 warning("krylov: trivial null space corrupted; set pflg=1 or eps1>%e",eps1); |
|
|
200 endif |
|
|
201 |
|
3211
|
202 endfunction |
