Actual source code: test4.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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;
26: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
27: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
28: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian plus diagonal, n=%D\n\n",n);
29: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: Compute the operator matrix for the 1-D Laplacian
31: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
33: MatCreate(PETSC_COMM_WORLD,&A);
34: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
35: MatSetFromOptions(A);
36: MatSetUp(A);
38: MatCreate(PETSC_COMM_WORLD,&B);
39: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);
40: MatSetFromOptions(B);
41: MatSetUp(B);
43: MatCreate(PETSC_COMM_WORLD,&C);
44: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
45: MatSetFromOptions(C);
46: MatSetUp(C);
48: MatCreate(PETSC_COMM_WORLD,&D);
49: MatSetSizes(D,PETSC_DECIDE,PETSC_DECIDE,n,n);
50: MatSetFromOptions(D);
51: MatSetUp(D);
53: MatGetOwnershipRange(A,&Istart,&Iend);
54: for (i=Istart;i<Iend;i++) {
55: MatSetValue(A,i,i,2.0,INSERT_VALUES);
56: if (i>0) {
57: MatSetValue(A,i,i-1,-1.0,INSERT_VALUES);
58: MatSetValue(B,i,i,(PetscScalar)i,INSERT_VALUES);
59: } else {
60: MatSetValue(B,i,i,-1.0,INSERT_VALUES);
61: }
62: if (i<n-1) {
63: MatSetValue(A,i,i+1,-1.0,INSERT_VALUES);
64: }
65: MatSetValue(C,i,n-i-1,1.0,INSERT_VALUES);
66: MatSetValue(D,i,i,i*.1,INSERT_VALUES);
67: if (i==0) {
68: MatSetValue(D,0,n-1,1.0,INSERT_VALUES);
69: }
70: if (i==n-1) {
71: MatSetValue(D,n-1,0,1.0,INSERT_VALUES);
72: }
73: }
75: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
77: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
79: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
81: MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY);
82: MatAssemblyEnd(D,MAT_FINAL_ASSEMBLY);
83: MatCreateVecs(A,&v,&w);
84: VecSet(v,1.0);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Create the spectral transformation object
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: STCreate(PETSC_COMM_WORLD,&st);
90: mat[0] = A;
91: mat[1] = B;
92: mat[2] = C;
93: mat[3] = D;
94: STSetMatrices(st,4,mat);
95: STGetKSP(st,&ksp);
96: KSPSetTolerances(ksp,100*PETSC_MACHINE_EPSILON,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
97: STSetFromOptions(st);
98: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: Apply the transformed operator for several ST's
100: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: /* shift, sigma=0.0 */
102: STSetUp(st);
103: STGetType(st,&type);
104: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
105: for (i=0;i<4;i++) {
106: STMatMult(st,i,v,w);
107: PetscPrintf(PETSC_COMM_WORLD,"k= %D\n",i);
108: VecView(w,NULL);
109: }
110: STMatSolve(st,v,w);
111: PetscPrintf(PETSC_COMM_WORLD,"solve\n");
112: VecView(w,NULL);
114: /* shift, sigma=0.1 */
115: sigma = 0.1;
116: STSetShift(st,sigma);
117: STGetShift(st,&sigma);
118: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
119: for (i=0;i<4;i++) {
120: STMatMult(st,i,v,w);
121: PetscPrintf(PETSC_COMM_WORLD,"k= %D\n",i);
122: VecView(w,NULL);
123: }
124: STMatSolve(st,v,w);
125: PetscPrintf(PETSC_COMM_WORLD,"solve\n");
126: VecView(w,NULL);
128: /* sinvert, sigma=0.1 */
129: STPostSolve(st);
130: STSetType(st,STSINVERT);
131: STGetType(st,&type);
132: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
133: STGetShift(st,&sigma);
134: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
135: for (i=0;i<4;i++) {
136: STMatMult(st,i,v,w);
137: PetscPrintf(PETSC_COMM_WORLD,"k= %D\n",i);
138: VecView(w,NULL);
139: }
140: STMatSolve(st,v,w);
141: PetscPrintf(PETSC_COMM_WORLD,"solve\n");
142: VecView(w,NULL);
144: /* sinvert, sigma=-0.5 */
145: sigma = -0.5;
146: STSetShift(st,sigma);
147: STGetShift(st,&sigma);
148: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
149: for (i=0;i<4;i++) {
150: STMatMult(st,i,v,w);
151: PetscPrintf(PETSC_COMM_WORLD,"k= %D\n",i);
152: VecView(w,NULL);
153: }
154: STMatSolve(st,v,w);
155: PetscPrintf(PETSC_COMM_WORLD,"solve\n");
156: VecView(w,NULL);
157: STDestroy(&st);
158: MatDestroy(&A);
159: MatDestroy(&B);
160: MatDestroy(&C);
161: MatDestroy(&D);
162: VecDestroy(&v);
163: VecDestroy(&w);
164: SlepcFinalize();
165: return ierr;
166: }
168: /*TEST
170: test:
171: suffix: 1
172: args: -st_transform -st_matmode {{copy shell}}
173: output_file: output/test4_1.out
174: requires: !single
176: test:
177: suffix: 2
178: args: -st_matmode {{copy shell}}
179: output_file: output/test4_2.out
181: TEST*/