Actual source code: ex120.c
petsc-3.6.4 2016-04-12
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*);
11: PetscInt main(PetscInt argc,char **args)
12: {
13: Mat A,A_dense,B;
14: Vec *evecs;
15: PetscBool flg,TestZHEEV=PETSC_TRUE,TestZHEEVX=PETSC_FALSE,TestZHEGV=PETSC_FALSE,TestZHEGVX=PETSC_FALSE;
17: PetscBool isSymmetric;
18: PetscScalar sigma,*arrayA,*arrayB,*evecs_array=NULL,*work;
19: PetscReal *evals,*rwork;
20: PetscMPIInt size;
21: PetscInt m,i,j,nevs,il,iu,cklvl=2;
22: PetscReal vl,vu,abstol=1.e-8;
23: PetscBLASInt *iwork,*ifail,lwork,lierr,bn;
24: PetscReal tols[2];
25: PetscInt nzeros[2],nz;
26: PetscReal ratio;
27: PetscScalar v,none = -1.0,sigma2,pfive = 0.5,*xa;
28: PetscRandom rctx;
29: PetscReal h2,sigma1 = 100.0;
30: PetscInt dim,Ii,J,Istart,Iend,n = 6,its,use_random,one=1;
32: PetscInitialize(&argc,&args,(char*)0,help);
33: #if !defined(PETSC_USE_COMPLEX)
34: SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex numbers");
35: #endif
36: MPI_Comm_size(PETSC_COMM_WORLD,&size);
37: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
39: PetscOptionsHasName(NULL, "-test_zheevx", &flg);
40: if (flg) {
41: TestZHEEV = PETSC_FALSE;
42: TestZHEEVX = PETSC_TRUE;
43: }
44: PetscOptionsHasName(NULL, "-test_zhegv", &flg);
45: if (flg) {
46: TestZHEEV = PETSC_FALSE;
47: TestZHEGV = PETSC_TRUE;
48: }
49: PetscOptionsHasName(NULL, "-test_zhegvx", &flg);
50: if (flg) {
51: TestZHEEV = PETSC_FALSE;
52: TestZHEGVX = PETSC_TRUE;
53: }
55: PetscOptionsGetReal(NULL,"-sigma1",&sigma1,NULL);
56: PetscOptionsGetInt(NULL,"-n",&n,NULL);
57: dim = n*n;
59: MatCreate(PETSC_COMM_SELF,&A);
60: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
61: MatSetType(A,MATSEQDENSE);
62: MatSetFromOptions(A);
64: PetscOptionsHasName(NULL,"-norandom",&flg);
65: if (flg) use_random = 0;
66: else use_random = 1;
67: if (use_random) {
68: PetscRandomCreate(PETSC_COMM_SELF,&rctx);
69: PetscRandomSetFromOptions(rctx);
70: PetscRandomSetInterval(rctx,0.0,PETSC_i);
71: } else {
72: sigma2 = 10.0*PETSC_i;
73: }
74: h2 = 1.0/((n+1)*(n+1));
75: for (Ii=0; Ii<dim; Ii++) {
76: v = -1.0; i = Ii/n; j = Ii - i*n;
77: if (i>0) {
78: J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
79: }
80: if (i<n-1) {
81: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
82: }
83: if (j>0) {
84: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
85: }
86: if (j<n-1) {
87: J = Ii+1; 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; J = 1;
96: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
97: v = -sigma2*h2;
98: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
99: if (use_random) {PetscRandomDestroy(&rctx);}
100: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
102: m = n = dim;
104: /* Check whether A is symmetric */
105: PetscOptionsHasName(NULL, "-check_symmetry", &flg);
106: if (flg) {
107: Mat Trans;
108: MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
109: MatEqual(A, Trans, &isSymmetric);
110: if (!isSymmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A must be symmetric");
111: MatDestroy(&Trans);
112: }
114: /* Convert aij matrix to MatSeqDense for LAPACK */
115: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
116: if (flg) {
117: MatDuplicate(A,MAT_COPY_VALUES,&A_dense);
118: } else {
119: MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);
120: }
122: MatCreate(PETSC_COMM_SELF,&B);
123: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
124: MatSetType(B,MATSEQDENSE);
125: MatSetFromOptions(B);
126: v = 1.0;
127: for (Ii=0; Ii<dim; Ii++) {
128: MatSetValues(B,1,&Ii,1,&Ii,&v,ADD_VALUES);
129: }
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: printf(" LAPACKsyev: compute all %d 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; iu=m;
148: }
149: if (TestZHEEVX) {
150: il = 1;
151: PetscBLASIntCast((0.2*m),&iu);
152: printf(" LAPACKsyevx: compute %d to %d-th eigensolutions...\n",il,iu);
153: PetscMalloc1(m*n+1,&evecs_array);
154: PetscMalloc1(7*n+1,&rwork);
155: PetscMalloc1(5*n+1,&iwork);
156: PetscMalloc1(n+1,&ifail);
158: /* in the case "I", vl and vu are not referenced */
159: vl = 0.0; vu = 8.0;
160: LAPACKsyevx_("V","I","U",&bn,arrayA,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&n,work,&lwork,rwork,iwork,ifail,&lierr);
161: PetscFree(iwork);
162: PetscFree(ifail);
163: PetscFree(rwork);
164: }
165: if (TestZHEGV) {
166: printf(" LAPACKsygv: compute all %d eigensolutions...\n",m);
167: PetscMalloc1(3*n+1,&rwork);
168: MatDenseGetArray(B,&arrayB);
169: LAPACKsygv_(&one,"V","U",&bn,arrayA,&bn,arrayB,&bn,evals,work,&lwork,rwork,&lierr);
170: evecs_array = arrayA;
171: nevs = m;
172: il = 1; iu=m;
173: MatDenseRestoreArray(B,&arrayB);
174: PetscFree(rwork);
175: }
176: if (TestZHEGVX) {
177: il = 1;
178: PetscBLASIntCast((0.2*m),&iu);
179: printf(" LAPACKsygv: compute %d to %d-th eigensolutions...\n",il,iu);
180: PetscMalloc1(m*n+1,&evecs_array);
181: PetscMalloc1(6*n+1,&iwork);
182: ifail = iwork + 5*n;
183: PetscMalloc1(7*n+1,&rwork);
184: MatDenseGetArray(B,&arrayB);
185: vl = 0.0; vu = 8.0;
186: LAPACKsygvx_(&one,"V","I","U",&bn,arrayA,&bn,arrayB,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&n,work,&lwork,rwork,iwork,ifail,&lierr);
187: MatDenseRestoreArray(B,&arrayB);
188: PetscFree(iwork);
189: PetscFree(rwork);
190: }
191: MatDenseRestoreArray(A_dense,&arrayA);
192: if (nevs <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);
194: /* View evals */
195: PetscOptionsHasName(NULL, "-eig_view", &flg);
196: if (flg) {
197: printf(" %d evals: \n",nevs);
198: for (i=0; i<nevs; i++) printf("%d %g\n",i+il,(double)evals[i]);
199: }
201: /* Check residuals and orthogonality */
202: PetscMalloc1(nevs+1,&evecs);
203: for (i=0; i<nevs; i++) {
204: VecCreate(PETSC_COMM_SELF,&evecs[i]);
205: VecSetSizes(evecs[i],PETSC_DECIDE,n);
206: VecSetFromOptions(evecs[i]);
207: VecPlaceArray(evecs[i],evecs_array+i*n);
208: }
210: tols[0] = 1.e-8; tols[1] = 1.e-8;
211: CkEigenSolutions(cklvl,A,il-1,iu-1,evals,evecs,tols);
212: for (i=0; i<nevs; i++) { VecDestroy(&evecs[i]);}
213: PetscFree(evecs);
215: /* Free work space. */
216: if (TestZHEEVX || TestZHEGVX) {
217: PetscFree(evecs_array);
218: }
219: PetscFree(evals);
220: PetscFree(work);
221: MatDestroy(&A_dense);
222: MatDestroy(&A);
223: MatDestroy(&B);
224: PetscFinalize();
225: return 0;
226: }
227: /*------------------------------------------------
228: Check the accuracy of the eigen solution
229: ----------------------------------------------- */
230: /*
231: input:
232: cklvl - check level:
233: 1: check residual
234: 2: 1 and check B-orthogonality locally
235: A - matrix
236: il,iu - lower and upper index bound of eigenvalues
237: eval, evec - eigenvalues and eigenvectors stored in this process
238: tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
239: tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
240: */
241: #undef DEBUG_CkEigenSolutions
244: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal *eval,Vec *evec,PetscReal *tols)
245: {
246: PetscInt ierr,i,j,nev;
247: Vec vt1,vt2; /* tmp vectors */
248: PetscReal norm,tmp,norm_max,dot_max,rdot;
249: PetscScalar dot;
252: nev = iu - il;
253: if (nev <= 0) return(0);
255: /*VecView(evec[0],PETSC_VIEWER_STDOUT_SELF); */
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: /*printf("ck %d-th\n",i); */
264: VecCopy(evec[i], vt1);
265: for (j=il; j<iu; j++) {
266: VecDot(evec[j],vt1,&dot);
267: if (j == i) {
268: rdot = PetscAbsScalar(dot - 1.0);
269: } else {
270: rdot = PetscAbsScalar(dot);
271: }
272: if (rdot > dot_max) dot_max = rdot;
273: #if defined(DEBUG_CkEigenSolutions)
274: if (rdot > tols[1]) {
275: VecNorm(evec[i],NORM_INFINITY,&norm);
276: PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %g, norm: %d\n",i,j,(double)dot,(double)norm);
277: }
278: #endif
279: }
280: }
281: PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max);
283: case 1:
284: norm_max = 0.0;
285: for (i = il; i< iu; i++) {
286: MatMult(A, evec[i], vt1);
287: VecCopy(evec[i], vt2);
288: tmp = -eval[i];
289: VecAXPY(vt1,tmp,vt2);
290: VecNorm(vt1, NORM_INFINITY, &norm);
291: norm = PetscAbsScalar(norm);
292: if (norm > norm_max) norm_max = norm;
293: #if defined(DEBUG_CkEigenSolutions)
294: /* sniff, and bark if necessary */
295: if (norm > tols[0]) {
296: printf(" residual violation: %d, resi: %g\n",i, norm);
297: }
298: #endif
299: }
300: PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max);
301: break;
302: default:
303: PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
304: }
305: VecDestroy(&vt2);
306: VecDestroy(&vt1);
307: return(0);
308: }