Actual source code: ex116.c

  1: static char help[] = "Test LAPACK routine DSYEV() or DSYEVX(). \n\
  2: Reads PETSc matrix A \n\
  3: then computes selected eigenvalues, and optionally, eigenvectors of \n\
  4: a real generalized symmetric-definite eigenproblem \n\
  5:  A*x = lambda*x \n\
  6: Input parameters include\n\
  7:   -f <input_file> : file to load\n\
  8: e.g. ./ex116 -f $DATAFILESPATH/matrices/small  \n\n";

 10: #include <petscmat.h>
 11: #include <petscblaslapack.h>

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

 15: int main(int argc, char **args)
 16: {
 17:   Mat           A, A_dense;
 18:   Vec          *evecs;
 19:   PetscViewer   fd;                          /* viewer */
 20:   char          file[1][PETSC_MAX_PATH_LEN]; /* input file name */
 21:   PetscBool     flg, TestSYEVX = PETSC_TRUE;
 22:   PetscBool     isSymmetric;
 23:   PetscScalar  *arrayA, *evecs_array, *work, *evals;
 24:   PetscMPIInt   size;
 25:   PetscInt      m, n, i, cklvl = 2;
 26:   PetscBLASInt  nevs, il, iu, in;
 27:   PetscReal     vl, vu, abstol = 1.e-8;
 28:   PetscBLASInt *iwork, *ifail, lwork, lierr, bn;
 29:   PetscReal     tols[2];

 31:   PetscFunctionBeginUser;
 32:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 33:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 34:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

 36:   PetscCall(PetscOptionsHasName(NULL, NULL, "-test_syev", &flg));
 37:   if (flg) TestSYEVX = PETSC_FALSE;

 39:   /* Determine files from which we read the two matrices */
 40:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file[0], sizeof(file[0]), &flg));

 42:   /* Load matrix A */
 43:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[0], FILE_MODE_READ, &fd));
 44:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 45:   PetscCall(MatSetType(A, MATSEQAIJ));
 46:   PetscCall(MatLoad(A, fd));
 47:   PetscCall(PetscViewerDestroy(&fd));
 48:   PetscCall(MatGetSize(A, &m, &n));

 50:   /* Check whether A is symmetric */
 51:   PetscCall(PetscOptionsHasName(NULL, NULL, "-check_symmetry", &flg));
 52:   if (flg) {
 53:     Mat Trans;
 54:     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Trans));
 55:     PetscCall(MatEqual(A, Trans, &isSymmetric));
 56:     PetscCheck(isSymmetric, PETSC_COMM_SELF, PETSC_ERR_USER, "A must be symmetric");
 57:     PetscCall(MatDestroy(&Trans));
 58:   }

 60:   /* Solve eigenvalue problem: A_dense*x = lambda*B*x */
 61:   /*==================================================*/
 62:   /* Convert aij matrix to MatSeqDense for LAPACK */
 63:   PetscCall(MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &A_dense));

 65:   PetscCall(PetscBLASIntCast(8 * n, &lwork));
 66:   PetscCall(PetscBLASIntCast(n, &bn));
 67:   PetscCall(PetscMalloc1(n, &evals));
 68:   PetscCall(PetscMalloc1(lwork, &work));
 69:   PetscCall(MatDenseGetArray(A_dense, &arrayA));

 71:   if (!TestSYEVX) { /* test syev() */
 72:     PetscCall(PetscPrintf(PETSC_COMM_SELF, " LAPACKsyev: compute all %" PetscInt_FMT " eigensolutions...\n", m));
 73:     LAPACKsyev_("V", "U", &bn, arrayA, &bn, evals, work, &lwork, &lierr);
 74:     evecs_array = arrayA;
 75:     PetscCall(PetscBLASIntCast(m, &nevs));
 76:     il = 1;
 77:     PetscCall(PetscBLASIntCast(m, &iu));
 78:   } else { /* test syevx()  */
 79:     il = 1;
 80:     PetscCall(PetscBLASIntCast(0.2 * m, &iu));
 81:     PetscCall(PetscBLASIntCast(n, &in));
 82:     PetscCall(PetscPrintf(PETSC_COMM_SELF, " LAPACKsyevx: compute %" PetscBLASInt_FMT " to %" PetscBLASInt_FMT "-th eigensolutions...\n", il, iu));
 83:     PetscCall(PetscMalloc1(m * n + 1, &evecs_array));
 84:     PetscCall(PetscMalloc1(6 * n + 1, &iwork));
 85:     ifail = iwork + 5 * n;

 87:     /* in the case "I", vl and vu are not referenced */
 88:     vl = 0.0;
 89:     vu = 8.0;
 90:     LAPACKsyevx_("V", "I", "U", &bn, arrayA, &bn, &vl, &vu, &il, &iu, &abstol, &nevs, evals, evecs_array, &in, work, &lwork, iwork, ifail, &lierr);
 91:     PetscCall(PetscFree(iwork));
 92:   }
 93:   PetscCall(MatDenseRestoreArray(A_dense, &arrayA));
 94:   PetscCheck(nevs > 0, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "nev=%" PetscBLASInt_FMT ", no eigensolution has found", nevs);

 96:   /* View eigenvalues */
 97:   PetscCall(PetscOptionsHasName(NULL, NULL, "-eig_view", &flg));
 98:   if (flg) {
 99:     PetscCall(PetscPrintf(PETSC_COMM_SELF, " %" PetscBLASInt_FMT " evals: \n", nevs));
100:     for (i = 0; i < nevs; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT "  %g\n", (PetscInt)(i + il), (double)evals[i]));
101:   }

