Actual source code: test6.c
slepc-3.20.2 2024-03-15
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,(char*)0,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));
37: PetscCall(MatSetUp(A));
39: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
40: for (i=Istart;i<Iend;i++) {
41: if (i>0) PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
42: if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
43: PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
44: }
45: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
46: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
47: PetscCall(MatCreateVecs(A,&v,&w));
48: PetscCall(VecSet(v,1.0));
50: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: Create the spectral transformation object
52: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: PetscCall(STCreate(PETSC_COMM_WORLD,&st));
55: mat[0] = A;
56: PetscCall(STSetMatrices(st,1,mat));
57: PetscCall(STSetType(st,STPRECOND));
58: PetscCall(STGetKSP(st,&ksp));
59: PetscCall(KSPSetType(ksp,KSPPREONLY));
60: PetscCall(STSetFromOptions(st));
62: /* set up */
63: /* - the transform flag is necessary so that A-sigma*I is built explicitly */
64: PetscCall(STSetTransform(st,PETSC_TRUE));
65: /* - the ksphasmat flag is necessary when using STApply(), otherwise can only use PCApply() */
66: PetscCall(STPrecondSetKSPHasMat(st,PETSC_TRUE));
67: /* no need to call STSetUp() explicitly */
68: PetscCall(STSetUp(st));
70: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: Apply the preconditioner
72: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: /* default shift */
75: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With default shift\n"));
76: PetscCall(STApply(st,v,w));
77: PetscCall(VecView(w,NULL));
79: /* change shift */
80: sigma = 0.1;
81: PetscCall(STSetShift(st,sigma));
82: PetscCall(STGetShift(st,&sigma));
83: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
84: PetscCall(STApply(st,v,w));
85: PetscCall(VecView(w,NULL));
86: PetscCall(STPostSolve(st)); /* undo changes if inplace */
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Test a user-provided preconditioner matrix
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: PetscCall(MatCreate(PETSC_COMM_WORLD,&P));
93: PetscCall(MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,n,n));
94: PetscCall(MatSetFromOptions(P));
95: PetscCall(MatSetUp(P));
97: PetscCall(MatGetOwnershipRange(P,&Istart,&Iend));
98: for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(P,i,i,2.0,INSERT_VALUES));
99: if (Istart==0) {
100: PetscCall(MatSetValue(P,1,0,-1.0,INSERT_VALUES));
101: PetscCall(MatSetValue(P,0,1,-1.0,INSERT_VALUES));
102: }
103: PetscCall(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY));
104: PetscCall(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY));
106: /* apply new preconditioner */
107: PetscCall(STSetPreconditionerMat(st,P));
108: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With user-provided matrix\n"));
109: PetscCall(STApply(st,v,w));
110: PetscCall(VecView(w,NULL));
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Test a user-provided preconditioner in split form
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: PetscCall(STGetMatMode(st,&matmode));
117: if (matmode==ST_MATMODE_COPY) {
118: PetscCall(STSetPreconditionerMat(st,NULL));
119: mat[0] = P;
120: PetscCall(STSetSplitPreconditioner(st,1,mat,SAME_NONZERO_PATTERN));
122: /* apply new preconditioner */
123: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With split preconditioner\n"));
124: PetscCall(STApply(st,v,w));
125: PetscCall(VecView(w,NULL));
126: }
128: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: Clean up
130: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: PetscCall(STDestroy(&st));
132: PetscCall(MatDestroy(&A));
133: PetscCall(MatDestroy(&P));
134: PetscCall(VecDestroy(&v));
135: PetscCall(VecDestroy(&w));
136: PetscCall(SlepcFinalize());
137: return 0;
138: }
140: /*TEST
142: test:
143: suffix: 1
144: args: -st_matmode {{copy inplace shell}separate output}
145: requires: !single
147: TEST*/