a confusion about QR factorization

[Powergon]

Dear fellows, I just registered in this forum b/z I found so many nice guys and helpful solutions in here.

I have a question and hope someone could give a solution.

I am looking for an orthogonal basis for {u: Au=0}, where u is in R^n, A is a matrix m*n (m<<n), with rank m.
One way is to do QR factorization for A^T (transpose of A), the last n-m columns of Q will be the orthogonal basis for the null space.
That is, if written Q=[Q1,Q2]. Q2 is what I need.

QR factorization can be done with LAPCK subroutine DGEQRF, then the orthogonal Q matrix can be found with DORMQR.

However, after done this, I tested that when apply A to one column of Q2, it's not zeros (as large as ~1.e-2). So, Q2 is not exactly the basis for the null space.

It really confuses me...could someone suggest the possible reasons?

Thanks a lot!

-Powergon

[Powergon]

Please, anyone here can help?

The 5*3 matrix I am testing
A=[9 5 3
4 5 0
8 4 5
9 6 4
5 6 7
];

I call dgeqrf first; then to get the 2*2 Q2, call dormqr, where I set SIDE="L", TRANS="N", the matrix variable
C=[0 0
0 0
0 0
1 0
0 1
];

Q2 is not even orthogonal !

Thanks a lot!

[Powergon]

Again, if i set, in the call of dormqr, C as the 5*5 identity matrix, on exit, C will be overwritten as Q*C which is Q.

And this Q is orthogonal. If i chose the last 2 columns as Q2, then Q2 is orthogonal. This is expected.

But whats wired is, if I set C as the last 2 columns of an identity matrix, the overwritten C, which is Q*C should directly equals to Q2, isn't it?

And still, A^T * Q2 is not zeros. (Q2 is expected to be the null basis)

[Julien Langou]

I am confused by your email, you start with A m-by-n with m<=n and rank m. Fine. Then you give a matrix A which 5-by-3, rank 3. Which case are you? Julien.

[Powergon]

I am confused by your email, you start with A m-by-n with m<=n and rank m. Fine. Then you give a matrix A which 5-by-3, rank 3. Which case are you? Julien.

Thank you so much, Julien, for the reply. And sorry for the confusion.
Actually the second A i gave is the transpose of the first one.
Actually I have figured out my problem: the stupid index order problem!

I have a C program, but am calling lapack subroutine in Fortan. C has row-major order while Fortran follows column-major roder..that causes the mess..

Forget about this problem, can I　ask you a new question here, Julien?

I need to find the smallest eigenvalue for a very LARGE (say, 4000*4000) symmetric positive definite real matrix.
I am using DSYEVR. it works fine for small and medium scale matrix. But for this large matrix, it generates "internal error" (returned info=4>0).The vector
storing eigenvalues was returned empty..

BTW, i am using threaded mkl.

So, what's the possible reason for this problem? I am not an expert, sincerely hope you or someone could give a kind explanation.

Thank you so much!

[Julien Langou]

Unfortunately, this looks a lot like a bug on our side ...
So can you please try to link not with Intel MKL but with reference LAPACK from netlib? And link with reference BLAS so as to be on the safe side ...
http://www.netlib.org/lapack/#_lapack_version_3_4_2
Also you can try to increase the workspaces just in case.
But when a code works on small sizes and crash on larger sizes this is never a good sign ...

( I am wondering if this could not be an integer overflow problem.)

[Powergon]

Unfortunately, this looks a lot like a bug on our side ...
So can you please try to link not with Intel MKL but with reference LAPACK from netlib? And link with reference BLAS so as to be on the safe side ...
http://www.netlib.org/lapack/#_lapack_version_3_4_2
Also you can try to increase the workspaces just in case.
But when a code works on small sizes and crash on larger sizes this is never a good sign ...

( I am wondering if this could not be an integer overflow problem.)

Thanks, Julien, I traced the error and found that the entries in my matrix (whose eigenvaules are needed) are NAN (so, it seems at least DSYEVR is innocent..).
The computing of my matrix involves the call of DORMQR, which computes Q*C, where Q is from QR factorization by DGEQRF and C is another large matrix that will be replaced by Q*C.
I tested, if I set C as identity matrix, it works well; but if I set C as another matrix, then it starts to give me NAN entries in matrix C. All before this, my matrix is "good".
So, what's the actual reason for it?

As you suggested, I tried to increase the dimension of workspace (though i was using the optimal value returned by a workspace query) and also, considering that DORMQR might not be thread-safe, i tried
sequential MKL. Still, the above problem arises.

Thank you very much!

[Powergon]

Unfortunately, this looks a lot like a bug on our side ...
So can you please try to link not with Intel MKL but with reference LAPACK from netlib? And link with reference BLAS so as to be on the safe side ...
http://www.netlib.org/lapack/#_lapack_version_3_4_2
Also you can try to increase the workspaces just in case.
But when a code works on small sizes and crash on larger sizes this is never a good sign ...

( I am wondering if this could not be an integer overflow problem.)

To simplify my question, why DORMQR(...A,....C1....) works while DORMQR(.....A.....C2...) not (gives NAN entries in C2 on exit), with all other variables the same ? What's the essential difference in C1 and C2 could cause such a problem?
Thank you very much!

[Julien Langou]

This is really weird ... I am not clear on why the reason. The entry C1 or C2 is essentially not relevant. It should always work or never work.
Can you please check for NaN in all your entries when you input DORMQR?

[Powergon]

This is really weird ... I am not clear on why the reason. The entry C1 or C2 is essentially not relevant. It should always work or never work.
Can you please check for NaN in all your entries when you input DORMQR?

Thank you very much, Julien..ye, it's my fault. I traced the entries and
figured out the problem. At present, LAPCK subroutines are innocent....

Really appreciate your patient help, Julien. And thanks a lot for your contributive work on LAPCK.!