[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Comparing between a PETSc matrix and a standard fortran array in compressed row format

Hi,

I have installed the latest ver of PETSc to make sure it's not due to any bug. I only read in the matrix into PETSc just before the eqn is to be solved. Another way which I've done is to convert the matrix to a CSR format (verified to be correct) and use MatCreateSeqAIJWithArrays to read in the matrix. Unfortunately, the answer is still the same. Is there any other recommendation?

Thanks

Matthew Knepley wrote: 
If

  1) The matrix and rhs are identical

  2) The matrix is nonsingular

then the solution will be identical. Thus, the code has a bug somewhere.
I suggest

  1) Read in  the other matrix and rhs

  2) Compare both matrices and rhs

  3) Solve both systems

  4) Check that the answer match

    Matt

On 8/5/07, Ben Tay <zonexo@xxxxxxxxx> wrote:
Well, here's what I got:

0 KSP preconditioned resid norm 1.498836640002e+04 true resid norm
1.531062859147e+03 ||Ae||/||Ax|| 9.947622397959e-01
1 KSP preconditioned resid norm 4.832736106865e-08 true resid norm
2.037181386899e-09 ||Ae||/||Ax|| 1.323597595746e-12

Looks like everything's ok.... right?

Any other suggestion?

Thanks.

Barry Smith wrote:
 Run with -ksp_type gmres -pc_type lu -ksp_monitor_true_residual -ksp_truemonitor

  At the bad iteration verify if the true residual is very small.

   Barry


On Sun, 5 Aug 2007, Ben Tay wrote:


Well, ya, they're the same ... that's why I don't know what's wrong. Is there
any way to check on singularity? The iteration no. is only 8 for AMG and LU
direct solver works. ..


Thanks

Barry Smith wrote:

  At the "bad step" are the two matrices and right hand sides the same?
The the resulting solution vectors are different (a lot)? Is there any
reason to think the matrix is (nearly) singular at that point?

   Barry



On Sun, 5 Aug 2007, Ben Tay wrote:



Hi,

Now I'm comparing at every time step. The A matrix and rhs of the NSPCG
and
PETSc matrix are exactly the same, NORM is zero. The solutions given by
both
solvers are very similar for e.g. from 1-315 time step, and varying slowly
since it is a oscillating body problem.  Then strangely, at time=316,
PETSc's
answer suddenly differs by a great difference. e.g. p=0.3 to 50. I must
also
add that the flowfield at that time step still "looks ok". However, the
large
pressure change of the pressure poisson eqn caused the computation of the
lift/drag coefficient to be erroneous. Moreover, after that, the pressure
slowly "returns back" to the same answer such that of the NSPCG solver
after a
few time steps ... then after a while the sudden change happens again. In
the
end, I got a cl/cd vs time graph with lots of spikes.

I hope that the thing can be solved because the NSPCG solver is much
slower
and it is not as stable.

Thanks

Barry Smith wrote:


   Ben,

    Are you comparing the two matrices at timestep 316? Also compare the
two right hand sides at 316. Are they similar, totally different??

   Barry


On Sat, 4 Aug 2007, Ben Tay wrote:



Hi,

After using MatAXPY, I used MatGetRowMaxAbs, VecAb and VecMax to get
the
row
with max value, get the absolute of that value and print out the
location/value. I managed to find some minor mistakes and corrected
them
so
that the 2 matrix are identical. Unfortunately, the same thing still
happens!

I've checked for convergence and the iteration no. is 1 and 8 for
direct
solving and HYPRE's AMG respectively.  So there's no convergence
problem.
I've
also tried to change to the same pc/ksp as that of NSPCG, which is
jacobi
presconditioning + KSPBCGS but the same thing happens.

Any idea what's happening?

Thanks

Barry Smith wrote:


  Unless the matrix values are huge then this is a big difference.
All
you
can do is print the difference matrix with ASCII MatView and
look for large entries. This will tell you the locations of the
trouble
entries.

   Barry


On Sat, 4 Aug 2007, Ben Tay wrote:



Hi,

