newmat4.cpp

非常好用的用C编写的矩阵类,可在不同编译器下编译使用.

// version with no conversion
void GeneralMatrix::Eq(const GeneralMatrix& X)
{
   GeneralMatrix* gmx = (GeneralMatrix*)&X;
   if (gmx!=this)
   {
      REPORT
      if (store)
      {
         MONITOR_REAL_DELETE("Free (operator=)",storage,store)
         REPORT delete [] store; storage = 0; store = 0;
      }
      GetMatrix(gmx);
   }
   else { REPORT }
   Protect();
}

// version to work with operator<<
void GeneralMatrix::Eq(const BaseMatrix& X, MatrixType mt, bool ldok)
{
   REPORT
   if (ldok) mt.SetDataLossOK();
   Eq(X, mt);
}

void GeneralMatrix::Eq2(const BaseMatrix& X, MatrixType mt)
// a cut down version of Eq for use with += etc.
// we know BaseMatrix points to two GeneralMatrix objects,
// the first being this (may be the same).
// we know tag_val has been set correctly in each.
{
   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(mt);
   if (gmx!=this) { REPORT GetMatrix(gmx); }  // simplify GetMatrix ?
   else { REPORT }
   Protect();
}

void GeneralMatrix::inject(const GeneralMatrix& X)
// copy stored values of X; otherwise leave els of *this unchanged
{
   REPORT
   Tracer tr("inject");
   if (nrows_val != X.nrows_val || ncols_val != X.ncols_val)
      Throw(IncompatibleDimensionsException());
   MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
   MatrixRow mrx(this, LoadOnEntry+StoreOnExit+DirectPart);
   int i=nrows_val;
   while (i--) { mrx.Inject(mr); mrx.Next(); mr.Next(); }
}

// ************* checking for data loss during conversion *******************/

bool Compare(const MatrixType& source, MatrixType& destination)
{
   if (!destination) { destination=source; return true; }
   if (destination==source) return true;
   if (!destination.DataLossOK && !(destination>=source))
      Throw(ProgramException("Illegal Conversion", source, destination));
   return false;
}

// ************* Make a copy of a matrix on the heap *********************/

GeneralMatrix* Matrix::Image() const
{
   REPORT
   GeneralMatrix* gm = new Matrix(*this); MatrixErrorNoSpace(gm);
   return gm;
}

GeneralMatrix* SquareMatrix::Image() const
{
   REPORT
   GeneralMatrix* gm = new SquareMatrix(*this); MatrixErrorNoSpace(gm);
   return gm;
}

GeneralMatrix* SymmetricMatrix::Image() const
{
   REPORT
   GeneralMatrix* gm = new SymmetricMatrix(*this); MatrixErrorNoSpace(gm);
   return gm;
}

GeneralMatrix* UpperTriangularMatrix::Image() const
{
   REPORT
   GeneralMatrix* gm = new UpperTriangularMatrix(*this);
   MatrixErrorNoSpace(gm); return gm;
}

GeneralMatrix* LowerTriangularMatrix::Image() const
{
   REPORT
   GeneralMatrix* gm = new LowerTriangularMatrix(*this);
   MatrixErrorNoSpace(gm); return gm;
}

GeneralMatrix* DiagonalMatrix::Image() const
{
   REPORT
   GeneralMatrix* gm = new DiagonalMatrix(*this); MatrixErrorNoSpace(gm);
   return gm;
}

GeneralMatrix* RowVector::Image() const
{
   REPORT
   GeneralMatrix* gm = new RowVector(*this); MatrixErrorNoSpace(gm);
   return gm;
}

GeneralMatrix* ColumnVector::Image() const
{
   REPORT
   GeneralMatrix* gm = new ColumnVector(*this); MatrixErrorNoSpace(gm);
   return gm;
}

GeneralMatrix* nricMatrix::Image() const
{
   REPORT
   GeneralMatrix* gm = new nricMatrix(*this); MatrixErrorNoSpace(gm);
   return gm;
}

GeneralMatrix* IdentityMatrix::Image() const
{
   REPORT
   GeneralMatrix* gm = new IdentityMatrix(*this); MatrixErrorNoSpace(gm);
   return gm;
}

GeneralMatrix* CroutMatrix::Image() const
{
   REPORT
   GeneralMatrix* gm = new CroutMatrix(*this); MatrixErrorNoSpace(gm);
   return gm;
}

GeneralMatrix* GeneralMatrix::Image() const
{
   bool dummy = true;
   if (dummy)                                   // get rid of warning message
      Throw(InternalException("Cannot apply Image to this matrix type"));
   return 0;
}


// *********************** nricMatrix routines *****************************/

