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 ST with shell matrices.\n\n";
12 :
13 : #include <slepcst.h>
14 :
15 : static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y);
16 : static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y);
17 : #if defined(PETSC_USE_COMPLEX)
18 : static PetscErrorCode MatMultHermitianTranspose_Shell(Mat S,Vec x,Vec y);
19 : #endif
20 : static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag);
21 : static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M);
22 :
23 2 : static PetscErrorCode MyShellMatCreate(Mat *A,Mat *M)
24 : {
25 2 : MPI_Comm comm;
26 2 : PetscInt n;
27 :
28 2 : PetscFunctionBeginUser;
29 2 : PetscCall(MatGetSize(*A,&n,NULL));
30 2 : PetscCall(PetscObjectGetComm((PetscObject)*A,&comm));
31 2 : PetscCall(MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,n,n,A,M));
32 2 : PetscCall(MatShellSetOperation(*M,MATOP_MULT,(void(*)(void))MatMult_Shell));
33 2 : PetscCall(MatShellSetOperation(*M,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell));
34 : #if defined(PETSC_USE_COMPLEX)
35 : PetscCall(MatShellSetOperation(*M,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_Shell));
36 : #endif
37 2 : PetscCall(MatShellSetOperation(*M,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell));
38 2 : PetscCall(MatShellSetOperation(*M,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Shell));
39 2 : PetscFunctionReturn(PETSC_SUCCESS);
40 : }
41 :
42 2 : int main(int argc,char **argv)
43 : {
44 2 : Mat A,S,mat[1];
45 2 : ST st;
46 2 : Vec v,w;
47 2 : STType type;
48 2 : KSP ksp;
49 2 : PC pc;
50 2 : PetscScalar sigma;
51 2 : PetscInt n=10,i,Istart,Iend;
52 :
53 2 : PetscFunctionBeginUser;
54 2 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
55 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
56 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian with shell matrices, n=%" PetscInt_FMT "\n\n",n));
57 :
58 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59 : Compute the operator matrix for the 1-D Laplacian
60 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61 :
62 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
63 2 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
64 2 : PetscCall(MatSetFromOptions(A));
65 :
66 2 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
67 22 : for (i=Istart;i<Iend;i++) {
68 20 : if (i>0) PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
69 20 : if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
70 20 : PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
71 : }
72 2 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
73 2 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
74 :
75 : /* create the shell version of A */
76 2 : PetscCall(MyShellMatCreate(&A,&S));
77 :
78 : /* work vectors */
79 2 : PetscCall(MatCreateVecs(A,&v,&w));
80 2 : PetscCall(VecSet(v,1.0));
81 :
82 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83 : Create the spectral transformation object
84 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85 :
86 2 : PetscCall(STCreate(PETSC_COMM_WORLD,&st));
87 2 : mat[0] = S;
88 2 : PetscCall(STSetMatrices(st,1,mat));
89 2 : PetscCall(STSetTransform(st,PETSC_TRUE));
90 2 : PetscCall(STSetFromOptions(st));
91 :
92 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93 : Apply the transformed operator for several ST's
94 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95 :
96 : /* shift, sigma=0.0 */
97 2 : PetscCall(STSetUp(st));
98 2 : PetscCall(STGetType(st,&type));
99 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type));
100 2 : PetscCall(STApply(st,v,w));
101 2 : PetscCall(VecView(w,NULL));
102 2 : PetscCall(STApplyHermitianTranspose(st,v,w));
103 2 : PetscCall(VecView(w,NULL));
104 :
105 : /* shift, sigma=0.1 */
106 2 : sigma = 0.1;
107 2 : PetscCall(STSetShift(st,sigma));
108 2 : PetscCall(STGetShift(st,&sigma));
109 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
110 2 : PetscCall(STApply(st,v,w));
111 2 : PetscCall(VecView(w,NULL));
112 :
113 : /* sinvert, sigma=0.1 */
114 2 : PetscCall(STPostSolve(st)); /* undo changes if inplace */
115 2 : PetscCall(STSetType(st,STSINVERT));
116 2 : PetscCall(STGetKSP(st,&ksp));
117 2 : PetscCall(KSPSetType(ksp,KSPGMRES));
118 2 : PetscCall(KSPGetPC(ksp,&pc));
119 2 : PetscCall(PCSetType(pc,PCJACOBI));
120 2 : PetscCall(STGetType(st,&type));
121 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type));
122 2 : PetscCall(STGetShift(st,&sigma));
123 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
124 2 : PetscCall(STApply(st,v,w));
125 2 : PetscCall(VecView(w,NULL));
126 :
127 : /* sinvert, sigma=-0.5 */
128 2 : sigma = -0.5;
129 2 : PetscCall(STSetShift(st,sigma));
130 2 : PetscCall(STGetShift(st,&sigma));
131 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
132 2 : PetscCall(STApply(st,v,w));
133 2 : PetscCall(VecView(w,NULL));
134 :
135 2 : PetscCall(STDestroy(&st));
136 2 : PetscCall(MatDestroy(&A));
137 2 : PetscCall(MatDestroy(&S));
138 2 : PetscCall(VecDestroy(&v));
139 2 : PetscCall(VecDestroy(&w));
140 2 : PetscCall(SlepcFinalize());
141 : return 0;
142 : }
143 :
144 24 : static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y)
145 : {
146 24 : Mat *A;
147 :
148 24 : PetscFunctionBeginUser;
149 24 : PetscCall(MatShellGetContext(S,&A));
150 24 : PetscCall(MatMult(*A,x,y));
151 24 : PetscFunctionReturn(PETSC_SUCCESS);
152 : }
153 :
154 2 : static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y)
155 : {
156 2 : Mat *A;
157 :
158 2 : PetscFunctionBeginUser;
159 2 : PetscCall(MatShellGetContext(S,&A));
160 2 : PetscCall(MatMultTranspose(*A,x,y));
161 2 : PetscFunctionReturn(PETSC_SUCCESS);
162 : }
163 :
164 : #if defined(PETSC_USE_COMPLEX)
165 : static PetscErrorCode MatMultHermitianTranspose_Shell(Mat S,Vec x,Vec y)
166 : {
167 : Mat *A;
168 :
169 : PetscFunctionBeginUser;
170 : PetscCall(MatShellGetContext(S,&A));
171 : PetscCall(MatMultHermitianTranspose(*A,x,y));
172 : PetscFunctionReturn(PETSC_SUCCESS);
173 : }
174 : #endif
175 :
176 4 : static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag)
177 : {
178 4 : Mat *A;
179 :
180 4 : PetscFunctionBeginUser;
181 4 : PetscCall(MatShellGetContext(S,&A));
182 4 : PetscCall(MatGetDiagonal(*A,diag));
183 4 : PetscFunctionReturn(PETSC_SUCCESS);
184 : }
185 :
186 0 : static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M)
187 : {
188 0 : Mat *A;
189 :
190 0 : PetscFunctionBeginUser;
191 0 : PetscCall(MatShellGetContext(S,&A));
192 0 : PetscCall(MyShellMatCreate(A,M));
193 0 : PetscFunctionReturn(PETSC_SUCCESS);
194 : }
195 :
196 : /*TEST
197 :
198 : test:
199 : suffix: 1
200 : args: -st_matmode {{inplace shell}}
201 : requires: !single
202 :
203 : TEST*/
|
