Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
7036 Discussions

## Inversion of matrix of dimension greater than 2^15.5

Beginner
1,092 Views
Hi, MKL matrix inversions are based on LAPACK and store matrices in vectors. Their dimension is thus limited by the integer size, namely sqrt(2^31). Is it possible to inverse a matrix greater than that?

PS. in animal breeding, large covariance matrix inversion is often used. Before, we are using the relationship expectation and its inversion has an O(n) algorithm. Now we are using realized relationship due to the advances of molecular genetics. Then we have to inverse such matrices with brute force.

At the moment, matirx dim=40k, dgetrf+dgetri take about 560min user time on Intel x5550. This is acceptable, so is their accuracy. If using dpotrf+dpotri, the time can be further reduced to 275 min.
13 Replies
Beginner
1,092 Views
Anyhow, using block matrix inversion can circumvent this. But a one step function is preferred.
Valued Contributor II
1,092 Views

I've checked IPP's Matrix Processing API:

- There are lots of functions to calculate a matrix inverse, like:

IppStatus ippmInvert_m_32f(
...
Ipp32u widthHeight,
... );

- A declaration for 'Ipp32u', is as follows:

typedef unsigned int Ipp32u;

- A max value for 'unsigned int' is as follows:

0xFFFFFFFF(base16) = 4294967295(Base10) = ((2^32) - 1)(Base10)

- And, there is another declaration:

#define IPP_MAX_32U ( 0xFFFFFFFF )

- A maximum size of a matrix could be 4294967295 x 4294967295

Moderator
1,092 Views
one note should be added specifically for IPP vector math function implementation - all of these functions are highly optimized for 2x2, 3x3... 6x6 matrixes. So, execution the ippmInvertr_m_32f function for the big inputs will take very long time...
Beginner
1,092 Views

Intel MKL uses instead of int variable MKL_INT. For MKL ILP64 integer type is MKL_INT64.

Beginner
1,092 Views
>PS. in animal breeding, large covariance matrix inversion is often used. Before, we are using the relationship expectation and its >inversion has an O(n) algorithm. Now we are using realized relationship due to the advances of molecular genetics. Then we have >to inverse such matrices with brute force.

Is the matrix sparse and symmetric?

Olaf
Honored Contributor III
1,092 Views
> Is the matrix sparse and symmetric?

It is symmetric since multiplication is commutative. Often, elements that should be zero show up with small non-zero values because of noise in the data from which the covariance matrix is computed.
Beginner
1,092 Views

>It is symmetric since multiplication is commutative. Often, elements that should be zero show up with small non-zero values >because of noise in the data from which the covariance matrix is computed.

If it is symmtric and sparse you might use other options than the LAPACK routines. We are computing diagonal elements of the inverse in A in an animal breeding applications with millions of equations. Send me an email and we can discuss it offline.

olaf.schenk@unibas.ch
Beginner
1,092 Views
With the realized relationship, the matrix is not sparse, all values are in [0, 1]. But the matrix is symmetric. The matrix dimension is acturally the number of ID in a population, ie number of animals/plant individuals. So relationship between ID A & B is same to relationship of B&A.

I also realized that large matrix needs really `super' computers, e.g., even using half storage, 1M id need 3.7T memory.
Beginner
1,092 Views

You read my message? There it is written that you can work with matrixes of dimension greater than 2^15.5 on one computer, using Intel MKL.

Beginner
1,092 Views
I wrote a sample program as below
```[cpp]#include
#include
using namespace std;
int main(int argc, char *argv[])
{
MKL_INT c(4294967295), inc(1), i;
double a;

for(i=0; i=.5;
cout< Then I compiled with the following command:icpc -DMKL_ILP64 -mkl t.cppWhen running the binary, I got a segment fault error. My ulimit stack size is 6G. So I think the BLAS core is not ready for int64.```
Beginner
1,092 Views
MKL_INTc(4294967295) ???
Dynamic selection of storage is necessary.
Beginner
1,092 Views
Quoting qtl
With the realized relationship, the matrix is not sparse, all values are in [0, 1]. But the matrix is symmetric.
Intel MKL is able to process the symmetric matrixes in the packed form. It saves storage. But these algorithms in Intel MKL are ineffective: http://software.intel.com/en-us/forums/showthread.php?t=76595&o=d&s=lr
Moderator
1,092 Views
Please add the all ILP64 libraries explicitly. See the MKL Linker adviser here.