Problem in Big input within dormqr - need help

Open discussion regarding features, bugs, issues, vendors, etc.

Problem in Big input within dormqr - need help

Postby Mat » Sat Sep 15, 2007 9:15 am

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
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Mat » Sat Sep 15, 2007 4:19 pm

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
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Mat » Sun Sep 16, 2007 2:16 pm

: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
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Mat » Sun Sep 16, 2007 3:41 pm

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
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Mat » Mon Sep 17, 2007 10:06 am

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
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julien Langou » Mon Sep 17, 2007 11:53 am

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);
}
Julien Langou
 
Posts: 735
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Postby Mat » Mon Sep 17, 2007 12:21 pm

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)
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julien Langou » Mon Sep 17, 2007 12:33 pm

Remove the print of tau[1]. tau is of size 1 in this case.
Julien Langou
 
Posts: 735
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Postby Mat » Mon Sep 17, 2007 3:56 pm

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;
}
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julien Langou » Mon Sep 17, 2007 5:17 pm

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);
}
Julien Langou
 
Posts: 735
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Postby Mat » Mon Sep 17, 2007 5:25 pm

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...
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julien Langou » Mon Sep 17, 2007 5:45 pm

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.
Julien Langou
 
Posts: 735
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Postby Mat » Mon Sep 17, 2007 6:08 pm

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
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julien Langou » Mon Sep 17, 2007 6:34 pm

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.
Julien Langou
 
Posts: 735
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Problem in Big input within dormqr - need help

Postby bravegag » Mon Apr 16, 2012 11:28 am

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
bravegag
 
Posts: 12
Joined: Sat Mar 24, 2012 3:23 pm


Return to User Discussion

Who is online

Users browsing this forum: Google [Bot] and 1 guest