// ----------------------------------------------------------------------
//
// la_pv_test -- test all features of the LinearAlgebra - PhysicsVectors
//	            bridge
//
// History:
//   25-Jun-1998  KLS   Initial draft created
//
// ----------------------------------------------------------------------

#ifdef __GNUG__
#include "ZMutility/iostream"
using std::cout;

int main() {
  cout << "LinearAlgebra bridge header currently unavailable for GCC.\n"
       << "Testing nothing.\n";
  return 0;
}

#else
#include "PhysicsVectors/LinAlg-PhysVec.h"


#include "PhysicsVectors/PlaneVector.h"
#include "PhysicsVectors/SpaceVector.h"
#include "PhysicsVectors/LorentzVector.h"
#include "PhysicsVectors/Rotation.h"
#include "PhysicsVectors/LorentzTransformation.h"

#include "LinearAlgebra/Matrix.h"

ZM_USING_NAMESPACE( zmpv )

#include "ZMutility/cmath"
using CMATH_NAMESPACE::fabs;
using CMATH_NAMESPACE::sqrt;

#include "ZMutility/iostream"
using std::cout;
using std::endl;

#include "Exceptions/ZMexHandler.h"

// Parts of this test:

void pv2la_constructors();
void pv2la_multiplication();
void la2pv_constructors();
void la2pv_multiplication();
void specialtycases();
void dotproducts();
void sizetests();


bool problems( false );
const Float8 PI =  3.1415926535897931;
const Float8 tolerance = 1.E-14;


int main() {

  ZMxPhysicsVectors::setHandler ( ZMexIgnoreAlways() );

  cout << "\n *************************************"
       << "\n ***** Testing PhysVec to LinAlg *****"
       << "\n *************************************\n";

  pv2la_constructors();
  pv2la_multiplication();

  cout << "\n *************************************"
       << "\n ***** Testing LinAlg to PhysVec *****"
       << "\n *************************************\n";

  la2pv_constructors();
  la2pv_multiplication();

  cout << "\n *************************************"
       << "\n *****   Miscellaneous testing   *****"
       << "\n *************************************\n";

  specialtycases();
  dotproducts();
  sizetests();

  if (problems) {
    cout << "\n ?????????????????? -- End of tests:  PROBLEMS DETECTED \n\n";
    return -1;
  } else {
    cout << "\nEnd of tests -- All automated checking passed.\n\n";
    return 0;
  }

} /* End of la_pv_test main() */



bool check( MatrixD & m, Float8 * ptr) {
  for (int i=0; i < m.rows(); ++i) {
    for (int j=0; j < m.columns(); ++j) {
      if ( fabs( m(i,j) - ptr[i*m.columns()+j] ) > tolerance ) {
         cout << "Element (" << i << "," << j
              << ") doesn't match by "
              << fabs( m(i,j) - ptr[i*m.columns()+j] ) << endl;
         return false;
      }
    }
  }
  return true;
}

Float8 * multiply( MatrixD & m, Float8 * ptr, int cols ) {

  Float8 * result = new Float8[m.rows() * cols];
  for (int a=0; a < m.rows() * cols; ++a) { result[a] = 0; }

  for ( int i=0; i < m.rows(); ++i ) {
    for ( int j=0; j < cols; ++j ) {
       for ( int k=0; k < m.columns(); ++k ) {
          result[i*cols+j] += m(i,k) * ptr[k*cols+j];
       }
    }
  }
  return result;
}


Float8 * multiply( Float8 * ptr, MatrixD & m, int rows ) {

  Float8 * result = new Float8[rows * m.columns()];
  for (int a=0; a < rows * m.columns(); ++a) { result[a] = 0; }

  for ( int i=0; i < rows; ++ i ) {
    for ( int j=0; j < m.columns(); ++j ) {
      for ( int k=0; k < m.rows(); ++k ) {
        result[i*m.columns()+j] += ptr[i*m.rows()+k] * m(k,j);
      }
    }
  }
  return result;
}


