Actual source code: test25.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[] = "Test for DSPEP and DSNEP.\n\n";

     13: #include <slepcds.h>

     15: #define NMAT 5

     17: int main(int argc,char **argv)
     18: {
     19:   DS             ds;
     20:   FN             f[NMAT],qfun;
     21:   SlepcSC        sc;
     22:   PetscScalar    *A,*wr,*wi,*X,*y,*r,numer[NMAT],alpha;
     23:   PetscReal      c[10] = { 0.6, 1.3, 1.3, 0.1, 0.1, 1.2, 1.0, 1.0, 1.2, 1.0 };
     24:   PetscReal      tol,radius=1.5,re,im,nrm;
     25:   PetscInt       i,j,ii,jj,II,k,m=3,n,ld,nev,nfun,d,*inside;
     26:   PetscViewer    viewer;
     27:   PetscBool      verbose,isnep=PETSC_FALSE;
     28:   RG             rg;
     29:   DSMatType      mat[5]={DS_MAT_E0,DS_MAT_E1,DS_MAT_E2,DS_MAT_E3,DS_MAT_E4};
     30: #if !defined(PETSC_USE_COMPLEX)
     31:   PetscScalar    *yi,*ri,alphai=0.0,t;
     32: #endif

     34:   PetscFunctionBeginUser;
     35:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     36:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
     37:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-isnep",&isnep,NULL));
     38:   n = m*m;
     39:   k = 10;
     40:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nButterfly problem, n=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",n,m));
     41:   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
     42:   PetscCall(PetscOptionsGetReal(NULL,NULL,"-radius",&radius,NULL));

     44:   /* Create DS object */
     45:   PetscCall(DSCreate(PETSC_COMM_WORLD,&ds));
     46:   tol  = 1000*n*PETSC_MACHINE_EPSILON;
     47:   if (isnep) {
     48:     PetscCall(DSSetType(ds,DSNEP));
     49:     PetscCall(DSSetMethod(ds,1));
     50:     PetscCall(DSNEPSetRefine(ds,tol,PETSC_DECIDE));
     51:   } else PetscCall(DSSetType(ds,DSPEP));
     52:   PetscCall(DSSetFromOptions(ds));

     54:   /* Set functions (prior to DSAllocate) f_i=x^i */
     55:   if (isnep) {
     56:     numer[0] = 1.0;
     57:     for (j=1;j<NMAT;j++) numer[j] = 0.0;
     58:     for (i=0;i<NMAT;i++) {
     59:       PetscCall(FNCreate(PETSC_COMM_WORLD,&f[i]));
     60:       PetscCall(FNSetType(f[i],FNRATIONAL));
     61:       PetscCall(FNRationalSetNumerator(f[i],i+1,numer));
     62:     }
     63:     PetscCall(DSNEPSetFN(ds,NMAT,f));
     64:   } else PetscCall(DSPEPSetDegree(ds,NMAT-1));

     66:   /* Set dimensions */
     67:   ld = n+2;  /* test leading dimension larger than n */
     68:   PetscCall(DSAllocate(ds,ld));
     69:   PetscCall(DSSetDimensions(ds,n,0,0));

     71:   /* Set region (used only in method=1) */
     72:   PetscCall(RGCreate(PETSC_COMM_WORLD,&rg));
     73:   PetscCall(RGSetType(rg,RGELLIPSE));
     74:   PetscCall(RGEllipseSetParameters(rg,1.5,radius,.5));
     75:   PetscCall(RGSetFromOptions(rg));
     76:   if (isnep) PetscCall(DSNEPSetRG(ds,rg));

     78:   /* Set up viewer */
     79:   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
     80:   PetscCall(DSViewFromOptions(ds,NULL,"-ds_view"));
     81:   if (verbose) {
     82:     PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
     83:     /* Show info about functions */
     84:     if (isnep) {
     85:       PetscCall(DSNEPGetNumFN(ds,&nfun));
     86:       for (i=0;i<nfun;i++) {
     87:         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Function %" PetscInt_FMT ":\n",i));
     88:         PetscCall(DSNEPGetFN(ds,i,&qfun));
     89:         PetscCall(FNView(qfun,NULL));
     90:       }
     91:     }
     92:   }

     94:   /* Fill matrices */
     95:   /* A0 */
     96:   PetscCall(DSGetArray(ds,DS_MAT_E0,&A));
     97:   for (II=0;II<n;II++) {
     98:     i = II/m; j = II-i*m;
     99:     A[II+II*ld] = 4.0*c[0]/6.0+4.0*c[1]/6.0;
    100:     if (j>0) A[II+(II-1)*ld] = c[0]/6.0;
    101:     if (j<m-1) A[II+ld*(II+1)] = c[0]/6.0;
    102:     if (i>0) A[II+ld*(II-m)] = c[1]/6.0;
    103:     if (i<m-1) A[II+ld*(II+m)] = c[1]/6.0;
    104:   }
    105:   PetscCall(DSRestoreArray(ds,DS_MAT_E0,&A));

    107:   /* A1 */
    108:   PetscCall(DSGetArray(ds,DS_MAT_E1,&A));
    109:   for (II=0;II<n;II++) {
    110:     i = II/m; j = II-i*m;
    111:     if (j>0) A[II+ld*(II-1)] = c[2];
    112:     if (j<m-1) A[II+ld*(II+1)] = -c[2];
    113:     if (i>0) A[II+ld*(II-m)] = c[3];
    114:     if (i<m-1) A[II+ld*(II+m)] = -c[3];
    115:   }
    116:   PetscCall(DSRestoreArray(ds,DS_MAT_E1,&A));

    118:   /* A2 */
    119:   PetscCall(DSGetArray(ds,DS_MAT_E2,&A));
    120:   for (II=0;II<n;II++) {
    121:     i = II/m; j = II-i*m;
    122:     A[II+ld*II] = -2.0*c[4]-2.0*c[5];
    123:     if (j>0) A[II+ld*(II-1)] = c[4];
    124:     if (j<m-1) A[II+ld*(II+1)] = c[4];
    125:     if (i>0) A[II+ld*(II-m)] = c[5];
    126:     if (i<m-1) A[II+ld*(II+m)] = c[5];
    127:   }
    128:   PetscCall(DSRestoreArray(ds,DS_MAT_E2,&A));

    130:   /* A3 */
    131:   PetscCall(DSGetArray(ds,DS_MAT_E3,&A));
    132:   for (II=0;II<n;II++) {
    133:     i = II/m; j = II-i*m;
    134:     if (j>0) A[II+ld*(II-1)] = c[6];
    135:     if (j<m-1) A[II+ld*(II+1)] = -c[6];
    136:     if (i>0) A[II+ld*(II-m)] = c[7];
    137:     if (i<m-1) A[II+ld*(II+m)] = -c[7];
    138:   }
    139:   PetscCall(DSRestoreArray(ds,DS_MAT_E3,&A));

    141:   /* A4 */
    142:   PetscCall(DSGetArray(ds,DS_MAT_E4,&A));
    143:   for (II=0;II<n;II++) {
    144:     i = II/m; j = II-i*m;
    145:     A[II+ld*II] = 2.0*c[8]+2.0*c[9];
    146:     if (j>0) A[II+ld*(II-1)] = -c[8];
    147:     if (j<m-1) A[II+ld*(II+1)] = -c[8];
    148:     if (i>0) A[II+ld*(II-m)] = -c[9];
    149:     if (i<m-1) A[II+ld*(II+m)] = -c[9];
    150:   }
    151:   PetscCall(DSRestoreArray(ds,DS_MAT_E4,&A));

    153:   if (verbose) {
    154:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Initial - - - - - - - - -\n"));
    155:     PetscCall(DSView(ds,viewer));
    156:   }

    158:   /* Solve */
    159:   if (isnep) PetscCall(DSNEPGetMinimality(ds,&d));
    160:   else PetscCall(DSPEPGetDegree(ds,&d));
    161:   PetscCall(PetscCalloc3(n*d,&wr,n*d,&wi,n*d,&inside));
    162:   PetscCall(DSGetSlepcSC(ds,&sc));
    163:   sc->comparison    = SlepcCompareLargestMagnitude;
    164:   sc->comparisonctx = NULL;
    165:   sc->map           = NULL;
    166:   sc->mapobj        = NULL;
    167:   PetscCall(DSSolve(ds,wr,wi));
    168:   PetscCall(DSSort(ds,wr,wi,NULL,NULL,NULL));

    170:   if (verbose) {
    171:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"After solve - - - - - - - - -\n"));
    172:     PetscCall(DSView(ds,viewer));
    173:   }
    174:   if (isnep) {
    175:     PetscCall(DSGetDimensions(ds,NULL,NULL,NULL,&nev));
    176:     for (i=0;i<nev;i++) inside[i] = i;
    177:   } else {
    178:     PetscCall(RGCheckInside(rg,d*n,wr,wi,inside));
    179:     nev = 0;
    180:     for (i=0;i<d*n;i++) if (inside[i]>0) inside[nev++] = i;
    181:   }

    183:   /* Print computed eigenvalues */
    184:   PetscCall(PetscMalloc2(ld,&y,ld,&r));
    185: #if !defined(PETSC_USE_COMPLEX)
    186:   PetscCall(PetscMalloc2(ld,&yi,ld,&ri));
    187: #endif
    188:   PetscCall(DSVectors(ds,DS_MAT_X,NULL,NULL));
    189:   PetscCall(DSGetArray(ds,DS_MAT_X,&X));
    190:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed eigenvalues in the region: %" PetscInt_FMT "\n",nev));
    191:   for (i=0;i<nev;i++) {
    192: #if defined(PETSC_USE_COMPLEX)
    193:     re = PetscRealPart(wr[inside[i]]);
    194:     im = PetscImaginaryPart(wr[inside[i]]);
    195: #else
    196:     re = wr[inside[i]];
    197:     im = wi[inside[i]];
    198: #endif
    199:     PetscCall(PetscArrayzero(r,n));
    200: #if !defined(PETSC_USE_COMPLEX)
    201:     PetscCall(PetscArrayzero(ri,n));
    202: #endif
    203:     /* Residual */
    204:     alpha = 1.0;
    205:     for (k=0;k<NMAT;k++) {
    206:       PetscCall(DSGetArray(ds,mat[k],&A));
    207:       for (ii=0;ii<n;ii++) {
    208:         y[ii] = 0.0;
    209:         for (jj=0;jj<n;jj++) y[ii] += A[jj*ld+ii]*X[inside[i]*ld+jj];
    210:       }
    211: #if !defined(PETSC_USE_COMPLEX)
    212:       for (ii=0;ii<n;ii++) {
    213:         yi[ii] = 0.0;
    214:         for (jj=0;jj<n;jj++) yi[ii] += A[jj*ld+ii]*X[inside[i+1]*ld+jj];
    215:       }
    216: #endif
    217:       PetscCall(DSRestoreArray(ds,mat[k],&A));
    218:       if (isnep) PetscCall(FNEvaluateFunction(f[k],wr[inside[i]],&alpha));
    219:       for (ii=0;ii<n;ii++) r[ii] += alpha*y[ii];
    220: #if !defined(PETSC_USE_COMPLEX)
    221:       for (ii=0;ii<n;ii++) r[ii]  -= alphai*yi[ii];
    222:       for (ii=0;ii<n;ii++) ri[ii] += alpha*yi[ii]+alphai*y[ii];
    223: #endif
    224:       if (!isnep) {
    225: #if defined(PETSC_USE_COMPLEX)
    226:         alpha *= wr[inside[i]];
    227: #else
    228:         t      = alpha;
    229:         alpha  = alpha*re-alphai*im;
    230:         alphai = alphai*re+t*im;
    231: #endif
    232:       }
    233:     }
    234:     nrm = 0.0;
    235:     for (k=0;k<n;k++) {
    236: #if !defined(PETSC_USE_COMPLEX)
    237:       nrm += r[k]*r[k]+ri[k]*ri[k];
    238: #else
    239:       nrm += PetscRealPart(r[k]*PetscConj(r[k]));
    240: #endif
    241:     }
    242:     nrm = PetscSqrtReal(nrm);
    243:     if (nrm/SlepcAbsEigenvalue(wr[inside[i]],wi[inside[i]])>tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the residual norm of the %" PetscInt_FMT "-th computed eigenpair %g\n",i,(double)nrm));
    244:     if (PetscAbs(im)<1e-10) PetscCall(PetscViewerASCIIPrintf(viewer,"  %.5f\n",(double)re));
    245:     else PetscCall(PetscViewerASCIIPrintf(viewer,"  %.5f%+.5fi\n",(double)re,(double)im));
    246: #if !defined(PETSC_USE_COMPLEX)
    247:     if (im!=0.0) i++;
    248:     if (PetscAbs(im)<1e-10) PetscCall(PetscViewerASCIIPrintf(viewer,"  %.5f\n",(double)re));
    249:     else PetscCall(PetscViewerASCIIPrintf(viewer,"  %.5f%+.5fi\n",(double)re,(double)-im));
    250: #endif
    251:   }
    252:   PetscCall(DSRestoreArray(ds,DS_MAT_X,&X));
    253:   PetscCall(PetscFree3(wr,wi,inside));
    254:   PetscCall(PetscFree2(y,r));
    255: #if !defined(PETSC_USE_COMPLEX)
    256:   PetscCall(PetscFree2(yi,ri));
    257: #endif
    258:   if (isnep) {
    259:     for (i=0;i<NMAT;i++) PetscCall(FNDestroy(&f[i]));
    260:   }
    261:   PetscCall(DSDestroy(&ds));
    262:   PetscCall(RGDestroy(&rg));
    263:   PetscCall(SlepcFinalize());
    264:   return 0;
    265: }

    267: /*TEST

    269:    testset:
    270:       filter: sed -e "s/[+-]\([0-9]\.[0-9]*i\)/+-\\1/" | sed -e "s/56808/56807/" | sed -e "s/34719/34720/"
    271:       output_file: output/test25_1.out
    272:       test:
    273:          suffix: 1
    274:       test:
    275:          suffix: 2
    276:          args: -isnep
    277:          requires: complex !single

    279: TEST*/