Actual source code: ex99.c

petsc-3.3-p6 2013-02-11
  1: static char help[] = "Test LAPACK routine DSYGV() or DSYGVX(). \n\
  2: Reads PETSc matrix A and B (or create B=I), \n\
  3: then computes selected eigenvalues, and optionally, eigenvectors of \n\
  4: a real generalized symmetric-definite eigenproblem \n\
  5:  A*x = lambda*B*x \n\
  6: Input parameters include\n\
  7:   -f0 <input_file> : first file to load (small system)\n\
  8:   -fA <input_file> -fB <input_file>: second files to load (larger system) \n\
  9: e.g. ./ex99 -f0 $D/small -fA $D/Eigdftb/dftb_bin/diamond_xxs_A -fB $D/Eigdftb/dftb_bin/diamond_xxs_B -mat_getrow_uppertriangular,\n\
 10:      where $D = /home/petsc/datafiles/matrices/Eigdftb/dftb_bin\n\n";

 12: #include <petscmat.h>
 13: #include <../src/mat/impls/sbaij/seq/sbaij.h>
 14: #include <petscblaslapack.h>

 16: extern PetscErrorCode CkEigenSolutions(PetscInt*,Mat*,PetscReal*,Vec*,PetscInt*,PetscInt*,PetscReal*);

 20: PetscInt main(PetscInt argc,char **args)
 21: {
 22:   Mat            A,B,A_dense,B_dense,mats[2],A_sp;
 23:   Vec            *evecs;
 24:   PetscViewer    fd;                /* viewer */
 25:   char           file[3][PETSC_MAX_PATH_LEN];     /* input file name */
 26:   PetscBool      flg,flgA=PETSC_FALSE,flgB=PETSC_FALSE,TestSYGVX=PETSC_TRUE;
 28:   PetscBool      preload=PETSC_TRUE,isSymmetric;
 29:   PetscScalar    sigma,one=1.0,*arrayA,*arrayB,*evecs_array,*work,*evals;
 30:   PetscMPIInt    size;
 31:   PetscInt       m,n,i,j,nevs,il,iu;
 32:   PetscLogStage  stages[2];
 33:   PetscReal      vl,vu,abstol=1.e-8;
 34:   PetscBLASInt   *iwork,*ifail,lone=1,lwork,lierr,bn;
 35:   PetscInt       ievbd_loc[2],offset=0,cklvl=2;
 36:   PetscReal      tols[2];
 37:   Mat_SeqSBAIJ   *sbaij;
 38:   PetscScalar    *aa;
 39:   PetscInt       *ai,*aj;
 40:   PetscInt       nzeros[2],nz;
 41:   PetscReal      ratio;
 42: 
 43:   PetscInitialize(&argc,&args,(char *)0,help);
 44:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 45:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
 46:   PetscLogStageRegister("EigSolve",&stages[0]);
 47:   PetscLogStageRegister("EigCheck",&stages[1]);

 49:   /* Determine files from which we read the two matrices */
 50:   PetscOptionsGetString(PETSC_NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
 51:   if (!flg) {
 52:     PetscOptionsGetString(PETSC_NULL,"-fA",file[0],PETSC_MAX_PATH_LEN,&flgA);
 53:     if (!flgA) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -fA or -fB options");
 54:     PetscOptionsGetString(PETSC_NULL,"-fB",file[1],PETSC_MAX_PATH_LEN,&flgB);
 55:     preload = PETSC_FALSE;
 56:   } else {
 57:     PetscOptionsGetString(PETSC_NULL,"-fA",file[1],PETSC_MAX_PATH_LEN,&flgA);
 58:     if (!flgA) {preload = PETSC_FALSE;} /* don't bother with second system */
 59:     PetscOptionsGetString(PETSC_NULL,"-fB",file[2],PETSC_MAX_PATH_LEN,&flgB);
 60:   }

 62:   PetscPreLoadBegin(preload,"Load system");
 63:     /* Load matrices */
 64:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);
 65:     MatCreate(PETSC_COMM_WORLD,&A);
 66:     MatSetType(A,MATSBAIJ);
 67:     MatLoad(A,fd);
 68:     PetscViewerDestroy(&fd);
 69:     MatGetSize(A,&m,&n);
 70:     if ((flgB && PetscPreLoadIt) || (flgB && !preload)){
 71:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt+1],FILE_MODE_READ,&fd);
 72:       MatCreate(PETSC_COMM_WORLD,&B);
 73:       MatSetType(B,MATSBAIJ);
 74:       MatLoad(B,fd);
 75:       PetscViewerDestroy(&fd);
 76:     } else { /* create B=I */
 77:       MatCreate(PETSC_COMM_WORLD,&B);
 78:       MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
 79:       MatSetType(B,MATSEQSBAIJ);
 80:       MatSetFromOptions(B);
 81:       for (i=0; i<m; i++) {
 82:         MatSetValues(B,1,&i,1,&i,&one,INSERT_VALUES);
 83:       }
 84:       MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 85:       MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 86:     }
 87: 
 88:     /* Add a shift to A */
 89:     PetscOptionsGetScalar(PETSC_NULL,"-mat_sigma",&sigma,&flg);
 90:     if(flg) {
 91:       MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN); /* A <- sigma*B + A */
 92:     }

 94:     /* Check whether A is symmetric */
 95:     PetscOptionsHasName(PETSC_NULL, "-check_symmetry", &flg);
 96:     if (flg) {
 97:       Mat Trans;
 98:       MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
 99:       MatEqual(A, Trans, &isSymmetric);
100:       if (!isSymmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A must be symmetric");
101:       MatDestroy(&Trans);
102:       if (flgB && PetscPreLoadIt){
103:         MatTranspose(B,MAT_INITIAL_MATRIX, &Trans);
104:         MatEqual(B, Trans, &isSymmetric);
105:         if (!isSymmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"B must be symmetric");
106:         MatDestroy(&Trans);
107:       }
108:     }

110:     /* View small entries of A */
111:     PetscOptionsHasName(PETSC_NULL, "-Asp_view", &flg);
112:     if (flg){
113:       MatCreate(PETSC_COMM_SELF,&A_sp);
114:       MatSetSizes(A_sp,PETSC_DECIDE,PETSC_DECIDE,m,n);
115:       MatSetType(A_sp,MATSEQSBAIJ);

117:       tols[0] = 1.e-6, tols[1] = 1.e-9;
118:       sbaij = (Mat_SeqSBAIJ*)A->data;
119:       ai    = sbaij->i;
120:       aj    = sbaij->j;
121:       aa    = sbaij->a;
122:       nzeros[0] = nzeros[1] = 0;
123:       for (i=0; i<m; i++) {
124:         nz = ai[i+1] - ai[i];
125:         for (j=0; j<nz; j++){
126:           if (PetscAbsScalar(*aa)<tols[0]) {
127:             MatSetValues(A_sp,1,&i,1,aj,aa,INSERT_VALUES);
128:             nzeros[0]++;
129:           }
130:           if (PetscAbsScalar(*aa)<tols[1]) nzeros[1]++;
131:           aa++; aj++;
132:         }
133:       }
134:       MatAssemblyBegin(A_sp,MAT_FINAL_ASSEMBLY);
135:       MatAssemblyEnd(A_sp,MAT_FINAL_ASSEMBLY);

137:       MatDestroy(&A_sp);

139:       ratio = (PetscReal)nzeros[0]/sbaij->nz;
140:       PetscPrintf(PETSC_COMM_SELF," %d matrix entries < %e, ratio %G of %d nonzeros\n",nzeros[0],tols[0],ratio,sbaij->nz);
141:       PetscPrintf(PETSC_COMM_SELF," %d matrix entries < %e\n",nzeros[1],tols[1]);
142:     }

144:     /* Convert aij matrix to MatSeqDense for LAPACK */
145:     PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
146:     if (!flg) {
147:       MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);
148:     }
149:     PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
150:     if (!flg) {MatConvert(B,MATSEQDENSE,MAT_INITIAL_MATRIX,&B_dense);}

