Actual source code: ex118.c

  1: static char help[] = "Test LAPACK routine DSTEBZ() and DTEIN().  \n\n";

  3: #include <petscmat.h>
  4: #include <petscblaslapack.h>

  6: extern PetscErrorCode CkEigenSolutions(PetscInt, Mat, PetscInt, PetscInt, PetscScalar *, Vec *, PetscReal *);

  8: int main(int argc, char **args)
  9: {
 10: #if defined(PETSC_USE_COMPLEX) || defined(PETSC_MISSING_LAPACK_STEBZ) || defined(PETSC_MISSING_LAPACK_STEIN)
 12:   PetscInitialize(&argc, &args, (char *)0, help);
 13:   SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, "This example requires LAPACK routines dstebz and stien and real numbers");
 14: #else
 15:   PetscReal   *work, tols[2];
 16:   PetscInt     i, j;
 17:   PetscBLASInt n, il = 1, iu = 5, *iblock, *isplit, *iwork, nevs, *ifail, cklvl = 2;
 18:   PetscMPIInt  size;
 19:   PetscBool    flg;
 20:   Vec         *evecs;
 21:   PetscScalar *evecs_array, *D, *E, *evals;
 22:   Mat          T;
 23:   PetscReal    vl = 0.0, vu = 4.0, tol = 1000 * PETSC_MACHINE_EPSILON;
 24:   PetscBLASInt nsplit, info;

 27:   PetscInitialize(&argc, &args, (char *)0, help);
 28:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 31:   n    = 100;
 32:   nevs = iu - il;
 33:   PetscMalloc1(3 * n + 1, &D);
 34:   E     = D + n;
 35:   evals = E + n;
 36:   PetscMalloc1(5 * n + 1, &work);
 37:   PetscMalloc1(3 * n + 1, &iwork);
 38:   PetscMalloc1(3 * n + 1, &iblock);
 39:   isplit = iblock + n;

 41:   /* Set symmetric tridiagonal matrix */
 42:   for (i = 0; i < n; i++) {
 43:     D[i] = 2.0;
 44:     E[i] = 1.0;
 45:   }

 47:   /* Solve eigenvalue problem: A*evec = eval*evec */
 48:   PetscPrintf(PETSC_COMM_SELF, " LAPACKstebz_: compute %d eigenvalues...\n", nevs);
 49:   LAPACKstebz_("I", "E", &n, &vl, &vu, &il, &iu, &tol, (PetscReal *)D, (PetscReal *)E, &nevs, &nsplit, (PetscReal *)evals, iblock, isplit, work, iwork, &info);

 52:   PetscPrintf(PETSC_COMM_SELF, " LAPACKstein_: compute %d found eigenvectors...\n", nevs);
 53:   PetscMalloc1(n * nevs, &evecs_array);
 54:   PetscMalloc1(nevs, &ifail);
 55:   LAPACKstein_(&n, (PetscReal *)D, (PetscReal *)E, &nevs, (PetscReal *)evals, iblock, isplit, evecs_array, &n, work, iwork, ifail, &info);
 57:   /* View evals */
 58:   PetscOptionsHasName(NULL, NULL, "-eig_view", &flg);
 59:   if (flg) {
 60:     PetscPrintf(PETSC_COMM_SELF, " %d evals: \n", nevs);
 61:     for (i = 0; i < nevs; i++) PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT "  %g\n", i, (double)evals[i]);
 62:   }

 64:   /* Check residuals and orthogonality */
 65:   MatCreate(PETSC_COMM_SELF, &T);
 66:   MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, n, n);
 67:   MatSetType(T, MATSBAIJ);
 68:   MatSetFromOptions(T);
 69:   MatSetUp(T);
 70:   for (i = 0; i < n; i++) {
 71:     MatSetValues(T, 1, &i, 1, &i, &D[i], INSERT_VALUES);
 72:     if (i != n - 1) {
 73:       j = i + 1;
 74:       MatSetValues(T, 1, &i, 1, &j, &E[i], INSERT_VALUES);
 75:     }
 76:   }
 77:   MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY);
 78:   MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY);

 80:   PetscMalloc1(nevs + 1, &evecs);
 81:   for (i = 0; i < nevs; i++) {
 82:     VecCreate(PETSC_COMM_SELF, &evecs[i]);
 83:     VecSetSizes(evecs[i], PETSC_DECIDE, n);
 84:     VecSetFromOptions(evecs[i]);
 85:     VecPlaceArray(evecs[i], evecs_array + i * n);
 86:   }

 88:   tols[0] = 1.e-8;
 89:   tols[1] = 1.e-8;
 90:   CkEigenSolutions(cklvl, T, il - 1, iu - 1, evals, evecs, tols);

 92:   for (i = 0; i < nevs; i++) VecResetArray(evecs[i]);

 94:   /* free space */

 96:   MatDestroy(&T);

 98:   for (i = 0; i < nevs; i++) VecDestroy(&evecs[i]);
 99:   PetscFree(evecs);
