Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : static char help[] = "Test ST with four matrices.\n\n";
12 :
13 : #include <slepcst.h>
14 :
15 4 : int main(int argc,char **argv)
16 : {
17 4 : Mat A,B,C,D,mat[4];
18 4 : ST st;
19 4 : KSP ksp;
20 4 : Vec v,w;
21 4 : STType type;
22 4 : PetscScalar sigma;
23 4 : PetscInt n=10,i,Istart,Iend;
24 :
25 4 : PetscFunctionBeginUser;
26 4 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
27 4 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
28 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nTest ST with four matrices, n=%" PetscInt_FMT "\n\n",n));
29 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30 : Compute the operator matrices
31 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
32 :
33 4 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
34 4 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
35 4 : PetscCall(MatSetFromOptions(A));
36 :
37 4 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
38 4 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n));
39 4 : PetscCall(MatSetFromOptions(B));
40 :
41 4 : PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
42 4 : PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
43 4 : PetscCall(MatSetFromOptions(C));
44 :
45 4 : PetscCall(MatCreate(PETSC_COMM_WORLD,&D));
46 4 : PetscCall(MatSetSizes(D,PETSC_DECIDE,PETSC_DECIDE,n,n));
47 4 : PetscCall(MatSetFromOptions(D));
48 :
49 4 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
50 44 : for (i=Istart;i<Iend;i++) {
51 40 : PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
52 40 : if (i>0) {
53 36 : PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
54 36 : PetscCall(MatSetValue(B,i,i,(PetscScalar)i,INSERT_VALUES));
55 4 : } else PetscCall(MatSetValue(B,i,i,-1.0,INSERT_VALUES));
56 40 : if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
57 40 : PetscCall(MatSetValue(C,i,n-i-1,1.0,INSERT_VALUES));
58 40 : PetscCall(MatSetValue(D,i,i,i*.1,INSERT_VALUES));
59 40 : if (i==0) PetscCall(MatSetValue(D,0,n-1,1.0,INSERT_VALUES));
60 40 : if (i==n-1) PetscCall(MatSetValue(D,n-1,0,1.0,INSERT_VALUES));
61 : }
62 :
63 4 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
64 4 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
65 4 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
66 4 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
67 4 : PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
68 4 : PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
69 4 : PetscCall(MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY));
70 4 : PetscCall(MatAssemblyEnd(D,MAT_FINAL_ASSEMBLY));
71 4 : PetscCall(MatCreateVecs(A,&v,&w));
72 4 : PetscCall(VecSet(v,1.0));
73 :
74 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75 : Create the spectral transformation object
76 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77 :
78 4 : PetscCall(STCreate(PETSC_COMM_WORLD,&st));
79 4 : mat[0] = A;
80 4 : mat[1] = B;
81 4 : mat[2] = C;
82 4 : mat[3] = D;
83 4 : PetscCall(STSetMatrices(st,4,mat));
84 4 : PetscCall(STGetKSP(st,&ksp));
85 4 : PetscCall(KSPSetTolerances(ksp,100*PETSC_MACHINE_EPSILON,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
86 4 : PetscCall(STSetFromOptions(st));
87 :
88 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89 : Apply the transformed operator for several ST's
90 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91 :
92 : /* shift, sigma=0.0 */
93 4 : PetscCall(STSetUp(st));
94 4 : PetscCall(STGetType(st,&type));
95 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type));
96 20 : for (i=0;i<4;i++) {
97 16 : PetscCall(STMatMult(st,i,v,w));
98 16 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"k= %" PetscInt_FMT "\n",i));
99 16 : PetscCall(VecView(w,NULL));
100 : }
101 4 : PetscCall(STMatSolve(st,v,w));
102 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"solve\n"));
103 4 : PetscCall(VecView(w,NULL));
104 :
105 : /* shift, sigma=0.1 */
106 4 : sigma = 0.1;
107 4 : PetscCall(STSetShift(st,sigma));
108 4 : PetscCall(STGetShift(st,&sigma));
109 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
110 20 : for (i=0;i<4;i++) {
111 16 : PetscCall(STMatMult(st,i,v,w));
112 16 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"k= %" PetscInt_FMT "\n",i));
113 16 : PetscCall(VecView(w,NULL));
114 : }
115 4 : PetscCall(STMatSolve(st,v,w));
116 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"solve\n"));
117 4 : PetscCall(VecView(w,NULL));
118 :
119 : /* sinvert, sigma=0.1 */
120 4 : PetscCall(STPostSolve(st));
121 4 : PetscCall(STSetType(st,STSINVERT));
122 4 : PetscCall(STGetType(st,&type));
123 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type));
124 4 : PetscCall(STGetShift(st,&sigma));
125 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
126 20 : for (i=0;i<4;i++) {
127 16 : PetscCall(STMatMult(st,i,v,w));
128 16 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"k= %" PetscInt_FMT "\n",i));
129 16 : PetscCall(VecView(w,NULL));
130 : }
131 4 : PetscCall(STMatSolve(st,v,w));
132 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"solve\n"));
133 4 : PetscCall(VecView(w,NULL));
134 :
135 : /* sinvert, sigma=-0.5 */
136 4 : sigma = -0.5;
137 4 : PetscCall(STSetShift(st,sigma));
138 4 : PetscCall(STGetShift(st,&sigma));
139 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
140 20 : for (i=0;i<4;i++) {
141 16 : PetscCall(STMatMult(st,i,v,w));
142 16 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"k= %" PetscInt_FMT "\n",i));
143 16 : PetscCall(VecView(w,NULL));
144 : }
145 4 : PetscCall(STMatSolve(st,v,w));
146 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"solve\n"));
147 4 : PetscCall(VecView(w,NULL));
148 :
149 4 : PetscCall(STDestroy(&st));
150 4 : PetscCall(MatDestroy(&A));
151 4 : PetscCall(MatDestroy(&B));
152 4 : PetscCall(MatDestroy(&C));
153 4 : PetscCall(MatDestroy(&D));
154 4 : PetscCall(VecDestroy(&v));
155 4 : PetscCall(VecDestroy(&w));
156 4 : PetscCall(SlepcFinalize());
157 : return 0;
158 : }
159 :
160 : /*TEST
161 :
162 : test:
163 : suffix: 1
164 : args: -st_transform -st_matmode {{copy shell}}
165 : requires: !single
166 :
167 : test:
168 : suffix: 2
169 : args: -st_matmode {{copy shell}}
170 :
171 : TEST*/
|
