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: }