Actual source code: ex99.c
petsc-3.4.5 2014-06-29
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;
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(NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
51: if (!flg) {
52: PetscOptionsGetString(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(NULL,"-fB",file[1],PETSC_MAX_PATH_LEN,&flgB);
55: preload = PETSC_FALSE;
56: } else {
57: PetscOptionsGetString(NULL,"-fA",file[1],PETSC_MAX_PATH_LEN,&flgA);
58: if (!flgA) preload = PETSC_FALSE; /* don't bother with second system */
59: PetscOptionsGetString(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: }
88: /* Add a shift to A */
89: PetscOptionsGetScalar(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(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(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: PetscBLASIntCast(8*n,&lwork);
155: PetscBLASIntCast(n,&bn);
156: PetscMalloc(n*sizeof(PetscScalar),&evals);
157: PetscMalloc(lwork*sizeof(PetscScalar),&work);
158: MatDenseGetArray(A_dense,&arrayA);
159: MatDenseGetArray(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;
168: PetscBLASIntCast(.6*m,&iu);
169: PetscMalloc((m*n+1)*sizeof(PetscScalar),&evecs_array);
170: PetscMalloc((6*n+1)*sizeof(PetscBLASInt),&iwork);
171: ifail = iwork + 5*n;
172: if (PetscPreLoadIt) {PetscLogStagePush(stages[0]);}
173: /* in the case "I", vl and vu are not referenced */
174: 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);
175: if (PetscPreLoadIt) PetscLogStagePop();
176: PetscFree(iwork);
177: }
178: MatDenseRestoreArray(A,&arrayA);
179: MatDenseRestoreArray(B,&arrayB);
181: if (nevs <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);
182: /* View evals */
183: PetscOptionsHasName(NULL, "-eig_view", &flg);
184: if (flg) {
185: printf(" %d evals: \n",nevs);
186: for (i=0; i<nevs; i++) printf("%d %G\n",i+il,evals[i]);
187: }
189: /* Check residuals and orthogonality */
190: if (PetscPreLoadIt) {
191: mats[0] = A; mats[1] = B;
192: one = (PetscInt)one;
193: PetscMalloc((nevs+1)*sizeof(Vec),&evecs);
194: for (i=0; i<nevs; i++) {
195: VecCreate(PETSC_COMM_SELF,&evecs[i]);
196: VecSetSizes(evecs[i],PETSC_DECIDE,n);
197: VecSetFromOptions(evecs[i]);
198: VecPlaceArray(evecs[i],evecs_array+i*n);
199: }
201: ievbd_loc[0] = 0; ievbd_loc[1] = nevs-1;
202: tols[0] = 1.e-8; tols[1] = 1.e-8;
204: PetscLogStagePush(stages[1]);
205: CkEigenSolutions(&cklvl,mats,evals,evecs,ievbd_loc,&offset,tols);
206: PetscLogStagePop();
207: for (i=0; i<nevs; i++) { VecDestroy(&evecs[i]);}
208: PetscFree(evecs);
209: }
211: /* Free work space. */
212: if (TestSYGVX) {PetscFree(evecs_array);}
214: PetscFree(evals);
215: PetscFree(work);
217: MatDestroy(&A_dense);
218: MatDestroy(&B_dense);
219: MatDestroy(&B);
220: MatDestroy(&A);
222: PetscPreLoadEnd();
223: PetscFinalize();
224: return 0;
225: }
226: /*------------------------------------------------
227: Check the accuracy of the eigen solution
228: ----------------------------------------------- */
229: /*
230: input:
231: cklvl - check level:
232: 1: check residual
233: 2: 1 and check B-orthogonality locally
234: mats - matrix pencil
235: eval, evec - eigenvalues and eigenvectors stored in this process
236: ievbd_loc - local eigenvalue bounds, see eigc()
237: offset - see eigc()
238: tols[0] - reporting tol_res: || A evec[i] - eval[i] B evec[i]||
239: tols[1] - reporting tol_orth: evec[i] B evec[j] - delta_ij
240: */
241: #undef DEBUG_CkEigenSolutions
244: PetscErrorCode CkEigenSolutions(PetscInt *fcklvl,Mat *mats,PetscReal *eval,Vec *evec,PetscInt *ievbd_loc,PetscInt *offset,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: #if defined(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: #if defined(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: }
300: PetscPrintf(PETSC_COMM_SELF," max_resi: %G\n", norm_max);
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: }