/*
 *
 *  February 21, 2000
 *  Incorporated ISOcxx compatibility
 *  - W. E. Brown
 *
 */


#ifdef ISOCXX__ISOCXX
#  include <ISOcxx/ISOcxx.h>
#endif


#include "LinearAlgebra/Matrix.h"
#include "LinearAlgebra/MLD33.h"
#include "LinearAlgebra/MLC33.h"   // Needed only for cloneComplex
#include "LinearAlgebra/MLDGen.h"  // For the arithmetic function table.
#include "ZMtools/pretendToUse.h"    // Needed for pretendToUse()

#ifdef __BORLAND_CPP__
#include <string.h>
#include <iomanip.h>
#else
#include <cstring>
#include <ZMutility/iomanip>
using namespace std;
#endif // __BORLAND_CPP__

// ===
// === Begin MLD33.icc file ==============================
// === This is temporarily put here for diagnostics ======
// ===

/* inline */ MLD* MLD33::EmptyClone() const
{
  return new MLD33();
}

// ??? REMOVE /* inline */ const Type_info& MLD33::typeID()
// ??? REMOVE {
// ??? REMOVE    return typeid( *this );
// ??? REMOVE }

/* inline */ int MLD33::vtblIndex() const
{
  return FPCL_INDEX_D_M33;
}

// ===
// === End MLD33.icc file ================================
// === This is temporarily put here for diagnostics ======
// ===

// ===
// === Declaration function ==============================
// ===

MatrixD& MatrixD::declareM33()
{
  // ??? REMOVE if( typeid( *_pml ) != typeid( MLD33 ) ) {
  // ??? REMOVE   if( --( _pml->_rc ) == 0 ) delete _pml;
  // ??? REMOVE   _pml = new MLD33();
  // ??? REMOVE }
  // ??? REMOVE return *this;
  // For a NxN data model, there is no difference from .makeMNN()
  return makeM33();
}


MatrixD& MatrixD::makeM33()
{
  static MLD33* pNewMLD;

  if( typeid( *_pml ) != typeid( MLD33 ) ) {
    if( ( _pml->_r == 3 ) && ( _pml->_c == 3 ) ) {
      pNewMLD = new MLD33();

      pNewMLD->_m[0][0] = _pml->value(0,0);
      pNewMLD->_m[0][1] = _pml->value(0,1);
      pNewMLD->_m[0][2] = _pml->value(0,2);

      pNewMLD->_m[1][0] = _pml->value(1,0);
      pNewMLD->_m[1][1] = _pml->value(1,1);
      pNewMLD->_m[1][2] = _pml->value(1,2);

      pNewMLD->_m[2][0] = _pml->value(2,0);
      pNewMLD->_m[2][1] = _pml->value(2,1);
      pNewMLD->_m[2][2] = _pml->value(2,2);

      if( --( _pml->_rc ) == 0 ) delete _pml;
      _pml = pNewMLD;
    }
    else { 
      ZMthrow ( ZMxlaDimensionMismatch(
                "MatrixD::makeM33 "
                "Dimensions are not correct") );
    }
  }
  return *this;
}

// ===
// === Two-way virtual table entries =====================
// ===

MLD* SumM33M33( const MLD* argx, const MLD* argy )
{
  static MLD33*  pNew;
  static Float8* px;
  static Float8* py;
  static Float8* pz;

  pNew = new MLD33( 0.0 );
  pz = pNew->_pdi;
  px = argx->_pdi;
  py = argy->_pdi;
  while( pz <= pNew->_pdf ) *(pz++) = *(px++) + *(py++);

  return pNew;
}

MLD* DiffM33M33( const MLD* argx, const MLD* argy )
{
  static MLD33*  pNew;
  static Float8* px;
  static Float8* py;
  static Float8* pz;

  pNew = new MLD33( 0.0 );
  pz = pNew->_pdi;
  px = argx->_pdi;
  py = argy->_pdi;
  while( pz <= pNew->_pdf ) *(pz++) = *(px++) - *(py++);

  return pNew;
}

