Actual source code: ex116.c

  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 $DATAFILESPATH/matrices/small  \n\n";

 10: #include <petscmat.h>
 11: #include <petscblaslapack.h>

 13: extern PetscErrorCode CkEigenSolutions(PetscInt,Mat,PetscInt,PetscInt,PetscReal*,Vec*,PetscReal*);

 15: int main(int argc,char **args)
 16: {
 17:   Mat            A,A_dense;
 18:   Vec            *evecs;
 19:   PetscViewer    fd;                /* viewer */
 20:   char           file[1][PETSC_MAX_PATH_LEN];     /* input file name */
 21:   PetscBool      flg,TestSYEVX=PETSC_TRUE;
 22:   PetscBool      isSymmetric;
 23:   PetscScalar    *arrayA,*evecs_array,*work,*evals;
 24:   PetscMPIInt    size;
 25:   PetscInt       m,n,i,cklvl=2;
 26:   PetscBLASInt   nevs,il,iu,in;
 27:   PetscReal      vl,vu,abstol=1.e-8;
 28:   PetscBLASInt   *iwork,*ifail,lwork,lierr,bn;
 29:   PetscReal      tols[2];

 31:   PetscInitialize(&argc,&args,(char*)0,help);
 32:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 35:   PetscOptionsHasName(NULL,NULL, "-test_syev", &flg);
 36:   if (flg) {
 37:     TestSYEVX = PETSC_FALSE;
 38:   }

 40:   /* Determine files from which we read the two matrices */
 41:   PetscOptionsGetString(NULL,NULL,"-f",file[0],sizeof(file[0]),&flg);

 43:   /* Load matrix A */
 44:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);
 45:   MatCreate(PETSC_COMM_WORLD,&A);
 46:   MatSetType(A,MATSEQAIJ);
 47:   MatLoad(A,fd);
 48:   PetscViewerDestroy(&fd);
 49:   MatGetSize(A,&m,&n);

 51:   /* Check whether A is symmetric */
 52:   PetscOptionsHasName(NULL,NULL, "-check_symmetry", &flg);
 53:   if (flg) {
 54:     Mat Trans;
 55:     MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
 56:     MatEqual(A, Trans, &isSymmetric);
 58:     MatDestroy(&Trans);
 59:   }

 61:   /* Solve eigenvalue problem: A_dense*x = lambda*B*x */
 62:   /*==================================================*/
 63:   /* Convert aij matrix to MatSeqDense for LAPACK */
 64:   MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);

 66:   PetscBLASIntCast(8*n,&lwork);
 67:   PetscBLASIntCast(n,&bn);
 68:   PetscMalloc1(n,&evals);
 69:   PetscMalloc1(lwork,&work);
 70:   MatDenseGetArray(A_dense,&arrayA);

 72:   if (!TestSYEVX) { /* test syev() */
 73:     PetscPrintf(PETSC_COMM_SELF," LAPACKsyev: compute all %" PetscInt_FMT " eigensolutions...\n",m);
 74:     LAPACKsyev_("V","U",&bn,arrayA,&bn,evals,work,&lwork,&lierr);
 75:     evecs_array = arrayA;
 76:     PetscBLASIntCast(m,&nevs);
 77:     il          = 1;
 78:     PetscBLASIntCast(m,&iu);
 79:   } else { /* test syevx()  */
 80:     il   = 1;
 81:     PetscBLASIntCast(0.2*m,&iu);
 82:     PetscBLASIntCast(n,&in);
 83:     PetscPrintf(PETSC_COMM_SELF," LAPACKsyevx: compute %" PetscBLASInt_FMT " to %" PetscBLASInt_FMT "-th eigensolutions...\n",il,iu);
 84:     PetscMalloc1(m*n+1,&evecs_array);
 85:     PetscMalloc1(6*n+1,&iwork);
 86:     ifail = iwork + 5*n;

 88:     /* in the case "I", vl and vu are not referenced */
 89:     vl = 0.0; vu = 8.0;
 90:     LAPACKsyevx_("V","I","U",&bn,arrayA,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&in,work,&lwork,iwork,ifail,&lierr);
 91:     PetscFree(iwork);
 92:   }
 93:   MatDenseRestoreArray(A_dense,&arrayA);

 96:   /* View eigenvalues */
 97:   PetscOptionsHasName(NULL,NULL, "-eig_view", &flg);
 98:   if (flg) {
 99:     PetscPrintf(PETSC_COMM_SELF," %" PetscBLASInt_FMT " evals: \n",nevs);
100:     for (i=0; i<nevs; i++) PetscPrintf(PETSC_COMM_SELF,"%" PetscInt_FMT "  %g\n",(PetscInt)(i+il),(double)evals[i]);
101:   }

