Actual source code: test25.c

slepc-3.17.2 2022-08-09
Report Typos and Errors
  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:   SlepcInitialize(&argc,&argv,(char*)0,help);
 35:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 36:   PetscOptionsGetBool(NULL,NULL,"-isnep",&isnep,NULL);
 37:   n = m*m;
 38:   k = 10;
 39:   PetscPrintf(PETSC_COMM_WORLD,"\nButterfly problem, n=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",n,m);
 40:   PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
 41:   PetscOptionsGetReal(NULL,NULL,"-radius",&radius,NULL);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

266: /*TEST

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

278: TEST*/