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[] = "A quadratic eigenproblem defined using shell matrices.\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = number of grid subdivisions in x and y dimensions.\n\n";
14 :
15 : #include <slepcpep.h>
16 :
17 : /*
18 : User-defined routines
19 : */
20 : PetscErrorCode MatMult_Laplacian2D(Mat A,Vec x,Vec y);
21 : PetscErrorCode MatGetDiagonal_Laplacian2D(Mat A,Vec diag);
22 : PetscErrorCode MatMult_Zero(Mat A,Vec x,Vec y);
23 : PetscErrorCode MatGetDiagonal_Zero(Mat A,Vec diag);
24 : PetscErrorCode MatMult_Identity(Mat A,Vec x,Vec y);
25 : PetscErrorCode MatGetDiagonal_Identity(Mat A,Vec diag);
26 :
27 3 : int main(int argc,char **argv)
28 : {
29 3 : Mat M,C,K,A[3]; /* problem matrices */
30 3 : PEP pep; /* polynomial eigenproblem solver context */
31 3 : PEPType type;
32 3 : PetscInt N,n=10,nev;
33 3 : PetscMPIInt size;
34 3 : PetscBool terse;
35 3 : ST st;
36 :
37 3 : PetscFunctionBeginUser;
38 3 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
39 3 : PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
40 3 : PetscCheck(size==1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only");
41 :
42 3 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
43 3 : N = n*n;
44 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nQuadratic Eigenproblem with shell matrices, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,n));
45 :
46 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47 : Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
48 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49 :
50 : /* K is the 2-D Laplacian */
51 3 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,&n,&K));
52 3 : PetscCall(MatShellSetOperation(K,MATOP_MULT,(void(*)(void))MatMult_Laplacian2D));
53 3 : PetscCall(MatShellSetOperation(K,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_Laplacian2D));
54 3 : PetscCall(MatShellSetOperation(K,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Laplacian2D));
55 :
56 : /* C is the zero matrix */
57 3 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,NULL,&C));
58 3 : PetscCall(MatShellSetOperation(C,MATOP_MULT,(void(*)(void))MatMult_Zero));
59 3 : PetscCall(MatShellSetOperation(C,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_Zero));
60 3 : PetscCall(MatShellSetOperation(C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Zero));
61 :
62 : /* M is the identity matrix */
63 3 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,NULL,&M));
64 3 : PetscCall(MatShellSetOperation(M,MATOP_MULT,(void(*)(void))MatMult_Identity));
65 3 : PetscCall(MatShellSetOperation(M,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_Identity));
66 3 : PetscCall(MatShellSetOperation(M,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Identity));
67 :
68 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69 : Create the eigensolver and set various options
70 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71 :
72 : /*
73 : Create eigensolver context
74 : */
75 3 : PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
76 :
77 : /*
78 : Set matrices and problem type
79 : */
80 3 : A[0] = K; A[1] = C; A[2] = M;
81 3 : PetscCall(PEPSetOperators(pep,3,A));
82 3 : PetscCall(PEPGetST(pep,&st));
83 3 : PetscCall(STSetMatMode(st,ST_MATMODE_SHELL));
84 :
85 : /*
86 : Set solver parameters at runtime
87 : */
88 3 : PetscCall(PEPSetFromOptions(pep));
89 :
90 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91 : Solve the eigensystem
92 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93 :
94 3 : PetscCall(PEPSolve(pep));
95 :
96 : /*
97 : Optional: Get some information from the solver and display it
98 : */
99 3 : PetscCall(PEPGetType(pep,&type));
100 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
101 3 : PetscCall(PEPGetDimensions(pep,&nev,NULL,NULL));
102 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
103 :
104 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105 : Display solution and clean up
106 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107 :
108 : /* show detailed info unless -terse option is given by user */
109 3 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
110 3 : if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_RELATIVE,NULL));
111 : else {
112 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
113 0 : PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
114 0 : PetscCall(PEPErrorView(pep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
115 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
116 : }
117 3 : PetscCall(PEPDestroy(&pep));
118 3 : PetscCall(MatDestroy(&M));
119 3 : PetscCall(MatDestroy(&C));
120 3 : PetscCall(MatDestroy(&K));
121 3 : PetscCall(SlepcFinalize());
122 : return 0;
123 : }
124 :
125 : /*
126 : Compute the matrix vector multiplication y<---T*x where T is a nx by nx
127 : tridiagonal matrix with DD on the diagonal, DL on the subdiagonal, and
128 : DU on the superdiagonal.
129 : */
130 4940 : static void tv(int nx,const PetscScalar *x,PetscScalar *y)
131 : {
132 4940 : PetscScalar dd,dl,du;
133 4940 : int j;
134 :
135 4940 : dd = 4.0;
136 4940 : dl = -1.0;
137 4940 : du = -1.0;
138 :
139 4940 : y[0] = dd*x[0] + du*x[1];
140 44460 : for (j=1;j<nx-1;j++)
141 39520 : y[j] = dl*x[j-1] + dd*x[j] + du*x[j+1];
142 4940 : y[nx-1] = dl*x[nx-2] + dd*x[nx-1];
143 4940 : }
144 :
145 : /*
146 : Matrix-vector product subroutine for the 2D Laplacian.
147 :
148 : The matrix used is the 2 dimensional discrete Laplacian on unit square with
149 : zero Dirichlet boundary condition.
150 :
151 : Computes y <-- A*x, where A is the block tridiagonal matrix
152 :
153 : | T -I |
154 : |-I T -I |
155 : A = | -I T |
156 : | ... -I|
157 : | -I T|
158 :
159 : The subroutine TV is called to compute y<--T*x.
160 : */
161 494 : PetscErrorCode MatMult_Laplacian2D(Mat A,Vec x,Vec y)
162 : {
163 494 : void *ctx;
164 494 : int nx,lo,i,j;
165 494 : const PetscScalar *px;
166 494 : PetscScalar *py;
167 :
168 494 : PetscFunctionBeginUser;
169 494 : PetscCall(MatShellGetContext(A,&ctx));
170 494 : nx = *(int*)ctx;
171 494 : PetscCall(VecGetArrayRead(x,&px));
172 494 : PetscCall(VecGetArray(y,&py));
173 :
174 494 : tv(nx,&px[0],&py[0]);
175 5928 : for (i=0;i<nx;i++) py[i] -= px[nx+i];
176 :
177 4446 : for (j=2;j<nx;j++) {
178 3952 : lo = (j-1)*nx;
179 3952 : tv(nx,&px[lo],&py[lo]);
180 47424 : for (i=0;i<nx;i++) py[lo+i] -= px[lo-nx+i] + px[lo+nx+i];
181 : }
182 :
183 494 : lo = (nx-1)*nx;
184 494 : tv(nx,&px[lo],&py[lo]);
185 5928 : for (i=0;i<nx;i++) py[lo+i] -= px[lo-nx+i];
186 :
187 494 : PetscCall(VecRestoreArrayRead(x,&px));
188 494 : PetscCall(VecRestoreArray(y,&py));
189 494 : PetscFunctionReturn(PETSC_SUCCESS);
190 : }
191 :
192 2 : PetscErrorCode MatGetDiagonal_Laplacian2D(Mat A,Vec diag)
193 : {
194 2 : PetscFunctionBeginUser;
195 2 : PetscCall(VecSet(diag,4.0));
196 2 : PetscFunctionReturn(PETSC_SUCCESS);
197 : }
198 :
199 : /*
200 : Matrix-vector product subroutine for the Null matrix.
201 : */
202 494 : PetscErrorCode MatMult_Zero(Mat A,Vec x,Vec y)
203 : {
204 494 : PetscFunctionBeginUser;
205 494 : PetscCall(VecSet(y,0.0));
206 494 : PetscFunctionReturn(PETSC_SUCCESS);
207 : }
208 :
209 2 : PetscErrorCode MatGetDiagonal_Zero(Mat A,Vec diag)
210 : {
211 2 : PetscFunctionBeginUser;
212 2 : PetscCall(VecSet(diag,0.0));
213 2 : PetscFunctionReturn(PETSC_SUCCESS);
214 : }
215 :
216 : /*
217 : Matrix-vector product subroutine for the Identity matrix.
218 : */
219 301 : PetscErrorCode MatMult_Identity(Mat A,Vec x,Vec y)
220 : {
221 301 : PetscFunctionBeginUser;
222 301 : PetscCall(VecCopy(x,y));
223 301 : PetscFunctionReturn(PETSC_SUCCESS);
224 : }
225 :
226 3 : PetscErrorCode MatGetDiagonal_Identity(Mat A,Vec diag)
227 : {
228 3 : PetscFunctionBeginUser;
229 3 : PetscCall(VecSet(diag,1.0));
230 3 : PetscFunctionReturn(PETSC_SUCCESS);
231 : }
232 :
233 : /*TEST
234 :
235 : test:
236 : suffix: 1
237 : args: -pep_type {{toar qarnoldi linear}} -pep_nev 4 -terse
238 : filter: grep -v Solution | sed -e "s/2.7996[1-8]i/2.79964i/g" | sed -e "s/2.7570[5-9]i/2.75708i/g" | sed -e "s/0.00000-2.79964i, 0.00000+2.79964i/0.00000+2.79964i, 0.00000-2.79964i/" | sed -e "s/0.00000-2.75708i, 0.00000+2.75708i/0.00000+2.75708i, 0.00000-2.75708i/"
239 :
240 : TEST*/
|