void pv2la_constructors() {

  cout << "\n{**** Testing Construction of matrices out of Vectors ****}\n";

  PlaneVector p1( 1.4, 0.9 );
  SpaceVector v1( 1.0, 6.5, 9.1 );
  LorentzVector w1( 0.6, -.4, 1.8, Tcomponent(15) );

  Rotation r1( PI/6, -PI/2, PI);

  LorentzBoost b1( .5, .7, .2);
  LorentzTransformation lt1( b1, r1 );

  Float8 floatArray[4];

  MatrixD m1 = mat( p1 );
  floatArray[0] = p1.x(); floatArray[1] = p1.y();
  if ( check( m1, (Float8 *)(&floatArray) ) ) {
     cout << "OK: Matrix m1 is same as PlaneVector p1 \n";
  } else {
     cout << "??? Matrix m1 isn't same as PlaneVector p1 \n";
     cout << m1 << endl << p1 << endl;
     problems = true;
  }

  MatrixD m2 = mat( v1 );
  floatArray[0] = v1.x(); floatArray[1] = v1.y();
  floatArray[2] = v1.z();
  if ( check( m2, (Float8 *)(&floatArray) ) ) {
    cout << "OK: Matrix m2 is same as SpaceVector v1 \n";
  } else {
    cout << "??? Matrix m2 isn't same as SpaceVector v1 \n";
    cout << m2 << endl << v1 << endl;
    problems = true;
  }

  MatrixD m3 = mat( w1 );
  floatArray[0] = w1.x(); floatArray[1] = w1.y();
  floatArray[2] = w1.z(); floatArray[3] = w1.t();
  if ( check( m3, (Float8 *)(&floatArray) ) ) {
    cout << "OK: Matrix m3 is same as LorentzVector w1 \n";
  } else {
    cout << "??? Matrix m3 isn't same as LorentzVector w1 \n";
    cout << m3 << endl << w1 << endl;
    problems = true;
  }

  MatrixD m4 = mat( r1 );
  ZMpvRep3x3 rep3 = r1.rep3x3();
  if ( check( m4, (Float8 *)(&rep3) ) ) {
    cout << "OK: Matrix m4 is same as Rotation r1 \n";
  } else {
    cout << "??? Matrix m4 isn't same as Rotation r1 \n";
    cout << m4 << endl << r1 << endl;
    problems = true;
  }

  MatrixD m5 = mat( lt1 );
  ZMpvRep4x4 rep4 = lt1.rep4x4();
  if ( check( m5, (Float8 *)(&rep4) ) ) {
    cout << "OK: Matrix m5 is same as LorentzTransformation lt1 \n";
  } else {
    cout << "??? Matrix m5 isn't same as LorentzTransformation lt1 \n";
    cout << m5 << endl << lt1 << endl;
    problems = true;
  }

} // pv2la_constructors()


