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;
28: PetscInitialize(&argc, &args, (char *)0, help);
29: MPI_Comm_size(PETSC_COMM_WORLD, &size);
32: PetscOptionsHasName(NULL, NULL, "-test_zheevx", &flg);
33: if (flg) {
34: TestZHEEV = PETSC_FALSE;
35: TestZHEEVX = PETSC_TRUE;
36: }
37: PetscOptionsHasName(NULL, NULL, "-test_zhegv", &flg);
38: if (flg) {
39: TestZHEEV = PETSC_FALSE;
40: TestZHEGV = PETSC_TRUE;
41: }
42: PetscOptionsHasName(NULL, NULL, "-test_zhegvx", &flg);
43: if (flg) {
44: TestZHEEV = PETSC_FALSE;
45: TestZHEGVX = PETSC_TRUE;
46: }
48: PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL);
49: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
50: dim = n * n;
52: MatCreate(PETSC_COMM_SELF, &A);
53: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim);
54: MatSetType(A, MATSEQDENSE);
55: MatSetFromOptions(A);
56: MatSetUp(A);
58: PetscOptionsHasName(NULL, NULL, "-norandom", &flg);
59: if (flg) use_random = 0;
60: else use_random = 1;
61: if (use_random) {
62: PetscRandomCreate(PETSC_COMM_SELF, &rctx);
63: PetscRandomSetFromOptions(rctx);
64: 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: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
76: }
77: if (i < n - 1) {
78: J = Ii + n;
79: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
80: }
81: if (j > 0) {
82: J = Ii - 1;
83: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
84: }
85: if (j < n - 1) {
86: J = Ii + 1;
87: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
88: }
89: if (use_random) PetscRandomGetValue(rctx, &sigma2);
90: v = 4.0 - sigma1 * h2;
91: 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: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
98: v = -sigma2 * h2;
99: MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES);
100: if (use_random) PetscRandomDestroy(&rctx);
101: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
102: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
103: m = n = dim;
105: /* Check whether A is symmetric */
106: PetscOptionsHasName(NULL, NULL, "-check_symmetry", &flg);
107: if (flg) {
108: Mat Trans;
109: MatTranspose(A, MAT_INITIAL_MATRIX, &Trans);
110: MatEqual(A, Trans, &isSymmetric);
112: MatDestroy(&Trans);
113: }
115: /* Convert aij matrix to MatSeqDense for LAPACK */
116: PetscObjectTypeCompare((PetscObject)A, MATSEQDENSE, &flg);
117: if (flg) {
118: MatDuplicate(A, MAT_COPY_VALUES, &A_dense);
119: } else {
120: MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &A_dense);
121: }
123: MatCreate(PETSC_COMM_SELF, &B);
124: MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, dim, dim);
125: MatSetType(B, MATSEQDENSE);
126: MatSetFromOptions(B);
127: MatSetUp(B);
128: v = 1.0;
129: for (Ii = 0; Ii < dim; Ii++) MatSetValues(B, 1, &Ii, 1, &Ii, &v, ADD_VALUES);
131: /* Solve standard eigenvalue problem: A*x = lambda*x */
132: /*===================================================*/
133: PetscBLASIntCast(2 * n, &lwork);
134: PetscBLASIntCast(n, &bn);
135: PetscMalloc1(n, &evals);
136: PetscMalloc1(lwork, &work);
137: MatDenseGetArray(A_dense, &arrayA);
139: if (TestZHEEV) { /* test zheev() */
140: PetscPrintf(PETSC_COMM_WORLD, " LAPACKsyev: compute all %" PetscInt_FMT " eigensolutions...\n", m);
141: PetscMalloc1(3 * n - 2, &rwork);
142: LAPACKsyev_("V", "U", &bn, arrayA, &bn, evals, work, &lwork, rwork, &lierr);
143: PetscFree(rwork);
145: evecs_array = arrayA;
146: nevs = m;
147: il = 1;
148: iu = m;
149: }
150: if (TestZHEEVX) {
151: il = 1;
152: PetscBLASIntCast((0.2 * m), &iu);
153: PetscPrintf(PETSC_COMM_WORLD, " LAPACKsyevx: compute %d to %d-th eigensolutions...\n", il, iu);
154: PetscMalloc1(m * n + 1, &evecs_array);
155: PetscMalloc1(7 * n + 1, &rwork);
156: PetscMalloc1(5 * n + 1, &iwork);
157: PetscMalloc1(n + 1, &ifail);
159: /* in the case "I", vl and vu are not referenced */
160: vl = 0.0;
161: vu = 8.0;
162: 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: PetscFree(iwork);
165: PetscFree(ifail);
166: PetscFree(rwork);
167: }
168: if (TestZHEGV) {
169: PetscPrintf(PETSC_COMM_WORLD, " LAPACKsygv: compute all %" PetscInt_FMT " eigensolutions...\n", m);
170: PetscMalloc1(3 * n + 1, &rwork);
171: 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: MatDenseRestoreArray(B, &arrayB);
178: PetscFree(rwork);
179: }
180: if (TestZHEGVX) {
181: il = 1;
182: PetscBLASIntCast((0.2 * m), &iu);
183: PetscPrintf(PETSC_COMM_WORLD, " LAPACKsygv: compute %d to %d-th eigensolutions...\n", il, iu);
184: PetscMalloc1(m * n + 1, &evecs_array);
185: PetscMalloc1(6 * n + 1, &iwork);
186: ifail = iwork + 5 * n;
187: PetscMalloc1(7 * n + 1, &rwork);
188: MatDenseGetArray(B, &arrayB);
189: vl = 0.0;
190: vu = 8.0;
191: 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: MatDenseRestoreArray(B, &arrayB);
194: PetscFree(iwork);
195: PetscFree(rwork);
196: }
197: MatDenseRestoreArray(A_dense, &arrayA);
200: /* View evals */
201: PetscOptionsHasName(NULL, NULL, "-eig_view", &flg);
202: if (flg) {
203: PetscPrintf(PETSC_COMM_WORLD, " %d evals: \n", nevs);
204: for (i = 0; i < nevs; i++) PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " %g\n", i + il, (double)evals[i]);
205: }
207: /* Check residuals and orthogonality */
208: PetscMalloc1(nevs + 1, &evecs);
209: for (i = 0; i < nevs; i++) {
210: VecCreate(PETSC_COMM_SELF, &evecs[i]);
211: VecSetSizes(evecs[i], PETSC_DECIDE, n);
212: VecSetFromOptions(evecs[i]);
213: VecPlaceArray(evecs[i], evecs_array + i * n);
214: }
216: tols[0] = PETSC_SQRT_MACHINE_EPSILON;
217: tols[1] = PETSC_SQRT_MACHINE_EPSILON;
218: CkEigenSolutions(cklvl, A, il - 1, iu - 1, evals, evecs, tols);
219: for (i = 0; i < nevs; i++) VecDestroy(&evecs[i]);
220: PetscFree(evecs);
222: /* Free work space. */
223: if (TestZHEEVX || TestZHEGVX) PetscFree(evecs_array);
224: PetscFree(evals);
225: PetscFree(work);
226: MatDestroy(&A_dense);
227: MatDestroy(&A);
228: MatDestroy(&B);
229: 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: nev = iu - il;
254: if (nev <= 0) return 0;
256: VecDuplicate(evec[0], &vt1);
257: VecDuplicate(evec[0], &vt2);
259: switch (cklvl) {
260: case 2:
261: dot_max = 0.0;
262: for (i = il; i < iu; i++) {
263: VecCopy(evec[i], vt1);
264: for (j = il; j < iu; j++) {
265: VecDot(evec[j], vt1, &dot);
266: if (j == i) {
267: rdot = PetscAbsScalar(dot - (PetscScalar)1.0);
268: } else {
269: rdot = PetscAbsScalar(dot);
270: }
271: if (rdot > dot_max) dot_max = rdot;
272: if (rdot > tols[1]) {
273: VecNorm(evec[i], NORM_INFINITY, &norm);
274: PetscPrintf(PETSC_COMM_SELF, "|delta(%" PetscInt_FMT ",%" PetscInt_FMT ")|: %g, norm: %g\n", i, j, (double)rdot, (double)norm);
275: }
276: }
277: }
278: PetscPrintf(PETSC_COMM_SELF, " max|(x_j^T*x_i) - delta_ji|: %g\n", (double)dot_max);
280: case 1:
281: norm_max = 0.0;
282: for (i = il; i < iu; i++) {
283: MatMult(A, evec[i], vt1);
284: VecCopy(evec[i], vt2);
285: tmp = -eval[i];
286: VecAXPY(vt1, tmp, vt2);
287: VecNorm(vt1, NORM_INFINITY, &norm);
288: norm = PetscAbs(norm);
289: if (norm > norm_max) norm_max = norm;
290: /* sniff, and bark if necessary */
291: if (norm > tols[0]) PetscPrintf(PETSC_COMM_WORLD, " residual violation: %" PetscInt_FMT ", resi: %g\n", i, (double)norm);
292: }
293: PetscPrintf(PETSC_COMM_SELF, " max_resi: %g\n", (double)norm_max);
294: break;
295: default:
296: PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%" PetscInt_FMT " is not supported \n", cklvl);
297: }
298: VecDestroy(&vt2);
299: VecDestroy(&vt1);
300: return 0;
301: }
303: /*TEST
305: build:
306: requires: complex
308: test:
310: test:
311: suffix: 2
312: args: -test_zheevx
314: test:
315: suffix: 3
316: args: -test_zhegv
318: test:
319: suffix: 4
320: args: -test_zhegvx
322: TEST*/