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