Actual source code: ex116.c
petsc-3.4.2 2013-07-02
1: static char help[] = "Test LAPACK routine DSYEV() or DSYEVX(). \n\
2: Reads PETSc matrix A \n\
3: then computes selected eigenvalues, and optionally, eigenvectors of \n\
4: a real generalized symmetric-definite eigenproblem \n\
5: A*x = lambda*x \n\
6: Input parameters include\n\
7: -f <input_file> : file to load\n\
8: e.g. ./ex116 -f /home/petsc/datafiles/matrices/small \n\n";
10: #include <petscmat.h>
11: #include <petscblaslapack.h>
13: extern PetscErrorCode CkEigenSolutions(PetscInt,Mat,PetscInt,PetscInt,PetscReal*,Vec*,PetscReal*);
17: PetscInt main(PetscInt argc,char **args)
18: {
19: Mat A,A_dense;
20: Vec *evecs;
21: PetscViewer fd; /* viewer */
22: char file[1][PETSC_MAX_PATH_LEN]; /* input file name */
23: PetscBool flg,TestSYEVX=PETSC_TRUE;
25: PetscBool isSymmetric;
26: PetscScalar *arrayA,*evecs_array,*work,*evals;
27: PetscMPIInt size;
28: PetscInt m,n,i,nevs,il,iu,cklvl=2;
29: PetscReal vl,vu,abstol=1.e-8;
30: PetscBLASInt *iwork,*ifail,lwork,lierr,bn;
31: PetscReal tols[2];
33: PetscInitialize(&argc,&args,(char*)0,help);
34: MPI_Comm_size(PETSC_COMM_WORLD,&size);
35: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
37: PetscOptionsHasName(NULL, "-test_syev", &flg);
38: if (flg) {
39: TestSYEVX = PETSC_FALSE;
40: }
42: /* Determine files from which we read the two matrices */
43: PetscOptionsGetString(NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
45: /* Load matrix A */
46: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);
47: MatCreate(PETSC_COMM_WORLD,&A);
48: MatSetType(A,MATSEQAIJ);
49: MatLoad(A,fd);
50: PetscViewerDestroy(&fd);
51: MatGetSize(A,&m,&n);
53: /* Check whether A is symmetric */
54: PetscOptionsHasName(NULL, "-check_symmetry", &flg);
55: if (flg) {
56: Mat Trans;
57: MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
58: MatEqual(A, Trans, &isSymmetric);
59: if (!isSymmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A must be symmetric");
60: MatDestroy(&Trans);
61: }
63: /* Solve eigenvalue problem: A_dense*x = lambda*B*x */
64: /*==================================================*/
65: /* Convert aij matrix to MatSeqDense for LAPACK */
66: MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);
68: PetscBLASIntCast(8*n,&lwork);
69: PetscBLASIntCast(n,&bn);
70: PetscMalloc(n*sizeof(PetscScalar),&evals);
71: PetscMalloc(lwork*sizeof(PetscScalar),&work);
72: MatDenseGetArray(A_dense,&arrayA);
74: if (!TestSYEVX) { /* test syev() */
75: printf(" LAPACKsyev: compute all %d eigensolutions...\n",m);
76: LAPACKsyev_("V","U",&bn,arrayA,&bn,evals,work,&lwork,&lierr);
77: evecs_array = arrayA;
78: PetscBLASIntCast(m,&nevs);
79: il = 1;
80: PetscBLASIntCast(m,&iu);
81: } else { /* test syevx() */
82: il = 1;
83: PetscBLASIntCast((0.2*m,&iu));
84: printf(" LAPACKsyevx: compute %d to %d-th eigensolutions...\n",il,iu);
85: PetscMalloc((m*n+1)*sizeof(PetscScalar),&evecs_array);
86: PetscMalloc((6*n+1)*sizeof(PetscBLASInt),&iwork);
87: ifail = iwork + 5*n;
89: /* in the case "I", vl and vu are not referenced */
90: vl = 0.0; vu = 8.0;
91: LAPACKsyevx_("V","I","U",&bn,arrayA,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&n,work,&lwork,iwork,ifail,&lierr);
92: PetscFree(iwork);
93: }
94: MatDenseRestoreArray(A_dense,&arrayA);
95: if (nevs <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);
97: /* View eigenvalues */
98: PetscOptionsHasName(NULL, "-eig_view", &flg);
99: if (flg) {
100: printf(" %d evals: \n",nevs);
101: for (i=0; i<nevs; i++) printf("%d %G\n",i+il,evals[i]);
102: }
104: /* Check residuals and orthogonality */
105: PetscMalloc((nevs+1)*sizeof(Vec),&evecs);
106: for (i=0; i<nevs; i++) {
107: VecCreate(PETSC_COMM_SELF,&evecs[i]);
108: VecSetSizes(evecs[i],PETSC_DECIDE,n);
109: VecSetFromOptions(evecs[i]);
110: VecPlaceArray(evecs[i],evecs_array+i*n);
111: }
113: tols[0] = 1.e-8; tols[1] = 1.e-8;
114: CkEigenSolutions(cklvl,A,il-1,iu-1,evals,evecs,tols);
116: /* Free work space. */
117: for (i=0; i<nevs; i++) { VecDestroy(&evecs[i]);}
118: PetscFree(evecs);
119: MatDestroy(&A_dense);
120: PetscFree(work);
121: if (TestSYEVX) {PetscFree(evecs_array);}
123: /* Compute SVD: A_dense = U*SIGMA*transpose(V),
124: JOBU=JOBV='S': the first min(m,n) columns of U and V are returned in the arrayU and arrayV; */
125: /*==============================================================================================*/
126: {
127: /* Convert aij matrix to MatSeqDense for LAPACK */
128: PetscScalar *arrayU,*arrayVT,*arrayErr,alpha=1.0,beta=-1.0;
129: Mat Err;
130: PetscBLASInt minMN,maxMN;
131: PetscInt j;
132: PetscReal norm;
134: MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);
136: minMN = PetscMin(m,n);
137: maxMN = PetscMax(m,n);
138: lwork = 5*minMN + maxMN;
139: PetscMalloc4(m*minMN,PetscScalar,&arrayU,m*minMN,PetscScalar,&arrayVT,m*minMN,PetscScalar,&arrayErr,lwork,PetscScalar,&work);
141: /* Create matrix Err for checking error */
142: MatCreate(PETSC_COMM_WORLD,&Err);
143: MatSetSizes(Err,PETSC_DECIDE,PETSC_DECIDE,m,minMN);
144: MatSetType(Err,MATSEQDENSE);
145: MatSeqDenseSetPreallocation(Err,(PetscScalar*)arrayErr);
147: /* Save A to arrayErr for checking accuracy later. arrayA will be destroyed by LAPACKgesvd_() */
148: MatDenseGetArray(A_dense,&arrayA);
149: PetscMemcpy(arrayErr,arrayA,sizeof(PetscScalar)*m*minMN);
151: /* Compute A = U*SIGMA*VT */
152: LAPACKgesvd_("S","S",&m,&n,arrayA,&m,evals,arrayU,&minMN,arrayVT,&minMN,work,&lwork,&lierr);
153: MatDenseRestoreArray(A_dense,&arrayA);
154: if (!lierr) {
155: printf(" 1st 10 of %d singular values: \n",minMN);
156: for (i=0; i<10; i++) printf("%d %G\n",i,evals[i]);
157: } else {
158: printf("LAPACKgesvd_ fails!");
159: }
161: /* Check Err = (U*Sigma*V^T - A) using BLASgemm() */
162: /* U = U*Sigma */
163: for (j=0; j<minMN; j++) { /* U[:,j] = sigma[j]*U[:,j] */
164: for (i=0; i<m; i++) arrayU[j*m+i] *= evals[j];
165: }
166: /* Err = U*VT - A = alpha*U*VT + beta*Err */
167: BLASgemm_("N","N",&m,&minMN,&minMN,&alpha,arrayU,&m,arrayVT,&minMN,&beta,arrayErr,&m);
168: /*printf(" Err:\n");*/
169: /*MatView(Err,PETSC_VIEWER_STDOUT_SELF);*/
170: MatNorm(Err,NORM_FROBENIUS,&norm);
171: printf(" || U*Sigma*VT - A || = %G\n",norm);
173: PetscFree4(arrayU,arrayVT,arrayErr,work);
174: PetscFree(evals);
175: MatDestroy(&A_dense);
176: MatDestroy(&Err);
177: }
179: MatDestroy(&A);
180: PetscFinalize();
181: return 0;
182: }
183: /*------------------------------------------------
184: Check the accuracy of the eigen solution
185: ----------------------------------------------- */
186: /*
187: input:
188: cklvl - check level:
189: 1: check residual
190: 2: 1 and check B-orthogonality locally
191: A - matrix
192: il,iu - lower and upper index bound of eigenvalues
193: eval, evec - eigenvalues and eigenvectors stored in this process
194: tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
195: tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
196: */
197: #undef DEBUG_CkEigenSolutions
200: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal *eval,Vec *evec,PetscReal *tols)
201: {
202: PetscInt ierr,i,j,nev;
203: Vec vt1,vt2; /* tmp vectors */
204: PetscReal norm,tmp,dot,norm_max,dot_max;
207: nev = iu - il;
208: if (nev <= 0) return(0);
210: /*VecView(evec[0],PETSC_VIEWER_STDOUT_WORLD);*/
211: VecDuplicate(evec[0],&vt1);
212: VecDuplicate(evec[0],&vt2);
214: switch (cklvl) {
215: case 2:
216: dot_max = 0.0;
217: for (i = il; i<iu; i++) {
218: /*printf("ck %d-th\n",i);*/
219: VecCopy(evec[i], vt1);
220: for (j=il; j<iu; j++) {
221: VecDot(evec[j],vt1,&dot);
222: if (j == i) {
223: dot = PetscAbsScalar(dot - 1.0);
224: } else {
225: dot = PetscAbsScalar(dot);
226: }
227: if (dot > dot_max) dot_max = dot;
228: #if defined(DEBUG_CkEigenSolutions)
229: if (dot > tols[1]) {
230: VecNorm(evec[i],NORM_INFINITY,&norm);
231: PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %G, norm: %G\n",i,j,dot,norm);
232: }
233: #endif
234: }
235: }
236: PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %G\n",dot_max);
238: case 1:
239: norm_max = 0.0;
240: for (i = il; i< iu; i++) {
241: MatMult(A, evec[i], vt1);
242: VecCopy(evec[i], vt2);
243: tmp = -eval[i];
244: VecAXPY(vt1,tmp,vt2);
245: VecNorm(vt1, NORM_INFINITY, &norm);
246: norm = PetscAbsScalar(norm);
247: if (norm > norm_max) norm_max = norm;
248: #if defined(DEBUG_CkEigenSolutions)
249: /* sniff, and bark if necessary */
250: if (norm > tols[0]) {
251: printf(" residual violation: %d, resi: %g\n",i, norm);
252: }
253: #endif
254: }
255: PetscPrintf(PETSC_COMM_SELF," max_resi: %G\n", norm_max);
256: break;
257: default:
258: PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
259: }
260: VecDestroy(&vt2);
261: VecDestroy(&vt1);
262: return(0);
263: }