# | Line 193 | Line 193 | bool SquareMatrix<Real, Dim>::jacobi(const SquareMatri | |
---|---|---|
193 | Vector<Real, Dim> b, z; | |
194 | ||
195 | // initialize | |
196 | < | for (ip=0; ip<N; ip++) |
197 | < | { |
198 | < | for (iq=0; iq<N; iq++) v(ip, iq) = 0.0; |
199 | < | v(ip, ip) = 1.0; |
196 | > | for (ip=0; ip<N; ip++) { |
197 | > | for (iq=0; iq<N; iq++) |
198 | > | v(ip, iq) = 0.0; |
199 | > | v(ip, ip) = 1.0; |
200 | } | |
201 | < | for (ip=0; ip<N; ip++) |
202 | < | { |
203 | < | b(ip) = w(ip) = a(ip, ip); |
204 | < | z(ip) = 0.0; |
201 | > | |
202 | > | for (ip=0; ip<N; ip++) { |
203 | > | b(ip) = w(ip) = a(ip, ip); |
204 | > | z(ip) = 0.0; |
205 | } | |
206 | ||
207 | // begin rotation sequence | |
208 | < | for (i=0; i<MAX_ROTATIONS; i++) |
209 | < | { |
210 | < | sm = 0.0; |
211 | < | for (ip=0; ip<2; ip++) |
212 | < | { |
213 | < | for (iq=ip+1; iq<N; iq++) sm += fabs(a(ip, iq)); |
214 | < | } |
215 | < | if (sm == 0.0) break; |
208 | > | for (i=0; i<MAX_ROTATIONS; i++) { |
209 | > | sm = 0.0; |
210 | > | for (ip=0; ip<2; ip++) { |
211 | > | for (iq=ip+1; iq<N; iq++) |
212 | > | sm += fabs(a(ip, iq)); |
213 | > | } |
214 | > | |
215 | > | if (sm == 0.0) |
216 | > | break; |
217 | ||
218 | < | if (i < 4) tresh = 0.2*sm/(9); |
219 | < | else tresh = 0.0; |
218 | > | if (i < 4) |
219 | > | tresh = 0.2*sm/(9); |
220 | > | else |
221 | > | tresh = 0.0; |
222 | ||
223 | < | for (ip=0; ip<2; ip++) |
224 | < | { |
225 | < | for (iq=ip+1; iq<N; iq++) |
226 | < | { |
227 | < | g = 100.0*fabs(a(ip, iq)); |
228 | < | if (i > 4 && (fabs(w(ip))+g) == fabs(w(ip)) |
229 | < | && (fabs(w(iq))+g) == fabs(w(iq))) |
230 | < | { |
231 | < | a(ip, iq) = 0.0; |
232 | < | } |
233 | < | else if (fabs(a(ip, iq)) > tresh) |
234 | < | { |
235 | < | h = w(iq) - w(ip); |
233 | < | if ( (fabs(h)+g) == fabs(h)) t = (a(ip, iq)) / h; |
234 | < | else |
235 | < | { |
236 | < | theta = 0.5*h / (a(ip, iq)); |
237 | < | t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta)); |
238 | < | if (theta < 0.0) t = -t; |
239 | < | } |
240 | < | c = 1.0 / sqrt(1+t*t); |
241 | < | s = t*c; |
242 | < | tau = s/(1.0+c); |
243 | < | h = t*a(ip, iq); |
244 | < | z(ip) -= h; |
245 | < | z(iq) += h; |
246 | < | w(ip) -= h; |
247 | < | w(iq) += h; |
248 | < | a(ip, iq)=0.0; |
249 | < | for (j=0;j<ip-1;j++) |
250 | < | { |
251 | < | ROT(a,j,ip,j,iq); |
252 | < | } |
253 | < | for (j=ip+1;j<iq-1;j++) |
254 | < | { |
255 | < | ROT(a,ip,j,j,iq); |
256 | < | } |
257 | < | for (j=iq+1; j<N; j++) |
258 | < | { |
259 | < | ROT(a,ip,j,iq,j); |
260 | < | } |
261 | < | for (j=0; j<N; j++) |
262 | < | { |
263 | < | ROT(v,j,ip,j,iq); |
264 | < | } |
265 | < | } |
266 | < | } |
267 | < | } |
223 | > | for (ip=0; ip<2; ip++) { |
224 | > | for (iq=ip+1; iq<N; iq++) { |
225 | > | g = 100.0*fabs(a(ip, iq)); |
226 | > | if (i > 4 && (fabs(w(ip))+g) == fabs(w(ip)) |
227 | > | && (fabs(w(iq))+g) == fabs(w(iq))) { |
228 | > | a(ip, iq) = 0.0; |
229 | > | } else if (fabs(a(ip, iq)) > tresh) { |
230 | > | h = w(iq) - w(ip); |
231 | > | if ( (fabs(h)+g) == fabs(h)) { |
232 | > | t = (a(ip, iq)) / h; |
233 | > | } else { |
234 | > | theta = 0.5*h / (a(ip, iq)); |
235 | > | t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta)); |
236 | ||
237 | < | for (ip=0; ip<N; ip++) |
238 | < | { |
239 | < | b(ip) += z(ip); |
272 | < | w(ip) = b(ip); |
273 | < | z(ip) = 0.0; |
274 | < | } |
275 | < | } |
237 | > | if (theta < 0.0) |
238 | > | t = -t; |
239 | > | } |
240 | ||
241 | + | c = 1.0 / sqrt(1+t*t); |
242 | + | s = t*c; |
243 | + | tau = s/(1.0+c); |
244 | + | h = t*a(ip, iq); |
245 | + | z(ip) -= h; |
246 | + | z(iq) += h; |
247 | + | w(ip) -= h; |
248 | + | w(iq) += h; |
249 | + | a(ip, iq)=0.0; |
250 | + | |
251 | + | for (j=0;j<ip-1;j++) |
252 | + | ROT(a,j,ip,j,iq); |
253 | + | |
254 | + | for (j=ip+1;j<iq-1;j++) |
255 | + | ROT(a,ip,j,j,iq); |
256 | + | |
257 | + | for (j=iq+1; j<N; j++) |
258 | + | ROT(a,ip,j,iq,j); |
259 | + | for (j=0; j<N; j++) |
260 | + | ROT(v,j,ip,j,iq); |
261 | + | } |
262 | + | } |
263 | + | }//for (ip=0; ip<2; ip++) |
264 | + | |
265 | + | for (ip=0; ip<N; ip++) { |
266 | + | b(ip) += z(ip); |
267 | + | w(ip) = b(ip); |
268 | + | z(ip) = 0.0; |
269 | + | } |
270 | + | |
271 | + | } // end for (i=0; i<MAX_ROTATIONS; i++) |
272 | + | |
273 | if ( i >= MAX_ROTATIONS ) | |
274 | < | return false; |
274 | > | return false; |
275 | ||
276 | // sort eigenfunctions | |
277 | < | for (j=0; j<N; j++) |
278 | < | { |
279 | < | k = j; |
280 | < | tmp = w(k); |
281 | < | for (i=j; i<N; i++) |
282 | < | { |
283 | < | if (w(i) >= tmp) |
284 | < | { |
285 | < | k = i; |
286 | < | tmp = w(k); |
287 | < | } |
288 | < | } |
289 | < | if (k != j) |
290 | < | { |
291 | < | w(k) = w(j); |
292 | < | w(j) = tmp; |
293 | < | for (i=0; i<N; i++) |
294 | < | { |
295 | < | tmp = v(i, j); |
300 | < | v(i, j) = v(i, k); |
301 | < | v(i, k) = tmp; |
302 | < | } |
303 | < | } |
277 | > | for (j=0; j<N; j++) { |
278 | > | k = j; |
279 | > | tmp = w(k); |
280 | > | for (i=j; i<N; i++) { |
281 | > | if (w(i) >= tmp) { |
282 | > | k = i; |
283 | > | tmp = w(k); |
284 | > | } |
285 | > | } |
286 | > | |
287 | > | if (k != j) { |
288 | > | w(k) = w(j); |
289 | > | w(j) = tmp; |
290 | > | for (i=0; i<N; i++) { |
291 | > | tmp = v(i, j); |
292 | > | v(i, j) = v(i, k); |
293 | > | v(i, k) = tmp; |
294 | > | } |
295 | > | } |
296 | } | |
297 | ||
298 | // insure eigenvector consistency (i.e., Jacobi can compute | |
# | Line 309 | Line 301 | bool SquareMatrix<Real, Dim>::jacobi(const SquareMatri | |
301 | // hyperstreamline/other stuff. We will select the most | |
302 | // positive eigenvector. | |
303 | int numPos; | |
304 | < | for (j=0; j<N; j++) |
305 | < | { |
306 | < | for (numPos=0, i=0; i<N; i++) if ( v(i, j) >= 0.0 ) numPos++; |
315 | < | if ( numPos < 2 ) for(i=0; i<N; i++) v(i, j) *= -1.0; |
304 | > | for (j=0; j<N; j++) { |
305 | > | for (numPos=0, i=0; i<N; i++) if ( v(i, j) >= 0.0 ) numPos++; |
306 | > | if ( numPos < 2 ) for(i=0; i<N; i++) v(i, j) *= -1.0; |
307 | } | |
308 | ||
309 | return true; |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |