Actual source code: test3.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 two matrices.\n\n";
13: #include <slepcst.h>
15: int main(int argc,char **argv)
16: {
17: Mat A,B,M,mat[2];
18: ST st;
19: Vec v,w;
20: STType type;
21: PetscScalar sigma,tau;
22: PetscInt n=10,i,Istart,Iend;
25: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
26: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
27: 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: MatGetOwnershipRange(A,&Istart,&Iend);
44: for (i=Istart;i<Iend;i++) {
45: MatSetValue(A,i,i,2.0,INSERT_VALUES);
46: if (i>0) {
47: MatSetValue(A,i,i-1,-1.0,INSERT_VALUES);
48: MatSetValue(B,i,i,(PetscScalar)i,INSERT_VALUES);
49: } else {
50: MatSetValue(B,i,i,-1.0,INSERT_VALUES);
51: }
52: if (i<n-1) { MatSetValue(A,i,i+1,-1.0,INSERT_VALUES); }
53: }
54: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
55: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
56: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
57: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
58: MatCreateVecs(A,&v,&w);
59: VecSet(v,1.0);
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Create the spectral transformation object
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: STCreate(PETSC_COMM_WORLD,&st);
66: mat[0] = A;
67: mat[1] = B;
68: STSetMatrices(st,2,mat);
69: STSetTransform(st,PETSC_TRUE);
70: STSetFromOptions(st);
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Apply the transformed operator for several ST's
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: /* shift, sigma=0.0 */
77: STSetUp(st);
78: STGetType(st,&type);
79: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
80: STApply(st,v,w);
81: VecView(w,NULL);
83: /* shift, sigma=0.1 */
84: sigma = 0.1;
85: STSetShift(st,sigma);
86: STGetShift(st,&sigma);
87: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
88: STApply(st,v,w);
89: VecView(w,NULL);
91: /* sinvert, sigma=0.1 */
92: STPostSolve(st); /* undo changes if inplace */
93: STSetType(st,STSINVERT);
94: STGetType(st,&type);
95: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
96: STGetShift(st,&sigma);
97: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
98: STApply(st,v,w);
99: VecView(w,NULL);
101: /* sinvert, sigma=-0.5 */
102: sigma = -0.5;
103: STSetShift(st,sigma);
104: STGetShift(st,&sigma);
105: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
106: STApply(st,v,w);
107: VecView(w,NULL);
109: /* cayley, sigma=-0.5, tau=-0.5 (equal to sigma by default) */
110: STPostSolve(st); /* undo changes if inplace */
111: STSetType(st,STCAYLEY);
112: STSetUp(st);
113: STGetType(st,&type);
114: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
115: STGetShift(st,&sigma);
116: STCayleyGetAntishift(st,&tau);
117: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g, antishift=%g\n",(double)PetscRealPart(sigma),(double)PetscRealPart(tau));
118: STApply(st,v,w);
119: VecView(w,NULL);
121: /* cayley, sigma=1.1, tau=1.1 (still equal to sigma) */
122: sigma = 1.1;
123: STSetShift(st,sigma);
124: STGetShift(st,&sigma);
125: STCayleyGetAntishift(st,&tau);
126: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g, antishift=%g\n",(double)PetscRealPart(sigma),(double)PetscRealPart(tau));
127: STApply(st,v,w);
128: VecView(w,NULL);
130: /* cayley, sigma=1.1, tau=-1.0 */
131: tau = -1.0;
132: STCayleySetAntishift(st,tau);
133: STGetShift(st,&sigma);
134: STCayleyGetAntishift(st,&tau);
135: PetscPrintf(PETSC_COMM_WORLD,"With shift=%g, antishift=%g\n",(double)PetscRealPart(sigma),(double)PetscRealPart(tau));
136: STApply(st,v,w);
137: VecView(w,NULL);
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Check inner product matrix in Cayley
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142: STGetBilinearForm(st,&M);
143: MatMult(M,v,w);
144: VecView(w,NULL);
146: STDestroy(&st);
147: MatDestroy(&A);
148: MatDestroy(&B);
149: MatDestroy(&M);
150: VecDestroy(&v);
151: VecDestroy(&w);
152: SlepcFinalize();
153: return ierr;
154: }
156: /*TEST
158: test:
159: suffix: 1
160: args: -st_matmode {{copy inplace shell}}
161: output_file: output/test3_1.out
162: requires: !single
164: TEST*/