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>
9: /* tql1.f -- translated by f2c (version of 25 March 1992 12:58:56).
10: By Barry Smith on March 27, 1994.
11: Eispack routine to determine eigenvalues of symmetric
12: tridiagonal matrix
14: Note that this routine always uses real numbers (not complex) even if the underlying
15: matrix is Hermitian. This is because the Lanczos process applied to Hermitian matrices
16: always produces a real, symmetric tridiagonal matrix.
17: */
19: static PetscReal LINPACKcgpthy(PetscReal *a, PetscReal *b)
20: {
21: /* System generated locals */
22: PetscReal d__1, d__2, d__3;
24: /* Local variables */
25: PetscReal p, r, s, t, u;
27: /* FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW */
29: /* Computing MAX */
30: d__1 = PetscAbsReal(*a);
31: d__2 = PetscAbsReal(*b);
32: p = PetscMax(d__1, d__2);
33: if (!p) goto L20;
34: /* Computing MIN */
35: d__2 = PetscAbsReal(*a);
36: d__3 = PetscAbsReal(*b);
37: /* Computing 2nd power */
38: d__1 = PetscMin(d__2, d__3) / p;
39: r = d__1 * d__1;
40: L10:
41: t = r + 4.;
42: if (t == 4.) goto L20;
43: s = r / t;
44: u = s * 2. + 1.;
45: p = u * p;
46: /* Computing 2nd power */
47: d__1 = s / u;
48: r = d__1 * d__1 * r;
49: goto L10;
50: L20:
51: return p;
52: } /* cgpthy_ */
54: static PetscErrorCode LINPACKcgtql1(PetscInt *n, PetscReal *d, PetscReal *e, PetscInt *ierr)
55: {
56: /* System generated locals */
57: PetscInt i__1, i__2;
58: PetscReal d__1, d__2, c_b10 = 1.0;
60: /* Local variables */
61: PetscReal c, f, g, h;
62: PetscInt i, j, l, m;
63: PetscReal p, r, s, c2, c3 = 0.0;
64: PetscInt l1, l2;
65: PetscReal s2 = 0.0;
66: PetscInt ii;
67: PetscReal dl1, el1;
68: PetscInt mml;
69: PetscReal tst1, tst2;
71: /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL1, */
72: /* NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND */
73: /* WILKINSON. */
74: /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). */
76: /* THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC */
77: /* TRIDIAGONAL MATRIX BY THE QL METHOD. */
79: /* ON INPUT */
81: /* N IS THE ORDER OF THE MATRIX. */
83: /* D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
85: /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
86: /* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. */
88: /* ON OUTPUT */
90: /* D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN */
91: /* ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND */
92: /* ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE */
93: /* THE SMALLEST EIGENVALUES. */
95: /* E HAS BEEN DESTROYED. */
97: /* IERR IS SET TO */
98: /* ZERO FOR NORMAL RETURN, */
99: /* J IF THE J-TH EIGENVALUE HAS NOT BEEN */
100: /* DETERMINED AFTER 30 ITERATIONS. */
102: /* CALLS CGPTHY FOR DSQRT(A*A + B*B) . */
104: /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
105: /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
106: */
108: /* THIS VERSION DATED AUGUST 1983. */
110: /* ------------------------------------------------------------------
111: */
112: PetscReal ds;
114: PetscFunctionBegin;
115: --e;
116: --d;
118: *ierr = 0;
119: if (*n == 1) goto L1001;
121: i__1 = *n;
122: for (i = 2; i <= i__1; ++i) e[i - 1] = e[i];
124: f = 0.;
125: tst1 = 0.;
126: e[*n] = 0.;
128: i__1 = *n;
129: for (l = 1; l <= i__1; ++l) {
130: j = 0;
131: d__1 = d[l];
132: d__2 = e[l];
133: h = PetscAbsReal(d__1) + PetscAbsReal(d__2);
134: if (tst1 < h) tst1 = h;
135: /* .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
136: i__2 = *n;
137: for (m = l; m <= i__2; ++m) {
138: d__1 = e[m];
139: tst2 = tst1 + PetscAbsReal(d__1);
140: if (tst2 == tst1) goto L120;
141: /* .......... E(N) IS ALWAYS ZERO,SO THERE IS NO EXIT */
142: /* THROUGH THE BOTTOM OF THE LOOP .......... */
143: }
144: L120:
145: if (m == l) goto L210;
146: L130:
147: if (j == 30) goto L1000;
148: ++j;
149: /* .......... FORM SHIFT .......... */
150: l1 = l + 1;
151: l2 = l1 + 1;
152: g = d[l];
153: p = (d[l1] - g) / (e[l] * 2.);
154: r = LINPACKcgpthy(&p, &c_b10);
155: ds = 1.0;
156: if (p < 0.0) ds = -1.0;
157: d[l] = e[l] / (p + ds * r);
158: d[l1] = e[l] * (p + ds * r);
159: dl1 = d[l1];
160: h = g - d[l];
161: if (l2 > *n) goto L145;
163: i__2 = *n;
164: for (i = l2; i <= i__2; ++i) d[i] -= h;
166: L145:
167: f += h;
168: /* .......... QL TRANSFORMATION .......... */
169: p = d[m];
170: c = 1.;
171: c2 = c;
172: el1 = e[l1];
173: s = 0.;
174: mml = m - l;
175: /* .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
176: i__2 = mml;
177: for (ii = 1; ii <= i__2; ++ii) {
178: c3 = c2;
179: c2 = c;
180: s2 = s;
181: i = m - ii;
182: g = c * e[i];
183: h = c * p;
184: r = LINPACKcgpthy(&p, &e[i]);
185: e[i + 1] = s * r;
186: s = e[i] / r;
187: c = p / r;
188: p = c * d[i] - s * g;
189: d[i + 1] = h + s * (c * g + s * d[i]);
190: }
192: p = -s * s2 * c3 * el1 * e[l] / dl1;
193: e[l] = s * p;
194: d[l] = c * p;
195: d__1 = e[l];
196: tst2 = tst1 + PetscAbsReal(d__1);
197: if (tst2 > tst1) goto L130;
198: L210:
199: p = d[l] + f;
200: /* .......... ORDER EIGENVALUES .......... */
201: if (l == 1) goto L250;
202: /* .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... */
203: i__2 = l;
204: for (ii = 2; ii <= i__2; ++ii) {
205: i = l + 2 - ii;
206: if (p >= d[i - 1]) goto L270;
207: d[i] = d[i - 1];
208: }
210: L250:
211: i = 1;
212: L270:
213: d[i] = p;
214: }
216: goto L1001;
217: /* .......... SET ERROR -- NO CONVERGENCE TO AN */
218: /* EIGENVALUE AFTER 30 ITERATIONS .......... */
219: L1000:
220: *ierr = l;
221: L1001:
222: PetscFunctionReturn(PETSC_SUCCESS);
223: } /* cgtql1_ */
225: PetscErrorCode KSPComputeEigenvalues_CG(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig)
226: {
227: KSP_CG *cgP = (KSP_CG *)ksp->data;
228: PetscScalar *d, *e;
229: PetscReal *ee;
230: PetscInt j, n = ksp->its;
232: PetscFunctionBegin;
233: PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
234: *neig = n;
236: PetscCall(PetscArrayzero(c, nmax));
237: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
238: d = cgP->d;
239: e = cgP->e;
240: ee = cgP->ee;
242: /* copy tridiagonal matrix to work space */
243: for (j = 0; j < n; j++) {
244: r[j] = PetscRealPart(d[j]);
245: ee[j] = PetscRealPart(e[j]);
246: }
248: PetscCall(LINPACKcgtql1(&n, r, ee, &j));
249: PetscCheck(j == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error from tql1(); eispack eigenvalue routine");
250: PetscCall(PetscSortReal(n, r));
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP ksp, PetscReal *emax, PetscReal *emin)
255: {
256: KSP_CG *cgP = (KSP_CG *)ksp->data;
257: PetscScalar *d, *e;
258: PetscReal *dd, *ee;
259: PetscInt j, n = ksp->its;
261: PetscFunctionBegin;
262: if (!n) {
263: *emax = *emin = 1.0;
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
266: d = cgP->d;
267: e = cgP->e;
268: dd = cgP->dd;
269: ee = cgP->ee;
271: /* copy tridiagonal matrix to work space */
272: for (j = 0; j < n; j++) {
273: dd[j] = PetscRealPart(d[j]);
274: ee[j] = PetscRealPart(e[j]);
275: }
277: PetscCall(LINPACKcgtql1(&n, dd, ee, &j));
278: PetscCheck(j == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error from tql1(); eispack eigenvalue routine");
279: *emin = dd[0];
280: *emax = dd[n - 1];
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }