Actual source code: cgeig.c
2: /*
3: Code for calculating extreme eigenvalues via the Lanczo method
4: running with CG. Note this only works for symmetric real and Hermitian
5: matrices (not complex matrices that are symmetric).
6: */
7: #include <../src/ksp/ksp/impls/cg/cgimpl.h>
8: static PetscErrorCode LINPACKcgtql1(PetscInt*,PetscReal*,PetscReal*,PetscInt*);
10: PetscErrorCode KSPComputeEigenvalues_CG(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig)
11: {
12: KSP_CG *cgP = (KSP_CG*)ksp->data;
13: PetscScalar *d,*e;
14: PetscReal *ee;
15: PetscInt j,n = ksp->its;
18: *neig = n;
20: PetscArrayzero(c,nmax);
21: if (!n) {
22: return 0;
23: }
24: d = cgP->d; e = cgP->e; ee = cgP->ee;
26: /* copy tridiagonal matrix to work space */
27: for (j=0; j<n; j++) {
28: r[j] = PetscRealPart(d[j]);
29: ee[j] = PetscRealPart(e[j]);
30: }
32: LINPACKcgtql1(&n,r,ee,&j);
34: PetscSortReal(n,r);
35: return 0;
36: }
38: PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP ksp,PetscReal *emax,PetscReal *emin)
39: {
40: KSP_CG *cgP = (KSP_CG*)ksp->data;
41: PetscScalar *d,*e;
42: PetscReal *dd,*ee;
43: PetscInt j,n = ksp->its;
45: if (!n) {
46: *emax = *emin = 1.0;
47: return 0;
48: }
49: d = cgP->d; e = cgP->e; dd = cgP->dd; ee = cgP->ee;
51: /* copy tridiagonal matrix to work space */
52: for (j=0; j<n; j++) {
53: dd[j] = PetscRealPart(d[j]);
54: ee[j] = PetscRealPart(e[j]);
55: }
57: LINPACKcgtql1(&n,dd,ee,&j);
59: *emin = dd[0]; *emax = dd[n-1];
60: return 0;
61: }
63: /* tql1.f -- translated by f2c (version of 25 March 1992 12:58:56).
64: By Barry Smith on March 27, 1994.
65: Eispack routine to determine eigenvalues of symmetric
66: tridiagonal matrix
68: Note that this routine always uses real numbers (not complex) even if the underlying
69: matrix is Hermitian. This is because the Lanczos process applied to Hermitian matrices
70: always produces a real, symmetric tridiagonal matrix.
71: */
73: static PetscReal LINPACKcgpthy(PetscReal*,PetscReal*);
75: static PetscErrorCode LINPACKcgtql1(PetscInt *n,PetscReal *d,PetscReal *e,PetscInt *ierr)
76: {
77: /* System generated locals */
78: PetscInt i__1,i__2;
79: PetscReal d__1,d__2,c_b10 = 1.0;
81: /* Local variables */
82: PetscReal c,f,g,h;
83: PetscInt i,j,l,m;
84: PetscReal p,r,s,c2,c3 = 0.0;
85: PetscInt l1,l2;
86: PetscReal s2 = 0.0;
87: PetscInt ii;
88: PetscReal dl1,el1;
89: PetscInt mml;
90: PetscReal tst1,tst2;
92: /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL1, */
93: /* NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND */
94: /* WILKINSON. */
95: /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). */
97: /* THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC */
98: /* TRIDIAGONAL MATRIX BY THE QL METHOD. */
100: /* ON INPUT */
102: /* N IS THE ORDER OF THE MATRIX. */
104: /* D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
106: /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
107: /* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. */
109: /* ON OUTPUT */
111: /* D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN */
112: /* ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND */
113: /* ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE */
114: /* THE SMALLEST EIGENVALUES. */
116: /* E HAS BEEN DESTROYED. */
118: /* IERR IS SET TO */
119: /* ZERO FOR NORMAL RETURN, */
120: /* J IF THE J-TH EIGENVALUE HAS NOT BEEN */
121: /* DETERMINED AFTER 30 ITERATIONS. */
123: /* CALLS CGPTHY FOR DSQRT(A*A + B*B) . */
125: /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
126: /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
127: */
129: /* THIS VERSION DATED AUGUST 1983. */
131: /* ------------------------------------------------------------------
132: */
133: PetscReal ds;
135: --e;
136: --d;
138: *0;
139: if (*n == 1) goto L1001;
141: i__1 = *n;
142: for (i = 2; i <= i__1; ++i) e[i - 1] = e[i];
144: f = 0.;
145: tst1 = 0.;
146: e[*n] = 0.;
148: i__1 = *n;
149: for (l = 1; l <= i__1; ++l) {
150: j = 0;
151: d__1 = d[l];
152: d__2 = e[l];
153: h = PetscAbsReal(d__1) + PetscAbsReal(d__2);
154: if (tst1 < h) tst1 = h;
155: /* .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
156: i__2 = *n;
157: for (m = l; m <= i__2; ++m) {
158: d__1 = e[m];
159: tst2 = tst1 + PetscAbsReal(d__1);
160: if (tst2 == tst1) goto L120;
161: /* .......... E(N) IS ALWAYS ZERO,SO THERE IS NO EXIT */
162: /* THROUGH THE BOTTOM OF THE LOOP .......... */
163: }
164: L120:
165: if (m == l) goto L210;
166: L130:
167: if (j == 30) goto L1000;
168: ++j;
169: /* .......... FORM SHIFT .......... */
170: l1 = l + 1;
171: l2 = l1 + 1;
172: g = d[l];
173: p = (d[l1] - g) / (e[l] * 2.);
174: r = LINPACKcgpthy(&p,&c_b10);
175: ds = 1.0; if (p < 0.0) ds = -1.0;
176: d[l] = e[l] / (p + ds*r);
177: d[l1] = e[l] * (p + ds*r);
178: dl1 = d[l1];
179: h = g - d[l];
180: if (l2 > *n) goto L145;
182: i__2 = *n;
183: for (i = l2; i <= i__2; ++i) d[i] -= h;
185: L145:
186: f += h;
187: /* .......... QL TRANSFORMATION .......... */
188: p = d[m];
189: c = 1.;
190: c2 = c;
191: el1 = e[l1];
192: s = 0.;
193: mml = m - l;
194: /* .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
195: i__2 = mml;
196: for (ii = 1; ii <= i__2; ++ii) {
197: c3 = c2;
198: c2 = c;
199: s2 = s;
200: i = m - ii;
201: g = c * e[i];
202: h = c * p;
203: r = LINPACKcgpthy(&p,&e[i]);
204: e[i + 1] = s * r;
205: s = e[i] / r;
206: c = p / r;
207: p = c * d[i] - s * g;
208: d[i + 1] = h + s * (c * g + s * d[i]);
209: }
211: p = -s * s2 * c3 * el1 * e[l] / dl1;
212: e[l] = s * p;
213: d[l] = c * p;
214: d__1 = e[l];
215: tst2 = tst1 + PetscAbsReal(d__1);
216: if (tst2 > tst1) goto L130;
217: L210:
218: p = d[l] + f;
219: /* .......... ORDER EIGENVALUES .......... */
220: if (l == 1) goto L250;
221: /* .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... */
222: i__2 = l;
223: for (ii = 2; ii <= i__2; ++ii) {
224: i = l + 2 - ii;
225: if (p >= d[i - 1]) goto L270;
226: d[i] = d[i - 1];
227: }
229: L250:
230: i = 1;
231: L270:
232: d[i] = p;
233: }
235: goto L1001;
236: /* .......... SET ERROR -- NO CONVERGENCE TO AN */
237: /* EIGENVALUE AFTER 30 ITERATIONS .......... */
238: L1000:
239: *l;
240: L1001:
241: return 0;
242: } /* cgtql1_ */
244: static PetscReal LINPACKcgpthy(PetscReal *a,PetscReal *b)
245: {
246: /* System generated locals */
247: PetscReal ret_val,d__1,d__2,d__3;
249: /* Local variables */
250: PetscReal p,r,s,t,u;
252: /* FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW */
254: /* Computing MAX */
255: d__1 = PetscAbsReal(*a);
256: d__2 = PetscAbsReal(*b);
257: p = PetscMax(d__1,d__2);
258: if (!p) goto L20;
259: /* Computing MIN */
260: d__2 = PetscAbsReal(*a);
261: d__3 = PetscAbsReal(*b);
262: /* Computing 2nd power */
263: d__1 = PetscMin(d__2,d__3) / p;
264: r = d__1 * d__1;
265: L10:
266: t = r + 4.;
267: if (t == 4.) goto L20;
268: s = r / t;
269: u = s * 2. + 1.;
270: p = u * p;
271: /* Computing 2nd power */
272: d__1 = s / u;
273: r = d__1 * d__1 * r;
274: goto L10;
275: L20:
276: ret_val = p;
277: return ret_val;
278: } /* cgpthy_ */