Actual source code: ex120.c
petsc-3.10.5 2019-03-28
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;
15: PetscBool isSymmetric;
16: PetscScalar *arrayA,*arrayB,*evecs_array=NULL,*work;
17: PetscReal *evals,*rwork;
18: PetscMPIInt size;
19: PetscInt m,i,j,cklvl=2;
20: PetscReal vl,vu,abstol=1.e-8;
21: PetscBLASInt nn,nevs,il,iu,*iwork,*ifail,lwork,lierr,bn,one=1;
22: PetscReal tols[2];
23: PetscScalar v,sigma2;
24: PetscRandom rctx;
25: PetscReal h2,sigma1 = 100.0;
26: PetscInt dim,Ii,J,n = 6,use_random;
28: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
29: MPI_Comm_size(PETSC_COMM_WORLD,&size);
30: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
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; i = Ii/n; j = Ii - i*n;
71: if (i>0) {
72: J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
73: }
74: if (i<n-1) {
75: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
76: }
77: if (j>0) {
78: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
79: }
80: if (j<n-1) {
81: J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
82: }
83: if (use_random) {PetscRandomGetValue(rctx,&sigma2);}
84: v = 4.0 - sigma1*h2;
85: MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
86: }
87: /* make A complex Hermitian */
88: v = sigma2*h2;
89: Ii = 0; J = 1;
90: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
91: v = -sigma2*h2;
92: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
93: if (use_random) {PetscRandomDestroy(&rctx);}
94: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
95: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
96: m = n = dim;
98: /* Check whether A is symmetric */
99: PetscOptionsHasName(NULL,NULL, "-check_symmetry", &flg);
100: if (flg) {
101: Mat Trans;
102: MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
103: MatEqual(A, Trans, &isSymmetric);
104: if (!isSymmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A must be symmetric");
105: MatDestroy(&Trans);
106: }
108: /* Convert aij matrix to MatSeqDense for LAPACK */
109: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
110: if (flg) {
111: MatDuplicate(A,MAT_COPY_VALUES,&A_dense);
112: } else {
113: MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);
114: }
116: MatCreate(PETSC_COMM_SELF,&B);
117: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
118: MatSetType(B,MATSEQDENSE);
119: MatSetFromOptions(B);
120: MatSetUp(B);
121: v = 1.0;
122: for (Ii=0; Ii<dim; Ii++) {
123: MatSetValues(B,1,&Ii,1,&Ii,&v,ADD_VALUES);
124: }
126: /* Solve standard eigenvalue problem: A*x = lambda*x */
127: /*===================================================*/
128: PetscBLASIntCast(2*n,&lwork);
129: PetscBLASIntCast(n,&bn);
130: PetscMalloc1(n,&evals);
131: PetscMalloc1(lwork,&work);
132: MatDenseGetArray(A_dense,&arrayA);
134: if (TestZHEEV) { /* test zheev() */
135: PetscPrintf(PETSC_COMM_WORLD," LAPACKsyev: compute all %D eigensolutions...\n",m);
136: PetscMalloc1(3*n-2,&rwork);
137: LAPACKsyev_("V","U",&bn,arrayA,&bn,evals,work,&lwork,rwork,&lierr);
138: PetscFree(rwork);
140: evecs_array = arrayA;
141: nevs = m;
142: il =1; iu=m;
143: }
144: if (TestZHEEVX) {
145: il = 1;
146: PetscBLASIntCast((0.2*m),&iu);
147: PetscPrintf(PETSC_COMM_WORLD," LAPACKsyevx: compute %d to %d-th eigensolutions...\n",il,iu);
148: PetscMalloc1(m*n+1,&evecs_array);
149: PetscMalloc1(7*n+1,&rwork);
150: PetscMalloc1(5*n+1,&iwork);
151: PetscMalloc1(n+1,&ifail);
153: /* in the case "I", vl and vu are not referenced */
154: vl = 0.0; vu = 8.0;
155: PetscBLASIntCast(n,&nn);
156: LAPACKsyevx_("V","I","U",&bn,arrayA,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&nn,work,&lwork,rwork,iwork,ifail,&lierr);
157: PetscFree(iwork);
158: PetscFree(ifail);
159: PetscFree(rwork);
160: }
161: if (TestZHEGV) {
162: PetscPrintf(PETSC_COMM_WORLD," LAPACKsygv: compute all %D eigensolutions...\n",m);
163: PetscMalloc1(3*n+1,&rwork);
164: MatDenseGetArray(B,&arrayB);
165: LAPACKsygv_(&one,"V","U",&bn,arrayA,&bn,arrayB,&bn,evals,work,&lwork,rwork,&lierr);
166: evecs_array = arrayA;
167: nevs = m;
168: il = 1; iu=m;
169: MatDenseRestoreArray(B,&arrayB);
170: PetscFree(rwork);
171: }
172: if (TestZHEGVX) {
173: il = 1;
174: PetscBLASIntCast((0.2*m),&iu);
175: PetscPrintf(PETSC_COMM_WORLD," LAPACKsygv: compute %d to %d-th eigensolutions...\n",il,iu);
176: PetscMalloc1(m*n+1,&evecs_array);
177: PetscMalloc1(6*n+1,&iwork);
178: ifail = iwork + 5*n;
179: PetscMalloc1(7*n+1,&rwork);
180: MatDenseGetArray(B,&arrayB);
181: vl = 0.0; vu = 8.0;
182: PetscBLASIntCast(n,&nn);
183: 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);
184: MatDenseRestoreArray(B,&arrayB);
185: PetscFree(iwork);
186: PetscFree(rwork);
187: }
188: MatDenseRestoreArray(A_dense,&arrayA);
189: if (nevs <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);
191: /* View evals */
192: PetscOptionsHasName(NULL,NULL, "-eig_view", &flg);
193: if (flg) {
194: PetscPrintf(PETSC_COMM_WORLD," %d evals: \n",nevs);
195: for (i=0; i<nevs; i++) {PetscPrintf(PETSC_COMM_WORLD,"%D %g\n",i+il,(double)evals[i]);}
196: }
198: /* Check residuals and orthogonality */
199: PetscMalloc1(nevs+1,&evecs);
200: for (i=0; i<nevs; i++) {
201: VecCreate(PETSC_COMM_SELF,&evecs[i]);
202: VecSetSizes(evecs[i],PETSC_DECIDE,n);
203: VecSetFromOptions(evecs[i]);
204: VecPlaceArray(evecs[i],evecs_array+i*n);
205: }
207: tols[0] = PETSC_SQRT_MACHINE_EPSILON; tols[1] = PETSC_SQRT_MACHINE_EPSILON;
208: CkEigenSolutions(cklvl,A,il-1,iu-1,evals,evecs,tols);
209: for (i=0; i<nevs; i++) { VecDestroy(&evecs[i]);}
210: PetscFree(evecs);
212: /* Free work space. */
213: if (TestZHEEVX || TestZHEGVX) {
214: PetscFree(evecs_array);
215: }
216: PetscFree(evals);
217: PetscFree(work);
218: MatDestroy(&A_dense);
219: MatDestroy(&A);
220: MatDestroy(&B);
221: PetscFinalize();
222: return ierr;
223: }
224: /*------------------------------------------------
225: Check the accuracy of the eigen solution
226: ----------------------------------------------- */
227: /*
228: input:
229: cklvl - check level:
230: 1: check residual
231: 2: 1 and check B-orthogonality locally
232: A - matrix
233: il,iu - lower and upper index bound of eigenvalues
234: eval, evec - eigenvalues and eigenvectors stored in this process
235: tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
236: tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
237: */
238: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal *eval,Vec *evec,PetscReal *tols)
239: {
240: PetscInt ierr,i,j,nev;
241: Vec vt1,vt2; /* tmp vectors */
242: PetscReal norm,tmp,norm_max,dot_max,rdot;
243: PetscScalar dot;
246: nev = iu - il;
247: if (nev <= 0) return(0);
249: VecDuplicate(evec[0],&vt1);
250: VecDuplicate(evec[0],&vt2);
252: switch (cklvl) {
253: case 2:
254: dot_max = 0.0;
255: for (i = il; i<iu; i++) {
256: VecCopy(evec[i], vt1);
257: for (j=il; j<iu; j++) {
258: VecDot(evec[j],vt1,&dot);
259: if (j == i) {
260: rdot = PetscAbsScalar(dot - (PetscScalar)1.0);
261: } else {
262: rdot = PetscAbsScalar(dot);
263: }
264: if (rdot > dot_max) dot_max = rdot;
265: if (rdot > tols[1]) {
266: VecNorm(evec[i],NORM_INFINITY,&norm);
267: PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %g, norm: %d\n",i,j,(double)rdot,(double)norm);
268: }
269: }
270: }
271: PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max);
273: case 1:
274: norm_max = 0.0;
275: for (i = il; i< iu; i++) {
276: MatMult(A, evec[i], vt1);
277: VecCopy(evec[i], vt2);
278: tmp = -eval[i];
279: VecAXPY(vt1,tmp,vt2);
280: VecNorm(vt1, NORM_INFINITY, &norm);
281: norm = PetscAbs(norm);
282: if (norm > norm_max) norm_max = norm;
283: /* sniff, and bark if necessary */
284: if (norm > tols[0]) {
285: PetscPrintf(PETSC_COMM_WORLD," residual violation: %d, resi: %g\n",i, norm);
286: }
287: }
288: PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max);
289: break;
290: default:
291: PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
292: }
293: VecDestroy(&vt2);
294: VecDestroy(&vt1);
295: return(0);
296: }
299: /*TEST
301: build:
302: requires: complex
304: test:
306: test:
307: suffix: 2
308: args: -test_zheevx
310: test:
311: suffix: 3
312: args: -test_zhegv
314: test:
315: suffix: 4
316: args: -test_zhegvx
318: TEST*/