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)
 11:   PetscFunctionBeginUser;
 12:   PetscCall(PetscInitialize(&argc, &args, NULL, 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;

 26:   PetscFunctionBeginUser;
 27:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 28:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 29:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

 31:   n    = 100;
 32:   nevs = iu - il;
 33:   PetscCall(PetscMalloc1(3 * n + 1, &D));
 34:   E     = D + n;
 35:   evals = E + n;
 36:   PetscCall(PetscMalloc1(5 * n + 1, &work));
 37:   PetscCall(PetscMalloc1(3 * n + 1, &iwork));
 38:   PetscCall(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:   PetscCall(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);
 50:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_USER, "LAPACKstebz_ fails. info %d", info);

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

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

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

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

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

 94:   /* free space */

 96:   PetscCall(MatDestroy(&T));

 98:   for (i = 0; i < nevs; i++) PetscCall(VecDestroy(&evecs[i]));
 99:   PetscCall(PetscFree(evecs));
100:   PetscCall(PetscFree(D));
101:   PetscCall(PetscFree(work));
102:   PetscCall(PetscFree(iwork));
103:   PetscCall(PetscFree(iblock));
104:   PetscCall(PetscFree(evecs_array));
105:   PetscCall(PetscFree(ifail));
106:   PetscCall(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:   PetscFunctionBegin;
134:   nev = iu - il;
135:   if (nev <= 0) PetscFunctionReturn(PETSC_SUCCESS);

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

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

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

183:   PetscCall(VecDestroy(&vt2));
184:   PetscCall(VecDestroy(&vt1));
185:   PetscFunctionReturn(PETSC_SUCCESS);
186: }