Actual source code: cgeig.c

petsc-3.14.0 2020-09-29
Report Typos and Errors
```
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;
16:   PetscInt       j,n = cgP->ned;

19:   if (nmax < n) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues");
20:   *neig = n;

22:   PetscArrayzero(c,nmax);
23:   if (!n) {
24:     return(0);
25:   }
26:   d = cgP->d; e = cgP->e; ee = cgP->ee;

28:   /* copy tridiagonal matrix to work space */
29:   for (j=0; j<n; j++) {
30:     r[j]  = PetscRealPart(d[j]);
31:     ee[j] = PetscRealPart(e[j]);
32:   }

34:   LINPACKcgtql1(&n,r,ee,&j);
35:   if (j != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error from tql1(); eispack eigenvalue routine");
36:   PetscSortReal(n,r);
37:   return(0);
38: }

40: PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP ksp,PetscReal *emax,PetscReal *emin)
41: {
42:   KSP_CG      *cgP = (KSP_CG*)ksp->data;
43:   PetscScalar *d,*e;
44:   PetscReal   *dd,*ee;
45:   PetscInt    j,n = cgP->ned;

48:   if (!n) {
49:     *emax = *emin = 1.0;
50:     return(0);
51:   }
52:   d = cgP->d; e = cgP->e; dd = cgP->dd; ee = cgP->ee;

54:   /* copy tridiagonal matrix to work space */
55:   for (j=0; j<n; j++) {
56:     dd[j] = PetscRealPart(d[j]);
57:     ee[j] = PetscRealPart(e[j]);
58:   }

60:   LINPACKcgtql1(&n,dd,ee,&j);
61:   if (j != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error from tql1(); eispack eigenvalue routine");
62:   *emin = dd[0]; *emax = dd[n-1];
63:   return(0);
64: }

66: /* tql1.f -- translated by f2c (version of 25 March 1992  12:58:56).
67:    By Barry Smith on March 27, 1994.
68:    Eispack routine to determine eigenvalues of symmetric
69:    tridiagonal matrix

71:   Note that this routine always uses real numbers (not complex) even if the underlying
72:   matrix is Hermitian. This is because the Lanczos process applied to Hermitian matrices
73:   always produces a real, symmetric tridiagonal matrix.
74: */

76: static PetscReal LINPACKcgpthy(PetscReal*,PetscReal*);

78: static PetscErrorCode LINPACKcgtql1(PetscInt *n,PetscReal *d,PetscReal *e,PetscInt *ierr)
79: {
80:   /* System generated locals */
81:   PetscInt  i__1,i__2;
82:   PetscReal d__1,d__2,c_b10 = 1.0;

84:   /* Local variables */
85:   PetscReal c,f,g,h;
86:   PetscInt  i,j,l,m;
87:   PetscReal p,r,s,c2,c3 = 0.0;
88:   PetscInt  l1,l2;
89:   PetscReal s2 = 0.0;
90:   PetscInt  ii;
91:   PetscReal dl1,el1;
92:   PetscInt  mml;
93:   PetscReal tst1,tst2;

95: /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL1, */
96: /*     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND */
97: /*     WILKINSON. */
98: /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). */

100: /*     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC */
101: /*     TRIDIAGONAL MATRIX BY THE QL METHOD. */

103: /*     ON INPUT */

105: /*        N IS THE ORDER OF THE MATRIX. */

107: /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */

109: /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
110: /*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */

112: /*      ON OUTPUT */

114: /*        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN */
115: /*          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND */
116: /*          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE */
117: /*          THE SMALLEST EIGENVALUES. */

119: /*        E HAS BEEN DESTROYED. */

121: /*        IERR IS SET TO */
122: /*          ZERO       FOR NORMAL RETURN, */
123: /*          J          IF THE J-TH EIGENVALUE HAS NOT BEEN */
124: /*                     DETERMINED AFTER 30 ITERATIONS. */

126: /*     CALLS CGPTHY FOR  DSQRT(A*A + B*B) . */

128: /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
129: /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
130: */

132: /*     THIS VERSION DATED AUGUST 1983. */

134: /*     ------------------------------------------------------------------
135: */
136:   PetscReal ds;

139:   --e;
140:   --d;

142:   *0;
143:   if (*n == 1) goto L1001;

146:   i__1 = *n;
147:   for (i = 2; i <= i__1; ++i) e[i - 1] = e[i];

149:   f     = 0.;
150:   tst1  = 0.;
151:   e[*n] = 0.;

153:   i__1 = *n;
154:   for (l = 1; l <= i__1; ++l) {
155:     j = 0;
156:     d__1 = d[l];
157:     d__2 = e[l];
158:     h = PetscAbsReal(d__1) + PetscAbsReal(d__2);
159:     if (tst1 < h) tst1 = h;
160: /*     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
161:     i__2 = *n;
162:     for (m = l; m <= i__2; ++m) {
163:       d__1 = e[m];
164:       tst2 = tst1 + PetscAbsReal(d__1);
165:       if (tst2 == tst1) goto L120;
166: /*     .......... E(N) IS ALWAYS ZERO,SO THERE IS NO EXIT */
167: /*                THROUGH THE BOTTOM OF THE LOOP .......... */
168:     }
169: L120:
170:     if (m == l) goto L210;
171: L130:
172:     if (j == 30) goto L1000;
173:     ++j;
174: /*     .......... FORM SHIFT .......... */
175:     l1    = l + 1;
176:     l2    = l1 + 1;
177:     g     = d[l];
178:     p     = (d[l1] - g) / (e[l] * 2.);
179:     r     = LINPACKcgpthy(&p,&c_b10);
180:     ds    = 1.0; if (p < 0.0) ds = -1.0;
181:     d[l]  = e[l] / (p + ds*r);
182:     d[l1] = e[l] * (p + ds*r);
183:     dl1   = d[l1];
184:     h     = g - d[l];
185:     if (l2 > *n) goto L145;

187:     i__2 = *n;
188:     for (i = l2; i <= i__2; ++i) d[i] -= h;

190: L145:
191:     f += h;
192: /*     .......... QL TRANSFORMATION .......... */
193:     p   = d[m];
194:     c   = 1.;
195:     c2  = c;
196:     el1 = e[l1];
197:     s   = 0.;
198:     mml = m - l;
199: /*     .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
200:     i__2 = mml;
201:     for (ii = 1; ii <= i__2; ++ii) {
202:       c3       = c2;
203:       c2       = c;
204:       s2       = s;
205:       i        = m - ii;
206:       g        = c * e[i];
207:       h        = c * p;
208:       r        = LINPACKcgpthy(&p,&e[i]);
209:       e[i + 1] = s * r;
210:       s        = e[i] / r;
211:       c        = p / r;
212:       p        = c * d[i] - s * g;
213:       d[i + 1] = h + s * (c * g + s * d[i]);
214:     }

216:     p    = -s * s2 * c3 * el1 * e[l] / dl1;
217:     e[l] = s * p;
218:     d[l] = c * p;
219:     d__1 = e[l];
220:     tst2 = tst1 + PetscAbsReal(d__1);
221:     if (tst2 > tst1) goto L130;
222: L210:
223:     p = d[l] + f;
224: /*     .......... ORDER EIGENVALUES .......... */
225:     if (l == 1) goto L250;
226: /*     .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... */
227:     i__2 = l;
228:     for (ii = 2; ii <= i__2; ++ii) {
229:       i = l + 2 - ii;
230:       if (p >= d[i - 1]) goto L270;
231:       d[i] = d[i - 1];
232:     }

234: L250:
235:     i = 1;
236: L270:
237:     d[i] = p;
238:   }

240:   goto L1001;
241: /*     .......... SET ERROR -- NO CONVERGENCE TO AN */
242: /*                EIGENVALUE AFTER 30 ITERATIONS .......... */
243: L1000:
244:   *l;
245: L1001:
246:   return(0);
247: } /* cgtql1_ */

249: static PetscReal LINPACKcgpthy(PetscReal *a,PetscReal *b)
250: {
251:   /* System generated locals */
252:   PetscReal ret_val,d__1,d__2,d__3;

254:   /* Local variables */
255:   PetscReal p,r,s,t,u;

258: /*     FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW */

261: /* Computing MAX */
262:   d__1 = PetscAbsReal(*a);
263:   d__2 = PetscAbsReal(*b);
264:   p    = PetscMax(d__1,d__2);
265:   if (!p) goto L20;
266: /* Computing MIN */
267:   d__2 = PetscAbsReal(*a);
268:   d__3 = PetscAbsReal(*b);
269: /* Computing 2nd power */
270:   d__1 = PetscMin(d__2,d__3) / p;
271:   r    = d__1 * d__1;
272: L10:
273:   t = r + 4.;
274:   if (t == 4.) goto L20;
275:   s = r / t;
276:   u = s * 2. + 1.;
277:   p = u * p;
278: /* Computing 2nd power */
279:   d__1 = s / u;
280:   r    = d__1 * d__1 * r;
281:   goto L10;
282: L20:
283:   ret_val = p;
284:   PetscFunctionReturn(ret_val);
285: } /* cgpthy_ */

```