Actual source code: ex99.c
petsc-3.4.2 2013-07-02

      1: static char help[] = "Test LAPACK routine DSYGV() or DSYGVX(). \n\
      2: Reads PETSc matrix A and B (or create B=I), \n\
      3: then computes selected eigenvalues, and optionally, eigenvectors of \n\
      4: a real generalized symmetric-definite eigenproblem \n\
      5:  A*x = lambda*B*x \n\
      6: Input parameters include\n\
      7:   -f0 <input_file> : first file to load (small system)\n\
      8:   -fA <input_file> -fB <input_file>: second files to load (larger system) \n\
      9: e.g. ./ex99 -f0 $D/small -fA $D/Eigdftb/dftb_bin/diamond_xxs_A -fB $D/Eigdftb/dftb_bin/diamond_xxs_B -mat_getrow_uppertriangular,\n\
     10:      where $D = /home/petsc/datafiles/matrices/Eigdftb/dftb_bin\n\n";

     12: #include <petscmat.h>
     13: #include <../src/mat/impls/sbaij/seq/sbaij.h>
     14: #include <petscblaslapack.h>

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

     20: PetscInt main(PetscInt argc,char **args)
     21: {
     22:   Mat            A,B,A_dense,B_dense,mats[2],A_sp;
     23:   Vec            *evecs;
     24:   PetscViewer    fd;                /* viewer */
     25:   char           file[3][PETSC_MAX_PATH_LEN];     /* input file name */
     26:   PetscBool      flg,flgA=PETSC_FALSE,flgB=PETSC_FALSE,TestSYGVX=PETSC_TRUE;
     28:   PetscBool      preload=PETSC_TRUE,isSymmetric;
     29:   PetscScalar    sigma,one=1.0,*arrayA,*arrayB,*evecs_array,*work,*evals;
     30:   PetscMPIInt    size;
     31:   PetscInt       m,n,i,j,nevs,il,iu;
     32:   PetscLogStage  stages[2];
     33:   PetscReal      vl,vu,abstol=1.e-8;
     34:   PetscBLASInt   *iwork,*ifail,lone=1,lwork,lierr,bn;
     35:   PetscInt       ievbd_loc[2],offset=0,cklvl=2;
     36:   PetscReal      tols[2];
     37:   Mat_SeqSBAIJ   *sbaij;
     38:   PetscScalar    *aa;
     39:   PetscInt       *ai,*aj;
     40:   PetscInt       nzeros[2],nz;
     41:   PetscReal      ratio;

     43:   PetscInitialize(&argc,&args,(char*)0,help);
     44:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     45:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
     46:   PetscLogStageRegister("EigSolve",&stages[0]);
     47:   PetscLogStageRegister("EigCheck",&stages[1]);

     49:   /* Determine files from which we read the two matrices */
     50:   PetscOptionsGetString(NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
     51:   if (!flg) {
     52:     PetscOptionsGetString(NULL,"-fA",file[0],PETSC_MAX_PATH_LEN,&flgA);
     53:     if (!flgA) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -fA or -fB options");
     54:     PetscOptionsGetString(NULL,"-fB",file[1],PETSC_MAX_PATH_LEN,&flgB);
     55:     preload = PETSC_FALSE;
     56:   } else {
     57:     PetscOptionsGetString(NULL,"-fA",file[1],PETSC_MAX_PATH_LEN,&flgA);
     58:     if (!flgA) preload = PETSC_FALSE; /* don't bother with second system */
     59:     PetscOptionsGetString(NULL,"-fB",file[2],PETSC_MAX_PATH_LEN,&flgB);
     60:   }

     62:   PetscPreLoadBegin(preload,"Load system");
     63:   /* Load matrices */
     64:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);
     65:   MatCreate(PETSC_COMM_WORLD,&A);
     66:   MatSetType(A,MATSBAIJ);
     67:   MatLoad(A,fd);
     68:   PetscViewerDestroy(&fd);
     69:   MatGetSize(A,&m,&n);
     70:   if ((flgB && PetscPreLoadIt) || (flgB && !preload)) {
     71:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt+1],FILE_MODE_READ,&fd);
     72:     MatCreate(PETSC_COMM_WORLD,&B);
     73:     MatSetType(B,MATSBAIJ);
     74:     MatLoad(B,fd);
     75:     PetscViewerDestroy(&fd);
     76:   } else {   /* create B=I */
     77:     MatCreate(PETSC_COMM_WORLD,&B);
     78:     MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
     79:     MatSetType(B,MATSEQSBAIJ);
     80:     MatSetFromOptions(B);
     81:     for (i=0; i<m; i++) {
     82:       MatSetValues(B,1,&i,1,&i,&one,INSERT_VALUES);
     83:     }
     84:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
     85:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
     86:   }

     88:   /* Add a shift to A */
     89:   PetscOptionsGetScalar(NULL,"-mat_sigma",&sigma,&flg);
     90:   if (flg) {
     91:     MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN);   /* A <- sigma*B + A */
     92:   }

     94:   /* Check whether A is symmetric */
     95:   PetscOptionsHasName(NULL, "-check_symmetry", &flg);
     96:   if (flg) {
     97:     Mat Trans;
     98:     MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
     99:     MatEqual(A, Trans, &isSymmetric);
    100:     if (!isSymmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A must be symmetric");
    101:     MatDestroy(&Trans);
    102:     if (flgB && PetscPreLoadIt) {
    103:       MatTranspose(B,MAT_INITIAL_MATRIX, &Trans);
    104:       MatEqual(B, Trans, &isSymmetric);
    105:       if (!isSymmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"B must be symmetric");
    106:       MatDestroy(&Trans);
    107:     }
    108:   }

    110:   /* View small entries of A */
    111:   PetscOptionsHasName(NULL, "-Asp_view", &flg);
    112:   if (flg) {
    113:     MatCreate(PETSC_COMM_SELF,&A_sp);
    114:     MatSetSizes(A_sp,PETSC_DECIDE,PETSC_DECIDE,m,n);
    115:     MatSetType(A_sp,MATSEQSBAIJ);

    117:     tols[0]   = 1.e-6, tols[1] = 1.e-9;
    118:     sbaij     = (Mat_SeqSBAIJ*)A->data;
    119:     ai        = sbaij->i;
    120:     aj        = sbaij->j;
    121:     aa        = sbaij->a;
    122:     nzeros[0] = nzeros[1] = 0;
    123:     for (i=0; i<m; i++) {
    124:       nz = ai[i+1] - ai[i];
    125:       for (j=0; j<nz; j++) {
    126:         if (PetscAbsScalar(*aa)<tols[0]) {
    127:           MatSetValues(A_sp,1,&i,1,aj,aa,INSERT_VALUES);
    128:           nzeros[0]++;
    129:         }
    130:         if (PetscAbsScalar(*aa)<tols[1]) nzeros[1]++;
    131:         aa++; aj++;
    132:       }
    133:     }
    134:     MatAssemblyBegin(A_sp,MAT_FINAL_ASSEMBLY);
    135:     MatAssemblyEnd(A_sp,MAT_FINAL_ASSEMBLY);

    137:     MatDestroy(&A_sp);

    139:     ratio = (PetscReal)nzeros[0]/sbaij->nz;
    140:     PetscPrintf(PETSC_COMM_SELF," %d matrix entries < %e, ratio %G of %d nonzeros\n",nzeros[0],tols[0],ratio,sbaij->nz);
    141:     PetscPrintf(PETSC_COMM_SELF," %d matrix entries < %e\n",nzeros[1],tols[1]);
    142:   }

    144:   /* Convert aij matrix to MatSeqDense for LAPACK */
    145:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
    146:   if (!flg) {
    147:     MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);
    148:   }
    149:   PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
    150:   if (!flg) {MatConvert(B,MATSEQDENSE,MAT_INITIAL_MATRIX,&B_dense);}

    152:   /* Solve eigenvalue problem: A*x = lambda*B*x */
    153:   /*============================================*/
    154:   PetscBLASIntCast(8*n,&lwork);
    155:   PetscBLASIntCast(n,&bn);
    156:   PetscMalloc(n*sizeof(PetscScalar),&evals);
    157:   PetscMalloc(lwork*sizeof(PetscScalar),&work);
    158:   MatDenseGetArray(A_dense,&arrayA);
    159:   MatDenseGetArray(B_dense,&arrayB);

    161:   if (!TestSYGVX) {   /* test sygv()  */
    162:     evecs_array = arrayA;
    163:     LAPACKsygv_(&lone,"V","U",&bn,arrayA,&bn,arrayB,&bn,evals,work,&lwork,&lierr);
    164:     nevs = m;
    165:     il   =1;
    166:   } else {   /* test sygvx()  */
    167:     il    = 1;
    168:     PetscBLASIntCast(.6*m,&iu);
    169:     PetscMalloc((m*n+1)*sizeof(PetscScalar),&evecs_array);
    170:     PetscMalloc((6*n+1)*sizeof(PetscBLASInt),&iwork);
    171:     ifail = iwork + 5*n;
    172:     if (PetscPreLoadIt) {PetscLogStagePush(stages[0]);}
    173:     /* in the case "I", vl and vu are not referenced */
    174:     LAPACKsygvx_(&lone,"V","I","U",&bn,arrayA,&bn,arrayB,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&n,work,&lwork,iwork,ifail,&lierr);
    175:     if (PetscPreLoadIt) PetscLogStagePop();
    176:     PetscFree(iwork);
    177:   }
    178:   MatDenseRestoreArray(A,&arrayA);
    179:   MatDenseRestoreArray(B,&arrayB);

    181:   if (nevs <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs);
    182:   /* View evals */
    183:   PetscOptionsHasName(NULL, "-eig_view", &flg);
    184:   if (flg) {
    185:     printf(" %d evals: \n",nevs);
    186:     for (i=0; i<nevs; i++) printf("%d  %G\n",i+il,evals[i]);
    187:   }

    189:   /* Check residuals and orthogonality */
    190:   if (PetscPreLoadIt) {
    191:     mats[0] = A; mats[1] = B;
    192:     one     = (PetscInt)one;
    193:     PetscMalloc((nevs+1)*sizeof(Vec),&evecs);
    194:     for (i=0; i<nevs; i++) {
    195:       VecCreate(PETSC_COMM_SELF,&evecs[i]);
    196:       VecSetSizes(evecs[i],PETSC_DECIDE,n);
    197:       VecSetFromOptions(evecs[i]);
    198:       VecPlaceArray(evecs[i],evecs_array+i*n);
    199:     }

    201:     ievbd_loc[0] = 0; ievbd_loc[1] = nevs-1;
    202:     tols[0]      = 1.e-8;  tols[1] = 1.e-8;

    204:     PetscLogStagePush(stages[1]);
    205:     CkEigenSolutions(&cklvl,mats,evals,evecs,ievbd_loc,&offset,tols);
    206:     PetscLogStagePop();
    207:     for (i=0; i<nevs; i++) { VecDestroy(&evecs[i]);}
    208:     PetscFree(evecs);
    209:   }

    211:   /* Free work space. */
    212:   if (TestSYGVX) {PetscFree(evecs_array);}

    214:   PetscFree(evals);
    215:   PetscFree(work);

    217:   MatDestroy(&A_dense);
    218:   MatDestroy(&B_dense);
    219:   MatDestroy(&B);
    220:   MatDestroy(&A);

    222:   PetscPreLoadEnd();
    223:   PetscFinalize();
    224:   return 0;
    225: }
    226: /*------------------------------------------------
    227:   Check the accuracy of the eigen solution
    228:   ----------------------------------------------- */
    229: /*
    230:   input:
    231:      cklvl      - check level:
    232:                     1: check residual
    233:                     2: 1 and check B-orthogonality locally
    234:      mats       - matrix pencil
    235:      eval, evec - eigenvalues and eigenvectors stored in this process
    236:      ievbd_loc  - local eigenvalue bounds, see eigc()
    237:      offset     - see eigc()
    238:      tols[0]    - reporting tol_res: || A evec[i] - eval[i] B evec[i]||
    239:      tols[1]    - reporting tol_orth: evec[i] B evec[j] - delta_ij
    240: */
    241: #undef DEBUG_CkEigenSolutions
    244: PetscErrorCode CkEigenSolutions(PetscInt *fcklvl,Mat *mats,PetscReal *eval,Vec *evec,PetscInt *ievbd_loc,PetscInt *offset,PetscReal *tols)
    245: {
    246:   PetscInt  ierr,cklvl=*fcklvl,nev_loc,i,j;
    247:   Mat       A=mats[0], B=mats[1];
    248:   Vec       vt1,vt2;    /* tmp vectors */
    249:   PetscReal norm,tmp,dot,norm_max,dot_max;

    252:   nev_loc = ievbd_loc[1] - ievbd_loc[0];
    253:   if (nev_loc == 0) return(0);

    255:   nev_loc += (*offset);
    256:   VecDuplicate(evec[*offset],&vt1);
    257:   VecDuplicate(evec[*offset],&vt2);

    259:   switch (cklvl) {
    260:   case 2:
    261:     dot_max = 0.0;
    262:     for (i = *offset; i<nev_loc; i++) {
    263:       MatMult(B, evec[i], vt1);
    264:       for (j=i; j<nev_loc; j++) {
    265:         VecDot(evec[j],vt1,&dot);
    266:         if (j == i) {
    267:           dot = PetscAbsScalar(dot - 1.0);
    268:         } else {
    269:           dot = PetscAbsScalar(dot);
    270:         }
    271:         if (dot > dot_max) dot_max = dot;
    272: #if defined(DEBUG_CkEigenSolutions)
    273:         if (dot > tols[1]) {
    274:           VecNorm(evec[i],NORM_INFINITY,&norm);
    275:           PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %G, norm: %G\n",i,j,dot,norm);
    276:         }
    277: #endif
    278:       } /* for (j=i; j<nev_loc; j++) */
    279:     }
    280:     PetscPrintf(PETSC_COMM_SELF,"    max|(x_j*B*x_i) - delta_ji|: %G\n",dot_max);

    282:   case 1:
    283:     norm_max = 0.0;
    284:     for (i = *offset; i< nev_loc; i++) {
    285:       MatMult(A, evec[i], vt1);
    286:       MatMult(B, evec[i], vt2);
    287:       tmp  = -eval[i];
    288:       VecAXPY(vt1,tmp,vt2);
    289:       VecNorm(vt1, NORM_INFINITY, &norm);
    290:       norm = PetscAbsScalar(norm);
    291:       if (norm > norm_max) norm_max = norm;
    292: #if defined(DEBUG_CkEigenSolutions)
    293:       /* sniff, and bark if necessary */
    294:       if (norm > tols[0]) {
    295:         printf("  residual violation: %d, resi: %g\n",i, norm);
    296:       }
    297: #endif
    298:     }

    300:     PetscPrintf(PETSC_COMM_SELF,"    max_resi:                    %G\n", norm_max);

    302:     break;
    303:   default:
    304:     PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
    305:   }
    306:   VecDestroy(&vt2);
    307:   VecDestroy(&vt1);
    308:   return(0);
    309: }