Actual source code: ex120.c

  1: static char help[] = "Test LAPACK routine ZHEEV, ZHEEVX, ZHEGV and ZHEGVX. \n\
  2: ZHEEV computes all eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix A. \n\n";

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

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

  9: int main(int argc, char **args)
 10: {
 11:   Mat          A, A_dense, B;
 12:   Vec         *evecs;
 13:   PetscBool    flg, TestZHEEV = PETSC_TRUE, TestZHEEVX = PETSC_FALSE, TestZHEGV = PETSC_FALSE, TestZHEGVX = PETSC_FALSE;
 14:   PetscBool    isSymmetric;
 15:   PetscScalar *arrayA, *arrayB, *evecs_array = NULL, *work;
 16:   PetscReal   *evals, *rwork;
 17:   PetscMPIInt  size;
 18:   PetscInt     m, i, j, cklvl = 2;
 19:   PetscReal    vl, vu, abstol = 1.e-8;
 20:   PetscBLASInt nn, nevs, il, iu, *iwork, *ifail, lwork, lierr, bn, one = 1;
 21:   PetscReal    tols[2];
 22:   PetscScalar  v, sigma2;
 23:   PetscRandom  rctx;
 24:   PetscReal    h2, sigma1 = 100.0;
 25:   PetscInt     dim, Ii, J, n = 6, use_random;

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

 32:   PetscCall(PetscOptionsHasName(NULL, NULL, "-test_zheevx", &flg));
 33:   if (flg) {
 34:     TestZHEEV  = PETSC_FALSE;
 35:     TestZHEEVX = PETSC_TRUE;
 36:   }
 37:   PetscCall(PetscOptionsHasName(NULL, NULL, "-test_zhegv", &flg));
 38:   if (flg) {
 39:     TestZHEEV = PETSC_FALSE;
 40:     TestZHEGV = PETSC_TRUE;
 41:   }
 42:   PetscCall(PetscOptionsHasName(NULL, NULL, "-test_zhegvx", &flg));
 43:   if (flg) {
 44:     TestZHEEV  = PETSC_FALSE;
 45:     TestZHEGVX = PETSC_TRUE;
 46:   }

 48:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
 49:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 50:   dim = n * n;

 52:   PetscCall(MatCreate(PETSC_COMM_SELF, &A));
 53:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim));
 54:   PetscCall(MatSetType(A, MATSEQDENSE));
 55:   PetscCall(MatSetFromOptions(A));
 56:   PetscCall(MatSetUp(A));

 58:   PetscCall(PetscOptionsHasName(NULL, NULL, "-norandom", &flg));
 59:   if (flg) use_random = 0;
 60:   else use_random = 1;
 61:   if (use_random) {
 62:     PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rctx));
 63:     PetscCall(PetscRandomSetFromOptions(rctx));
 64:     PetscCall(PetscRandomSetInterval(rctx, 0.0, PETSC_i));
 65:   } else {
 66:     sigma2 = 10.0 * PETSC_i;
 67:   }
 68:   h2 = 1.0 / ((n + 1) * (n + 1));
 69:   for (Ii = 0; Ii < dim; Ii++) {
 70:     v = -1.0;
 71:     i = Ii / n;
 72:     j = Ii - i * n;
 73:     if (i > 0) {
 74:       J = Ii - n;
 75:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 76:     }
 77:     if (i < n - 1) {
 78:       J = Ii + n;
 79:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 80:     }
 81:     if (j > 0) {
 82:       J = Ii - 1;
 83:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 84:     }
 85:     if (j < n - 1) {
 86:       J = Ii + 1;
 87:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 88:     }
 89:     if (use_random) PetscCall(PetscRandomGetValue(rctx, &sigma2));
 90:     v = 4.0 - sigma1 * h2;
 91:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
 92:   }
 93:   /* make A complex Hermitian */
 94:   v  = sigma2 * h2;
 95:   Ii = 0;
 96:   J  = 1;
 97:   PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 98:   v = -sigma2 * h2;
 99:   PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES));
100:   if (use_random) PetscCall(PetscRandomDestroy(&rctx));
101:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
102:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
103:   m = n = dim;

105:   /* Check whether A is symmetric */
106:   PetscCall(PetscOptionsHasName(NULL, NULL, "-check_symmetry", &flg));
107:   if (flg) {
108:     Mat Trans;
109:     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Trans));
110:     PetscCall(MatEqual(A, Trans, &isSymmetric));
111:     PetscCheck(isSymmetric, PETSC_COMM_SELF, PETSC_ERR_USER, "A must be symmetric");
112:     PetscCall(MatDestroy(&Trans));
113:   }

115:   /* Convert aij matrix to MatSeqDense for LAPACK */
116:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQDENSE, &flg));
117:   if (flg) {
118:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A_dense));
119:   } else {
120:     PetscCall(MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &A_dense));
121:   }

