Mercurial > hg > octave-lyh

comparison scripts/linear-algebra/krylov.m @ 8506:bc982528de11

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
comment style fixes
author John W. Eaton <jwe@octave.org>
date Tue, 13 Jan 2009 11:56:00 -0500
parents df9519e9990c
children eb63fbe60fab
comparison
equal deleted inserted replaced
8505:124dd27bedae 8506:bc982528de11
66 endif 66 endif
67 67
68 if (nargin < 3 || nargin > 5) 68 if (nargin < 3 || nargin > 5)
69 print_usage (); 69 print_usage ();
70 elseif (nargin < 5) 70 elseif (nargin < 5)
71 pflg = 0; # default permutation flag 71 ## Default permutation flag.
72 pflg = 0;
72 endif 73 endif
73 74
74 if(nargin < 4) 75 if(nargin < 4)
75 eps1 = defeps; # default tolerance parameter 76 ## Default tolerance parameter.
77 eps1 = defeps;
76 endif 78 endif
77 79
78 if (isempty (eps1)) 80 if (isempty (eps1))
79 eps1 = defeps; 81 eps1 = defeps;
80 endif 82 endif
94 error ("krylov: third argument must be a scalar integer"); 96 error ("krylov: third argument must be a scalar integer");
95 endif 97 endif
96 98
97 Vnrm = norm (V, Inf); 99 Vnrm = norm (V, Inf);
98 100
99 ## check for trivial solution 101 ## check for trivial solution.
100 if (Vnrm == 0) 102 if (Vnrm == 0)
101 Uret = []; 103 Uret = [];
102 H = []; 104 H = [];
103 nu = 0; 105 nu = 0;
104 return; 106 return;
105 endif 107 endif
106 108
107 # identify trivial null space 109 ## Identify trivial null space.
108 abm = max (abs ([A, V]')); 110 abm = max (abs ([A, V]'));
109 zidx = find (abm == 0); 111 zidx = find (abm == 0);
110 112
111 # set up vector of pivot points 113 ## Set up vector of pivot points.
112 pivot_vec = 1:na; 114 pivot_vec = 1:na;
113 115
114 iter = 0; 116 iter = 0;
115 alpha = []; 117 alpha = [];
116 nh = 0; 118 nh = 0;
117 while (length(alpha) < na) && (columns(V) > 0) && (iter < k) 119 while (length(alpha) < na) && (columns(V) > 0) && (iter < k)
118 iter++; 120 iter++;
119 121
120 ## get orthogonal basis of V 122 ## Get orthogonal basis of V.
121 jj = 1; 123 jj = 1;
122 while (jj <= columns (V) && length (alpha) < na) 124 while (jj <= columns (V) && length (alpha) < na)
123 ## index of next Householder reflection 125 ## Index of next Householder reflection.
124 nu = length(alpha)+1; 126 nu = length(alpha)+1;
125 127
126 short_pv = pivot_vec(nu:na); 128 short_pv = pivot_vec(nu:na);
127 q = V(:,jj); 129 q = V(:,jj);
128 short_q = q(short_pv); 130 short_q = q(short_pv);
129 131
130 if (norm (short_q) < eps1) 132 if (norm (short_q) < eps1)
131 ## insignificant column; delete 133 ## Insignificant column; delete.
132 nv = columns (V); 134 nv = columns (V);
133 if (jj != nv) 135 if (jj != nv)
134 [V(:,jj), V(:,nv)] = swap (V(:,jj), V(:,nv)); 136 [V(:,jj), V(:,nv)] = swap (V(:,jj), V(:,nv));
135 ## FIXME -- H columns should be swapped too. Not done 137 ## FIXME -- H columns should be swapped too. Not done
136 ## since Block Hessenberg structure is lost anyway. 138 ## since Block Hessenberg structure is lost anyway.
137 endif 139 endif
138 V = V(:,1:(nv-1)); 140 V = V(:,1:(nv-1));
139 ## one less reflection 141 ## One less reflection.
140 nu--; 142 nu--;
141 else 143 else
142 ## new householder reflection 144 ## New householder reflection.
143 if (pflg) 145 if (pflg)
144 ## locate max magnitude element in short_q 146 ## Locate max magnitude element in short_q.
145 asq = abs (short_q); 147 asq = abs (short_q);
146 maxv = max (asq); 148 maxv = max (asq);
147 maxidx = find (asq == maxv, 1); 149 maxidx = find (asq == maxv, 1);
148 pivot_idx = short_pv(maxidx); 150 pivot_idx = short_pv(maxidx);
149 151
150 ## see if need to change the pivot list 152 ## See if need to change the pivot list.
151 if (pivot_idx != pivot_vec(nu)) 153 if (pivot_idx != pivot_vec(nu))
152 swapidx = maxidx + (nu-1); 154 swapidx = maxidx + (nu-1);
153 [pivot_vec(nu), pivot_vec(swapidx)] = ... 155 [pivot_vec(nu), pivot_vec(swapidx)] = ...
154 swap (pivot_vec(nu), pivot_vec(swapidx)); 156 swap (pivot_vec(nu), pivot_vec(swapidx));
155 endif 157 endif
156 endif 158 endif
157 159
158 ## isolate portion of vector for reflection 160 ## Isolate portion of vector for reflection.
159 idx = pivot_vec(nu:na); 161 idx = pivot_vec(nu:na);
160 jdx = pivot_vec(1:nu); 162 jdx = pivot_vec(1:nu);
161 163
162 [hv, av, z] = housh (q(idx), 1, 0); 164 [hv, av, z] = housh (q(idx), 1, 0);
163 alpha(nu) = av; 165 alpha(nu) = av;
164 U(idx,nu) = hv; 166 U(idx,nu) = hv;
165 167
166 # reduce V per the reflection 168 ## Reduce V per the reflection.
167 V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:)); 169 V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:));
168 if(iter > 1) 170 if(iter > 1)
169 ## FIXME -- not done correctly for block case 171 ## FIXME -- not done correctly for block case.
170 H(nu,nu-1) = V(pivot_vec(nu),jj); 172 H(nu,nu-1) = V(pivot_vec(nu),jj);
171 endif 173 endif
172 174
173 ## advance to next column of V 175 ## Advance to next column of V.
174 jj++; 176 jj++;
175 endif 177 endif
176 endwhile 178 endwhile
177 179
178 ## check for oversize V (due to full rank) 180 ## Check for oversize V (due to full rank).
179 if ((columns (V) > na) && (length (alpha) == na)) 181 if ((columns (V) > na) && (length (alpha) == na))
180 ## trim to size 182 ## Trim to size.
181 V = V(:,1:na); 183 V = V(:,1:na);
182 elseif (columns(V) > na) 184 elseif (columns(V) > na)
183 krylov_V = V 185 krylov_V = V
184 krylov_na = na 186 krylov_na = na
185 krylov_length_alpha = length (alpha) 187 krylov_length_alpha = length (alpha)
186 error ("This case should never happen; submit a bug report"); 188 error ("This case should never happen; submit a bug report");
187 endif 189 endif
188 190
189 if (columns (V) > 0) 191 if (columns (V) > 0)
190 ## construct next Q and multiply 192 ## Construct next Q and multiply.
191 Q = zeros (size (V)); 193 Q = zeros (size (V));
192 for kk = 1:columns (Q) 194 for kk = 1:columns (Q)
193 Q(pivot_vec(nu-columns(Q)+kk),kk) = 1; 195 Q(pivot_vec(nu-columns(Q)+kk),kk) = 1;
194 endfor 196 endfor
195 197
196 ## apply Householder reflections 198 ## Apply Householder reflections.
197 for ii = nu:-1:1 199 for ii = nu:-1:1
198 idx = pivot_vec(ii:na); 200 idx = pivot_vec(ii:na);
199 hv = U(idx,ii); 201 hv = U(idx,ii);
200 av = alpha(ii); 202 av = alpha(ii);
201 Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:)); 203 Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:));
202 endfor 204 endfor
203 endif 205 endif
204 206
205 ## multiply to get new vector; 207 ## Multiply to get new vector.
206 V = A*Q; 208 V = A*Q;
207 ## project off of previous vectors 209 ## Project off of previous vectors.
208 nu = length (alpha); 210 nu = length (alpha);
209 for i = 1:nu 211 for i = 1:nu
210 hv = U(:,i); 212 hv = U(:,i);
211 av = alpha(i); 213 av = alpha(i);
212 V = V - av*hv*(hv'*V); 214 V = V - av*hv*(hv'*V);
213 H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:); 215 H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:);
214 endfor 216 endfor
215 217
216 endwhile 218 endwhile
217 219
218 ## Back out complete U matrix 220 ## Back out complete U matrix.
219 ## back out U matrix ; 221 ## back out U matrix.
220 j1 = columns (U); 222 j1 = columns (U);
221 for i = j1:-1:1; 223 for i = j1:-1:1;
222 idx = pivot_vec(i:na); 224 idx = pivot_vec(i:na);
223 hv = U(idx,i); 225 hv = U(idx,i);
224 av = alpha(i); 226 av = alpha(i);