Inversion of matrix of dimension greater than 2^15.5

[qtl]

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.

[qtl]

Anyhow, using block matrix inversion can circumvent this. But a one step function is preferred.

[SergeyKostrov]

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

[Gennady_F_Intel]

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...

[yuriisig]

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

[basel]

>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

[mecej4]

> 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.

[basel]

>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

[qtl]

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.

[yuriisig]

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.

[qtl]

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.cpp
When 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.

[yuriisig]

MKL_INTc(4294967295) ???
Dynamic selection of storage is necessary.

[yuriisig]

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

[Gennady_F_Intel]

Please add the all ILP64 libraries explicitly. See the MKL Linker adviser here.