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) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A_dense));
118: else PetscCall(MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &A_dense));
120: PetscCall(MatCreate(PETSC_COMM_SELF, &B));
121: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, dim, dim));
122: PetscCall(MatSetType(B, MATSEQDENSE));
123: PetscCall(MatSetFromOptions(B));
124: PetscCall(MatSetUp(B));
125: v = 1.0;
126: for (Ii = 0; Ii < dim; Ii++) PetscCall(MatSetValues(B, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
128: /* Solve standard eigenvalue problem: A*x = lambda*x */
129: /*===================================================*/
130: PetscCall(PetscBLASIntCast(2 * n, &lwork));
131: PetscCall(PetscBLASIntCast(n, &bn));
132: PetscCall(PetscMalloc1(n, &evals));
133: PetscCall(PetscMalloc1(lwork, &work));
134: PetscCall(MatDenseGetArray(A_dense, &arrayA));
136: if (TestZHEEV) { /* test zheev() */
137: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsyev: compute all %" PetscInt_FMT " eigensolutions...\n", m));
138: PetscCall(PetscMalloc1(3 * n - 2, &rwork));
139: PetscCallBLAS("LAPACKsyev_", LAPACKsyev_("V", "U", &bn, arrayA, &bn, evals, work, &lwork, rwork, &lierr));
140: PetscCall(PetscFree(rwork));
142: evecs_array = arrayA;
143: PetscCall(PetscBLASIntCast(m, &nevs));
144: il = 1;
145: PetscCall(PetscBLASIntCast(m, &iu));
146: }
147: if (TestZHEEVX) {
148: il = 1;
149: PetscCall(PetscBLASIntCast(0.2 * m, &iu));
150: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsyevx: compute %d to %d-th eigensolutions...\n", il, iu));
151: PetscCall(PetscMalloc1(m * n + 1, &evecs_array));
152: PetscCall(PetscMalloc1(7 * n + 1, &rwork));
153: PetscCall(PetscMalloc1(5 * n + 1, &iwork));
154: PetscCall(PetscMalloc1(n + 1, &ifail));
156: /* in the case "I", vl and vu are not referenced */
157: vl = 0.0;
158: vu = 8.0;
159: PetscCall(PetscBLASIntCast(n, &nn));
160: PetscCallBLAS("LAPACKsyevx_", LAPACKsyevx_("V", "I", "U", &bn, arrayA, &bn, &vl, &vu, &il, &iu, &abstol, &nevs, evals, evecs_array, &nn, work, &lwork, rwork, iwork, ifail, &lierr));
161: PetscCall(PetscFree(iwork));
162: PetscCall(PetscFree(ifail));
163: PetscCall(PetscFree(rwork));
164: }
165: if (TestZHEGV) {
166: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsygv: compute all %" PetscInt_FMT " eigensolutions...\n", m));
167: PetscCall(PetscMalloc1(3 * n + 1, &rwork));
168: PetscCall(MatDenseGetArray(B, &arrayB));
169: PetscCallBLAS("LAPACKsygv_", LAPACKsygv_(&one, "V", "U", &bn, arrayA, &bn, arrayB, &bn, evals, work, &lwork, rwork, &lierr));
170: evecs_array = arrayA;
171: PetscCall(PetscBLASIntCast(m, &nevs));
172: il = 1;
173: PetscCall(PetscBLASIntCast(m, &iu));
174: PetscCall(MatDenseRestoreArray(B, &arrayB));
175: PetscCall(PetscFree(rwork));
176: }
177: if (TestZHEGVX) {
178: il = 1;
179: PetscCall(PetscBLASIntCast(0.2 * m, &iu));
180: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsygv: compute %d to %d-th eigensolutions...\n", il, iu));
181: PetscCall(PetscMalloc1(m * n + 1, &evecs_array));
182: PetscCall(PetscMalloc1(6 * n + 1, &iwork));
183: ifail = iwork + 5 * n;
184: PetscCall(PetscMalloc1(7 * n + 1, &rwork));
185: PetscCall(MatDenseGetArray(B, &arrayB));
186: vl = 0.0;
187: vu = 8.0;
188: PetscCall(PetscBLASIntCast(n, &nn));
189: PetscCallBLAS("LAPACKsygvx_", 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));
190: PetscCall(MatDenseRestoreArray(B, &arrayB));
191: PetscCall(PetscFree(iwork));
192: PetscCall(PetscFree(rwork));
193: }
194: PetscCall(MatDenseRestoreArray(A_dense, &arrayA));
195: PetscCheck(nevs > 0, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);
197: /* View evals */
198: PetscCall(PetscOptionsHasName(NULL, NULL, "-eig_view", &flg));
199: if (flg) {
200: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %d evals: \n", nevs));
201: for (i = 0; i < nevs; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " %g\n", i + il, (double)evals[i]));
202: }
204: /* Check residuals and orthogonality */
205: PetscCall(PetscMalloc1(nevs + 1, &evecs));
206: for (i = 0; i < nevs; i++) {
207: PetscCall(VecCreate(PETSC_COMM_SELF, &evecs[i]));
208: PetscCall(VecSetSizes(evecs[i], PETSC_DECIDE, n));
209: PetscCall(VecSetFromOptions(evecs[i]));
210: PetscCall(VecPlaceArray(evecs[i], evecs_array + i * n));
211: }
213: tols[0] = PETSC_SQRT_MACHINE_EPSILON;
214: tols[1] = PETSC_SQRT_MACHINE_EPSILON;
215: PetscCall(CkEigenSolutions(cklvl, A, il - 1, iu - 1, evals, evecs, tols));
216: for (i = 0; i < nevs; i++) PetscCall(VecDestroy(&evecs[i]));
217: PetscCall(PetscFree(evecs));
219: /* Free work space. */
220: if (TestZHEEVX || TestZHEGVX) PetscCall(PetscFree(evecs_array));
221: PetscCall(PetscFree(evals));
222: PetscCall(PetscFree(work));
223: PetscCall(MatDestroy(&A_dense));
224: PetscCall(MatDestroy(&A));
225: PetscCall(MatDestroy(&B));
226: PetscCall(PetscFinalize());
227: return 0;
228: }
229: /*------------------------------------------------
230: Check the accuracy of the eigen solution
231: ----------------------------------------------- */
232: /*
233: input:
234: cklvl - check level:
235: 1: check residual
236: 2: 1 and check B-orthogonality locally
237: A - matrix
238: il,iu - lower and upper index bound of eigenvalues
239: eval, evec - eigenvalues and eigenvectors stored in this process
240: tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
241: tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
242: */
243: PetscErrorCode CkEigenSolutions(PetscInt cklvl, Mat A, PetscInt il, PetscInt iu, PetscReal *eval, Vec *evec, PetscReal *tols)
244: {
245: PetscInt i, j, nev;
246: Vec vt1, vt2; /* tmp vectors */
247: PetscReal norm, tmp, norm_max, dot_max, rdot;
248: PetscScalar dot;
250: PetscFunctionBegin;
251: nev = iu - il;
252: if (nev <= 0) PetscFunctionReturn(PETSC_SUCCESS);
254: PetscCall(VecDuplicate(evec[0], &vt1));
255: PetscCall(VecDuplicate(evec[0], &vt2));
257: switch (cklvl) {
258: case 2:
259: dot_max = 0.0;
260: for (i = il; i < iu; i++) {
261: PetscCall(VecCopy(evec[i], vt1));
262: for (j = il; j < iu; j++) {
263: PetscCall(VecDot(evec[j], vt1, &dot));
264: if (j == i) {
265: rdot = PetscAbsScalar(dot - (PetscScalar)1.0);
266: } else {
267: rdot = PetscAbsScalar(dot);
268: }
269: if (rdot > dot_max) dot_max = rdot;
270: if (rdot > tols[1]) {
271: PetscCall(VecNorm(evec[i], NORM_INFINITY, &norm));
272: PetscCall(PetscPrintf(PETSC_COMM_SELF, "|delta(%" PetscInt_FMT ",%" PetscInt_FMT ")|: %g, norm: %g\n", i, j, (double)rdot, (double)norm));
273: }
274: }
275: }
276: PetscCall(PetscPrintf(PETSC_COMM_SELF, " max|(x_j^T*x_i) - delta_ji|: %g\n", (double)dot_max)); /* fall through */
277: case 1:
278: norm_max = 0.0;
279: for (i = il; i < iu; i++) {
280: PetscCall(MatMult(A, evec[i], vt1));
281: PetscCall(VecCopy(evec[i], vt2));
282: tmp = -eval[i];
283: PetscCall(VecAXPY(vt1, tmp, vt2));
284: PetscCall(VecNorm(vt1, NORM_INFINITY, &norm));
285: norm = PetscAbs(norm);
286: if (norm > norm_max) norm_max = norm;
287: /* sniff, and bark if necessary */
288: if (norm > tols[0]) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual violation: %" PetscInt_FMT ", resi: %g\n", i, (double)norm));
289: }
290: PetscCall(PetscPrintf(PETSC_COMM_SELF, " max_resi: %g\n", (double)norm_max));
291: break;
292: default:
293: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%" PetscInt_FMT " is not supported \n", cklvl));
294: }
295: PetscCall(VecDestroy(&vt2));
296: PetscCall(VecDestroy(&vt1));
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /*TEST
302: build:
303: requires: complex
305: test:
307: test:
308: suffix: 2
309: args: -test_zheevx
311: test:
312: suffix: 3
313: args: -test_zhegv
315: test:
316: suffix: 4
317: args: -test_zhegvx
319: TEST*/