Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- dgelss does not work

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Hai

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

01-17-2013
08:48 AM

135 Views

dgelss does not work

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;

}

Link Copied

5 Replies

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

01-17-2013
11:38 AM

135 Views

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?

Hai

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

01-17-2013
12:39 PM

135 Views

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?

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

01-17-2013
01:20 PM

135 Views

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

01-17-2013
01:21 PM

135 Views

Duplicate Message: Please Delete!

Hai

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

01-17-2013
02:30 PM

135 Views

Yeah, it works. I know this trick now.

thanks mecej4!

mecej4 wrote:

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.