mattest.cpp

用C++写的矩阵和矢量运算库

// $Id: mattest.cpp,v 1.32 2002/11/29 18:56:10 hkuiper Exp $

// CwMtx matrix and vector math library
// Copyright (C) 1999-2000  Harry Kuiper, Will DeVore(template conversion)

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include <iostream>
#include <stdlib.h>
#include <cwmtx.h>

// included for a simple test of namespaces
#include <list>

using namespace std;
using namespace CwMtx;

#ifndef M_PI
#define M_PI 3.1415926535897932384626433832795
#endif

// Test matrices with structured elements.
void test_struct_matrix()
{
  cout << endl << "Start test of matrices with structured elements" << endl;

  typedef CWTSMat<CWSquareMatrix, 2> CWSMat2;
  typedef CWTMat<CWTMatrix<CWSMat2>, 2, 3> CWMat23SMat2;
  typedef CWTMat<CWTMatrix<CWSMat2>, 3, 2> CWMat32SMat2;
  typedef CWTMat<CWTMatrix<CWSMat2>, 2, 2> CWMat22SMat2;

  cout << endl << "CWTMatrix<T>" << endl << endl;

  CWSMat2 smtx;
  cout << "smtx = CWTUnity<CWSMat2>() "
       << (smtx = CWTUnity<CWSMat2>()) << endl;
  cout << "smtx *= 3.0 "
       << (smtx *= 3.0) << endl;

  CWMat23SMat2 mtx23Smtx;
  mtx23Smtx.Fill(smtx);
  cout << "mtx23Smtx.Fill(smtx) "
       << mtx23Smtx << endl;

  CWMat32SMat2 mtx32Smtx;
  mtx32Smtx.Fill(smtx);
  cout << "mtx32Smtx.Fill(smtx) "
       << mtx32Smtx << endl;

  CWMat22SMat2 mtx22Smtx;
  cout << "mtx22Smtx = mtx23Smtx*mtx32Smtx "
       << (mtx22Smtx = mtx23Smtx*mtx32Smtx) << endl;
  cout << "mtx22Smtx = mtx22Smtx + mtx22Smtx "
       << (mtx22Smtx = mtx22Smtx + mtx22Smtx) << endl;
  cout << "mtx22Smtx += mtx22Smtx "
       << (mtx22Smtx += mtx22Smtx) << endl;
  cout << "mtx22Smtx = mtx22Smtx - mtx22Smtx*smtx "
       << (mtx22Smtx = mtx22Smtx - mtx22Smtx*smtx) << endl;
  cout << "mtx22Smtx -= mtx22Smtx*smtx "
       << (mtx22Smtx -= mtx22Smtx*smtx) << endl;
  cout << "mtx22Smtx = -mtx22Smtx "
       << (mtx22Smtx = -mtx22Smtx) << endl;
  cout << "mtx22Smtx = mtx22Smtx*smtx "
       << (mtx22Smtx = mtx22Smtx*smtx) << endl;
  cout << "mtx22Smtx = smtx*mtx22Smtx "
       << (mtx22Smtx = smtx*mtx22Smtx) << endl;
  cout << "mtx22Smtx *= smtx "
       << (mtx22Smtx *= smtx) << endl;
  cout << "mtx22Smtx = mtx22Smtx/smtx "
       << (mtx22Smtx = mtx22Smtx/smtx) << endl;
  cout << "mtx22Smtx /= smtx "
       << (mtx22Smtx /= smtx) << endl;

  typedef CWTSMat<CWTSquareMatrix<CWSMat2>, 2> CWSMat2SMat2;

  cout << endl << "CWTSquareMatrix<T>" << endl << endl;

  CWSMat2SMat2 smtxSmtx1, smtxSmtx2, smtxSmtx3;
  cout << "smtxSmtx1 = CWTUnity<CWSMat2SMat2>() "
       << (smtxSmtx1 = CWTUnity<CWSMat2SMat2>()) << endl;
  cout << "smtxSmtx2 = CWTUnity<CWSMat2SMat2>() "
       << (smtxSmtx2 = CWTUnity<CWSMat2SMat2>()) << endl;
  cout << "smtxSmtx1 = smtxSmtx1*smtx "
       << (smtxSmtx1 = smtxSmtx1*smtx) << endl;
  cout << "smtxSmtx1 = smtx*smtxSmtx1 "
       << (smtxSmtx1 = smtx*smtxSmtx1) << endl;
  cout << "smtxSmtx1 *= smtx "
       << (smtxSmtx1 *= smtx) << endl;
  cout << "smtxSmtx1 = smtxSmtx1/smtx "
       << (smtxSmtx1 = smtxSmtx1/smtx) << endl;
  cout << "smtxSmtx1 /= smtx "
       << (smtxSmtx1 /= smtx) << endl;

  cout << "smtxSmtx1 *= smtx "
       << (smtxSmtx1 *= smtx) << endl;

  cout << "smtxSmtx3 = smtxSmtx1 + smtxSmtx2 "
       << (smtxSmtx3 = smtxSmtx1 + smtxSmtx2) << endl;
  cout << "smtxSmtx3 += smtxSmtx1 "
       << (smtxSmtx3 += smtxSmtx1) << endl;
  cout << "smtxSmtx3 = smtxSmtx1 - smtxSmtx2 "
       << (smtxSmtx3 = smtxSmtx1 - smtxSmtx2) << endl;
  cout << "smtxSmtx3 -= smtxSmtx1 "
       << (smtxSmtx3 -= smtxSmtx1) << endl;
  cout << "smtxSmtx3 = -smtxSmtx1 "
       << (smtxSmtx3 = -smtxSmtx1) << endl;
  cout << "smtxSmtx3 = smtxSmtx1*smtxSmtx2 "
       << (smtxSmtx3 = smtxSmtx1*smtxSmtx2) << endl;
  cout << "smtxSmtx3 *= smtxSmtx1 "
       << (smtxSmtx3 *= smtxSmtx1) << endl;
  cout << "smtxSmtx3 = smtxSmtx2/smtxSmtx1 "
       << (smtxSmtx3 = smtxSmtx2/smtxSmtx1) << endl;
  cout << "smtxSmtx3 /= smtxSmtx1 "
       << (smtxSmtx3 /= smtxSmtx1) << endl;

  typedef CWTVec<CWTVector<CWSMat2>, 3> CWVec3SMat2;

  cout << endl << "CWTVector<T>" << endl << endl;

  CWSMat2 smtx2 = smtx;
  CWVec3SMat2 vecSmtx1, vecSmtx2, vecSmtx3;
  vecSmtx1.Fill(smtx2);
  cout << "vecSmtx1.Fill(smtx2) "
       << vecSmtx1 << endl;
  cout << "vecSmtx2 = vecSmtx1 "
       << (vecSmtx2 = vecSmtx1) << endl;
  cout << "vecSmtx3 = CWTZero<CWVec3SMat2>() "
       << (vecSmtx3 = CWTZero<CWVec3SMat2>()) << endl;

  cout << "vecSmtx3 = vecSmtx1 + vecSmtx2 "
       << (vecSmtx3 = vecSmtx1 + vecSmtx2) << endl;
  cout << "vecSmtx3 += vecSmtx1 "
       << (vecSmtx3 += vecSmtx1) << endl;
  cout << "vecSmtx3 = vecSmtx1 - vecSmtx2 "
       << (vecSmtx3 = vecSmtx1 - vecSmtx2) << endl;
  cout << "vecSmtx3 -= vecSmtx1 "
       << (vecSmtx3 -= vecSmtx1) << endl;
  cout << "vecSmtx3 = -vecSmtx1 "
       << (vecSmtx3 = -vecSmtx1) << endl;
  cout << "vecSmtx3 = vecSmtx3*smtx2 "
       << (vecSmtx3 = vecSmtx3*smtx2) << endl;
  cout << "vecSmtx3 = smtx2*vecSmtx3 "
       << (vecSmtx3 = smtx2*vecSmtx3) << endl;
  cout << "vecSmtx3 *= smtx2 "
       << (vecSmtx3 *= smtx2) << endl;
  cout << "smtx2 = vecSmtx1*vecSmtx2 "
       << (smtx2 = vecSmtx1*vecSmtx2) << endl;
  cout << "vecSmtx3 = vecSmtx3/smtx2 "
       << (vecSmtx3 = vecSmtx3/smtx2) << endl;
  cout << "vecSmtx3 /= smtx2 "
       << (vecSmtx3 /= smtx2) << endl;

  int f;
  cout << "f = vecSmtx1 != vecSmtx2 "
       << (f = vecSmtx1 != vecSmtx2) << endl;

  typedef CWTSpaceVector<CWSMat2> CWSvecSMat2;

  cout << endl << "CWTSpaceVector<T>" << endl << endl;

  CWSvecSMat2 svecSmtx1, svecSmtx2, svecSmtx3;
  cout << "svecSmtx1 = CWTZero<CWSvecSMat2>() "
       << (svecSmtx1 = CWTZero<CWSvecSMat2>()) << endl;
  svecSmtx2.Fill(smtx);
  cout << "svecSmtx2.Fill(smtx) "
       << svecSmtx2 << endl;

  cout << "svecSmtx3 = svecSmtx1 + svecSmtx2 "
       << (svecSmtx3 = svecSmtx1 + svecSmtx2) << endl;
  cout << "svecSmtx3 += svecSmtx2 "
       << (svecSmtx3 += svecSmtx2) << endl;
  cout << "svecSmtx3 = svecSmtx1 - svecSmtx2 "
       << (svecSmtx3 = svecSmtx1 - svecSmtx2) << endl;
  cout << "svecSmtx3 -= svecSmtx2 "
       << (svecSmtx3 -= svecSmtx2) << endl;
  cout << "svecSmtx3 = -svecSmtx2 "
       << (svecSmtx3 = -svecSmtx2) << endl;
  cout << "svecSmtx3 = svecSmtx2*smtx "
       << (svecSmtx3 = svecSmtx2*smtx) << endl;
  cout << "svecSmtx3 = smtx*svecSmtx2 "
       << (svecSmtx3 = smtx*svecSmtx2) << endl;
  cout << "svecSmtx3 *= smtx "
       << (svecSmtx3 *= smtx) << endl;
  cout << "smtx2 = svecSmtx2*svecSmtx2 "
       << (smtx2 = svecSmtx2*svecSmtx2) << endl;

  cout << "svecSmtx1[0] = smtx "
       << (svecSmtx1[0] = smtx) << endl;
  cout << "svecSmtx3 = svecSmtx1%svecSmtx2 "
       << (svecSmtx3 = svecSmtx1%svecSmtx2) << endl;
  cout << "svecSmtx3 %= svecSmtx2 "
       << (svecSmtx3 %= svecSmtx2) << endl;

  cout << "svecSmtx3 = svecSmtx2/smtx "
       << (svecSmtx3 = svecSmtx2/smtx) << endl;
  cout << "svecSmtx3 /= smtx "
       << (svecSmtx3 /= smtx) << endl;

  typedef CWTQuaternion<CWSMat2> CWQuatSMat2;

  cout << endl << "CWTQuaternion<T>" << endl << endl;

  CWQuatSMat2 qtnSmtx1, qtnSmtx2, qtnSmtx3;
  cout << "qtnSmtx1 = CWTZero<CWQuatSMat2>() "
       << (qtnSmtx1 = CWTZero<CWQuatSMat2>()) << endl;
  qtnSmtx1.Fill(CWTUnity<CWSMat2>());
  cout << "qtnSmtx1.Fill(CWTUnity<CWSMat2>()) "
       << qtnSmtx1  << endl;
  cout << "qtnSmtx2 = qtnSmtx1*smtx2 "
       << (qtnSmtx2 = qtnSmtx1*smtx2) << endl;
  cout << "qtnSmtx3 = qtnSmtx1 "
       << (qtnSmtx3 = qtnSmtx1) << endl;

  cout << "qtnSmtx2 = qtnSmtx1 + qtnSmtx1 "
       << (qtnSmtx2 = qtnSmtx1 + qtnSmtx1) << endl;
  cout << "qtnSmtx3 += qtnSmtx1 "
       << (qtnSmtx3 += qtnSmtx1) << endl;
  cout << "qtnSmtx2 = qtnSmtx1 - qtnSmtx1 "
       << (qtnSmtx2 = qtnSmtx1 - qtnSmtx1) << endl;
  cout << "qtnSmtx3 -= qtnSmtx1 "
       << (qtnSmtx3 -= qtnSmtx1) << endl;
  cout << "qtnSmtx3 = -qtnSmtx1 "
       << (qtnSmtx3 = -qtnSmtx1) << endl;
  cout << "qtnSmtx2 = smtx2*qtnSmtx1 "
       << (qtnSmtx2 = smtx2*qtnSmtx1) << endl;
  cout << "qtnSmtx3 *= qtnSmtx1 "
       << (qtnSmtx3 *= qtnSmtx1) << endl;
  cout << "qtnSmtx3 = qtnSmtx1*conj(qtnSmtx1) "
       << (qtnSmtx3 = qtnSmtx1*conj(qtnSmtx1)) << endl;
  cout << "qtnSmtx3 = inv(qtnSmtx1) "
       << (qtnSmtx3 = inv(qtnSmtx1)) << endl;
  cout << "qtnSmtx3 = qtnSmtx1/smtx2 "
       << (qtnSmtx3 = qtnSmtx1/smtx2) << endl;
  cout << "qtnSmtx3 /= smtx2 "
       << (qtnSmtx3 /= smtx2) << endl;
  cout << "qtnSmtx3 = qtnSmtx1/qtnSmtx2 "
       << (qtnSmtx3 = qtnSmtx1/qtnSmtx2) << endl;
  cout << "qtnSmtx3 /= qtnSmtx1 "
       << (qtnSmtx3 /= qtnSmtx1) << endl;

  cout << endl << "End test of matrices with structured elements" << endl;
}

