Comparing between a PETSc matrix and a standard fortran array in compressed row format
Matthew Knepley
knepley at gmail.com
Sun Aug 5 20:12:45 CDT 2007
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 at gmail.com> 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.
> >>>>>>>>>>
> >>>>>>>>>>
> >>>>>>>>>>
> >>>>>>>>>>
> >>>>>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>
> >>>>>>>
> >>>>>>
> >>>>>>
> >>>>>
> >>>>>
> >>>>
> >>>>
> >>>
> >>>
> >>
> >
> >
> >
>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which
their experiments lead.
-- Norbert Wiener
More information about the petsc-users
mailing list