Comparing between a PETSc matrix and a standard fortran array in compressed row format
Ben Tay
zonexo at gmail.com
Sun Aug 5 20:04:37 CDT 2007
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