103:   /* Check residuals and orthogonality */
104:   PetscMalloc1(nevs+1,&evecs);
105:   for (i=0; i<nevs; i++) {
106:     VecCreate(PETSC_COMM_SELF,&evecs[i]);
107:     VecSetSizes(evecs[i],PETSC_DECIDE,n);
108:     VecSetFromOptions(evecs[i]);
109:     VecPlaceArray(evecs[i],evecs_array+i*n);
110:   }

112:   tols[0] = tols[1] = PETSC_SQRT_MACHINE_EPSILON;
113:   CkEigenSolutions(cklvl,A,il-1,iu-1,evals,evecs,tols);

115:   /* Free work space. */
116:   for (i=0; i<nevs; i++) VecDestroy(&evecs[i]);
117:   PetscFree(evecs);
118:   MatDestroy(&A_dense);
119:   PetscFree(work);
120:   if (TestSYEVX) PetscFree(evecs_array);

122:   /* Compute SVD: A_dense = U*SIGMA*transpose(V),
123:      JOBU=JOBV='S':  the first min(m,n) columns of U and V are returned in the arrayU and arrayV; */
124:   /*==============================================================================================*/
125:   {
126:     /* Convert aij matrix to MatSeqDense for LAPACK */
127:     PetscScalar  *arrayU,*arrayVT,*arrayErr,alpha=1.0,beta=-1.0;
128:     Mat          Err;
129:     PetscBLASInt minMN,maxMN,im,in;
130:     PetscInt     j;
131:     PetscReal    norm;

133:     MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);

135:     minMN = PetscMin(m,n);
136:     maxMN = PetscMax(m,n);
137:     lwork = 5*minMN + maxMN;
138:     PetscMalloc4(m*minMN,&arrayU,m*minMN,&arrayVT,m*minMN,&arrayErr,lwork,&work);

140:     /* Create matrix Err for checking error */
141:     MatCreate(PETSC_COMM_WORLD,&Err);
142:     MatSetSizes(Err,PETSC_DECIDE,PETSC_DECIDE,m,minMN);
143:     MatSetType(Err,MATSEQDENSE);
144:     MatSeqDenseSetPreallocation(Err,(PetscScalar*)arrayErr);

146:     /* Save A to arrayErr for checking accuracy later. arrayA will be destroyed by LAPACKgesvd_() */
147:     MatDenseGetArray(A_dense,&arrayA);
148:     PetscArraycpy(arrayErr,arrayA,m*minMN);

150:     PetscBLASIntCast(m,&im);
151:     PetscBLASIntCast(n,&in);
152:     /* Compute A = U*SIGMA*VT */
153:     LAPACKgesvd_("S","S",&im,&in,arrayA,&im,evals,arrayU,&minMN,arrayVT,&minMN,work,&lwork,&lierr);
154:     MatDenseRestoreArray(A_dense,&arrayA);
155:     if (!lierr) {
156:       PetscPrintf(PETSC_COMM_SELF," 1st 10 of %" PetscBLASInt_FMT " singular values: \n",minMN);
157:       for (i=0; i<10; i++) PetscPrintf(PETSC_COMM_SELF,"%" PetscInt_FMT "  %g\n",i,(double)evals[i]);
158:     } else {
159:       PetscPrintf(PETSC_COMM_SELF,"LAPACKgesvd_ fails!");
160:     }

