Problem in Big input within dormqr - need help

[Mat]

Hello,

first i want to describe what i want to solve:

I have a matrix A and want to make QR Factorization: A= QR;
Then i want to solve the Q^T * unit_vector. vec = Q^T * vec;

I have a loop around this operation in which the A matrix which is an
double array changes as well as its dimension.

Now the problem:
For small input matrix A everythink works correct even if the number of loops within for-loop is very big > 1000.
For big input matrices A everything works correkt if the number of loop iterations is small < 30. But if set the iterations to e.g.: 35 than my programms get a segmentation fault within a blas routine. If I check this with valgrind i get several times the following output :

Code: Select all

Use of uninitialised value of size 4
==7388==    at 0x46BE6DC: (within /usr/lib/sse2/libf77blas.so.3.0)


After valgrind returns i have about 16000 errors which cause my seg-fault i think.

I think that this could be because a parameter is not set correct within my for-loop which is neede by the blas routine. Perhaps anybody would see if i made a mistake?


I'm programming in C++ so my code looks like this:

Code: Select all

for(...)
{
       double* work, *tau;
       int one = 1, info, lwork, lda,k;
     
        //m and n are set correctly!!!

      k       = std::min(m, n);
   lda    = std::max(m, 1);
   lwork    = std::max(n, 1);
   int max = std::max(1, lwork);
   *work = new double[max];
   *tau = new double[k];
   
   memset(*work, 0, max * sizeof(double));
   memset(*tau, 0, k * sizeof(double));

   dgeqrf_(&m,
         &n,
         A,
         &m,
         tau,
         &work,
         &lwork,
         &info);

        char * SIDE = "L";
        char * TRANS = "T";
   dormqr_(SIDE,
         TRANS,
         &m,
         &one,
         &k,   
         A,
         &lda,
         tau,
         vec,
         &lda,
         work,
         &lwork,
         &info);

     delete [] work;
     delete [] tau;
}


The strange thing is, if is reset the tau array between dgeqrf and dormqr with the following line everything works correct:

memset(tau, 0, k * sizeof(double));

But this way the output is incorrect of course.

Thus i suspect anything is incorrect with the tau array. Can anybode help me with that or say what exactly the tau array is?

For help i'm really thankful.

[Mat]

Ok, now i transfered my problem into an exterior mini-programm without the loop solving only once my mathematical problem.

I have exactly the input which causes the error, but i don't know why.

The problem is, the input is quite large, the Matrix A has dimension 66*135.

I would like to know whether the tau array can be given to dormqr without doing anything with it as it is returned from dgeqrf?

As is can see all input-parameters are correct.

This is very strange...

Perhaps anybode has an idea what i could test next?

[Mat]

:cry:
I cannot find the bug.

Now i have the following input matrix:

Code: Select all

double* A = new double[135 * 66];


for(int i = 0; i < 66; i++)
{
   for (int j = 0; j < 135; j++)
      A[i*135 + j] = 10.0 + i * 2.3 - j* 5.4;
}

It's just a full filled matrix.
Now it is obvious that something with the lapack parameters is incorrect. But what?

Even my elementary reflectors (tau) looks correct:

Code: Select all

1.0021  1.16911  1.03951  1.05207  1.0095  1.01703  1.02418  1.03819  1.01235  1.01501  1.00746  1.00521  1.01367  1.01908  1.04091  1.03324  1.0743  1.02194  1.02594  1.03217  1.01987  1.00419  1.04184  1.02427  1.0284  1.01624  1.01455  1.01376  1.03884  1.08709  1.13553  1.00554  1.00855  1.12581  1.01849  1.00019  1.00075  1.00776  1.04228  1.03543  1.12985  1.07531  1.02637  1.13092  1.04228  1.08012  1.04777  1.08982  1.07744  1.01998  1.10473  1.11752  1.08096  1.05638  1.09937  1.05999  1.0009  1.05829  1.06111  1.19898  1.20038  1.03701  1.09602  1.04784  1.2104  1.03766

The other parameters are:
m: 135
n: 66
k: 66
lda: 135
lwork: 66

[Mat]

Ok,
i coded now the mathematical operation with dlarf this way and it works correct.

