Actual source code: test6.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 STPRECOND operations.\n\n";
13: #include <slepcst.h>
15: int main(int argc,char **argv)
16: {
17: Mat A,P,mat[1];
18: ST st;
19: KSP ksp;
20: Vec v,w;
21: PetscScalar sigma;
22: PetscInt n=10,i,Istart,Iend;
23: STMatMode matmode;
25: PetscFunctionBeginUser;
26: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
27: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
28: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nPreconditioner for 1-D Laplacian, n=%" PetscInt_FMT "\n\n",n));
30: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31: Compute the operator matrix for the 1-D Laplacian
32: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
34: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
35: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
36: PetscCall(MatSetFromOptions(A));
38: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
39: for (i=Istart;i<Iend;i++) {
40: if (i>0) PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
41: if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
42: PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
43: }
44: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
45: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
46: PetscCall(MatCreateVecs(A,&v,&w));
47: PetscCall(VecSet(v,1.0));
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Create the spectral transformation object
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53: PetscCall(STCreate(PETSC_COMM_WORLD,&st));
54: mat[0] = A;
55: PetscCall(STSetMatrices(st,1,mat));
56: PetscCall(STSetType(st,STPRECOND));
57: PetscCall(STGetKSP(st,&ksp));
58: PetscCall(KSPSetType(ksp,KSPPREONLY));
59: PetscCall(STSetFromOptions(st));
61: /* set up */
62: /* - the transform flag is necessary so that A-sigma*I is built explicitly */
63: PetscCall(STSetTransform(st,PETSC_TRUE));
64: /* - the ksphasmat flag is necessary when using STApply(), otherwise can only use PCApply() */
65: PetscCall(STPrecondSetKSPHasMat(st,PETSC_TRUE));
66: /* no need to call STSetUp() explicitly */
67: PetscCall(STSetUp(st));
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Apply the preconditioner
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: /* default shift */
74: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With default shift\n"));
75: PetscCall(STApply(st,v,w));
76: PetscCall(VecView(w,NULL));
78: /* change shift */
79: sigma = 0.1;
80: PetscCall(STSetShift(st,sigma));
81: PetscCall(STGetShift(st,&sigma));
82: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
83: PetscCall(STApply(st,v,w));
84: PetscCall(VecView(w,NULL));
85: PetscCall(STPostSolve(st)); /* undo changes if inplace */
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Test a user-provided preconditioner matrix
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: PetscCall(MatCreate(PETSC_COMM_WORLD,&P));
92: PetscCall(MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,n,n));
93: PetscCall(MatSetFromOptions(P));
95: PetscCall(MatGetOwnershipRange(P,&Istart,&Iend));
96: for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(P,i,i,2.0,INSERT_VALUES));
97: if (Istart==0) {
98: PetscCall(MatSetValue(P,1,0,-1.0,INSERT_VALUES));
99: PetscCall(MatSetValue(P,0,1,-1.0,INSERT_VALUES));
100: }
101: PetscCall(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY));
102: PetscCall(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY));
104: /* apply new preconditioner */
105: PetscCall(STSetPreconditionerMat(st,P));
106: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With user-provided matrix\n"));
107: PetscCall(STApply(st,v,w));
108: PetscCall(VecView(w,NULL));
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Test a user-provided preconditioner in split form
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: PetscCall(STGetMatMode(st,&matmode));
115: if (matmode==ST_MATMODE_COPY) {
116: PetscCall(STSetPreconditionerMat(st,NULL));
117: mat[0] = P;
118: PetscCall(STSetSplitPreconditioner(st,1,mat,SAME_NONZERO_PATTERN));
120: /* apply new preconditioner */
121: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With split preconditioner\n"));
122: PetscCall(STApply(st,v,w));
123: PetscCall(VecView(w,NULL));
124: }
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Clean up
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: PetscCall(STDestroy(&st));
130: PetscCall(MatDestroy(&A));
131: PetscCall(MatDestroy(&P));
132: PetscCall(VecDestroy(&v));
133: PetscCall(VecDestroy(&w));
134: PetscCall(SlepcFinalize());
135: return 0;
136: }
138: /*TEST
140: test:
141: suffix: 1
142: args: -st_matmode {{copy inplace shell}separate output}
143: requires: !single
145: TEST*/