162:     /* Check Err = (U*Sigma*V^T - A) using BLASgemm() */
163:     /* U = U*Sigma */
164:     for (j=0; j<minMN; j++) { /* U[:,j] = sigma[j]*U[:,j] */
165:       for (i=0; i<m; i++) arrayU[j*m+i] *= evals[j];
166:     }
167:     /* Err = U*VT - A = alpha*U*VT + beta*Err */
168:     BLASgemm_("N","N",&im,&minMN,&minMN,&alpha,arrayU,&im,arrayVT,&minMN,&beta,arrayErr,&im);
169:     MatNorm(Err,NORM_FROBENIUS,&norm);
170:     PetscPrintf(PETSC_COMM_SELF," || U*Sigma*VT - A || = %g\n",(double)norm);

172:     PetscFree4(arrayU,arrayVT,arrayErr,work);
173:     PetscFree(evals);
174:     MatDestroy(&A_dense);
175:     MatDestroy(&Err);
176:   }

178:   MatDestroy(&A);
179:   PetscFinalize();
180:   return 0;
181: }
182: /*------------------------------------------------
183:   Check the accuracy of the eigen solution
184:   ----------------------------------------------- */
185: /*
186:   input:
187:      cklvl      - check level:
188:                     1: check residual
189:                     2: 1 and check B-orthogonality locally
190:      A          - matrix
191:      il,iu      - lower and upper index bound of eigenvalues
192:      eval, evec - eigenvalues and eigenvectors stored in this process
193:      tols[0]    - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
194:      tols[1]    - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
195: */
196: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal *eval,Vec *evec,PetscReal *tols)
197: {
198:   PetscInt  i,j,nev;
199:   Vec       vt1,vt2;    /* tmp vectors */
200:   PetscReal norm,tmp,dot,norm_max,dot_max;

202:   nev = iu - il;
203:   if (nev <= 0) return 0;

205:   /*VecView(evec[0],PETSC_VIEWER_STDOUT_WORLD);*/
206:   VecDuplicate(evec[0],&vt1);
207:   VecDuplicate(evec[0],&vt2);

209:   switch (cklvl) {
210:   case 2:
211:     dot_max = 0.0;
212:     for (i = il; i<iu; i++) {
213:       VecCopy(evec[i], vt1);
214:       for (j=il; j<iu; j++) {
215:         VecDot(evec[j],vt1,&dot);
216:         if (j == i) {
217:           dot = PetscAbsScalar(dot - 1);
218:         } else {
219:           dot = PetscAbsScalar(dot);
220:         }
221:         if (dot > dot_max) dot_max = dot;
222:         if (dot > tols[1]) {
223:           VecNorm(evec[i],NORM_INFINITY,&norm);
224:           PetscPrintf(PETSC_COMM_SELF,"|delta(%" PetscInt_FMT ",%" PetscInt_FMT ")|: %g, norm: %g\n",i,j,(double)dot,(double)norm);
225:         }
226:       }
227:     }
228:     PetscPrintf(PETSC_COMM_SELF,"    max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max);

230:   case 1:
231:     norm_max = 0.0;
232:     for (i = il; i< iu; i++) {
233:       MatMult(A, evec[i], vt1);
234:       VecCopy(evec[i], vt2);
235:       tmp  = -eval[i];
236:       VecAXPY(vt1,tmp,vt2);
237:       VecNorm(vt1, NORM_INFINITY, &norm);
238:       norm = PetscAbsScalar(norm);
239:       if (norm > norm_max) norm_max = norm;
240:       /* sniff, and bark if necessary */
241:       if (norm > tols[0]) {
242:         PetscPrintf(PETSC_COMM_SELF,"  residual violation: %" PetscInt_FMT ", resi: %g\n",i, (double)norm);
243:       }
244:     }
245:     PetscPrintf(PETSC_COMM_SELF,"    max_resi:                    %g\n", (double)norm_max);
246:     break;
247:   default:
248:     PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%" PetscInt_FMT " is not supported \n",cklvl);
249:   }
250:   VecDestroy(&vt2);
251:   VecDestroy(&vt1);
252:   return 0;
253: }

255: /*TEST

257:    build:
258:       requires: !complex

260:    test:
261:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
262:       args: -f ${DATAFILESPATH}/matrices/small
263:       output_file: output/ex116_1.out

265:    test:
266:       suffix: 2
267:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
268:       args: -f ${DATAFILESPATH}/matrices/small -test_syev -check_symmetry

270: TEST*/