Comparing between a PETSc matrix and a standard fortran array in compressed row format
Barry Smith
bsmith at mcs.anl.gov
Sun Aug 12 21:12:50 CDT 2007
On Mon, 13 Aug 2007, Ben Tay wrote:
> 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?
Ben,
Are you running with -info and looking at what the line search does?
(Are you using SNES?) does the line search fail? You'll need to see
what NSPC does at this point.
*) same rhs
*) same matrix
*) same solution (?) to some good number of digits
NOW what does NSPCG do to update the Newton step and how and why is
that update different then what PETSC does?
Barry
>
> 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 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.
> > > > > > > > > > > > >
> > > > > > > > > > > > >
> > > > > > > > > > > > >
> > > > > > > > > > > > >
> > > > > > > > > > > > >
> > > > > > > > > > > > >
> > > > > > > > > > > >
> > > > > > > > > > >
> > > > > > > > > >
> > > > > > > > >
> > > > > > > >
> > > > > > >
> > > > > >
> > > >
> > > >
> > >
> >
> >
> >
>
>
More information about the petsc-users
mailing list