Matrix Inversions

[0p3n-s0urc3r3r]

How would I go about determining the inverse of a matrix with the existing API? It seems that DGETRI is not implemented, as suggested by users on Stack Overflow. Are there plans to implement DGETRI?

[giaf]

Hi, dgetri may be implemented at some point but there are no concrete plans for it yet.

Using the current API, once you have an LU factorization (using dgetrf) you can always compute the inverse by solving with an identity matrix B * B_inv = L * U * B_inv = I that gives B_inv = U \ (L \ I) where the backslash means triangular solve with matrix right hand side dtrsm (or dgetrs to compute both solutions at once). The downside is that the flop count is higher than with a dedicated dgetri routine.

[0p3n-s0urc3r3r]

Thanks, @giaf!

[0p3n-s0urc3r3r]

Hi @giaf, I am running into some issues with dgetrf. I have not fully diagnosed the issue, but I suspect it has something to do with panel and column formats.

I have defined the following in my build:

• TARGET_ARMV7A_ARM_CORTEX_A9
• LA_HIGH_PERFORMANCE
• MF_PANELMAJ

While attempting to get the LU factorization of a variety of matrices, I noticed that only 4x4 matrices show expected results.

Here is a minimal example:

#define MATRIX_SIZE_A (4)
static INT aiPIVOT_BUFFER[MATRIX_SIZE_A ];
static DOUBLE adMEMORY_BUFFER_A[(MATRIX_SIZE_A * MATRIX_SIZE_A )];

static const DOUBLE DGETRF_EXAMPLE_A[MATRIX_SIZE_A][MATRIX_SIZE_A] = {
   { 32.7880,  3.9799, 57.3508, 49.3076},
   {  1.7340, 13.3143, 69.3625, 16.5685},
   { 43.3401, 90.6032, 89.8322, 16.4636},
   { 45.2754, 77.2201, 39.4537, 83.0328}
};

ULONG ulMatrixSize_ = MATRIX_SIZE_A;
CHAR cN = 'n';
INT iInfo = 0;

struct blasfeo_dmat stMatrixA;
blasfeo_create_dmat(ulMatrixSize_, ulMatrixSize_, &stMatrixA, (void*)(&adMEMORY_BUFFER_A[0]));
blasfeo_dgese(ulMatrixSize_, ulMatrixSize_, 0.0, &stMatrixA, 0, 0);

for(UINT uiIndexI = 0; uiIndexI < ulMatrixSize_; uiIndexI++)
{
   for(UINT uiIndexJ = 0; uiIndexJ < ulMatrixSize_; uiIndexJ++)
   {
      BLASFEO_DMATEL(&stMatrixA, uiIndexJ, uiIndexI) = DGETRF_EXAMPLE_A[uiIndexJ][uiIndexI];
   }
}

dgetrf_((INT*)(&ulMatrixSize_), (INT*)(&ulMatrixSize_), stMatrixA.pA, (INT*)(&ulMatrixSize_), &aiPIVOT_BUFFER[0], &iInfo);

I have tested a series of matrices, as shown below.

N = 4: Input:

32.788  3.980   57.351  49.308
1.734   13.314  69.362  16.569
43.340  90.603  89.832  16.464
45.275  77.220  39.454  83.033

Result:

45.275  77.220  39.454  83.033
0.724   -51.942 28.779  -10.824
0.038   -0.199  73.590  11.230
0.957   -0.321  0.833   -75.853

Expected:

45.275  77.220  39.454  83.033
0.724   -51.942 28.779  -10.824
0.038   -0.199  73.590  11.230
0.957   -0.321  0.833   -75.853

N = 3: Input:

32.788  3.980   57.351
1.734   13.314  69.362
43.340  90.603  89.832

Result:

43.340  -10.072 13.865
0.757   -0.342  69.362
0.040   57.351  89.832

Expected:

43.340  90.603  89.832
0.757   -64.564 -10.610
0.040   -0.150  64.176

N = 5: Input:

32.788  3.980   57.351  49.308  53.087
1.734   13.314  69.362  16.569  87.710
43.340  90.603  89.832  16.464  17.398
45.275  77.220  39.454  83.033  68.406
82.716  75.683  56.646  67.629  45.323

Result:

45.275  0.957   0.728   90.224  52.420
0.038   57.351  0.252   -0.147  87.166
0.724   88.407  16.569  0.262   41.525
0.088   -0.319  38.819  17.398  0.842
82.716  75.683  56.646  67.629  45.323

Expected:

45.275  77.220  39.454  83.033  68.406
0.038   10.357  67.851  13.389  85.090
0.724   -5.015  369.080 56.327  430.311
1.827   -6.314  1.119   -62.557 -23.880
0.957   1.611   -0.155  1.213   -89.468

If you need more information, I can do my best to provide it. I hope that you can help me determine if something is awry in the implementation of dgetrf_() with my configuration.

[0p3n-s0urc3r3r]

Since I suspected that this issue is actually related to the Panel Format, I switched my configuration to Column-major, and the problem seems to have gone away.

I suppose this could be a bug in the implementation for dgetrf_ in either the memory access for the panel-major format.

[giaf]

Hi,

the issue is the fact that you are mixing different BLASFEO interfaces (APIs).

If you use the standard BLAS and LAPACK API provided by BLASFEO (i.e. dgetrf_), then you should not use the BLASFEO matrix structure blasfeo_dmat (and in particular, not attempt to use its internal memory in the .pA struct member, as the actual data layout in memory in implementation-dependent), but simply pass the pointer to the array of doubles (in column-major order, so the transposed of DGETRF_EXAMPLE_A) as you would do with any other standard BLAS and LAPACK implementation.

If instead you decide to use BLASFEO's own API (i.e. blasfeo_dgetrf_rp), then you have to pass a pointer to the BLASFEO matrix struct blasfeo_dmat, as you would do with any other BLASFEO routine.

The internals of the BLASFEO structs like blasfeo_dmat should never be accessed directly in any case, as they are implementation-dependent; the provided functions and macros should be used instead. For example, a different implementation of the macro BLASFEO_DMATEL exists for each different data layout in memory, so you can safely access the matrix elements even without knowing the internal matrix format.

[0p3n-s0urc3r3r]

That should have been more obvious to me! Thank you very much for the information.