152:     /* Solve eigenvalue problem: A*x = lambda*B*x */
153:     /*============================================*/
154:     lwork = PetscBLASIntCast(8*n);
155:     bn    = PetscBLASIntCast(n);
156:     PetscMalloc(n*sizeof(PetscScalar),&evals);
157:     PetscMalloc(lwork*sizeof(PetscScalar),&work);
158:     MatGetArray(A_dense,&arrayA);
159:     MatGetArray(B_dense,&arrayB);

161:     if (!TestSYGVX){ /* test sygv()  */
162:       evecs_array = arrayA;
163:       LAPACKsygv_(&lone,"V","U",&bn,arrayA,&bn,arrayB,&bn,evals,work,&lwork,&lierr);
164:       nevs = m;
165:       il=1;
166:     } else { /* test sygvx()  */
167:       il = 1; iu=PetscBLASIntCast(.6*m); /* request 1 to 60%m evalues */
168:       PetscMalloc((m*n+1)*sizeof(PetscScalar),&evecs_array);
169:       PetscMalloc((6*n+1)*sizeof(PetscBLASInt),&iwork);
170:       ifail = iwork + 5*n;
171:       if(PetscPreLoadIt){PetscLogStagePush(stages[0]);}
172:       /* in the case "I", vl and vu are not referenced */
173:       LAPACKsygvx_(&lone,"V","I","U",&bn,arrayA,&bn,arrayB,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&n,work,&lwork,iwork,ifail,&lierr);
174:       if(PetscPreLoadIt){PetscLogStagePop();}
175:       PetscFree(iwork);
176:     }
177:     MatRestoreArray(A,&arrayA);
178:     MatRestoreArray(B,&arrayB);

180:     if (nevs <= 0 ) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);
181:     /* View evals */
182:     PetscOptionsHasName(PETSC_NULL, "-eig_view", &flg);
183:     if (flg){
184:       printf(" %d evals: \n",nevs);
185:       for (i=0; i<nevs; i++) printf("%d  %G\n",i+il,evals[i]);
186:     }

