1 /* These tests taken from https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip
2 Examples in CSP11/Algorithms/MINRESQLP/minresQLP.m comments */
3 static char help[] = "Tests MINRES-QLP.\n\n";
4
5 #include <petscksp.h>
6
Get2DStencil(PetscInt i,PetscInt j,PetscInt n,PetscInt idxs[])7 static PetscErrorCode Get2DStencil(PetscInt i, PetscInt j, PetscInt n, PetscInt idxs[])
8 {
9 PetscInt k = 0;
10
11 PetscFunctionBeginUser;
12 for (k = 0; k < 9; k++) idxs[k] = -1;
13 k = 0;
14 for (PetscInt k1 = -1; k1 <= 1; k1++)
15 for (PetscInt k2 = -1; k2 <= 1; k2++)
16 if (i + k1 >= 0 && i + k1 < n && j + k2 >= 0 && j + k2 < n) idxs[k++] = n * (i + k1) + (j + k2);
17 PetscFunctionReturn(PETSC_SUCCESS);
18 }
19
main(int argc,char ** args)20 int main(int argc, char **args)
21 {
22 Vec x, b;
23 Mat A, P;
24 KSP ksp;
25 PC pc;
26 PetscInt testcase = 0, m, nnz, pnnz;
27 PetscMPIInt rank;
28 PCType pctype;
29 PetscBool flg;
30 PetscReal radius = 0.0;
31 PetscReal rtol = PETSC_DEFAULT;
32
33 PetscFunctionBeginUser;
34 PetscCall(PetscInitialize(&argc, &args, NULL, help));
35 PetscCall(PetscOptionsGetInt(NULL, NULL, "-testcase", &testcase, NULL));
36 switch (testcase) {
37 case 1:
38 m = 100;
39 nnz = 3;
40 pnnz = 1;
41 pctype = PCMAT;
42 rtol = 1.e-10;
43 break;
44 case 2:
45 m = 2500;
46 nnz = 9;
47 pnnz = 0;
48 pctype = PCNONE;
49 rtol = 1.e-5;
50 radius = 1.e2;
51 break;
52 case 3:
53 m = 21;
54 nnz = 1;
55 pnnz = 0;
56 pctype = PCNONE;
57 radius = 1.e2;
58 break;
59 default: /* Example 7.1 in https://stanford.edu/group/SOL/software/minresqlp/MINRESQLP-SISC-2011.pdf */
60 m = 4;
61 nnz = 3;
62 pnnz = 1;
63 pctype = PCMAT;
64 break;
65 }
66
67 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
68 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m));
69 PetscCall(MatSetFromOptions(A));
70 PetscCall(MatMPIAIJSetPreallocation(A, nnz, NULL, nnz, NULL));
71 PetscCall(MatSeqAIJSetPreallocation(A, nnz, NULL));
72 PetscCall(MatSetUp(A));
73 PetscCall(MatCreate(PETSC_COMM_WORLD, &P));
74 PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, m, m));
75 PetscCall(MatSetFromOptions(P));
76 PetscCall(MatMPIAIJSetPreallocation(P, pnnz, NULL, pnnz, NULL));
77 PetscCall(MatSeqAIJSetPreallocation(P, pnnz, NULL));
78 PetscCall(MatSetUp(P));
79
80 PetscCall(MatCreateVecs(A, &x, &b));
81
82 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
83 /* dummy assemble */
84 if (rank == 0) {
85 PetscScalar *vals;
86 PetscInt *cols;
87 PetscInt row;
88 PetscCall(PetscMalloc2(nnz, &cols, nnz, &vals));
89 switch (testcase) {
90 case 1:
91 vals[0] = -2.0;
92 vals[1] = 4.0;
93 vals[2] = -2.0;
94 for (row = 0; row < m; row++) {
95 cols[0] = row - 1;
96 cols[1] = row;
97 cols[2] = row == m - 1 ? -1 : row + 1;
98 PetscCall(MatSetValues(A, 1, &row, 3, cols, vals, INSERT_VALUES));
99 PetscCall(MatSetValue(P, row, row, 1.0 / 4.0, INSERT_VALUES));
100 }
101 break;
102 case 2:
103 for (PetscInt i = 0; i < 9; i++) vals[i] = 1.0;
104 for (PetscInt i = 0; i < 50; i++) {
105 for (PetscInt j = 0; j < 50; j++) {
106 PetscInt row = i * 50 + j;
107 PetscCall(Get2DStencil(i, j, 50, cols));
108 PetscCall(MatSetValues(A, 1, &row, 9, cols, vals, INSERT_VALUES));
109 }
110 }
111 break;
112 case 3:
113 for (row = 0; row < m; row++) {
114 if (row == 10) continue;
115 vals[0] = row - 10.0;
116 PetscCall(MatSetValue(A, row, row, vals[0], INSERT_VALUES));
117 }
118 break;
119 default:
120 vals[0] = vals[1] = vals[2] = 1.0;
121 row = 0;
122 cols[0] = 0;
123 cols[1] = 1;
124 PetscCall(MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES));
125 PetscCall(VecSetValue(b, row, 6.0, INSERT_VALUES));
126 PetscCall(MatSetValue(P, row, row, PetscSqr(0.84201), INSERT_VALUES));
127 row = 1;
128 cols[0] = 0;
129 cols[1] = 1;
130 cols[2] = 2;
131 PetscCall(MatSetValues(A, 1, &row, 3, cols, vals, INSERT_VALUES));
132 PetscCall(VecSetValue(b, row, 9.0, INSERT_VALUES));
133 PetscCall(MatSetValue(P, row, row, PetscSqr(0.81228), INSERT_VALUES));
134 row = 2;
135 cols[0] = 1;
136 cols[1] = 3;
137 PetscCall(MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES));
138 PetscCall(VecSetValue(b, row, 6.0, INSERT_VALUES));
139 PetscCall(MatSetValue(P, row, row, PetscSqr(0.30957), INSERT_VALUES));
140 row = 3;
141 cols[0] = 2;
142 PetscCall(MatSetValues(A, 1, &row, 1, cols, vals, INSERT_VALUES));
143 PetscCall(VecSetValue(b, row, 3.0, INSERT_VALUES));
144 PetscCall(MatSetValue(P, row, row, PetscSqr(3.2303), INSERT_VALUES));
145 break;
146 }
147 PetscCall(PetscFree2(cols, vals));
148 }
149 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
150 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
151 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
152 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
153
154 /* rhs */
155 switch (testcase) {
156 case 1:
157 case 2:
158 PetscCall(VecSet(x, 1.0));
159 PetscCall(MatMult(A, x, b));
160 break;
161 case 3:
162 PetscCall(VecSet(b, 1.0));
163 break;
164 default:
165 PetscCall(VecAssemblyBegin(b));
166 PetscCall(VecAssemblyEnd(b));
167 break;
168 }
169
170 /* Create linear solver context */
171 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
172 PetscCall(KSPSetOperators(ksp, A, P));
173 PetscCall(KSPSetTolerances(ksp, rtol, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
174 PetscCall(KSPSetType(ksp, KSPMINRES));
175 PetscCall(KSPMINRESSetUseQLP(ksp, PETSC_TRUE));
176 if (radius > 0.0) PetscCall(KSPMINRESSetRadius(ksp, radius));
177 PetscCall(KSPGetPC(ksp, &pc));
178 PetscCall(PCSetType(pc, pctype));
179 PetscCall(KSPSetFromOptions(ksp));
180 PetscCall(KSPSetUp(ksp));
181
182 /* print info */
183 PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPMINRES, &flg));
184 if (flg) {
185 PetscCall(KSPMINRESGetUseQLP(ksp, &flg));
186 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Using KSPMINRES%s\n", flg ? "-QLP" : ""));
187 } else {
188 KSPType ksptype;
189 PetscCall(KSPGetType(ksp, &ksptype));
190 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Using %s\n", ksptype));
191 }
192
193 /* solve */
194 PetscCall(KSPSolve(ksp, b, x));
195
196 /* test reuse */
197 PetscCall(KSPSolve(ksp, b, x));
198
199 /* Free work space. */
200 PetscCall(KSPDestroy(&ksp));
201 PetscCall(VecDestroy(&x));
202 PetscCall(VecDestroy(&b));
203 PetscCall(MatDestroy(&A));
204 PetscCall(MatDestroy(&P));
205
206 PetscCall(PetscFinalize());
207 return 0;
208 }
209
210 /*TEST
211
212 test:
213 suffix: qlp_sisc
214 args: -ksp_converged_reason -ksp_minres_monitor -ksp_view_solution
215
216 test:
217 suffix: qlp_sisc_none
218 args: -ksp_converged_reason -ksp_minres_monitor -ksp_view_solution -pc_type none
219
220 test:
221 suffix: qlp_1
222 args: -ksp_converged_reason -testcase 1 -ksp_minres_monitor
223 filter: sed -e "s/49 3/49 1/g" -e "s/50 3/50 1/g" -e "s/CONVERGED_HAPPY_BREAKDOWN/CONVERGED_RTOL/g"
224
225 test:
226 suffix: qlp_2
227 args: -ksp_converged_reason -testcase 2 -ksp_minres_monitor
228
229 test:
230 suffix: qlp_3
231 args: -ksp_converged_reason -testcase 3 -ksp_minres_monitor
232 filter: sed -e "s/24 2/24 6/g" -e "s/50 3/50 1/g" -e "s/CONVERGED_RTOL/CONVERGED_STEP_LENGTH/g"
233
234 TEST*/
235