Actual source code: test2.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2013, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Test ST with one matrix.\n\n";
24: #include <slepcst.h>
28: int main(int argc,char **argv)
29: {
30: Mat A,B,mat[1];
31: ST st;
32: Vec v,w;
33: STType type;
34: PetscScalar value[3],sigma,tau;
35: PetscInt n=10,i,Istart,Iend,col[3];
36: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
39: SlepcInitialize(&argc,&argv,(char*)0,help);
40: PetscOptionsGetInt(NULL,"-n",&n,NULL);
41: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian, n=%D\n\n",n);
43: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: Compute the operator matrix for the 1-D Laplacian
45: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
47: MatCreate(PETSC_COMM_WORLD,&A);
48: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
49: MatSetFromOptions(A);
50: MatSetUp(A);
52: MatGetOwnershipRange(A,&Istart,&Iend);
53: if (Istart==0) FirstBlock=PETSC_TRUE;
54: if (Iend==n) LastBlock=PETSC_TRUE;
55: value[0]=-1.0; value[1]=2.0; value[2]=-1.0;
56: for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
57: col[0]=i-1; col[1]=i; col[2]=i+1;
58: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
59: }
60: if (LastBlock) {
61: i=n-1; col[0]=n-2; col[1]=n-1;
62: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
63: }
64: if (FirstBlock) {
65: i=0; col[0]=0; col[1]=1; value[0]=2.0; value[1]=-1.0;
66: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
67: }
69: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
70: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
71: MatGetVecs(A,&v,&w);
72: VecSet(v,1.0);
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Create the spectral transformation object
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: STCreate(PETSC_COMM_WORLD,&st);
79: mat[0] = A;
80: STSetOperators(st,1,mat);
81: STSetFromOptions(st);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Apply the transformed operator for several ST's
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: /* shift, sigma=0.0 */
88: STSetUp(st);
89: STGetType(st,&type);
90: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
91: STApply(st,v,w);
92: VecView(w,NULL);
94: /* shift, sigma=0.1 */
95: sigma = 0.1;
96: STSetShift(st,sigma);
97: STGetShift(st,&sigma);
98: PetscPrintf(PETSC_COMM_WORLD,"With shift=%G\n",PetscRealPart(sigma));
99: STApply(st,v,w);
100: VecView(w,NULL);
102: /* sinvert, sigma=0.1 */
103: STPostSolve(st); /* undo changes if inplace */
104: STSetType(st,STSINVERT);
105: STGetType(st,&type);
106: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
107: STGetShift(st,&sigma);
108: PetscPrintf(PETSC_COMM_WORLD,"With shift=%G\n",PetscRealPart(sigma));
109: STApply(st,v,w);
110: VecView(w,NULL);
112: /* sinvert, sigma=-0.5 */
113: sigma = -0.5;
114: STSetShift(st,sigma);
115: STGetShift(st,&sigma);
116: PetscPrintf(PETSC_COMM_WORLD,"With shift=%G\n",PetscRealPart(sigma));
117: STApply(st,v,w);
118: VecView(w,NULL);
120: /* cayley, sigma=-0.5, tau=-0.5 (equal to sigma by default) */
121: STPostSolve(st); /* undo changes if inplace */
122: STSetType(st,STCAYLEY);
123: STSetUp(st);
124: STGetType(st,&type);
125: PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
126: STGetShift(st,&sigma);
127: STCayleyGetAntishift(st,&tau);
128: PetscPrintf(PETSC_COMM_WORLD,"With shift=%G, antishift=%G\n",PetscRealPart(sigma),PetscRealPart(tau));
129: STApply(st,v,w);
130: VecView(w,NULL);
132: /* cayley, sigma=1.1, tau=1.1 (still equal to sigma) */
133: sigma = 1.1;
134: STSetShift(st,sigma);
135: STGetShift(st,&sigma);
136: STCayleyGetAntishift(st,&tau);
137: PetscPrintf(PETSC_COMM_WORLD,"With shift=%G, antishift=%G\n",PetscRealPart(sigma),PetscRealPart(tau));
138: STApply(st,v,w);
139: VecView(w,NULL);
141: /* cayley, sigma=1.1, tau=-1.0 */
142: tau = -1.0;
143: STCayleySetAntishift(st,tau);
144: STGetShift(st,&sigma);
145: STCayleyGetAntishift(st,&tau);
146: PetscPrintf(PETSC_COMM_WORLD,"With shift=%G, antishift=%G\n",PetscRealPart(sigma),PetscRealPart(tau));
147: STApply(st,v,w);
148: VecView(w,NULL);
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Check inner product matrix in Cayley
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: STGetBilinearForm(st,&B);
154: MatMult(B,v,w);
155: VecView(w,NULL);
157: STDestroy(&st);
158: MatDestroy(&A);
159: MatDestroy(&B);
160: VecDestroy(&v);
161: VecDestroy(&w);
162: SlepcFinalize();
163: return 0;
164: }