188:     /* Check residuals and orthogonality */
189:     if(PetscPreLoadIt){
190:       mats[0] = A; mats[1] = B;
191:       one = (PetscInt)one;
192:       PetscMalloc((nevs+1)*sizeof(Vec),&evecs);
193:       for (i=0; i<nevs; i++){
194:         VecCreate(PETSC_COMM_SELF,&evecs[i]);
195:         VecSetSizes(evecs[i],PETSC_DECIDE,n);
196:         VecSetFromOptions(evecs[i]);
197:         VecPlaceArray(evecs[i],evecs_array+i*n);
198:       }
199: 
200:       ievbd_loc[0] = 0; ievbd_loc[1] = nevs-1;
201:       tols[0] = 1.e-8;  tols[1] = 1.e-8;
202:       PetscLogStagePush(stages[1]);
203:       CkEigenSolutions(&cklvl,mats,evals,evecs,ievbd_loc,&offset,tols);
204:       PetscLogStagePop();
205:       for (i=0; i<nevs; i++){ VecDestroy(&evecs[i]);}
206:       PetscFree(evecs);
207:     }
208: 
209:     /* Free work space. */
210:     if (TestSYGVX){PetscFree(evecs_array);}
211: 
212:     PetscFree(evals);
213:     PetscFree(work);

215:     MatDestroy(&A_dense);
216:     MatDestroy(&B_dense);
217:     MatDestroy(&B);
218:     MatDestroy(&A);

220:   PetscPreLoadEnd();
221:   PetscFinalize();
222:   return 0;
223: }
224: /*------------------------------------------------
225:   Check the accuracy of the eigen solution
226:   ----------------------------------------------- */
227: /*
228:   input: 
229:      cklvl      - check level: 
230:                     1: check residual
231:                     2: 1 and check B-orthogonality locally 
232:      mats       - matrix pencil
233:      eval, evec - eigenvalues and eigenvectors stored in this process
234:      ievbd_loc  - local eigenvalue bounds, see eigc()
235:      offset     - see eigc()
236:      tols[0]    - reporting tol_res: || A evec[i] - eval[i] B evec[i]||
237:      tols[1]    - reporting tol_orth: evec[i] B evec[j] - delta_ij
238: */
239: #undef DEBUG_CkEigenSolutions
242: PetscErrorCode CkEigenSolutions(PetscInt *fcklvl,Mat *mats,
243:                    PetscReal *eval,Vec *evec,PetscInt *ievbd_loc,PetscInt *offset, 
244:                    PetscReal *tols)
245: {
246:   PetscInt     ierr,cklvl=*fcklvl,nev_loc,i,j;
247:   Mat          A=mats[0], B=mats[1];
248:   Vec          vt1,vt2; /* tmp vectors */
249:   PetscReal    norm,tmp,dot,norm_max,dot_max;

252:   nev_loc = ievbd_loc[1] - ievbd_loc[0];
253:   if (nev_loc == 0) return(0);

255:   nev_loc += (*offset);
256:   VecDuplicate(evec[*offset],&vt1);
257:   VecDuplicate(evec[*offset],&vt2);

259:   switch (cklvl){
260:   case 2:
261:     dot_max = 0.0;
262:     for (i = *offset; i<nev_loc; i++){
263:       MatMult(B, evec[i], vt1);
264:       for (j=i; j<nev_loc; j++){
265:         VecDot(evec[j],vt1,&dot);
266:         if (j == i){
267:           dot = PetscAbsScalar(dot - 1.0);
268:         } else {
269:           dot = PetscAbsScalar(dot);
270:         }
271:         if (dot > dot_max) dot_max = dot;
272: #ifdef DEBUG_CkEigenSolutions
273:         if (dot > tols[1] ) {
274:           VecNorm(evec[i],NORM_INFINITY,&norm);
275:           PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %G, norm: %G\n",i,j,dot,norm);
276:         }
277: #endif
278:       } /* for (j=i; j<nev_loc; j++) */
279:     }
280:     PetscPrintf(PETSC_COMM_SELF,"    max|(x_j*B*x_i) - delta_ji|: %G\n",dot_max);

282:   case 1:
283:     norm_max = 0.0;
284:     for (i = *offset; i< nev_loc; i++){
285:       MatMult(A, evec[i], vt1);
286:       MatMult(B, evec[i], vt2);
287:       tmp  = -eval[i];
288:       VecAXPY(vt1,tmp,vt2);
289:       VecNorm(vt1, NORM_INFINITY, &norm);
290:       norm = PetscAbsScalar(norm);
291:       if (norm > norm_max) norm_max = norm;
292: #ifdef DEBUG_CkEigenSolutions
293:       /* sniff, and bark if necessary */
294:       if (norm > tols[0]){
295:         printf( "  residual violation: %d, resi: %g\n",i, norm);
296:       }
297: #endif
298:     }
299: 
300:       PetscPrintf(PETSC_COMM_SELF,"    max_resi:                    %G\n", norm_max);
301: 
302:    break;
303:   default:
304:     PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
305:   }
306:   VecDestroy(&vt2);
307:   VecDestroy(&vt1);
308:   return(0);
309: }