In-place Class Functions

[Andrew Morgan]

Hi

I have a matrix class:

template<class Type> class TMatrix
{
//...
TMatrix& T();
//...
}

The transposition function above does not execute in place, ie it
returns a new object. This allows expressions like:

TMatrix<single> A,L,D;
L = ...
D = ...
A = L.T()*D*L;

L.T() is not L and therefore it makes sense to return a new object. If
L.T() were executed in-place, the final product D*L would in fact be
D*L.T() which is incorrect.

However, there are situations where the transpose is required in an
assignment, and

L = L.T();

seems clumsy and inefficient with all the constructors and destructors
getting in on the action (unless good return by value optimizations are
being generated).

Should I do something like:

TMatrix& T(bool inplace = false);

which does not break:
A = L.T()*D*L;

but allows:
L.T(true);

and even:
A = L.T(true)*D*L.T();

as a pathological example. The difference between this expression and
the first version being that L is ends in a transposed state (first
transpose in-place, the second not).

The code may be something like:
template<class Type> TMatrix<Type>& TMatrix<Type>::T(bool inplace)
{
if(inplace)
{
// Transpose in-place.
//...
return *this;
}
/* Create temporary and transpose in-place. Admittedly, a bit lazy
since the transpose code could be re-done below, but on the other hand
one does not want two code fragments doing the same job. */
return TMatrix<Type>(this).T(true);
}

Since a temporary has to be created, maybe
L = L.T();
is not so inefficient after all?

Any comments or suggestions?

Thanks
Andrew

[Julián Albo]

Andrew Morgan escribió:

template<class Type> class TMatrix
{
//...
TMatrix& T();
//...
}

The transposition function above does not execute in place, ie it
returns a new object. This allows expressions like:

TMatrix<single> A,L,D;
L = ...
D = ...
A = L.T()*D*L;

L.T() is not L and therefore it makes sense to return a new object. If
L.T() were executed in-place, the final product D*L would in fact be
D*L.T() which is incorrect.

However, there are situations where the transpose is required in an
assignment, and

L = L.T();

Make T return a TMatrix_T object that holds a reference to the TMatrix.
Define a TMatrix::operator = (const TMatrix_T &) and an operator *
(const TMatrix_T &, const TMatrix &), and additional operators if
neccessary.

Regards.

[lilburne]

Andrew Morgan wrote:

Hi

I have a matrix class:

template<class Type> class TMatrix
{
//...
TMatrix& T();
//...
}

The transposition function above does not execute in place, ie it
returns a new object. This allows expressions like:

TMatrix<single> A,L,D;
L = ...
D = ...
A = L.T()*D*L;

L.T() is not L and therefore it makes sense to return a new object. If
L.T() were executed in-place, the final product D*L would in fact be
D*L.T() which is incorrect.

However, there are situations where the transpose is required in an
assignment, and

L = L.T();

seems clumsy and inefficient with all the constructors and destructors
getting in on the action (unless good return by value optimizations are
being generated).

Should I do something like:

TMatrix& T(bool inplace = false);

which does not break:
A = L.T()*D*L;

but allows:
L.T(true);

and even:
A = L.T(true)*D*L.T();

as a pathological example. The difference between this expression and
the first version being that L is ends in a transposed state (first
transpose in-place, the second not).

The code may be something like:
template<class Type> TMatrix<Type>& TMatrix<Type>::T(bool inplace)
{
if(inplace)
{
// Transpose in-place.
//...
return *this;
}
/* Create temporary and transpose in-place. Admittedly, a bit lazy
since the transpose code could be re-done below, but on the other hand
one does not want two code fragments doing the same job. */
return TMatrix<Type>(this).T(true);
}

Since a temporary has to be created, maybe
L = L.T();
is not so inefficient after all?

Any comments or suggestions?

It seems to me that you are getting into a kerfuffle trying
to provide a unified syntax for doing different things.

T(TMatrix& result) const;
// assign 'result' the transposition of this matrix

TMatrix& T();
// transpose matrix this matrix
anyone reading the code won't then have to wonder what the
signifigance of 'true' is.

[Andrew Morgan]

lilburne wrote:

It seems to me that you are getting into a kerfuffle trying to provide a
unified syntax for doing different things. True.
T(TMatrix& result) const;
// assign 'result' the transposition of this matrix

TMatrix& T();
// transpose matrix this matrix
anyone reading the code won't then have to wonder what the signifigance
of 'true' is.

I now have:

template <class Type> class TMatrix
{
//...
void T();
void PLU();
void UFSub(TMatrix &x);
void BSub(TMatrix &y);
//...
}
template<class Type> TMatrix<Type> T(const TMatrix<Type> &A)
{
// Create temporary, transpose and return.
TMatrix<Type> B(A);
B.T();
return B;
};

The void return types in the class is meant to convey the idea that the
operation are done in-place, while the function T (which does not alter
its argument - this would be catastrophic) returns a new transposed object.

All in-expression transpositions are handled by the function:

A = B + (C * T(D)) + E;

[BTW, is A = B + C * T(D) + E; evaluated left to right, right to left (I
suspect) or according to mathematical precedence where * would be
evaluated first?]

The reason why the in-place/copy issues has come up is because the
function PLU (Partially pivoted LU decomposition) is done in-place.

Solving a finite element-type problem:
[f] = [K][x]

where [K] is a 1000 node 3D problem (and hence 6,000x6,000 = 36,000,000
elements = 275Mb! - worst case obviously) is done in-place simply
because I may not have a spare 275M to create a copy!

Re-arranging:
[x] = [f]/[K] (or strictly, [x] = INV[K][f])

Again, [x] may be "largish" (6,000 elements = 47K) so [f] can be
decomposed in-place into [x].

Coding would look like this:
f = ...
K = ...
f /= K; // [f] is now [x] - the solution.

which wakes up:

template<class Type> TMatrix<Type>& TMatrix<Type>::operator/=(const
TMatrix &K)
{
// Division is strictly inversion and pre-multiplication.
/* Which version???
TMatrix A(K);
A.PLU();
A.UFSub(*this);
A.BSub(*this);
return *this;
*/
// Factorize K to LU by Partial LU.
K.PLU();
// Solve [y] = INV[L][f] by Unit Forward Substitution.
K.UFSub(*this);
// Solve [x] = INV[u][y] by Back Substitution.
K.BSub(*this);
return *this;
};

As you can see from the code, my dilemma is:
* do I make &K const, in which case I must make a temporary A (275M
free?), or do a leave &K non-const (ie developer please note that
operator/= is entitled to modify non-const K, and similarly, non-const
&x and &y in UFSub() and BSub() respectively)?

Thus, f /= K decomposes both f and K in-place.

If I were the only user of the code I can live with this since I know
its behaviour, but it is not necessarily good class design. But then,
good design sometimes runs foul of practical problems. Comments?

Thanks for the posts.

Regards
Andrew

[lilburne]

Andrew Morgan wrote:

lilburne wrote:
Thus, f /= K decomposes both f and K in-place.

If I were the only user of the code I can live with this since I know
its behaviour, but it is not necessarily good class design. But then,
good design sometimes runs foul of practical problems. Comments?

I looked at our own matrix class and unsurprisingly none of
the methods return a matrix by value. Now our matrix is
mainly used for geometric transformation so they never get
beyond 4x4 in size and are usually contained within a class
that provides geometric transforms. However, the design of
the class is such that a user of it never gets any
surprises. We provide some methods so that users can do
Vector * Transform or Point * Transform, but mostly we'll
have methods which look like:

Matrix& Transpose();
Matrix& Transpose(Matrix& copy);
Matrix& Postmultiple(const Matrix&);
Matrix& Postmultiple(const Matrix&, Matrix& result) const;
Matrix& Premultiple(const Matrix);

obviously this forces users to write their code in a
functional style rather than by using operator shorthands,
the reasoning being that the class user is never left in any
doubt about whether something is going to occur in place or not.

So we sidestep the issue of temporaries by saying

"if you require a temporary then you actively create one
yourself, you don't rely on the compiler doing it for
you."

this is a style of coding which we carry through with all
our data structures.

[Andrew Morgan]

lilburne wrote:

