Actual source code: test1.c
slepc-3.16.2 2022-02-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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);
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Apply the transformed operator for several ST's
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: /* shift, sigma=0.0 */
93: STSetUp(st);
94: STGetType(st,&type);
95: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
96: STApply(st,v,w);
97: VecView(w,NULL);
98: STApplyTranspose(st,v,w);
99: VecView(w,NULL);
101: /* shift, sigma=0.1 */
102: sigma = 0.1;
103: STSetShift(st,sigma);
104: STGetShift(st,&sigma);
105: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
106: STApply(st,v,w);
107: VecView(w,NULL);
109: /* sinvert, sigma=0.1 */
110: STPostSolve(st); /* undo changes if inplace */
111: STSetType(st,STSINVERT);
112: STGetKSP(st,&ksp);
113: KSPSetType(ksp,KSPGMRES);
114: KSPGetPC(ksp,&pc);
115: PCSetType(pc,PCJACOBI);
116: STGetType(st,&type);
117: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
118: STGetShift(st,&sigma);
119: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
120: STApply(st,v,w);
121: VecView(w,NULL);
123: /* sinvert, sigma=-0.5 */
124: sigma = -0.5;
125: STSetShift(st,sigma);
126: STGetShift(st,&sigma);
127: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
128: STApply(st,v,w);
129: VecView(w,NULL);
131: STDestroy(&st);
132: MatDestroy(&A);
133: MatDestroy(&S);
134: VecDestroy(&v);
135: VecDestroy(&w);
136: SlepcFinalize();
137: return ierr;
138: }
140: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y)
141: {
142: PetscErrorCode ierr;
143: Mat *A;
146: MatShellGetContext(S,&A);
147: MatMult(*A,x,y);
148: return(0);
149: }
151: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y)
152: {
153: PetscErrorCode ierr;
154: Mat *A;
157: MatShellGetContext(S,&A);
158: MatMultTranspose(*A,x,y);
159: return(0);
160: }
162: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag)
163: {
164: PetscErrorCode ierr;
165: Mat *A;
168: MatShellGetContext(S,&A);
169: MatGetDiagonal(*A,diag);
170: return(0);
171: }
173: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M)
174: {
176: Mat *A;
179: MatShellGetContext(S,&A);
180: MyShellMatCreate(A,M);
181: return(0);
182: }
184: /*TEST
186: test:
187: suffix: 1
188: args: -st_matmode {{inplace shell}}
189: output_file: output/test1_1.out
190: requires: !single
192: TEST*/