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

      1: static char help[] = "Test LAPACK routine DSTEBZ() and DTEIN().  \n\n";

      3: #include <petscmat.h>
      4: #include <petscblaslapack.h>

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

     10: int main(int argc,char **args)
     11: {
     12: #if !defined(PETSC_USE_COMPLEX)
     14:   PetscReal      *work,tols[2];
     15:   PetscInt       i,j;
     16:   PetscBLASInt   n,il=1,iu=5,*iblock,*isplit,*iwork,nevs,*ifail,cklvl=2;
     17:   PetscMPIInt    size;
     18:   PetscBool      flg;
     19:   Vec            *evecs;
     20:   PetscScalar    *evecs_array,*D,*E,*evals;
     21:   Mat            T;
     22: #if !defined(PETSC_MISSING_LAPACK_DSTEBZ)
     23:   PetscReal vl=0.0,vu=4.0,tol=1.e-10;
     24:   PetscBLASInt  nsplit,info;
     25: #endif
     26: #endif

     28:   PetscInitialize(&argc,&args,(char*)0,help);
     29: #if defined(PETSC_USE_COMPLEX)
     30:   SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
     31: #else
     32:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     33:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");

     35:   n      = 100;
     36:   nevs   = iu - il;
     37:   PetscMalloc((3*n+1)*sizeof(PetscScalar),&D);
     38:   E      = D + n;
     39:   evals  = E + n;
     40:   PetscMalloc((5*n+1)*sizeof(PetscReal),&work);
     41:   PetscMalloc((3*n+1)*sizeof(PetscBLASInt),&iwork);
     42:   PetscMalloc((3*n+1)*sizeof(PetscBLASInt),&iblock);
     43:   isplit = iblock + n;

     45:   /* Set symmetric tridiagonal matrix */
     46:   for (i=0; i<n; i++) {
     47:     D[i] = 2.0;
     48:     E[i] = 1.0;
     49:   }

     51:   /* Solve eigenvalue problem: A*evec = eval*evec */
     52: #if defined(PETSC_MISSING_LAPACK_STEBZ)
     53:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
     54: #else
     55:   printf(" LAPACKstebz_: compute %d eigenvalues...\n",nevs);
     56:   LAPACKstebz_("I","E",&n,&vl,&vu,&il,&iu,&tol,(PetscReal*)D,(PetscReal*)E,&nevs,&nsplit,(PetscReal*)evals,iblock,isplit,work,iwork,&info);
     57:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstebz_ fails. info %d",info);
     58: #endif

     60:   printf(" LAPACKstein_: compute %d found eigenvectors...\n",nevs);
     61:   PetscMalloc(n*nevs*sizeof(PetscScalar),&evecs_array);
     62:   PetscMalloc(nevs*sizeof(PetscBLASInt),&ifail);
     63: #if defined(PETSC_MISSING_LAPACK_STEIN)
     64:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
     65: #else
     66:   LAPACKstein_(&n,(PetscReal*)D,(PetscReal*)E,&nevs,(PetscReal*)evals,iblock,isplit,evecs_array,&n,work,iwork,ifail,&info);
     67:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstein_ fails. info %d",info);
     68: #endif
     69:   /* View evals */
     70:   PetscOptionsHasName(NULL, "-eig_view", &flg);
     71:   if (flg) {
     72:     PetscPrintf(PETSC_COMM_SELF," %d evals: \n",nevs);
     73:     for (i=0; i<nevs; i++) PetscPrintf(PETSC_COMM_SELF,"%d  %G\n",i,evals[i]);
     74:   }

     76:   /* Check residuals and orthogonality */
     77:   MatCreate(PETSC_COMM_SELF,&T);
     78:   MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,n,n);
     79:   MatSetType(T,MATSBAIJ);
     80:   MatSetFromOptions(T);
     81:   MatSetUp(T);
     82:   for (i=0; i<n; i++) {
     83:     MatSetValues(T,1,&i,1,&i,&D[i],INSERT_VALUES);
     84:     if (i != n-1) {
     85:       j    = i+1;
     86:       MatSetValues(T,1,&i,1,&j,&E[i],INSERT_VALUES);
     87:     }
     88:   }
     89:   MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY);
     90:   MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY);

     92:   PetscMalloc((nevs+1)*sizeof(Vec),&evecs);
     93:   for (i=0; i<nevs; i++) {
     94:     VecCreate(PETSC_COMM_SELF,&evecs[i]);
     95:     VecSetSizes(evecs[i],PETSC_DECIDE,n);
     96:     VecSetFromOptions(evecs[i]);
     97:     VecPlaceArray(evecs[i],evecs_array+i*n);
     98:   }

    100:   tols[0] = 1.e-8;  tols[1] = 1.e-8;
    101:   CkEigenSolutions(cklvl,T,il-1,iu-1,evals,evecs,tols);

    103:   for (i=0; i<nevs; i++) {
    104:     VecResetArray(evecs[i]);
    105:   }

    107:   /* free space */

    109:   MatDestroy(&T);

    111:   for (i=0; i<nevs; i++) { VecDestroy(&evecs[i]);}
    112:   PetscFree(evecs);
    113:   PetscFree(D);
    114:   PetscFree(work);
    115:   PetscFree(iwork);
    116:   PetscFree(iblock);
    117:   PetscFree(evecs_array);
    118:   PetscFree(ifail);
    119:   PetscFinalize();
    120: #endif
    121:   return 0;
    122: }
    123: /*------------------------------------------------
    124:   Check the accuracy of the eigen solution
    125:   ----------------------------------------------- */
    126: /*
    127:   input:
    128:      cklvl      - check level:
    129:                     1: check residual
    130:                     2: 1 and check B-orthogonality locally
    131:      A          - matrix
    132:      il,iu      - lower and upper index bound of eigenvalues
    133:      eval, evec - eigenvalues and eigenvectors stored in this process
    134:      tols[0]    - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
    135:      tols[1]    - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
    136: */
    137: #undef DEBUG_CkEigenSolutions
    140: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscScalar *eval,Vec *evec,PetscReal *tols)
    141: {
    142:   PetscInt    ierr,i,j,nev;
    143:   Vec         vt1,vt2;  /* tmp vectors */
    144:   PetscReal   norm,norm_max;
    145:   PetscScalar dot,tmp;
    146:   PetscReal   dot_max;

    149:   nev = iu - il;
    150:   if (nev <= 0) return(0);

    152:   VecDuplicate(evec[0],&vt1);
    153:   VecDuplicate(evec[0],&vt2);

    155:   switch (cklvl) {
    156:   case 2:
    157:     dot_max = 0.0;
    158:     for (i = il; i<iu; i++) {
    159:       VecCopy(evec[i], vt1);
    160:       for (j=il; j<iu; j++) {
    161:         VecDot(evec[j],vt1,&dot);
    162:         if (j == i) {
    163:           dot = PetscAbsScalar(dot - 1.0);
    164:         } else {
    165:           dot = PetscAbsScalar(dot);
    166:         }
    167:         if (PetscAbsScalar(dot) > dot_max) dot_max = PetscAbsScalar(dot);
    168: #if defined(DEBUG_CkEigenSolutions)
    169:         if (dot > tols[1]) {
    170:           VecNorm(evec[i],NORM_INFINITY,&norm);
    171:           PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %G, norm: %G\n",i,j,dot,norm);
    172:         }
    173: #endif
    174:       }
    175:     }
    176:     PetscPrintf(PETSC_COMM_SELF,"    max|(x_j^T*x_i) - delta_ji|: %G\n",dot_max);

    178:   case 1:
    179:     norm_max = 0.0;
    180:     for (i = il; i< iu; i++) {
    181:       MatMult(A, evec[i], vt1);
    182:       VecCopy(evec[i], vt2);
    183:       tmp  = -eval[i];
    184:       VecAXPY(vt1,tmp,vt2);
    185:       VecNorm(vt1, NORM_INFINITY, &norm);
    186:       norm = PetscAbsScalar(norm);
    187:       if (norm > norm_max) norm_max = norm;
    188: #if defined(DEBUG_CkEigenSolutions)
    189:       /* sniff, and bark if necessary */
    190:       if (norm > tols[0]) {
    191:         printf("  residual violation: %d, resi: %g\n",i, norm);
    192:       }
    193: #endif
    194:     }
    195:     PetscPrintf(PETSC_COMM_SELF,"    max_resi:                    %G\n", norm_max);
    196:     break;
    197:   default:
    198:     PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
    199:   }

    201:   VecDestroy(&vt2);
    202:   VecDestroy(&vt1);
    203:   return(0);
    204: }