MLD* ProdM33M33( const MLD* argx, const MLD* argy )
{
  static MLD33* pNew;
  static const MLD33* x;
  static const MLD33* y;

#ifndef NO_DYNAMIC_CAST
  x = dynamic_cast<const MLD33*>( argx );
  y = dynamic_cast<const MLD33*>( argy );
#endif
#ifdef NO_DYNAMIC_CAST
  x = (MLD33*)( argx );
  y = (MLD33*)( argy );
#endif

  pNew = new MLD33( 0.0 );

  pNew->_m[0][0] = x->_m[0][0] * y->_m[0][0]
                 + x->_m[0][1] * y->_m[1][0]
                 + x->_m[0][2] * y->_m[2][0];

  pNew->_m[0][1] = x->_m[0][0] * y->_m[0][1]
                 + x->_m[0][1] * y->_m[1][1]
                 + x->_m[0][2] * y->_m[2][1];

  pNew->_m[0][2] = x->_m[0][0] * y->_m[0][2]
                 + x->_m[0][1] * y->_m[1][2]
                 + x->_m[0][2] * y->_m[2][2];

  pNew->_m[1][0] = x->_m[1][0] * y->_m[0][0]
                 + x->_m[1][1] * y->_m[1][0]
                 + x->_m[1][2] * y->_m[2][0];

  pNew->_m[1][1] = x->_m[1][0] * y->_m[0][1]
                 + x->_m[1][1] * y->_m[1][1]
                 + x->_m[1][2] * y->_m[2][1];

  pNew->_m[1][2] = x->_m[1][0] * y->_m[0][2]
                 + x->_m[1][1] * y->_m[1][2]
                 + x->_m[1][2] * y->_m[2][2];

  pNew->_m[2][0] = x->_m[2][0] * y->_m[0][0]
                 + x->_m[2][1] * y->_m[1][0]
                 + x->_m[2][2] * y->_m[2][0];

  pNew->_m[2][1] = x->_m[2][0] * y->_m[0][1]
                 + x->_m[2][1] * y->_m[1][1]
                 + x->_m[2][2] * y->_m[2][1];

  pNew->_m[2][2] = x->_m[2][0] * y->_m[0][2]
                 + x->_m[2][1] * y->_m[1][2]
                 + x->_m[2][2] * y->_m[2][2];

  return pNew;
}


MLD* SumGenM33( const MLD* argx, const MLD* argy )
{
  static MLDGeneric* pNew;
  static Float8* px;
  static Float8* py;
  static Float8* pz;

  if( ( argx->_r != 3 ) ||
      ( argx->_c != 3 )    )
  {
    return 0;
  }

  pNew = new MLDGeneric( 3, 3, 0.0 );

  pz = pNew->_pdi;
  px = argx->_pdi;
  py = argy->_pdi;
  while( pz <= pNew->_pdf ) *(pz++) = *(px++) + *(py++);

  return pNew;
}


MLD* SumM33Gen( const MLD* argy, const MLD* argx )
{
  // This is a direct copy of SumGenM33
  static MLDGeneric* pNew;
  static Float8* px;
  static Float8* py;
  static Float8* pz;

  if( ( argx->_r != 3 ) ||
      ( argx->_c != 3 )    )
  {
    return 0;
  }

  pNew = new MLDGeneric( 3, 3, 0.0 );

  pz = pNew->_pdi;
  px = argx->_pdi;
  py = argy->_pdi;
  while( pz <= pNew->_pdf ) *(pz++) = *(px++) + *(py++);

  return pNew;
}

