topic Quote:Meysam J. wrote: in Intel® oneAPI Math Kernel Library

SSOR preconditioner using RCI

Meysam_J_ — Tue, 26 Nov 2013 14:36:23 GMT

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

Quote:Meysam J. wrote:

Meysam_J_ — Tue, 26 Nov 2013 14:38:00 GMT

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

Hi Guys agian,

Meysam_J_ — Wed, 27 Nov 2013 10:25:34 GMT

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

Can you take a look at the

Zhang_Z_Intel — Wed, 27 Nov 2013 21:17:27 GMT

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

Hi Meysam,

Sergey_P_Intel2 — Thu, 28 Nov 2013 09:59:49 GMT

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

Thank mecej4,

Meysam_J_ — Thu, 28 Nov 2013 13:25:32 GMT

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

Hey Zhang Z,

Meysam_J_ — Thu, 28 Nov 2013 13:30:22 GMT

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

Despite MKL not having a

mecej4 — Thu, 28 Nov 2013 13:34:34 GMT

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.

Quote:Meysam J. wrote:

mecej4 — Thu, 28 Nov 2013 13:41:00 GMT

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.