void pv2la_multiplication() {

  cout << "\n{**** Testing multiplication of Vectors by matrices ****}\n";

  Float8  twoNums[2] =  { 3.0, 4.5 };
  Float8 threNums[3] =  { 1.0, -6.7, 8.1 };
  Float8 fourNums[4] =  { 9.1, 6.2, 5.2, -13 };
  Float8 nineNums[9] =  {  0.5         ,  0.5         , sqrt(2.0)/2,
                          -0.5         , -0.5         , sqrt(2.0)/2,
                           sqrt(2.0)/2 , -sqrt(2.0)/2 ,     0.0
                        };
  Float8 sxtnNums[16] = { 1, 0, 0, 0,
                          0, 1, 0, 0,
                          0, 0, 1, 0,
                          0, 0, 0, 1
                        };

  Float8 floatArray[4];


  MatrixD m1x2(1, 2, twoNums);
  PlaneVector p1( 1.4, 0.9 );

  MatrixD m1 = m1x2 * p1;
  floatArray[0] = p1.x(); floatArray[1] = p1.y();
  if ( check(m1, multiply( m1x2, (Float8 *)(&floatArray), 1) ) ) {
    cout << "OK: Result matrix is correct for multiplication of a matrix"
         << " by a PlaneVector" << endl;
  } else {
    cout << "??? Result matrix incorrect for multiplication of a matrix"
         << " by a PlaneVector" << endl;
    cout << m1 << endl << m1x2 << endl << p1 << endl;
    problems = true;
  }



  MatrixD m1x3(1, 3, threNums);
  SpaceVector v1( 1.0, 6.5, 9.1 );

  MatrixD m2 = m1x3 * v1;
  floatArray[0] = v1.x(); floatArray[1] = v1.y();
  floatArray[2] = v1.z();
  if ( check(m2, multiply(m1x3, (Float8 *)(&floatArray), 1) ) ) {
    cout << "OK: Result matrix is correct for multiplication of a matrix"
         << " by a SpaceVector" << endl;
  } else {
    cout << "??? Result matrix incorrect for multiplication of a matrix"
         << " by a SpaceVector" << endl;
    cout << m1 << endl << m1x3 << endl << v1 << endl;
    problems = true;
  }



  MatrixD m1x4(1, 4, fourNums);
  LorentzVector w1( 0.6, -.4, 1.8, Tcomponent(15) );

  MatrixD m3 = m1x4 * w1;
  floatArray[0] = w1.x(); floatArray[1] = w1.y();
  floatArray[2] = w1.z(); floatArray[3] = w1.t();
  if ( check(m3, multiply(m1x4, (Float8 *)(&floatArray), 1) ) ) {
    cout << "OK: Result matrix is correct for multiplication of a matrix"
         << " by a LorentzVector" << endl;
  } else {
    cout << "??? Result matrix incorrect for multiplication of a matrix"
         << " by a LorentzVector" << endl;
    cout << m3 << endl << m1x4 << endl << w1 << endl;
    problems = true;
  }



  MatrixD m3x3(3, 3, nineNums);
  Rotation r1( PI/6, -PI/2, PI);

  MatrixD m4 = m3x3 * r1;
  ZMpvRep3x3 rep3 = r1.rep3x3();
  if ( check(m4, multiply(m3x3, (Float8 *)(&rep3), 3) ) ) {
    cout << "OK: Result matrix is correct for multiplication of a matrix"
         << " by a Rotation" << endl;
  } else {
    cout << "??? Result matrix incorrect for multiplication of a matrix"
         << " by a Rotation" << endl;
    cout << m4 << endl << m3x3 << endl << r1 << endl;
    problems = true;
  }



  MatrixD m4x4(4, 4, sxtnNums);
  LorentzBoost b1( .5, .7, .2);
  LorentzTransformation lt1( b1, r1 );

  MatrixD m5 = m4x4 * lt1;
  ZMpvRep4x4 rep4 = lt1.rep4x4();
  if ( check(m5, multiply(m4x4, (Float8 *)(&rep4), 4) ) ) {
    cout << "OK: Result matrix is correct for multiplication of a matrix"
         << " by a LorentzTransformation" << endl;
  } else {
    cout << "??? Result matrix incorrect for multiplication of a matrix"
         << " by a LorentzTransformation" << endl;
    cout << m5 << endl << m4x4 << endl << lt1 << endl;
    problems = true;
  }

} // pv2la_multiplication()