MLD* DiffGenM33( const MLD* argx, const MLD* argy )
{
  static MLDGeneric* pNew;
  static Float8* px;
  static Float8* py;
  static Float8* pz;

  if( ( argx->_r != 3 ) ||
      ( argx->_c != 3 )    )
  {
    return 0;
  }

  pNew = new MLDGeneric( 3, 3, 0.0 );

  pz = pNew->_pdi;
  px = argx->_pdi;
  py = argy->_pdi;
  while( pz <= pNew->_pdf ) *(pz++) = *(px++) - *(py++);

  return pNew;
}


MLD* DiffM33Gen( const MLD* argy, const MLD* argx )
{
  // This is a direct copy of DiffGenM33
  // except interchanging y and x at the end.
  static MLDGeneric* pNew;
  static Float8* px;
  static Float8* py;
  static Float8* pz;

  if( ( argx->_r != 3 ) ||
      ( argx->_c != 3 )    )
  {
    return 0;
  }

  pNew = new MLDGeneric( 3, 3, 0.0 );

  pz = pNew->_pdi;
  px = argx->_pdi;
  py = argy->_pdi;
  while( pz <= pNew->_pdf ) *(pz++) = *(py++) - *(px++);

  return pNew;
}

MLD* ProdGenM33( const MLD* argx, const MLD* argy )
{
  static MLDGeneric* pNew;
  static const MLDGeneric* x;
  static const MLD33* y;
  static int i;
  static int r;

  if( argx->_c != 3 )
  {
    return 0;
  }

  r = argx->_r;
#ifndef NO_DYNAMIC_CAST
  x = dynamic_cast<const MLDGeneric*>( argx );
  y = dynamic_cast<const MLD33*>( argy );
#endif
#ifdef NO_DYNAMIC_CAST
  x = (MLDGeneric*)( argx );
  y = (MLD33*)( argy );
#endif

  pNew = new MLDGeneric( r, 3, 0.0 );

  for( i = 0; i < r; i++ )  {
    pNew->_m[i][0] = x->_m[i][0] * y->_m[0][0]
                   + x->_m[i][1] * y->_m[1][0]
                   + x->_m[i][2] * y->_m[2][0];
    pNew->_m[i][1] = x->_m[i][0] * y->_m[0][1]
                   + x->_m[i][1] * y->_m[1][1]
                   + x->_m[i][2] * y->_m[2][1];
    pNew->_m[i][2] = x->_m[i][0] * y->_m[0][2]
                   + x->_m[i][1] * y->_m[1][2]
                   + x->_m[i][2] * y->_m[2][2];
  }

  return pNew;
}

MLD* ProdM33Gen( const MLD* argy, const MLD* argx )
{
  static MLDGeneric* pNew;
  static const MLDGeneric* x;
  static const MLD33* y;
  static int j;
  static int c;

  if( argx->_r!= 3 )
  {
    return 0;
  }

  c = argx->_c;
#ifndef NO_DYNAMIC_CAST
  x = dynamic_cast<const MLDGeneric*>( argx );
  y = dynamic_cast<const MLD33*>( argy );
#endif
#ifdef NO_DYNAMIC_CAST
  x = (MLDGeneric*)( argx );
  y = (MLD33*)( argy );
#endif

  pNew = new MLDGeneric( 3, c, 0.0 );

  for( j = 0; j < c; j++ )  {
    pNew->_m[0][j] = y->_m[0][0] * x->_m[0][j]
                   + y->_m[0][1] * x->_m[1][j]
                   + y->_m[0][2] * x->_m[2][j];

    pNew->_m[1][j] = y->_m[1][0] * x->_m[0][j]
                   + y->_m[1][1] * x->_m[1][j]
                   + y->_m[1][2] * x->_m[2][j];

    pNew->_m[2][j] = y->_m[2][0] * x->_m[0][j]
                   + y->_m[2][1] * x->_m[1][j]
                   + y->_m[2][2] * x->_m[2][j];
  }

  return pNew;
}


