Community
cancel
Showing results for 
Search instead for 
Did you mean: 
qtl
Beginner
157 Views

Inversion of matrix of dimension greater than 2^15.5

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.
0 Kudos
13 Replies
qtl
Beginner
157 Views

Anyhow, using block matrix inversion can circumvent this. But a one step function is preferred.
SergeyKostrov
Valued Contributor II
157 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

Gennady_F_Intel
Moderator
157 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...
yuriisig
Beginner
157 Views

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

basel
Beginner
157 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
mecej4
Black Belt
157 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.
basel
Beginner
157 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
qtl
Beginner
157 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.
yuriisig
Beginner
157 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.

qtl
Beginner
157 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.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
Beginner
157 Views

MKL_INTc(4294967295) ???
Dynamic selection of storage is necessary.
yuriisig
Beginner
157 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
Gennady_F_Intel
Moderator
157 Views

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