Actual source code: ex2.c

  1: static char help[] = "Tests repeated solving linear system on 2 by 2 matrix provided by MUMPS developer, Dec 17, 2012.\n\n";
  2: /*
  3: We have investigated the problem further, and we have
  4: been able to reproduce it and obtain an erroneous
  5: solution with an even smaller, 2x2, matrix:
  6:     [1 2]
  7:     [2 3]
  8: and a right-hand side vector with all ones (1,1)
  9: The correct solution is the vector (-1,1), in both solves.

 11: mpiexec -n 2 ./ex2 -ksp_type preonly -pc_type lu  -pc_factor_mat_solver_type mumps  -mat_mumps_icntl_7 6 -mat_mumps_cntl_1 0.99

 13: With this combination of options, I get off-diagonal pivots during the
 14: factorization, which is the cause of the problem (different isol_loc
 15: returned in the second solve, whereas, as I understand it, Petsc expects
 16: isol_loc not to change between successive solves).
 17: */

 19: #include <petscksp.h>

 21: int main(int argc, char **args)
 22: {
 23:   Mat         C;
 24:   PetscInt    N = 2, rowidx, colidx;
 25:   Vec         u, b, r;
 26:   KSP         ksp;
 27:   PetscReal   norm;
 28:   PetscMPIInt rank, size;
 29:   PetscScalar v;

 31:   PetscFunctionBeginUser;
 32:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 33:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 34:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 36:   /* create stiffness matrix C = [1 2; 2 3] */
 37:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 38:   PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N));
 39:   PetscCall(MatSetFromOptions(C));
 40:   PetscCall(MatSetUp(C));
 41:   if (rank == 0) {
 42:     rowidx = 0;
 43:     colidx = 0;
 44:     v      = 1.0;
 45:     PetscCall(MatSetValues(C, 1, &rowidx, 1, &colidx, &v, INSERT_VALUES));
 46:     rowidx = 0;
 47:     colidx = 1;
 48:     v      = 2.0;
 49:     PetscCall(MatSetValues(C, 1, &rowidx, 1, &colidx, &v, INSERT_VALUES));

 51:     rowidx = 1;
 52:     colidx = 0;
 53:     v      = 2.0;
 54:     PetscCall(MatSetValues(C, 1, &rowidx, 1, &colidx, &v, INSERT_VALUES));
 55:     rowidx = 1;
 56:     colidx = 1;
 57:     v      = 3.0;
 58:     PetscCall(MatSetValues(C, 1, &rowidx, 1, &colidx, &v, INSERT_VALUES));
 59:   }
 60:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 61:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

 63:   /* create right-hand side and solution */
 64:   PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
 65:   PetscCall(VecSetSizes(u, PETSC_DECIDE, N));
 66:   PetscCall(VecSetFromOptions(u));
 67:   PetscCall(VecDuplicate(u, &b));
 68:   PetscCall(VecDuplicate(u, &r));
 69:   PetscCall(VecSet(u, 0.0));
 70:   PetscCall(VecSet(b, 1.0));

 72:   /* solve linear system C*u = b */
 73:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 74:   PetscCall(KSPSetOperators(ksp, C, C));
 75:   PetscCall(KSPSetFromOptions(ksp));
 76:   PetscCall(KSPSolve(ksp, b, u));

 78:   /* check residual r = C*u - b */
 79:   PetscCall(MatMult(C, u, r));
 80:   PetscCall(VecAXPY(r, -1.0, b));
 81:   PetscCall(VecNorm(r, NORM_2, &norm));
 82:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "|| C*u - b|| = %g\n", (double)norm));

 84:   /* solve C^T*u = b twice */
 85:   PetscCall(KSPSolveTranspose(ksp, b, u));
 86:   /* check residual r = C^T*u - b */
 87:   PetscCall(MatMultTranspose(C, u, r));
 88:   PetscCall(VecAXPY(r, -1.0, b));
 89:   PetscCall(VecNorm(r, NORM_2, &norm));
 90:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "|| C^T*u - b|| =  %g\n", (double)norm));

 92:   PetscCall(KSPSolveTranspose(ksp, b, u));
 93:   PetscCall(MatMultTranspose(C, u, r));
 94:   PetscCall(VecAXPY(r, -1.0, b));
 95:   PetscCall(VecNorm(r, NORM_2, &norm));
 96:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "|| C^T*u - b|| =  %g\n", (double)norm));

 98:   /* solve C*u = b again */
 99:   PetscCall(KSPSolve(ksp, b, u));
100:   PetscCall(MatMult(C, u, r));
101:   PetscCall(VecAXPY(r, -1.0, b));
102:   PetscCall(VecNorm(r, NORM_2, &norm));
103:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "|| C*u - b|| = %g\n", (double)norm));

105:   PetscCall(KSPDestroy(&ksp));
106:   PetscCall(VecDestroy(&u));
107:   PetscCall(VecDestroy(&r));
108:   PetscCall(VecDestroy(&b));
109:   PetscCall(MatDestroy(&C));
110:   PetscCall(PetscFinalize());
111:   return 0;
112: }