Actual source code: test4.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 ST with four matrices.\n\n";
13: #include <slepcst.h>
15: int main(int argc,char **argv)
16: {
17: Mat A,B,C,D,mat[4];
18: ST st;
19: KSP ksp;
20: Vec v,w;
21: STType type;
22: PetscScalar sigma;
23: PetscInt n=10,i,Istart,Iend;
25: PetscFunctionBeginUser;
26: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
27: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
28: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nTest ST with four matrices, n=%" PetscInt_FMT "\n\n",n));
29: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: Compute the operator matrices
31: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
33: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
34: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
35: PetscCall(MatSetFromOptions(A));
36: PetscCall(MatSetUp(A));
38: PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
39: PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n));
40: PetscCall(MatSetFromOptions(B));
41: PetscCall(MatSetUp(B));
43: PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
44: PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
45: PetscCall(MatSetFromOptions(C));
46: PetscCall(MatSetUp(C));
48: PetscCall(MatCreate(PETSC_COMM_WORLD,&D));
49: PetscCall(MatSetSizes(D,PETSC_DECIDE,PETSC_DECIDE,n,n));
50: PetscCall(MatSetFromOptions(D));
51: PetscCall(MatSetUp(D));
53: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
54: for (i=Istart;i<Iend;i++) {
55: PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
56: if (i>0) {
57: PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
58: PetscCall(MatSetValue(B,i,i,(PetscScalar)i,INSERT_VALUES));
59: } else PetscCall(MatSetValue(B,i,i,-1.0,INSERT_VALUES));
60: if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
61: PetscCall(MatSetValue(C,i,n-i-1,1.0,INSERT_VALUES));
62: PetscCall(MatSetValue(D,i,i,i*.1,INSERT_VALUES));
63: if (i==0) PetscCall(MatSetValue(D,0,n-1,1.0,INSERT_VALUES));
64: if (i==n-1) PetscCall(MatSetValue(D,n-1,0,1.0,INSERT_VALUES));
65: }
67: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
68: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
69: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
70: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
71: PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
72: PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
73: PetscCall(MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY));
74: PetscCall(MatAssemblyEnd(D,MAT_FINAL_ASSEMBLY));
75: PetscCall(MatCreateVecs(A,&v,&w));
76: PetscCall(VecSet(v,1.0));
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Create the spectral transformation object
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82: PetscCall(STCreate(PETSC_COMM_WORLD,&st));
83: mat[0] = A;
84: mat[1] = B;
85: mat[2] = C;
86: mat[3] = D;
87: PetscCall(STSetMatrices(st,4,mat));
88: PetscCall(STGetKSP(st,&ksp));
89: PetscCall(KSPSetTolerances(ksp,100*PETSC_MACHINE_EPSILON,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
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: for (i=0;i<4;i++) {
101: PetscCall(STMatMult(st,i,v,w));
102: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"k= %" PetscInt_FMT "\n",i));
103: PetscCall(VecView(w,NULL));
104: }
105: PetscCall(STMatSolve(st,v,w));
106: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"solve\n"));
107: PetscCall(VecView(w,NULL));
109: /* shift, sigma=0.1 */
110: sigma = 0.1;
111: PetscCall(STSetShift(st,sigma));
112: PetscCall(STGetShift(st,&sigma));
113: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
114: for (i=0;i<4;i++) {
115: PetscCall(STMatMult(st,i,v,w));
116: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"k= %" PetscInt_FMT "\n",i));
117: PetscCall(VecView(w,NULL));
118: }
119: PetscCall(STMatSolve(st,v,w));
120: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"solve\n"));
121: PetscCall(VecView(w,NULL));
123: /* sinvert, sigma=0.1 */
124: PetscCall(STPostSolve(st));
125: PetscCall(STSetType(st,STSINVERT));
126: PetscCall(STGetType(st,&type));
127: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type));
128: PetscCall(STGetShift(st,&sigma));
129: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
130: for (i=0;i<4;i++) {
131: PetscCall(STMatMult(st,i,v,w));
132: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"k= %" PetscInt_FMT "\n",i));
133: PetscCall(VecView(w,NULL));
134: }
135: PetscCall(STMatSolve(st,v,w));
136: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"solve\n"));
137: PetscCall(VecView(w,NULL));
139: /* sinvert, sigma=-0.5 */
140: sigma = -0.5;
141: PetscCall(STSetShift(st,sigma));
142: PetscCall(STGetShift(st,&sigma));
143: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
144: for (i=0;i<4;i++) {
145: PetscCall(STMatMult(st,i,v,w));
146: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"k= %" PetscInt_FMT "\n",i));
147: PetscCall(VecView(w,NULL));
148: }
149: PetscCall(STMatSolve(st,v,w));
150: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"solve\n"));
151: PetscCall(VecView(w,NULL));
153: PetscCall(STDestroy(&st));
154: PetscCall(MatDestroy(&A));
155: PetscCall(MatDestroy(&B));
156: PetscCall(MatDestroy(&C));
157: PetscCall(MatDestroy(&D));
158: PetscCall(VecDestroy(&v));
159: PetscCall(VecDestroy(&w));
160: PetscCall(SlepcFinalize());
161: return 0;
162: }
164: /*TEST
166: test:
167: suffix: 1
168: args: -st_transform -st_matmode {{copy shell}}
169: output_file: output/test4_1.out
170: requires: !single
172: test:
173: suffix: 2
174: args: -st_matmode {{copy shell}}
175: output_file: output/test4_2.out
177: TEST*/