Actual source code: test1.c
slepc-3.6.1 2015-09-03
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2015, 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[] = "Test ST with shell matrices.\n\n";
24: #include <slepcst.h>
26: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag);
27: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y);
28: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y);
29: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M);
33: static PetscErrorCode MyShellMatCreate(Mat *A,Mat *M)
34: {
36: MPI_Comm comm;
37: PetscInt n;
40: MatGetSize(*A,&n,NULL);
41: PetscObjectGetComm((PetscObject)*A,&comm);
42: MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,n,n,A,M);
43: MatSetFromOptions(*M);
44: MatShellSetOperation(*M,MATOP_MULT,(void(*)())MatMult_Shell);
45: MatShellSetOperation(*M,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_Shell);
46: MatShellSetOperation(*M,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Shell);
47: MatShellSetOperation(*M,MATOP_DUPLICATE,(void(*)())MatDuplicate_Shell);
48: return(0);
49: }
53: int main(int argc,char **argv)
54: {
55: Mat A,S,mat[1];
56: ST st;
57: Vec v,w;
58: STType type;
59: KSP ksp;
60: PC pc;
61: PetscScalar sigma;
62: PetscInt n=10,i,Istart,Iend;
65: SlepcInitialize(&argc,&argv,(char*)0,help);
66: PetscOptionsGetInt(NULL,"-n",&n,NULL);
67: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian with shell matrices, n=%D\n\n",n);
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Compute the operator matrix for the 1-D Laplacian
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: MatCreate(PETSC_COMM_WORLD,&A);
74: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
75: MatSetFromOptions(A);
76: MatSetUp(A);
78: MatGetOwnershipRange(A,&Istart,&Iend);
79: for (i=Istart;i<Iend;i++) {
80: if (i>0) { MatSetValue(A,i,i-1,-1.0,INSERT_VALUES); }
81: if (i<n-1) { MatSetValue(A,i,i+1,-1.0,INSERT_VALUES); }
82: MatSetValue(A,i,i,2.0,INSERT_VALUES);
83: }
84: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
85: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
87: /* create the shell version of A */
88: MyShellMatCreate(&A,&S);
90: /* work vectors */
91: MatCreateVecs(A,&v,&w);
92: VecSet(v,1.0);
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Create the spectral transformation object
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: STCreate(PETSC_COMM_WORLD,&st);
99: mat[0] = S;
100: STSetOperators(st,1,mat);
101: STSetTransform(st,PETSC_TRUE);
102: STSetFromOptions(st);
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Apply the transformed operator for several ST's
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: /* shift, sigma=0.0 */
110: STSetUp(st);
111: STGetType(st,&type);
112: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
113: STApply(st,v,w);
114: VecView(w,NULL);
116: /* shift, sigma=0.1 */
117: sigma = 0.1;
118: STSetShift(st,sigma);
119: STGetShift(st,&sigma);
120: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
121: STApply(st,v,w);
122: VecView(w,NULL);
124: /* sinvert, sigma=0.1 */
125: STPostSolve(st); /* undo changes if inplace */
126: STSetType(st,STSINVERT);
127: STGetKSP(st,&ksp);
128: KSPSetType(ksp,KSPGMRES);
129: KSPGetPC(ksp,&pc);
130: PCSetType(pc,PCJACOBI);
131: STGetType(st,&type);
132: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
133: STGetShift(st,&sigma);
134: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
135: STApply(st,v,w);
136: VecView(w,NULL);
138: /* sinvert, sigma=-0.5 */
139: sigma = -0.5;
140: STSetShift(st,sigma);
141: STGetShift(st,&sigma);
142: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
143: STApply(st,v,w);
144: VecView(w,NULL);
146: STDestroy(&st);
147: MatDestroy(&A);
148: MatDestroy(&S);
149: VecDestroy(&v);
150: VecDestroy(&w);
151: SlepcFinalize();
152: return 0;
153: }
157: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y)
158: {
159: PetscErrorCode ierr;
160: Mat *A;
163: MatShellGetContext(S,&A);
164: MatMult(*A,x,y);
165: return(0);
166: }
170: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y)
171: {
172: PetscErrorCode ierr;
173: Mat *A;
176: MatShellGetContext(S,&A);
177: MatMultTranspose(*A,x,y);
178: return(0);
179: }
183: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag)
184: {
185: PetscErrorCode ierr;
186: Mat *A;
189: MatShellGetContext(S,&A);
190: MatGetDiagonal(*A,diag);
191: return(0);
192: }
196: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M)
197: {
199: Mat *A;
202: MatShellGetContext(S,&A);
203: MyShellMatCreate(A,M);
204: return(0);
205: }