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 the NLEIGS solver with shell matrices.\n\n"
12 : "This is based on ex27.\n"
13 : "The command line options are:\n"
14 : " -n <n>, where <n> = matrix dimension.\n"
15 : " -split <0/1>, to select the split form in the problem definition (enabled by default).\n";
16 :
17 : /*
18 : Solve T(lambda)x=0 using NLEIGS solver
19 : with T(lambda) = -D+sqrt(lambda)*I
20 : where D is the Laplacian operator in 1 dimension
21 : and with the interpolation interval [.01,16]
22 : */
23 :
24 : #include <slepcnep.h>
25 :
26 : /* User-defined routines */
27 : PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
28 : PetscErrorCode ComputeSingularities(NEP,PetscInt*,PetscScalar*,void*);
29 : PetscErrorCode MatMult_A0(Mat,Vec,Vec);
30 : PetscErrorCode MatGetDiagonal_A0(Mat,Vec);
31 : PetscErrorCode MatDuplicate_A0(Mat,MatDuplicateOption,Mat*);
32 : PetscErrorCode MatMult_A1(Mat,Vec,Vec);
33 : PetscErrorCode MatGetDiagonal_A1(Mat,Vec);
34 : PetscErrorCode MatDuplicate_A1(Mat,MatDuplicateOption,Mat*);
35 : PetscErrorCode MatMult_F(Mat,Vec,Vec);
36 : PetscErrorCode MatGetDiagonal_F(Mat,Vec);
37 : PetscErrorCode MatDuplicate_F(Mat,MatDuplicateOption,Mat*);
38 : PetscErrorCode MatDestroy_F(Mat);
39 :
40 : typedef struct {
41 : PetscScalar t; /* square root of lambda */
42 : } MatCtx;
43 :
44 6 : int main(int argc,char **argv)
45 : {
46 6 : NEP nep;
47 6 : KSP *ksp;
48 6 : PC pc;
49 6 : Mat F,A[2];
50 6 : NEPType type;
51 6 : PetscInt i,n=100,nev,its,nsolve;
52 6 : PetscReal keep,tol=PETSC_SQRT_MACHINE_EPSILON/10;
53 6 : RG rg;
54 6 : FN f[2];
55 6 : PetscBool terse,flg,lock,split=PETSC_TRUE;
56 6 : PetscScalar coeffs;
57 6 : MatCtx *ctx;
58 :
59 6 : PetscFunctionBeginUser;
60 6 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
61 6 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
62 6 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL));
63 9 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSquare root eigenproblem, n=%" PetscInt_FMT "%s\n\n",n,split?" (in split form)":""));
64 :
65 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66 : Create NEP context, configure NLEIGS with appropriate options
67 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68 :
69 6 : PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
70 6 : PetscCall(NEPSetType(nep,NEPNLEIGS));
71 6 : PetscCall(NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,NULL));
72 6 : PetscCall(NEPGetRG(nep,&rg));
73 6 : PetscCall(RGSetType(rg,RGINTERVAL));
74 : #if defined(PETSC_USE_COMPLEX)
75 6 : PetscCall(RGIntervalSetEndpoints(rg,0.01,16.0,-0.001,0.001));
76 : #else
77 : PetscCall(RGIntervalSetEndpoints(rg,0.01,16.0,0,0));
78 : #endif
79 6 : PetscCall(NEPSetTarget(nep,1.1));
80 6 : PetscCall(NEPNLEIGSGetKSPs(nep,&nsolve,&ksp));
81 12 : for (i=0;i<nsolve;i++) {
82 6 : PetscCall(KSPSetType(ksp[i],KSPBICG));
83 6 : PetscCall(KSPGetPC(ksp[i],&pc));
84 6 : PetscCall(PCSetType(pc,PCJACOBI));
85 6 : PetscCall(KSPSetTolerances(ksp[i],tol,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
86 : }
87 :
88 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89 : Define the nonlinear problem
90 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91 :
92 6 : if (split) {
93 : /* Create matrix A0 (tridiagonal) */
94 3 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,NULL,&A[0]));
95 3 : PetscCall(MatShellSetOperation(A[0],MATOP_MULT,(void(*)(void))MatMult_A0));
96 3 : PetscCall(MatShellSetOperation(A[0],MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A0));
97 3 : PetscCall(MatShellSetOperation(A[0],MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A0));
98 3 : PetscCall(MatShellSetOperation(A[0],MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A0));
99 :
100 : /* Create matrix A0 (identity) */
101 3 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,NULL,&A[1]));
102 3 : PetscCall(MatShellSetOperation(A[1],MATOP_MULT,(void(*)(void))MatMult_A1));
103 3 : PetscCall(MatShellSetOperation(A[1],MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A1));
104 3 : PetscCall(MatShellSetOperation(A[1],MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A1));
105 3 : PetscCall(MatShellSetOperation(A[1],MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A1));
106 :
107 : /* Define functions for the split form */
108 3 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f[0]));
109 3 : PetscCall(FNSetType(f[0],FNRATIONAL));
110 3 : coeffs = 1.0;
111 3 : PetscCall(FNRationalSetNumerator(f[0],1,&coeffs));
112 3 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f[1]));
113 3 : PetscCall(FNSetType(f[1],FNSQRT));
114 3 : PetscCall(NEPSetSplitOperator(nep,2,A,f,SUBSET_NONZERO_PATTERN));
115 : } else {
116 : /* Callback form: create shell matrix for F=A0+sqrt(lambda)*A1 */
117 3 : PetscCall(PetscNew(&ctx));
118 3 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctx,&F));
119 3 : PetscCall(MatShellSetOperation(F,MATOP_MULT,(void(*)(void))MatMult_F));
120 3 : PetscCall(MatShellSetOperation(F,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_F));
121 3 : PetscCall(MatShellSetOperation(F,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_F));
122 3 : PetscCall(MatShellSetOperation(F,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_F));
123 3 : PetscCall(MatShellSetOperation(F,MATOP_DESTROY,(void(*)(void))MatDestroy_F));
124 : /* Set Function evaluation routine */
125 3 : PetscCall(NEPSetFunction(nep,F,F,FormFunction,NULL));
126 : }
127 :
128 : /* Set solver parameters at runtime */
129 6 : PetscCall(NEPSetFromOptions(nep));
130 :
131 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132 : Solve the eigensystem
133 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134 6 : PetscCall(NEPSolve(nep));
135 6 : PetscCall(NEPGetType(nep,&type));
136 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type));
137 6 : PetscCall(NEPGetDimensions(nep,&nev,NULL,NULL));
138 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
139 6 : PetscCall(PetscObjectTypeCompare((PetscObject)nep,NEPNLEIGS,&flg));
140 6 : if (flg) {
141 6 : PetscCall(NEPNLEIGSGetRestart(nep,&keep));
142 6 : PetscCall(NEPNLEIGSGetLocking(nep,&lock));
143 6 : PetscCall(NEPNLEIGSGetInterpolation(nep,&tol,&its));
144 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Restart factor is %3.2f",(double)keep));
145 6 : if (lock) PetscCall(PetscPrintf(PETSC_COMM_WORLD," (locking activated)"));
146 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n Divided diferences with tol=%6.2g maxit=%" PetscInt_FMT "\n",(double)tol,its));
147 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
148 : }
149 :
150 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151 : Display solution and clean up
152 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153 :
154 : /* show detailed info unless -terse option is given by user */
155 6 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
156 6 : if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
157 : else {
158 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
159 0 : PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
160 0 : PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
161 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
162 : }
163 6 : PetscCall(NEPDestroy(&nep));
164 6 : if (split) {
165 3 : PetscCall(MatDestroy(&A[0]));
166 3 : PetscCall(MatDestroy(&A[1]));
167 3 : PetscCall(FNDestroy(&f[0]));
168 3 : PetscCall(FNDestroy(&f[1]));
169 3 : } else PetscCall(MatDestroy(&F));
170 6 : PetscCall(SlepcFinalize());
171 : return 0;
172 : }
173 :
174 : /*
175 : FormFunction - Computes Function matrix T(lambda)
176 : */
177 75 : PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
178 : {
179 75 : MatCtx *ctxF;
180 :
181 75 : PetscFunctionBeginUser;
182 75 : PetscCall(MatShellGetContext(fun,&ctxF));
183 75 : ctxF->t = PetscSqrtScalar(lambda);
184 75 : PetscFunctionReturn(PETSC_SUCCESS);
185 : }
186 :
187 : /*
188 : ComputeSingularities - Computes maxnp points (at most) in the complex plane where
189 : the function T(.) is not analytic.
190 :
191 : In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
192 : */
193 6 : PetscErrorCode ComputeSingularities(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *pt)
194 : {
195 6 : PetscReal h;
196 6 : PetscInt i;
197 :
198 6 : PetscFunctionBeginUser;
199 6 : h = 12.0/(*maxnp-1);
200 6 : xi[0] = -1e-6; xi[*maxnp-1] = -1e+6;
201 59994 : for (i=1;i<*maxnp-1;i++) xi[i] = -PetscPowReal(10,-6+h*i);
202 6 : PetscFunctionReturn(PETSC_SUCCESS);
203 : }
204 :
205 : /* -------------------------------- A0 ----------------------------------- */
206 :
207 16931 : PetscErrorCode MatMult_A0(Mat A,Vec x,Vec y)
208 : {
209 16931 : PetscInt i,n;
210 16931 : PetscMPIInt rank,size,next,prev;
211 16931 : const PetscScalar *px;
212 16931 : PetscScalar *py,upper=0.0,lower=0.0;
213 16931 : MPI_Comm comm;
214 :
215 16931 : PetscFunctionBeginUser;
216 16931 : PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
217 16931 : PetscCallMPI(MPI_Comm_size(comm,&size));
218 16931 : PetscCallMPI(MPI_Comm_rank(comm,&rank));
219 16931 : next = rank==size-1? MPI_PROC_NULL: rank+1;
220 16931 : prev = rank==0? MPI_PROC_NULL: rank-1;
221 :
222 16931 : PetscCall(VecGetArrayRead(x,&px));
223 16931 : PetscCall(VecGetArray(y,&py));
224 16931 : PetscCall(VecGetLocalSize(x,&n));
225 :
226 16931 : PetscCallMPI(MPI_Sendrecv(px,1,MPIU_SCALAR,prev,0,&lower,1,MPIU_SCALAR,next,0,comm,MPI_STATUS_IGNORE));
227 16931 : PetscCallMPI(MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,next,0,&upper,1,MPIU_SCALAR,prev,0,comm,MPI_STATUS_IGNORE));
228 :
229 16931 : py[0] = upper-2.0*px[0]+px[1];
230 1112069 : for (i=1;i<n-1;i++) py[i] = px[i-1]-2.0*px[i]+px[i+1];
231 16931 : py[n-1] = px[n-2]-2.0*px[n-1]+lower;
232 16931 : PetscCall(VecRestoreArrayRead(x,&px));
233 16931 : PetscCall(VecRestoreArray(y,&py));
234 16931 : PetscFunctionReturn(PETSC_SUCCESS);
235 : }
236 :
237 3 : PetscErrorCode MatGetDiagonal_A0(Mat A,Vec diag)
238 : {
239 3 : PetscFunctionBeginUser;
240 3 : PetscCall(VecSet(diag,-2.0));
241 3 : PetscFunctionReturn(PETSC_SUCCESS);
242 : }
243 :
244 9 : PetscErrorCode MatDuplicate_A0(Mat A,MatDuplicateOption op,Mat *B)
245 : {
246 9 : PetscInt m,n,M,N;
247 9 : MPI_Comm comm;
248 :
249 9 : PetscFunctionBegin;
250 9 : PetscCall(MatGetSize(A,&M,&N));
251 9 : PetscCall(MatGetLocalSize(A,&m,&n));
252 9 : PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
253 9 : PetscCall(MatCreateShell(comm,m,n,M,N,NULL,B));
254 9 : PetscCall(MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_A0));
255 9 : PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A0));
256 9 : PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A0));
257 9 : PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A0));
258 9 : PetscFunctionReturn(PETSC_SUCCESS);
259 : }
260 :
261 : /* -------------------------------- A1 ----------------------------------- */
262 :
263 16931 : PetscErrorCode MatMult_A1(Mat A,Vec x,Vec y)
264 : {
265 16931 : PetscFunctionBeginUser;
266 16931 : PetscCall(VecCopy(x,y));
267 16931 : PetscFunctionReturn(PETSC_SUCCESS);
268 : }
269 :
270 3 : PetscErrorCode MatGetDiagonal_A1(Mat A,Vec diag)
271 : {
272 3 : PetscFunctionBeginUser;
273 3 : PetscCall(VecSet(diag,1.0));
274 3 : PetscFunctionReturn(PETSC_SUCCESS);
275 : }
276 :
277 3 : PetscErrorCode MatDuplicate_A1(Mat A,MatDuplicateOption op,Mat *B)
278 : {
279 3 : PetscInt m,n,M,N;
280 3 : MPI_Comm comm;
281 :
282 3 : PetscFunctionBegin;
283 3 : PetscCall(MatGetSize(A,&M,&N));
284 3 : PetscCall(MatGetLocalSize(A,&m,&n));
285 3 : PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
286 3 : PetscCall(MatCreateShell(comm,m,n,M,N,NULL,B));
287 3 : PetscCall(MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_A1));
288 3 : PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A1));
289 3 : PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A1));
290 3 : PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A1));
291 3 : PetscFunctionReturn(PETSC_SUCCESS);
292 : }
293 :
294 : /* -------------------------------- F ----------------------------------- */
295 :
296 370908 : PetscErrorCode MatMult_F(Mat A,Vec x,Vec y)
297 : {
298 370908 : PetscInt i,n;
299 370908 : PetscMPIInt rank,size,next,prev;
300 370908 : const PetscScalar *px;
301 370908 : PetscScalar *py,d,upper=0.0,lower=0.0;
302 370908 : MatCtx *ctx;
303 370908 : MPI_Comm comm;
304 :
305 370908 : PetscFunctionBeginUser;
306 370908 : PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
307 370908 : PetscCallMPI(MPI_Comm_size(comm,&size));
308 370908 : PetscCallMPI(MPI_Comm_rank(comm,&rank));
309 370908 : next = rank==size-1? MPI_PROC_NULL: rank+1;
310 370908 : prev = rank==0? MPI_PROC_NULL: rank-1;
311 :
312 370908 : PetscCall(MatShellGetContext(A,&ctx));
313 370908 : PetscCall(VecGetArrayRead(x,&px));
314 370908 : PetscCall(VecGetArray(y,&py));
315 370908 : PetscCall(VecGetLocalSize(x,&n));
316 :
317 370908 : PetscCallMPI(MPI_Sendrecv(px,1,MPIU_SCALAR,prev,0,&lower,1,MPIU_SCALAR,next,0,comm,MPI_STATUS_IGNORE));
318 370908 : PetscCallMPI(MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,next,0,&upper,1,MPIU_SCALAR,prev,0,comm,MPI_STATUS_IGNORE));
319 :
320 370908 : d = -2.0+ctx->t;
321 370908 : py[0] = upper+d*px[0]+px[1];
322 24624692 : for (i=1;i<n-1;i++) py[i] = px[i-1]+d*px[i]+px[i+1];
323 370908 : py[n-1] = px[n-2]+d*px[n-1]+lower;
324 370908 : PetscCall(VecRestoreArrayRead(x,&px));
325 370908 : PetscCall(VecRestoreArray(y,&py));
326 370908 : PetscFunctionReturn(PETSC_SUCCESS);
327 : }
328 :
329 66 : PetscErrorCode MatGetDiagonal_F(Mat A,Vec diag)
330 : {
331 66 : MatCtx *ctx;
332 :
333 66 : PetscFunctionBeginUser;
334 66 : PetscCall(MatShellGetContext(A,&ctx));
335 66 : PetscCall(VecSet(diag,-2.0+ctx->t));
336 66 : PetscFunctionReturn(PETSC_SUCCESS);
337 : }
338 :
339 66 : PetscErrorCode MatDuplicate_F(Mat A,MatDuplicateOption op,Mat *B)
340 : {
341 66 : MatCtx *actx,*bctx;
342 66 : PetscInt m,n,M,N;
343 66 : MPI_Comm comm;
344 :
345 66 : PetscFunctionBegin;
346 66 : PetscCall(MatShellGetContext(A,&actx));
347 66 : PetscCall(MatGetSize(A,&M,&N));
348 66 : PetscCall(MatGetLocalSize(A,&m,&n));
349 66 : PetscCall(PetscNew(&bctx));
350 66 : bctx->t = actx->t;
351 66 : PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
352 66 : PetscCall(MatCreateShell(comm,m,n,M,N,(void*)bctx,B));
353 66 : PetscCall(MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_F));
354 66 : PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_F));
355 66 : PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_F));
356 66 : PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_F));
357 66 : PetscCall(MatShellSetOperation(*B,MATOP_DESTROY,(void(*)(void))MatDestroy_F));
358 66 : PetscFunctionReturn(PETSC_SUCCESS);
359 : }
360 :
361 69 : PetscErrorCode MatDestroy_F(Mat A)
362 : {
363 69 : MatCtx *ctx;
364 :
365 69 : PetscFunctionBegin;
366 69 : PetscCall(MatShellGetContext(A,&ctx));
367 69 : PetscCall(PetscFree(ctx));
368 69 : PetscFunctionReturn(PETSC_SUCCESS);
369 : }
370 :
371 : /*TEST
372 :
373 : testset:
374 : nsize: {{1 2}}
375 : args: -nep_nev 3 -nep_tol 1e-8 -terse
376 : filter: sed -e "s/[+-]0\.0*i//g"
377 : requires: !single
378 : test:
379 : suffix: 1
380 : args: -nep_nleigs_locking 0 -nep_nleigs_interpolation_degree 90 -nep_nleigs_interpolation_tol 1e-8 -nep_nleigs_restart 0.4
381 : test:
382 : suffix: 2
383 : args: -split 0
384 :
385 : TEST*/
|
