Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : static char help[] = "Test for DSPEP and DSNEP.\n\n";
12 :
13 : #include <slepcds.h>
14 :
15 : #define NMAT 5
16 :
17 1 : int main(int argc,char **argv)
18 : {
19 1 : DS ds;
20 1 : FN f[NMAT],qfun;
21 1 : SlepcSC sc;
22 1 : PetscScalar *A,*wr,*wi,*X,*y,*r,numer[NMAT],alpha;
23 1 : PetscReal c[10] = { 0.6, 1.3, 1.3, 0.1, 0.1, 1.2, 1.0, 1.0, 1.2, 1.0 };
24 1 : PetscReal tol,radius=1.5,re,im,nrm;
25 1 : PetscInt i,j,ii,jj,II,k,m=3,n,ld,nev,nfun,d,*inside;
26 1 : PetscViewer viewer;
27 1 : PetscBool verbose,isnep=PETSC_FALSE;
28 1 : RG rg;
29 1 : 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 1 : PetscScalar *yi,*ri,alphai=0.0,t;
32 : #endif
33 :
34 1 : PetscFunctionBeginUser;
35 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
36 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
37 1 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-isnep",&isnep,NULL));
38 1 : n = m*m;
39 1 : k = 10;
40 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nButterfly problem, n=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",n,m));
41 1 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
42 1 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-radius",&radius,NULL));
43 :
44 : /* Create DS object */
45 1 : PetscCall(DSCreate(PETSC_COMM_WORLD,&ds));
46 1 : tol = 1000*n*PETSC_MACHINE_EPSILON;
47 1 : if (isnep) {
48 0 : PetscCall(DSSetType(ds,DSNEP));
49 0 : PetscCall(DSSetMethod(ds,1));
50 0 : PetscCall(DSNEPSetRefine(ds,tol,PETSC_DECIDE));
51 1 : } else PetscCall(DSSetType(ds,DSPEP));
52 1 : PetscCall(DSSetFromOptions(ds));
53 :
54 : /* Set functions (prior to DSAllocate) f_i=x^i */
55 1 : if (isnep) {
56 0 : numer[0] = 1.0;
57 0 : for (j=1;j<NMAT;j++) numer[j] = 0.0;
58 0 : for (i=0;i<NMAT;i++) {
59 0 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f[i]));
60 0 : PetscCall(FNSetType(f[i],FNRATIONAL));
61 0 : PetscCall(FNRationalSetNumerator(f[i],i+1,numer));
62 : }
63 0 : PetscCall(DSNEPSetFN(ds,NMAT,f));
64 1 : } else PetscCall(DSPEPSetDegree(ds,NMAT-1));
65 :
66 : /* Set dimensions */
67 1 : ld = n+2; /* test leading dimension larger than n */
68 1 : PetscCall(DSAllocate(ds,ld));
69 1 : PetscCall(DSSetDimensions(ds,n,0,0));
70 :
71 : /* Set region (used only in method=1) */
72 1 : PetscCall(RGCreate(PETSC_COMM_WORLD,&rg));
73 1 : PetscCall(RGSetType(rg,RGELLIPSE));
74 1 : PetscCall(RGEllipseSetParameters(rg,1.5,radius,.5));
75 1 : PetscCall(RGSetFromOptions(rg));
76 1 : if (isnep) PetscCall(DSNEPSetRG(ds,rg));
77 :
78 : /* Set up viewer */
79 1 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
80 1 : PetscCall(DSViewFromOptions(ds,NULL,"-ds_view"));
81 1 : if (verbose) {
82 0 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
83 : /* Show info about functions */
84 0 : if (isnep) {
85 0 : PetscCall(DSNEPGetNumFN(ds,&nfun));
86 0 : for (i=0;i<nfun;i++) {
87 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Function %" PetscInt_FMT ":\n",i));
88 0 : PetscCall(DSNEPGetFN(ds,i,&qfun));
89 0 : PetscCall(FNView(qfun,NULL));
90 : }
91 : }
92 : }
93 :
94 : /* Fill matrices */
95 : /* A0 */
96 1 : PetscCall(DSGetArray(ds,DS_MAT_E0,&A));
97 10 : for (II=0;II<n;II++) {
98 9 : i = II/m; j = II-i*m;
99 9 : A[II+II*ld] = 4.0*c[0]/6.0+4.0*c[1]/6.0;
100 9 : if (j>0) A[II+(II-1)*ld] = c[0]/6.0;
101 9 : if (j<m-1) A[II+ld*(II+1)] = c[0]/6.0;
102 9 : if (i>0) A[II+ld*(II-m)] = c[1]/6.0;
103 9 : if (i<m-1) A[II+ld*(II+m)] = c[1]/6.0;
104 : }
105 1 : PetscCall(DSRestoreArray(ds,DS_MAT_E0,&A));
106 :
107 : /* A1 */
108 1 : PetscCall(DSGetArray(ds,DS_MAT_E1,&A));
109 10 : for (II=0;II<n;II++) {
110 9 : i = II/m; j = II-i*m;
111 9 : if (j>0) A[II+ld*(II-1)] = c[2];
112 9 : if (j<m-1) A[II+ld*(II+1)] = -c[2];
113 9 : if (i>0) A[II+ld*(II-m)] = c[3];
114 9 : if (i<m-1) A[II+ld*(II+m)] = -c[3];
115 : }
116 1 : PetscCall(DSRestoreArray(ds,DS_MAT_E1,&A));
117 :
118 : /* A2 */
119 1 : PetscCall(DSGetArray(ds,DS_MAT_E2,&A));
120 10 : for (II=0;II<n;II++) {
121 9 : i = II/m; j = II-i*m;
122 9 : A[II+ld*II] = -2.0*c[4]-2.0*c[5];
123 9 : if (j>0) A[II+ld*(II-1)] = c[4];
124 9 : if (j<m-1) A[II+ld*(II+1)] = c[4];
125 9 : if (i>0) A[II+ld*(II-m)] = c[5];
126 9 : if (i<m-1) A[II+ld*(II+m)] = c[5];
127 : }
128 1 : PetscCall(DSRestoreArray(ds,DS_MAT_E2,&A));
129 :
130 : /* A3 */
131 1 : PetscCall(DSGetArray(ds,DS_MAT_E3,&A));
132 10 : for (II=0;II<n;II++) {
133 9 : i = II/m; j = II-i*m;
134 9 : if (j>0) A[II+ld*(II-1)] = c[6];
135 9 : if (j<m-1) A[II+ld*(II+1)] = -c[6];
136 9 : if (i>0) A[II+ld*(II-m)] = c[7];
137 9 : if (i<m-1) A[II+ld*(II+m)] = -c[7];
138 : }
139 1 : PetscCall(DSRestoreArray(ds,DS_MAT_E3,&A));
140 :
141 : /* A4 */
142 1 : PetscCall(DSGetArray(ds,DS_MAT_E4,&A));
143 10 : for (II=0;II<n;II++) {
144 9 : i = II/m; j = II-i*m;
145 9 : A[II+ld*II] = 2.0*c[8]+2.0*c[9];
146 9 : if (j>0) A[II+ld*(II-1)] = -c[8];
147 9 : if (j<m-1) A[II+ld*(II+1)] = -c[8];
148 9 : if (i>0) A[II+ld*(II-m)] = -c[9];
149 9 : if (i<m-1) A[II+ld*(II+m)] = -c[9];
150 : }
151 1 : PetscCall(DSRestoreArray(ds,DS_MAT_E4,&A));
152 :
153 1 : if (verbose) {
154 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Initial - - - - - - - - -\n"));
155 0 : PetscCall(DSView(ds,viewer));
156 : }
157 :
158 : /* Solve */
159 1 : if (isnep) PetscCall(DSNEPGetMinimality(ds,&d));
160 1 : else PetscCall(DSPEPGetDegree(ds,&d));
161 1 : PetscCall(PetscCalloc3(n*d,&wr,n*d,&wi,n*d,&inside));
162 1 : PetscCall(DSGetSlepcSC(ds,&sc));
163 1 : sc->comparison = SlepcCompareLargestMagnitude;
164 1 : sc->comparisonctx = NULL;
165 1 : sc->map = NULL;
166 1 : sc->mapobj = NULL;
167 1 : PetscCall(DSSolve(ds,wr,wi));
168 1 : PetscCall(DSSort(ds,wr,wi,NULL,NULL,NULL));
169 :
170 1 : if (verbose) {
171 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"After solve - - - - - - - - -\n"));
172 0 : PetscCall(DSView(ds,viewer));
173 : }
174 1 : if (isnep) {
175 0 : PetscCall(DSGetDimensions(ds,NULL,NULL,NULL,&nev));
176 0 : for (i=0;i<nev;i++) inside[i] = i;
177 : } else {
178 1 : PetscCall(RGCheckInside(rg,d*n,wr,wi,inside));
179 1 : nev = 0;
180 37 : for (i=0;i<d*n;i++) if (inside[i]>0) inside[nev++] = i;
181 : }
182 :
183 : /* Print computed eigenvalues */
184 1 : PetscCall(PetscMalloc2(ld,&y,ld,&r));
185 : #if !defined(PETSC_USE_COMPLEX)
186 1 : PetscCall(PetscMalloc2(ld,&yi,ld,&ri));
187 : #endif
188 1 : PetscCall(DSVectors(ds,DS_MAT_X,NULL,NULL));
189 1 : PetscCall(DSGetArray(ds,DS_MAT_X,&X));
190 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed eigenvalues in the region: %" PetscInt_FMT "\n",nev));
191 9 : 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 8 : re = wr[inside[i]];
197 8 : im = wi[inside[i]];
198 : #endif
199 8 : PetscCall(PetscArrayzero(r,n));
200 : #if !defined(PETSC_USE_COMPLEX)
201 8 : PetscCall(PetscArrayzero(ri,n));
202 : #endif
203 : /* Residual */
204 8 : alpha = 1.0;
205 48 : for (k=0;k<NMAT;k++) {
206 40 : PetscCall(DSGetArray(ds,mat[k],&A));
207 400 : for (ii=0;ii<n;ii++) {
208 360 : y[ii] = 0.0;
209 3600 : 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 400 : for (ii=0;ii<n;ii++) {
213 360 : yi[ii] = 0.0;
214 3600 : for (jj=0;jj<n;jj++) yi[ii] += A[jj*ld+ii]*X[inside[i+1]*ld+jj];
215 : }
216 : #endif
217 40 : PetscCall(DSRestoreArray(ds,mat[k],&A));
218 40 : if (isnep) PetscCall(FNEvaluateFunction(f[k],wr[inside[i]],&alpha));
219 400 : for (ii=0;ii<n;ii++) r[ii] += alpha*y[ii];
220 : #if !defined(PETSC_USE_COMPLEX)
221 400 : for (ii=0;ii<n;ii++) r[ii] -= alphai*yi[ii];
222 400 : for (ii=0;ii<n;ii++) ri[ii] += alpha*yi[ii]+alphai*y[ii];
223 : #endif
224 40 : if (!isnep) {
225 : #if defined(PETSC_USE_COMPLEX)
226 : alpha *= wr[inside[i]];
227 : #else
228 40 : t = alpha;
229 40 : alpha = alpha*re-alphai*im;
230 40 : alphai = alphai*re+t*im;
231 : #endif
232 : }
233 : }
234 : nrm = 0.0;
235 80 : for (k=0;k<n;k++) {
236 : #if !defined(PETSC_USE_COMPLEX)
237 72 : nrm += r[k]*r[k]+ri[k]*ri[k];
238 : #else
239 : nrm += PetscRealPart(r[k]*PetscConj(r[k]));
240 : #endif
241 : }
242 8 : nrm = PetscSqrtReal(nrm);
243 8 : 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 8 : if (PetscAbs(im)<1e-10) PetscCall(PetscViewerASCIIPrintf(viewer," %.5f\n",(double)re));
245 8 : else PetscCall(PetscViewerASCIIPrintf(viewer," %.5f%+.5fi\n",(double)re,(double)im));
246 : #if !defined(PETSC_USE_COMPLEX)
247 8 : if (im!=0.0) i++;
248 8 : if (PetscAbs(im)<1e-10) PetscCall(PetscViewerASCIIPrintf(viewer," %.5f\n",(double)re));
249 8 : else PetscCall(PetscViewerASCIIPrintf(viewer," %.5f%+.5fi\n",(double)re,(double)-im));
250 : #endif
251 : }
252 1 : PetscCall(DSRestoreArray(ds,DS_MAT_X,&X));
253 1 : PetscCall(PetscFree3(wr,wi,inside));
254 1 : PetscCall(PetscFree2(y,r));
255 : #if !defined(PETSC_USE_COMPLEX)
256 1 : PetscCall(PetscFree2(yi,ri));
257 : #endif
258 1 : if (isnep) {
259 0 : for (i=0;i<NMAT;i++) PetscCall(FNDestroy(&f[i]));
260 : }
261 1 : PetscCall(DSDestroy(&ds));
262 1 : PetscCall(RGDestroy(&rg));
263 1 : PetscCall(SlepcFinalize());
264 : return 0;
265 : }
266 :
267 : /*TEST
268 :
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
278 :
279 : TEST*/
|
