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

[Ben Tay]

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