Sparse Matrix Storage

[pourmatin85]

Dear guys;

I'm trying to use PARDISO routine insteade of ?GETRS to use the advantage of sparse matrix. However, I couldn't find any MKL routine for converting a matrix to the proper sparse storage.

Is there any library that I can use for sparse storage, or I should write my own code for it?

Thanks in advance
Hossein

[ArturGuzik]

Quoting - pourmatin85

Dear guys;

I'm trying to use PARDISO routine insteade of ?GETRS to use the advantage of sparse matrix. However, I couldn't find any MKL routine for converting a matrix to the proper sparse storage.

Is there any library that I can use for sparse storage, or I should write my own code for it?

Thanks in advance
Hossein

there is. see description of mkl_dcsrcoo routine. also this thread and Gennady's response for additional info.

A.

[pourmatin85]

Quoting - ArturGuzik

there is. see description of mkl_dcsrcoo routine. also this thread and Gennady's response for additional info.

A.

Thanks alot for your quick reply. I found mkl_dcsrcoo and learned to use it. However, I just found out that my actual problem is with the matrix at the beginning of the execution. It means that for a heavy input data (for example 15000 elements) the code would allocate a, say, 20000 by 20000 matrix, which will led to memory stack error.

So, seems to me I have to generate CSR at the very first beginning of the code so i would never use such a huge matrix.

Is there any routine which can help me in this case?

Regards
Hossein

[Gennady_F_Intel]

Hossein,
It seems that your input matrix is not sparse matrix but the dense one.
The number of non-zeros elements for sparse matrix is pretty small - several percentage.
So that is the reason why the sparse Solvers work with very big system and 20000x20000 is the typical sizes.
What is the number of non-zeros elements in your input matrix and what is the size of this matrix?
-Gennady

[pourmatin85]

Quoting - Gennady Fedorov (Intel)

Hossein,
It seems that your input matrix is not sparse matrix but the dense one.
The number of non-zeros elements for sparse matrix is pretty small - several percentage.
So that is the reason why the sparse Solvers work with very big system and 20000x20000 is the typical sizes.
What is the number of non-zeros elements in your input matrix and what is the size of this matrix?
-Gennady

I can assure you that my matrix is a sparse matrix, for the number of non-zero elements are 2% to 5% of total elements.
However, I'm very surperised about what you said about 20000x20000 matrixes, is it really a typical size to you!!
anyway, is there any routine that could avoid allocation of big matrixes. If not, what can I do with the memory stack errors?

[ArturGuzik]

Quoting - pourmatin85

I can assure you that my matrix is a sparse matrix, for the number of non-zero elements are 2% to 5% of total elements.
However, I'm very surperised about what you said about 20000x20000 matrixes, is it really a typical size to you!!
anyway, is there any routine that could avoid allocation of big matrixes. If not, what can I do with the memory stack errors?

Hossein,

before trying to suggest something, a few things (additional info) need to be clarified:

(1) what architecture (IA32/x64) and what system (win/linux) are you on?
(2) what is the 'nature' of your problem to solve? Is it FEM or something else?
(3) how do you allocate matrix (COMMON block, ALLOCATABLE statement) and how do you pass the matrix to subroutines, if at all?

Comments:
(1) stack (for windows exe) can be increased.
(2) matrices passed incorrectly require compiler (sometimes) to make temp copy
(3) you can use heap arrays switch
(4) anyway you should try to avoid creating rectangular matrix just for the sake of making a sparse (compressed) version at the next step.

A.

[pourmatin85]

I have a 64 achitecture intel processor (T7700) and vista (32-bit). However, my code works in win32 platform (becuase of of problems that I had with including mkl for 64x). Anyway my program is FEM and I allocate my stiffness matrix and I pass it through a madule.

[ArturGuzik]

Quoting - pourmatin85

I have a 64 achitecture intel processor (T7700) and vista (32-bit). However, my code works in win32 platform (becuase of of problems that I had with including mkl for 64x). Anyway my program is FEM and I allocate my stiffness matrix and I pass it through a madule.

OK then. You not only have a spare but also a symmetric (stiffness) matrix. For IA32 with 20,000 x 20,000 (double?) real matrix you just reach the (theoretical) 2GB (you may try boot with 3GB option) per process limits. So basically you have two options at hand:

(1) fix problem once for all and create sparse matrix (avoiding the huge matrix) only. In case of FEM is not that difficult, as you need something like:

(a) loop over all nodes
(b) for each node determine column and row indexes
(c) write routine for assembling the matrix (symm part only).

Try to find something (source code) on Internet as this was done in FEM for ages.

(2) switch to x64 system (what problems did you have????, MKL works very well on x64 system) and try to go along existing path of assembling matrix and then using mkl_xxx routine to create compressed format matrix.


A.

[pourmatin85]

Ok, here is my problem about x64 system:
when I start to build my project in x64, this error appears:
Error 1 fatal error LNK1112: module machine type 'IA64' conflicts with target machine type 'x64' mkl_intel_ilp64.lib(_dgetrf.obj)

