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

[Barry Smith]

   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.
> > > > > 
> > > > > 
> > > > > 
> > > > >             
> > > >         
> > >     
> > 
> > 
> >   
> 
>