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: PetscInitialize(&argc,&args,(char*)0,help);
28: MPI_Comm_size(PETSC_COMM_WORLD,&size);
31: PetscOptionsHasName(NULL,NULL, "-test_zheevx", &flg);
32: if (flg) {
33: TestZHEEV = PETSC_FALSE;
34: TestZHEEVX = PETSC_TRUE;
35: }
36: PetscOptionsHasName(NULL,NULL, "-test_zhegv", &flg);
37: if (flg) {
38: TestZHEEV = PETSC_FALSE;
39: TestZHEGV = PETSC_TRUE;
40: }
41: PetscOptionsHasName(NULL,NULL, "-test_zhegvx", &flg);
42: if (flg) {
43: TestZHEEV = PETSC_FALSE;
44: TestZHEGVX = PETSC_TRUE;
45: }
47: PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);
48: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
49: dim = n*n;
51: MatCreate(PETSC_COMM_SELF,&A);
52: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
53: MatSetType(A,MATSEQDENSE);
54: MatSetFromOptions(A);
55: MatSetUp(A);
57: PetscOptionsHasName(NULL,NULL,"-norandom",&flg);
58: if (flg) use_random = 0;
59: else use_random = 1;
60: if (use_random) {
61: PetscRandomCreate(PETSC_COMM_SELF,&rctx);
62: PetscRandomSetFromOptions(rctx);
63: PetscRandomSetInterval(rctx,0.0,PETSC_i);
64: } else {
65: sigma2 = 10.0*PETSC_i;
66: }
67: h2 = 1.0/((n+1)*(n+1));
68: for (Ii=0; Ii<dim; Ii++) {
69: v = -1.0; i = Ii/n; j = Ii - i*n;
70: if (i>0) {
71: J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
72: }
73: if (i<n-1) {
74: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
75: }
76: if (j>0) {
77: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
78: }
79: if (j<n-1) {
80: J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
81: }
82: if (use_random) PetscRandomGetValue(rctx,&sigma2);
83: v = 4.0 - sigma1*h2;
84: MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
85: }
86: /* make A complex Hermitian */
87: v = sigma2*h2;
88: Ii = 0; J = 1;
89: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
90: v = -sigma2*h2;
91: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
92: if (use_random) PetscRandomDestroy(&rctx);
93: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
94: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
95: m = n = dim;
97: /* Check whether A is symmetric */
98: PetscOptionsHasName(NULL,NULL, "-check_symmetry", &flg);
99: if (flg) {
100: Mat Trans;
101: MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
102: MatEqual(A, Trans, &isSymmetric);
104: MatDestroy(&Trans);
105: }
107: /* Convert aij matrix to MatSeqDense for LAPACK */
108: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
109: if (flg) {
110: MatDuplicate(A,MAT_COPY_VALUES,&A_dense);
111: } else {
112: MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);
113: }
115: MatCreate(PETSC_COMM_SELF,&B);
116: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
117: MatSetType(B,MATSEQDENSE);
118: MatSetFromOptions(B);
119: MatSetUp(B);
120: v = 1.0;
121: for (Ii=0; Ii<dim; Ii++) {
122: MatSetValues(B,1,&Ii,1,&Ii,&v,ADD_VALUES);
123: }
125: /* Solve standard eigenvalue problem: A*x = lambda*x */
126: /*===================================================*/
127: PetscBLASIntCast(2*n,&lwork);
128: PetscBLASIntCast(n,&bn);
129: PetscMalloc1(n,&evals);
130: PetscMalloc1(lwork,&work);
131: MatDenseGetArray(A_dense,&arrayA);
133: if (TestZHEEV) { /* test zheev() */
134: PetscPrintf(PETSC_COMM_WORLD," LAPACKsyev: compute all %" PetscInt_FMT " eigensolutions...\n",m);
135: PetscMalloc1(3*n-2,&rwork);
136: LAPACKsyev_("V","U",&bn,arrayA,&bn,evals,work,&lwork,rwork,&lierr);
137: PetscFree(rwork);
139: evecs_array = arrayA;
140: nevs = m;
141: il =1; iu=m;
142: }
143: if (TestZHEEVX) {
144: il = 1;
145: PetscBLASIntCast((0.2*m),&iu);
146: PetscPrintf(PETSC_COMM_WORLD," LAPACKsyevx: compute %d to %d-th eigensolutions...\n",il,iu);
147: PetscMalloc1(m*n+1,&evecs_array);
148: PetscMalloc1(7*n+1,&rwork);
149: PetscMalloc1(5*n+1,&iwork);
150: PetscMalloc1(n+1,&ifail);
152: /* in the case "I", vl and vu are not referenced */
153: vl = 0.0; vu = 8.0;
154: PetscBLASIntCast(n,&nn);
155: LAPACKsyevx_("V","I","U",&bn,arrayA,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&nn,work,&lwork,rwork,iwork,ifail,&lierr);
156: PetscFree(iwork);
157: PetscFree(ifail);
158: PetscFree(rwork);
159: }
160: if (TestZHEGV) {
161: PetscPrintf(PETSC_COMM_WORLD," LAPACKsygv: compute all %" PetscInt_FMT " eigensolutions...\n",m);
162: PetscMalloc1(3*n+1,&rwork);
163: MatDenseGetArray(B,&arrayB);
164: LAPACKsygv_(&one,"V","U",&bn,arrayA,&bn,arrayB,&bn,evals,work,&lwork,rwork,&lierr);
165: evecs_array = arrayA;
166: nevs = m;
167: il = 1; iu=m;
168: MatDenseRestoreArray(B,&arrayB);
169: PetscFree(rwork);
170: }
171: if (TestZHEGVX) {
172: il = 1;
173: PetscBLASIntCast((0.2*m),&iu);
174: PetscPrintf(PETSC_COMM_WORLD," LAPACKsygv: compute %d to %d-th eigensolutions...\n",il,iu);
175: PetscMalloc1(m*n+1,&evecs_array);
176: PetscMalloc1(6*n+1,&iwork);
177: ifail = iwork + 5*n;
178: PetscMalloc1(7*n+1,&rwork);
179: MatDenseGetArray(B,&arrayB);
180: vl = 0.0; vu = 8.0;
181: PetscBLASIntCast(n,&nn);
182: 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);
183: MatDenseRestoreArray(B,&arrayB);
184: PetscFree(iwork);
185: PetscFree(rwork);
186: }
187: MatDenseRestoreArray(A_dense,&arrayA);
190: /* View evals */
191: PetscOptionsHasName(NULL,NULL, "-eig_view", &flg);
192: if (flg) {
193: PetscPrintf(PETSC_COMM_WORLD," %d evals: \n",nevs);
194: for (i=0; i<nevs; i++) PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT " %g\n",i+il,(double)evals[i]);
195: }
197: /* Check residuals and orthogonality */
198: PetscMalloc1(nevs+1,&evecs);
199: for (i=0; i<nevs; i++) {
200: VecCreate(PETSC_COMM_SELF,&evecs[i]);
201: VecSetSizes(evecs[i],PETSC_DECIDE,n);
202: VecSetFromOptions(evecs[i]);
203: VecPlaceArray(evecs[i],evecs_array+i*n);
204: }
206: tols[0] = PETSC_SQRT_MACHINE_EPSILON; tols[1] = PETSC_SQRT_MACHINE_EPSILON;
207: CkEigenSolutions(cklvl,A,il-1,iu-1,evals,evecs,tols);
208: for (i=0; i<nevs; i++) VecDestroy(&evecs[i]);
209: PetscFree(evecs);
211: /* Free work space. */
212: if (TestZHEEVX || TestZHEGVX) {
213: PetscFree(evecs_array);
214: }
215: PetscFree(evals);
216: PetscFree(work);
217: MatDestroy(&A_dense);
218: MatDestroy(&A);
219: MatDestroy(&B);
220: PetscFinalize();
221: return 0;
222: }
223: /*------------------------------------------------
224: Check the accuracy of the eigen solution
225: ----------------------------------------------- */
226: /*
227: input:
228: cklvl - check level:
229: 1: check residual
230: 2: 1 and check B-orthogonality locally
231: A - matrix
232: il,iu - lower and upper index bound of eigenvalues
233: eval, evec - eigenvalues and eigenvectors stored in this process
234: tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
235: tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
236: */
237: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal *eval,Vec *evec,PetscReal *tols)
238: {
239: PetscInt i,j,nev;
240: Vec vt1,vt2; /* tmp vectors */
241: PetscReal norm,tmp,norm_max,dot_max,rdot;
242: PetscScalar dot;
244: nev = iu - il;
245: if (nev <= 0) return 0;
247: VecDuplicate(evec[0],&vt1);
248: VecDuplicate(evec[0],&vt2);
250: switch (cklvl) {
251: case 2:
252: dot_max = 0.0;
253: for (i = il; i<iu; i++) {
254: VecCopy(evec[i], vt1);
255: for (j=il; j<iu; j++) {
256: VecDot(evec[j],vt1,&dot);
257: if (j == i) {
258: rdot = PetscAbsScalar(dot - (PetscScalar)1.0);
259: } else {
260: rdot = PetscAbsScalar(dot);
261: }
262: if (rdot > dot_max) dot_max = rdot;
263: if (rdot > tols[1]) {
264: VecNorm(evec[i],NORM_INFINITY,&norm);
265: PetscPrintf(PETSC_COMM_SELF,"|delta(%" PetscInt_FMT ",%" PetscInt_FMT ")|: %g, norm: %g\n",i,j,(double)rdot,(double)norm);
266: }
267: }
268: }
269: PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max);
271: case 1:
272: norm_max = 0.0;
273: for (i = il; i< iu; i++) {
274: MatMult(A, evec[i], vt1);
275: VecCopy(evec[i], vt2);
276: tmp = -eval[i];
277: VecAXPY(vt1,tmp,vt2);
278: VecNorm(vt1, NORM_INFINITY, &norm);
279: norm = PetscAbs(norm);
280: if (norm > norm_max) norm_max = norm;
281: /* sniff, and bark if necessary */
282: if (norm > tols[0]) {
283: PetscPrintf(PETSC_COMM_WORLD," residual violation: %" PetscInt_FMT ", resi: %g\n",i, norm);
284: }
285: }
286: PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max);
287: break;
288: default:
289: PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%" PetscInt_FMT " is not supported \n",cklvl);
290: }
291: VecDestroy(&vt2);
292: VecDestroy(&vt1);
293: return 0;
294: }
296: /*TEST
298: build:
299: requires: complex
301: test:
303: test:
304: suffix: 2
305: args: -test_zheevx
307: test:
308: suffix: 3
309: args: -test_zhegv
311: test:
312: suffix: 4
313: args: -test_zhegvx
315: TEST*/