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