void nricMatrix::MakeRowPointer()
{
   REPORT
   if (nrows_val > 0)
   {
      row_pointer = new Real* [nrows_val]; MatrixErrorNoSpace(row_pointer);
      Real* s = Store() - 1; int i = nrows_val; Real** rp = row_pointer;
      if (i) for (;;) { *rp++ = s; if (!(--i)) break; s+=ncols_val; }
   }
   else row_pointer = 0;
}

void nricMatrix::DeleteRowPointer()
   { REPORT if (nrows_val) delete [] row_pointer; }

void GeneralMatrix::CheckStore() const
{
   REPORT
   if (!store)
      Throw(ProgramException("NRIC accessing matrix with unset dimensions"));
}


// *************************** CleanUp routines *****************************/

void GeneralMatrix::cleanup()
{
   // set matrix dimensions to zero, delete storage
   REPORT
   if (store && storage)
   {
      MONITOR_REAL_DELETE("Free (cleanup)    ",storage,store)
      REPORT delete [] store;
   }
   store=0; storage=0; nrows_val=0; ncols_val=0; tag_val = -1;
}

void nricMatrix::cleanup()
   { REPORT DeleteRowPointer(); GeneralMatrix::cleanup(); }

void nricMatrix::MiniCleanUp()
   { REPORT DeleteRowPointer(); GeneralMatrix::MiniCleanUp(); }

void RowVector::cleanup()
   { REPORT GeneralMatrix::cleanup(); nrows_val=1; }

void ColumnVector::cleanup()
   { REPORT GeneralMatrix::cleanup(); ncols_val=1; }

void CroutMatrix::cleanup()
{
   REPORT
   if (nrows_val) delete [] indx;
   GeneralMatrix::cleanup();
}

void CroutMatrix::MiniCleanUp()
{
   REPORT
   if (nrows_val) delete [] indx;
   GeneralMatrix::MiniCleanUp();
}

void BandLUMatrix::cleanup()
{
   REPORT
   if (nrows_val) delete [] indx;
   if (storage2) delete [] store2;
   GeneralMatrix::cleanup();
}

void BandLUMatrix::MiniCleanUp()
{
   REPORT
   if (nrows_val) delete [] indx;
   if (storage2) delete [] store2;
   GeneralMatrix::MiniCleanUp();
}

// ************************ simple integer array class ***********************

// construct a new array of length xn. Check that xn is non-negative and
// that space is available

SimpleIntArray::SimpleIntArray(int xn) : n(xn)
{
   if (n < 0) Throw(Logic_error("invalid array length"));
   else if (n == 0) { REPORT  a = 0; }
   else { REPORT  a = new int [n]; if (!a) Throw(Bad_alloc()); }
}

// destroy an array - return its space to memory

SimpleIntArray::~SimpleIntArray() { REPORT  if (a) delete [] a; }

// access an element of an array; return a "reference" so elements
// can be modified.
// check index is within range
// in this array class the index runs from 0 to n-1

int& SimpleIntArray::operator[](int i)
{
   REPORT
   if (i < 0 || i >= n) Throw(Logic_error("array index out of range"));
   return a[i];
}

// same thing again but for arrays declared constant so we can't
// modify its elements

int SimpleIntArray::operator[](int i) const
{
   REPORT
   if (i < 0 || i >= n) Throw(Logic_error("array index out of range"));
   return a[i];
}

// set all the elements equal to a given value

void SimpleIntArray::operator=(int ai)
   { REPORT  for (int i = 0; i < n; i++) a[i] = ai; }

// set the elements equal to those of another array.
// now allow length to be changed

void SimpleIntArray::operator=(const SimpleIntArray& b)
{
   REPORT
   if (b.n != n) resize(b.n);
   for (int i = 0; i < n; i++) a[i] = b.a[i];
}

// construct a new array equal to an existing array
// check that space is available

SimpleIntArray::SimpleIntArray(const SimpleIntArray& b) : Janitor(), n(b.n)
{
   if (n == 0) { REPORT  a = 0; }
   else
   {
      REPORT
      a = new int [n]; if (!a) Throw(Bad_alloc());
      for (int i = 0; i < n; i++) a[i] = b.a[i];
   }
}

// change the size of an array; optionally copy data from old array to
// new array

