Actual source code: ex13.c
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";
5: #include <petscksp.h>
7: static PetscErrorCode Get2DStencil(PetscInt i, PetscInt j, PetscInt n, PetscInt idxs[])
8: {
9: PetscInt k = 0;
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: }
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;
33: PetscFunctionBeginUser;
34: PetscCall(PetscInitialize(&argc, &args, (char *)0, 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: }
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));
80: PetscCall(MatCreateVecs(A, &x, &b));
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));
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: }
170: /* Create linear solver context */
171: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
172: PetscCall(KSPSetOperators(ksp, A, P));
173: PetscCall(KSPSetTolerances(ksp, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
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));
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: }
193: /* solve */
194: PetscCall(KSPSolve(ksp, b, x));
196: /* test reuse */
197: PetscCall(KSPSolve(ksp, b, x));
199: /* Free work space. */
200: PetscCall(KSPDestroy(&ksp));
201: PetscCall(VecDestroy(&x));
202: PetscCall(VecDestroy(&b));
203: PetscCall(MatDestroy(&A));
204: PetscCall(MatDestroy(&P));
206: PetscCall(PetscFinalize());
207: return 0;
208: }
210: /*TEST
212: test:
213: suffix: qlp_sisc
214: args: -ksp_converged_reason -ksp_minres_monitor -ksp_view_solution
216: test:
217: suffix: qlp_sisc_none
218: args: -ksp_converged_reason -ksp_minres_monitor -ksp_view_solution -pc_type none
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"
225: test:
226: suffix: qlp_2
227: args: -ksp_converged_reason -testcase 2 -ksp_minres_monitor
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"
234: TEST*/