void la2pv_constructors() {

  cout << "\n{**** Testing Construction of Vectors from matrices ****}\n";

  Float8 twoNums[2] = { 3.0, 4.5 };
  MatrixD m1 (2, 1, twoNums);
  MatrixD m1a(1, 2, twoNums);

  Float8 threeNums[3] = { 1.0, -6.7, 8.1 };
  MatrixD m2 (3, 1, threeNums);
  MatrixD m2a(1, 3, threeNums);

  Float8 fourNums[4] = { 9.1, 6.2, 5.2, -13 };
  MatrixD m3 (4, 1, fourNums);
  MatrixD m3a(1, 4, fourNums);

  Float8 nineNums[9] = {  0.5        ,  0.5        , sqrt(2.0)/2,
                         -0.5        , -0.5        , sqrt(2.0)/2,
                          sqrt(2.0)/2, -sqrt(2.0)/2, 0.0
                       };
  MatrixD m4(3, 3, nineNums);

  Float8 sixteenNums[16] = { 1, 0, 0, 0,
                             0, 1, 0, 0,
                             0, 0, 1, 0,
                             0, 0, 0, 1
                           };
  MatrixD m5(4, 4, sixteenNums);

  Float8 floatArray[4];


  PlaneVector p1 = pv( m1 );
  floatArray[0] = p1.x(); floatArray[1] = p1.y();
  if ( check( m1, (Float8 *)(&floatArray) ) ) {
    cout << "OK: PlaneVector p1 is same as Matrix m1 \n";
  } else {
    cout << "??? PlaneVector p1 isn't same as Matrix m1 \n";
    cout << p1 << endl << m1 << endl;
    problems = true;
  }

  PlaneVector p1a = pv( m1a );
  floatArray[0] = p1a.x(); floatArray[1] = p1a.y();
  if ( check( m1a, (Float8 *)(&floatArray) ) ) {
    cout << "OK: PlaneVector p1a is same as Matrix m1a \n";
  } else {
    cout << "??? PlaneVector p1a isn't same as Matrix m1a \n";
    cout << p1a << endl << m1a << endl;
    problems = true;
  }


  SpaceVector v1 = sv( m2 );
  floatArray[0] = v1.x(); floatArray[1] = v1.y();
  floatArray[2] = v1.z();
  if ( check( m2, (Float8 *)(&floatArray) ) ) {
    cout << "OK: SpaceVector v1 is same as Matrix m2 \n";
  } else {
    cout << "??? SpaceVector v1 isn't same as Matrix m2 \n";
    cout << v1 << endl << m2 << endl;
    problems = true;
  }


  SpaceVector v1a = sv( m2a );
  floatArray[0] = v1a.x(); floatArray[1] = v1a.y(); floatArray[2] = v1a.z();
  if ( check( m2a, (Float8 *)(&floatArray) ) ) {
    cout << "OK: SpaceVector v1a is same as Matrix m2a \n";
  } else {
    cout << "??? SpaceVector v1a isn't same as Matrix m2a \n";
    cout << v1a << endl << m2a << endl;
    problems = true;
  }


  LorentzVector w1 = lv( m3 );
  floatArray[0] = w1.x(); floatArray[1] = w1.y();
  floatArray[2] = w1.z(); floatArray[3] = w1.t();
  if ( check( m3, (Float8 *)(&floatArray) ) ) {
    cout << "OK: LorentzVector w1 is same as Matrix m3 \n";
  } else {
    cout << "??? LorentzVector w1 isn't same as Matrix m3 \n";
    cout << w1 << endl << m3 << endl;
    problems = true;
  }

  LorentzVector w1a = lv( m3a );
  floatArray[0] = w1a.x(); floatArray[1] = w1a.y();
  floatArray[2] = w1a.z(); floatArray[3] = w1a.t();
  if ( check( m3a, (Float8 *)(&floatArray) ) ) {
    cout << "OK: LorentzVector w1a is same as Matrix m3a \n";
  } else {
    cout << "??? LorentzVector w1a isn't same as Matrix m3a \n";
    cout << w1a << endl << m3a << endl;
    problems = true;
  }

  Rotation r1 = rot( m4 );
  ZMpvRep3x3 rep3 = r1.rep3x3();
  if ( check( m4, (Float8 *)(&rep3) ) ) {
    cout << "OK: Rotation r1 is same as Matrix m4 \n";
  } else {
    cout << "??? Rotation r1 isn't same as Matrix m4 \n";
    cout << r1 << endl << m4 << endl;
    problems = true;
  }

  LorentzTransformation lt1 = lt( m5 );
  ZMpvRep4x4 rep4 = lt1.rep4x4();
  if ( check( m5, (Float8 *)(&rep4) ) ) {
    cout << "OK: LorentzTransformation lt1 is same as Matrix m5 \n";
  } else {
    cout << "??? LorentzTransformation lt1 isn't same as Matrix m5 \n";
    cout << lt1 << endl << m5 << endl;
    problems = true;
  }

} // la2pv_constructors()


