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

Sun Aug 5 20:07:52 CDT 2007

  Looks ok to me. Did you compare the right hand side and solution for the two 
codes at that iteration?

On Mon, 6 Aug 2007, Ben Tay 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.
> > > > > > > > > > > 
> > > > > > > > > > > 
> > > > > > > > > > > 
> > > > > > > > > > >                                                             
> > > > > > > > > >                                                   
> > > > > > > > >                                         
> > > > > > > >                                 
> > > > > > >                         
> > > > > >                   
> > > > >             
> > > >         
> > >     
> > 
> > 
> >   
> 
> 