123:   PetscCall(MatCreate(PETSC_COMM_SELF, &B));
124:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, dim, dim));
125:   PetscCall(MatSetType(B, MATSEQDENSE));
126:   PetscCall(MatSetFromOptions(B));
127:   PetscCall(MatSetUp(B));
128:   v = 1.0;
129:   for (Ii = 0; Ii < dim; Ii++) PetscCall(MatSetValues(B, 1, &Ii, 1, &Ii, &v, ADD_VALUES));

131:   /* Solve standard eigenvalue problem: A*x = lambda*x */
132:   /*===================================================*/
133:   PetscCall(PetscBLASIntCast(2 * n, &lwork));
134:   PetscCall(PetscBLASIntCast(n, &bn));
135:   PetscCall(PetscMalloc1(n, &evals));
136:   PetscCall(PetscMalloc1(lwork, &work));
137:   PetscCall(MatDenseGetArray(A_dense, &arrayA));

139:   if (TestZHEEV) { /* test zheev() */
140:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsyev: compute all %" PetscInt_FMT " eigensolutions...\n", m));
141:     PetscCall(PetscMalloc1(3 * n - 2, &rwork));
142:     LAPACKsyev_("V", "U", &bn, arrayA, &bn, evals, work, &lwork, rwork, &lierr);
143:     PetscCall(PetscFree(rwork));

145:     evecs_array = arrayA;
146:     nevs        = m;
147:     il          = 1;
148:     iu          = m;
149:   }
150:   if (TestZHEEVX) {
151:     il = 1;
152:     PetscCall(PetscBLASIntCast(0.2 * m, &iu));
153:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsyevx: compute %d to %d-th eigensolutions...\n", il, iu));
154:     PetscCall(PetscMalloc1(m * n + 1, &evecs_array));
155:     PetscCall(PetscMalloc1(7 * n + 1, &rwork));
156:     PetscCall(PetscMalloc1(5 * n + 1, &iwork));
157:     PetscCall(PetscMalloc1(n + 1, &ifail));

159:     /* in the case "I", vl and vu are not referenced */
160:     vl = 0.0;
161:     vu = 8.0;
162:     PetscCall(PetscBLASIntCast(n, &nn));
163:     LAPACKsyevx_("V", "I", "U", &bn, arrayA, &bn, &vl, &vu, &il, &iu, &abstol, &nevs, evals, evecs_array, &nn, work, &lwork, rwork, iwork, ifail, &lierr);
164:     PetscCall(PetscFree(iwork));
165:     PetscCall(PetscFree(ifail));
166:     PetscCall(PetscFree(rwork));
167:   }
168:   if (TestZHEGV) {
169:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsygv: compute all %" PetscInt_FMT " eigensolutions...\n", m));
170:     PetscCall(PetscMalloc1(3 * n + 1, &rwork));
171:     PetscCall(MatDenseGetArray(B, &arrayB));
172:     LAPACKsygv_(&one, "V", "U", &bn, arrayA, &bn, arrayB, &bn, evals, work, &lwork, rwork, &lierr);
173:     evecs_array = arrayA;
174:     nevs        = m;
175:     il          = 1;
176:     iu          = m;
177:     PetscCall(MatDenseRestoreArray(B, &arrayB));
178:     PetscCall(PetscFree(rwork));
179:   }
180:   if (TestZHEGVX) {
181:     il = 1;
182:     PetscCall(PetscBLASIntCast(0.2 * m, &iu));
183:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsygv: compute %d to %d-th eigensolutions...\n", il, iu));
184:     PetscCall(PetscMalloc1(m * n + 1, &evecs_array));
185:     PetscCall(PetscMalloc1(6 * n + 1, &iwork));
186:     ifail = iwork + 5 * n;
187:     PetscCall(PetscMalloc1(7 * n + 1, &rwork));
188:     PetscCall(MatDenseGetArray(B, &arrayB));
189:     vl = 0.0;
190:     vu = 8.0;
191:     PetscCall(PetscBLASIntCast(n, &nn));
192:     LAPACKsygvx_(&one, "V", "I", "U", &bn, arrayA, &bn, arrayB, &bn, &vl, &vu, &il, &iu, &abstol, &nevs, evals, evecs_array, &nn, work, &lwork, rwork, iwork, ifail, &lierr);
193:     PetscCall(MatDenseRestoreArray(B, &arrayB));
194:     PetscCall(PetscFree(iwork));
195:     PetscCall(PetscFree(rwork));
196:   }
197:   PetscCall(MatDenseRestoreArray(A_dense, &arrayA));
198:   PetscCheck(nevs > 0, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);

200:   /* View evals */
201:   PetscCall(PetscOptionsHasName(NULL, NULL, "-eig_view", &flg));
202:   if (flg) {
203:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %d evals: \n", nevs));
204:     for (i = 0; i < nevs; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "  %g\n", i + il, (double)evals[i]));
205:   }

