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[] = "Solves the same eigenproblem as in example ex2, but using a shell matrix. "
12 : "The problem is a standard symmetric eigenproblem corresponding to the 2-D Laplacian operator.\n\n"
13 : "The command line options are:\n"
14 : " -n <n>, where <n> = number of grid subdivisions in both x and y dimensions.\n\n";
15 :
16 : #include <slepceps.h>
17 :
18 : /*
19 : User-defined routines
20 : */
21 : PetscErrorCode MatMult_Laplacian2D(Mat A,Vec x,Vec y);
22 : PetscErrorCode MatGetDiagonal_Laplacian2D(Mat A,Vec diag);
23 :
24 1 : int main(int argc,char **argv)
25 : {
26 1 : Mat A; /* operator matrix */
27 1 : EPS eps; /* eigenproblem solver context */
28 1 : EPSType type;
29 1 : PetscMPIInt size;
30 1 : PetscInt N,n=10,nev;
31 1 : PetscBool terse;
32 :
33 1 : PetscFunctionBeginUser;
34 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
35 1 : PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
36 1 : PetscCheck(size==1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only");
37 :
38 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
39 1 : N = n*n;
40 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n2-D Laplacian Eigenproblem (matrix-free version), N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,n));
41 :
42 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43 : Create the operator matrix that defines the eigensystem, Ax=kx
44 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
45 :
46 1 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,&n,&A));
47 1 : PetscCall(MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_Laplacian2D));
48 1 : PetscCall(MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_Laplacian2D));
49 1 : PetscCall(MatShellSetOperation(A,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Laplacian2D));
50 :
51 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52 : Create the eigensolver and set various options
53 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54 :
55 : /*
56 : Create eigensolver context
57 : */
58 1 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
59 :
60 : /*
61 : Set operators. In this case, it is a standard eigenvalue problem
62 : */
63 1 : PetscCall(EPSSetOperators(eps,A,NULL));
64 1 : PetscCall(EPSSetProblemType(eps,EPS_HEP));
65 :
66 : /*
67 : Set solver parameters at runtime
68 : */
69 1 : PetscCall(EPSSetFromOptions(eps));
70 :
71 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72 : Solve the eigensystem
73 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74 :
75 1 : PetscCall(EPSSolve(eps));
76 :
77 : /*
78 : Optional: Get some information from the solver and display it
79 : */
80 1 : PetscCall(EPSGetType(eps,&type));
81 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
82 1 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
83 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
84 :
85 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86 : Display solution and clean up
87 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88 :
89 : /* show detailed info unless -terse option is given by user */
90 1 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
91 1 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
92 : else {
93 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
94 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
95 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
96 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
97 : }
98 1 : PetscCall(EPSDestroy(&eps));
99 1 : PetscCall(MatDestroy(&A));
100 1 : PetscCall(SlepcFinalize());
101 : return 0;
102 : }
103 :
104 : /*
105 : Compute the matrix vector multiplication y<---T*x where T is a nx by nx
106 : tridiagonal matrix with DD on the diagonal, DL on the subdiagonal, and
107 : DU on the superdiagonal.
108 : */
109 25128 : static void tv(int nx,const PetscScalar *x,PetscScalar *y)
110 : {
111 25128 : PetscScalar dd,dl,du;
112 25128 : int j;
113 :
114 25128 : dd = 4.0;
115 25128 : dl = -1.0;
116 25128 : du = -1.0;
117 :
118 25128 : y[0] = dd*x[0] + du*x[1];
119 1784088 : for (j=1;j<nx-1;j++)
120 1758960 : y[j] = dl*x[j-1] + dd*x[j] + du*x[j+1];
121 25128 : y[nx-1] = dl*x[nx-2] + dd*x[nx-1];
122 25128 : }
123 :
124 : /*
125 : Matrix-vector product subroutine for the 2D Laplacian.
126 :
127 : The matrix used is the 2 dimensional discrete Laplacian on unit square with
128 : zero Dirichlet boundary condition.
129 :
130 : Computes y <-- A*x, where A is the block tridiagonal matrix
131 :
132 : | T -I |
133 : |-I T -I |
134 : A = | -I T |
135 : | ... -I|
136 : | -I T|
137 :
138 : The subroutine TV is called to compute y<--T*x.
139 : */
140 349 : PetscErrorCode MatMult_Laplacian2D(Mat A,Vec x,Vec y)
141 : {
142 349 : void *ctx;
143 349 : int nx,lo,i,j;
144 349 : const PetscScalar *px;
145 349 : PetscScalar *py;
146 :
147 349 : PetscFunctionBeginUser;
148 349 : PetscCall(MatShellGetContext(A,&ctx));
149 349 : nx = *(int*)ctx;
150 349 : PetscCall(VecGetArrayRead(x,&px));
151 349 : PetscCall(VecGetArray(y,&py));
152 :
153 349 : tv(nx,&px[0],&py[0]);
154 25826 : for (i=0;i<nx;i++) py[i] -= px[nx+i];
155 :
156 24779 : for (j=2;j<nx;j++) {
157 24430 : lo = (j-1)*nx;
158 24430 : tv(nx,&px[lo],&py[lo]);
159 1807820 : for (i=0;i<nx;i++) py[lo+i] -= px[lo-nx+i] + px[lo+nx+i];
160 : }
161 :
162 349 : lo = (nx-1)*nx;
163 349 : tv(nx,&px[lo],&py[lo]);
164 25826 : for (i=0;i<nx;i++) py[lo+i] -= px[lo-nx+i];
165 :
166 349 : PetscCall(VecRestoreArrayRead(x,&px));
167 349 : PetscCall(VecRestoreArray(y,&py));
168 349 : PetscFunctionReturn(PETSC_SUCCESS);
169 : }
170 :
171 0 : PetscErrorCode MatGetDiagonal_Laplacian2D(Mat A,Vec diag)
172 : {
173 0 : PetscFunctionBeginUser;
174 0 : PetscCall(VecSet(diag,4.0));
175 0 : PetscFunctionReturn(PETSC_SUCCESS);
176 : }
177 :
178 : /*TEST
179 :
180 : test:
181 : suffix: 1
182 : args: -n 72 -eps_nev 4 -eps_ncv 20 -terse
183 : requires: !single
184 :
185 : TEST*/
|
