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