1 | /*************************************************************************
|
---|
2 | Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
|
---|
3 |
|
---|
4 | Contributors:
|
---|
5 | * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
|
---|
6 | pseudocode.
|
---|
7 |
|
---|
8 | See subroutines comments for additional copyrights.
|
---|
9 |
|
---|
10 | >>> SOURCE LICENSE >>>
|
---|
11 | This program is free software; you can redistribute it and/or modify
|
---|
12 | it under the terms of the GNU General Public License as published by
|
---|
13 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
14 | License, or (at your option) any later version.
|
---|
15 |
|
---|
16 | This program is distributed in the hope that it will be useful,
|
---|
17 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
18 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
19 | GNU General Public License for more details.
|
---|
20 |
|
---|
21 | A copy of the GNU General Public License is available at
|
---|
22 | http://www.fsf.org/licensing/licenses
|
---|
23 |
|
---|
24 | >>> END OF LICENSE >>>
|
---|
25 | *************************************************************************/
|
---|
26 |
|
---|
27 | using System;
|
---|
28 |
|
---|
29 | namespace alglib
|
---|
30 | {
|
---|
31 | public class hsschur
|
---|
32 | {
|
---|
33 | /*************************************************************************
|
---|
34 | Subroutine performing the Schur decomposition of a matrix in upper
|
---|
35 | Hessenberg form using the QR algorithm with multiple shifts.
|
---|
36 |
|
---|
37 | The source matrix H is represented as S'*H*S = T, where H - matrix in
|
---|
38 | upper Hessenberg form, S - orthogonal matrix (Schur vectors), T - upper
|
---|
39 | quasi-triangular matrix (with blocks of sizes 1x1 and 2x2 on the main
|
---|
40 | diagonal).
|
---|
41 |
|
---|
42 | Input parameters:
|
---|
43 | H - matrix to be decomposed.
|
---|
44 | Array whose indexes range within [1..N, 1..N].
|
---|
45 | N - size of H, N>=0.
|
---|
46 |
|
---|
47 |
|
---|
48 | Output parameters:
|
---|
49 | H contains the matrix T.
|
---|
50 | Array whose indexes range within [1..N, 1..N].
|
---|
51 | All elements below the blocks on the main diagonal are equal
|
---|
52 | to 0.
|
---|
53 | S - contains Schur vectors.
|
---|
54 | Array whose indexes range within [1..N, 1..N].
|
---|
55 |
|
---|
56 | Note 1:
|
---|
57 | The block structure of matrix T could be easily recognized: since all
|
---|
58 | the elements below the blocks are zeros, the elements a[i+1,i] which
|
---|
59 | are equal to 0 show the block border.
|
---|
60 |
|
---|
61 | Note 2:
|
---|
62 | the algorithm performance depends on the value of the internal
|
---|
63 | parameter NS of InternalSchurDecomposition subroutine which defines
|
---|
64 | the number of shifts in the QR algorithm (analog of the block width
|
---|
65 | in block matrix algorithms in linear algebra). If you require maximum
|
---|
66 | performance on your machine, it is recommended to adjust this
|
---|
67 | parameter manually.
|
---|
68 |
|
---|
69 | Result:
|
---|
70 | True, if the algorithm has converged and the parameters H and S contain
|
---|
71 | the result.
|
---|
72 | False, if the algorithm has not converged.
|
---|
73 |
|
---|
74 | Algorithm implemented on the basis of subroutine DHSEQR (LAPACK 3.0 library).
|
---|
75 | *************************************************************************/
|
---|
76 | public static bool upperhessenbergschurdecomposition(ref double[,] h,
|
---|
77 | int n,
|
---|
78 | ref double[,] s)
|
---|
79 | {
|
---|
80 | bool result = new bool();
|
---|
81 | double[] wi = new double[0];
|
---|
82 | double[] wr = new double[0];
|
---|
83 | int info = 0;
|
---|
84 |
|
---|
85 | internalschurdecomposition(ref h, n, 1, 2, ref wr, ref wi, ref s, ref info);
|
---|
86 | result = info==0;
|
---|
87 | return result;
|
---|
88 | }
|
---|
89 |
|
---|
90 |
|
---|
91 | public static void internalschurdecomposition(ref double[,] h,
|
---|
92 | int n,
|
---|
93 | int tneeded,
|
---|
94 | int zneeded,
|
---|
95 | ref double[] wr,
|
---|
96 | ref double[] wi,
|
---|
97 | ref double[,] z,
|
---|
98 | ref int info)
|
---|
99 | {
|
---|
100 | double[] work = new double[0];
|
---|
101 | int i = 0;
|
---|
102 | int i1 = 0;
|
---|
103 | int i2 = 0;
|
---|
104 | int ierr = 0;
|
---|
105 | int ii = 0;
|
---|
106 | int itemp = 0;
|
---|
107 | int itn = 0;
|
---|
108 | int its = 0;
|
---|
109 | int j = 0;
|
---|
110 | int k = 0;
|
---|
111 | int l = 0;
|
---|
112 | int maxb = 0;
|
---|
113 | int nr = 0;
|
---|
114 | int ns = 0;
|
---|
115 | int nv = 0;
|
---|
116 | double absw = 0;
|
---|
117 | double ovfl = 0;
|
---|
118 | double smlnum = 0;
|
---|
119 | double tau = 0;
|
---|
120 | double temp = 0;
|
---|
121 | double tst1 = 0;
|
---|
122 | double ulp = 0;
|
---|
123 | double unfl = 0;
|
---|
124 | double[,] s = new double[0,0];
|
---|
125 | double[] v = new double[0];
|
---|
126 | double[] vv = new double[0];
|
---|
127 | double[] workc1 = new double[0];
|
---|
128 | double[] works1 = new double[0];
|
---|
129 | double[] workv3 = new double[0];
|
---|
130 | double[] tmpwr = new double[0];
|
---|
131 | double[] tmpwi = new double[0];
|
---|
132 | bool initz = new bool();
|
---|
133 | bool wantt = new bool();
|
---|
134 | bool wantz = new bool();
|
---|
135 | double cnst = 0;
|
---|
136 | bool failflag = new bool();
|
---|
137 | int p1 = 0;
|
---|
138 | int p2 = 0;
|
---|
139 | double vt = 0;
|
---|
140 | int i_ = 0;
|
---|
141 | int i1_ = 0;
|
---|
142 |
|
---|
143 |
|
---|
144 | //
|
---|
145 | // Set the order of the multi-shift QR algorithm to be used.
|
---|
146 | // If you want to tune algorithm, change this values
|
---|
147 | //
|
---|
148 | ns = 12;
|
---|
149 | maxb = 50;
|
---|
150 |
|
---|
151 | //
|
---|
152 | // Now 2 < NS <= MAXB < NH.
|
---|
153 | //
|
---|
154 | maxb = Math.Max(3, maxb);
|
---|
155 | ns = Math.Min(maxb, ns);
|
---|
156 |
|
---|
157 | //
|
---|
158 | // Initialize
|
---|
159 | //
|
---|
160 | cnst = 1.5;
|
---|
161 | work = new double[Math.Max(n, 1)+1];
|
---|
162 | s = new double[ns+1, ns+1];
|
---|
163 | v = new double[ns+1+1];
|
---|
164 | vv = new double[ns+1+1];
|
---|
165 | wr = new double[Math.Max(n, 1)+1];
|
---|
166 | wi = new double[Math.Max(n, 1)+1];
|
---|
167 | workc1 = new double[1+1];
|
---|
168 | works1 = new double[1+1];
|
---|
169 | workv3 = new double[3+1];
|
---|
170 | tmpwr = new double[Math.Max(n, 1)+1];
|
---|
171 | tmpwi = new double[Math.Max(n, 1)+1];
|
---|
172 | System.Diagnostics.Debug.Assert(n>=0, "InternalSchurDecomposition: incorrect N!");
|
---|
173 | System.Diagnostics.Debug.Assert(tneeded==0 | tneeded==1, "InternalSchurDecomposition: incorrect TNeeded!");
|
---|
174 | System.Diagnostics.Debug.Assert(zneeded==0 | zneeded==1 | zneeded==2, "InternalSchurDecomposition: incorrect ZNeeded!");
|
---|
175 | wantt = tneeded==1;
|
---|
176 | initz = zneeded==2;
|
---|
177 | wantz = zneeded!=0;
|
---|
178 | info = 0;
|
---|
179 |
|
---|
180 | //
|
---|
181 | // Initialize Z, if necessary
|
---|
182 | //
|
---|
183 | if( initz )
|
---|
184 | {
|
---|
185 | z = new double[n+1, n+1];
|
---|
186 | for(i=1; i<=n; i++)
|
---|
187 | {
|
---|
188 | for(j=1; j<=n; j++)
|
---|
189 | {
|
---|
190 | if( i==j )
|
---|
191 | {
|
---|
192 | z[i,j] = 1;
|
---|
193 | }
|
---|
194 | else
|
---|
195 | {
|
---|
196 | z[i,j] = 0;
|
---|
197 | }
|
---|
198 | }
|
---|
199 | }
|
---|
200 | }
|
---|
201 |
|
---|
202 | //
|
---|
203 | // Quick return if possible
|
---|
204 | //
|
---|
205 | if( n==0 )
|
---|
206 | {
|
---|
207 | return;
|
---|
208 | }
|
---|
209 | if( n==1 )
|
---|
210 | {
|
---|
211 | wr[1] = h[1,1];
|
---|
212 | wi[1] = 0;
|
---|
213 | return;
|
---|
214 | }
|
---|
215 |
|
---|
216 | //
|
---|
217 | // Set rows and columns 1 to N to zero below the first
|
---|
218 | // subdiagonal.
|
---|
219 | //
|
---|
220 | for(j=1; j<=n-2; j++)
|
---|
221 | {
|
---|
222 | for(i=j+2; i<=n; i++)
|
---|
223 | {
|
---|
224 | h[i,j] = 0;
|
---|
225 | }
|
---|
226 | }
|
---|
227 |
|
---|
228 | //
|
---|
229 | // Test if N is sufficiently small
|
---|
230 | //
|
---|
231 | if( ns<=2 | ns>n | maxb>=n )
|
---|
232 | {
|
---|
233 |
|
---|
234 | //
|
---|
235 | // Use the standard double-shift algorithm
|
---|
236 | //
|
---|
237 | internalauxschur(wantt, wantz, n, 1, n, ref h, ref wr, ref wi, 1, n, ref z, ref work, ref workv3, ref workc1, ref works1, ref info);
|
---|
238 |
|
---|
239 | //
|
---|
240 | // fill entries under diagonal blocks of T with zeros
|
---|
241 | //
|
---|
242 | if( wantt )
|
---|
243 | {
|
---|
244 | j = 1;
|
---|
245 | while( j<=n )
|
---|
246 | {
|
---|
247 | if( (double)(wi[j])==(double)(0) )
|
---|
248 | {
|
---|
249 | for(i=j+1; i<=n; i++)
|
---|
250 | {
|
---|
251 | h[i,j] = 0;
|
---|
252 | }
|
---|
253 | j = j+1;
|
---|
254 | }
|
---|
255 | else
|
---|
256 | {
|
---|
257 | for(i=j+2; i<=n; i++)
|
---|
258 | {
|
---|
259 | h[i,j] = 0;
|
---|
260 | h[i,j+1] = 0;
|
---|
261 | }
|
---|
262 | j = j+2;
|
---|
263 | }
|
---|
264 | }
|
---|
265 | }
|
---|
266 | return;
|
---|
267 | }
|
---|
268 | unfl = AP.Math.MinRealNumber;
|
---|
269 | ovfl = 1/unfl;
|
---|
270 | ulp = 2*AP.Math.MachineEpsilon;
|
---|
271 | smlnum = unfl*(n/ulp);
|
---|
272 |
|
---|
273 | //
|
---|
274 | // I1 and I2 are the indices of the first row and last column of H
|
---|
275 | // to which transformations must be applied. If eigenvalues only are
|
---|
276 | // being computed, I1 and I2 are set inside the main loop.
|
---|
277 | //
|
---|
278 | i1 = 1;
|
---|
279 | i2 = n;
|
---|
280 |
|
---|
281 | //
|
---|
282 | // ITN is the total number of multiple-shift QR iterations allowed.
|
---|
283 | //
|
---|
284 | itn = 30*n;
|
---|
285 |
|
---|
286 | //
|
---|
287 | // The main loop begins here. I is the loop index and decreases from
|
---|
288 | // IHI to ILO in steps of at most MAXB. Each iteration of the loop
|
---|
289 | // works with the active submatrix in rows and columns L to I.
|
---|
290 | // Eigenvalues I+1 to IHI have already converged. Either L = ILO or
|
---|
291 | // H(L,L-1) is negligible so that the matrix splits.
|
---|
292 | //
|
---|
293 | i = n;
|
---|
294 | while( true )
|
---|
295 | {
|
---|
296 | l = 1;
|
---|
297 | if( i<1 )
|
---|
298 | {
|
---|
299 |
|
---|
300 | //
|
---|
301 | // fill entries under diagonal blocks of T with zeros
|
---|
302 | //
|
---|
303 | if( wantt )
|
---|
304 | {
|
---|
305 | j = 1;
|
---|
306 | while( j<=n )
|
---|
307 | {
|
---|
308 | if( (double)(wi[j])==(double)(0) )
|
---|
309 | {
|
---|
310 | for(i=j+1; i<=n; i++)
|
---|
311 | {
|
---|
312 | h[i,j] = 0;
|
---|
313 | }
|
---|
314 | j = j+1;
|
---|
315 | }
|
---|
316 | else
|
---|
317 | {
|
---|
318 | for(i=j+2; i<=n; i++)
|
---|
319 | {
|
---|
320 | h[i,j] = 0;
|
---|
321 | h[i,j+1] = 0;
|
---|
322 | }
|
---|
323 | j = j+2;
|
---|
324 | }
|
---|
325 | }
|
---|
326 | }
|
---|
327 |
|
---|
328 | //
|
---|
329 | // Exit
|
---|
330 | //
|
---|
331 | return;
|
---|
332 | }
|
---|
333 |
|
---|
334 | //
|
---|
335 | // Perform multiple-shift QR iterations on rows and columns ILO to I
|
---|
336 | // until a submatrix of order at most MAXB splits off at the bottom
|
---|
337 | // because a subdiagonal element has become negligible.
|
---|
338 | //
|
---|
339 | failflag = true;
|
---|
340 | for(its=0; its<=itn; its++)
|
---|
341 | {
|
---|
342 |
|
---|
343 | //
|
---|
344 | // Look for a single small subdiagonal element.
|
---|
345 | //
|
---|
346 | for(k=i; k>=l+1; k--)
|
---|
347 | {
|
---|
348 | tst1 = Math.Abs(h[k-1,k-1])+Math.Abs(h[k,k]);
|
---|
349 | if( (double)(tst1)==(double)(0) )
|
---|
350 | {
|
---|
351 | tst1 = blas.upperhessenberg1norm(ref h, l, i, l, i, ref work);
|
---|
352 | }
|
---|
353 | if( (double)(Math.Abs(h[k,k-1]))<=(double)(Math.Max(ulp*tst1, smlnum)) )
|
---|
354 | {
|
---|
355 | break;
|
---|
356 | }
|
---|
357 | }
|
---|
358 | l = k;
|
---|
359 | if( l>1 )
|
---|
360 | {
|
---|
361 |
|
---|
362 | //
|
---|
363 | // H(L,L-1) is negligible.
|
---|
364 | //
|
---|
365 | h[l,l-1] = 0;
|
---|
366 | }
|
---|
367 |
|
---|
368 | //
|
---|
369 | // Exit from loop if a submatrix of order <= MAXB has split off.
|
---|
370 | //
|
---|
371 | if( l>=i-maxb+1 )
|
---|
372 | {
|
---|
373 | failflag = false;
|
---|
374 | break;
|
---|
375 | }
|
---|
376 |
|
---|
377 | //
|
---|
378 | // Now the active submatrix is in rows and columns L to I. If
|
---|
379 | // eigenvalues only are being computed, only the active submatrix
|
---|
380 | // need be transformed.
|
---|
381 | //
|
---|
382 | if( its==20 | its==30 )
|
---|
383 | {
|
---|
384 |
|
---|
385 | //
|
---|
386 | // Exceptional shifts.
|
---|
387 | //
|
---|
388 | for(ii=i-ns+1; ii<=i; ii++)
|
---|
389 | {
|
---|
390 | wr[ii] = cnst*(Math.Abs(h[ii,ii-1])+Math.Abs(h[ii,ii]));
|
---|
391 | wi[ii] = 0;
|
---|
392 | }
|
---|
393 | }
|
---|
394 | else
|
---|
395 | {
|
---|
396 |
|
---|
397 | //
|
---|
398 | // Use eigenvalues of trailing submatrix of order NS as shifts.
|
---|
399 | //
|
---|
400 | blas.copymatrix(ref h, i-ns+1, i, i-ns+1, i, ref s, 1, ns, 1, ns);
|
---|
401 | internalauxschur(false, false, ns, 1, ns, ref s, ref tmpwr, ref tmpwi, 1, ns, ref z, ref work, ref workv3, ref workc1, ref works1, ref ierr);
|
---|
402 | for(p1=1; p1<=ns; p1++)
|
---|
403 | {
|
---|
404 | wr[i-ns+p1] = tmpwr[p1];
|
---|
405 | wi[i-ns+p1] = tmpwi[p1];
|
---|
406 | }
|
---|
407 | if( ierr>0 )
|
---|
408 | {
|
---|
409 |
|
---|
410 | //
|
---|
411 | // If DLAHQR failed to compute all NS eigenvalues, use the
|
---|
412 | // unconverged diagonal elements as the remaining shifts.
|
---|
413 | //
|
---|
414 | for(ii=1; ii<=ierr; ii++)
|
---|
415 | {
|
---|
416 | wr[i-ns+ii] = s[ii,ii];
|
---|
417 | wi[i-ns+ii] = 0;
|
---|
418 | }
|
---|
419 | }
|
---|
420 | }
|
---|
421 |
|
---|
422 | //
|
---|
423 | // Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns))
|
---|
424 | // where G is the Hessenberg submatrix H(L:I,L:I) and w is
|
---|
425 | // the vector of shifts (stored in WR and WI). The result is
|
---|
426 | // stored in the local array V.
|
---|
427 | //
|
---|
428 | v[1] = 1;
|
---|
429 | for(ii=2; ii<=ns+1; ii++)
|
---|
430 | {
|
---|
431 | v[ii] = 0;
|
---|
432 | }
|
---|
433 | nv = 1;
|
---|
434 | for(j=i-ns+1; j<=i; j++)
|
---|
435 | {
|
---|
436 | if( (double)(wi[j])>=(double)(0) )
|
---|
437 | {
|
---|
438 | if( (double)(wi[j])==(double)(0) )
|
---|
439 | {
|
---|
440 |
|
---|
441 | //
|
---|
442 | // real shift
|
---|
443 | //
|
---|
444 | p1 = nv+1;
|
---|
445 | for(i_=1; i_<=p1;i_++)
|
---|
446 | {
|
---|
447 | vv[i_] = v[i_];
|
---|
448 | }
|
---|
449 | blas.matrixvectormultiply(ref h, l, l+nv, l, l+nv-1, false, ref vv, 1, nv, 1.0, ref v, 1, nv+1, -wr[j]);
|
---|
450 | nv = nv+1;
|
---|
451 | }
|
---|
452 | else
|
---|
453 | {
|
---|
454 | if( (double)(wi[j])>(double)(0) )
|
---|
455 | {
|
---|
456 |
|
---|
457 | //
|
---|
458 | // complex conjugate pair of shifts
|
---|
459 | //
|
---|
460 | p1 = nv+1;
|
---|
461 | for(i_=1; i_<=p1;i_++)
|
---|
462 | {
|
---|
463 | vv[i_] = v[i_];
|
---|
464 | }
|
---|
465 | blas.matrixvectormultiply(ref h, l, l+nv, l, l+nv-1, false, ref v, 1, nv, 1.0, ref vv, 1, nv+1, -(2*wr[j]));
|
---|
466 | itemp = blas.vectoridxabsmax(ref vv, 1, nv+1);
|
---|
467 | temp = 1/Math.Max(Math.Abs(vv[itemp]), smlnum);
|
---|
468 | p1 = nv+1;
|
---|
469 | for(i_=1; i_<=p1;i_++)
|
---|
470 | {
|
---|
471 | vv[i_] = temp*vv[i_];
|
---|
472 | }
|
---|
473 | absw = blas.pythag2(wr[j], wi[j]);
|
---|
474 | temp = temp*absw*absw;
|
---|
475 | blas.matrixvectormultiply(ref h, l, l+nv+1, l, l+nv, false, ref vv, 1, nv+1, 1.0, ref v, 1, nv+2, temp);
|
---|
476 | nv = nv+2;
|
---|
477 | }
|
---|
478 | }
|
---|
479 |
|
---|
480 | //
|
---|
481 | // Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero,
|
---|
482 | // reset it to the unit vector.
|
---|
483 | //
|
---|
484 | itemp = blas.vectoridxabsmax(ref v, 1, nv);
|
---|
485 | temp = Math.Abs(v[itemp]);
|
---|
486 | if( (double)(temp)==(double)(0) )
|
---|
487 | {
|
---|
488 | v[1] = 1;
|
---|
489 | for(ii=2; ii<=nv; ii++)
|
---|
490 | {
|
---|
491 | v[ii] = 0;
|
---|
492 | }
|
---|
493 | }
|
---|
494 | else
|
---|
495 | {
|
---|
496 | temp = Math.Max(temp, smlnum);
|
---|
497 | vt = 1/temp;
|
---|
498 | for(i_=1; i_<=nv;i_++)
|
---|
499 | {
|
---|
500 | v[i_] = vt*v[i_];
|
---|
501 | }
|
---|
502 | }
|
---|
503 | }
|
---|
504 | }
|
---|
505 |
|
---|
506 | //
|
---|
507 | // Multiple-shift QR step
|
---|
508 | //
|
---|
509 | for(k=l; k<=i-1; k++)
|
---|
510 | {
|
---|
511 |
|
---|
512 | //
|
---|
513 | // The first iteration of this loop determines a reflection G
|
---|
514 | // from the vector V and applies it from left and right to H,
|
---|
515 | // thus creating a nonzero bulge below the subdiagonal.
|
---|
516 | //
|
---|
517 | // Each subsequent iteration determines a reflection G to
|
---|
518 | // restore the Hessenberg form in the (K-1)th column, and thus
|
---|
519 | // chases the bulge one step toward the bottom of the active
|
---|
520 | // submatrix. NR is the order of G.
|
---|
521 | //
|
---|
522 | nr = Math.Min(ns+1, i-k+1);
|
---|
523 | if( k>l )
|
---|
524 | {
|
---|
525 | p1 = k-1;
|
---|
526 | p2 = k+nr-1;
|
---|
527 | i1_ = (k) - (1);
|
---|
528 | for(i_=1; i_<=nr;i_++)
|
---|
529 | {
|
---|
530 | v[i_] = h[i_+i1_,p1];
|
---|
531 | }
|
---|
532 | }
|
---|
533 | reflections.generatereflection(ref v, nr, ref tau);
|
---|
534 | if( k>l )
|
---|
535 | {
|
---|
536 | h[k,k-1] = v[1];
|
---|
537 | for(ii=k+1; ii<=i; ii++)
|
---|
538 | {
|
---|
539 | h[ii,k-1] = 0;
|
---|
540 | }
|
---|
541 | }
|
---|
542 | v[1] = 1;
|
---|
543 |
|
---|
544 | //
|
---|
545 | // Apply G from the left to transform the rows of the matrix in
|
---|
546 | // columns K to I2.
|
---|
547 | //
|
---|
548 | reflections.applyreflectionfromtheleft(ref h, tau, ref v, k, k+nr-1, k, i2, ref work);
|
---|
549 |
|
---|
550 | //
|
---|
551 | // Apply G from the right to transform the columns of the
|
---|
552 | // matrix in rows I1 to min(K+NR,I).
|
---|
553 | //
|
---|
554 | reflections.applyreflectionfromtheright(ref h, tau, ref v, i1, Math.Min(k+nr, i), k, k+nr-1, ref work);
|
---|
555 | if( wantz )
|
---|
556 | {
|
---|
557 |
|
---|
558 | //
|
---|
559 | // Accumulate transformations in the matrix Z
|
---|
560 | //
|
---|
561 | reflections.applyreflectionfromtheright(ref z, tau, ref v, 1, n, k, k+nr-1, ref work);
|
---|
562 | }
|
---|
563 | }
|
---|
564 | }
|
---|
565 |
|
---|
566 | //
|
---|
567 | // Failure to converge in remaining number of iterations
|
---|
568 | //
|
---|
569 | if( failflag )
|
---|
570 | {
|
---|
571 | info = i;
|
---|
572 | return;
|
---|
573 | }
|
---|
574 |
|
---|
575 | //
|
---|
576 | // A submatrix of order <= MAXB in rows and columns L to I has split
|
---|
577 | // off. Use the double-shift QR algorithm to handle it.
|
---|
578 | //
|
---|
579 | internalauxschur(wantt, wantz, n, l, i, ref h, ref wr, ref wi, 1, n, ref z, ref work, ref workv3, ref workc1, ref works1, ref info);
|
---|
580 | if( info>0 )
|
---|
581 | {
|
---|
582 | return;
|
---|
583 | }
|
---|
584 |
|
---|
585 | //
|
---|
586 | // Decrement number of remaining iterations, and return to start of
|
---|
587 | // the main loop with a new value of I.
|
---|
588 | //
|
---|
589 | itn = itn-its;
|
---|
590 | i = l-1;
|
---|
591 | }
|
---|
592 | }
|
---|
593 |
|
---|
594 |
|
---|
595 | private static void internalauxschur(bool wantt,
|
---|
596 | bool wantz,
|
---|
597 | int n,
|
---|
598 | int ilo,
|
---|
599 | int ihi,
|
---|
600 | ref double[,] h,
|
---|
601 | ref double[] wr,
|
---|
602 | ref double[] wi,
|
---|
603 | int iloz,
|
---|
604 | int ihiz,
|
---|
605 | ref double[,] z,
|
---|
606 | ref double[] work,
|
---|
607 | ref double[] workv3,
|
---|
608 | ref double[] workc1,
|
---|
609 | ref double[] works1,
|
---|
610 | ref int info)
|
---|
611 | {
|
---|
612 | int i = 0;
|
---|
613 | int i1 = 0;
|
---|
614 | int i2 = 0;
|
---|
615 | int itn = 0;
|
---|
616 | int its = 0;
|
---|
617 | int j = 0;
|
---|
618 | int k = 0;
|
---|
619 | int l = 0;
|
---|
620 | int m = 0;
|
---|
621 | int nh = 0;
|
---|
622 | int nr = 0;
|
---|
623 | int nz = 0;
|
---|
624 | double ave = 0;
|
---|
625 | double cs = 0;
|
---|
626 | double disc = 0;
|
---|
627 | double h00 = 0;
|
---|
628 | double h10 = 0;
|
---|
629 | double h11 = 0;
|
---|
630 | double h12 = 0;
|
---|
631 | double h21 = 0;
|
---|
632 | double h22 = 0;
|
---|
633 | double h33 = 0;
|
---|
634 | double h33s = 0;
|
---|
635 | double h43h34 = 0;
|
---|
636 | double h44 = 0;
|
---|
637 | double h44s = 0;
|
---|
638 | double ovfl = 0;
|
---|
639 | double s = 0;
|
---|
640 | double smlnum = 0;
|
---|
641 | double sn = 0;
|
---|
642 | double sum = 0;
|
---|
643 | double t1 = 0;
|
---|
644 | double t2 = 0;
|
---|
645 | double t3 = 0;
|
---|
646 | double tst1 = 0;
|
---|
647 | double unfl = 0;
|
---|
648 | double v1 = 0;
|
---|
649 | double v2 = 0;
|
---|
650 | double v3 = 0;
|
---|
651 | bool failflag = new bool();
|
---|
652 | double dat1 = 0;
|
---|
653 | double dat2 = 0;
|
---|
654 | int p1 = 0;
|
---|
655 | double him1im1 = 0;
|
---|
656 | double him1i = 0;
|
---|
657 | double hiim1 = 0;
|
---|
658 | double hii = 0;
|
---|
659 | double wrim1 = 0;
|
---|
660 | double wri = 0;
|
---|
661 | double wiim1 = 0;
|
---|
662 | double wii = 0;
|
---|
663 | double ulp = 0;
|
---|
664 |
|
---|
665 | info = 0;
|
---|
666 | dat1 = 0.75;
|
---|
667 | dat2 = -0.4375;
|
---|
668 | ulp = AP.Math.MachineEpsilon;
|
---|
669 |
|
---|
670 | //
|
---|
671 | // Quick return if possible
|
---|
672 | //
|
---|
673 | if( n==0 )
|
---|
674 | {
|
---|
675 | return;
|
---|
676 | }
|
---|
677 | if( ilo==ihi )
|
---|
678 | {
|
---|
679 | wr[ilo] = h[ilo,ilo];
|
---|
680 | wi[ilo] = 0;
|
---|
681 | return;
|
---|
682 | }
|
---|
683 | nh = ihi-ilo+1;
|
---|
684 | nz = ihiz-iloz+1;
|
---|
685 |
|
---|
686 | //
|
---|
687 | // Set machine-dependent constants for the stopping criterion.
|
---|
688 | // If norm(H) <= sqrt(OVFL), overflow should not occur.
|
---|
689 | //
|
---|
690 | unfl = AP.Math.MinRealNumber;
|
---|
691 | ovfl = 1/unfl;
|
---|
692 | smlnum = unfl*(nh/ulp);
|
---|
693 |
|
---|
694 | //
|
---|
695 | // I1 and I2 are the indices of the first row and last column of H
|
---|
696 | // to which transformations must be applied. If eigenvalues only are
|
---|
697 | // being computed, I1 and I2 are set inside the main loop.
|
---|
698 | //
|
---|
699 | i1 = 1;
|
---|
700 | i2 = n;
|
---|
701 |
|
---|
702 | //
|
---|
703 | // ITN is the total number of QR iterations allowed.
|
---|
704 | //
|
---|
705 | itn = 30*nh;
|
---|
706 |
|
---|
707 | //
|
---|
708 | // The main loop begins here. I is the loop index and decreases from
|
---|
709 | // IHI to ILO in steps of 1 or 2. Each iteration of the loop works
|
---|
710 | // with the active submatrix in rows and columns L to I.
|
---|
711 | // Eigenvalues I+1 to IHI have already converged. Either L = ILO or
|
---|
712 | // H(L,L-1) is negligible so that the matrix splits.
|
---|
713 | //
|
---|
714 | i = ihi;
|
---|
715 | while( true )
|
---|
716 | {
|
---|
717 | l = ilo;
|
---|
718 | if( i<ilo )
|
---|
719 | {
|
---|
720 | return;
|
---|
721 | }
|
---|
722 |
|
---|
723 | //
|
---|
724 | // Perform QR iterations on rows and columns ILO to I until a
|
---|
725 | // submatrix of order 1 or 2 splits off at the bottom because a
|
---|
726 | // subdiagonal element has become negligible.
|
---|
727 | //
|
---|
728 | failflag = true;
|
---|
729 | for(its=0; its<=itn; its++)
|
---|
730 | {
|
---|
731 |
|
---|
732 | //
|
---|
733 | // Look for a single small subdiagonal element.
|
---|
734 | //
|
---|
735 | for(k=i; k>=l+1; k--)
|
---|
736 | {
|
---|
737 | tst1 = Math.Abs(h[k-1,k-1])+Math.Abs(h[k,k]);
|
---|
738 | if( (double)(tst1)==(double)(0) )
|
---|
739 | {
|
---|
740 | tst1 = blas.upperhessenberg1norm(ref h, l, i, l, i, ref work);
|
---|
741 | }
|
---|
742 | if( (double)(Math.Abs(h[k,k-1]))<=(double)(Math.Max(ulp*tst1, smlnum)) )
|
---|
743 | {
|
---|
744 | break;
|
---|
745 | }
|
---|
746 | }
|
---|
747 | l = k;
|
---|
748 | if( l>ilo )
|
---|
749 | {
|
---|
750 |
|
---|
751 | //
|
---|
752 | // H(L,L-1) is negligible
|
---|
753 | //
|
---|
754 | h[l,l-1] = 0;
|
---|
755 | }
|
---|
756 |
|
---|
757 | //
|
---|
758 | // Exit from loop if a submatrix of order 1 or 2 has split off.
|
---|
759 | //
|
---|
760 | if( l>=i-1 )
|
---|
761 | {
|
---|
762 | failflag = false;
|
---|
763 | break;
|
---|
764 | }
|
---|
765 |
|
---|
766 | //
|
---|
767 | // Now the active submatrix is in rows and columns L to I. If
|
---|
768 | // eigenvalues only are being computed, only the active submatrix
|
---|
769 | // need be transformed.
|
---|
770 | //
|
---|
771 | if( its==10 | its==20 )
|
---|
772 | {
|
---|
773 |
|
---|
774 | //
|
---|
775 | // Exceptional shift.
|
---|
776 | //
|
---|
777 | s = Math.Abs(h[i,i-1])+Math.Abs(h[i-1,i-2]);
|
---|
778 | h44 = dat1*s+h[i,i];
|
---|
779 | h33 = h44;
|
---|
780 | h43h34 = dat2*s*s;
|
---|
781 | }
|
---|
782 | else
|
---|
783 | {
|
---|
784 |
|
---|
785 | //
|
---|
786 | // Prepare to use Francis' double shift
|
---|
787 | // (i.e. 2nd degree generalized Rayleigh quotient)
|
---|
788 | //
|
---|
789 | h44 = h[i,i];
|
---|
790 | h33 = h[i-1,i-1];
|
---|
791 | h43h34 = h[i,i-1]*h[i-1,i];
|
---|
792 | s = h[i-1,i-2]*h[i-1,i-2];
|
---|
793 | disc = (h33-h44)*0.5;
|
---|
794 | disc = disc*disc+h43h34;
|
---|
795 | if( (double)(disc)>(double)(0) )
|
---|
796 | {
|
---|
797 |
|
---|
798 | //
|
---|
799 | // Real roots: use Wilkinson's shift twice
|
---|
800 | //
|
---|
801 | disc = Math.Sqrt(disc);
|
---|
802 | ave = 0.5*(h33+h44);
|
---|
803 | if( (double)(Math.Abs(h33)-Math.Abs(h44))>(double)(0) )
|
---|
804 | {
|
---|
805 | h33 = h33*h44-h43h34;
|
---|
806 | h44 = h33/(extschursign(disc, ave)+ave);
|
---|
807 | }
|
---|
808 | else
|
---|
809 | {
|
---|
810 | h44 = extschursign(disc, ave)+ave;
|
---|
811 | }
|
---|
812 | h33 = h44;
|
---|
813 | h43h34 = 0;
|
---|
814 | }
|
---|
815 | }
|
---|
816 |
|
---|
817 | //
|
---|
818 | // Look for two consecutive small subdiagonal elements.
|
---|
819 | //
|
---|
820 | for(m=i-2; m>=l; m--)
|
---|
821 | {
|
---|
822 |
|
---|
823 | //
|
---|
824 | // Determine the effect of starting the double-shift QR
|
---|
825 | // iteration at row M, and see if this would make H(M,M-1)
|
---|
826 | // negligible.
|
---|
827 | //
|
---|
828 | h11 = h[m,m];
|
---|
829 | h22 = h[m+1,m+1];
|
---|
830 | h21 = h[m+1,m];
|
---|
831 | h12 = h[m,m+1];
|
---|
832 | h44s = h44-h11;
|
---|
833 | h33s = h33-h11;
|
---|
834 | v1 = (h33s*h44s-h43h34)/h21+h12;
|
---|
835 | v2 = h22-h11-h33s-h44s;
|
---|
836 | v3 = h[m+2,m+1];
|
---|
837 | s = Math.Abs(v1)+Math.Abs(v2)+Math.Abs(v3);
|
---|
838 | v1 = v1/s;
|
---|
839 | v2 = v2/s;
|
---|
840 | v3 = v3/s;
|
---|
841 | workv3[1] = v1;
|
---|
842 | workv3[2] = v2;
|
---|
843 | workv3[3] = v3;
|
---|
844 | if( m==l )
|
---|
845 | {
|
---|
846 | break;
|
---|
847 | }
|
---|
848 | h00 = h[m-1,m-1];
|
---|
849 | h10 = h[m,m-1];
|
---|
850 | tst1 = Math.Abs(v1)*(Math.Abs(h00)+Math.Abs(h11)+Math.Abs(h22));
|
---|
851 | if( (double)(Math.Abs(h10)*(Math.Abs(v2)+Math.Abs(v3)))<=(double)(ulp*tst1) )
|
---|
852 | {
|
---|
853 | break;
|
---|
854 | }
|
---|
855 | }
|
---|
856 |
|
---|
857 | //
|
---|
858 | // Double-shift QR step
|
---|
859 | //
|
---|
860 | for(k=m; k<=i-1; k++)
|
---|
861 | {
|
---|
862 |
|
---|
863 | //
|
---|
864 | // The first iteration of this loop determines a reflection G
|
---|
865 | // from the vector V and applies it from left and right to H,
|
---|
866 | // thus creating a nonzero bulge below the subdiagonal.
|
---|
867 | //
|
---|
868 | // Each subsequent iteration determines a reflection G to
|
---|
869 | // restore the Hessenberg form in the (K-1)th column, and thus
|
---|
870 | // chases the bulge one step toward the bottom of the active
|
---|
871 | // submatrix. NR is the order of G.
|
---|
872 | //
|
---|
873 | nr = Math.Min(3, i-k+1);
|
---|
874 | if( k>m )
|
---|
875 | {
|
---|
876 | for(p1=1; p1<=nr; p1++)
|
---|
877 | {
|
---|
878 | workv3[p1] = h[k+p1-1,k-1];
|
---|
879 | }
|
---|
880 | }
|
---|
881 | reflections.generatereflection(ref workv3, nr, ref t1);
|
---|
882 | if( k>m )
|
---|
883 | {
|
---|
884 | h[k,k-1] = workv3[1];
|
---|
885 | h[k+1,k-1] = 0;
|
---|
886 | if( k<i-1 )
|
---|
887 | {
|
---|
888 | h[k+2,k-1] = 0;
|
---|
889 | }
|
---|
890 | }
|
---|
891 | else
|
---|
892 | {
|
---|
893 | if( m>l )
|
---|
894 | {
|
---|
895 | h[k,k-1] = -h[k,k-1];
|
---|
896 | }
|
---|
897 | }
|
---|
898 | v2 = workv3[2];
|
---|
899 | t2 = t1*v2;
|
---|
900 | if( nr==3 )
|
---|
901 | {
|
---|
902 | v3 = workv3[3];
|
---|
903 | t3 = t1*v3;
|
---|
904 |
|
---|
905 | //
|
---|
906 | // Apply G from the left to transform the rows of the matrix
|
---|
907 | // in columns K to I2.
|
---|
908 | //
|
---|
909 | for(j=k; j<=i2; j++)
|
---|
910 | {
|
---|
911 | sum = h[k,j]+v2*h[k+1,j]+v3*h[k+2,j];
|
---|
912 | h[k,j] = h[k,j]-sum*t1;
|
---|
913 | h[k+1,j] = h[k+1,j]-sum*t2;
|
---|
914 | h[k+2,j] = h[k+2,j]-sum*t3;
|
---|
915 | }
|
---|
916 |
|
---|
917 | //
|
---|
918 | // Apply G from the right to transform the columns of the
|
---|
919 | // matrix in rows I1 to min(K+3,I).
|
---|
920 | //
|
---|
921 | for(j=i1; j<=Math.Min(k+3, i); j++)
|
---|
922 | {
|
---|
923 | sum = h[j,k]+v2*h[j,k+1]+v3*h[j,k+2];
|
---|
924 | h[j,k] = h[j,k]-sum*t1;
|
---|
925 | h[j,k+1] = h[j,k+1]-sum*t2;
|
---|
926 | h[j,k+2] = h[j,k+2]-sum*t3;
|
---|
927 | }
|
---|
928 | if( wantz )
|
---|
929 | {
|
---|
930 |
|
---|
931 | //
|
---|
932 | // Accumulate transformations in the matrix Z
|
---|
933 | //
|
---|
934 | for(j=iloz; j<=ihiz; j++)
|
---|
935 | {
|
---|
936 | sum = z[j,k]+v2*z[j,k+1]+v3*z[j,k+2];
|
---|
937 | z[j,k] = z[j,k]-sum*t1;
|
---|
938 | z[j,k+1] = z[j,k+1]-sum*t2;
|
---|
939 | z[j,k+2] = z[j,k+2]-sum*t3;
|
---|
940 | }
|
---|
941 | }
|
---|
942 | }
|
---|
943 | else
|
---|
944 | {
|
---|
945 | if( nr==2 )
|
---|
946 | {
|
---|
947 |
|
---|
948 | //
|
---|
949 | // Apply G from the left to transform the rows of the matrix
|
---|
950 | // in columns K to I2.
|
---|
951 | //
|
---|
952 | for(j=k; j<=i2; j++)
|
---|
953 | {
|
---|
954 | sum = h[k,j]+v2*h[k+1,j];
|
---|
955 | h[k,j] = h[k,j]-sum*t1;
|
---|
956 | h[k+1,j] = h[k+1,j]-sum*t2;
|
---|
957 | }
|
---|
958 |
|
---|
959 | //
|
---|
960 | // Apply G from the right to transform the columns of the
|
---|
961 | // matrix in rows I1 to min(K+3,I).
|
---|
962 | //
|
---|
963 | for(j=i1; j<=i; j++)
|
---|
964 | {
|
---|
965 | sum = h[j,k]+v2*h[j,k+1];
|
---|
966 | h[j,k] = h[j,k]-sum*t1;
|
---|
967 | h[j,k+1] = h[j,k+1]-sum*t2;
|
---|
968 | }
|
---|
969 | if( wantz )
|
---|
970 | {
|
---|
971 |
|
---|
972 | //
|
---|
973 | // Accumulate transformations in the matrix Z
|
---|
974 | //
|
---|
975 | for(j=iloz; j<=ihiz; j++)
|
---|
976 | {
|
---|
977 | sum = z[j,k]+v2*z[j,k+1];
|
---|
978 | z[j,k] = z[j,k]-sum*t1;
|
---|
979 | z[j,k+1] = z[j,k+1]-sum*t2;
|
---|
980 | }
|
---|
981 | }
|
---|
982 | }
|
---|
983 | }
|
---|
984 | }
|
---|
985 | }
|
---|
986 | if( failflag )
|
---|
987 | {
|
---|
988 |
|
---|
989 | //
|
---|
990 | // Failure to converge in remaining number of iterations
|
---|
991 | //
|
---|
992 | info = i;
|
---|
993 | return;
|
---|
994 | }
|
---|
995 | if( l==i )
|
---|
996 | {
|
---|
997 |
|
---|
998 | //
|
---|
999 | // H(I,I-1) is negligible: one eigenvalue has converged.
|
---|
1000 | //
|
---|
1001 | wr[i] = h[i,i];
|
---|
1002 | wi[i] = 0;
|
---|
1003 | }
|
---|
1004 | else
|
---|
1005 | {
|
---|
1006 | if( l==i-1 )
|
---|
1007 | {
|
---|
1008 |
|
---|
1009 | //
|
---|
1010 | // H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
|
---|
1011 | //
|
---|
1012 | // Transform the 2-by-2 submatrix to standard Schur form,
|
---|
1013 | // and compute and store the eigenvalues.
|
---|
1014 | //
|
---|
1015 | him1im1 = h[i-1,i-1];
|
---|
1016 | him1i = h[i-1,i];
|
---|
1017 | hiim1 = h[i,i-1];
|
---|
1018 | hii = h[i,i];
|
---|
1019 | aux2x2schur(ref him1im1, ref him1i, ref hiim1, ref hii, ref wrim1, ref wiim1, ref wri, ref wii, ref cs, ref sn);
|
---|
1020 | wr[i-1] = wrim1;
|
---|
1021 | wi[i-1] = wiim1;
|
---|
1022 | wr[i] = wri;
|
---|
1023 | wi[i] = wii;
|
---|
1024 | h[i-1,i-1] = him1im1;
|
---|
1025 | h[i-1,i] = him1i;
|
---|
1026 | h[i,i-1] = hiim1;
|
---|
1027 | h[i,i] = hii;
|
---|
1028 | if( wantt )
|
---|
1029 | {
|
---|
1030 |
|
---|
1031 | //
|
---|
1032 | // Apply the transformation to the rest of H.
|
---|
1033 | //
|
---|
1034 | if( i2>i )
|
---|
1035 | {
|
---|
1036 | workc1[1] = cs;
|
---|
1037 | works1[1] = sn;
|
---|
1038 | rotations.applyrotationsfromtheleft(true, i-1, i, i+1, i2, ref workc1, ref works1, ref h, ref work);
|
---|
1039 | }
|
---|
1040 | workc1[1] = cs;
|
---|
1041 | works1[1] = sn;
|
---|
1042 | rotations.applyrotationsfromtheright(true, i1, i-2, i-1, i, ref workc1, ref works1, ref h, ref work);
|
---|
1043 | }
|
---|
1044 | if( wantz )
|
---|
1045 | {
|
---|
1046 |
|
---|
1047 | //
|
---|
1048 | // Apply the transformation to Z.
|
---|
1049 | //
|
---|
1050 | workc1[1] = cs;
|
---|
1051 | works1[1] = sn;
|
---|
1052 | rotations.applyrotationsfromtheright(true, iloz, iloz+nz-1, i-1, i, ref workc1, ref works1, ref z, ref work);
|
---|
1053 | }
|
---|
1054 | }
|
---|
1055 | }
|
---|
1056 |
|
---|
1057 | //
|
---|
1058 | // Decrement number of remaining iterations, and return to start of
|
---|
1059 | // the main loop with new value of I.
|
---|
1060 | //
|
---|
1061 | itn = itn-its;
|
---|
1062 | i = l-1;
|
---|
1063 | }
|
---|
1064 | }
|
---|
1065 |
|
---|
1066 |
|
---|
1067 | private static void aux2x2schur(ref double a,
|
---|
1068 | ref double b,
|
---|
1069 | ref double c,
|
---|
1070 | ref double d,
|
---|
1071 | ref double rt1r,
|
---|
1072 | ref double rt1i,
|
---|
1073 | ref double rt2r,
|
---|
1074 | ref double rt2i,
|
---|
1075 | ref double cs,
|
---|
1076 | ref double sn)
|
---|
1077 | {
|
---|
1078 | double multpl = 0;
|
---|
1079 | double aa = 0;
|
---|
1080 | double bb = 0;
|
---|
1081 | double bcmax = 0;
|
---|
1082 | double bcmis = 0;
|
---|
1083 | double cc = 0;
|
---|
1084 | double cs1 = 0;
|
---|
1085 | double dd = 0;
|
---|
1086 | double eps = 0;
|
---|
1087 | double p = 0;
|
---|
1088 | double sab = 0;
|
---|
1089 | double sac = 0;
|
---|
1090 | double scl = 0;
|
---|
1091 | double sigma = 0;
|
---|
1092 | double sn1 = 0;
|
---|
1093 | double tau = 0;
|
---|
1094 | double temp = 0;
|
---|
1095 | double z = 0;
|
---|
1096 |
|
---|
1097 | multpl = 4.0;
|
---|
1098 | eps = AP.Math.MachineEpsilon;
|
---|
1099 | if( (double)(c)==(double)(0) )
|
---|
1100 | {
|
---|
1101 | cs = 1;
|
---|
1102 | sn = 0;
|
---|
1103 | }
|
---|
1104 | else
|
---|
1105 | {
|
---|
1106 | if( (double)(b)==(double)(0) )
|
---|
1107 | {
|
---|
1108 |
|
---|
1109 | //
|
---|
1110 | // Swap rows and columns
|
---|
1111 | //
|
---|
1112 | cs = 0;
|
---|
1113 | sn = 1;
|
---|
1114 | temp = d;
|
---|
1115 | d = a;
|
---|
1116 | a = temp;
|
---|
1117 | b = -c;
|
---|
1118 | c = 0;
|
---|
1119 | }
|
---|
1120 | else
|
---|
1121 | {
|
---|
1122 | if( (double)(a-d)==(double)(0) & extschursigntoone(b)!=extschursigntoone(c) )
|
---|
1123 | {
|
---|
1124 | cs = 1;
|
---|
1125 | sn = 0;
|
---|
1126 | }
|
---|
1127 | else
|
---|
1128 | {
|
---|
1129 | temp = a-d;
|
---|
1130 | p = 0.5*temp;
|
---|
1131 | bcmax = Math.Max(Math.Abs(b), Math.Abs(c));
|
---|
1132 | bcmis = Math.Min(Math.Abs(b), Math.Abs(c))*extschursigntoone(b)*extschursigntoone(c);
|
---|
1133 | scl = Math.Max(Math.Abs(p), bcmax);
|
---|
1134 | z = p/scl*p+bcmax/scl*bcmis;
|
---|
1135 |
|
---|
1136 | //
|
---|
1137 | // If Z is of the order of the machine accuracy, postpone the
|
---|
1138 | // decision on the nature of eigenvalues
|
---|
1139 | //
|
---|
1140 | if( (double)(z)>=(double)(multpl*eps) )
|
---|
1141 | {
|
---|
1142 |
|
---|
1143 | //
|
---|
1144 | // Real eigenvalues. Compute A and D.
|
---|
1145 | //
|
---|
1146 | z = p+extschursign(Math.Sqrt(scl)*Math.Sqrt(z), p);
|
---|
1147 | a = d+z;
|
---|
1148 | d = d-bcmax/z*bcmis;
|
---|
1149 |
|
---|
1150 | //
|
---|
1151 | // Compute B and the rotation matrix
|
---|
1152 | //
|
---|
1153 | tau = blas.pythag2(c, z);
|
---|
1154 | cs = z/tau;
|
---|
1155 | sn = c/tau;
|
---|
1156 | b = b-c;
|
---|
1157 | c = 0;
|
---|
1158 | }
|
---|
1159 | else
|
---|
1160 | {
|
---|
1161 |
|
---|
1162 | //
|
---|
1163 | // Complex eigenvalues, or real (almost) equal eigenvalues.
|
---|
1164 | // Make diagonal elements equal.
|
---|
1165 | //
|
---|
1166 | sigma = b+c;
|
---|
1167 | tau = blas.pythag2(sigma, temp);
|
---|
1168 | cs = Math.Sqrt(0.5*(1+Math.Abs(sigma)/tau));
|
---|
1169 | sn = -(p/(tau*cs)*extschursign(1, sigma));
|
---|
1170 |
|
---|
1171 | //
|
---|
1172 | // Compute [ AA BB ] = [ A B ] [ CS -SN ]
|
---|
1173 | // [ CC DD ] [ C D ] [ SN CS ]
|
---|
1174 | //
|
---|
1175 | aa = a*cs+b*sn;
|
---|
1176 | bb = -(a*sn)+b*cs;
|
---|
1177 | cc = c*cs+d*sn;
|
---|
1178 | dd = -(c*sn)+d*cs;
|
---|
1179 |
|
---|
1180 | //
|
---|
1181 | // Compute [ A B ] = [ CS SN ] [ AA BB ]
|
---|
1182 | // [ C D ] [-SN CS ] [ CC DD ]
|
---|
1183 | //
|
---|
1184 | a = aa*cs+cc*sn;
|
---|
1185 | b = bb*cs+dd*sn;
|
---|
1186 | c = -(aa*sn)+cc*cs;
|
---|
1187 | d = -(bb*sn)+dd*cs;
|
---|
1188 | temp = 0.5*(a+d);
|
---|
1189 | a = temp;
|
---|
1190 | d = temp;
|
---|
1191 | if( (double)(c)!=(double)(0) )
|
---|
1192 | {
|
---|
1193 | if( (double)(b)!=(double)(0) )
|
---|
1194 | {
|
---|
1195 | if( extschursigntoone(b)==extschursigntoone(c) )
|
---|
1196 | {
|
---|
1197 |
|
---|
1198 | //
|
---|
1199 | // Real eigenvalues: reduce to upper triangular form
|
---|
1200 | //
|
---|
1201 | sab = Math.Sqrt(Math.Abs(b));
|
---|
1202 | sac = Math.Sqrt(Math.Abs(c));
|
---|
1203 | p = extschursign(sab*sac, c);
|
---|
1204 | tau = 1/Math.Sqrt(Math.Abs(b+c));
|
---|
1205 | a = temp+p;
|
---|
1206 | d = temp-p;
|
---|
1207 | b = b-c;
|
---|
1208 | c = 0;
|
---|
1209 | cs1 = sab*tau;
|
---|
1210 | sn1 = sac*tau;
|
---|
1211 | temp = cs*cs1-sn*sn1;
|
---|
1212 | sn = cs*sn1+sn*cs1;
|
---|
1213 | cs = temp;
|
---|
1214 | }
|
---|
1215 | }
|
---|
1216 | else
|
---|
1217 | {
|
---|
1218 | b = -c;
|
---|
1219 | c = 0;
|
---|
1220 | temp = cs;
|
---|
1221 | cs = -sn;
|
---|
1222 | sn = temp;
|
---|
1223 | }
|
---|
1224 | }
|
---|
1225 | }
|
---|
1226 | }
|
---|
1227 | }
|
---|
1228 | }
|
---|
1229 |
|
---|
1230 | //
|
---|
1231 | // Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
|
---|
1232 | //
|
---|
1233 | rt1r = a;
|
---|
1234 | rt2r = d;
|
---|
1235 | if( (double)(c)==(double)(0) )
|
---|
1236 | {
|
---|
1237 | rt1i = 0;
|
---|
1238 | rt2i = 0;
|
---|
1239 | }
|
---|
1240 | else
|
---|
1241 | {
|
---|
1242 | rt1i = Math.Sqrt(Math.Abs(b))*Math.Sqrt(Math.Abs(c));
|
---|
1243 | rt2i = -rt1i;
|
---|
1244 | }
|
---|
1245 | }
|
---|
1246 |
|
---|
1247 |
|
---|
1248 | private static double extschursign(double a,
|
---|
1249 | double b)
|
---|
1250 | {
|
---|
1251 | double result = 0;
|
---|
1252 |
|
---|
1253 | if( (double)(b)>=(double)(0) )
|
---|
1254 | {
|
---|
1255 | result = Math.Abs(a);
|
---|
1256 | }
|
---|
1257 | else
|
---|
1258 | {
|
---|
1259 | result = -Math.Abs(a);
|
---|
1260 | }
|
---|
1261 | return result;
|
---|
1262 | }
|
---|
1263 |
|
---|
1264 |
|
---|
1265 | private static int extschursigntoone(double b)
|
---|
1266 | {
|
---|
1267 | int result = 0;
|
---|
1268 |
|
---|
1269 | if( (double)(b)>=(double)(0) )
|
---|
1270 | {
|
---|
1271 | result = 1;
|
---|
1272 | }
|
---|
1273 | else
|
---|
1274 | {
|
---|
1275 | result = -1;
|
---|
1276 | }
|
---|
1277 | return result;
|
---|
1278 | }
|
---|
1279 | }
|
---|
1280 | }
|
---|