Actual source code: ex116.c
petsc-3.3-p6 2013-02-11
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,flgA=PETSC_FALSE,flgB=PETSC_FALSE,TestSYEVX=PETSC_TRUE;
25: PetscBool isSymmetric;
26: PetscScalar sigma,*arrayA,*arrayB,*evecs_array,*work,*evals;
27: PetscMPIInt size;
28: PetscInt m,n,i,j,nevs,il,iu,cklvl=2;
29: PetscReal vl,vu,abstol=1.e-8;
30: PetscBLASInt *iwork,*ifail,lwork,lierr,bn;
31: PetscReal tols[2];
32: PetscInt nzeros[2],nz;
33: PetscReal ratio;
34:
35: PetscInitialize(&argc,&args,(char *)0,help);
36: MPI_Comm_size(PETSC_COMM_WORLD,&size);
37: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
39: PetscOptionsHasName(PETSC_NULL, "-test_syev", &flg);
40: if (flg){
41: TestSYEVX = PETSC_FALSE;
42: }
44: /* Determine files from which we read the two matrices */
45: PetscOptionsGetString(PETSC_NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
47: /* Load matrix A */
48: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);
49: MatCreate(PETSC_COMM_WORLD,&A);
50: MatSetType(A,MATSEQAIJ);
51: MatLoad(A,fd);
52: PetscViewerDestroy(&fd);
53: MatGetSize(A,&m,&n);
55: /* Check whether A is symmetric */
56: PetscOptionsHasName(PETSC_NULL, "-check_symmetry", &flg);
57: if (flg) {
58: Mat Trans;
59: MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
60: MatEqual(A, Trans, &isSymmetric);
61: if (!isSymmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A must be symmetric");
62: MatDestroy(&Trans);
63: }
65: /* Convert aij matrix to MatSeqDense for LAPACK */
66: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
67: if (!flg) {
68: MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);
69: }
71: /* Solve eigenvalue problem: A*x = lambda*B*x */
72: /*============================================*/
73: lwork = PetscBLASIntCast(8*n);
74: bn = PetscBLASIntCast(n);
75: PetscMalloc(n*sizeof(PetscScalar),&evals);
76: PetscMalloc(lwork*sizeof(PetscScalar),&work);
77: MatGetArray(A_dense,&arrayA);
79: if (!TestSYEVX){ /* test syev() */
80: printf(" LAPACKsyev: compute all %d eigensolutions...\n",m);
81: LAPACKsyev_("V","U",&bn,arrayA,&bn,evals,work,&lwork,&lierr);
82: evecs_array = arrayA;
83: nevs = PetscBLASIntCast(m);
84: il=1; iu=PetscBLASIntCast(m);
85: } else { /* test syevx() */
86: il = 1; iu=PetscBLASIntCast((0.2*m)); /* request 1 to 20%m evalues */
87: printf(" LAPACKsyevx: compute %d to %d-th eigensolutions...\n",il,iu);
88: PetscMalloc((m*n+1)*sizeof(PetscScalar),&evecs_array);
89: PetscMalloc((6*n+1)*sizeof(PetscBLASInt),&iwork);
90: ifail = iwork + 5*n;
91:
92: /* in the case "I", vl and vu are not referenced */
93: vl = 0.0; vu = 8.0;
94: LAPACKsyevx_("V","I","U",&bn,arrayA,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&n,work,&lwork,iwork,ifail,&lierr);
95: PetscFree(iwork);
96: }
97: MatRestoreArray(A,&arrayA);
98: if (nevs <= 0 ) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);
100: /* View evals */
101: PetscOptionsHasName(PETSC_NULL, "-eig_view", &flg);
102: if (flg){
103: printf(" %d evals: \n",nevs);
104: for (i=0; i<nevs; i++) printf("%d %G\n",i+il,evals[i]);
105: }
107: /* Check residuals and orthogonality */
108: PetscMalloc((nevs+1)*sizeof(Vec),&evecs);
109: for (i=0; i<nevs; i++){
110: VecCreate(PETSC_COMM_SELF,&evecs[i]);
111: VecSetSizes(evecs[i],PETSC_DECIDE,n);
112: VecSetFromOptions(evecs[i]);
113: VecPlaceArray(evecs[i],evecs_array+i*n);
114: }
115:
116: tols[0] = 1.e-8; tols[1] = 1.e-8;
117: CkEigenSolutions(cklvl,A,il-1,iu-1,evals,evecs,tols);
118: for (i=0; i<nevs; i++){ VecDestroy(&evecs[i]);}
119: PetscFree(evecs);
120:
121: /* Free work space. */
122: if (TestSYEVX){PetscFree(evecs_array);}
123:
124: PetscFree(evals);
125: PetscFree(work);
127: MatDestroy(&A_dense);
128: MatDestroy(&A);
129: PetscFinalize();
130: return 0;
131: }
132: /*------------------------------------------------
133: Check the accuracy of the eigen solution
134: ----------------------------------------------- */
135: /*
136: input:
137: cklvl - check level:
138: 1: check residual
139: 2: 1 and check B-orthogonality locally
140: A - matrix
141: il,iu - lower and upper index bound of eigenvalues
142: eval, evec - eigenvalues and eigenvectors stored in this process
143: tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
144: tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
145: */
146: #undef DEBUG_CkEigenSolutions
149: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal *eval,Vec *evec,PetscReal *tols)
150: {
151: PetscInt ierr,i,j,nev;
152: Vec vt1,vt2; /* tmp vectors */
153: PetscReal norm,tmp,dot,norm_max,dot_max;
156: nev = iu - il;
157: if (nev <= 0) return(0);
159: //VecView(evec[0],PETSC_VIEWER_STDOUT_WORLD);
160: VecDuplicate(evec[0],&vt1);
161: VecDuplicate(evec[0],&vt2);
163: switch (cklvl){
164: case 2:
165: dot_max = 0.0;
166: for (i = il; i<iu; i++){
167: //printf("ck %d-th\n",i);
168: VecCopy(evec[i], vt1);
169: for (j=il; j<iu; j++){
170: VecDot(evec[j],vt1,&dot);
171: if (j == i){
172: dot = PetscAbsScalar(dot - 1.0);
173: } else {
174: dot = PetscAbsScalar(dot);
175: }
176: if (dot > dot_max) dot_max = dot;
177: #ifdef DEBUG_CkEigenSolutions
178: if (dot > tols[1] ) {
179: VecNorm(evec[i],NORM_INFINITY,&norm);
180: PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %G, norm: %G\n",i,j,dot,norm);
181: }
182: #endif
183: }
184: }
185: PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %G\n",dot_max);
187: case 1:
188: norm_max = 0.0;
189: for (i = il; i< iu; i++){
190: MatMult(A, evec[i], vt1);
191: VecCopy(evec[i], vt2);
192: tmp = -eval[i];
193: VecAXPY(vt1,tmp,vt2);
194: VecNorm(vt1, NORM_INFINITY, &norm);
195: norm = PetscAbsScalar(norm);
196: if (norm > norm_max) norm_max = norm;
197: #ifdef DEBUG_CkEigenSolutions
198: /* sniff, and bark if necessary */
199: if (norm > tols[0]){
200: printf( " residual violation: %d, resi: %g\n",i, norm);
201: }
202: #endif
203: }
204: PetscPrintf(PETSC_COMM_SELF," max_resi: %G\n", norm_max);
205: break;
206: default:
207: PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
208: }
209: VecDestroy(&vt2);
210: VecDestroy(&vt1);
211: return(0);
212: }