Am I including wrong libraries?

[ArturGuzik]

Quoting - pourmatin85

Ok, here is my problem about x64 system:
when I start to build my project in x64, this error appears:
Error 1 fatal error LNK1112: module machine type 'IA64' conflicts with target machine type 'x64' mkl_intel_ilp64.lib(_dgetrf.obj)

Am I including wrong libraries?

yes, you do (use wrong libs). IA64 means Itanium processor. Your architecture is known as em64t or Intel64, so you need to point to something like C:Program Files (x86)IntelCompiler11.0�72fortranmklem64tlib.

Please take a look also at this thread (end of it).

A.

[pourmatin85]

Quoting - ArturGuzik

yes, you do (use wrong libs). IA64 means Itanium processor. Your architecture is known as em64t or Intel64, so you need to point to something like C:Program Files (x86)IntelCompiler11.0�72fortranmklem64tlib.

Please take a look also at this thread (end of it).

A.

Thanks alot A.
You were right. My problem with x64 platform was just solved. However, now I cannot use MKL_DCSRCOO correctly. When I call it, it just does not do anything and doesn't show any error either! (info=0)
For the record, I want to convert the coordinate format sparse that I've made into CSR so I can use the output as an input for PRADISO

Regards
Hossein

[ArturGuzik]

Quoting - pourmatin85

Thanks alot A.
You were right. My problem with x64 platform was just solved. However, now I cannot use MKL_DCSRCOO correctly. When I call it, it just does not do anything and doesn't show any error either! (info=0)
For the record, I want to convert the coordinate format sparse that I've made into CSR so I can use the output as an input for PRADISO

Regards
Hossein

This one I don't know. If the code is working on IA32 system, then the (very wild) guess would be that there is something wrong with the data type (say, integer/real size), (did you modify any default VS setting related to that?). Show snippet of the code, so guys from Intel can take a look and offer some advice.

A.

[pourmatin85]

Quoting - ArturGuzik

This one I don't know. If the code is working on IA32 system, then the (very wild) guess would be that there is something wrong with the data type (say, integer/real size), (did you modify any default VS setting related to that?). Show snippet of the code, so guys from Intel can take a look and offer some advice.

A.

No, it didn't work in win32 either. The problem is I cannot find any reliable resource about it and MKL manual hasn't mentioned it at all.
However, here is my code:
INTEGER INFO, JOB(8), IA(SIZE(GSTIFF,1)+1), NZZ, TEMP_ROW(SIZE(GSTIFF)), TEMP_COL(SIZE(GSTIFF))
INTEGER, ALLOCATABLE :: JA(:), STIFF_ROW(:), STIFF_COL(:)
REAL*8, ALLOCATABLE :: ACSR(:), STIFF_VALUE(:)
REAL*8 TEMP_VALUE(SIZE(GSTIFF))
.
.
.
I=1
DO INODE=1,SIZE(GSTIFF,1)
DO JNODE=1,SIZE(GSTIFF,1)
IF (ABS(GSTIFF(INODE,JNODE))>1.D-4) THEN
TEMP_VALUE(I)=GSTIFF(INODE,JNODE)
TEMP_ROW(I)=INODE
TEMP_COL(I)=JNODE
I=I+1
ENDIF
ENDDO
ENDDO
I=I-1
ALLOCATE(STIFF_VALUE(I), STIFF_ROW(I), STIFF_COL(I))
STIFF_VALUE=TEMP_VALUE(1:I)
STIFF_ROW=TEMP_ROW(1:I)
STIFF_COL=TEMP_COL(1:I)
JOB(5)=I
DATA (JOB(1)=1, JOB(2)=1, JOB(3)=1, JOB(6)=1 )
CALL MKL_DCSRCOO(JOB, SIZE(GSTIFF,1), ACSR, JA, IA, NNZ, STIFF_VALUE, STIFF_ROW, STIFF_COL, INFO)

[ArturGuzik]

Quoting - pourmatin85

No, it didn't work in win32 either. The problem is I cannot find any reliable resource about it and MKL manual hasn't mentioned it at all.
However, here is my code:
INTEGER INFO, JOB(8), IA(SIZE(GSTIFF,1)+1), NZZ, TEMP_ROW(SIZE(GSTIFF)), TEMP_COL(SIZE(GSTIFF))
INTEGER, ALLOCATABLE :: JA(:), STIFF_ROW(:), STIFF_COL(:)
REAL*8, ALLOCATABLE :: ACSR(:), STIFF_VALUE(:)
REAL*8 TEMP_VALUE(SIZE(GSTIFF))
.
.
.
I=1
DO INODE=1,SIZE(GSTIFF,1)
DO JNODE=1,SIZE(GSTIFF,1)
IF (ABS(GSTIFF(INODE,JNODE))>1.D-4) THEN
TEMP_VALUE(I)=GSTIFF(INODE,JNODE)
TEMP_ROW(I)=INODE
TEMP_COL(I)=JNODE
I=I+1
ENDIF
ENDDO
ENDDO
I=I-1
ALLOCATE(STIFF_VALUE(I), STIFF_ROW(I), STIFF_COL(I))
STIFF_VALUE=TEMP_VALUE(1:I)
STIFF_ROW=TEMP_ROW(1:I)
STIFF_COL=TEMP_COL(1:I)
JOB(5)=I
DATA (JOB(1)=1, JOB(2)=1, JOB(3)=1, JOB(6)=1 )
CALL MKL_DCSRCOO(JOB, SIZE(GSTIFF,1), ACSR, JA, IA, NNZ, STIFF_VALUE, STIFF_ROW, STIFF_COL, INFO)

