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: }