Iterative Solver Internal Settings

[josefin]

Hi,

I'm doing ice sheet modeling, and am modifying/adding stuff to some old version of a FlowSolver routine, (so where Au=b is solved, describing the Stokes equations).

I'm solving a dual problem, i.e. I'm solving A^T x=f, where A^T is the transpose of A. This works fine if I'm using a direct solver, but if I use an iterative solver (gcr, BigStabl) it immediately crashes (the norm in the first iteration is NaN). Is there any setting in the iterative solver that assumes that the matrix I'm solving with is the Solver % Matrix, and if so, how do I change it? This is my code:

A => Solver % Matrix
AT => AllocateMatrix() !Will be transponate of A
AT % Format = MATRIX_LIST

.. Some allocations...

!Transpose A
AT = A
AT = CRS_Transpose(AT)

!Allocate and set the right hand side
ALLOCATE(AT % RHS(SIZE(A % RHS)))
AT % RHS = f
...
!Set system matrix to AT, and values to x
Solver % Matrix => AT
A => Solver % Matrix
Solver % Variable % Values => x
...
Solver % Variable % Values=0._dp
UNorm = DefaultSolve()

Cheers,

Josefin Ahlkrona

[Franz Pichler]

HI there,

i am not quite sure but isnt the AT=A giving the pointer A (which is pointing to solver%matrix) into AT so that effectively AT is pointing to solver%matrix? And the stuff that you have allocated for AT is floating without a pointer?

best regards
Franz

[josefin]

Hi Franz,

You're right I probably pointed back and forth a big too much I could also just set AT = CRS_Transpose(A) directly (although I do need to allocate AT first)

I don't believe this is the problem though, since it works perfectly fine to solve with a direct solver :/

Best,

Josefin

[josefin]

It must have something to do with the CRS_Transpose, since if I exhange

AT = A
AT = CRS_Transpose(AT)

in the code I posted originally to just

AT = A

it works, but if I change it to

AT = A
AT = CRS_Transpose(AT)
AT = CRS_Transpose(AT)

It does not work. I also tried using CALL List_toCRSMatrix(AT) before transposing but it doesn't change. In the function CRS_Transpose routine the only suspcisious thing I can see is that it sets AT % Diag = 0.

Any suggestions anyone? The output I get when I run is:

IterSolver: Using iterative method: bicgstabl
CRS_IncompleteLU: ILU(0) (Real), Starting Factorization:
CRS_IncompleteLU: Allocated LU matrix of size: 3174080
CRS_IncompleteLU: Performing incomplete LU
CRS_IncompleteLU: ILU(0) (Real), NOF nonzeros: 3174081
CRS_IncompleteLU: ILU(0) (Real), filling (%) : 99
CRS_IncompleteLU: ILU(0) (Real), Factorization ready at (s): 0.02
ERROR:: RealBiCGStab(l): Breakdown error: sigma == NaN.

Cheers,

Josefin

[raback]

Hi Josefin

You could perhaps try to call CRS_SortMatrix just after transposing it. It sets the Diag properly, if that is the problem.

Is there any difference depending in the type of preconditioning (ILUn, none, diag, vanka, etc.). I might guess that some preconditioners take use of the diagonal entries, or assume matrix symmetry in topology, while the direct solver probably don't use these.

-Peter

[josefin]

Hej Peter,

Yey! Thanks a lot, using CRS_SortMatrix afterwards solved it!

Cheers,

Josefin