Why does LAPACKE_dormqr require double length of input array b?

[Wu__Wentao]

Using QR factorization to solve a linear least square problem requires three functions:
LAPACKE_dgeqrf for QR factorization matrix A = Q*R,
LAPACKE_dormqr for calculation of c = Q^T*b, and
dtrsm for solution of R*X = c.

The declaration of the second function, LAPACKE_dormqu, is

lapack_int LAPACKE_dormqr( int matrix_order, char side, char trans, lapack_int m, lapack_int n, lapack_int k, const double* a, lapack_int lda, const double* tau, double* b, lapack_int ldb); .

The question is array b must hold double elements than b should. Here is an example:

//Least_square_test.cpp

#include

#include

#include

#include

using namespace std;

void print(double*, int, int = 12);

int main(){

int matrix_order = LAPACK_COL_MAJOR;

lapack_int info, m(3), n(2), lda(3);

double a[6] = {0,1,2,1,1,1};

double tau[2] = {0.0,0.0};

info = LAPACKE_dgeqrf(matrix_order, m, n, a, lda, tau)); // It works well here.

lapack_int k(2), ldb(3);

const int dim_b = 6; //Note that it must double length of b!

double b[dim_b] = {0,1.1,2};

if(LAPACKE_dormqr(matrix_order, 'L', 'T', m, n, k, a, lda, tau, b, ldb)) // This will result right answer here!

cerr<<"Something wrong in Q'*B multiplication!"<

else

cout<<"Succeed in Q'*B multiplication!"<

cout<<"a = "<

print(a,6);

cout<<"tau = "<

print(tau,2);

cout<<"b = "<

print(b,dimb);

cout<

return 0;

}

void print(double *a, int N, int width){

for(int i = 0; i != N; ++i)

cout<<

cout<

}

I compile it with Intel icpc (ICC) 12.1.3 20120212 with least MKL in Ubuntu 10.04 platform, and it display the right answer b = -2.28079 0.0365148 -0.0816497.

But if I set dim_b = 3, the answer is incorrect: -2.28079 0.0879783 -0.0323584, and const array tau is also changed.

Any comments are welcome and thank you for your answers!

David Wu

[TimP]

You could study the open source and see how the arrays are used.

[Ying_H_Intel]

Hi David,

I haven't tried the code with exact icc 12.1.3 in Ubuntu 10.04 platform. But one of parameter of LAPACKE_dormqr(matrix_order, 'L', 'T', m, n, k, a, lda, tau, b, ldb)) seems be given wrong value. So you have to setting b=m*n=6.

Here m, n is not related to A, but the B in QT B ( or C in MKL manual:

m

INTEGER. The number of rows in the matrix C (m 0).

n

INTEGER. The number of columns in C (n 0).


In your test case, the B is 1 column. So if you change

const int dim_b =3; //Note that it must double length of b!
double b[dim_b] = {0,1.1,2};
n=1;

You will get correct value: 2.28079 0.0365148 -0.0816497.

Best Regards,
Ying