[petsc-users] Matrix missing diagonal entry

[Florian Lindner]

Am Dienstag, 7. Oktober 2014, 08:34:55 schrieb Barry Smith:
> 
>    PETSc has this “extra” requirement that you provide values for the entire diagonal of the matrix, even those locations with zero. So make sure that you use MatSetValues() to put a zero on each diagonal entry that has zero.

Thanks for that information. I came up with this scheme:

  _matCLU.assemble(MAT_FLUSH_ASSEMBLY);
  petsc::Vector zeros(_matCLU);   // construct a vector like the given matrix
  MatDiagonalSet(_matCLU.matrix, zeros.vector, ADD_VALUES);
  _matCLU.assemble(MAT_FINAL_ASSEMBLY);

(contains a bit of my petsc helper functions but should be clear I guess)

To set all unset entries to zero. Do you see any problem with that?

Best Regards,
Florian

> 
>   Also from below you are using ICC that is incomplete cholesky factorization. If you want to use a direct solver you should select -pc_type cholesky
> 
> 
>    Barry
> 
> On Oct 7, 2014, at 8:14 AM, Florian Lindner <mailinglists at xgm.de> wrote:
> 
> > Hello,
> > 
> > when I try to KSPSolve my matrix petsc prints out that a diagonal entry is missing:
> > 
> > KSPSetOperators(solver, _matCLU.matrix, _matCLU.matrix );
> > KSPSolve(solver, vin.vector, vp.vector);
> > 
> > And petsc is perfectly right about that specific entry being zero. So this may be more of a mathematical problem then pure petsc.
> > 
> > The algorithm that I try to reimplement does the following:
> > 
> >  lu(_matrixCLU, _pivotsCLU);  // Compute LU decomposition
> >  // Left out: Permutate in according to pivotsCLU
> >  forwardSubstitution(_matrixCLU, in, y); // CLU^-1 * in = y  (lower triangle of CLU)
> >  backSubstitution(_matrixCLU, y, p);     // CLU^-1 * y = p (upper triangle of CLU)
> >  multiply(_matrixA, p, out );            // out = A * p (not relevant here)
> > 
> > I want to replace to lu, forward- and backwardSubstituation by using a petsc KSP solver. Before doing any of these operations _matrixCLU and _matCLU are equal:
> > 
> > =============== CLU Before LU =======
> > ==== Dynamic Matrix of size 7 x 7 ====
> > 1.0000 0.3679 0.1353 0.3679 1.0000 0.0000 0.0000
> > 0.3679 1.0000 0.3679 0.1353 1.0000 1.0000 0.0000
> > 0.1353 0.3679 1.0000 0.3679 1.0000 1.0000 1.0000
> > 0.3679 0.1353 0.3679 1.0000 1.0000 0.0000 1.0000
> > 1.0000 1.0000 1.0000 1.0000 0.0000 0.0000 0.0000
> > 0.0000 1.0000 1.0000 0.0000 0.0000 0.0000 0.0000
> > 0.0000 0.0000 1.0000 1.0000 0.0000 0.0000 0.0000
> > 
> > Mat Object:CLU 1 MPI processes
> >  type: seqsbaij
> > row 0: (0, 1)  (1, 0.367879)  (2, 0.135335)  (3, 0.367879)  (4, 1)  (5, 0)  (6, 0) 
> > row 1: (1, 1)  (2, 0.367879)  (3, 0.135335)  (4, 1)  (5, 1)  (6, 0) 
> > row 2: (2, 1)  (3, 0.367879)  (4, 1)  (5, 1)  (6, 1) 
> > row 3: (3, 1)  (4, 1)  (5, 0)  (6, 1) 
> > row 4:
> > row 5:
> > row 6:
> > 
> > LU decomposition of matrixCLU works, no rank deficiency (no 0 on the main diagonal after LU and I tested for full rank of matrixCLU). Petsc complains with this matrix:
> > 
> > [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> > [0]PETSC ERROR: Object is in wrong state
> > [0]PETSC ERROR: Matrix is missing diagonal entry 4
> > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
> > [0]PETSC ERROR: Petsc Release Version 3.5.2, unknown 
> > [0]PETSC ERROR: ../build/debug/binprecice on a arch-linux2-c-debug named helium by lindnefn Mon Oct  6 16:26:49 2014
> > [0]PETSC ERROR: Configure options 
> > [0]PETSC ERROR: #1 MatICCFactorSymbolic_SeqSBAIJ() line 2430 in /data2/scratch/lindner/petsc/src/mat/impls/sbaij/seq/sbaijfact2.c
> > [0]PETSC ERROR: #2 MatICCFactorSymbolic() line 6266 in /data2/scratch/lindner/petsc/src/mat/interface/matrix.c
> > [0]PETSC ERROR: #3 PCSetup_ICC() line 20 in /data2/scratch/lindner/petsc/src/ksp/pc/impls/factor/icc/icc.c
> > [0]PETSC ERROR: #4 PCSetUp() line 902 in /data2/scratch/lindner/petsc/src/ksp/pc/interface/precon.c
> > [0]PETSC ERROR: #5 KSPSetUp() line 305 in /data2/scratch/lindner/petsc/src/ksp/ksp/interface/itfunc.c
> > [0]PETSC ERROR: #6 KSPSolve() line 417 in /data2/scratch/lindner/petsc/src/ksp/ksp/interface/itfunc.c
> > [0]PETSC ERROR: #7 map() line 696 in src/mapping/PetRadialBasisFctMapping.hpp
> >  numerical equality test failed: 0.00000000000000000000e+00 instead of 1.00000000000000000000e+00
> >  file: src/mapping/tests/PetRadialBasisFctMappingTest.cpp       line: 196
> >  statement: value==1.0
> >  difference:  1.00000000000000000000e+00
> > 
> > 
> > What's indeed perfectly right. But what is the best way to cope with that?
> > 
> > Thanks,
> > Florian
>