Community
cancel
Showing results for
Did you mean: Beginner
74 Views

## QR decomposition for tall matrices

Sir,

i use the routines *qwgrf and the *orgqr to generate the Q and R .

However if the matrix is tall ie more rows than coloumns these routines just dont work (my mistake ) also while using *orgqr

it requires the mtrix to passed as an arguement. given that for tall matrices Q will be square which matrix should be passed as an arguement.

The routine works just fine for square and wide but for tall i am not able to fiqure it out.

in the code attached the matrix size is 36 elemnts, 6*6 square, 9*4 tall, 4*9 wide.

arguents passed to the executable can be either of the three eg. ./a.out 9 4

1 Solution Black Belt
74 Views

There is no Lapack-E routine containing "qwgrf" in its name, and I am unable to decipher your second sentence.

Here is an example program that performs the QR decomposition of the 6 X 4 matrix in your example:

```#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>

int main(){
float A[]={-0.57, -1.28, -0.39,  0.25,
-1.93,  1.08, -0.31, -2.14,
2.30,  0.24,  0.40, -0.35,
-1.93,  0.64, -0.66,  0.08,
0.15,  0.30,  0.15, -2.13,
-0.02,  1.03, -1.43,  0.50};
float *tau;
int iret,m,n,stride,i,j;

m=6; n=4; stride=n; tau=(float *)malloc(min(m,n)*sizeof(float));

iret=LAPACKE_sgeqrf(LAPACK_ROW_MAJOR,m,n,A,stride,tau);
printf("SGEQRF return code = %d\n",iret);
iret=LAPACKE_sorgqr(LAPACK_ROW_MAJOR,m,n,n,A,stride,tau);
printf("SORGQR return code = %d\n",iret);
for(i=0; i<m; i++){
printf("\n"); for(j=0; j<n; j++)printf(" %8.4f",A[i*stride+j]);
}
}
```

The results are in agreement with those given by Matlab.

3 Replies Black Belt
75 Views

There is no Lapack-E routine containing "qwgrf" in its name, and I am unable to decipher your second sentence.

Here is an example program that performs the QR decomposition of the 6 X 4 matrix in your example:

```#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>

int main(){
float A[]={-0.57, -1.28, -0.39,  0.25,
-1.93,  1.08, -0.31, -2.14,
2.30,  0.24,  0.40, -0.35,
-1.93,  0.64, -0.66,  0.08,
0.15,  0.30,  0.15, -2.13,
-0.02,  1.03, -1.43,  0.50};
float *tau;
int iret,m,n,stride,i,j;

m=6; n=4; stride=n; tau=(float *)malloc(min(m,n)*sizeof(float));

iret=LAPACKE_sgeqrf(LAPACK_ROW_MAJOR,m,n,A,stride,tau);
printf("SGEQRF return code = %d\n",iret);
iret=LAPACKE_sorgqr(LAPACK_ROW_MAJOR,m,n,n,A,stride,tau);
printf("SORGQR return code = %d\n",iret);
for(i=0; i<m; i++){
printf("\n"); for(j=0; j<n; j++)printf(" %8.4f",A[i*stride+j]);
}
}
```

The results are in agreement with those given by Matlab. Employee
74 Views

Hi,

As mecej4 point out, the stride in your code is wrong.  It is stride of row (row-major). So it is N in your cases.

And the sgeqrf and sorgqr can work square, tall and wide.

for example, for M=9, N=4, the leading 4 column of Q are same as Matlab.

if (M>N)

K=N;

else
K=M;

tau = (float *)malloc(K*sizeof(float));

printf("yo2\n");

LAPACKE_sgeqrf( LAPACK_ROW_MAJOR, M, N, q_mat, N, tau );

LAPACKE_sorgqr( LAPACK_ROW_MAJOR, M, K, K, q_mat, N, tau );

printf("\nQ\n");
for(i=0; i<M; i++){
printf("\n"); for(j=0; j<N; j++)printf(" %8.4f",q_mat[i*N+j]);
}

If you'd like to get the whole of Q when M=9 and N=4.

Then you may need to a temp space A  (M, M).

A=(float *)malloc(M*M*sizeof(float));

LAPACKE_sgeqrf( LAPACK_ROW_MAJOR, M, N, q_mat, N, tau );

for(i=0; i<M; i++){   for(j=0; j<N; j++)  A[i*M+j]=q_mat[i*N+j];   }

LAPACKE_sorgqr( LAPACK_ROW_MAJOR, M, M, K, A, M, tau );

for(i=0; i<M; i++){    printf("\n"); for(j=0; j<M; j++) printf(" %6.3f",A[i*M+j]);
}

Best Regards,
Ying

yo2

6.8000  -6.0500  -0.4500   8.3200
-9.6700  -2.1100  -3.3000   2.5800
2.7100  -5.1400   5.6600   5.3600
-2.7000   4.3500  -7.2600   5.9700
-4.4400   0.2700  -7.1700   6.0800
8.2300   1.0800   9.0400   2.1400
-6.8700   5.6600   5.3600  -2.7000
4.3500  -7.2600   5.9700  -4.4400
0.2700  -7.1700   6.0800   4.3500

Q

-0.3875  -0.2654   0.3130  -0.3702
0.5510  -0.4131  -0.0310  -0.1207
-0.1544  -0.3071  -0.2144  -0.3834
0.1538   0.2491   0.3303  -0.4066
0.2530  -0.0983   0.3349  -0.3051
-0.4689   0.2987  -0.3680  -0.3813
0.3914   0.2347  -0.6283  -0.1845
-0.2479  -0.4199  -0.1563   0.4326
-0.0154  -0.5218  -0.2815  -0.2765

Case 2:

Q

-0.387 -0.265  0.313 -0.370 -0.139 -0.105  0.679 -0.154 -0.171
0.551 -0.413 -0.031 -0.121 -0.308  0.509 -0.001 -0.095 -0.384
-0.154 -0.307 -0.214 -0.383  0.115 -0.470 -0.488 -0.017 -0.464
0.154  0.249  0.330 -0.407 -0.426 -0.114 -0.125  0.650  0.085
0.253 -0.098  0.335 -0.305  0.810  0.168  0.034  0.175  0.068
-0.469  0.299 -0.368 -0.381  0.055  0.621 -0.084  0.085 -0.086
0.391  0.235 -0.628 -0.185  0.106 -0.275  0.499  0.148 -0.069
-0.248 -0.420 -0.156  0.433  0.094  0.082  0.150  0.696 -0.167
-0.015 -0.522 -0.281 -0.277 -0.106  0.007 -0.067 -0.005  0.746Press any key to
continue . . . Beginner
74 Views

Sir,

Thank You all for your prompt help.

:) 