103:   /* Check residuals and orthogonality */
104:   PetscCall(PetscMalloc1(nevs + 1, &evecs));
105:   for (i = 0; i < nevs; i++) {
106:     PetscCall(VecCreate(PETSC_COMM_SELF, &evecs[i]));
107:     PetscCall(VecSetSizes(evecs[i], PETSC_DECIDE, n));
108:     PetscCall(VecSetFromOptions(evecs[i]));
109:     PetscCall(VecPlaceArray(evecs[i], evecs_array + i * n));
110:   }

112:   tols[0] = tols[1] = PETSC_SQRT_MACHINE_EPSILON;
113:   PetscCall(CkEigenSolutions(cklvl, A, il - 1, iu - 1, evals, evecs, tols));

115:   /* Free work space. */
116:   for (i = 0; i < nevs; i++) PetscCall(VecDestroy(&evecs[i]));
117:   PetscCall(PetscFree(evecs));
118:   PetscCall(MatDestroy(&A_dense));
119:   PetscCall(PetscFree(work));
120:   if (TestSYEVX) PetscCall(PetscFree(evecs_array));

122:   /* Compute SVD: A_dense = U*SIGMA*transpose(V),
123:      JOBU=JOBV='S':  the first min(m,n) columns of U and V are returned in the arrayU and arrayV; */
124:   /*==============================================================================================*/
125:   {
126:     /* Convert aij matrix to MatSeqDense for LAPACK */
127:     PetscScalar *arrayU, *arrayVT, *arrayErr, alpha = 1.0, beta = -1.0;
128:     Mat          Err;
129:     PetscBLASInt minMN, maxMN, im, in;
130:     PetscInt     j;
131:     PetscReal    norm;

133:     PetscCall(MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &A_dense));

135:     PetscCall(PetscBLASIntCast(PetscMin(m, n), &minMN));
136:     PetscCall(PetscBLASIntCast(PetscMax(m, n), &maxMN));
137:     PetscCall(PetscBLASIntCast(5 * minMN + maxMN, &lwork));
138:     PetscCall(PetscMalloc4(m * minMN, &arrayU, m * minMN, &arrayVT, m * minMN, &arrayErr, lwork, &work));

140:     /* Create matrix Err for checking error */
141:     PetscCall(MatCreate(PETSC_COMM_WORLD, &Err));
142:     PetscCall(MatSetSizes(Err, PETSC_DECIDE, PETSC_DECIDE, m, minMN));
143:     PetscCall(MatSetType(Err, MATSEQDENSE));
144:     PetscCall(MatSeqDenseSetPreallocation(Err, (PetscScalar *)arrayErr));

146:     /* Save A to arrayErr for checking accuracy later. arrayA will be destroyed by LAPACKgesvd_() */
147:     PetscCall(MatDenseGetArray(A_dense, &arrayA));
148:     PetscCall(PetscArraycpy(arrayErr, arrayA, m * minMN));

150:     PetscCall(PetscBLASIntCast(m, &im));
151:     PetscCall(PetscBLASIntCast(n, &in));
152:     /* Compute A = U*SIGMA*VT */
153:     LAPACKgesvd_("S", "S", &im, &in, arrayA, &im, evals, arrayU, &minMN, arrayVT, &minMN, work, &lwork, &lierr);
154:     PetscCall(MatDenseRestoreArray(A_dense, &arrayA));
155:     if (!lierr) {
156:       PetscCall(PetscPrintf(PETSC_COMM_SELF, " 1st 10 of %" PetscBLASInt_FMT " singular values: \n", minMN));
157:       for (i = 0; i < 10; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT "  %g\n", i, (double)evals[i]));
158:     } else {
159:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "LAPACKgesvd_ fails!"));
160:     }

