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
11: PetscInt main(PetscInt argc,char **args)
12: {
13: Mat A,A_dense,B;
14: Vec *evecs;
15: PetscTruth flg,TestZHEEV=PETSC_TRUE,TestZHEEVX=PETSC_FALSE,TestZHEGV=PETSC_FALSE,TestZHEGVX=PETSC_FALSE;
17: PetscTruth isSymmetric;
18: PetscScalar sigma,*arrayA,*arrayB,*evecs_array=PETSC_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;
31:
32: PetscInitialize(&argc,&args,(char *)0,help);
33: #if !defined(PETSC_USE_COMPLEX)
34: SETERRQ(1,"This example requires complex numbers");
35: #endif
36: MPI_Comm_size(PETSC_COMM_WORLD,&size);
37: if (size != 1) SETERRQ(PETSC_ERR_SUP,"This is a uniprocessor example only!");
39: PetscOptionsHasName(PETSC_NULL, "-test_zheevx", &flg);
40: if (flg){
41: TestZHEEV = PETSC_FALSE;
42: TestZHEEVX = PETSC_TRUE;
43: }
44: PetscOptionsHasName(PETSC_NULL, "-test_zhegv", &flg);
45: if (flg){
46: TestZHEEV = PETSC_FALSE;
47: TestZHEGV= PETSC_TRUE;
48: }
49: PetscOptionsHasName(PETSC_NULL, "-test_zhegvx", &flg);
50: if (flg){
51: TestZHEEV = PETSC_FALSE;
52: TestZHEGVX = PETSC_TRUE;
53: }
55: PetscOptionsGetReal(PETSC_NULL,"-sigma1",&sigma1,PETSC_NULL);
56: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_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(PETSC_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: if (i<n-1) {
80: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
81: if (j>0) {
82: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
83: if (j<n-1) {
84: J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
85: if (use_random) {PetscRandomGetValue(rctx,&sigma2);}
86: v = 4.0 - sigma1*h2;
87: MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
88: }
89: /* make A complex Hermitian */
90: v = sigma2*h2;
91: Ii = 0; J = 1;
92: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
93: v = -sigma2*h2;
94: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
95: if (use_random) {PetscRandomDestroy(rctx);}
96: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
97: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
98: m = n = dim;
100: /* Check whether A is symmetric */
101: PetscOptionsHasName(PETSC_NULL, "-check_symmetry", &flg);
102: if (flg) {
103: Mat Trans;
104: MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
105: MatEqual(A, Trans, &isSymmetric);
106: if (!isSymmetric) SETERRQ(PETSC_ERR_USER,"A must be symmetric");
107: MatDestroy(Trans);
108: }
110: /* Convert aij matrix to MatSeqDense for LAPACK */
111: PetscTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
112: if (flg) {
113: MatDuplicate(A,MAT_COPY_VALUES,&A_dense);
114: } else {
115: MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);
116: }
118: MatCreate(PETSC_COMM_SELF,&B);
119: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
120: MatSetType(B,MATSEQDENSE);
121: MatSetFromOptions(B);
122: v = 1.0;
123: for (Ii=0; Ii<dim; Ii++) {
124: MatSetValues(B,1,&Ii,1,&Ii,&v,ADD_VALUES);
125: }
127: /* Solve standard eigenvalue problem: A*x = lambda*x */
128: /*===================================================*/
129: lwork = PetscBLASIntCast(2*n);
130: bn = PetscBLASIntCast(n);
131: PetscMalloc(n*sizeof(PetscReal),&evals);
132: PetscMalloc(lwork*sizeof(PetscScalar),&work);
133: MatGetArray(A_dense,&arrayA);
135: if (TestZHEEV){ /* test zheev() */
136: printf(" LAPACKsyev: compute all %d eigensolutions...\n",m);
137: PetscMalloc((3*n-2)*sizeof(PetscReal),&rwork);
138: LAPACKsyev_("V","U",&bn,arrayA,&bn,evals,work,&lwork,rwork,&lierr);
139: PetscFree(rwork);
140: evecs_array = arrayA;
141: nevs = m;
142: il=1; iu=m;
143: }
144: if (TestZHEEVX){
145: il = 1; iu=PetscBLASIntCast((0.2*m)); /* request 1 to 20%m evalues */
146: printf(" LAPACKsyevx: compute %d to %d-th eigensolutions...\n",il,iu);
147: PetscMalloc((m*n+1)*sizeof(PetscScalar),&evecs_array);
148: PetscMalloc((7*n+1)*sizeof(PetscReal),&rwork);
149: PetscMalloc((5*n+1)*sizeof(PetscBLASInt),&iwork);
150: PetscMalloc((n+1)*sizeof(PetscBLASInt),&ifail);
151:
152: /* in the case "I", vl and vu are not referenced */
153: vl = 0.0; vu = 8.0;
154: LAPACKsyevx_("V","I","U",&bn,arrayA,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&n,work,&lwork,rwork,iwork,ifail,&lierr);
155: PetscFree(iwork);
156: PetscFree(ifail);
157: PetscFree(rwork);
158: }
159: if (TestZHEGV){
160: printf(" LAPACKsygv: compute all %d eigensolutions...\n",m);
161: PetscMalloc((3*n+1)*sizeof(PetscReal),&rwork);
162: MatGetArray(B,&arrayB);
163: LAPACKsygv_(&one,"V","U",&bn,arrayA,&bn,arrayB,&bn,evals,work,&lwork,rwork,&lierr);
164: evecs_array = arrayA;
165: nevs = m;
166: il=1; iu=m;
167: MatRestoreArray(B,&arrayB);
168: PetscFree(rwork);
169: }
170: if (TestZHEGVX){
171: il = 1; iu=PetscBLASIntCast((0.2*m)); /* request 1 to 20%m evalues */
172: printf(" LAPACKsygv: compute %d to %d-th eigensolutions...\n",il,iu);
173: PetscMalloc((m*n+1)*sizeof(PetscScalar),&evecs_array);
174: PetscMalloc((6*n+1)*sizeof(PetscBLASInt),&iwork);
175: ifail = iwork + 5*n;
176: PetscMalloc((7*n+1)*sizeof(PetscReal),&rwork);
177: MatGetArray(B,&arrayB);
178: vl = 0.0; vu = 8.0;
179: 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);
180: MatRestoreArray(B,&arrayB);
181: PetscFree(iwork);
182: PetscFree(rwork);
183: }
184: MatRestoreArray(A_dense,&arrayA);
185: if (nevs <= 0 ) SETERRQ1(PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);
187: /* View evals */
188: PetscOptionsHasName(PETSC_NULL, "-eig_view", &flg);
189: if (flg){
190: printf(" %d evals: \n",nevs);
191: for (i=0; i<nevs; i++) printf("%d %G\n",i+il,evals[i]);
192: }
194: /* Check residuals and orthogonality */
195: PetscMalloc((nevs+1)*sizeof(Vec),&evecs);
196: for (i=0; i<nevs; i++){
197: VecCreate(PETSC_COMM_SELF,&evecs[i]);
198: VecSetSizes(evecs[i],PETSC_DECIDE,n);
199: VecSetFromOptions(evecs[i]);
200: VecPlaceArray(evecs[i],evecs_array+i*n);
201: }
202:
203: tols[0] = 1.e-8; tols[1] = 1.e-8;
204: CkEigenSolutions(cklvl,A,il-1,iu-1,evals,evecs,tols);
205: for (i=0; i<nevs; i++){ VecDestroy(evecs[i]);}
206: PetscFree(evecs);
207:
208: /* Free work space. */
209: if (TestZHEEVX || TestZHEGVX){
210: PetscFree(evecs_array);
211: }
212: PetscFree(evals);
213: PetscFree(work);
214: MatDestroy(A_dense);
215: MatDestroy(A);
216: MatDestroy(B);
217: PetscFinalize();
218: return 0;
219: }
220: /*------------------------------------------------
221: Check the accuracy of the eigen solution
222: ----------------------------------------------- */
223: /*
224: input:
225: cklvl - check level:
226: 1: check residual
227: 2: 1 and check B-orthogonality locally
228: A - matrix
229: il,iu - lower and upper index bound of eigenvalues
230: eval, evec - eigenvalues and eigenvectors stored in this process
231: tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
232: tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
233: */
234: #undef DEBUG_CkEigenSolutions
237: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal *eval,Vec *evec,PetscReal *tols)
238: {
239: PetscInt ierr,i,j,nev;
240: Vec vt1,vt2; /* tmp vectors */
241: PetscReal norm,tmp,norm_max,dot_max,rdot;
242: PetscScalar dot;
245: nev = iu - il;
246: if (nev <= 0) return(0);
248: //VecView(evec[0],PETSC_VIEWER_STDOUT_SELF);
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: //printf("ck %d-th\n",i);
257: VecCopy(evec[i], vt1);
258: for (j=il; j<iu; j++){
259: VecDot(evec[j],vt1,&dot);
260: if (j == i){
261: rdot = PetscAbsScalar(dot - 1.0);
262: } else {
263: rdot = PetscAbsScalar(dot);
264: }
265: if (rdot > dot_max) dot_max = rdot;
266: #ifdef DEBUG_CkEigenSolutions
267: if (rdot > tols[1] ) {
268: VecNorm(evec[i],NORM_INFINITY,&norm);
269: PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %G, norm: %G\n",i,j,dot,norm);
270: }
271: #endif
272: }
273: }
274: PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %G\n",dot_max);
276: case 1:
277: norm_max = 0.0;
278: for (i = il; i< iu; i++){
279: MatMult(A, evec[i], vt1);
280: VecCopy(evec[i], vt2);
281: tmp = -eval[i];
282: VecAXPY(vt1,tmp,vt2);
283: VecNorm(vt1, NORM_INFINITY, &norm);
284: norm = PetscAbsScalar(norm);
285: if (norm > norm_max) norm_max = norm;
286: #ifdef DEBUG_CkEigenSolutions
287: /* sniff, and bark if necessary */
288: if (norm > tols[0]){
289: printf( " residual violation: %d, resi: %g\n",i, norm);
290: }
291: #endif
292: }
293: PetscPrintf(PETSC_COMM_SELF," max_resi: %G\n", norm_max);
294: break;
295: default:
296: PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
297: }
298: VecDestroy(vt2);
299: VecDestroy(vt1);
300: return(0);
301: }