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*/