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