//
// Test base matrix
//
void test_matrix ( void )
{

  cout << endl << "Start Matrix test" << endl;
  //_CWMatrix<> mx; must use namespace as below
  CwMtx::CWTMatrix<> mx;

  //
  // Test Phase 1. Check the mechcanics of template.
  //

  // The only way to test the template is through
  // instantiation. Otherwise the compiler does not include the code
  // and hence doesn't compile anything.  So the first set of tests in
  // Test Phase 1 is merely just compilation, the numbers are
  // ficticous.
  cout << endl << "CWMatrix constructors" << endl;
  CWMatrix M0;
  CWMatrix M1( 2 , 2 );
  CWMatrix M2( 3 , 3 );
  CWMatrix M3( M1 );
  CWMatrix M6 = M2;

  cout << endl << "Dimension CWMatrix" << endl;
  M0.Dimension( 4 , 5 );
  cout << "M0 = " << M0 << endl;
  CWMatrix M4( M0 , 1 , 2 , 3 , 4 );
  M4.Fill(1);
  cout << "M4 = " << M4 << endl;

  cout << endl << "MapInto" << endl;
  M1.MapInto( M0 , 0 , 0 , 1 , 1 );
  M1.Fill(1);
  cout << "M1 = " << M1 << endl;
  cout << "M0 = " << M0 << endl;
  cout << endl << "GetStatus" << endl;
  int s = M0.GetStatus();
  cout << "Status = " << s << endl;
	
  cout << endl << "GetCols and rows" << endl;
  M0.GetCols();
  M0.GetRows();
  double d = *M4[0];
  cout << "d = " << d << endl;
	
  cout << endl << "Matrix + Matrix" << endl;
  M3.Fill(1);
  cout << "M1 = " << M1 << endl;
  cout << "M3 = " << M3 << endl;
  CWMatrix M5 = M1 + M3;
  cout << "M5 = " << M5 << endl;
  cout << endl << "Matrix - Matrix" << endl;
  // this found an underscore "_" missing in the template
  M5 = M1 - M3;
  cout << "M5 = " << M5 << endl;
  cout << endl << "-Matrix" << endl;
  M5 = -M3;
  cout << "M5 = " << M5 << endl;

  cout << endl << "Matrix * scalar" << endl;
  M2.Fill(2);
  M5 = M2 * 2.0;
  cout << "M5 = " << M5 << endl;
  cout << endl << "scalar * Matrix" << endl;
  M5 = 2.0 * M2;
  cout << "M5 = " << M5 << endl;
  cout << endl << "Matrix * Matrix" << endl;
  M6.Fill(3);
  M5 = M2 * M6;
  cout << "M5 = " << M5 << endl;
  cout << endl << "Matrix / scalar" << endl;
  // this found the infinite recursion occuring in the "operator /()"
  M5 = M2 / 3.0;
  cout << "M5 = " << M5 << endl;

  M6 -= M2;
  M6 += M2;
  M6 *= 3.0;
  M6 /= 4.0;

  if ( M6 == M2 )
    cout << "M6 equal M2" << endl;
  else
    cout << "M6 not equal M2" << endl;

  if ( M6 != M2 )
    cout << "M6 not equal M2" << endl;
  else
    cout << "M6 equal M2" << endl;

  // Can't do any of the Store* methods until we allocate.
  CWMatrix M7( 3 , 3 );
  M7.StoreProduct( M6 , M2 );
  M7.StoreTranspose( M6 );
  // this core dumps if M7 was not defined is size.
  M7.StoreAtPosition( 0 , 0 , M2 );
  M7.Fill( 2 );
  M6.InterchangeRows( 1 , 2 );
  M2.AddRowToRow( 1 , 2 , 2.0 );
  M3.MultiplyRow( 0 , 3.0 );

  //
  // Test Phase 2. Check the mathimatical validty of the matrix template
  //

  // Some checks are going to be done with 2x2 matrices and other 3x3 matrices.
  CWMatrix M10( 2 , 2 );
  CWMatrix M11( 2 , 2 );
  M10[0][0] = 2.5;
  M10[0][1] = 0.0;
  M10[1][0] = 0.0;
  M10[1][1] = 2.5;

  M11[0][0] = 2.5;
  M11[0][1] = 0.0;
  M11[1][0] = 0.0;
  M11[1][1] = 2.5;
  cout << "M10 = " << M10 << endl;
  cout << "M11 = " << M11 << endl;
  cout << "M10*M11 = " << M10*M11 << endl;

  //
  // Test Phase 3
  //

  // Can't be tested until the other sibling classes are converted to templates

  cout << endl << "Deallocate" << endl;
  M0.Deallocate();
  cout << "M0 = " << M0 << endl;
  cout << endl << "GetStatus" << endl;
  s = M0.GetStatus();
  cout << "Status = " << s << endl;

  cout << endl << "End Matrix test" << endl;

}

// Test vector

void test_vector( )
{

  cout << endl << "Start Vector test" << endl;

  //
  // Test Phase 1. Check the mechcanics of template.
  //
  CWVector v0( 3 );
  v0[0] = 1.0;
  v0[1] = 1.0;
  v0[2] = 1.0;
  cout << "v0 = " << v0 << endl;
	
  CWMatrix m0( 3 , 1 );
  CWVector v1( m0 );
  CWVector v2( v1 );
  CWVector v3( m0 , 1 , 2 , 3 );