lapackpp-devel

[Lapackpp-devel] bugs in operator= for LaIndex/LaGenMatComplex/LaGenMatDouble?

[Brian W. <bw...@cs...>]

Hello,

I think I've uncovered a few bugs.

In LaIndex::oeprator=(const LaIndex &i), the following definition
ignores the increment of the LaIndex it is copying:

        inline LaIndex& operator=(const LaIndex& i){
            start_=i.start_; inc_=1; end_=i.end_;
            return *this;}

I think inc_=1 should be inc_=i.inc_:

        inline LaIndex& operator=(const LaIndex& i){
            start_=i.start_; inc_=i.inc_; end_=i.end_;
            return *this;}

Also, LaGenMatComplex::operator=(COMPLEX s) and
LaGenMatDouble::operator=(double s) seem to be problematic since both
assign to the vectorDouble v without accounting for the active view
of the LaGenMatDouble object.  I believe this means that if we assign
(via operator=) through a submatrix view, the entire matrix will be
effected.

So, instead of doing

    v = s;

in both LaGenMatComplex::operator=(COMPLEX s) and
LaGenMatDouble::operator=(double s), we should instead have:

    int i,j;

    for (j=0; j<size(1); j++)
    {
        for (i=0; i<size(0); i++)
        {
            (*this)(i,j) = s;
        }
    }

and something similar for the complex case.  Note that this
mirrors operator+= (obviously, with += replaced by =).

Below is some example code (taken essentially from the Lapack++ 1.1
manual when it was still being developed by NIST/Tennessee/Oak Ridge).
I've included what I believe to be the correct output and also
the output I get before making the above changes.

By the way, this differs from the example in the manual only
in so far as assignments are replaced by references.
i.e., B = A(II,II) is replaced below by B.ref(A(II,II)).
If I'm reading that manual correctly, B = A(II,II) should
establish B as a reference to elements in A.  Instead it
appears to perform a (deep) copy.

Thanks,
Brian

int main(int argc, char **argv)
{

  LaGenMatDouble A(10,10), B, C;
  LaIndex II(1,9,2);
  LaIndex JJ(1,1,1);

  B.ref(A(II,II));   // B references the following elements of A:
                     //    (1, 1), (1, 3), (1, 5), (1, 7), (1, 9)
                     //    (3, 1), (3, 3), (3, 5), (3, 7), (3, 9)
                     //    (5, 1), (5, 3), (5, 5), (5, 7), (5, 9)
                     //    (7, 1), (7, 3), (7, 5), (7, 7), (7, 9)
                     //    (9, 1), (9, 3), (9, 5), (9, 7), (9, 9)
  B(2,3) = 3.1;      // Set A(5,7) = 3.1

  C.ref(B(LaIndex(2,2,4), JJ));
                     // C references the following elements of B:
                     //    (2, 1)
                     // C references the following elements of A:
                     //    (5, 3)

  C = 1.1;           // Set B(2,1) = A(5,3) = 1.1

  cout << A << endl;
  cout << B << endl;
  cout << C << endl;

  return 0;
}

Correct output:

0  0  0  0  0  0  0  0  0  0
7.94687e-311  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  1.1  0  0  0  3.1  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0

// Note the following matrix is 5x5
0  0  0  0  0
0  0  0  0  0
0  1.1  0  3.1  0
0  0  0  0  0
0  0  0  0  0

1.1


Incorrect output:

1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1

// The following matrix is 9x9 because the
// view (II=(1,9,2),II) is incorrectly copied to
// (I=(1,9,1),I) when the increment is set to 1
// irrespective of the LaIndex being copied.

// Also, we have assigned to the entire matrix,
// rather than simply the view specified in
//   C.ref(B(LaIndex(2,2,4), JJ));
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1
1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1  1.1

1.1

Re: [Lapackpp-devel] bugs in operator= for LaIndex/LaGenMatComplex/LaGenMatDouble?

[Christian S. <sti...@tu...>]

Dear Brian,

Brian White schrieb:
> In LaIndex::oeprator=(const LaIndex &i), the following definition
> ignores the increment of the LaIndex it is copying:
> 
>         inline LaIndex& operator=(const LaIndex& i){
>             start_=i.start_; inc_=1; end_=i.end_;
>             return *this;}
> 
> I think inc_=1 should be inc_=i.inc_:

Yes, that's a bug. Thanks for reporting this. It has been fixed in CVS now.

> Also, LaGenMatComplex::operator=(COMPLEX s) and
> LaGenMatDouble::operator=(double s) seem to be problematic since both
> assign to the vectorDouble v without accounting for the active view
> of the LaGenMatDouble object.  I believe this means that if we assign
> (via operator=) through a submatrix view, the entire matrix will be
> effected.

Yes, that's a bug, too. I've fixed this by first checking for unit 
stride and then either assigning to the vectorDouble object (*much* 
faster!) or iterating over all elements.

> Below is some example code (taken essentially from the Lapack++ 1.1
> manual when it was still being developed by NIST/Tennessee/Oak Ridge).
> I've included what I believe to be the correct output and also
> the output I get before making the above changes.

Thanks for adding this. It helped a lot to verify the correct fix.

