Actual source code: ex99.c
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 petscblaslapack.h
14: #include src/mat/impls/sbaij/seq/sbaij.h
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: PetscTruth flg,flgA=PETSC_FALSE,flgB=PETSC_FALSE,TestSYGVX=PETSC_TRUE;
28: PetscTruth 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,stages[2];
32: PetscReal vl,vu,abstol=1.e-8;
33: PetscBLASInt *iwork,*ifail,lone=1,lwork,lierr,bn;
34: PetscInt ievbd_loc[2],offset=0,cklvl=2;
35: PetscReal tols[2];
36: Mat_SeqSBAIJ *sbaij;
37: PetscScalar *aa;
38: PetscInt *ai,*aj;
39: PetscInt nzeros[2],nz;
40: PetscReal ratio;
41:
42: PetscInitialize(&argc,&args,(char *)0,help);
43: MPI_Comm_size(PETSC_COMM_WORLD,&size);
44: if (size != 1) SETERRQ(PETSC_ERR_SUP,"This is a uniprocessor example only!");
45: PetscLogStageRegister(&stages[0],"EigSolve");
46: PetscLogStageRegister(&stages[1],"EigCheck");
48: /* Determine files from which we read the two matrices */
49: PetscOptionsGetString(PETSC_NULL,"-f0",file[0],PETSC_MAX_PATH_LEN-1,&flg);
50: if (!flg) {
51: PetscOptionsGetString(PETSC_NULL,"-fA",file[0],PETSC_MAX_PATH_LEN-1,&flgA);
52: if (!flgA) SETERRQ(PETSC_ERR_USER,"Must indicate binary file with the -fA or -fB options");
53: PetscOptionsGetString(PETSC_NULL,"-fB",file[1],PETSC_MAX_PATH_LEN-1,&flgB);
54: preload = PETSC_FALSE;
55: } else {
56: PetscOptionsGetString(PETSC_NULL,"-fA",file[1],PETSC_MAX_PATH_LEN-1,&flgA);
57: if (!flgA) {preload = PETSC_FALSE;} /* don't bother with second system */
58: PetscOptionsGetString(PETSC_NULL,"-fB",file[2],PETSC_MAX_PATH_LEN-1,&flgB);
59: }
61: PreLoadBegin(preload,"Load system");
62: /* Load matrices */
63: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PreLoadIt],FILE_MODE_READ,&fd);
64: MatLoad(fd,MATSBAIJ,&A);
65: PetscViewerDestroy(fd);
66: MatGetSize(A,&m,&n);
67: if ((flgB && PreLoadIt) || (flgB && !preload)){
68: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PreLoadIt+1],FILE_MODE_READ,&fd);
69: MatLoad(fd,MATSBAIJ,&B);
70: PetscViewerDestroy(fd);
71: } else { /* create B=I */
72: MatCreate(PETSC_COMM_WORLD,&B);
73: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
74: MatSetType(B,MATSEQSBAIJ);
75: MatSetFromOptions(B);
76: for (i=0; i<m; i++) {
77: MatSetValues(B,1,&i,1,&i,&one,INSERT_VALUES);
78: }
79: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
81: }
82:
83: /* Add a shift to A */
84: PetscOptionsGetScalar(PETSC_NULL,"-mat_sigma",&sigma,&flg);
85: if(flg) {
86: MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN); /* A <- sigma*B + A */
87: }
89: /* Check whether A is symmetric */
90: PetscOptionsHasName(PETSC_NULL, "-check_symmetry", &flg);
91: if (flg) {
92: Mat Trans;
93: MatTranspose(A, &Trans);
94: MatEqual(A, Trans, &isSymmetric);
95: if (!isSymmetric) SETERRQ(PETSC_ERR_USER,"A must be symmetric");
96: MatDestroy(Trans);
97: if (flgB && PreLoadIt){
98: MatTranspose(B, &Trans);
99: MatEqual(B, Trans, &isSymmetric);
100: if (!isSymmetric) SETERRQ(PETSC_ERR_USER,"B must be symmetric");
101: MatDestroy(Trans);
102: }
103: }
105: /* View small entries of A */
106: PetscOptionsHasName(PETSC_NULL, "-Asp_view", &flg);
107: if (flg){
108: MatCreate(PETSC_COMM_SELF,&A_sp);
109: MatSetSizes(A_sp,PETSC_DECIDE,PETSC_DECIDE,m,n);
110: MatSetType(A_sp,MATSEQSBAIJ);
112: tols[0] = 1.e-6, tols[1] = 1.e-9;
113: sbaij = (Mat_SeqSBAIJ*)A->data;
114: ai = sbaij->i;
115: aj = sbaij->j;
116: aa = sbaij->a;
117: nzeros[0] = nzeros[1] = 0;
118: for (i=0; i<m; i++) {
119: nz = ai[i+1] - ai[i];
120: for (j=0; j<nz; j++){
121: if (PetscAbsScalar(*aa)<tols[0]) {
122: MatSetValues(A_sp,1,&i,1,aj,aa,INSERT_VALUES);
123: nzeros[0]++;
124: }
125: if (PetscAbsScalar(*aa)<tols[1]) nzeros[1]++;
126: aa++; aj++;
127: }
128: }
129: MatAssemblyBegin(A_sp,MAT_FINAL_ASSEMBLY);
130: MatAssemblyEnd(A_sp,MAT_FINAL_ASSEMBLY);
132: MatDestroy(A_sp);
134: ratio = (PetscReal)nzeros[0]/sbaij->nz;
135: PetscPrintf(PETSC_COMM_SELF," %d matrix entries < %e, ratio %G of %d nonzeros\n",nzeros[0],tols[0],ratio,sbaij->nz);
136: PetscPrintf(PETSC_COMM_SELF," %d matrix entries < %e\n",nzeros[1],tols[1]);
137: }
139: /* Convert aij matrix to MatSeqDense for LAPACK */
140: PetscTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
141: if (!flg) {
142: MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);
143: }
144: PetscTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
145: if (!flg) {MatConvert(B,MATSEQDENSE,MAT_INITIAL_MATRIX,&B_dense);}
147: /* Solve eigenvalue problem: A*x = lambda*B*x */
148: /*============================================*/
149: lwork = 8*n;
150: bn = (PetscBLASInt)n;
151: PetscMalloc(n*sizeof(PetscScalar),&evals);
152: PetscMalloc(lwork*sizeof(PetscScalar),&work);
153: MatGetArray(A_dense,&arrayA);
154: MatGetArray(B_dense,&arrayB);
156: if (!TestSYGVX){ /* test sygv() */
157: evecs_array = arrayA;
158: LAPACKsygv_(&lone,"V","U",&bn,arrayA,&bn,arrayB,&bn,evals,work,&lwork,&lierr);
159: nevs = m;
160: il=1;
161: } else { /* test sygvx() */
162: il = 1; iu=(PetscBLASInt)(.6*m); /* request 1 to 60%m evalues */
163: PetscMalloc((m*n+1)*sizeof(PetscScalar),&evecs_array);
164: PetscMalloc((6*n+1)*sizeof(PetscBLASInt),&iwork);
165: ifail = iwork + 5*n;
166: if(PreLoadIt){PetscLogStagePush(stages[0]);}
167: /* in the case "I", vl and vu are not referenced */
168: 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);
169: if(PreLoadIt){PetscLogStagePop();}
170: PetscFree(iwork);
171: }
172: MatRestoreArray(A,&arrayA);
173: MatRestoreArray(B,&arrayB);
175: if (nevs <= 0 ) SETERRQ1(PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);
176: /* View evals */
177: PetscOptionsHasName(PETSC_NULL, "-eig_view", &flg);
178: if (flg){
179: printf(" %d evals: \n",nevs);
180: for (i=0; i<nevs; i++) printf("%d %G\n",i+il,evals[i]);
181: }
183: /* Check residuals and orthogonality */
184: if(PreLoadIt){
185: mats[0] = A; mats[1] = B;
186: one = (PetscInt)one;
187: PetscMalloc((nevs+1)*sizeof(Vec),&evecs);
188: for (i=0; i<nevs; i++){
189: VecCreate(PETSC_COMM_SELF,&evecs[i]);
190: VecSetSizes(evecs[i],PETSC_DECIDE,n);
191: VecSetFromOptions(evecs[i]);
192: VecPlaceArray(evecs[i],evecs_array+i*n);
193: }
194:
195: ievbd_loc[0] = 0; ievbd_loc[1] = nevs-1;
196: tols[0] = 1.e-8; tols[1] = 1.e-8;
197: PetscLogStagePush(stages[1]);
198: CkEigenSolutions(&cklvl,mats,evals,evecs,ievbd_loc,&offset,tols);
199: PetscLogStagePop();
200: for (i=0; i<nevs; i++){ VecDestroy(evecs[i]);}
201: PetscFree(evecs);
202: }
203:
204: /* Free work space. */
205: if (TestSYGVX){PetscFree(evecs_array);}
206:
207: PetscFree(evals);
208: PetscFree(work);
210: MatDestroy(A_dense);
211: MatDestroy(B_dense);
212: MatDestroy(B);
213: MatDestroy(A);
215: PreLoadEnd();
216: PetscFinalize();
217: return 0;
218: }
219: /*------------------------------------------------
220: Check the accuracy of the eigen solution
221: ----------------------------------------------- */
222: /*
223: input:
224: cklvl - check level:
225: 1: check residual
226: 2: 1 and check B-orthogonality locally
227: mats - matrix pencil
228: eval, evec - eigenvalues and eigenvectors stored in this process
229: ievbd_loc - local eigenvalue bounds, see eigc()
230: offset - see eigc()
231: tols[0] - reporting tol_res: || A evec[i] - eval[i] B evec[i]||
232: tols[1] - reporting tol_orth: evec[i] B evec[j] - delta_ij
233: */
234: #undef DEBUG_CkEigenSolutions
237: PetscErrorCode CkEigenSolutions(PetscInt *fcklvl,Mat *mats,
238: PetscReal *eval,Vec *evec,PetscInt *ievbd_loc,PetscInt *offset,
239: PetscReal *tols)
240: {
241: PetscInt ierr,cklvl=*fcklvl,nev_loc,i,j;
242: Mat A=mats[0], B=mats[1];
243: Vec vt1,vt2; /* tmp vectors */
244: PetscReal norm,tmp,dot,norm_max,dot_max;
247: nev_loc = ievbd_loc[1] - ievbd_loc[0];
248: if (nev_loc == 0) return(0);
250: nev_loc += (*offset);
251: VecDuplicate(evec[*offset],&vt1);
252: VecDuplicate(evec[*offset],&vt2);
254: switch (cklvl){
255: case 2:
256: dot_max = 0.0;
257: for (i = *offset; i<nev_loc; i++){
258: MatMult(B, evec[i], vt1);
259: for (j=i; j<nev_loc; j++){
260: VecDot(evec[j],vt1,&dot);
261: if (j == i){
262: dot = PetscAbsScalar(dot - 1.0);
263: } else {
264: dot = PetscAbsScalar(dot);
265: }
266: if (dot > dot_max) dot_max = dot;
267: #ifdef DEBUG_CkEigenSolutions
268: if (dot > tols[1] ) {
269: VecNorm(evec[i],NORM_INFINITY,&norm);
270: PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %G, norm: %G\n",i,j,dot,norm);
271: }
272: #endif
273: } /* for (j=i; j<nev_loc; j++) */
274: }
275: PetscPrintf(PETSC_COMM_SELF," max|(x_j*B*x_i) - delta_ji|: %G\n",dot_max);
277: case 1:
278: norm_max = 0.0;
279: for (i = *offset; i< nev_loc; i++){
280: MatMult(A, evec[i], vt1);
281: MatMult(B, evec[i], vt2);
282: tmp = -eval[i];
283: VecAXPY(vt1,tmp,vt2);
284: VecNorm(vt1, NORM_INFINITY, &norm);
285: norm = PetscAbsScalar(norm);
286: if (norm > norm_max) norm_max = norm;
287: #ifdef DEBUG_CkEigenSolutions
288: /* sniff, and bark if necessary */
289: if (norm > tols[0]){
290: printf( " residual violation: %d, resi: %g\n",i, norm);
291: }
292: #endif
293: }
294:
295: PetscPrintf(PETSC_COMM_SELF," max_resi: %G\n", norm_max);
296:
297: break;
298: default:
299: PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
300: }
301: VecDestroy(vt2);
302: VecDestroy(vt1);
303: return(0);
304: }