162:     /* Check Err = (U*Sigma*V^T - A) using BLASgemm() */
163:     /* U = U*Sigma */
164:     for (j = 0; j < minMN; j++) { /* U[:,j] = sigma[j]*U[:,j] */
165:       for (i = 0; i < m; i++) arrayU[j * m + i] *= evals[j];
166:     }
167:     /* Err = U*VT - A = alpha*U*VT + beta*Err */
168:     BLASgemm_("N", "N", &im, &minMN, &minMN, &alpha, arrayU, &im, arrayVT, &minMN, &beta, arrayErr, &im);
169:     PetscCall(MatNorm(Err, NORM_FROBENIUS, &norm));
170:     PetscCall(PetscPrintf(PETSC_COMM_SELF, " || U*Sigma*VT - A || = %g\n", (double)norm));

172:     PetscCall(PetscFree4(arrayU, arrayVT, arrayErr, work));
173:     PetscCall(PetscFree(evals));
174:     PetscCall(MatDestroy(&A_dense));
175:     PetscCall(MatDestroy(&Err));
176:   }

178:   PetscCall(MatDestroy(&A));
179:   PetscCall(PetscFinalize());
180:   return 0;
181: }
182: /*------------------------------------------------
183:   Check the accuracy of the eigen solution
184:   ----------------------------------------------- */
185: /*
186:   input:
187:      cklvl      - check level:
188:                     1: check residual
189:                     2: 1 and check B-orthogonality locally
190:      A          - matrix
191:      il,iu      - lower and upper index bound of eigenvalues
192:      eval, evec - eigenvalues and eigenvectors stored in this process
193:      tols[0]    - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
194:      tols[1]    - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
195: */
196: PetscErrorCode CkEigenSolutions(PetscInt cklvl, Mat A, PetscInt il, PetscInt iu, PetscReal *eval, Vec *evec, PetscReal *tols)
197: {
198:   PetscInt  i, j, nev;
199:   Vec       vt1, vt2; /* tmp vectors */
200:   PetscReal norm, tmp, dot, norm_max, dot_max;

202:   PetscFunctionBegin;
203:   nev = iu - il;
204:   if (nev <= 0) PetscFunctionReturn(PETSC_SUCCESS);

206:   /*ierr = VecView(evec[0],PETSC_VIEWER_STDOUT_WORLD);*/
207:   PetscCall(VecDuplicate(evec[0], &vt1));
208:   PetscCall(VecDuplicate(evec[0], &vt2));

210:   switch (cklvl) {
211:   case 2:
212:     dot_max = 0.0;
213:     for (i = il; i < iu; i++) {
214:       PetscCall(VecCopy(evec[i], vt1));
215:       for (j = il; j < iu; j++) {
216:         PetscCall(VecDot(evec[j], vt1, &dot));
217:         if (j == i) {
218:           dot = PetscAbsScalar(dot - 1);
219:         } else {
220:           dot = PetscAbsScalar(dot);
221:         }
222:         if (dot > dot_max) dot_max = dot;
223:         if (dot > tols[1]) {
224:           PetscCall(VecNorm(evec[i], NORM_INFINITY, &norm));
225:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "|delta(%" PetscInt_FMT ",%" PetscInt_FMT ")|: %g, norm: %g\n", i, j, (double)dot, (double)norm));
226:         }
227:       }
228:     }
229:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "    max|(x_j^T*x_i) - delta_ji|: %g\n", (double)dot_max));

231:   case 1:
232:     norm_max = 0.0;
233:     for (i = il; i < iu; i++) {
234:       PetscCall(MatMult(A, evec[i], vt1));
235:       PetscCall(VecCopy(evec[i], vt2));
236:       tmp = -eval[i];
237:       PetscCall(VecAXPY(vt1, tmp, vt2));
238:       PetscCall(VecNorm(vt1, NORM_INFINITY, &norm));
239:       norm = PetscAbsScalar(norm);
240:       if (norm > norm_max) norm_max = norm;
241:       /* sniff, and bark if necessary */
242:       if (norm > tols[0]) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  residual violation: %" PetscInt_FMT ", resi: %g\n", i, (double)norm));
243:     }
244:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "    max_resi:                    %g\n", (double)norm_max));
245:     break;
246:   default:
247:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%" PetscInt_FMT " is not supported \n", cklvl));
248:   }
249:   PetscCall(VecDestroy(&vt2));
250:   PetscCall(VecDestroy(&vt1));
251:   PetscFunctionReturn(PETSC_SUCCESS);
252: }

254: /*TEST

256:    build:
257:       requires: !complex

259:    test:
260:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
261:       args: -f ${DATAFILESPATH}/matrices/small
262:       output_file: output/ex116_1.out

264:    test:
265:       suffix: 2
266:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
267:       args: -f ${DATAFILESPATH}/matrices/small -test_syev -check_symmetry

269: TEST*/