> By the way, this differs from the example in the manual only
> in so far as assignments are replaced by references.
> i.e., B = A(II,II) is replaced below by B.ref(A(II,II)).
> If I'm reading that manual correctly, B = A(II,II) should
> establish B as a reference to elements in A.  Instead it
> appears to perform a (deep) copy.

Yes, that's true and it's a feature :-) . Seriously, I've decided 
against that original lapackpp-1.1 approach (now also more clearly 
documented in the header) because of two reasons: 
LaGenMatDouble::LaGenMatDouble(const LaGenMatDouble &) [i.e. the copy 
constructor] was working as a by-value deep copy anyway. So having 
operator= [i.e. the copy assignment] work as a by-reference copy means 
that the copy constructor and the copy assignment were doing two 
different things, which is a bad thing. Secondly, many users of lapackpp 
got quite confused when they constructed a new LaGenMatDouble, assigned 
an initial value by the copy assignment (i.e. LaGenMatDouble B = A;) and 
to their surprise discovered that the newly created matrix won't behave 
as a standalone object but rather as a reference to some other object. 
If you want, you can say that the by-reference semantics were too 
difficult to grok for the average programmer. Therefore I changed the 
copy assignment to be the same as the copy constructor: a by-value deep 
copy. If you want a by-reference copy, use ref(). Those who use it are 
probably quite sure they know what they do.

Regards,

Christian

[Lapackpp-devel] bugs in reference counting in vd.{cc,h}?

[Brian W. <bw...@cs...>]

> Yes, that's true and it's a feature :-) .

Thanks for the explanation, Christian.

I may have another bug.  I can't pinpoint exactly what the problem is, but
I can manifest it repeatedly and I can fix it. :)

I've included a program called manifestMemProb.C, which is similar to the
last program I sent out.  It creates two different views of a matrix (plus
the original) and modifies them.  It does the same procedure on a
different matrix.  Somehow, the second matrix is being corrupted, as you
can see in the file manifestMemProb.incorrect.

Like I said, I was unable to track down the exact problem.  It may be
related to new'ing and delete'ing 0 byte chunks of memory?  I added a
check for n=0 in the following constructor:

VectorDouble::VectorDouble(int n=0)

In that case, I set ref_count = 0 and p->data = data = NULL:

 VectorDouble::VectorDouble(int n=0)
{
    assert(n>=0);
    p = new vrefDouble;
    p->sz = n;
    if ( n > 0 ) {
      p->data = data = new double[n];
      p->ref_count = 1;
    } else {
      p->data = data = NULL;
      p->ref_count = 0;
    }
}

I then only decremented ref_count if it was > 0.  And asserted explicitly
that p->data was non-NULL before delete'ing it.  However, even in the
original code this assertion was never hit.  At any rate, this seems to
have fixed the manifestation of the bug.  Sorry I can't provide a better
understanding of what was going on.

I ended up changing all of the constructors to do something similar.  I
also changed the VectorDouble destructor and ref (where ref_count is
decremented) as described above.  Also, in ref I don't increment the
ref_count unless the associated data is non-NULL.  In short then:
ref_counts are only > 0 when data != NULL.  Also, I never allocate or free
0-byte data.

I imagine other vectors will have the same type of problems, but I didn't
fix those.

I'm sending my modified vd.cc and vd.h.

By the way, please let me know if I'm somehow abusing the interface in my
example such that I am causing the seemingly incorrect behavior.

Thanks,
Brian

[Lapackpp-devel] Re: bugs in reference counting in vd.{cc,h}?

[Christian S. <sti...@tu...>]

Dear Brian,

Brian White schrieb:
> I may have another bug.  I can't pinpoint exactly what the problem is, but
> I can manifest it repeatedly and I can fix it. :)
> 
> I've included a program called manifestMemProb.C, which is similar to the
> last program I sent out.  It creates two different views of a matrix (plus
> the original) and modifies them.  It does the same procedure on a
> different matrix.  Somehow, the second matrix is being corrupted, as you
> can see in the file manifestMemProb.incorrect.

No, there is not a bug. Instead, it seems you didn't understand the 
implications of the LaGenMatDouble constructor that you used. The 
LaGenMatDouble(int, int) constructor clearly states: "Matrix elements 
are NOT initialized!" (see 
http://lapackpp.sourceforge.net/html/classLaGenMatDouble.html#z41_1 ) -- 
and that is precisely what you see here. If you want an all-zero matrix, 
you need to do

   LaGenMatDouble A(10,10);
   A = 0.0;

The reason for this is that the extra initialization (setting all 
elements to zero) can be an expensive operation (think of a 1e5x1e5 
matrix), and it really depends on the application whether this is really 
needed. Therefore the standard constructor will only allocate the 
memory, but the values inside that memory can be anything. It is just 
pure coincidence that the first matrix in your test program happened to 
contain "almost" zeros.

Lapackpp is explicitly intended for "high-performance linear algebra", 
and in that field you really need to think for yourself about which 
prerequisites you depend on and which you don't. The library will give 
you many choices, especially all the high-performance choices. 
Construction without initialization clearly is one of those choices.

Christian