topic I see that in your C++ caller in Intel® oneAPI Math Kernel Library

dgelss does not work

Hai — Thu, 17 Jan 2013 16:48:38 GMT

Hi,

I am working on function "dgelss" to simulate pinv in Matlab. After I called dgelss() in MKL (lapack), the output info will be 0 which means execution is successful. But the result matrix does not provide the correct value. Can someone help me figure out why?

MKL version: 10.2.2.025

int main(int argc, char *argv[])
{
   Init();
   int step_extended, step_orig;
   int unitMatrixSize = 1;
   IppStatus status;
   int m = 3, n = 3, nrhs = 1;
   float *matrixA = (float *)malloc(sizeof(float)*m*n);

    //Data Matrix, size is 3x3
   matrixA[0] =   11.0;
   matrixA[1] =   9.0;
   matrixA[2] =   -6.0;
   matrixA[3] =   3.0;
   matrixA[4] =   -3.0;
   matrixA[5] =   0.0;
   matrixA[6] =   -3.0;
   matrixA[7] =   0.0;
   matrixA[8] =   1.0;

    //identity matrix, size is 3x3
   float * unitMatrix = (float *)malloc(sizeof(float) * 9);
   for(int i = 0; i < 3; i++)
       for(int j = 0; j < 3; j++)
           unitMatrix[j*m+i] = 0.0;

   unitMatrix[0] = 1.0;
   unitMatrix[4] = 1.0;
   unitMatrix[8] = 1.0;

   int lda = m;
   int ldb = m;
   double work;
   int rank = -1, info, lwork=-1;
   float rcond = .01;
   int dim = min(m,n);
   float *S = (float *)malloc(sizeof(float) *dim);
   float *WORK = (float *)malloc(sizeof(float) *1);

   sgelss(&m, &n, &m, matrixA, &lda, unitMatrix,&ldb, S, &rcond, &rank, WORK, &lwork, &info);
   free(matrixA);
   free(unitMatrix);
   free(S);
   free(WORK);

   return TRUE;
}

Quote:But the result matrix

mecej4 — Thu, 17 Jan 2013 19:38:00 GMT

But the result matrix does not provide the correct value.

On what basis do you rest this claim? The matrix you gave is of full rank. Therefore, the minimum norm solution is the same as the result of Gaussian elimination, and the minimum residuals are zero. Are you aware that you are calling the Fortran77 routine and, therefore, the input matrices are expected to be in column major order?

The matrix, as entered in your program, is

11 3 -3

9 -3 0

-6 0 1

Is this the matrix you wanted, or its transpose?

thanks for your reply first.

Hai — Thu, 17 Jan 2013 20:39:33 GMT

thanks for your reply first.

What I want to do is to find the psuedo inverse of a given matrix. The size of the matrix is arbitrary. In my code, I use 3x3 as an example. For the matrix B in ?gelss, I used a identity matrix. Then the result should be the psuedo inverse of input matrix A. I did not notice about the column major order issue. But right now, I just cannot get even a look-good result.

For example, if I use Matlab to compute, if the input matrix A is:

[11 3 -3

9 -3 0

-6 0 1]

Then the result matrix B should be pinv(A)

     [0.5000    0.5000    1.5000
    1.5000    1.1667    4.5000
    3.0000    3.0000   10.0000]

thanks

mecej4 wrote:

Quote:

But the result matrix does not provide the correct value.

On what basis do you rest this claim? The matrix you gave is of full rank. Therefore, the minimum norm solution is the same as the result of Gaussian elimination, and the minimum residuals are zero. Are you aware that you are calling the Fortran77 routine and, therefore, the input matrices are expected to be in column major order?

The matrix, as entered in your program, is

11 3 -3

9 -3 0

-6 0 1

Is this the matrix you wanted, or its transpose?

I see that in your C++ caller

mecej4 — Thu, 17 Jan 2013 21:20:23 GMT

I see that in your C++ caller you set lwork = -1. That is the convention for asking Lapack to tell you how much workspace is required. Such a call does not do anything else. You have to take the value returned in lwork, allocate that much real workspace, and call Lapack again (with lwork unchanged) to do the actual calculation of the pseudoinverse.

I see that in your C++ caller

mecej4 — Thu, 17 Jan 2013 21:21:00 GMT

Duplicate Message: Please Delete!

Yeah, it works. I know this

Hai — Thu, 17 Jan 2013 22:30:00 GMT

Yeah, it works. I know this trick now.

thanks mecej4!

mecej4 wrote:

I see that in your C++ caller you set lwork = -1. That is the convention for asking Lapack to tell you how much workspace is required. Such a call does not do anything else. You have to take the value returned in lwork, allocate that much real workspace, and call Lapack again (with lwork unchanged) to do the actual calculation of the pseudoinverse.