LU-Factorization differs in Lapack and Matlab

Question

I'm computing a LU-Factorization with Lapack and compare the results. I use either dgetrf or dgetf2. Comparing the results to the corresponding function lu(...) in Matlab, they are identical as long as the the input matrix is square or has more columns than rows. For matrixes with more rows than columns, the values at the bottom right lu-matrix differ. For example, dgetrf and dgetf2 in Lapack for

A=[
    1  2  3;
    4  5  6;
    7  8  9;
   10 11 12] 

gives

lu = [
   10.0 11.0    12.0; 
    0.1  0.9     1.8;
    0.7  0.33333 0.0;
    0.4  0.66666 0.45989304812834225]

but the same calculated in Matlab results in

lu = [
   10.0 11.0    12.0; 
    0.1  0.9     1.8;
    0.7  0.33333 0.0;
    0.4  0.66666 0.0]

I have no clue, where difference in the last value of the last row comes from. The same effect is present, when having even more rows. I'm using Lapack in Swift, but I don't think it has anything to do with the programming language. Any clue? Your help is appreciated.

I've tried varying the input parameters of dgetrf and dgetf2, but still the same effect.

Answer 1

I wrote some MATLAB m-code to test this on my win PC running R2021a. I also wrote C-code to call the MATLAB LAPACK library directly for dgetrf( ) and dgetf2( ). Here is what I get for results:

>> lu_test
A =
     1     2     3
     4     5     6
     7     8     9
    10    11    12
LU =
   10.0000   11.0000   12.0000
    0.1000    0.9000    1.8000
    0.7000    0.3333    0.0000
    0.4000    0.6667    0.4599
LUdgetrf =
   10.0000   11.0000   12.0000
    0.1000    0.9000    1.8000
    0.7000    0.3333    0.0000
    0.4000    0.6667    0.4599
LUdgetf2 =
   10.0000   11.0000   12.0000
    0.1000    0.9000    1.8000
    0.7000    0.3333    0.0000
    0.4000    0.6667    0.4599

So the MATLAB lu( ) function produces the same result as calling the MATLAB LAPACK library functions dgetrf( ) and dgetf2( ) directly, and they all have that lower right corner spot filled. Maybe you can compare your code to my code below to help determine what is going wrong at your end.

The m-code:

A=[ 1  2  3;
    4  5  6;
    7  8  9;
    10 11 12] 
LU = lu(A)
LUdgetrf = dgetrf(A)
LUdgetf2 = dgetf2(A)

The C-code for dgetrf:

/* MATLAB mex routine, Calls LAPACK dgetrf routine for LU factoring     */
/* Programmer: James Tursa                                              */

/* Includes ----------------------------------------------------------- */

#include "mex.h"
#include "blas.h"
#include "lapack.h"

/* Gateway ------------------------------------------------------------ */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSignedIndex M, N, LDA, INFO, MN;
    mwSignedIndex *IPIV;
    mxArray *A;
    
    if( nlhs > 3 ) {
        mexErrMsgTxt("Too many outputs.");
    }
    if( nrhs != 1 || !mxIsDouble(prhs[0]) || mxIsSparse(prhs[0]) || mxIsComplex(prhs[0])
                  || mxGetNumberOfDimensions(prhs[0]) != 2 ) {
        mexErrMsgTxt("Need one full real double 2D matrix input.");
    }
    M = mxGetM(prhs[0]);
    N = mxGetN(prhs[0]);
    LDA = M;
    MN = M<N ? M : N;
    IPIV = mxMalloc(MN * sizeof(*IPIV));
    A = mxDuplicateArray(prhs[0]);
    dgetrf( &M, &N, mxGetData(A), &LDA, IPIV, &INFO );
    if( INFO ) {
        mexErrMsgTxt("dgetrf returned an error code.");
    }
    if( nrhs <= 1 ) {
        plhs[0] = A;
    } else if( nrhs == 2 ) {
        mexErrMsgTxt("2 outputs not yet implemented.");
    } else {
        mexErrMsgTxt("3 outputs not yet implemented.");
    }
    mxFree(IPIV);
}

The C-code for dgetf2:

/* MATLAB mex routine, Calls LAPACK dgetf2 routine for LU factoring     */
/* Programmer: James Tursa                                              */

/* Includes ----------------------------------------------------------- */

#include "mex.h"
#include "blas.h"
#include "lapack.h"

/* Gateway ------------------------------------------------------------ */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSignedIndex M, N, LDA, INFO, MN;
    mwSignedIndex *IPIV;
    mxArray *A;
    
    if( nlhs > 3 ) {
        mexErrMsgTxt("Too many outputs.");
    }
    if( nrhs != 1 || !mxIsDouble(prhs[0]) || mxIsSparse(prhs[0]) || mxIsComplex(prhs[0])
                  || mxGetNumberOfDimensions(prhs[0]) != 2 ) {
        mexErrMsgTxt("Need one full real double 2D matrix input.");
    }
    M = mxGetM(prhs[0]);
    N = mxGetN(prhs[0]);
    LDA = M;
    MN = M<N ? M : N;
    IPIV = mxMalloc(MN * sizeof(*IPIV));
    A = mxDuplicateArray(prhs[0]);
    dgetf2( &M, &N, mxGetData(A), &LDA, IPIV, &INFO );
    if( INFO ) {
        mexErrMsgTxt("dgetf2 returned an error code.");
    }
    if( nrhs <= 1 ) {
        plhs[0] = A;
    } else if( nrhs == 2 ) {
        mexErrMsgTxt("2 outputs not yet implemented.");
    } else {
        mexErrMsgTxt("3 outputs not yet implemented.");
    }
    mxFree(IPIV);
}

The mex compile routine:

function compile_lapack(varargin)
libdir = 'microsoft';
lib_blas = fullfile(matlabroot,'extern','lib',computer('arch'),libdir,'libmwblas.lib');
lib_lapack = fullfile(matlabroot,'extern','lib',computer('arch'),libdir,'libmwlapack.lib');
mex(varargin{:},lib_blas,lib_lapack,'-largeArrayDims');
end

To compile the mex routines:

>> compile_lapack('dgetf2.c')
Building with 'MinGW64 Compiler (C)'.
MEX completed successfully.

>> compile_lapack('dgetrf.c')
Building with 'MinGW64 Compiler (C)'.
MEX completed successfully.

UPDATE -----------------------------------------------------------

Online MATLAB R2023a gives these results:

>> LUA = lu(A)
LUA =
   10.0000   11.0000   12.0000
    0.1000    0.9000    1.8000
    0.7000    0.3333         0
    0.4000    0.6667         0
>> LUA(3,3)
ans =
   0
>> LUA(4,3)
ans =
   0

So there does appear to be a somewhat recent change in the LU algorithm used that gives an exact 0 for the U(3,3) and L(4,3) spots. Whereas my PC MATLAB R2021a gives:

>> LUA(3,3)
ans =
   6.9204e-16
>> LUA(4,3)
ans =
    0.4599

Yes, it looks funny, but both answers are "correct" in the sense that these residual values are within what one might expect from simple floating point numerical effects.

Feedback on the MATLAB Answers forum indicates that this change may have occurred in the R2022a release.