void SimpleIntArray::resize(int n1, bool keep)
{
   if (n1 == n) { REPORT  return; }
   else if (n1 == 0) { REPORT  n = 0; delete [] a; a = 0; }
   else if (n == 0)
   {
      REPORT
      a = new int [n1]; if (!a) Throw(Bad_alloc());
      n = n1;
      if (keep) operator=(0);
   }
   else
   {
      int* a1 = a;
      if (keep)
      {
         REPORT
         int i;
         a = new int [n1]; if (!a) Throw(Bad_alloc());
         if (n > n1) n = n1;
         else for (i = n; i < n1; i++) a[i] = 0;
         for (i = 0; i < n; i++) a[i] = a1[i];
         n = n1; delete [] a1;
      }
      else
      {
         REPORT  n = n1; delete [] a1;
         a = new int [n]; if (!a) Throw(Bad_alloc());
      }
   }
}

//************************** swap values ********************************

// swap values

void GeneralMatrix::swap(GeneralMatrix& gm)
{
   REPORT
   int t;
   t = tag_val; tag_val = gm.tag_val; gm.tag_val = t;
   t = nrows_val; nrows_val = gm.nrows_val; gm.nrows_val = t;
   t = ncols_val; ncols_val = gm.ncols_val; gm.ncols_val = t;
   t = storage; storage = gm.storage; gm.storage = t;
   Real* s = store; store = gm.store; gm.store = s;
}
   
void nricMatrix::swap(nricMatrix& gm)
{
   REPORT
   GeneralMatrix::swap((GeneralMatrix&)gm);
   Real** rp = row_pointer; row_pointer = gm.row_pointer; gm.row_pointer = rp;
}

void CroutMatrix::swap(CroutMatrix& gm)
{
   REPORT
   GeneralMatrix::swap((GeneralMatrix&)gm);
   int* i = indx; indx = gm.indx; gm.indx = i;
   bool b;
   b = d; d = gm.d; gm.d = b;
   b = sing; sing = gm.sing; gm.sing = b;
}

void BandMatrix::swap(BandMatrix& gm)
{
   REPORT
   GeneralMatrix::swap((GeneralMatrix&)gm);
   int i;
   i = lower_val; lower_val = gm.lower_val; gm.lower_val = i;
   i = upper_val; upper_val = gm.upper_val; gm.upper_val = i;
}

void SymmetricBandMatrix::swap(SymmetricBandMatrix& gm)
{
   REPORT
   GeneralMatrix::swap((GeneralMatrix&)gm);
   int i;
   i = lower_val; lower_val = gm.lower_val; gm.lower_val = i;
}

void BandLUMatrix::swap(BandLUMatrix& gm)
{
   REPORT
   GeneralMatrix::swap((GeneralMatrix&)gm);
   int* i = indx; indx = gm.indx; gm.indx = i;
   bool b;
   b = d; d = gm.d; gm.d = b;
   b = sing; sing = gm.sing; gm.sing = b;
   int m;
   m = storage2; storage2 = gm.storage2; gm.storage2 = m;
   m = m1; m1 = gm.m1; gm.m1 = m;
   m = m2; m2 = gm.m2; gm.m2 = m;
   Real* s = store2; store2 = gm.store2; gm.store2 = s;
}

void GenericMatrix::swap(GenericMatrix& x)
{
   REPORT
   GeneralMatrix* tgm = gm; gm = x.gm; x.gm = tgm;
}

// ********************** C subscript classes ****************************

RealStarStar::RealStarStar(Matrix& A)
{
   REPORT
   Tracer tr("RealStarStar");
   int n = A.ncols();
   int m = A.nrows();
   a = new Real*[m];
   MatrixErrorNoSpace(a);
   Real* d = A.data();
   for (int i = 0; i < m; ++i) a[i] = d + i * n;
} 

ConstRealStarStar::ConstRealStarStar(const Matrix& A)
{
   REPORT
   Tracer tr("ConstRealStarStar");
   int n = A.ncols();
   int m = A.nrows();
   a = new const Real*[m];
   MatrixErrorNoSpace(a);
   const Real* d = A.data();
   for (int i = 0; i < m; ++i) a[i] = d + i * n;
} 



#ifdef use_namespace
}
#endif


/// \fn GeneralMatrix::SimpleAddOK(const GeneralMatrix* gm)
/// Can we add two matrices with simple vector add.
/// SimpleAddOK shows when we can add two matrices by a simple vector add
/// and when we can add one matrix into another
///
/// *gm must be the same type as *this
/// - return 0 if simple add is OK
/// - return 1 if we can add into *gm only
/// - return 2 if we can add into *this only
/// - return 3 if we can't add either way
///
/// Also applies to subtract;
/// for SP this will still be valid if we swap 1 and 2
///
/// For types Matrix, DiagonalMatrix, UpperTriangularMatrix,
/// LowerTriangularMatrix, SymmetricMatrix etc return 0.
/// For band matrices we will need to check bandwidths.







///@}