Actual source code: test25.c
slepc-3.21.1 2024-04-26
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,(char*)0,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*/