Actual source code: test29.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */

     11: static char help[] = "Illustrates the computation of left eigenvectors for generalized eigenproblems.\n\n"
     12:   "The command line options are:\n"
     13:   "  -f1 <filename> -f2 <filename>, PETSc binary files containing A and B\n\n";

     15: #include <slepceps.h>

     17: /*
     18:    User-defined routines
     19: */
     20: PetscErrorCode ComputeResidualNorm(Mat,Mat,PetscBool,PetscScalar,PetscScalar,Vec,Vec,Vec*,PetscReal*);
     21: PetscErrorCode CheckBiorthogonality(Vec*,Vec*,PetscScalar*,PetscInt,Mat,PetscReal*);

     23: int main(int argc,char **argv)
     24: {
     25:   Mat            A,B;
     26:   EPS            eps;
     27:   EPSType        type;
     28:   PetscInt       i,nconv;
     29:   PetscBool      twosided,flg;
     30:   PetscReal      nrmr,nrml=0.0,re,im,lev;
     31:   PetscScalar    *kr,*ki;
     32:   Vec            t,*xr,*xi,*yr,*yi,*z;
     33:   char           filename[PETSC_MAX_PATH_LEN];
     34:   PetscViewer    viewer;

     36:   PetscFunctionBeginUser;
     37:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));

     39:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     40:         Load the matrices that define the eigensystem, Ax=kBx
     41:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     43:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n"));
     44:   PetscCall(PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg));
     45:   PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -f1 option");

     47: #if defined(PETSC_USE_COMPLEX)
     48:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n"));
     49: #else
     50:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n"));
     51: #endif
     52:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
     53:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
     54:   PetscCall(MatSetFromOptions(A));
     55:   PetscCall(MatLoad(A,viewer));
     56:   PetscCall(PetscViewerDestroy(&viewer));

     58:   PetscCall(PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg));
     59:   if (flg) {
     60:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
     61:     PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
     62:     PetscCall(MatSetFromOptions(B));
     63:     PetscCall(MatLoad(B,viewer));
     64:     PetscCall(PetscViewerDestroy(&viewer));
     65:   } else {
     66:     PetscCall(PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n"));
     67:     B = NULL;
     68:   }
     69:   PetscCall(MatCreateVecs(A,NULL,&t));

     71:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     72:                 Create the eigensolver and set various options
     73:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     75:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
     76:   PetscCall(EPSSetOperators(eps,A,B));

     78:   /* use a two-sided algorithm to compute left eigenvectors as well */
     79:   PetscCall(EPSSetTwoSided(eps,PETSC_TRUE));

     81:   /* allow user to change settings at run time */
     82:   PetscCall(EPSSetFromOptions(eps));
     83:   PetscCall(EPSGetTwoSided(eps,&twosided));

     85:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     86:                       Solve the eigensystem
     87:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     89:   PetscCall(EPSSolve(eps));

     91:   /*
     92:      Optional: Get some information from the solver and display it
     93:   */
     94:   PetscCall(EPSGetType(eps,&type));
     95:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));

     97:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     98:                     Display solution and clean up
     99:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    101:   /*
    102:      Get number of converged approximate eigenpairs
    103:   */
    104:   PetscCall(EPSGetConverged(eps,&nconv));
    105:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %" PetscInt_FMT "\n\n",nconv));
    106:   PetscCall(PetscMalloc2(nconv,&kr,nconv,&ki));
    107:   PetscCall(VecDuplicateVecs(t,3,&z));
    108:   PetscCall(VecDuplicateVecs(t,nconv,&xr));
    109:   PetscCall(VecDuplicateVecs(t,nconv,&xi));
    110:   if (twosided) {
    111:     PetscCall(VecDuplicateVecs(t,nconv,&yr));
    112:     PetscCall(VecDuplicateVecs(t,nconv,&yi));
    113:   }

    115:   if (nconv>0) {
    116:     /*
    117:        Display eigenvalues and relative errors
    118:     */
    119:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,
    120:          "           k            ||Ax-kBx||         ||y'A-y'Bk||\n"
    121:          "   ---------------- ------------------ ------------------\n"));

    123:     for (i=0;i<nconv;i++) {
    124:       /*
    125:         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
    126:         ki (imaginary part)
    127:       */
    128:       PetscCall(EPSGetEigenpair(eps,i,&kr[i],&ki[i],xr[i],xi[i]));
    129:       if (twosided) PetscCall(EPSGetLeftEigenvector(eps,i,yr[i],yi[i]));
    130:       /*
    131:          Compute the residual norms associated to each eigenpair
    132:       */
    133:       PetscCall(ComputeResidualNorm(A,B,PETSC_FALSE,kr[i],ki[i],xr[i],xi[i],z,&nrmr));
    134:       if (twosided) PetscCall(ComputeResidualNorm(A,B,PETSC_TRUE,kr[i],ki[i],yr[i],yi[i],z,&nrml));

    136: #if defined(PETSC_USE_COMPLEX)
    137:       re = PetscRealPart(kr[i]);
    138:       im = PetscImaginaryPart(kr[i]);
    139: #else
    140:       re = kr[i];
    141:       im = ki[i];
    142: #endif
    143:       if (im!=0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %8f%+8fi %12g       %12g\n",(double)re,(double)im,(double)nrmr,(double)nrml));
    144:       else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g       %12g\n",(double)re,(double)nrmr,(double)nrml));
    145:     }
    146:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
    147:     /*
    148:        Check bi-orthogonality of eigenvectors
    149:     */
    150:     if (twosided) {
    151:       PetscCall(CheckBiorthogonality(xr,yr,ki,nconv,B,&lev));
    152:       if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
    153:       else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
    154:     }
    155:   }

    157:   PetscCall(EPSDestroy(&eps));
    158:   PetscCall(MatDestroy(&A));
    159:   PetscCall(MatDestroy(&B));
    160:   PetscCall(VecDestroy(&t));
    161:   PetscCall(PetscFree2(kr,ki));
    162:   PetscCall(VecDestroyVecs(3,&z));
    163:   PetscCall(VecDestroyVecs(nconv,&xr));
    164:   PetscCall(VecDestroyVecs(nconv,&xi));
    165:   if (twosided) {
    166:     PetscCall(VecDestroyVecs(nconv,&yr));
    167:     PetscCall(VecDestroyVecs(nconv,&yi));
    168:   }
    169:   PetscCall(SlepcFinalize());
    170:   return 0;
    171: }

    173: /*
    174:    ComputeResidualNorm - Computes the norm of the residual vector
    175:    associated with an eigenpair.

    177:    Input Parameters:
    178:      trans - whether A' must be used instead of A
    179:      kr,ki - eigenvalue
    180:      xr,xi - eigenvector
    181:      z     - three work vectors (the second one not referenced in complex scalars)
    182: */
    183: PetscErrorCode ComputeResidualNorm(Mat A,Mat B,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
    184: {
    185:   Vec            u,w=NULL;
    186:   PetscScalar    alpha;
    187: #if !defined(PETSC_USE_COMPLEX)
    188:   Vec            v;
    189:   PetscReal      ni,nr;
    190: #endif
    191:   PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;

    193:   PetscFunctionBegin;
    194:   u = z[0];
    195:   if (B) w = z[2];

    197: #if !defined(PETSC_USE_COMPLEX)
    198:   v = z[1];
    199:   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
    200: #endif
    201:     PetscCall((*matmult)(A,xr,u));                          /* u=A*x */
    202:     if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
    203:       if (B) PetscCall((*matmult)(B,xr,w));             /* w=B*x */
    204:       else w = xr;
    205:       alpha = trans? -PetscConj(kr): -kr;
    206:       PetscCall(VecAXPY(u,alpha,w));                        /* u=A*x-k*B*x */
    207:     }
    208:     PetscCall(VecNorm(u,NORM_2,norm));
    209: #if !defined(PETSC_USE_COMPLEX)
    210:   } else {
    211:     PetscCall((*matmult)(A,xr,u));                          /* u=A*xr */
    212:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
    213:       if (B) PetscCall((*matmult)(B,xr,v));             /* v=B*xr */
    214:       else PetscCall(VecCopy(xr,v));
    215:       PetscCall(VecAXPY(u,-kr,v));                          /* u=A*xr-kr*B*xr */
    216:       if (B) PetscCall((*matmult)(B,xi,w));             /* w=B*xi */
    217:       else w = xi;
    218:       PetscCall(VecAXPY(u,trans?-ki:ki,w));                 /* u=A*xr-kr*B*xr+ki*B*xi */
    219:     }
    220:     PetscCall(VecNorm(u,NORM_2,&nr));
    221:     PetscCall((*matmult)(A,xi,u));                          /* u=A*xi */
    222:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
    223:       PetscCall(VecAXPY(u,-kr,w));                          /* u=A*xi-kr*B*xi */
    224:       PetscCall(VecAXPY(u,trans?ki:-ki,v));                 /* u=A*xi-kr*B*xi-ki*B*xr */
    225:     }
    226:     PetscCall(VecNorm(u,NORM_2,&ni));
    227:     *norm = SlepcAbsEigenvalue(nr,ni);
    228:   }
    229: #endif
    230:   PetscFunctionReturn(PETSC_SUCCESS);
    231: }

    233: /*
    234:    CheckOrthogonality - Similar to VecCheckOrthogonality but taking into account complex eigenvalues.

    236:    Input Parameters:
    237:      V - right eigenvectors
    238:      W - left eigenvectors
    239:      eigi - imaginary part of eigenvalues (only to check if the corresponding vector is complex)
    240:      n - number of V,W vectors
    241:      B - matrix defining the inner product
    242:    Output Parameter:
    243:      lev - level of bi-orthogonality
    244: */
    245: PetscErrorCode CheckBiorthogonality(Vec *V,Vec *W,PetscScalar *ki,PetscInt n,Mat B,PetscReal *lev)
    246: {
    247:   PetscInt       i,j;
    248:   PetscScalar    val1;
    249:   PetscBool      wcmplx;
    250:   Vec            wr;
    251: #if !defined(PETSC_USE_COMPLEX)
    252:   Vec            wi;
    253:   PetscScalar    val2,dotr,doti;
    254: #endif

    256:   PetscFunctionBegin;
    257:   if (n<=0) PetscFunctionReturn(PETSC_SUCCESS);
    258:   PetscCall(VecDuplicate(V[0],&wr));
    259: #if !defined(PETSC_USE_COMPLEX)
    260:   PetscCall(VecDuplicate(V[0],&wi));
    261: #endif
    262:   *lev = 0.0;
    263:   for (i=0;i<n;i++) {
    264:     wcmplx = PETSC_FALSE;
    265:     PetscCall(MatMult(B,W[i],wr));
    266: #if !defined(PETSC_USE_COMPLEX)
    267:     if (ki[i] != 0.0) {
    268:       PetscCall(MatMult(B,W[i+1],wi));
    269:       wcmplx = PETSC_TRUE;
    270:     }
    271: #endif
    272:     for (j=0;j<n;j++) {
    273:       if (j==i || (wcmplx && j==i+1)) continue;
    274:       PetscCall(VecDot(wr,V[j],&val1));
    275: #if !defined(PETSC_USE_COMPLEX)
    276:       dotr = val1;
    277:       if (ki[j] != 0.0) {  /* V complex */
    278:         if (wcmplx) {
    279:           PetscCall(VecDot(wi,V[j+1],&val2));
    280:           dotr = val1+val2;
    281:           PetscCall(VecDot(wi,V[j],&val1));
    282:           PetscCall(VecDot(wr,V[j+1],&val2));
    283:           doti = val1-val2;
    284:         } else {
    285:           PetscCall(VecDot(wr,V[j+1],&val1));
    286:           doti = -val1;
    287:         }
    288:       } else if (wcmplx) {
    289:         PetscCall(VecDot(wi,V[j],&doti));
    290:       } else doti = 0.0;
    291:       *lev = PetscMax(*lev,SlepcAbsEigenvalue(dotr,doti));
    292:       if (ki[j] != 0.0) j++;
    293: #else
    294:       *lev = PetscMax(*lev,PetscAbsScalar(val1));
    295: #endif
    296:     }
    297:     if (wcmplx) i++;
    298:   }
    299:   PetscCall(VecDestroy(&wr));
    300: #if !defined(PETSC_USE_COMPLEX)
    301:   PetscCall(VecDestroy(&wi));
    302: #endif
    303:   PetscFunctionReturn(PETSC_SUCCESS);
    304: }

    306: /*TEST

    308:    testset:
    309:       args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc -eps_nev 4 -st_type sinvert -eps_target -190000
    310:       filter: grep -v "method" | sed -e "s/[+-]0\.0*i//g" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
    311:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
    312:       test:
    313:          suffix: 1
    314:       test:
    315:          suffix: 1_cmplxvals
    316:          args: -eps_target -250000
    317:       test:
    318:          suffix: 1_rqi
    319:          args: -eps_type power -eps_power_shift_type rayleigh -eps_nev 2 -eps_target -2000
    320:       test:
    321:          suffix: 1_rqi_singular
    322:          args: -eps_type power -eps_power_shift_type rayleigh -eps_nev 1 -eps_target -195500

    324:    test:
    325:       suffix: 2
    326:       args: -f1 ${DATAFILESPATH}/matrices/complex/mhd1280a.petsc -f2 ${DATAFILESPATH}/matrices/complex/mhd1280b.petsc -eps_nev 6 -eps_tol 1e-11
    327:       filter: sed -e "s/-892/+892/" | sed -e "s/-759/+759/" | sed -e "s/-674/+674/" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
    328:       requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
    329:       timeoutfactor: 2

    331:    test:
    332:       suffix: 3
    333:       args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc -eps_nev 4 -st_type sinvert -eps_target -250000
    334:       filter: grep -v "method" | sed -e "s/[+-]0\.0*i//g" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g" | grep -v bi-orthog
    335:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)

    337: TEST*/