QR decomposition for tall matrices

[hemantp_p_]

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

[mecej4 · Accepted solution]

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.

[Ying_H_Intel]

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 . . .

[hemantp_p_]

Sir,

Thank You all for your prompt help.

:)