obviously this forces users to write their code in a functional style
rather than by using operator shorthands, the reasoning being that the
class user is never left in any doubt about whether something is going
to occur in place or not.

So we sidestep the issue of temporaries by saying

"if you require a temporary then you actively create one
yourself, you don't rely on the compiler doing it for
you."

this is a style of coding which we carry through with all our data
structures.

I agree that is a safe approach.

I am however, trying to capture the essence of object orientation which
urges user-defined types to behave like built-ins.

That is if:
double x,y,z;
x = y + z;
y /= z;

are valid for doubles, then
TMatrix<double> X,Y,Z;
X = Y + Z;
Y /= Z;

should also be (without surprises and bad side effects -- the
responsibility of the designer).

To this end, I have the following global templates:

//---------------------------------------------------------------------------
template<class TClass> const TClass operator+(const TClass &Left,const
TClass &Right)
{
return TClass(Left) += Right;
}
//---------------------------------------------------------------------------
template<class TClass> const TClass operator-(const TClass &Left,const
TClass &Right)
{
return TClass(Left) -= Right;
}
//---------------------------------------------------------------------------
template<class TClass> const TClass operator*(const TClass &Left,const
TClass &Right)
{
return TClass(Left) *= Right;
}
//---------------------------------------------------------------------------
template<class TClass> const TClass operator/(const TClass &Left,const
TClass &Right)
{
return TClass(Left) /= Right;
}
//---------------------------------------------------------------------------

These are always implemented as op= (per Scott Meyers). To support
matrix operations with these template, my dilemma is solved, since I
have had to redefined operator/= with a const parameter, and since the
parameter is decomposed internally, I must necessarily use a temporary:

template<class Type> TMatrix<Type>& TMatrix<Type>::operator/=(const
TMatrix &B)
{
// Division is strictly inversion and pre-multiplication.
TMatrix A(B);
A.PLU();
A.UFSub(*this);
A.BSub(*this);
return *this;
};

So, in this scenario I can match the behaviour of the built-ins (without
surprises), but I have not lost the ability to operations without
temporaries:

TMatrix f,K,x;
f = ...
K = ...
----------
x = f; // Contrived, but f not lost.
x /= K; // Uses temporary for K.
----------
OR
----------
f /= K; // Uses temporary for K, x not needed, but f modified
---------- // (no surprise)
OR
----------
K.Inverse(); // Very slow but OK.
x = K*f; // Solution, but non-intuitive if you did not spot the //
Inverse()
----------
OR
---------- // Fastest!
K.PLU(); // In-place decomposition - no surprise or temporary
K.UFSub(f);
K.BSub(f); // Solution in f
----------
OR
----------
x = f; // Contrived, but f not lost.
K.PLU(); // In-place decomposition - no surprise or temporary
K.UFSub(x);
K.BSub(x); // Solution in x.
----------

The last two could be wrapped in functions:

void Solve(TMatrix &f,const TMatrix &K);
TMatrix& Solve(const TMatrix &f,const TMatrix &K);

respectively.

Having given this a lot of thought recently, built-in behaviour can be
replicated by TMatrix operators which use temporaries. When performance
and memory storage are issues, use functions like x = Solve(f,K).

This strategy differs somewhat from yours in that when one wishes to use
the extremely convenient short-hand, one pays the price in temporaries,
but more traditional alternatives are available should the need arise.

Thanks for the input.

Regards
Andrew

[E. Robert Tisdale]

Andrew Morgan wrote:

I have a matrix class:

template<class Type> class TMatrix {
//...
TMatrix& T(void);
//...
};

The transposition function above does not execute in place, ie it
returns a new object.
No it doesn't. It returns a *reference* to a TMatrix object.
You probably just screwed up. Instead, you should define:

template<class Type> class TMatrix {
//...
TMatrix T(void);
//...
};
This allows expressions like:
TMatrix<single> A,L,D;
L = ...
D = ...
A = L.T()*D*L;

L.T() is not L and therefore it makes sense to return a new object.
If L.T() were executed in-place,
the final product D*L would in fact be D*L.T() which is incorrect.

However, there are situations where the transpose is required in an
assignment, and

