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 <mailingli...@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 >