Code: Select all

   char * side = "L";
   double* work_new = new double[n];
   double diag_safe;
   memset(work_new, 0, n * sizeof(double));
   double tau_i;
   int cols_dlarf;

   for(int i = 0; i < n; i++)
   {
      tau_i = tau[i];
      cols_dlarf = m - i;
      diag_safe = *A;
      *A  = 1;
      dlarf_(side, &cols_dlarf, &one, A, &one, &tau_i, vec, &m, work_new );
      *A = diag_safe;
      vec++;
      A += (m + 1);
   }


Can anybode tell me why my method with using dormqr doesn't work correctly? I would prefer avoiding this for-loop and using dormqr. But why the hell it doesn't work with the above code for huge problems?

[Mat]

ok - now i need really your help.

I now have reduced the problem to the minimum which still cuases 2 errors although i think my program works correctly.

Perhaps anybode could just test it within a main-method.

Code: Select all

double * A = new double[2*9];

A[0] = -16809.7
A[1] =  6.66667;
A[2] =   0;
A[3] =   160;
A[4] =   0;
A[5] =   6250;
A[6] =   0;
A[7] =   25.6;
A[8] =   3.33333;
A[9] =   3.33333;
A[10] =   -16809.7;
A[11] =   6.66667;
A[12] =   0;
A[13] =   160;
A[14] =   0;
A[15] =   6250;
A[16] =   25.6;
A[17] =   0;
 
        int m = 9;
        int n = 2;

        int info = 0;
   int k       = std::min(m, n);
   int lda    = std::max(m, 1);
   int lwork    = std::max(n, 1);
   int max = std::max(1, lwork);
   double *work = new double[max];
   double *tau = new double[k];
   
   memset(work, 0, max * sizeof(double));
   memset(tau, 0, k * sizeof(double));

dgeqrf_(&m,
         &n,
         A,
         &lda,
         tau,
         work,
         &lwork,
         &info);

        double* vec = new double[m];
        memset(vec, 0, m * sizeof(double));
        vec[2] = 1.0;

        char * side = "L";
        char * TR    = "T";
        int one = 1;
        dormqr_(side,
                        TR,
                        &m,
                        &one,
                        &k,
                        A,
                        &lda,
                        tau,
                        vec,
                        &lda,
                        work,
                        &lwork,
                        &info);


This causes 2 errors, but i don't know why - it would be very nice if anybode could help!

The output of dgeqrf_ produces tau as following:
tau =
1.93727 1.93727

[Julien Langou]

The C code below seems to work fine. (Valgrind do not complain at least.)
What is your problem exactly?
(NB if you only have one vector on which you want to apply the Householder
reflection, then it makes sense to call DLARF directly.)

Code: Select all

int main(){

    double A[2*9];

    int m = 9;
    int n = 2;

    int info = 0;
    int k = n;          /* k = min(m,n);       */
    int lda = m;        /* lda = max(m,1);     */
    int lwork = n;      /* lwork = max(n,1);   */
    int max = lwork;    /* max = max(lwork,1); */

    double *work;
    double *tau;

    char *side = "L";
    char *TR    = "T";
    int one = 1;

    double *vec;

    A[0] = -16809.7;
    A[1] =  6.66667;
    A[2] =   0;
    A[3] =   160;
    A[4] =   0;
    A[5] =   6250;
    A[6] =   0;
    A[7] =   25.6;
    A[8] =   3.33333;
    A[9] =   3.33333;
    A[10] =   -16809.7;
    A[11] =   6.66667;
    A[12] =   0;
    A[13] =   160;
    A[14] =   0;
    A[15] =   6250;
    A[16] =   25.6;
    A[17] =   0;

    work = (double *) malloc( max * sizeof( double ) );
    tau  = (double *) malloc( k * sizeof( double ) );
    vec  = (double *) malloc( m * sizeof( double ) );


    memset(work, 0, max * sizeof(double));
    memset(tau, 0, k * sizeof(double));

    dgeqrf_(&m, &n, A, &lda, tau, work, &lwork, &info);
    printf("tau[0] = %f tau[1] = %f\n",tau[0],tau[1]);

    memset(vec, 0, m * sizeof(double));
    vec[2] = 1.0;

    dormqr_(side, TR, &m, &one, &k, A, &lda, tau, vec, &lda, work, &lwork, &info);

    free(vec);
    free(tau);
    free(work);
}


[Mat]

Sorry, i posted the wrong code which causes the error:

The following code produces this error:

Code: Select all

nvalid read of size 8
at 0x80489DB: main (in /home/mat/Desktop/test/test_exec)
Address 0x599F758 is 0 bytes after a block of size 8 alloc'd
at 0x4021620: malloc (vg_replace_malloc.c:149)
by 0x804896B: main (in /home/mat/Desktop/test/test_exec)
tau[0] = 1.000512 tau[1] = 0.000000