void la2pv_multiplication() {

    cout << "\n{**** Testing multiplication of matrices by Vectors ****}\n";

  Float8  twoNums[2] =  { 4.5, 0.1 };
  Float8 threNums[3] =  { 9.0, -6.2, 3.8 };
  Float8 fourNums[4] =  { 1.0, 2.1, 0.9, -5 };
  Float8 nineNums[9] =  {  0.5     ,  0.5     , 0.707107,
                          -0.5     , -0.5     , 0.707107,
                           0.707107, -0.707107, 0
                        };
  Float8 sxtnNums[16] = { 1, 0, 0, 0,
                          0, 1, 0, 0,
                          0, 0, 1, 0,
                          0, 0, 0, 1
                        };

  Float8 floatArray[4];


  MatrixD m2x1(2, 1, twoNums);
  PlaneVector p1( 0.5, 1.6 );

  MatrixD m1 = p1 * m2x1;
  floatArray[0] = p1.x(); floatArray[1] = p1.y();
  if ( check(m1, multiply((Float8 *)(&floatArray), m2x1, 1) ) ) {
    cout << "OK: Result matrix is correct for multiplication of a"
         << " PlaneVector by a matrix" << endl;
  } else {
    cout << "??? Result matrix incorrect for multiplication of a"
         << " PlaneVector by a matrix" << endl;
    cout << m1 << endl << m2x1 << endl << p1 << endl;
    problems = true;
  }



  MatrixD m3x1(3, 1, threNums);
  SpaceVector v1( 9.4, 1.6, 6.9 );

  MatrixD m2 = v1 * m3x1;
  floatArray[0] = v1.x(); floatArray[1] = v1.y();
  floatArray[2] = v1.z();
  if ( check(m2, multiply((Float8 *)(&floatArray), m3x1, 1) ) ) {
    cout << "OK: Result matrix is correct for multiplication of a"
         << " SpaceVector by a matrix" << endl;
  } else {
    cout << "??? Result matrix incorrect for multiplication of a"
         << " SpaceVector by a matrix" << endl;
    cout << m1 << endl << m3x1 << endl << v1 << endl;
    problems = true;
  }



  MatrixD m4x1(4, 1, fourNums);
  LorentzVector w1( 0.6, -.4, 1.8, Tcomponent(-3) );

  MatrixD m3 = w1 * m4x1;
  floatArray[0] = w1.x(); floatArray[1] = w1.y();
  floatArray[2] = w1.z(); floatArray[3] = w1.t();
  if ( check(m3, multiply((Float8 *)(&floatArray), m4x1, 1) ) ) {
    cout << "OK: Result matrix is correct for multiplication of a"
         << " LorentzVector by a matrix" << endl;
  } else {
    cout << "??? Result matrix incorrect for multiplication of a"
         << " LorentzVector by a matrix" << endl;
    cout << m3 << endl << m4x1 << endl << w1 << endl;
    problems = true;
  }



  MatrixD m3x3(3, 3, nineNums);
  Rotation r1( PI, -PI/6, PI/5);

  MatrixD m4 = r1 * m3x3;
  ZMpvRep3x3 rep3 = r1.rep3x3();
  if ( check(m4, multiply((Float8 *)(&rep3), m3x3, 3) ) ) {
    cout << "OK: Result matrix is correct for multiplication of a"
         << " Rotation by a matrix" << endl;
  } else {
    cout << "??? Result matrix incorrect for multiplication of a"
         << " Rotation by a matrix" << endl;
    cout << m4 << endl << m3x3 << endl << r1 << endl;
    problems = true;
  }



  MatrixD m4x4(4, 4, sxtnNums);
  LorentzBoost b1( .6, -.7, .1);
  LorentzTransformation lt1( b1, r1 );

  MatrixD m5 = lt1 * m4x4;
  ZMpvRep4x4 rep4 = lt1.rep4x4();
  if ( check(m5, multiply((Float8 *)(&rep4), m4x4, 4) ) ) {
    cout << "OK: Result matrix is correct for multiplication of a"
         << " LorentzTransformation by a matrix" << endl;
  } else {
    cout << "??? Result matrix incorrect for multiplication of a"
         << " LorentzTransformation by a matrix" << endl;
    cout << m5 << endl << m4x4 << endl << lt1 << endl;
    problems = true;
  }

} // la2pv_multiplication()


