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