Code: Select all

int main(){

    double A[1*9];

    int m = 9;
    int n = 1;

    int info = 0;
    int k = n;          /* k = min(m,n);       */
    int lda = m;        /* lda = max(m,1);     */
    int lwork = n;      /* lwork = max(n,1);   */
    int max = lwork;    /* max = max(lwork,1); */

    double *work;
    double *tau;

    char *side = "L";
    char *TR    = "T";
    int one = 1;

    double *vec;

    A[0] = 36.5714;
    A[1] = 36.5714;
    A[2] = 36.5714;
    A[3] = 36.5714;
    A[4] = 146.286;
    A[5] = 3.33333;
    A[6] = -66901.5;
    A[7] = 6.66667;
    A[8] = 25000;



    work = (double *) malloc( max * sizeof( double ) );
    tau  = (double *) malloc( k * sizeof( double ) );
    vec  = (double *) malloc( m * sizeof( double ) );


    memset(work, 0, max * sizeof(double));
    memset(tau, 0, k * sizeof(double));

    dgeqrf_(&m, &n, A, &lda, tau, work, &lwork, &info);
    printf("tau[0] = %f tau[1] = %f\n",tau[0],tau[1]);

    memset(vec, 0, m * sizeof(double));
    vec[2] = 1.0;

    dormqr_(side, TR, &m, &one, &k, A, &lda, tau, vec, &lda, work, &lwork, &info);

    free(vec);
    free(tau);
    free(work);
}



in my main program this code produces the following error:

Code: Select all

Conditional jump or move depends on uninitialised value(s)
==12999==    at 0x473427D: ATL_daxpy (in /usr/lib/sse2/libf77blas.so.3.0)


[Julien Langou]

Remove the print of tau[1]. tau is of size 1 in this case.

[Mat]

Thanks,

i can reproduce the previous examples with no errors.

But what i wanted is to reduce my matrix A because it is very big. Now i have a file with all values (8910) of the matrix. And the single file produces the same error as my program.

It would be very kind if anybody could test the following code once more.
The problem is, that i would attach the file with the entries but there is no option here.

I attached the matrix file to another forum under:
http://www.prad.de/board/thread.php?pos ... post281633
at the bottom there is an attachement named test.txt
You can download it.


The code is:

Code: Select all

#include <iostream>
#include "Lapack.h"

int main(){

    int m = 135;
    int n = 66;

    double * A = new double[n*m];


    int info = 0;
    int k = n;          /* k = min(m,n);       */
    int lda = m;        /* lda = max(m,1);     */
    int lwork = n;      /* lwork = max(n,1);   */
    int max = lwork;    /* max = max(lwork,1); */

    double *work;
    double *tau;

     char *side = "L";
     char *TR    = "T";
     int one = 1;

    double *vec;

   float val;
   char    line[128];
        //NOTICE: ENTER PATH to file
   FILE *f=fopen("/path_to_file","r");
   for(int i = 0; i< 8910 ;i++)
   {
      fgets(line,128,f);

      
      sscanf(line, "%g\n", &val);
      A[i] = val;
      std::cout << val << std::endl;
   }
   fclose(f);


    work = new double[max];
    tau = new double[k];
    vec = new double[m];

       
    memset(work, 0, max * sizeof(double));
    memset(tau, 0, k * sizeof(double));

    dgeqrf_(&m, &n, A, &lda, tau, work, &lwork, &info);

    memset(vec, 0, m * sizeof(double));
    vec[2] = 1.0;

    dormqr_(side, TR, &m, &one, &k, A, &lda, tau, vec, &lda, work, &lwork, &info);

   
   //Free heap storage
   delete [] A,
   delete [] vec;
   delete [] work;
   delete [] tau;
}


[Julien Langou]

The following code works fine for me (under valgrind as well) ...
It is your code, I have just removed the C++ syntax.
Can you try it?

Code: Select all

