Actual source code: ex27.c
slepc-3.7.4 2017-05-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2013, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Simple nonlinear eigenproblem using the NLEIGS solver.\n\n"
23: "The command line options are:\n"
24: " -n <n>, where <n> = matrix dimension.\n"
25: " -split <0/1>, to select the split form in the problem definition (enabled by default)\n";
28: /*
29: Solve T(lambda)x=0 using NLEIGS solver
30: with T(lambda) = -D+sqrt(lambda)*I
31: where D is the Laplacian operator in 1 dimension
32: and with the interpolation interval [.01,16]
33: */
35: #include <slepcnep.h>
37: /*
38: User-defined routines
39: */
40: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
41: PetscErrorCode ComputeSingularities(NEP,PetscInt*,PetscScalar*,void*);
45: int main(int argc,char **argv)
46: {
47: NEP nep; /* nonlinear eigensolver context */
48: Mat F,A[2];
49: NEPType type;
50: PetscInt n=100,nev,Istart,Iend,i;
52: PetscBool split=PETSC_TRUE;
53: RG rg;
54: FN f[2];
55: PetscBool terse;
56: PetscScalar coeffs;
58: SlepcInitialize(&argc,&argv,(char*)0,help);
59: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
60: PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL);
61: PetscPrintf(PETSC_COMM_WORLD,"\nSquare root eigenproblem, n=%D%s\n\n",n,split?" (in split form)":"");
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Create nonlinear eigensolver context
65: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: NEPCreate(PETSC_COMM_WORLD,&nep);
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Select the NLEIGS solver and set required options for it
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: NEPSetType(nep,NEPNLEIGS);
74: NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,NULL);
75: NEPGetRG(nep,&rg);
76: RGSetType(rg,RGINTERVAL);
77: #if defined(PETSC_USE_COMPLEX)
78: RGIntervalSetEndpoints(rg,0.01,16.0,-0.001,0.001);
79: #else
80: RGIntervalSetEndpoints(rg,0.01,16.0,0,0);
81: #endif
82: NEPSetTarget(nep,1.1);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Define the nonlinear problem
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87:
88: if (split) {
89: /*
90: Create matrices for the split form
91: */
92: MatCreate(PETSC_COMM_WORLD,&A[0]);
93: MatSetSizes(A[0],PETSC_DECIDE,PETSC_DECIDE,n,n);
94: MatSetFromOptions(A[0]);
95: MatSetUp(A[0]);
96: MatGetOwnershipRange(A[0],&Istart,&Iend);
97: for (i=Istart;i<Iend;i++) {
98: if (i>0) { MatSetValue(A[0],i,i-1,1.0,INSERT_VALUES); }
99: if (i<n-1) { MatSetValue(A[0],i,i+1,1.0,INSERT_VALUES); }
100: MatSetValue(A[0],i,i,-2.0,INSERT_VALUES);
101: }
102: MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY);
103: MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);
105: MatCreate(PETSC_COMM_WORLD,&A[1]);
106: MatSetSizes(A[1],PETSC_DECIDE,PETSC_DECIDE,n,n);
107: MatSetFromOptions(A[1]);
108: MatSetUp(A[1]);
109: MatAssemblyBegin(A[1],MAT_FINAL_ASSEMBLY);
110: MatAssemblyEnd(A[1],MAT_FINAL_ASSEMBLY);
111: MatShift(A[1],1.0);
113: /*
114: Define funcions for the split form
115: */
116: FNCreate(PETSC_COMM_WORLD,&f[0]);
117: FNSetType(f[0],FNRATIONAL);
118: coeffs = 1.0;
119: FNRationalSetNumerator(f[0],1,&coeffs);
120: FNCreate(PETSC_COMM_WORLD,&f[1]);
121: FNSetType(f[1],FNSQRT);
122: NEPSetSplitOperator(nep,2,A,f,SUBSET_NONZERO_PATTERN);
124: } else {
125: /*
126: Callback form: create matrix and set Function evaluation routine
127: */
128: MatCreate(PETSC_COMM_WORLD,&F);
129: MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n);
130: MatSetFromOptions(F);
131: MatSeqAIJSetPreallocation(F,3,NULL);
132: MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
133: MatSetUp(F);
134: NEPSetFunction(nep,F,F,FormFunction,NULL);
135: }
137: /*
138: Set solver parameters at runtime
139: */
140: NEPSetFromOptions(nep);
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Solve the eigensystem
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: NEPSolve(nep);
146: NEPGetType(nep,&type);
147: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
148: NEPGetDimensions(nep,&nev,NULL,NULL);
149: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Display solution and clean up
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155: /* show detailed info unless -terse option is given by user */
156: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
157: if (terse) {
158: NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
159: } else {
160: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
161: NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
162: NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
163: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
164: }
165: NEPDestroy(&nep);
166: if (split) {
167: MatDestroy(&A[0]);
168: MatDestroy(&A[1]);
169: FNDestroy(&f[0]);
170: FNDestroy(&f[1]);
171: } else {
172: MatDestroy(&F);
173: }
174: SlepcFinalize();
175: return ierr;
176: }
178: /* ------------------------------------------------------------------- */
181: /*
182: FormFunction - Computes Function matrix T(lambda)
183: */
184: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
185: {
187: PetscInt i,n,col[3],Istart,Iend;
188: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
189: PetscScalar value[3],t;
192: /*
193: Compute Function entries and insert into matrix
194: */
195: t = PetscSqrtScalar(lambda);
196: MatGetSize(fun,&n,NULL);
197: MatGetOwnershipRange(fun,&Istart,&Iend);
198: if (Istart==0) FirstBlock=PETSC_TRUE;
199: if (Iend==n) LastBlock=PETSC_TRUE;
200: value[0]=1.0; value[1]=t-2.0; value[2]=1.0;
201: for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
202: col[0]=i-1; col[1]=i; col[2]=i+1;
203: MatSetValues(fun,1,&i,3,col,value,INSERT_VALUES);
204: }
205: if (LastBlock) {
206: i=n-1; col[0]=n-2; col[1]=n-1;
207: MatSetValues(fun,1,&i,2,col,value,INSERT_VALUES);
208: }
209: if (FirstBlock) {
210: i=0; col[0]=0; col[1]=1; value[0]=t-2.0; value[1]=1.0;
211: MatSetValues(fun,1,&i,2,col,value,INSERT_VALUES);
212: }
214: /*
215: Assemble matrix
216: */
217: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
218: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
219: if (fun != B) {
220: MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
221: MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
222: }
223: return(0);
224: }
228: /*
229: ComputeSingularities - Computes maxnp points (at most) in the complex plane where
230: the function T(.) is not analytic.
232: In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
233: */
234: PetscErrorCode ComputeSingularities(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *pt)
235: {
236: PetscReal h;
237: PetscInt i;
240: h = 12.0/(*maxnp-1);
241: xi[0] = -1e-6; xi[*maxnp-1] = -1e+6;
242: for (i=1;i<*maxnp-1;i++) xi[i] = -PetscPowReal(10,-6+h*i);
243: return(0);
244: }