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: }