Comments (6)
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.
from blasfeo.
Thanks, @giaf!
from blasfeo.
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.
from blasfeo.
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.
from blasfeo.
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.
from blasfeo.
That should have been more obvious to me! Thank you very much for the information.
from blasfeo.
Related Issues (20)
- Need to link against math library? HOT 1
- Linker error: SHF_MERGE section size (456) must be a multiple of sh_entsize (32) HOT 4
- Tests fail to build: libblasfeo.so: undefined reference to kernel_dpack_buffer_fn HOT 1
- Are there routines for matrix norms? HOT 2
- Incorrect documentation for dtrmm in blasfeo_d_blasfeo_api.h? HOT 3
- blasfeo_dtrmm_rltn not implemented HOT 1
- Missing symbols kernel_dpack_buffer_* in the shared library HOT 2
- Tests fail: error: undefined symbol: blasfeo_sgemm HOT 1
- When can we use parameter as both input and output? HOT 1
- What are m, n, k in dgemm routines? HOT 1
- blasfeo_target.h:1:0: error: unterminated #ifndef HOT 2
- MacOS M2 compiling issue HOT 2
- Calling certain triangular matrices routines leads to `undefined symbol` error HOT 2
- `blasfeo_dtrmm_rlnn` accesses invalid memory with offset on the lower triangular matrix HOT 1
- How can I build BLASFEO and HPIPM for microcontrollers? HOT 1
- what is panel format? HOT 2
- blasfeo_strmm_rlnn fails on HASWELL
- sgemm_nn fails for k > 8 on HASWELL HOT 1
- BLASFEO API explanation for Riccati recursion in HPIPM
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from blasfeo.