[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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

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.