SSOR preconditioner using RCI

[Meysam_J_]

Hey Guys,

I have a problem with impelementing SSOR precondtioner using RCI and hope you can help me with that.

I have written a cg solver using RCI for a symmetric matrix that is stored in CSR format 3 array version. I have also implemented the jacobi preconditioner and it works fine. Now I want to implement SSOR but unfortunately I am lost somehow. There is an example of SSOR in MKL example directory, but in that example the SSOR is implemented in a loop which I don't know why!

Here is the SSOR that I am aware of.

Let's say M is the preconditioner which works as follows:

M rho = r   (eq.1)

where r is the current residual ( stored in &tmp[2n] in RCI) and rho is the modified residual (shall be stored in &tmp[3n] in RCI)

M is defined as follow [SAAD 2002]

M = 1/(w(2-w)) (D+wL) D^-1 (D+wU)

where

w = relaxation parameter 0<w<2

U = is the upper part of the matrix A

L = the lower part of the matrix A

D = diagonal of the matrix A

so A = L+D+U and since A is symmetric L= U'

To solve the eq. 1 the following is suggested 

Q = 1/(w(2-w)) (D+wL) D^-1

G = (D+wU)

1. Solve Q z = r

2. Solve G rho = z

To perform the above steps we need to triangular solvers. Now my question, what function from mkl should I use to this end? Is it necessary to build Q and G explicitly or the mkl functions can do so?

I would be happy if somebody can give me some hints here.

With best regards.

Meysam

[Meysam_J_]

Meysam J. wrote:

Hey Guys,

I have a problem with impelementing SSOR precondtioner using RCI and hope you can help me with that.

I have written a cg solver using RCI for a symmetric matrix that is stored in CSR format 3 array version. I have also implemented the jacobi preconditioner and it works fine. Now I want to implement SSOR but unfortunately I am lost somehow. There is an example of SSOR in MKL example directory, but in that example the SSOR is implemented in a loop which I don't know why!

Here is the SSOR that I am aware of.

Let's say M is the preconditioner which works as follows:

M rho = r   (eq.1)

where r is the current residual ( stored in &tmp[2n] in RCI) and rho is the modified residual (shall be stored in &tmp[3n] in RCI)

M is defined as follow [SAAD 2002]

M = 1/(w(2-w)) (D+wL) D^-1 (D+wU)

where

w = relaxation parameter 0<w<2

U = is the upper part of the matrix A

L = the lower part of the matrix A

D = diagonal of the matrix A

so A = L+D+U and since A is symmetric L= U'

To solve the eq. 1 the following is suggested 

Q = 1/(w(2-w)) (D+wL) D^-1

G = (D+wU)

1. Solve Q z = r

2. Solve G rho = z

To perform the above steps we need to triangular solvers. Now my question, what function from mkl should I use to this end? Is it necessary to build Q and G explicitly or the mkl functions can do so?

I would be happy if somebody can give me some hints here.

With best regards.

Meysam

Meysam

[Meysam_J_]

Hi Guys agian,

Let me ask my question as the following:

Is there any function in mkl to solve the following equation system:

(D+wL) x = y

where:

D = diagonal of a CSR matrix

L = Lower part of a CSR matrix

w = scalar

thanks, Meysam

[Zhang_Z_Intel]

Can you take a look at the MKL sample code $MKLROOT//examples/solverc/source/cg_ssor_precon_c.c, which shows one implementation of SSOR?

To solve the triangular sparse matrix, you can use 'mkl_?csrsv'. See the description of this function here: http://software.intel.com/en-us/node/468586

[Sergey_P_Intel2]

Hi Meysam,

Intel MKL SparseBLAS has functionality - *csrsv - which can solve triangular sparse matrices that should be built explicitly in the application. There are no functions which can solve equation in the form:

(D+wL)x=y

Regards, Sergey

[Meysam_J_]

Thank mecej4,

This was exactly what I was looking for.  I think in your code I don't even need to define ic and jc since the resulting matrix C has the same structure as A.

Thank you so much for the help.

Cheers,Meysam

[Meysam_J_]

Hey Zhang Z,

I am already aware of that example. But there is a "for loop" around the preconditioner part which I cannot understand it (niter_ssor!). I wonder if you can point me to some literatures  where this routine has been explained.

Regards,

Meysam

[mecej4]

Despite MKL not having a routine to do the C = diag(A) + w B operation, this is easy to implement in Fortran or C, as long as there are no missing diagonal entries in the CSR representations of A and B. I have attached a small example where A and B have different numbers of non-zero entries, and are square matrices. The code does not need to be changed if the matrices are triangular, since the CSC representation is not bothered by where the non-zero elements are located.

If there are missing diagonal entries, the task is more complicated, but it is clear that it can be done using existing MKL routines augmented with similar code to the example that I have provided.

[mecej4]

Meysam J. wrote:

...the resulting matrix C has the same structure as A B.

Note the correction A -> B. It is possible that A has different off-diagonal entries than B.