Comparing between a PETSc matrix and a standard fortran array in compressed row format
Barry Smith
bsmith at mcs.anl.gov
Sun Aug 5 10:58:58 CDT 2007
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.
> > > > > > > > >
> > > > > > > > >
> > > > > > > > >
> > > > > > > > >
> > > > > > > >
> > > > > > >
> > > > > >
> > > > >
> > > >
> > >
> >
> >
> >
>
>
More information about the petsc-users
mailing list