dgeqrf(LAPACK) not working properly in c++

Question

I would grateful if you could help me to resolve the following issue. I have stumbled upon the wall and i can't really tell what is the problem here. I would like to use the LAPACK functions dgeqrf, dormqr and dtrtrs under c++ to solve the homogeneous system (15 times 15). But i tried to see if i can solve much simpler system for the first hand. What i am posting here is my simple usage of the dgeqrf.

#include "lapacke.h"
#include "lapacke_config.h"
#include "lapacke_mangling.h"
#include "lapacke_utils.h"
using namespace std;

int main()
{

int     INFO=3;
double WORK[3];
int LWORK=3;
double TAU[3];
int     LDA = 3;
int     LDB = 3;
int     N = 3;


double det;

double  A[9] =
{
    1, 1, 0,
    0, 1, 0,
    0, 1, 1
};

// end of declarations
LAPACK_dgeqrf(&N,&N,A,&LDA,TAU,WORK,&LWORK,&INFO);
for(int i=0;i<N*N;i++){
cout<<A[i]<<",";
}
cout<<endl;
if (INFO!=0){
return 0;
}
det=pow(-1,N-1);
 for(int i=0;i<N;i++){
det = det * A[i*N+i];
}
cout<<det<<endl;

What i do not understand is why i am getting R-matrix which is not correct? The printout gives:

 -1.41421,0.414214,0,-0.707107,0.707107,0,-0.707107,0.707107,1,

If i check it against the other online matrix calculators, it is wrong. I understand that upper triangle together with the diagonal should be my R-matrix. And what is more surprising i obtain the determinant which is correct (i checked it against more simple examples). R-matrix should be unique, right? If i want to solve the 15 times 15 problem, i should have this first step resolved...i hope someone could explain it to me. I would like to say upfront that i tried to use all three subroutines to solve the simple 3 times 3 system but it gave meaningless results. I hope that once i got this one right, i will understand the issue with dormqr and drtrs as well.

Cheers!

score 3 · Accepted Answer · edited Feb 09 '17 at 15:10

3

This is a row-major vs. column major issue it seems. In C/C++ we usually have matrices row-major, but in FORTRAN (in which the lapack you use is probably implemented) uses a column-major format.

Try changing rows and columns in A:

double A[9] = { 1, 0, 0, 1, 1, 1, 0, 0, 1 };

Considering the matrix you posted as a column-major matrix and plugging it into an online QR-calculator (http://www.bluebit.gr/matrix-calculator/) I get the R-matrix you post, with the exception that the first row has its sign inverted, which means that the first Q-vector just points in the opposite direction.

$\mathbf{A} = \begin{pmatrix} 1.000 & 0.000 & 0.000 \\ 1.000 & 1.000 & 1.000 \\ 0.000 & 0.000 & 1.000 \end{pmatrix} \ \ \ \ \mathbf{R} = \begin{pmatrix} 1.414 & 0.707 & 0.707 \\ 0.000 & 0.707 & 0.707 \\ 0.000 & 0.000 & 1.000 \\ \end{pmatrix}$

edited Feb 09 '17 at 15:10

Community

1
1

answered Nov 13 '15 at 13:02

Banan

434
6
12

Hi Banan, thanks, this solves my problem! I got the second row inverted, when i checked against the wolfram alpha calculator. I see that the output is also column-major. Thanks again! – OwlBite Nov 13 '15 at 13:32
@OwlBite No problem, these things can be tricky. I've been working with stuff like this forever, and still get it wrong sometimes ;) – Banan Nov 13 '15 at 13:35

dgeqrf(LAPACK) not working properly in c++

1 Answers1