need help with arguments of dormqr_

[Mat]

Hello,
i realized that my arguments within my dormqr_ function cause some memory problems...
i don't know what the number of elementary reflectors is?
In dormqr.f is:

Code: Select all

 K       (input) INTEGER
*          The number of elementary reflectors whose product defines
*          the matrix Q.
*          If SIDE = 'L', M >= K >= 0;
*          if SIDE = 'R', N >= K >= 0.
*


i have L but how great should k be? 0? M? or something between?
It would be very helpful because tau has to be as great as k! If i have a tau computed from previous calculation with several values - do i have to resize tau then to the size of k? But how great is k?

for example if i have a matrix like this:
A =
1 3
0 -2
9 1

is it right that k = 2 ? is k in every case the number of columns if side = L ?

but my tau computed by dgeqrf is still as great as k is??

Thanks a lot for help
Greets

[sven]

Dear Greets,

If you used DGEQRF to perform a QR factorization of an m by n matrix A, then k should be set to min(m,n). You do not need to resize tau.

Best wishes,

Sven Hammarling.

[Mat]

thanks
the code produces correct results then BUT valgrind complies with this - I really tested it and it is the dormqr_ function!

Code: Select all

p_col =
1
0
0
1
==19442== Use of uninitialised value of size 8
==19442==    at 0x1C2B66EB: ATL_dJIK0x0x0TN0x0x0_aX_bX (in /usr/lib/atlas/libbla s.so.3.0)
==19442==
==19442== Use of uninitialised value of size 8
==19442==    at 0x1C2D9838: (within /usr/lib/atlas/libblas.so.3.0)
p_col =
-0.99449
-0.104828
0
1



if i set k to 0 (which only produces wrong results) these errors don't appear

EDIT:
here is the exampel...you can easily compile and look:

Code: Select all

#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <math.h>
#include <iomanip>
#include <vector>


extern "C"
{
    //! The QR decomposition of the input matrix A
    /*!
        Computes a QR factorization of a real M-by-N matrix A:
        A = Q * R
        \see http://www.netlib.org/lapack/double/dgeqrf.f

        \param M  The number of rows of the matrix A.  M >= 0.
        \param N  The number of columns of the matrix A.  N >= 0.
        \param A  The M-by-N matrix A.
        \param LDA The leading dimension of the array A.  LDA >= max(1,M).
        \param TAU The scalar factors of the elementary reflectors
        \param WORK  Workspace/Output
        \param LWORK The dimension of the array WORK.
        \param INFO  Function return value
    */
    void    dgeqrf_(int*       M,
                    int*       N,
                    double     A[],
                    int*       LDA,
                    double     TAU[],
                    double     WORK[],
                    int*       LWORK,
                    int*       INFO);
   
   
    //! Overwrites the matrix A with Q and C where Q is matrix defined as the product of k elementary reflectors
    /*!
        Computing:  Q^T * e_k_hat with SIDE = L and TRANS = T
        \see http://www.netlib.org/lapack/double/dormqr.f
       
        \param SIDE 'L': apply Q or Q**T from the Left; 'R': apply Q or Q**T from the Right.
        \param TRANS 'N':  No transpose, apply Q; 'T':  Transpose, apply Q**T.
        \param M The number of rows of the matrix C. M >= 0.
        \param N The number of columns of the matrix C. N >= 0.
        \param K The number of elementary reflectors.
        \param A The input matrix
        \param LDA The leading dimension of the array A.
        \param TAU TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGEQRF.
        \param C On entry, the M-by-N matrix C. On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
        \param LDC The leading dimension of the array C. LDC >= max(1,M).
        \param WORK Workspace/Output
        \param LWORK The dimension of the array WORK.
        \param INFO  Function return value
     */
    void    dormqr_(char*      SIDE,
                    char*      TRANS,
                    int*       M,
                    int*       N,
                    int*       K,
                    double     A[],
                    int*       LDA,
                    double     TAU[],
                    double     C[],
                    int*       LDC,
                    double     WORK[],
                    int*       LWORK,
                    int*       INFO);
   
   
    //! Solves a triangular system A * X = B
    /*!
        Solves a triangular system of the form
        A * X = B  or  A**T * X = B
        \see http://www.netlib.org/lapack/double/dtrtrs.f
     
        \param UPLO 'U':  A is upper triangular;'L':  A is lower triangular.
        \param TRANS  Specifies the form of the system of equations:
        \param DIAG 'N':  A is non-unit triangular; 'U':  A is unit triangular.
        \param N The order of the matrix A.  N >= 0.
        \param NRHS The number of columns of the matrix B.  NRHS >= 0.
        \param A The input matrix
        \param LDA The leading dimension of the array A.  LDA >= max(1,N).
        \param B The right hand side matrix B.
        \param LDB The leading dimension of the array B.  LDB >= max(1,N).
        \param INFO Function return value
     */
    void    dtrtrs_(char*      UPLO,
                    char*      TRANS,
                    char*      DIAG,
                    int*       N,
                    int*       NHRS,
                    double     A[],
                    int*       LDA,
                    double     B[],
                    int*       LDB,
                    int*       INFO);
}