int MLD33::initOps()
{
  SumVtbl[ FPCL_INDEX_D_M33 ][ FPCL_INDEX_D_M33 ] = SumM33M33;
  SumVtbl[ FPCL_INDEX_D_GEN ][ FPCL_INDEX_D_M33 ] = SumGenM33;
  SumVtbl[ FPCL_INDEX_D_M33 ][ FPCL_INDEX_D_GEN ] = SumM33Gen;

  DiffVtbl[ FPCL_INDEX_D_M33 ][ FPCL_INDEX_D_M33 ] = DiffM33M33;
  DiffVtbl[ FPCL_INDEX_D_GEN ][ FPCL_INDEX_D_M33 ] = DiffGenM33;
  DiffVtbl[ FPCL_INDEX_D_M33 ][ FPCL_INDEX_D_GEN ] = DiffM33Gen;

  ProdVtbl[ FPCL_INDEX_D_M33 ][ FPCL_INDEX_D_M33 ] = ProdM33M33;
  ProdVtbl[ FPCL_INDEX_D_GEN ][ FPCL_INDEX_D_M33 ] = ProdGenM33;
  ProdVtbl[ FPCL_INDEX_D_M33 ][ FPCL_INDEX_D_GEN ] = ProdM33Gen;

  return 1;
}



// ===
// === MLD33   constructor ===============================
// ===

MLD33::MLD33( const Float8& initval )
: MLD( 3, 3 )
{
  // Initialize the entries into the two-way virtual function table
  // This should be done once and once only.
  static int opsDone = initOps();

  // pretend to use the variable in order to avoid unused variable
  // warnings at compile time
  pretendToUse(opsDone);  
  
             _m[0][1] = _m[0][2] = 0.0;
  _m[1][0]            = _m[1][2] = 0.0;
  _m[2][0] = _m[2][1]            = 0.0;

  _m[0][0] = _m[1][1] = _m[2][2] = initval;

  _pdi = _m[0];
  _pdf = &( _m[2][2] );
  _size = 9*sizeof( Float8 );
}


MLD33::MLD33( const Float8* initval )
: MLD( 3, 3 )
{
  // Initialize the entries into the two-way virtual function table
  // This should be done once and once only.
  static int opsDone = initOps();

  // pretend to use the variable in order to avoid unused variable
  // warnings at compile time
  pretendToUse(opsDone);  
  
  _pdi = _m[0];
  _pdf = &( _m[2][2] );
  _size = 9*sizeof( Float8 );

  _m[0][0] = initval[0]; // ??? Use memcpy here.
  _m[0][1] = initval[1];
  _m[0][2] = initval[2];
  _m[1][0] = initval[3];
  _m[1][1] = initval[4];
  _m[1][2] = initval[5];
  _m[2][0] = initval[6];
  _m[2][1] = initval[7];
  _m[2][2] = initval[8];

}


void MLD33::loadFrom( const MLD* x )
{
  static Float8* p;
  if( ( _r == x->_r ) && ( _c == x->_c ) ) {
    p = _pdi;

    *(p++) = x->value(0,0);
    *(p++) = x->value(0,1);
    *(p++) = x->value(0,2);

    *(p++) = x->value(1,0);
    *(p++) = x->value(1,1);
    *(p++) = x->value(1,2);

    *(p++) = x->value(2,0);
    *(p++) = x->value(2,1);
    *p++   = x->value(2,2);
  }

  // Dimensions were not correct.
  else { 
    ZMthrow ( ZMxlaDimensionMismatch(
              "MLD33::loadFrom( MLD* ) "
              "Dimensions do not match") );
  }
}


bool MLD33::setValue( int i, int j, const Float8& x )
{
  // Normally, the bounds check should be done
  // at a higher level.
  #if defined(ZM_MLD_DEBUG) || defined(ZM_MLD_TESTS)
  if( ( i < 0 ) || ( i > 2 ) ||
      ( j < 0 ) || ( j > 2 )
    ) { 
    ZMthrow ( ZMxlaInvalidIndex(
              "MLD33::setValue "
              "Indices out of bounds") );
    return 1;
  };
  #endif

  _m[i][j] = x;
  return 0;
}