I realised that the matrix storage format of NSPCG is actually
ITPACK
storage
format. Anyway, I tried to insert the values into a PETSc matrix,
use
MatAXPY
with the scalar a = -1 and then use MatNorm on the output matrix.

Using NORM_1, NORM_FROBENIUS and NORM_INFINITY
<../Vec/NORM_FROBENIUS.html#NORM_FROBENIUS>, I got ~543, 3194 and
556
respectively. Does this mean that the matrix are different? How
"much"
different does that means?

So what is the next step? Is there anyway to find out the location
of
the
value which is different?

This is my comparison subroutine:

subroutine matrix_compare

   !compare big_A & PETSc matrix

   integer :: k,kk,II,JJ,ierr

   Vec    xx_test,b_rhs_test
     Mat    A_mat_test

   PetscScalar aa

   PetscReal nrm

   aa=-1.

   call
MatCreateSeqAIJ(PETSC_COMM_SELF,total_k,total_k,10,PETSC_NULL_INTEGER,A_mat_test,ierr)

   call VecCreateSeq(PETSC_COMM_SELF,total_k,b_rhs_test,ierr)

   call VecDuplicate(b_rhs_test,xx_test,ierr)

   call VecAssemblyBegin(b_rhs_test,ierr)

   call VecAssemblyEnd(b_rhs_test,ierr)

   call VecAssemblyBegin(xx_test,ierr)

   call VecAssemblyEnd(xx_test,ierr)

       do k=1,total_k
         do kk=1,10

           II=k-1

           JJ=int_A(k,kk)-1
                 call
MatSetValues(A_mat_test,1,II,1,JJ,big_A(k,kk),INSERT_VALUES,ierr)
call VecSetValue(b_rhs_test,II,q_p(k),INSERT_VALUES,ierr)
             end do

   end do

!    call MatCopy(A_mat,A_mat_test,DIFFERENT_NONZERO_PATTERN
,ierr)

   call MatAssemblyBegin(A_mat_test,MAT_FINAL_ASSEMBLY,ierr)
       call MatAssemblyEnd(A_mat_test,MAT_FINAL_ASSEMBLY,ierr)

   call MatAXPY(A_mat_test,aa,A_mat,SAME_NONZERO_PATTERN,ierr)

   call MatNorm(A_mat_test,NORM_INFINITY,nrm,ierr)

   print *, nrm

   end subroutine matrix_compare

Thanks!



Barry Smith wrote:


  You can use MatCreateSeqAIJWithArrays() to "convert" the NSPCG
matrix
into
PETSc format and then MatAXPY() to difference them and then
MatNorm() to
see
how large the result is.   Make sure the PETSc or hypre solver
is
always
converging. Run with
-ksp_converged_reason and or -ksp_monitor. My guess is that the
matrix
is
becoming very ill-conditioned so the solvers with the default
options
are
failing.

    Barry


On Fri, 3 Aug 2007, Ben Tay wrote:



Hi,

I used 2 packages to solve my poisson eqn, which is part of my
NS
unsteady
solver. One is the NSPCG solver package which uses the
compressed
row
format
to store the A matrix. The other is PETSc. I found that using
both
solvers
gave me similar answers for a number of time step. However,
after
that,
the
answer will suddenly change drastically for the PETSc case.
This
does
not
happen for the NSPCG solver.

For e.g. time step 1-315, oscillating airfoil case, pressure
changes
smoothly,
similar answers in both cases

at time=316, pressure changes from -3.22 to -3.2 for NSPCG,
but
pressure
changes from -3.21 to -60.2 for PETSc

This happens when I use HYPRE's AMG or PETSc's direct solver
LU.

I have been trying to find out what's the cause and I can't
find
the
answer in
debugging. I would like to compare the values of the matrix of
the
2
different
solvers and see if there's any difference. However, NSPCG's
matrix
is
in
compressed row format while PETSc's one is just an address and
it
can't be
viewed easily. Moreover, it's a big matrix so it's not
possible to
check
by
inspection. I'm thinking of subtracting one matrix by the
other
and
find
if
it's zero. What's the best way to solve this problem? Btw, I'm
using
fortran
and there's no mpi

Thank you.