100:   PetscFree(D);
101:   PetscFree(work);
102:   PetscFree(iwork);
103:   PetscFree(iblock);
104:   PetscFree(evecs_array);
105:   PetscFree(ifail);
106:   PetscFinalize();
107:   return 0;
108: #endif
109: }
110: /*------------------------------------------------
111:   Check the accuracy of the eigen solution
112:   ----------------------------------------------- */
113: /*
114:   input:
115:      cklvl      - check level:
116:                     1: check residual
117:                     2: 1 and check B-orthogonality locally
118:      A          - matrix
119:      il,iu      - lower and upper index bound of eigenvalues
120:      eval, evec - eigenvalues and eigenvectors stored in this process
121:      tols[0]    - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
122:      tols[1]    - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
123: */
124: #undef DEBUG_CkEigenSolutions
125: PetscErrorCode CkEigenSolutions(PetscInt cklvl, Mat A, PetscInt il, PetscInt iu, PetscScalar *eval, Vec *evec, PetscReal *tols)
126: {
127:   PetscInt    ierr, i, j, nev;
128:   Vec         vt1, vt2; /* tmp vectors */
129:   PetscReal   norm, norm_max;
130:   PetscScalar dot, tmp;
131:   PetscReal   dot_max;

133:   nev = iu - il;
134:   if (nev <= 0) return 0;

136:   VecDuplicate(evec[0], &vt1);
137:   VecDuplicate(evec[0], &vt2);

139:   switch (cklvl) {
140:   case 2:
141:     dot_max = 0.0;
142:     for (i = il; i < iu; i++) {
143:       VecCopy(evec[i], vt1);
144:       for (j = il; j < iu; j++) {
145:         VecDot(evec[j], vt1, &dot);
146:         if (j == i) {
147:           dot = PetscAbsScalar(dot - (PetscScalar)1.0);
148:         } else {
149:           dot = PetscAbsScalar(dot);
150:         }
151:         if (PetscAbsScalar(dot) > dot_max) dot_max = PetscAbsScalar(dot);
152: #if defined(DEBUG_CkEigenSolutions)
153:         if (dot > tols[1]) {
154:           VecNorm(evec[i], NORM_INFINITY, &norm);
155:           PetscPrintf(PETSC_COMM_SELF, "|delta(%d,%d)|: %g, norm: %d\n", i, j, (double)dot, (double)norm);
156:         }
157: #endif
158:       }
159:     }
160:     PetscPrintf(PETSC_COMM_SELF, "    max|(x_j^T*x_i) - delta_ji|: %g\n", (double)dot_max);

162:   case 1:
163:     norm_max = 0.0;
164:     for (i = il; i < iu; i++) {
165:       MatMult(A, evec[i], vt1);
166:       VecCopy(evec[i], vt2);
167:       tmp = -eval[i];
168:       VecAXPY(vt1, tmp, vt2);
169:       VecNorm(vt1, NORM_INFINITY, &norm);
170:       norm = PetscAbsReal(norm);
171:       if (norm > norm_max) norm_max = norm;
172: #if defined(DEBUG_CkEigenSolutions)
173:       if (norm > tols[0]) PetscPrintf(PETSC_COMM_SELF, "  residual violation: %d, resi: %g\n", i, norm);
174: #endif
175:     }
176:     PetscPrintf(PETSC_COMM_SELF, "    max_resi:                    %g\n", (double)norm_max);
177:     break;
178:   default:
179:     PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%d is not supported \n", cklvl);
180:   }

182:   VecDestroy(&vt2);
183:   VecDestroy(&vt1);
184:   return 0;
185: }