Float8 MLD33::value( int i, int j ) const
{
  return _m[i][j];
}


MLD* MLD33::Clone() const
{
  static MLD*     q;
  // ??? REMOVE static Float8* py;
  // ??? REMOVE static Float8* px;
  q = new MLD33();
  // ??? REMOVE py = q->_pdi;
  // ??? REMOVE px =    _pdi;
  // ??? REMOVE while( px <= _pdf ) *(py++) = *(px++);  // ??? Use memcpy here.
  ::memcpy( q->_pdi, _pdi, _size );
  return q;
}


MLC* MLD33::cloneComplex() const
{
  static MLC* q;
  static Complex8* pf;
  static Float8* pi;
  static int i;

  q = new MLC33();
  pi = _pdi;
  pf = q->_pdi;
  for( i = 0; i < 9; i++, pf++, pi++ ) {
    *pf = *pi;
  }

  return q;
}


// ===
// === Utility functions =================================
// ===

Float8 MLD33::determinant() const
{
  return ( _m[0][0]*( _m[1][1]*_m[2][2] - _m[1][2]*_m[2][1] )
         + _m[0][1]*( _m[1][2]*_m[2][0] - _m[1][0]*_m[2][2] )
         + _m[0][2]*( _m[1][0]*_m[2][1] - _m[1][1]*_m[2][0] )
         );
}

Float8 MLD33::trace() const
{
  return ( _m[0][0] + _m[1][1] + _m[2][2] );
}

MLD* MLD33::transpose() const
{
  static MLD33* ret;
  ret = new MLD33;
  ret->_m[0][0] = _m[0][0];
  ret->_m[0][1] = _m[1][0];
  ret->_m[0][2] = _m[2][0];

  ret->_m[1][0] = _m[0][1];
  ret->_m[1][1] = _m[1][1];
  ret->_m[1][2] = _m[2][1];

  ret->_m[2][0] = _m[0][2];
  ret->_m[2][1] = _m[1][2];
  ret->_m[2][2] = _m[2][2];
  return ret;
}


// ===
// === Inversion, solving linear equations ===============
// ===

MLD* MLD33::inverse() const
{
  static MLD33* ret;
  static Float8 det;

  // Calculation of determinant done in place rather than
  // calling MLD33::determinant(). Is this stupid?
  det =  _m[0][0]*( _m[1][1]*_m[2][2] - _m[1][2]*_m[2][1] )
       + _m[0][1]*( _m[1][2]*_m[2][0] - _m[1][0]*_m[2][2] )
       + _m[0][2]*( _m[1][0]*_m[2][1] - _m[1][1]*_m[2][0] );

  // ??? Should the test be made: fabs(det) <= epsilon ???
  if( det == 0.0 ) { 
    ZMthrow ( ZMxlaSingularMatrix(
              "MLD33::inverse() "
              "Matrix is singular") );
    return 0;
  }

  ret = new MLD33;

  ret->_m[0][0] = ( _m[1][1]*_m[2][2] - _m[1][2]*_m[2][1] ) / det;
  ret->_m[1][0] = ( _m[1][2]*_m[2][0] - _m[1][0]*_m[2][2] ) / det;
  ret->_m[2][0] = ( _m[1][0]*_m[2][1] - _m[1][1]*_m[2][0] ) / det;

  ret->_m[0][1] = ( _m[2][1]*_m[0][2] - _m[2][2]*_m[0][1] ) / det;
  ret->_m[1][1] = ( _m[2][2]*_m[0][0] - _m[2][0]*_m[0][2] ) / det;
  ret->_m[2][1] = ( _m[2][0]*_m[0][1] - _m[2][1]*_m[0][0] ) / det;

  ret->_m[0][2] = ( _m[0][1]*_m[1][2] - _m[0][2]*_m[1][1] ) / det;
  ret->_m[1][2] = ( _m[0][2]*_m[1][0] - _m[0][0]*_m[1][2] ) / det;
  ret->_m[2][2] = ( _m[0][0]*_m[1][1] - _m[0][1]*_m[1][0] ) / det;

  return ret;
}



