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";
3
4 #include <petscmat.h>
5 #include <petscblaslapack.h>
6
7 extern PetscErrorCode CkEigenSolutions(PetscInt, Mat, PetscInt, PetscInt, PetscReal *, Vec *, PetscReal *);
8
main(int argc,char ** args)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;
26
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!");
31
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 }
47
48 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
49 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
50 dim = n * n;
51
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));
57
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;
104
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 }
114
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 }
122
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));
130
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));
138
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));
144
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));
158
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);
199
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 }
206
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 }
215
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));
221
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 */
CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal * eval,Vec * evec,PetscReal * tols)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;
252
253 PetscFunctionBegin;
254 nev = iu - il;
255 if (nev <= 0) PetscFunctionReturn(PETSC_SUCCESS);
256
257 PetscCall(VecDuplicate(evec[0], &vt1));
258 PetscCall(VecDuplicate(evec[0], &vt2));
259
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));
280
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 }
303
304 /*TEST
305
306 build:
307 requires: complex
308
309 test:
310
311 test:
312 suffix: 2
313 args: -test_zheevx
314
315 test:
316 suffix: 3
317 args: -test_zhegv
318
319 test:
320 suffix: 4
321 args: -test_zhegvx
322
323 TEST*/
324