void specialtycases() {

  cout << "\n{**** Testing the special cases of Rotations and "
       << "LorentzTransformations **** }\n";

  RotationX rx1( .62 );
  RotationY ry1( -.5 );
  RotationZ rz1( 1.57);
  ZMpvRep3x3 rep3;

  MatrixD m1 = mat( rx1 );
  rep3 = rx1.rep3x3();
  if ( check(m1,(Float8 *)(&rep3)) ) {
    cout << "OK: Result matrix is same as RotationX" << endl;
  } else {
    cout << "??? Result matrix isn't same as RotationX" << endl;
    cout << m1 << endl << rx1 << endl;
    problems = true;
  }

  MatrixD m2 = mat( ry1 );
  rep3 = ry1.rep3x3();
  if ( check(m2,(Float8 *)(&rep3)) ) {
    cout << "OK: Result matrix is same as RotationY" << endl;
  } else {
    cout << "??? Result matrix isn't same as RotationY" << endl;
    cout << m2 << endl << ry1 << endl;
    problems = true;
  }

  MatrixD m3 = mat( rz1 );
  rep3 = rz1.rep3x3();
  if ( check(m3,(Float8 *)(&rep3)) ) {
    cout << "OK: Result matrix is same as RotationZ" << endl;
  } else {
    cout << "??? Result matrix isn't same as RotationZ" << endl;
    cout << m3 << endl << rz1 << endl;
    problems = true;
  }

  LorentzBoost   b1( .5, .72, .1 );
  LorentzBoostX bx1( .61 );
  LorentzBoostY by1( -.2 );
  LorentzBoostZ bz1( .99 );
  ZMpvRep4x4 rep4;

  MatrixD m4 = mat( b1 );
  rep4 = b1.rep4x4();
  if ( check(m4,(Float8 *)(&rep4)) ) {
    cout << "OK: Result matrix is same as LorentzBoost" << endl;
  } else {
    cout << "??? Result matrix isn't same as LorentzBoost" << endl;
    cout << m4 << endl << b1 << endl;
    problems = true;
  }

  MatrixD m5 = mat( bx1 );
  rep4 = bx1.rep4x4();
  if ( check(m5,(Float8 *)(&rep4)) ) {
    cout << "OK: Result matrix is same as LorentzBoostX" << endl;
  } else {
    cout << "??? Result matrix isn't same as LorentzBoostX" << endl;
    cout << m5 << endl << bx1 << endl;
    problems = true;
  }

  MatrixD m6 = mat( by1 );
  rep4 = by1.rep4x4();
  if ( check(m6,(Float8 *)(&rep4)) ) {
    cout << "OK: Result matrix is same as LorentzBoostY" << endl;
  } else {
    cout << "??? Result matrix isn't same as LorentzBoostY" << endl;
    cout << m6 << endl << by1 << endl;
    problems = true;
  }

  MatrixD m7 = mat( bz1 );
  rep4 = bz1.rep4x4();
  if ( check(m7,(Float8 *)(&rep4)) ) {
    cout << "OK: Result matrix is same as LorentzBoostZ" << endl;
  } else {
    cout << "??? Result matrix isn't same as LorentzBoostZ" << endl;
    cout << m7 << endl << bz1 << endl;
    problems = true;
  }


} // specialtycases()


