sorted compressed sparse row matrix

[moortgatgmail_com]

Hi, I've been struggling for days with the following: in my finite element routine a huge sparse matrix is constructed in triplet/coordinate form. I used to solve this system with the old umfpack linear solver. Now I want to switch to the Pardiso solver included with the MKL. However, Pardiso needs the input in the SORTED compressed sparse row (CSR) format. I have implemented themkl_dcsrcoo converter, which seems to work fine. Unfortunately, this routine doesn't sort the column indices though. The fastest way to achieve this is apparently a linear bucket sort, which is equivalent to the transpose of the sparse matrix. I can't find anything like this in the MKL though. There are simple and fast routines available in CSparse for matlab/c, but I am unable to translate those to fortran.

Any help on how to program this, or pointers to existing routines in MKL that I could use would be GREATLY appreciated.

[Gennady_F_Intel]

Hi,
1.MKL contains such routine internally but this routine is not declared externally.
2. As a temporarily workaround for sorting the column indexes you can usemkl_dcsradd(...) or mkl_dcsrmultcsr(...) routines ( see the description of these functions into Sparse BLAS Level 2 and Level 3 Routines.)
Please pay attention on parameter "sort" which specifies the type of reordering.
namely:
sort: If this parameter is not set (default), the routine does not perform reordering. If sort=1, the routine arranges the column indices ja for each row in the increasing order and reorders the corresponding values of the matrix A in the array a. If sort=2, the routine arranges the column indices jb for each row in the increasing order and reorders the corresponding values of the matrix B in the array b.If sort=3, the routine performs reordering for both input matrices A and B.
3. if you are really interested in this functionality, I would recommend you submit the issue against MKL to Premier support( https://premier.intel.com/ ) as a Feature Reguest.
-Gennady

[moortgatgmail_com]

Thanks so much. I wish I had checked back sooner (I thought I'd probably get an email notification; probably something I should have set up). Certainly, it would be useful if the sort-option was an option of mkl_dcsr_coo as well. In the meantime, I had written a fortran routine to do a fast bucket sort on the output sparse matrix, based on the transpose algorithm as implemented in CSparse. I will just paste it below in case it might be useful to others browsing the forum with related issues.

--Joachim

axout = 0

ww = 0; axtmp = 0; aiout = 0; ajout = 0; ajperm=0; ajtmp =0; aitmp=0

job(1) = 1 ! This setting indicates that we want to convert a coordinate/triplet based sparse

! matrix to a compressed CSR based format (however, unfortunately with unsorted column indices).

job(2:3) = 1 ! Both input and output matrices are 1-based (not 0-bases as in c).

job(6)=1 ! Somehow the routine only seems to work when proceeding in 2 steps; this step only finds the compressed AI

! Next we will use the supplied conversion routine provided in the Intel Math Kernel Libaray (MKL) with syntax:

! call mkl_dcsrcoo(job, n, acsr, ja, ia, nnz, acoo, rowind, colind, info)

call mkl_dcsrcoo (job, n, AXout, ajout,aiout,nnz, AX, Ai,Aj,info)

job(6)=2 ! Now, find AX and AJ in CRS form [this shouldn't change the arrays actually and can probably be omitted; as this routine

! is only called once at time = 0, we'll leave it in for now.]

call mkl_dcsrcoo (job, n, AXout, ajout,aiout,nnz, AX, Ai,Aj,info)

! The remainer of this routine is to sort the column indeces by transposing the sparse matrix (and back).

! Traverse the array with column numbers once to find

! the nr of non-zero entries in each column (rather than rows).

aitmp = 1

DO i = 1, nnz

ww(ajout(i)) = ww(ajout(i)) + 1

enddo

! Now it's trivial to find the compressed AI list for the A_transpose from the cumulative sumes of ww:

do i = 2, n+1

aitmp(i) = aitmp(i-1) + ww(i-1)

enddo

! Set the work-space ww equal to AI. In tranposing A, each entry of ww will be incremented

! each time it's encountered.

ww = aitmp

! Transpose once

PI = 0

DO i = 1, n ! Traverse each row i

DO l = aitmp(i), aitmp(i+1)-1 ! l runs through all the columns in row i

pi = ww(ajout(l)) ! next we will place the row index i in column j at position pi = w[..]

ww(ajout(l)) = ww(ajout(l)) + 1 ! which is post-incremented to prepare for the next entry to be inserted into column w[..]

IF(PI > NE) PRINT*, "Error in pardiso solver. PI> NE:", PI, NE

ajtmp(pi) = i

axtmp(pi) = Ax(l)

ENDDO

Enddo

! Transpose back [should be re-written as subroutine for ellegance ]

PI = 0

ww = aiout

DO i = 1, n

DO l = aiout(i), aiout(i+1)-1

pi = ww(ajtmp(l))

ww(ajtmp(l)) = ww(ajtmp(l)) + 1

ajout(pi) = i

axout(pi) = Axtmp(l)

ENDDO

Enddo

! The structure of the sparse matrix A will not change during this simulation

! [it only depends on the geometry/mesh used, which defines which edge-faces are at what positions in A]

! Therefore, once we know the original order of row- and column indeces in AI, AJ, and AX, and the

! CSR sorted arrays, we can obtain a permutation array that can take array elements from the list as

! provided in the main finite-element matrix contruction loop, above, and permute/ put them at the required,

! sorted positions to speed up the Pardiso solver:

! Find permutation array:

DO i = 1, n

DO l = aiout(i), aiout(i+1)-1

thiscol = ajout(l)

DO k = aiout(i), aiout(i+1)-1

IF(thiscol == AJ(k)) ajperm(l) = k

ENDDO

ENDDO

Enddo