Actual source code: ex2.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      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;
     26:   PetscInt       N = 2,rowidx,colidx;
     27:   Vec            u,b,r;
     28:   KSP            ksp;
     29:   PetscReal      norm;
     30:   PetscMPIInt    rank,size;
     31:   PetscScalar    v;

     33:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
     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) {
     43:     rowidx = 0; colidx = 0; v = 1.0;
     44:     MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
     45:     rowidx = 0; colidx = 1; v = 2.0;
     46:     MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);

     48:     rowidx = 1; colidx = 0; v = 2.0;
     49:     MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
     50:     rowidx = 1; colidx = 1; v = 3.0;
     51:     MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
     52:   }
     53:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
     54:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

     56:   /* create right hand side and solution */
     57:   VecCreate(PETSC_COMM_WORLD,&u);
     58:   VecSetSizes(u,PETSC_DECIDE,N);
     59:   VecSetFromOptions(u);
     60:   VecDuplicate(u,&b);
     61:   VecDuplicate(u,&r);
     62:   VecSet(u,0.0);
     63:   VecSet(b,1.0);

     65:   /* solve linear system C*u = b */
     66:   KSPCreate(PETSC_COMM_WORLD,&ksp);
     67:   KSPSetOperators(ksp,C,C);
     68:   KSPSetFromOptions(ksp);
     69:   KSPSolve(ksp,b,u);

     71:   /* check residual r = C*u - b */
     72:   MatMult(C,u,r);
     73:   VecAXPY(r,-1.0,b);
     74:   VecNorm(r,NORM_2,&norm);
     75:   PetscPrintf(PETSC_COMM_WORLD,"|| C*u - b|| = %g\n",(double)norm);

     77:   /* solve C^T*u = b twice */
     78:   KSPSolveTranspose(ksp,b,u);
     79:   /* check residual r = C^T*u - b */
     80:   MatMultTranspose(C,u,r);
     81:   VecAXPY(r,-1.0,b);
     82:   VecNorm(r,NORM_2,&norm);
     83:   PetscPrintf(PETSC_COMM_WORLD,"|| C^T*u - b|| =  %g\n",(double)norm);

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

     91:   /* solve C*u = b again */
     92:   KSPSolve(ksp,b,u);
     93:   MatMult(C,u,r);
     94:   VecAXPY(r,-1.0,b);
     95:   VecNorm(r,NORM_2,&norm);
     96:   PetscPrintf(PETSC_COMM_WORLD,"|| C*u - b|| = %g\n",(double)norm);

     98:   KSPDestroy(&ksp);
     99:   VecDestroy(&u);
    100:   VecDestroy(&r);
    101:   VecDestroy(&b);
    102:   MatDestroy(&C);
    103:   PetscFinalize();
    104:   return ierr;
    105: }