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 and split preconditioner.\n\n";
12 :
13 : #include <slepcst.h>
14 :
15 2 : int main(int argc,char **argv)
16 : {
17 2 : Mat A,B,C,D,Pa,Pb,Pc,Pd,Pmat,mat[4];
18 2 : ST st;
19 2 : KSP ksp;
20 2 : PC pc;
21 2 : Vec v,w;
22 2 : STType type;
23 2 : PetscScalar sigma;
24 2 : PetscInt n=10,i,Istart,Iend;
25 :
26 2 : PetscFunctionBeginUser;
27 2 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
28 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
29 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nTest ST with four matrices, n=%" PetscInt_FMT "\n\n",n));
30 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31 : Compute the operator matrices
32 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
33 :
34 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
35 2 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
36 2 : PetscCall(MatSetFromOptions(A));
37 :
38 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
39 2 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n));
40 2 : PetscCall(MatSetFromOptions(B));
41 :
42 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
43 2 : PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
44 2 : PetscCall(MatSetFromOptions(C));
45 :
46 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&D));
47 2 : PetscCall(MatSetSizes(D,PETSC_DECIDE,PETSC_DECIDE,n,n));
48 2 : PetscCall(MatSetFromOptions(D));
49 :
50 2 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
51 22 : for (i=Istart;i<Iend;i++) {
52 20 : PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
53 20 : if (i>0) {
54 18 : PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
55 18 : PetscCall(MatSetValue(B,i,i,(PetscScalar)i,INSERT_VALUES));
56 2 : } else PetscCall(MatSetValue(B,i,i,-1.0,INSERT_VALUES));
57 20 : if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
58 20 : PetscCall(MatSetValue(C,i,n-i-1,1.0,INSERT_VALUES));
59 20 : PetscCall(MatSetValue(D,i,i,i*.1,INSERT_VALUES));
60 20 : if (i==0) PetscCall(MatSetValue(D,0,n-1,1.0,INSERT_VALUES));
61 20 : if (i==n-1) PetscCall(MatSetValue(D,n-1,0,1.0,INSERT_VALUES));
62 : }
63 :
64 2 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
65 2 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
66 2 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
67 2 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
68 2 : PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
69 2 : PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
70 2 : PetscCall(MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY));
71 2 : PetscCall(MatAssemblyEnd(D,MAT_FINAL_ASSEMBLY));
72 2 : PetscCall(MatCreateVecs(A,&v,&w));
73 2 : PetscCall(VecSet(v,1.0));
74 :
75 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76 : Compute the split preconditioner matrices (four diagonals)
77 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78 :
79 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&Pa));
80 2 : PetscCall(MatSetSizes(Pa,PETSC_DECIDE,PETSC_DECIDE,n,n));
81 2 : PetscCall(MatSetFromOptions(Pa));
82 :
83 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&Pb));
84 2 : PetscCall(MatSetSizes(Pb,PETSC_DECIDE,PETSC_DECIDE,n,n));
85 2 : PetscCall(MatSetFromOptions(Pb));
86 :
87 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&Pc));
88 2 : PetscCall(MatSetSizes(Pc,PETSC_DECIDE,PETSC_DECIDE,n,n));
89 2 : PetscCall(MatSetFromOptions(Pc));
90 :
91 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&Pd));
92 2 : PetscCall(MatSetSizes(Pd,PETSC_DECIDE,PETSC_DECIDE,n,n));
93 2 : PetscCall(MatSetFromOptions(Pd));
94 :
95 2 : PetscCall(MatGetOwnershipRange(Pa,&Istart,&Iend));
96 22 : for (i=Istart;i<Iend;i++) {
97 20 : PetscCall(MatSetValue(Pa,i,i,2.0,INSERT_VALUES));
98 20 : if (i>0) PetscCall(MatSetValue(Pb,i,i,(PetscScalar)i,INSERT_VALUES));
99 2 : else PetscCall(MatSetValue(Pb,i,i,-1.0,INSERT_VALUES));
100 20 : PetscCall(MatSetValue(Pd,i,i,i*.1,INSERT_VALUES));
101 : }
102 :
103 2 : PetscCall(MatAssemblyBegin(Pa,MAT_FINAL_ASSEMBLY));
104 2 : PetscCall(MatAssemblyEnd(Pa,MAT_FINAL_ASSEMBLY));
105 2 : PetscCall(MatAssemblyBegin(Pb,MAT_FINAL_ASSEMBLY));
106 2 : PetscCall(MatAssemblyEnd(Pb,MAT_FINAL_ASSEMBLY));
107 2 : PetscCall(MatAssemblyBegin(Pc,MAT_FINAL_ASSEMBLY));
108 2 : PetscCall(MatAssemblyEnd(Pc,MAT_FINAL_ASSEMBLY));
109 2 : PetscCall(MatAssemblyBegin(Pd,MAT_FINAL_ASSEMBLY));
110 2 : PetscCall(MatAssemblyEnd(Pd,MAT_FINAL_ASSEMBLY));
111 :
112 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113 : Create the spectral transformation object
114 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115 :
116 2 : PetscCall(STCreate(PETSC_COMM_WORLD,&st));
117 2 : mat[0] = A;
118 2 : mat[1] = B;
119 2 : mat[2] = C;
120 2 : mat[3] = D;
121 2 : PetscCall(STSetMatrices(st,4,mat));
122 2 : mat[0] = Pa;
123 2 : mat[1] = Pb;
124 2 : mat[2] = Pc;
125 2 : mat[3] = Pd;
126 2 : PetscCall(STSetSplitPreconditioner(st,4,mat,SUBSET_NONZERO_PATTERN));
127 2 : PetscCall(STGetKSP(st,&ksp));
128 2 : PetscCall(KSPSetTolerances(ksp,100*PETSC_MACHINE_EPSILON,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
129 2 : PetscCall(STSetTransform(st,PETSC_TRUE));
130 2 : PetscCall(STSetFromOptions(st));
131 2 : PetscCall(STGetKSP(st,&ksp));
132 2 : PetscCall(KSPGetPC(ksp,&pc));
133 :
134 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135 : Apply the operator
136 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137 :
138 : /* sigma=0.0 */
139 2 : PetscCall(STSetUp(st));
140 2 : PetscCall(STGetType(st,&type));
141 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type));
142 2 : PetscCall(PCGetOperators(pc,NULL,&Pmat));
143 2 : PetscCall(MatView(Pmat,NULL));
144 2 : PetscCall(STMatSolve(st,v,w));
145 2 : PetscCall(VecView(w,NULL));
146 :
147 : /* sigma=0.1 */
148 2 : sigma = 0.1;
149 2 : PetscCall(STSetShift(st,sigma));
150 2 : PetscCall(STGetShift(st,&sigma));
151 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
152 2 : PetscCall(PCGetOperators(pc,NULL,&Pmat));
153 2 : PetscCall(MatView(Pmat,NULL));
154 2 : PetscCall(STMatSolve(st,v,w));
155 2 : PetscCall(VecView(w,NULL));
156 :
157 2 : PetscCall(STDestroy(&st));
158 2 : PetscCall(MatDestroy(&A));
159 2 : PetscCall(MatDestroy(&B));
160 2 : PetscCall(MatDestroy(&C));
161 2 : PetscCall(MatDestroy(&D));
162 2 : PetscCall(MatDestroy(&Pa));
163 2 : PetscCall(MatDestroy(&Pb));
164 2 : PetscCall(MatDestroy(&Pc));
165 2 : PetscCall(MatDestroy(&Pd));
166 2 : PetscCall(VecDestroy(&v));
167 2 : PetscCall(VecDestroy(&w));
168 2 : PetscCall(SlepcFinalize());
169 : return 0;
170 : }
171 :
172 : /*TEST
173 :
174 : test:
175 : suffix: 1
176 : args: -st_type {{shift sinvert}separate output} -st_pc_type jacobi
177 : requires: !single
178 :
179 : TEST*/
|
