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*/