Actual source code: test25.c
slepc-3.18.0 2022-10-01
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
35: SlepcInitialize(&argc,&argv,(char*)0,help);
36: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
37: PetscOptionsGetBool(NULL,NULL,"-isnep",&isnep,NULL);
38: n = m*m;
39: k = 10;
40: PetscPrintf(PETSC_COMM_WORLD,"\nButterfly problem, n=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",n,m);
41: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
42: PetscOptionsGetReal(NULL,NULL,"-radius",&radius,NULL);
44: /* Create DS object */
45: DSCreate(PETSC_COMM_WORLD,&ds);
46: tol = 1000*n*PETSC_MACHINE_EPSILON;
47: if (isnep) {
48: DSSetType(ds,DSNEP);
49: DSSetMethod(ds,1);
50: DSNEPSetRefine(ds,tol,PETSC_DECIDE);
51: } else DSSetType(ds,DSPEP);
52: 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: FNCreate(PETSC_COMM_WORLD,&f[i]);
60: FNSetType(f[i],FNRATIONAL);
61: FNRationalSetNumerator(f[i],i+1,numer);
62: }
63: DSNEPSetFN(ds,NMAT,f);
64: } else DSPEPSetDegree(ds,NMAT-1);
66: /* Set dimensions */
67: ld = n+2; /* test leading dimension larger than n */
68: DSAllocate(ds,ld);
69: DSSetDimensions(ds,n,0,0);
71: /* Set region (used only in method=1) */
72: RGCreate(PETSC_COMM_WORLD,&rg);
73: RGSetType(rg,RGELLIPSE);
74: RGEllipseSetParameters(rg,1.5,radius,.5);
75: RGSetFromOptions(rg);
76: if (isnep) DSNEPSetRG(ds,rg);
78: /* Set up viewer */
79: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer);
80: DSViewFromOptions(ds,NULL,"-ds_view");
81: if (verbose) {
82: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
83: /* Show info about functions */
84: if (isnep) {
85: DSNEPGetNumFN(ds,&nfun);
86: for (i=0;i<nfun;i++) {
87: PetscPrintf(PETSC_COMM_WORLD,"Function %" PetscInt_FMT ":\n",i);
88: DSNEPGetFN(ds,i,&qfun);
89: FNView(qfun,NULL);
90: }
91: }
92: }
94: /* Fill matrices */
95: /* A0 */
96: 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: DSRestoreArray(ds,DS_MAT_E0,&A);
107: /* A1 */
108: 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: DSRestoreArray(ds,DS_MAT_E1,&A);
118: /* A2 */
119: 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: DSRestoreArray(ds,DS_MAT_E2,&A);
130: /* A3 */
131: 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: DSRestoreArray(ds,DS_MAT_E3,&A);
141: /* A4 */
142: 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: DSRestoreArray(ds,DS_MAT_E4,&A);
153: if (verbose) {
154: PetscPrintf(PETSC_COMM_WORLD,"Initial - - - - - - - - -\n");
155: DSView(ds,viewer);
156: }
158: /* Solve */
159: if (isnep) DSNEPGetMinimality(ds,&d);
160: else DSPEPGetDegree(ds,&d);
161: PetscCalloc3(n*d,&wr,n*d,&wi,n*d,&inside);
162: DSGetSlepcSC(ds,&sc);
163: sc->comparison = SlepcCompareLargestMagnitude;
164: sc->comparisonctx = NULL;
165: sc->map = NULL;
166: sc->mapobj = NULL;
167: DSSolve(ds,wr,wi);
168: DSSort(ds,wr,wi,NULL,NULL,NULL);
170: if (verbose) {
171: PetscPrintf(PETSC_COMM_WORLD,"After solve - - - - - - - - -\n");
172: DSView(ds,viewer);
173: }
174: if (isnep) {
175: DSGetDimensions(ds,NULL,NULL,NULL,&nev);
176: for (i=0;i<nev;i++) inside[i] = i;
177: } else {
178: 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: PetscMalloc2(ld,&y,ld,&r);
185: #if !defined(PETSC_USE_COMPLEX)
186: PetscMalloc2(ld,&yi,ld,&ri);
187: #endif
188: DSVectors(ds,DS_MAT_X,NULL,NULL);
189: DSGetArray(ds,DS_MAT_X,&X);
190: 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: PetscArrayzero(r,n);
200: #if !defined(PETSC_USE_COMPLEX)
201: PetscArrayzero(ri,n);
202: #endif
203: /* Residual */
204: alpha = 1.0;
205: for (k=0;k<NMAT;k++) {
206: 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: DSRestoreArray(ds,mat[k],&A);
218: if (isnep) 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) 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) PetscViewerASCIIPrintf(viewer," %.5f\n",(double)re);
245: else 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) PetscViewerASCIIPrintf(viewer," %.5f\n",(double)re);
249: else PetscViewerASCIIPrintf(viewer," %.5f%+.5fi\n",(double)re,(double)-im);
250: #endif
251: }
252: DSRestoreArray(ds,DS_MAT_X,&X);
253: PetscFree3(wr,wi,inside);
254: PetscFree2(y,r);
255: #if !defined(PETSC_USE_COMPLEX)
256: PetscFree2(yi,ri);
257: #endif
258: if (isnep) {
259: for (i=0;i<NMAT;i++) FNDestroy(&f[i]);
260: }
261: DSDestroy(&ds);
262: RGDestroy(&rg);
263: 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*/