207:   /* Check residuals and orthogonality */
208:   PetscCall(PetscMalloc1(nevs + 1, &evecs));
209:   for (i = 0; i < nevs; i++) {
210:     PetscCall(VecCreate(PETSC_COMM_SELF, &evecs[i]));
211:     PetscCall(VecSetSizes(evecs[i], PETSC_DECIDE, n));
212:     PetscCall(VecSetFromOptions(evecs[i]));
213:     PetscCall(VecPlaceArray(evecs[i], evecs_array + i * n));
214:   }

216:   tols[0] = PETSC_SQRT_MACHINE_EPSILON;
217:   tols[1] = PETSC_SQRT_MACHINE_EPSILON;
218:   PetscCall(CkEigenSolutions(cklvl, A, il - 1, iu - 1, evals, evecs, tols));
219:   for (i = 0; i < nevs; i++) PetscCall(VecDestroy(&evecs[i]));
220:   PetscCall(PetscFree(evecs));

222:   /* Free work space. */
223:   if (TestZHEEVX || TestZHEGVX) PetscCall(PetscFree(evecs_array));
224:   PetscCall(PetscFree(evals));
225:   PetscCall(PetscFree(work));
226:   PetscCall(MatDestroy(&A_dense));
227:   PetscCall(MatDestroy(&A));
228:   PetscCall(MatDestroy(&B));
229:   PetscCall(PetscFinalize());
230:   return 0;
231: }
232: /*------------------------------------------------
233:   Check the accuracy of the eigen solution
234:   ----------------------------------------------- */
235: /*
236:   input:
237:      cklvl      - check level:
238:                     1: check residual
239:                     2: 1 and check B-orthogonality locally
240:      A          - matrix
241:      il,iu      - lower and upper index bound of eigenvalues
242:      eval, evec - eigenvalues and eigenvectors stored in this process
243:      tols[0]    - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
244:      tols[1]    - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
245: */
246: PetscErrorCode CkEigenSolutions(PetscInt cklvl, Mat A, PetscInt il, PetscInt iu, PetscReal *eval, Vec *evec, PetscReal *tols)
247: {
248:   PetscInt    i, j, nev;
249:   Vec         vt1, vt2; /* tmp vectors */
250:   PetscReal   norm, tmp, norm_max, dot_max, rdot;
251:   PetscScalar dot;

253:   PetscFunctionBegin;
254:   nev = iu - il;
255:   if (nev <= 0) PetscFunctionReturn(PETSC_SUCCESS);

257:   PetscCall(VecDuplicate(evec[0], &vt1));
258:   PetscCall(VecDuplicate(evec[0], &vt2));

260:   switch (cklvl) {
261:   case 2:
262:     dot_max = 0.0;
263:     for (i = il; i < iu; i++) {
264:       PetscCall(VecCopy(evec[i], vt1));
265:       for (j = il; j < iu; j++) {
266:         PetscCall(VecDot(evec[j], vt1, &dot));
267:         if (j == i) {
268:           rdot = PetscAbsScalar(dot - (PetscScalar)1.0);
269:         } else {
270:           rdot = PetscAbsScalar(dot);
271:         }
272:         if (rdot > dot_max) dot_max = rdot;
273:         if (rdot > tols[1]) {
274:           PetscCall(VecNorm(evec[i], NORM_INFINITY, &norm));
275:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "|delta(%" PetscInt_FMT ",%" PetscInt_FMT ")|: %g, norm: %g\n", i, j, (double)rdot, (double)norm));
276:         }
277:       }
278:     }
279:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "    max|(x_j^T*x_i) - delta_ji|: %g\n", (double)dot_max));

281:   case 1:
282:     norm_max = 0.0;
283:     for (i = il; i < iu; i++) {
284:       PetscCall(MatMult(A, evec[i], vt1));
285:       PetscCall(VecCopy(evec[i], vt2));
286:       tmp = -eval[i];
287:       PetscCall(VecAXPY(vt1, tmp, vt2));
288:       PetscCall(VecNorm(vt1, NORM_INFINITY, &norm));
289:       norm = PetscAbs(norm);
290:       if (norm > norm_max) norm_max = norm;
291:       /* sniff, and bark if necessary */
292:       if (norm > tols[0]) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  residual violation: %" PetscInt_FMT ", resi: %g\n", i, (double)norm));
293:     }
294:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "    max_resi:                    %g\n", (double)norm_max));
295:     break;
296:   default:
297:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%" PetscInt_FMT " is not supported \n", cklvl));
298:   }
299:   PetscCall(VecDestroy(&vt2));
300:   PetscCall(VecDestroy(&vt1));
301:   PetscFunctionReturn(PETSC_SUCCESS);
302: }

304: /*TEST

306:    build:
307:       requires: complex

309:    test:

311:    test:
312:       suffix: 2
313:       args: -test_zheevx

315:    test:
316:       suffix: 3
317:       args: -test_zhegv

319:    test:
320:       suffix: 4
321:       args: -test_zhegvx

323: TEST*/