void dotproducts() {

  cout << "\n{**** Testing equality of matrix products and dot products ****}"
       << endl;

  PlaneVector   p1( 1.6, 9.2 );
  SpaceVector   v1( 4.5, 1.2, 9.8 );
  LorentzVector w1( 5.2, 4.8, -6.1, Tcomponent(1.2) );

  MatrixD m1 = p1 * mat(p1);
  Float8 fl1 = p1.dot(p1);

  if ( check(m1,(Float8 *)(&fl1)) ) {
    cout << "OK: The matrix product of p1 yields the same as the dot product\n";
  } else {
    cout << "??? The matrix product of p1 is different than the dot product\n";
    cout << m1 << endl << fl1 << endl;
    problems = true;
  }


  MatrixD m2 = v1 * mat(v1);
  Float8 fl2 = v1.dot(v1);

  if ( check(m2,(Float8 *)(&fl2)) ) {
    cout << "OK: The matrix product of v1 yields the same as the dot product\n";
  } else {
    cout << "??? The matrix product of v1 is different than the dot product\n";
    cout << m2 << endl << fl2 << endl;
    problems = true;
  }

  MatrixD m3 = w1 * mat(w1);
  Float8 fl3 = w1.euclideanNorm2();

  if ( check(m3,(Float8 *)(&fl3)) ) {
    cout << "OK: The matrix product of w1 yields the same as the components "
         << "squared\n";
  } else {
    cout << "??? The matrix product of w1 is different than the components "
         << "squared\n";
    cout << m3 << endl << fl3 << endl;
    problems = true;
  }

} // dotproducts()


void sizetests() {

  cout << "\n{**** Testing to ensure proper exceptions are thrown ****}\n";

  Float8 twoNums[2] = { 3.0, 4.5 };
  MatrixD m1 (2, 1, twoNums);
  MatrixD m1a(1, 2, twoNums);

  Float8 threeNums[3] = { 1.0, -6.7, 8.1 };
  MatrixD m2 (3, 1, threeNums);
  MatrixD m2a(1, 3, threeNums);

  Float8 fourNums[4] = { 9.1, 6.2, 5.2, -13 };
  MatrixD m3 (4, 1, fourNums);
  MatrixD m3a(1, 4, fourNums);

  Float8 nineNums[9] = {  0.5     ,  0.5     , 0.707107,
                         -0.5     , -0.5     , 0.707107,
                          0.707107, -0.707107, 0
                       };
  MatrixD m4(3, 3, nineNums);

  Float8 sixteenNums[16] = { 1, 0, 0, 0,
                             0, 1, 0, 0,
                             0, 0, 1, 0,
                             0, 0, 0, 1
                           };
  MatrixD m5(4, 4, sixteenNums);

  cout << "Note: Expect only one exception per class: " << endl << endl;

  cout << "Attempting to incorrectly make a PlaneVector:" << endl;

  PlaneVector p1 = pv( m1 );  // Should be fine
  PlaneVector p2 = pv( m1a);  // Should be fine
  PlaneVector p3 = pv( m2 );  // nope
  cout << "Problem PlaneVector is now " << p3 << endl << endl;

  cout << "Attempting to incorrectly make a SpaceVector:" << endl;

  SpaceVector v1 = sv( m2 );  // Should be fine
  SpaceVector v2 = sv( m2a);  // Should be fine
  SpaceVector v3 = sv( m3 );  // nope
  cout << "Problem SpaceVector is now " << v3 << endl << endl;

  cout << "Attempting to incorrectly make a LorentzVector:" << endl;

  LorentzVector w1 = lv( m3 );  // Should be fine
  LorentzVector w2 = lv( m3a);  // Should be fine
  LorentzVector w3 = lv( m4 );  // nope
  cout << "Problem LorentzVector is now " << w3 << endl << endl;

  cout << "Attempting to incorrectly make a Rotation:" << endl;

  Rotation r1 = rot( m4 );  // Should be fine
  Rotation r2 = rot( m5 );  // nope
  cout << "Problem Rotation is now " << r2 << endl << endl;

  cout << "Attempting to incorrectly make a LorentzTransformation:" << endl;

  LorentzTransformation lt1 = lt( m5 ); // Should be fine
  LorentzTransformation lt2 = lt( m4 ); // nope
  cout << "Problem LorentzTransform is now: " << lt2 << endl << endl;

} // sizetests()

#endif // __GNUG__