L = L.T();

seems clumsy and inefficient with all the constructors and destructors
getting in on the action (unless good return by value optimizations are
being generated).

Should I do something like:

TMatrix& T(bool inplace = false);

which does not break:
A = L.T()*D*L;

but allows:
L.T(true);

and even:
A = L.T(true)*D*L.T();

as a pathological example. The difference between this expression and
the first version being that L is ends in a transposed state (first
transpose in-place, the second not).

The code may be something like:
template<class Type> TMatrix<Type>& TMatrix<Type>::T(bool inplace) {
if(inplace) {
// Transpose in-place.
//...
return *this;
};
/* Create temporary and transpose in-place. Admittedly, a bit lazy
since the transpose code could be re-done below, but on the other hand
one does not want two code fragments doing the same job. */
return TMatrix<Type>(this).T(true);
}

Since a temporary has to be created, maybe
L = L.T();
is not so inefficient after all?

Any comments or suggestions?

Take a look at
The C++ Scalar, Vector, Matrix and Tensor class Library

http://www.netwood.net/~edwin/svmtl/

Remember that the transpose is applied *only* to
*Square* matrix types.

[Andrew Morgan]

E. Robert Tisdale wrote:

Andrew Morgan wrote:

I have a matrix class:

template<class Type> class TMatrix {
//...
TMatrix& T(void);
//...
};

The transposition function above does not execute in place, ie it
returns a new object.

No it doesn't. It returns a *reference* to a TMatrix object.
You probably just screwed up. Instead, you should define:

template<class Type> class TMatrix {
//...
TMatrix T(void);
//...
};

Part of the discussion was simply to debate this point. I have
concluded that member functions that modify their own data should return
void to prevent them being coerced into expressions.

template<class Type> class TMatrix {
//...
void T(); // Transpose.
void PLU(); // Partial Pivot LU Decomposition.
///...
};

and a simple function returning the transpose:

TMatrix T(const TMatrix&);

which can safely be used in expressions:

A = T(B) * D * B;

Take a look at
The C++ Scalar, Vector, Matrix and Tensor class Library

http://www.netwood.net/~edwin/svmtl/
Still chewing...
Remember that the transpose is applied *only* to
*Square* matrix types.

Huh? I make very heavy use of xT*x (dot product) as the work horse of
many of the matrix inner loops using optimized assembler. A significant
application of the class processes huge non-square matrices.

[E. Robert Tisdale]

Andrew Morgan wrote:

Part of the discussion was simply to debate this point.
I have concluded that member functions that modify their own data
should return void to prevent them being coerced into expressions.

template<class Type> class TMatrix {
//...
void T(); // Transpose.
void PLU(); // Partial Pivot LU Decomposition.
///...
};
This is a *very* bad idea.
In place operations *should* return a reference to the modified object
so that it *can* be used in an expression.

For example:

TMatrix<Type>& TMatrix<Type>::operator=(
const TMatrix<Type> rhs) {
/* . . . */
return *this; }

and a simple function returning the transpose:

TMatrix T(const TMatrix&);

which can safely be used in expressions:

A = T(B) * D * B;
No.

Use operator* for element-by-element multiplication.
Use operator/ for element-by-element division.

Take a look at
The C++ Scalar, Vector, Matrix and Tensor class Library

http://www.netwood.net/~edwin/svmtl/

Still chewing...

Remember that the transpose is applied *only* to
*Square* matrix types.

Huh? I make very heavy use of xT*x (dot product) as the work horse
of many of the matrix inner loops using optimized assembler.
A significant application of the class
processes huge non-square matrices.

No!

The matrix-matrix dot product

C = A.dot(B); // C <-- AB^T

is a very special function.

You are obviously a Fortran programmer.
In C++, a matrix is a collection of *row* vectors --
*not* column vectors.

Anyway, this is the wrong place to discuss this.

Try the Object Oriented Numerics Web Page:

http://www.oonumerics.org/oon/

and subscribe to the mailing list.

Meanwhile, study
The C++ Scalar, Vector, Matrix and Tensor class library
standard proposal at

http://www.netwood.net/~edwin/svmtl/