// ===
// === Eigen-analysis ====================================
// ===

MatrixEigenData MLD33::eigen() const
{
  MatrixEigenData ret( 3 );
  MatrixC z( 3, 1 );
  MatrixD c( 4, 1 );

  // Solve characteristic polynomial ...
  c( 0, 0 ) = determinant();
  {
  MatrixD temp( 3, _pdi );
  c( 1, 0 ) = - ( temp.cofactor(0,0) +
                  temp.cofactor(1,1) +
                  temp.cofactor(2,2) );
  } // ... temp goes out of scope
  c( 2, 0 ) = trace();
  c( 3, 0 ) = -1.0;
  MatrixEigenData::_polySolve( c, z );

  // Load eigenvalues ...
  ret._values( 0, 0 ) = z( 0, 0 );
  ret._values( 1, 1 ) = z( 1, 0 );
  ret._values( 2, 2 ) = z( 2, 0 );

  // Solve and load eigenvectors ...
  int i, j, ii, jj, iii, jjj, k;
  Float8 amp;
  MatrixC v( 2, 1 );
  MatrixC aux( 3, 3 );
  aux.declareM33();
  MatrixC x( 2, 1 );
  MatrixC u( 2, 2 );
  u.declareM22();

  for( k = 0; k < 3; k++ )
  {
    // Construct auxiliary matrix and find
    // largest cofactor

    for( i = 0; i < 3; i++ ) for( j = 0; j < 3; j++ ) {
      aux( i, j ) = _m[i][j];
      if( i == j ) aux( i, j ) -= z( k, 0 );
    }

    amp = 0.0;
    ii = 0;
    jj = 0;
    for( i = 0; i < 3; i++ ) for( j = 0; j < 3; j++ ) {
      if( abs( aux.cofactor( i, j ) ) > amp ) {
        ii = i;
        jj = j;
      }
    }


    // Set up and solve linear equation for eigenvectors ...
    iii = 0;
    for( i = 0; i < 3; i++ ) {
      if( i != ii ) {
        v( iii, 0 ) = - aux( i, jj );
        iii++;
      }
    }

    iii = 0;
    for( i = 0; i < 3; i++ ) {
      if( i != ii ) {
        jjj = 0;
        for( j = 0; j < 3; j++ ) {
          if( j != jj ) {
            u( iii, jjj ) = aux( i, j );
            jjj++;
          }
        }
        iii++;
      }
    }

    // ... Solution step ...
    x = u.inverse() * v;


    // Load eigenvectors ...
    iii = 0;
    for( i = 0; i < 3; i++ ) {
       if( i != jj )  {
         ret._vectors( i, k ) = x( iii, 0 );
         iii++;
       }
       else {
         ret._vectors( i, k ) = (Complex8)1.0;
       }
    }

    // Normalize ...
    amp = norm( Complex8(ret._vectors( 0, k )) )
        + norm( Complex8(ret._vectors( 1, k )) )
        + norm( Complex8(ret._vectors( 2, k )) );
    amp = ::sqrt( amp );
    ret._vectors( 0, k ) /= amp;
    ret._vectors( 1, k ) /= amp;
    ret._vectors( 2, k ) /= amp;

  }

  return ret;
}



// ===
// === Row, column, and element operations ===============
// ===

void MLD33::switchRows( int i, int j )
{
  static Float8 temp;
  temp     = _m[i][0];
  _m[i][0] = _m[j][0];
  _m[j][0] = temp;
  temp     = _m[i][1];
  _m[i][1] = _m[j][1];
  _m[j][1] = temp;
  temp     = _m[i][2];
  _m[i][2] = _m[j][2];
  _m[j][2] = temp;
}


