Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- Fast solving of block-tridiagonal matrices

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Mikhail_Matrosov

Novice

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-07-2011
10:11 AM

152 Views

Fast solving of block-tridiagonal matrices

I'm trying to solve SLE of form Ax=b with matrix A being large, sparse, symmetric, positive defined and block-tridiagonal. With exactly the same structure as one arising from five-point finite difference approximation, used for solving Poisson equation. I.e.:

D1 L2

L2 D2 L3

L3 D3 L4

L4 D4 L5 ...

where Di is symmetric tridiagonal and Li is diagonal. When solving Poisson equation, Li = -I with I denoting unit matrix and Di is of form

4 -1

-1 4 -1

-1 4 -1 ...

I tried direct solver (DSS) and iterative solver (PCG) with different preconditioners. I was able to compute solution, but I'm looking for possible speed-ups. I've solved Poisson eqation of the same size and solution was found much faster. E.g., 0.1 seconds against 15 seconds for PCG and 325 seconds for DSS.

The question is, are there some special routines for dealing with matrices of such a structure in MKL? I understand, that Poisson-case is optimised for matrix with given values, while I have arbitrary values in A. But still, are there some alternatives for PCG?

Link Copied

12 Replies

Alexander_K_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-07-2011
07:24 PM

152 Views

The DSS is direct solver, PCG is iterative and Poisson solver is direct solver for matrix with form you wrote in topic above. It's clear that for some matrix each of these 3 solvers spend different time to compute solution and this computational depend on condition number, matrix structure, matrix value and etc... It difficult to say for general matrix which solver is best for it, but if you provide to us matrix you want to solve we could recommended you one of these solvers. For matrix from you topic the best solver is Poisson solver.

With best regards,

Alexander Kalinkin

Mikhail_Matrosov

Novice

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-07-2011
11:59 PM

152 Views

12.7972 -9.4535 0 -2.3438 0 0 0 0 0 0 0 0

-9.4535 39.8238 -26.2896 0 -3.0808 0 0 0 0 0 0 0

0 -26.2896 31.1540 0 0 -3.8644 0 0 0 0 0 0

-2.3438 0 0 103.7057 -33.4790 0 -66.8830 0 0 0 0 0

0 -3.0808 0 -33.4790 66.0553 -11.2585 0 -17.2370 0 0 0 0

0 0 -3.8644 0 -11.2585 93.8882 0 0 -77.7652 0 0 0

0 0 0 -66.8830 0 0 87.0312 -13.1618 0 -5.9864 0 0

0 0 0 0 -17.2370 0 -13.1618 41.9111 -5.2094 0 -5.3030 0

0 0 0 0 0 -77.7652 0 -5.2094 88.3185 0 0 -4.3439

0 0 0 0 0 0 -5.9864 0 0 24.3179 -17.3315 0

0 0 0 0 0 0 0 -5.3030 0 -17.3315 27.9133 -4.2788

0 0 0 0 0 0 0 0 -4.3439 0 -4.2788 9.6227

For my experiments, PCG with incomplete Cholesky factorization as preconditioner works much faster than DSS. Is there a way to use Poisson solver for this matrix? As I understand, values for matrix for Poisson solver are predefined (4, -1, 0, -1, etc...), while I have arbitrary values (with mentioned constraints).

IDZ_A_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-08-2011
01:17 AM

152 Views

Hi,

Have you try PCG with poisson solver as preconditioner?

With best regards,

Alexander Kalinkin

Mikhail_Matrosov

Novice

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-08-2011
01:31 AM

152 Views

f(x, y) stands for b in SLE of the form Ax=b. And u stands for x. However, coefficients of matrix A are determined from finite approximation of second derivative. So how it is possible to provide my own coefficients for A?..

Alexander_K_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-08-2011
01:43 AM

152 Views

Maybe I wrote a bit unclear... Any iterative process for solving Ax=f could be representing in next form:

B(x^(k+1)-x^k)=tau_k*(Ax^k-f), or

x^(k+1)-x^k=B^(-1)*(tau_k*(Ax^k-f)) where B is preconditioner. So you could use PCG with Poisson solver as preconditioner instead of PCG with incomplete Cholesky preconditioner (when you want to multiply preconditioner on some vector just call poisson solver with this array)

With best regards,

Alexander Kalinkin

Mikhail_Matrosov

Novice

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-08-2011
01:54 AM

152 Views

4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0

-1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0

0 -1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0

0 0 -1 4 0 0 0 -1 0 0 0 0 0 0 0 0

-1 0 0 0 4 -1 0 0 -1 0 0 0 0 0 0 0

0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0 0

0 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0

0 0 0 -1 0 0 -1 4 0 0 0 -1 0 0 0 0

0 0 0 0 -1 0 0 0 4 -1 0 0 -1 0 0 0

0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0

0 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0

0 0 0 0 0 0 0 -1 0 0 -1 4 0 0 0 -1

0 0 0 0 0 0 0 0 -1 0 0 0 4 -1 0 0

0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 -1 0

0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 -1

0 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4

Am I right?

Alexander_K_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-08-2011
02:35 AM

152 Views

With best regards,

Alexander Kalinkin

Mikhail_Matrosov

Novice

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-09-2011
02:33 AM

152 Views

So, back to my question, are there some special routines for dealing with matrices of such a structure in MKL?

Alexander_K_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-09-2011
08:47 AM

152 Views

Back to iss solver with Poisson precondition: what kind of boundary condition have you used for your solver? It's seems that for your matrix there is kernel vector (all elements of this vector are equal 1). So probably it's better to use Neumann boundary condition... I will try to check it later.

With best regards,

Alexander Kalinkin

Mikhail_Matrosov

Novice

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-09-2011
09:10 AM

152 Views

I checked convergence of Poisson preconditioner in MATLAB, I simply provided matrix mentioned above as preconditioner to MATLAB's PCG routine. I'm always testing all preconditioners in MATLAB first to swiftly estimate convergence rate - it works properly. As I understand, this is the same as to run Poisson solver with Dirichlet boundary conditions. Don't think Neumann will change anything. Moreover, solution for Neumann boundary conditions is determined up to a constant value - it's not acceptable in my case.

styc

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-11-2011
07:56 PM

152 Views

Mikhail_Matrosov

Novice

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-11-2011
10:09 PM

152 Views

The most narrow place in current PCG is applying preconditioner, not matrix-vector multiplication (time ratio is approx 6:1). To apply incomplete Cholesky preconditioner, two triangular systems are solved using MKL sparse BLAS routine mkl_ddiatrsv. It fits perfectly since Cholesky lower triangle has exactly three non-zero diagonals. Maybe applying the LDL^T factorization will be a bit faster than current LL^T (since L will be unit diagonal), but it is essentially the same.

Matrix-vector multiplication uses CSR and might be boosted using the same diagonal representation, but it will not benefit much on the overall time, since applying the preconditioner is the most expensive operation.

One more thing: applying Jacobi preconditioner is much cheaper, but convergence is too low, so Cholesky is now the best option.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.