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[] = "User-defined split preconditioner when solving a quadratic eigenproblem.\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
14 : " -m <m>, where <m> = number of grid subdivisions in y dimension.\n\n";
15 :
16 : #include <slepcpep.h>
17 :
18 3 : int main(int argc,char **argv)
19 : {
20 3 : Mat A[3],P[3]; /* problem matrices and split preconditioner matrices */
21 3 : PEP pep; /* polynomial eigenproblem solver context */
22 3 : ST st;
23 3 : PetscInt N,n=10,m,Istart,Iend,II,i,j;
24 3 : PetscBool flag,terse;
25 :
26 3 : PetscFunctionBeginUser;
27 3 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
28 :
29 3 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
30 3 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
31 3 : if (!flag) m=n;
32 3 : N = n*m;
33 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nQuadratic Eigenproblem, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
34 :
35 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36 : Compute the matrices for (k^2*A_2+k*A_1+A_0)x=0, and their approximations
37 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
38 :
39 : /* A[0] is the 2-D Laplacian */
40 3 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A[0]));
41 3 : PetscCall(MatSetSizes(A[0],PETSC_DECIDE,PETSC_DECIDE,N,N));
42 3 : PetscCall(MatSetFromOptions(A[0]));
43 3 : PetscCall(MatCreate(PETSC_COMM_WORLD,&P[0]));
44 3 : PetscCall(MatSetSizes(P[0],PETSC_DECIDE,PETSC_DECIDE,N,N));
45 3 : PetscCall(MatSetFromOptions(P[0]));
46 :
47 3 : PetscCall(MatGetOwnershipRange(A[0],&Istart,&Iend));
48 435 : for (II=Istart;II<Iend;II++) {
49 432 : i = II/n; j = II-i*n;
50 432 : if (i>0) PetscCall(MatSetValue(A[0],II,II-n,-1.0,INSERT_VALUES));
51 432 : if (i<m-1) PetscCall(MatSetValue(A[0],II,II+n,-1.0,INSERT_VALUES));
52 432 : if (j>0) PetscCall(MatSetValue(A[0],II,II-1,-1.0,INSERT_VALUES));
53 432 : if (j<n-1) PetscCall(MatSetValue(A[0],II,II+1,-1.0,INSERT_VALUES));
54 432 : PetscCall(MatSetValue(A[0],II,II,4.0,INSERT_VALUES));
55 432 : if (j>0) PetscCall(MatSetValue(P[0],II,II-1,-1.0,INSERT_VALUES));
56 432 : if (j<n-1) PetscCall(MatSetValue(P[0],II,II+1,-1.0,INSERT_VALUES));
57 432 : PetscCall(MatSetValue(P[0],II,II,4.0,INSERT_VALUES));
58 : }
59 3 : PetscCall(MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY));
60 3 : PetscCall(MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY));
61 3 : PetscCall(MatAssemblyBegin(P[0],MAT_FINAL_ASSEMBLY));
62 3 : PetscCall(MatAssemblyEnd(P[0],MAT_FINAL_ASSEMBLY));
63 :
64 : /* A[1] is the 1-D Laplacian on horizontal lines */
65 3 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A[1]));
66 3 : PetscCall(MatSetSizes(A[1],PETSC_DECIDE,PETSC_DECIDE,N,N));
67 3 : PetscCall(MatSetFromOptions(A[1]));
68 3 : PetscCall(MatCreate(PETSC_COMM_WORLD,&P[1]));
69 3 : PetscCall(MatSetSizes(P[1],PETSC_DECIDE,PETSC_DECIDE,N,N));
70 3 : PetscCall(MatSetFromOptions(P[1]));
71 :
72 3 : PetscCall(MatGetOwnershipRange(A[1],&Istart,&Iend));
73 435 : for (II=Istart;II<Iend;II++) {
74 432 : i = II/n; j = II-i*n;
75 432 : if (j>0) PetscCall(MatSetValue(A[1],II,II-1,-1.0,INSERT_VALUES));
76 432 : if (j<n-1) PetscCall(MatSetValue(A[1],II,II+1,-1.0,INSERT_VALUES));
77 432 : PetscCall(MatSetValue(A[1],II,II,2.0,INSERT_VALUES));
78 432 : PetscCall(MatSetValue(P[1],II,II,2.0,INSERT_VALUES));
79 : }
80 3 : PetscCall(MatAssemblyBegin(A[1],MAT_FINAL_ASSEMBLY));
81 3 : PetscCall(MatAssemblyEnd(A[1],MAT_FINAL_ASSEMBLY));
82 3 : PetscCall(MatAssemblyBegin(P[1],MAT_FINAL_ASSEMBLY));
83 3 : PetscCall(MatAssemblyEnd(P[1],MAT_FINAL_ASSEMBLY));
84 :
85 : /* A[2] is a diagonal matrix */
86 3 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A[2]));
87 3 : PetscCall(MatSetSizes(A[2],PETSC_DECIDE,PETSC_DECIDE,N,N));
88 3 : PetscCall(MatSetFromOptions(A[2]));
89 3 : PetscCall(MatCreate(PETSC_COMM_WORLD,&P[2]));
90 3 : PetscCall(MatSetSizes(P[2],PETSC_DECIDE,PETSC_DECIDE,N,N));
91 3 : PetscCall(MatSetFromOptions(P[2]));
92 :
93 3 : PetscCall(MatGetOwnershipRange(A[2],&Istart,&Iend));
94 435 : for (II=Istart;II<Iend;II++) {
95 432 : PetscCall(MatSetValue(A[2],II,II,(PetscReal)(II+1),INSERT_VALUES));
96 432 : PetscCall(MatSetValue(P[2],II,II,(PetscReal)(II+1),INSERT_VALUES));
97 : }
98 3 : PetscCall(MatAssemblyBegin(A[2],MAT_FINAL_ASSEMBLY));
99 3 : PetscCall(MatAssemblyEnd(A[2],MAT_FINAL_ASSEMBLY));
100 3 : PetscCall(MatAssemblyBegin(P[2],MAT_FINAL_ASSEMBLY));
101 3 : PetscCall(MatAssemblyEnd(P[2],MAT_FINAL_ASSEMBLY));
102 :
103 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104 : Create the eigensolver and set various options
105 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106 :
107 3 : PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
108 3 : PetscCall(PEPSetOperators(pep,3,A));
109 3 : PetscCall(PEPSetProblemType(pep,PEP_HERMITIAN));
110 :
111 3 : PetscCall(PEPGetST(pep,&st));
112 3 : PetscCall(STSetType(st,STSINVERT));
113 3 : PetscCall(STSetSplitPreconditioner(st,3,P,SUBSET_NONZERO_PATTERN));
114 :
115 3 : PetscCall(PEPSetTarget(pep,-2.0));
116 3 : PetscCall(PEPSetWhichEigenpairs(pep,PEP_TARGET_MAGNITUDE));
117 :
118 : /*
119 : Set solver parameters at runtime
120 : */
121 3 : PetscCall(PEPSetFromOptions(pep));
122 :
123 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124 : Solve the eigensystem, display solution and clean up
125 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126 :
127 3 : PetscCall(PEPSolve(pep));
128 : /* show detailed info unless -terse option is given by user */
129 3 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
130 3 : if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
131 : else {
132 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
133 0 : PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
134 0 : PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD));
135 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
136 : }
137 3 : PetscCall(PEPDestroy(&pep));
138 3 : PetscCall(MatDestroy(&A[0]));
139 3 : PetscCall(MatDestroy(&A[1]));
140 3 : PetscCall(MatDestroy(&A[2]));
141 3 : PetscCall(MatDestroy(&P[0]));
142 3 : PetscCall(MatDestroy(&P[1]));
143 3 : PetscCall(MatDestroy(&P[2]));
144 3 : PetscCall(SlepcFinalize());
145 : return 0;
146 : }
147 :
148 : /*TEST
149 :
150 : testset:
151 : args: -pep_nev 4 -pep_ncv 28 -n 12 -terse
152 : output_file: output/ex50_1.out
153 : requires: double
154 : test:
155 : suffix: 1
156 : args: -pep_type {{toar qarnoldi}}
157 : test:
158 : suffix: 1_linear
159 : args: -pep_type linear -pep_general
160 :
161 : testset:
162 : args: -pep_all -n 12 -pep_type ciss -rg_type ellipse -rg_ellipse_center -1+1.5i -rg_ellipse_radius .3 -terse
163 : output_file: output/ex50_2.out
164 : requires: complex double
165 : timeoutfactor: 2
166 : test:
167 : suffix: 2
168 : test:
169 : suffix: 2_par
170 : nsize: 2
171 : args: -pep_ciss_partitions 2
172 :
173 : TEST*/
|