int main(int argc, char *argv[]){
   int                             m;
    int n;
    int lda;
    int k;
    int info;
    int lwork;
    int one = 1;
                                   
    char                            *SIDE="L";
    char *TRANS="T"; //from left and transpose
    char *UPLO="U";
    char *NCHAR="N"; // upper triangular matrix not transposed
       
    std::vector< double >            work;
    std::vector< double >            tau;
    std::vector< double >            A_hat;
    std::vector< double >            e_k_hat;
    std::vector< double >            p_col;
                         
    std::vector< double >::iterator  iter;
 
 
       
    p_col.clear();
         //Create the identity vector.
    e_k_hat.push_back(1.0);
    e_k_hat.push_back(0.0);
    e_k_hat.push_back(0.0);

    A_hat.push_back(2.0);
    A_hat.push_back(0.0);
    A_hat.push_back(6.0);
    A_hat.push_back(0.0);
    A_hat.push_back(3.0);
    A_hat.push_back(1.0);
         
    p_col.push_back(1.0);
    p_col.push_back(0.0);
    p_col.push_back(0.0);
    p_col.push_back(1.0);
 
         
    m = 3;
    n = 2;
    tau.clear();
    work.clear();
    k = std::min(m,n);
    lda = std::max(m,1);
    lwork = std::max(n,1);           
    tau.resize(k);   
    work.resize(std::max(1,lwork));
 
 
         
         //QR decomposition -> Computing :   A_hat = QR
    dgeqrf_                     (&m,
                                  &n,
                                  &A_hat[0],
                                  &lda,
                                  &tau[0],
                                  &work[0],
                                  &lwork,
                                  &info);   
 
    k = std::min(m,n);
    std::cout << "p_col = " << std::endl;
    std::vector< double >::iterator it;
    for(it = p_col.begin(); it != p_col.end(); it++)
        std::cout << *it << std::endl;
         //Computing :   Q^T * e_k_hat
    dormqr_                     (SIDE,
                                 TRANS,
                                 &m,
                                 &n,
                                 &k,
                                 &A_hat[0],
                                 &lda,
                                 &tau[0],
                                 &e_k_hat[0],
                                 &lda,
                                 &work[0],
                                 &lwork,
                                 &info);
    std::cout << "p_col = " << std::endl;
    for(it = p_col.begin(); it != p_col.end(); it++)
        std::cout << *it << std::endl;
   
    return 0;
}



[Mat]

I have the bug.
In dormqr i have a wrong paramter for N! If it is greater than the number of columns of C than the function will write anywhere ....
i'm wondering why there is no exception call if the number of columns is greater than the real number of columns in the vector....? This is the reason why i had to spent over 60 hours of debugging work with 3 times reimplementation of the same code....very very very annoying!

[Julien Langou]

Hello Mat,

Sorry but this is the way it is.

If you are using LAPACK or BLAS and you provide wrong matrix-size, we have no way to figure this out at the LAPACK or BLAS level. So this is your reponsability to provide correct input parameters.

If this is too much of a burden for you then you need to have a higher level interface to LAPACK. For example Matlab. Then you will have nice error messages. Other higher level interfaces to LAPACK based on C or C++ are available around as well I think.

At our level there is nothing we can do except assuming you gave correct input parameters.

Julien