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 interface functions of DSNEP.\n\n";
12 :
13 : #include <slepcds.h>
14 :
15 1 : int main(int argc,char **argv)
16 : {
17 1 : DS ds;
18 1 : FN f1,f2,f3,funs[3];
19 1 : SlepcSC sc;
20 1 : PetscScalar *Id,*A,*B,*wr,*wi,*X,*W,coeffs[2],auxr,alpha;
21 1 : PetscReal tau=0.001,h,a=20,xi,re,im,nrm,aux;
22 1 : PetscInt i,j,ii,jj,k,n=10,ld,nev,nfun,midx,ip,rits,meth,spls;
23 1 : PetscViewer viewer;
24 1 : PetscBool verbose;
25 1 : RG rg;
26 1 : DSMatType mat[3]={DS_MAT_E0,DS_MAT_E1,DS_MAT_E2};
27 : #if defined(PETSC_USE_COMPLEX)
28 1 : PetscBool flg;
29 : #else
30 : PetscScalar auxi;
31 : #endif
32 :
33 1 : PetscFunctionBeginUser;
34 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
35 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
36 1 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-tau",&tau,NULL));
37 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Solve a Dense System of type NEP - dimension %" PetscInt_FMT ", tau=%g.\n",n,(double)tau));
38 1 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
39 :
40 : /* Create DS object and set options */
41 1 : PetscCall(DSCreate(PETSC_COMM_WORLD,&ds));
42 1 : PetscCall(DSSetType(ds,DSNEP));
43 : #if defined(PETSC_USE_COMPLEX)
44 1 : PetscCall(DSSetMethod(ds,1)); /* contour integral */
45 : #endif
46 1 : PetscCall(DSNEPGetRG(ds,&rg));
47 1 : PetscCall(RGSetType(rg,RGELLIPSE));
48 1 : PetscCall(DSNEPSetMinimality(ds,1));
49 1 : PetscCall(DSNEPSetIntegrationPoints(ds,16));
50 1 : PetscCall(DSNEPSetRefine(ds,PETSC_CURRENT,2));
51 1 : PetscCall(DSNEPSetSamplingSize(ds,25));
52 1 : PetscCall(DSSetFromOptions(ds));
53 :
54 : /* Print current options */
55 1 : PetscCall(DSGetMethod(ds,&meth));
56 : #if defined(PETSC_USE_COMPLEX)
57 1 : PetscCheck(meth==1,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"This example requires ds_method=1");
58 1 : PetscCall(RGIsTrivial(rg,&flg));
59 1 : PetscCheck(!flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must at least set the radius of the ellipse");
60 : #endif
61 :
62 1 : PetscCall(DSNEPGetMinimality(ds,&midx));
63 1 : PetscCall(DSNEPGetIntegrationPoints(ds,&ip));
64 1 : PetscCall(DSNEPGetRefine(ds,NULL,&rits));
65 1 : PetscCall(DSNEPGetSamplingSize(ds,&spls));
66 1 : if (meth==1) {
67 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Contour integral method with %" PetscInt_FMT " integration points, minimality index %" PetscInt_FMT ", and sampling size %" PetscInt_FMT "\n",ip,midx,spls));
68 1 : if (rits) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Doing %" PetscInt_FMT " iterations of Newton refinement\n",rits));
69 : }
70 :
71 : /* Set functions (prior to DSAllocate) */
72 1 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f1));
73 1 : PetscCall(FNSetType(f1,FNRATIONAL));
74 1 : coeffs[0] = -1.0; coeffs[1] = 0.0;
75 1 : PetscCall(FNRationalSetNumerator(f1,2,coeffs));
76 :
77 1 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f2));
78 1 : PetscCall(FNSetType(f2,FNRATIONAL));
79 1 : coeffs[0] = 1.0;
80 1 : PetscCall(FNRationalSetNumerator(f2,1,coeffs));
81 :
82 1 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f3));
83 1 : PetscCall(FNSetType(f3,FNEXP));
84 1 : PetscCall(FNSetScale(f3,-tau,1.0));
85 :
86 1 : funs[0] = f1;
87 1 : funs[1] = f2;
88 1 : funs[2] = f3;
89 1 : PetscCall(DSNEPSetFN(ds,3,funs));
90 :
91 : /* Set dimensions */
92 1 : ld = n;
93 1 : PetscCall(DSAllocate(ds,ld));
94 1 : PetscCall(DSSetDimensions(ds,n,0,0));
95 :
96 : /* Set up viewer */
97 1 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
98 1 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
99 1 : PetscCall(PetscViewerPopFormat(viewer));
100 1 : if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
101 :
102 : /* Fill matrices */
103 1 : PetscCall(DSGetArray(ds,DS_MAT_E0,&Id));
104 11 : for (i=0;i<n;i++) Id[i+i*ld]=1.0;
105 1 : PetscCall(DSRestoreArray(ds,DS_MAT_E0,&Id));
106 1 : h = PETSC_PI/(PetscReal)(n+1);
107 1 : PetscCall(DSGetArray(ds,DS_MAT_E1,&A));
108 11 : for (i=0;i<n;i++) A[i+i*ld]=-2.0/(h*h)+a;
109 10 : for (i=1;i<n;i++) {
110 9 : A[i+(i-1)*ld]=1.0/(h*h);
111 9 : A[(i-1)+i*ld]=1.0/(h*h);
112 : }
113 1 : PetscCall(DSRestoreArray(ds,DS_MAT_E1,&A));
114 1 : PetscCall(DSGetArray(ds,DS_MAT_E2,&B));
115 11 : for (i=0;i<n;i++) {
116 10 : xi = (i+1)*h;
117 10 : B[i+i*ld] = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
118 : }
119 1 : PetscCall(DSRestoreArray(ds,DS_MAT_E2,&B));
120 :
121 1 : if (verbose) {
122 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Initial - - - - - - - - -\n"));
123 0 : PetscCall(DSView(ds,viewer));
124 : }
125 :
126 : /* Solve */
127 1 : PetscCall(PetscCalloc2(n,&wr,n,&wi));
128 1 : PetscCall(DSGetSlepcSC(ds,&sc));
129 1 : sc->comparison = SlepcCompareLargestMagnitude;
130 1 : sc->comparisonctx = NULL;
131 1 : sc->map = NULL;
132 1 : sc->mapobj = NULL;
133 1 : PetscCall(DSSolve(ds,wr,wi));
134 1 : PetscCall(DSSort(ds,wr,wi,NULL,NULL,NULL));
135 :
136 1 : if (verbose) {
137 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"After solve - - - - - - - - -\n"));
138 0 : PetscCall(DSView(ds,viewer));
139 : }
140 1 : PetscCall(DSGetDimensions(ds,NULL,NULL,NULL,&nev));
141 :
142 : /* Print computed eigenvalues */
143 1 : PetscCall(DSNEPGetNumFN(ds,&nfun));
144 1 : PetscCall(PetscMalloc1(ld*ld,&W));
145 1 : PetscCall(DSVectors(ds,DS_MAT_X,NULL,NULL));
146 1 : PetscCall(DSGetArray(ds,DS_MAT_X,&X));
147 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed eigenvalues =\n"));
148 4 : for (i=0;i<nev;i++) {
149 : #if defined(PETSC_USE_COMPLEX)
150 3 : re = PetscRealPart(wr[i]);
151 3 : im = PetscImaginaryPart(wr[i]);
152 : #else
153 : re = wr[i];
154 : im = wi[i];
155 : #endif
156 : /* Residual */
157 3 : PetscCall(PetscArrayzero(W,ld*ld));
158 12 : for (k=0;k<nfun;k++) {
159 9 : PetscCall(FNEvaluateFunction(funs[k],wr[i],&alpha));
160 9 : PetscCall(DSGetArray(ds,mat[k],&A));
161 999 : for (jj=0;jj<n;jj++) for (ii=0;ii<n;ii++) W[jj*ld+ii] += alpha*A[jj*ld+ii];
162 9 : PetscCall(DSRestoreArray(ds,mat[k],&A));
163 : }
164 : nrm = 0.0;
165 33 : for (k=0;k<n;k++) {
166 : auxr = 0.0;
167 : #if !defined(PETSC_USE_COMPLEX)
168 : auxi = 0.0;
169 : #endif
170 330 : for (j=0;j<n;j++) {
171 300 : auxr += W[k+j*ld]*X[i*ld+j];
172 : #if !defined(PETSC_USE_COMPLEX)
173 : if (PetscAbs(wi[j])!=0.0) auxi += W[k+j*ld]*X[(i+1)*ld+j];
174 : #endif
175 : }
176 30 : aux = SlepcAbsEigenvalue(auxr,auxi);
177 30 : nrm += aux*aux;
178 : }
179 3 : nrm = PetscSqrtReal(nrm);
180 3 : if (nrm>1000*n*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the residual norm of the %" PetscInt_FMT "-th computed eigenpair %g\n",i,(double)nrm));
181 3 : if (PetscAbs(im)<1e-10) PetscCall(PetscViewerASCIIPrintf(viewer," %.5f\n",(double)re));
182 3 : else PetscCall(PetscViewerASCIIPrintf(viewer," %.5f%+.5fi\n",(double)re,(double)im));
183 : }
184 1 : PetscCall(PetscFree(W));
185 1 : PetscCall(DSRestoreArray(ds,DS_MAT_X,&X));
186 1 : PetscCall(DSRestoreArray(ds,DS_MAT_W,&W));
187 1 : PetscCall(PetscFree2(wr,wi));
188 1 : PetscCall(FNDestroy(&f1));
189 1 : PetscCall(FNDestroy(&f2));
190 1 : PetscCall(FNDestroy(&f3));
191 1 : PetscCall(DSDestroy(&ds));
192 1 : PetscCall(SlepcFinalize());
193 : return 0;
194 : }
195 :
196 : /*TEST
197 :
198 : testset:
199 : test:
200 : suffix: 1
201 : requires: !complex
202 : test:
203 : suffix: 2
204 : args: -ds_nep_rg_ellipse_radius 10
205 : filter: sed -e "s/[+-]0\.0*i//g" | sed -e "s/37411/37410/"
206 : requires: complex
207 :
208 : TEST*/
|