void MLD33::switchColumns( int i, int j )
{
  static Float8 temp;
    temp = _m[0][i];
    _m[0][i] = _m[0][j];
    _m[0][j] = temp;
    temp = _m[1][i];
    _m[1][i] = _m[1][j];
    _m[1][j] = temp;
    temp = _m[2][i];
    _m[2][i] = _m[2][j];
    _m[2][j] = temp;
}



// ===
// === Remaining MLDVirtF.h functions ====================
// ===

void MLD33::opPlusEqFlt  ( const Float8& x )
{
  _m[0][0] += x;
  _m[1][1] += x;
  _m[2][2] += x;
}

void MLD33::opSubEqFlt   ( const Float8& x )
{
  _m[0][0] -= x;
  _m[1][1] -= x;
  _m[2][2] -= x;
}

void MLD33::opMultEqMLD  ( const MLD* x )
{
  static Float8 temp[3];

  temp[0] = _m[0][0];
  temp[1] = _m[0][1];
  temp[2] = _m[0][2];
  (( _m[0][0] = temp[0] * x->value( 0, 0 ) )
             += temp[1] * x->value( 1, 0 ) )
             += temp[2] * x->value( 2, 0 );
  (( _m[0][1] = temp[0] * x->value( 0, 1 ) )
             += temp[1] * x->value( 1, 1 ) )
             += temp[2] * x->value( 2, 1 );
  (( _m[0][2] = temp[0] * x->value( 0, 2 ) )
             += temp[1] * x->value( 1, 2 ) )
             += temp[2] * x->value( 2, 2 );

  temp[0] = _m[1][0];
  temp[1] = _m[1][1];
  temp[2] = _m[1][2];
  (( _m[1][0] = temp[0] * x->value( 0, 0 ) )
             += temp[1] * x->value( 1, 0 ) )
             += temp[2] * x->value( 2, 0 );
  (( _m[1][1] = temp[0] * x->value( 0, 1 ) )
             += temp[1] * x->value( 1, 1 ) )
             += temp[2] * x->value( 2, 1 );
  (( _m[1][2] = temp[0] * x->value( 0, 2 ) )
             += temp[1] * x->value( 1, 2 ) )
             += temp[2] * x->value( 2, 2 );

  temp[0] = _m[2][0];
  temp[1] = _m[2][1];
  temp[2] = _m[2][2];
  (( _m[2][0] = temp[0] * x->value( 0, 0 ) )
             += temp[1] * x->value( 1, 0 ) )
             += temp[2] * x->value( 2, 0 );
  (( _m[2][1] = temp[0] * x->value( 0, 1 ) )
             += temp[1] * x->value( 1, 1 ) )
             += temp[2] * x->value( 2, 1 );
  (( _m[2][2] = temp[0] * x->value( 0, 2 ) )
             += temp[1] * x->value( 1, 2 ) )
             += temp[2] * x->value( 2, 2 );

}

MLD* MLD33::add( const Float8& x ) const
{
  static MLD33*   ret;

#ifdef NO_DYNAMIC_CAST
  ret = (MLD33*) ( this->Clone() );
#else
  ret = dynamic_cast<MLD33*>( this->Clone() );
#endif

  ret->_m[0][0] += x;
  ret->_m[1][1] += x;
  ret->_m[2][2] += x;
  return ret;
}

MLD* MLD33::subtract( const Float8& x ) const
{
  static MLD33*   ret;

#ifdef NO_DYNAMIC_CAST
  ret = (MLD33*) ( this->Clone() );
#else
  ret = dynamic_cast<MLD33*>( this->Clone() );
#endif

  ret->_m[0][0] -= x;
  ret->_m[1][1] -= x;
  ret->_m[2][2] -= x;
  return ret;
}

void MLD33::readFrom( istream& is )
{
#include "LinearAlgebra/MLDREADFROM.ins"
}