I will take a look, however, doesn't see how you get GSTIFF.

A.

[ArturGuzik]

Quoting - ArturGuzik

I will take a look, however, doesn't see how you get GSTIFF.

A.

I must say, that I'm VERY frequently wondering who writes/maintains the MKL docs. Really. Trying to decipher the docs and MKL_dcsrcoo i found that:

(1) seems to be bug (?) in Interface, correct me if I'm wrong, but where "m" is defined there? Should be n

[cpp]      INTERFACE
      subroutine mkl_dcsrcoo ( job, n, Acsr, AJR, AIR, nnz              &
     &Acoo, ir, jc, info)
      Integer            job(8)
      integer            n, nnz, info 
      integer 		     AJR(*), AIR(m+1), ir(*), jc(*)
      double precision   Acsr(*), Acoo(*)
      END
      END INTERFACE
[/cpp]

(2) it seems that in handling request is also a bug or rather documentation is out of date/touch with reality.

If JOB(6) = 0, which seems to get JA, IA, ACSR on output, it does nothing (INFO =0) but complains about input parameter 1 (JOB). However, if called twice, first time with JOB(6) = 1 and then JOB(6) = 2, then all is fine.

So, call routine twice, and then you'll get what you want.

A.

[Gennady_F_Intel]

Quoting - pourmatin85

No, it didn't work in win32 either. The problem is I cannot find any reliable resource about it and MKL manual hasn't mentioned it at all.
However, here is my code:
INTEGER INFO, JOB(8), IA(SIZE(GSTIFF,1)+1), NZZ, TEMP_ROW(SIZE(GSTIFF)), TEMP_COL(SIZE(GSTIFF))
INTEGER, ALLOCATABLE :: JA(:), STIFF_ROW(:), STIFF_COL(:)
REAL*8, ALLOCATABLE :: ACSR(:), STIFF_VALUE(:)
REAL*8 TEMP_VALUE(SIZE(GSTIFF))
.
.
.
I=1
DO INODE=1,SIZE(GSTIFF,1)
DO JNODE=1,SIZE(GSTIFF,1)
IF (ABS(GSTIFF(INODE,JNODE))>1.D-4) THEN
TEMP_VALUE(I)=GSTIFF(INODE,JNODE)
TEMP_ROW(I)=INODE
TEMP_COL(I)=JNODE
I=I+1
ENDIF
ENDDO
ENDDO
I=I-1
ALLOCATE(STIFF_VALUE(I), STIFF_ROW(I), STIFF_COL(I))
STIFF_VALUE=TEMP_VALUE(1:I)
STIFF_ROW=TEMP_ROW(1:I)
STIFF_COL=TEMP_COL(1:I)
JOB(5)=I
DATA (JOB(1)=1, JOB(2)=1, JOB(3)=1, JOB(6)=1 )
CALL MKL_DCSRCOO(JOB, SIZE(GSTIFF,1), ACSR, JA, IA, NNZ, STIFF_VALUE, STIFF_ROW, STIFF_COL, INFO)

Hossein,
"No, it didn't work in win32 either. The problem is I cannot find any reliable resource about it and MKL manual hasn't mentioned it at all."

There are two examples program for using MKL Sparse format converters for "C" and "Fortran".
You can find these files into directory: ..examplesspblassourceconverters_f.f and converters_c.c.
Both of these examples are used the following routines:
MKL_DDNSCSR, MKL_DCSRCOO, MKL_DCSRBSR, MKL_DCSRDIA,MKL_DCSRCSC andMKL_DCSRSKY
I hope it will help you.
--Gennady

[ArturGuzik]

Quoting - Gennady Fedorov (Intel)

There are two examples program for using MKL Sparse format converters for "C" and "Fortran".
You can find these files into directory: ..examplesspblassourceconverters_f.f and converters_c.c.
Both of these examples are used the following routines:
MKL_DDNSCSR, MKL_DCSRCOO, MKL_DCSRBSR, MKL_DCSRDIA,MKL_DCSRCSC andMKL_DCSRSKY
I hope it will help you.
--Gennady

Gennady,

what about Interface (include file)? Is that correct?

And what about documentation? I mean the job(6)=0? The examples you mentioned, are calling routine in a sequence (of tasks), and then all is fine (see my previous post).


A.