forums.silverfrost.com

[mecej4]

There are, as far as I know, no special routines in MKL for sparse matrices where the nonzeros are confined to recognizable blocks. From matrix theory, we know that many of the properties of matrices of real-number elements also apply to matrix (i.e., block) elements.

I have not had much experience with block matrices. Have you read the book "Matrix Computations" by Golub and Van Loan? They cover such topics very well.

People who work with splines have developed solvers for ABD (Almost Block Diagonal) matrices.

MKL includes a version of Pardiso. It also contains a DSS interface to Pardiso. I think that it is worth trying out Pardiso on some of your block-sparse matrices before looking for special block solvers. If you send me a couple of big (say, 10K X 10K) matrices in COO or CSR format, I'll try out Pardiso on them and compare with Laipe timings. That should help you decide if Pardiso will suffice for your purposes.

[DanRRight]

I suggest to use for the matrix A and vector B the same approach like it was used with dense matrix case and fill them just with the random numbers. The code in this case is almost the same too. This will create two squares on the major diagonal like was plotted above

[mecej4]

With four threads on my i7-2720QM laptop, I get:

[JohnCampbell]

Dan,

For your block matrix, I would expect that a skyline approach would be effective, as the operations count would be the same as for any block specific approach.

The disadvantage of a skyline solver is that they typically do not have any pivoting capability. The skyline solvers I am familiar with are for symmetric matrices, and apply artificial constraints to (near) zero diagonals, which is not a problem for the Finite Element matrices I am using.

[DanRRight]

Broke again my favorite glasses with UV coating (too fragile), but good that this pushed me to find lost 3 months ago (!) spare ones without coating. The coating specifically made for monitors, it clearly changes tint, that means it starts screening even some far blue. I think it is important to have that coating if we sit 10+ hours in front of larger and larger screens, some of which emit dangerous UV spectrum transforming it to blue which may leak enough to do harm over large time (not clear how much monitor front glass leaks UV, better to do clean experiment and measure the spectrum of LED and OLED monitors)

Mecej4, many thanks again, you have done great job and I owe you now even more. Have you noticed that while MKL dense matrix solver consumes almost 100% of CPU time in Task Manager the sparse one grabs only ~40% ? LAIPE2 dense case grabs even less, why it is slower. Sparse LAIPE2 though get closer to 100%

Anyway, here are preliminary test results (in seconds) for the block sparse matrix. MKL wins on larger size matrices >3000 but loses on smaller ones (processor I7-4770K 4.5Ghz)

[mecej4]

Dan, I think you will find it interesting to do a log-log plot of solution time against N, perhaps after removing data that took less than 0.1 s to run.

For full matrices, one expects O(N^3) for Gaussian elimination. Your Laipe2 results give about O(N^2.85), whereas MKL gives O(N^1.45). The difference in the exponents is striking. Perhaps that is not a fair comparison, because Pardiso has been asked to do equation reordering, but Laipe2 was not asked to do the same.

If your ultimate goal is to solve block matrix equations with N = 10^4 or more, Laipe2 may be too slow.

[DanRRight]

That scaling is indeed good.

[mecej4]

The value in iparm(2) controls reordering. You cannot turn reordering off in Pardiso -- you can choose which reordering algorithm to use.

When you factorise a banded matrix the factors have the same bandwidths as the original matrix, but the fill-in within the band is affected by reordering. You should probably try Metis reordering, and see if asking Laipe to reorder helps to increase its speed.

[JohnCampbell]

Dan,

If you are using MKL with reordering, you need to be careful defining the initial matrix, especially the % non-zero within the profile.
It is interesting that MKL will attempt re-ordering of a general sparse matrix on a per equation basis. I would have thought that to be unproductive, but the MKL results are impressive. (How does it manage the unknown increase in storage required?)
Lower % cpu can be an indicator that there is not good load sharing between the threads, which can be due to the specific sparsity characteristics of the test matrix. You may want to experiment with different %non-zero values within the profile.
If you can define the profile of the test matrices, I could run them in my profile solver and report how they compare.

[DanRRight]

John, Just take simple 2-block matrix as I plotted above first (matrix A in the source below). Later you can make 3 or more such blocks and then takes them of different size.
Here is the source for the Laipe VAG_8 I used in comparisons above. Compilation:

[DanRRight]

John,
The code above is a bit messy, so I will probably repeat how the matrix was set by simpler way

[JohnCampbell]

Dan,

Trying to review the solver times, I stumbled at the first hurdle, as you have used laipe$Decompose_VAG_8, while I can only perform laipe$Decompose_VSG_8 solutions.
With Variable-Bandwidth and Asymmetric Systems, I presume the solver first checks that you have provided for non-zero infill during reduction.
I haven't used MKL, so I presume non-zero infill is managed there also. Does the MKL sparse input format require this allowance ?
In the 32-bit days, allocating sufficient temporary space could be a significant issue ?

Are all your tests assuming non-symmetric matrices ?

[DanRRight]

Blocks inside the matrix are squares which go along the major diagonal, so the geometrically the shape of matrix is symmetric, but matrix elements all are different, so formally it is not symmetric (not equal to transposed as a definition of symmetric).

Dense 20000x20000 real*8 matrix contains 3GB of data. My block sparse one 1/2 of that

[DanRRight]

Here is the test with single precision using LAIPE block sparse solver. I substituted VAG_8 by VAG_4 and all arrays with real*4. It shows twice the speed of VAG_8 above

[mecej4]

If by "older LAIPE circa 1997 made by the Microsoft Fortran" you mean a library made for use with Fortran Powerstation (FPS) 4, to call routines in the library you need to use the STDCALL convention. In the caller, you may need to add declarations for this purpose, or use compiler options, depending on whether your code is compiled with IFort or FTN95.

Using the FPS-4 compiler, I ran the code that I posted in my post of March 22, 2017 after changing the subroutine names to LAIPE-1 names and providing interface blocks for the three LAIPE subroutines that are called. The run times were roughly the same as with Intel Fortran, because most of the run time is spent in the LAIPE-1 library.