Actual source code: ex2.c


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

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

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

 20: #include <petscksp.h>

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

 33:   PetscInitialize(&argc, &args, (char *)0, help);
 34:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 35:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

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

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

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

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

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

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

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

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

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