Actual source code: test1.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: static char help[] = "Test ST with shell matrices.\n\n";
13: #include <slepcst.h>
15: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag);
16: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y);
17: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y);
18: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M);
20: static PetscErrorCode MyShellMatCreate(Mat *A,Mat *M)
21: {
23: MPI_Comm comm;
24: PetscInt n;
27: MatGetSize(*A,&n,NULL);
28: PetscObjectGetComm((PetscObject)*A,&comm);
29: MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,n,n,A,M);
30: MatShellSetOperation(*M,MATOP_MULT,(void(*)(void))MatMult_Shell);
31: MatShellSetOperation(*M,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell);
32: MatShellSetOperation(*M,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell);
33: MatShellSetOperation(*M,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Shell);
34: return(0);
35: }
37: int main(int argc,char **argv)
38: {
39: Mat A,S,mat[1];
40: ST st;
41: Vec v,w;
42: STType type;
43: KSP ksp;
44: PC pc;
45: PetscScalar sigma;
46: PetscInt n=10,i,Istart,Iend;
49: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
50: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
51: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian with shell matrices, n=%D\n\n",n);
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Compute the operator matrix for the 1-D Laplacian
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: MatCreate(PETSC_COMM_WORLD,&A);
58: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
59: MatSetFromOptions(A);
60: MatSetUp(A);
62: MatGetOwnershipRange(A,&Istart,&Iend);
63: for (i=Istart;i<Iend;i++) {
64: if (i>0) { MatSetValue(A,i,i-1,-1.0,INSERT_VALUES); }
65: if (i<n-1) { MatSetValue(A,i,i+1,-1.0,INSERT_VALUES); }
66: MatSetValue(A,i,i,2.0,INSERT_VALUES);
67: }
68: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
71: /* create the shell version of A */
72: MyShellMatCreate(&A,&S);
74: /* work vectors */
75: MatCreateVecs(A,&v,&w);
76: VecSet(v,1.0);
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Create the spectral transformation object
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82: STCreate(PETSC_COMM_WORLD,&st);
83: mat[0] = S;
84: STSetMatrices(st,1,mat);
85: STSetTransform(st,PETSC_TRUE);
86: STSetFromOptions(st);
89: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: Apply the transformed operator for several ST's
91: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: /* shift, sigma=0.0 */
94: STSetUp(st);
95: STGetType(st,&type);
96: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
97: STApply(st,v,w);
98: VecView(w,NULL);
99: STApplyTranspose(st,v,w);
100: VecView(w,NULL);
102: /* shift, sigma=0.1 */
103: sigma = 0.1;
104: STSetShift(st,sigma);
105: STGetShift(st,&sigma);
106: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
107: STApply(st,v,w);
108: VecView(w,NULL);
110: /* sinvert, sigma=0.1 */
111: STPostSolve(st); /* undo changes if inplace */
112: STSetType(st,STSINVERT);
113: STGetKSP(st,&ksp);
114: KSPSetType(ksp,KSPGMRES);
115: KSPGetPC(ksp,&pc);
116: PCSetType(pc,PCJACOBI);
117: STGetType(st,&type);
118: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
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.5 */
125: sigma = -0.5;
126: STSetShift(st,sigma);
127: STGetShift(st,&sigma);
128: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
129: STApply(st,v,w);
130: VecView(w,NULL);
132: STDestroy(&st);
133: MatDestroy(&A);
134: MatDestroy(&S);
135: VecDestroy(&v);
136: VecDestroy(&w);
137: SlepcFinalize();
138: return ierr;
139: }
141: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y)
142: {
143: PetscErrorCode ierr;
144: Mat *A;
147: MatShellGetContext(S,&A);
148: MatMult(*A,x,y);
149: return(0);
150: }
152: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y)
153: {
154: PetscErrorCode ierr;
155: Mat *A;
158: MatShellGetContext(S,&A);
159: MatMultTranspose(*A,x,y);
160: return(0);
161: }
163: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag)
164: {
165: PetscErrorCode ierr;
166: Mat *A;
169: MatShellGetContext(S,&A);
170: MatGetDiagonal(*A,diag);
171: return(0);
172: }
174: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M)
175: {
177: Mat *A;
180: MatShellGetContext(S,&A);
181: MyShellMatCreate(A,M);
182: return(0);
183: }
185: /*TEST
187: test:
188: suffix: 1
189: args: -st_matmode {{inplace shell}}
190: output_file: output/test1_1.out
191: requires: !single
193: TEST*/