Actual source code: cgeig.c

      1: #define PETSCKSP_DLL

      3: /*
      4:       Code for calculating extreme eigenvalues via the Lanczo method
      5:    running with CG. Note this only works for symmetric real and Hermitian
      6:    matrices (not complex matrices that are symmetric).
      7: */
     8:  #include src/ksp/ksp/impls/cg/cgctx.h
      9: static PetscErrorCode LINPACKcgtql1(PetscInt*,PetscReal *,PetscReal *,PetscInt *);

     13: PetscErrorCode KSPComputeEigenvalues_CG(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig)
     14: {
     15:   KSP_CG         *cgP = (KSP_CG*)ksp->data;
     16:   PetscScalar    *d,*e;
     17:   PetscReal      *ee;
     19:   PetscInt       j,n = ksp->its;

     22:   if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues");
     23:   *neig = n;

     25:   PetscMemzero(c,nmax*sizeof(PetscReal));
     26:   if (!n) {
     27:     *r = 0.0;
     28:     return(0);
     29:   }
     30:   d = cgP->d; e = cgP->e; ee = cgP->ee;

     32:   /* copy tridiagonal matrix to work space */
     33:   for (j=0; j<n ; j++) {
     34:     r[j]  = PetscRealPart(d[j]);
     35:     ee[j] = PetscRealPart(e[j]);
     36:   }

     38:   LINPACKcgtql1(&n,r,ee,&j);
     39:   if (j != 0) SETERRQ(PETSC_ERR_LIB,"Error from tql1(); eispack eigenvalue routine");
     40:   PetscSortReal(n,r);
     41:   return(0);
     42: }

     46: PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP ksp,PetscReal *emax,PetscReal *emin)
     47: {
     48:   KSP_CG      *cgP = (KSP_CG*)ksp->data;
     49:   PetscScalar *d,*e;
     50:   PetscReal   *dd,*ee;
     51:   PetscInt    j,n = ksp->its;

     54:   if (!n) {
     55:     *emax = *emin = 1.0;
     56:     return(0);
     57:   }
     58:   d = cgP->d; e = cgP->e; dd = cgP->dd; ee = cgP->ee;

     60:   /* copy tridiagonal matrix to work space */
     61:   for (j=0; j<n ; j++) {
     62:     dd[j] = PetscRealPart(d[j]);
     63:     ee[j] = PetscRealPart(e[j]);
     64:   }

     66:   LINPACKcgtql1(&n,dd,ee,&j);
     67:   if (j != 0) SETERRQ(PETSC_ERR_LIB,"Error from tql1(); eispack eigenvalue routine");
     68:   *emin = dd[0]; *emax = dd[n-1];
     69:   return(0);
     70: }

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

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

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

     86: static PetscErrorCode LINPACKcgtql1(PetscInt *n,PetscReal *d,PetscReal *e,PetscInt *ierr)
     87: {
     88:     /* System generated locals */
     89:     PetscInt  i__1,i__2;
     90:     PetscReal d__1,d__2,c_b10 = 1.0;

     92:     /* Local variables */
     93:      PetscReal c,f,g,h;
     94:      PetscInt  i,j,l,m;
     95:      PetscReal p,r,s,c2,c3 = 0.0;
     96:      PetscInt  l1,l2;
     97:      PetscReal s2 = 0.0;
     98:      PetscInt  ii;
     99:      PetscReal dl1,el1;
    100:      PetscInt  mml;
    101:      PetscReal tst1,tst2;

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

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

    111: /*     ON INPUT */

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

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

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

    120: /*      ON OUTPUT */

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

    127: /*        E HAS BEEN DESTROYED. */

    129: /*        IERR IS SET TO */
    130: /*          ZERO       FOR NORMAL RETURN, */
    131: /*          J          IF THE J-TH EIGENVALUE HAS NOT BEEN */
    132: /*                     DETERMINED AFTER 30 ITERATIONS. */

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

    136: /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
    137: /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
    138: */

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

    142: /*     ------------------------------------------------------------------
    143: */
    144:     PetscReal ds;

    147:     --e;
    148:     --d;

    150:     *0;
    151:     if (*n == 1) {
    152:         goto L1001;
    153:     }

    155:     i__1 = *n;
    156:     for (i = 2; i <= i__1; ++i) {
    157:         e[i - 1] = e[i];
    158:     }

    160:     f = 0.;
    161:     tst1 = 0.;
    162:     e[*n] = 0.;

    164:     i__1 = *n;
    165:     for (l = 1; l <= i__1; ++l) {
    166:         j = 0;
    167:         h = (d__1 = d[l],PetscAbsReal(d__1)) + (d__2 = e[l],PetscAbsReal(d__2));
    168:         if (tst1 < h) {
    169:             tst1 = h;
    170:         }
    171: /*     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
    172:         i__2 = *n;
    173:         for (m = l; m <= i__2; ++m) {
    174:             tst2 = tst1 + (d__1 = e[m],PetscAbsReal(d__1));
    175:             if (tst2 == tst1) {
    176:                 goto L120;
    177:             }
    178: /*     .......... E(N) IS ALWAYS ZERO,SO THERE IS NO EXIT */
    179: /*                THROUGH THE BOTTOM OF THE LOOP .......... */
    180:         }
    181: L120:
    182:         if (m == l) {
    183:             goto L210;
    184:         }
    185: L130:
    186:         if (j == 30) {
    187:             goto L1000;
    188:         }
    189:         ++j;
    190: /*     .......... FORM SHIFT .......... */
    191:         l1 = l + 1;
    192:         l2 = l1 + 1;
    193:         g = d[l];
    194:         p = (d[l1] - g) / (e[l] * 2.);
    195:         r = LINPACKcgpthy(&p,&c_b10);
    196:         ds = 1.0; if (p < 0.0) ds = -1.0;
    197:         d[l] = e[l] / (p + ds*r);
    198:         d[l1] = e[l] * (p + ds*r);
    199:         dl1 = d[l1];
    200:         h = g - d[l];
    201:         if (l2 > *n) {
    202:             goto L145;
    203:         }

    205:         i__2 = *n;
    206:         for (i = l2; i <= i__2; ++i) {
    207:             d[i] -= h;
    208:         }

    210: L145:
    211:         f += h;
    212: /*     .......... QL TRANSFORMATION .......... */
    213:         p = d[m];
    214:         c = 1.;
    215:         c2 = c;
    216:         el1 = e[l1];
    217:         s = 0.;
    218:         mml = m - l;
    219: /*     .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
    220:         i__2 = mml;
    221:         for (ii = 1; ii <= i__2; ++ii) {
    222:             c3 = c2;
    223:             c2 = c;
    224:             s2 = s;
    225:             i = m - ii;
    226:             g = c * e[i];
    227:             h = c * p;
    228:             r = LINPACKcgpthy(&p,&e[i]);
    229:             e[i + 1] = s * r;
    230:             s = e[i] / r;
    231:             c = p / r;
    232:             p = c * d[i] - s * g;
    233:             d[i + 1] = h + s * (c * g + s * d[i]);
    234:         }
    235:
    236:         p = -s * s2 * c3 * el1 * e[l] / dl1;
    237:         e[l] = s * p;
    238:         d[l] = c * p;
    239:         tst2 = tst1 + (d__1 = e[l],PetscAbsReal(d__1));
    240:         if (tst2 > tst1) {
    241:             goto L130;
    242:         }
    243: L210:
    244:         p = d[l] + f;
    245: /*     .......... ORDER EIGENVALUES .......... */
    246:         if (l == 1) {
    247:             goto L250;
    248:         }
    249: /*     .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... */
    250:         i__2 = l;
    251:         for (ii = 2; ii <= i__2; ++ii) {
    252:             i = l + 2 - ii;
    253:             if (p >= d[i - 1]) {
    254:                 goto L270;
    255:             }
    256:             d[i] = d[i - 1];
    257:         }

    259: L250:
    260:         i = 1;
    261: L270:
    262:         d[i] = p;
    263:     }

    265:     goto L1001;
    266: /*     .......... SET ERROR -- NO CONVERGENCE TO AN */
    267: /*                EIGENVALUE AFTER 30 ITERATIONS .......... */
    268: L1000:
    269:     *l;
    270: L1001:
    271:     return(0);
    272: } /* cgtql1_ */

    276: static PetscReal LINPACKcgpthy(PetscReal *a,PetscReal *b)
    277: {
    278:     /* System generated locals */
    279:     PetscReal ret_val,d__1,d__2,d__3;
    280:
    281:     /* Local variables */
    282:     PetscReal p,r,s,t,u;

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

    288: /* Computing MAX */
    289:     d__1 = PetscAbsReal(*a),d__2 = PetscAbsReal(*b);
    290:     p = PetscMax(d__1,d__2);
    291:     if (!p) {
    292:         goto L20;
    293:     }
    294: /* Computing MIN */
    295:     d__2 = PetscAbsReal(*a),d__3 = PetscAbsReal(*b);
    296: /* Computing 2nd power */
    297:     d__1 = PetscMin(d__2,d__3) / p;
    298:     r = d__1 * d__1;
    299: L10:
    300:     t = r + 4.;
    301:     if (t == 4.) {
    302:         goto L20;
    303:     }
    304:     s = r / t;
    305:     u = s * 2. + 1.;
    306:     p = u * p;
    307: /* Computing 2nd power */
    308:     d__1 = s / u;
    309:     r = d__1 * d__1 * r;
    310:     goto L10;
    311: L20:
    312:     ret_val = p;
    313:     PetscFunctionReturn(ret_val);
    314: } /* cgpthy_ */