Thread Subject: Confused with QR decomposition and LAPACK

Hey!

I am trying to implement the qr() function (of Matlab) in C. In other
words, I am trying to make use of the LAPACK functions directly.

Matlab calls the "SGEQRF function" when we do: [Q,R] = qr(A)
So I'm trying to make use of the corresponding C-file (http://www.netlib.org/clapack/double/dgeqrf.c)
n my code. However, the results from C and Matlab don't match. I
believe I am not using the LAPACK function correctly. Could one of
you help me?

**** Here's my Matlab Code ******

x = [1, 3, 7, 9, 2, 6, 8, 0, 9, 12];
x = x(:);
mu = [mean(x); std(x)];
x = (x - mu(1))/mu(2);

n = 1;
V(:,n+1) = ones(length(x),1,class(x));
for j = n:-1:1
   V(:,j) = x.*V(:,j+1);
end

[Q,R] = qr(V);

**** Here's the C Code ****

int p;
/* 10 elements */
double x[] = {1, 3, 7, 9, 2, 6, 8, 0, 9, 12};
double mean, std;

double sum = 0.0;
for(p=0; p<10; p++){
  sum += x[p];
}
mean = sum/10;

double var = 0.0;
for(p=0; p<10; p++){
  var += powf((x[p]-mean),2);
}
std = sqrt(var/(10-1));

for(p=0; p<10; p++){
  x[p] = (x[p] - mean)/std;
}

int n = 1;
double V[10][n+1];
int i,j;
for(i=n; i>=0; i--) {
   if(i!=n) {
for(j=0; j<=10; j++) {
V[j][i] = x[j]*V[j][i+1];
}
   } else {
for(j=0; j<=10; j++) {
V[j][i] = 1;
}
   }
}

/* Converting 2D array to 1D array */
double A[10*(n+1)];
p = 0;
for(i=0; i<n1+1; i++) {
   for(j=0; j<10; j++) {
A[p++] = V[j][i];
   }
}

long M = 10, N = n1+1, LDA = MAX(1,M), lwork = N, info;
double tau[MIN(M,N)], work[lwork];
long info1 = dgeqrf_(&M, &N, A, &LDA, tau, work, &lwork, &info);

********************

Can you point out my error?

Thanks,
Nevine

Would someone be able to help me with the problem?

Thanks,
Nevine

Nevine Jacob skrev:
> Hey!
>
> I am trying to implement the qr() function (of Matlab) in C. In other
> words, I am trying to make use of the LAPACK functions directly.
>
> Matlab calls the "SGEQRF function" when we do: [Q,R] = qr(A)
> So I'm trying to make use of the corresponding C-file (http://www.netlib.org/clapack/double/dgeqrf.c)
> n my code. However, the results from C and Matlab don't match. I
> believe I am not using the LAPACK function correctly. Could one of
> you help me?

[code snipped]

> Can you point out my error?

People generally don't answer questions that are related to non-matlab
problems here. It is reasonable to assume your matlab results are
correct, and that the problem is with your lapack code. Not quite
on topic for comp.soft-sys.matlab, in other words.

Second, you don't give any indications about what matlab produces,
what your code produces, how different they are, nor exactly what
makes you think there is an error.

You COULD have posted the output of both matlab and your program
for, say, a 3x3 matrix, for people to see.

BEFORE you do that, though, check the following:

- Are the coefficients in the Q matrices equal, except for signs?
- How large are the absolute values of the differences?
- Do you call the correct version (double/float,real/complex) variant
  of LAPACK's QR routines?

You may get more response if you post in a newsgroup like
sci.maths.numerical.

Rune

> People generally don't answer questions that are related to
> non-matlab
> problems here. It is reasonable to assume your matlab results are
> correct, and that the problem is with your lapack code. Not quite
> on topic for comp.soft-sys.matlab, in other words.
>
> Second, you don't give any indications about what matlab produces,
> what your code produces, how different they are, nor exactly what
> makes you think there is an error.
>
> You COULD have posted the output of both matlab and your program
> for, say, a 3x3 matrix, for people to see.
>
> BEFORE you do that, though, check the following:
>
> - Are the coefficients in the Q matrices equal, except for signs?
> - How large are the absolute values of the differences?
> - Do you call the correct version (double/float,real/complex)
> variant
> of LAPACK's QR routines?
>
> You may get more response if you post in a newsgroup like
> sci.maths.numerical.
>
> Rune

Hey Rune,

Matlab uses the LAPACK to perform QR decomposition and I have seen a
lot of threads regarding the topic. So I believed that someone in the
group might know how to use the LAPACK functions. Using those
functions were not that direct and I am sure that the difference in
output was because of wrong usage rather than wrong function
implementation.

A few not-so-obvious pointers I got by reading the website:

1) Only pointers should be used as arguments for function.
2) A 2D array cannot be used. It needs to be converted to a
column-major 1D array.

I have incorporated both these changes but may be I am going wrong
somewhere else and I was hoping someone would point that out to me.

And by the way, the output of the LAPACK function would not give Q &
R matrices directly but can be calculated from the Output matrix. But
if I was doing things correctly, then the elements on and above the
diagonal of the Output matrix must be equal to R. But that is not the
result I am getting.

So, if you can help, please do. Also let me know if you need more
information.

Thanks,
Nevine

The code works. My bad. Will get back if I have any questions.

Thanks,
Nevine

Hi
I am having trouble recovering the q matrix from the output of LAPACK. R matrix is coming correctly, with sign change sometimes.

Could you give more pointers on how you constructed the q matrix?

Thanks

Hari

"Nevine Jacob" <nevinejacob1981@hotmail.com> wrote in message <ef42ac8.3@webcrossing.raydaftYaTP>...
> The code works. My bad. Will get back if I have any questions.
>
> Thanks,
> Nevine

You might find the Armadillo library (available at http://arma.sourceforge.net/) useful; it's basically a wrapper for LAPACK designed to simplify the syntax necessary to use LAPACK with C++. Among other things, it has an easy-to-use QR function wrapper (http://arma.sourceforge.net/docs.html#qr) and a MATLAB-to-Armadillo syntax conversion table is also provided (http://arma.sourceforge.net/docs.html#syntax), which is quite useful for converting operations originally written in MATLAB over to C++.

W Frane

"Hari Kannan" wrote in message <jgqolt$mdk$1@newscl01ah.mathworks.com>...
> Hi
> I am having trouble recovering the q matrix from the output of LAPACK. R matrix is coming correctly, with sign change sometimes.
>
> Could you give more pointers on how you constructed the q matrix?
>
> Thanks
>
> Hari
>
> "Nevine Jacob" <nevinejacob1981@hotmail.com> wrote in message <ef42ac8.3@webcrossing.raydaftYaTP>...
> > The code works. My bad. Will get back if I have any questions.
> >
> > Thanks,
> > Nevine