#include <stdio.h>
int main(){


    int m = 135;
    int n = 66;
    double A[m*n];

    int info = 0;
    int k = n;          /* k = min(m,n);       */
    int lda = m;        /* lda = max(m,1);     */
    int lwork = n;      /* lwork = max(n,1);   */
    int max = lwork;    /* max = max(lwork,1); */

    double *work;
    double *tau;

    char *side = "L";
    char *TR    = "T";
    int one = 1;
    int i;

    double *vec;

    float val;
    char line[128];
    FILE *f;

    f =fopen("text.txt","r");
    for( i = 0; i< 8910 ;i++)
    {
        fgets(line,128,f);
        sscanf(line, "%g\n", &val);
        A[i] = val;
    }
    fclose(f);

    work = (double *) malloc( max * sizeof( double ) );
    tau  = (double *) malloc( k * sizeof( double ) );
    vec  = (double *) malloc( m * sizeof( double ) );

    memset(work, 0, max * sizeof(double));
    memset(tau, 0, k * sizeof(double));

    dgeqrf_(&m, &n, A, &lda, tau, work, &lwork, &info);
    //printf("tau[0] = %f tau[1] = %f\n",tau[0],tau[1]);

    memset(vec, 0, m * sizeof(double));
    vec[2] = 1.0;

    dormqr_(side, TR, &m, &one, &k, A, &lda, tau, vec, &lda, work, &lwork, &info);

    free(vec);
    free(tau);
    free(work);
}


[Mat]

Thanks a lot, really - your help is unbelievale.

But i get the following message from valgrind as in all my other programms:

Code: Select all

Use of uninitialised value of size 4
==16500==    at 0x46BA2EC: (within /usr/lib/sse2/libf77blas.so.3.0)
==16500==
==16500== Use of uninitialised value of size 4
==16500==    at 0x46BE6DC: (within /usr/lib/sse2/libf77blas.so.3.0)
==16500==
==16500== Use of uninitialised value of size 4
==16500==    at 0x46BE31C: (within /usr/lib/sse2/libf77blas.so.3.0)
==16500==
==16500== ERROR SUMMARY: 139 errors from 3 contexts (suppressed: 27 from 1)
==16500== malloc/free: in use at exit: 1,902 bytes in 1 blocks.
==16500== malloc/free: 134 allocs, 133 frees, 41,814 bytes allocated.
==16500== For counts of detected errors, rerun with: -v
==16500== searching for pointers to 1 not-freed blocks.
==16500== checked 1,310,872 bytes.


I'm working under Kubuntu 7.04.
My lapack version is: 3.020000531a-6ubuntu3
My blas version is: Basic Linear Algebra Subroutines 3 1.2.8-ubuntu2

What can I do now?
This bug makes me crazy...

[Julien Langou]

Hello,
Can you try to install BLAS and LAPACK by yourself?
o Download them from http://www.netlib.org/lapack/lapack-lite-3.1.1.tgz
o Then untar,
o Edit the make.inc file with your compiler (there is a make.inc.example that should almost
do it)
o Type make blaslib then make lib
After this you'll have two .a files. (lapack and blas), a good idea is to rename as
libreflapack.a and librefblas.a
Link your application with these.
I am curious to see the result.
Julien.

[Mat]

YES!
That was it!

Now i have 0 errors as well as under valgrind.

I try to stay calm and not to explode, this took about 4 days of hard work!
Thanks a lot Julien for your patient help!

I try to get back under control.
In future i will use only the "original" staff from the developers sites.

Thank again. :D

[Julien Langou]

By the way, since your problem are failry big, you should definitely try to use an optimized
BLAS library. So replace librefblas.a by another BLAS. For example, Goto BLAS.
Julien.

[bravegag]

Hi,

First off all, very nice thread! it helped me solve an issue I had while using dlarf.

Is there something like a dormqr "lite" version? This is in the context of doing QR delete column block updates. Basically I need to apply a set of elementary reflectors compacted in matrix A starting from a column K and only a block close to a small band right below the diagonal, this small band is of width P (how many columns were deleted). See the code below which is an adaptation of http://www.maths.manchester.ac.uk/~cluc ... delcolsq.f in order to apply the elementary reflectors to an arbitrary matrix. My question is whether there exist a dormqr "lite" that could compute the same?

Code: Select all

      LAST = MIN( M-1, N )
*
      DO 10 J = K, LAST
*
         LENH = MIN( P+1, M-J+1 )
*
*        Apply H(J) to trailing matrix from left
*
         AJJ = A( J, J )
         A( J, J ) = ONE
         CALL DLARF( 'L', LENH, M-J, A( J, J ), 1, TAU( J-K+1 ),
     $               Q( J, J ), LDQ, WORK )
*
         A( J, J ) = AJJ
*
   10 CONTINUE


Many thanks in advance,
Best regards,
Giovanni