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 2 : int main(int argc,char **argv)
18 : {
19 2 : DS ds;
20 2 : FN f[NMAT],qfun;
21 2 : SlepcSC sc;
22 2 : PetscScalar *A,*wr,*wi,*X,*y,*r,numer[NMAT],alpha;
23 2 : PetscReal c[10] = { 0.6, 1.3, 1.3, 0.1, 0.1, 1.2, 1.0, 1.0, 1.2, 1.0 };
24 2 : PetscReal tol,radius=1.5,re,im,nrm;
25 2 : PetscInt i,j,ii,jj,II,k,m=3,n,ld,nev,nfun,d,*inside;
26 2 : PetscViewer viewer;
27 2 : PetscBool verbose,isnep=PETSC_FALSE;
28 2 : RG rg;
29 2 : 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
33 :
34 2 : PetscFunctionBeginUser;
35 2 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
36 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
37 2 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-isnep",&isnep,NULL));
38 2 : n = m*m;
39 2 : k = 10;
40 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nButterfly problem, n=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",n,m));
41 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
42 2 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-radius",&radius,NULL));
43 :
44 : /* Create DS object */
45 2 : PetscCall(DSCreate(PETSC_COMM_WORLD,&ds));
46 2 : tol = 1000*n*PETSC_MACHINE_EPSILON;
47 2 : if (isnep) {
48 1 : PetscCall(DSSetType(ds,DSNEP));
49 1 : PetscCall(DSSetMethod(ds,1));
50 1 : PetscCall(DSNEPSetRefine(ds,tol,PETSC_DECIDE));
51 1 : } else PetscCall(DSSetType(ds,DSPEP));
52 2 : PetscCall(DSSetFromOptions(ds));
53 :
54 : /* Set functions (prior to DSAllocate) f_i=x^i */
55 2 : if (isnep) {
56 1 : numer[0] = 1.0;
57 5 : for (j=1;j<NMAT;j++) numer[j] = 0.0;
58 6 : for (i=0;i<NMAT;i++) {
59 5 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f[i]));
60 5 : PetscCall(FNSetType(f[i],FNRATIONAL));
61 5 : PetscCall(FNRationalSetNumerator(f[i],i+1,numer));
62 : }
63 1 : PetscCall(DSNEPSetFN(ds,NMAT,f));
64 1 : } else PetscCall(DSPEPSetDegree(ds,NMAT-1));
65 :
66 : /* Set dimensions */
67 2 : ld = n+2; /* test leading dimension larger than n */
68 2 : PetscCall(DSAllocate(ds,ld));
69 2 : PetscCall(DSSetDimensions(ds,n,0,0));
70 :
71 : /* Set region (used only in method=1) */
72 2 : PetscCall(RGCreate(PETSC_COMM_WORLD,&rg));
73 2 : PetscCall(RGSetType(rg,RGELLIPSE));
74 2 : PetscCall(RGEllipseSetParameters(rg,1.5,radius,.5));
75 2 : PetscCall(RGSetFromOptions(rg));
76 2 : if (isnep) PetscCall(DSNEPSetRG(ds,rg));
77 :
78 : /* Set up viewer */
79 2 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
80 2 : PetscCall(DSViewFromOptions(ds,NULL,"-ds_view"));
81 2 : 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 2 : PetscCall(DSGetArray(ds,DS_MAT_E0,&A));
97 20 : for (II=0;II<n;II++) {
98 18 : i = II/m; j = II-i*m;
99 18 : A[II+II*ld] = 4.0*c[0]/6.0+4.0*c[1]/6.0;
100 18 : if (j>0) A[II+(II-1)*ld] = c[0]/6.0;
101 18 : if (j<m-1) A[II+ld*(II+1)] = c[0]/6.0;
102 18 : if (i>0) A[II+ld*(II-m)] = c[1]/6.0;
103 18 : if (i<m-1) A[II+ld*(II+m)] = c[1]/6.0;
104 : }
105 2 : PetscCall(DSRestoreArray(ds,DS_MAT_E0,&A));
106 :
107 : /* A1 */
108 2 : PetscCall(DSGetArray(ds,DS_MAT_E1,&A));
109 20 : for (II=0;II<n;II++) {
110 18 : i = II/m; j = II-i*m;
111 18 : if (j>0) A[II+ld*(II-1)] = c[2];
112 18 : if (j<m-1) A[II+ld*(II+1)] = -c[2];
113 18 : if (i>0) A[II+ld*(II-m)] = c[3];
114 18 : if (i<m-1) A[II+ld*(II+m)] = -c[3];
115 : }
116 2 : PetscCall(DSRestoreArray(ds,DS_MAT_E1,&A));
117 :
118 : /* A2 */
119 2 : PetscCall(DSGetArray(ds,DS_MAT_E2,&A));
120 20 : for (II=0;II<n;II++) {
121 18 : i = II/m; j = II-i*m;
122 18 : A[II+ld*II] = -2.0*c[4]-2.0*c[5];
123 18 : if (j>0) A[II+ld*(II-1)] = c[4];
124 18 : if (j<m-1) A[II+ld*(II+1)] = c[4];
125 18 : if (i>0) A[II+ld*(II-m)] = c[5];
126 18 : if (i<m-1) A[II+ld*(II+m)] = c[5];
127 : }
128 2 : PetscCall(DSRestoreArray(ds,DS_MAT_E2,&A));
129 :
130 : /* A3 */
131 2 : PetscCall(DSGetArray(ds,DS_MAT_E3,&A));
132 20 : for (II=0;II<n;II++) {
133 18 : i = II/m; j = II-i*m;
134 18 : if (j>0) A[II+ld*(II-1)] = c[6];
135 18 : if (j<m-1) A[II+ld*(II+1)] = -c[6];
136 18 : if (i>0) A[II+ld*(II-m)] = c[7];
137 18 : if (i<m-1) A[II+ld*(II+m)] = -c[7];
138 : }
139 2 : PetscCall(DSRestoreArray(ds,DS_MAT_E3,&A));
140 :
141 : /* A4 */
142 2 : PetscCall(DSGetArray(ds,DS_MAT_E4,&A));
143 20 : for (II=0;II<n;II++) {
144 18 : i = II/m; j = II-i*m;
145 18 : A[II+ld*II] = 2.0*c[8]+2.0*c[9];
146 18 : if (j>0) A[II+ld*(II-1)] = -c[8];
147 18 : if (j<m-1) A[II+ld*(II+1)] = -c[8];
148 18 : if (i>0) A[II+ld*(II-m)] = -c[9];
149 18 : if (i<m-1) A[II+ld*(II+m)] = -c[9];
150 : }
151 2 : PetscCall(DSRestoreArray(ds,DS_MAT_E4,&A));
152 :
153 2 : if (verbose) {
154 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Initial - - - - - - - - -\n"));
155 0 : PetscCall(DSView(ds,viewer));
156 : }
157 :
158 : /* Solve */
159 2 : if (isnep) PetscCall(DSNEPGetMinimality(ds,&d));
160 1 : else PetscCall(DSPEPGetDegree(ds,&d));
161 2 : PetscCall(PetscCalloc3(n*d,&wr,n*d,&wi,n*d,&inside));
162 2 : PetscCall(DSGetSlepcSC(ds,&sc));
163 2 : sc->comparison = SlepcCompareLargestMagnitude;
164 2 : sc->comparisonctx = NULL;
165 2 : sc->map = NULL;
166 2 : sc->mapobj = NULL;
167 2 : PetscCall(DSSolve(ds,wr,wi));
168 2 : PetscCall(DSSort(ds,wr,wi,NULL,NULL,NULL));
169 :
170 2 : if (verbose) {
171 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"After solve - - - - - - - - -\n"));
172 0 : PetscCall(DSView(ds,viewer));
173 : }
174 2 : if (isnep) {
175 1 : PetscCall(DSGetDimensions(ds,NULL,NULL,NULL,&nev));
176 17 : 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 2 : PetscCall(PetscMalloc2(ld,&y,ld,&r));
185 : #if !defined(PETSC_USE_COMPLEX)
186 : PetscCall(PetscMalloc2(ld,&yi,ld,&ri));
187 : #endif
188 2 : PetscCall(DSVectors(ds,DS_MAT_X,NULL,NULL));
189 2 : PetscCall(DSGetArray(ds,DS_MAT_X,&X));
190 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed eigenvalues in the region: %" PetscInt_FMT "\n",nev));
191 34 : for (i=0;i<nev;i++) {
192 : #if defined(PETSC_USE_COMPLEX)
193 32 : re = PetscRealPart(wr[inside[i]]);
194 32 : im = PetscImaginaryPart(wr[inside[i]]);
195 : #else
196 : re = wr[inside[i]];
197 : im = wi[inside[i]];
198 : #endif
199 32 : PetscCall(PetscArrayzero(r,n));
200 : #if !defined(PETSC_USE_COMPLEX)
201 : PetscCall(PetscArrayzero(ri,n));
202 : #endif
203 : /* Residual */
204 32 : alpha = 1.0;
205 192 : for (k=0;k<NMAT;k++) {
206 160 : PetscCall(DSGetArray(ds,mat[k],&A));
207 1600 : for (ii=0;ii<n;ii++) {
208 1440 : y[ii] = 0.0;
209 14400 : 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 160 : PetscCall(DSRestoreArray(ds,mat[k],&A));
218 160 : if (isnep) PetscCall(FNEvaluateFunction(f[k],wr[inside[i]],&alpha));
219 1600 : 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 160 : if (!isnep) {
225 : #if defined(PETSC_USE_COMPLEX)
226 80 : 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 320 : 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 288 : nrm += PetscRealPart(r[k]*PetscConj(r[k]));
240 : #endif
241 : }
242 32 : nrm = PetscSqrtReal(nrm);
243 32 : 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 32 : if (PetscAbs(im)<1e-10) PetscCall(PetscViewerASCIIPrintf(viewer," %.5f\n",(double)re));
245 32 : 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 2 : PetscCall(DSRestoreArray(ds,DS_MAT_X,&X));
253 2 : PetscCall(PetscFree3(wr,wi,inside));
254 2 : PetscCall(PetscFree2(y,r));
255 : #if !defined(PETSC_USE_COMPLEX)
256 : PetscCall(PetscFree2(yi,ri));
257 : #endif
258 2 : if (isnep) {
259 6 : for (i=0;i<NMAT;i++) PetscCall(FNDestroy(&f[i]));
260 : }
261 2 : PetscCall(DSDestroy(&ds));
262 2 : PetscCall(RGDestroy(&rg));
263 2 : 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*/
|
