newmat1.cxx

科学和工程计算中使用统计功能开发工具包

//$$ newmat1.cxx   Matrix type functions

// Copyright (C) 1991,2: R B Davies

//#define WANT_STREAM                     // to include stream package

#include "include.h"

#include "newmat.h"



/************************* MatrixType functions *****************************/


int MatrixType::BestFit() const
// given attributes find a matrix type that best fits them
{

   if (!(attribute & Valid)) return USX;

   if (attribute & (LUDeco + OneRow + OneCol))
   {
      if (attribute & OneCol) return CVX;
      if (attribute & OneRow) return RVX;
      if (attribute & LUDeco) return (attribute & Band) ? BCX : CtX;
   }

   if (attribute & Band)
   {
      if (attribute & Lower)
	 return (attribute & (Symmetric + Upper)) ? DgX : LBX;
      if (attribute & Upper)
	 return (attribute & (Symmetric + Lower)) ? DgX : UBX;
      if (attribute & Symmetric) return SBX;
      return BMX;
   }

   if (attribute & (Lower + Upper + Symmetric))
   {
      if (attribute & Lower)
	 return (attribute & (Symmetric + Upper)) ? DgX : LTX;
      if (attribute & Upper)
	 return (attribute & (Symmetric + Lower)) ? DgX : UTX;
      return SmX;
   }

   return RtX;   

}

Boolean MatrixType::operator>=(const MatrixType& t) const
{
   return ( attribute & (t.attribute | (OneRow + OneCol)) ) == attribute;
}

MatrixType MatrixType::operator+(const MatrixType& mt) const
{ return attribute & mt.attribute; }

MatrixType MatrixType::operator*(const MatrixType& mt) const
{
   if ( ((attribute & (Upper+Lower)) == (Upper+Lower))
      && ((mt.attribute & (Upper+Lower)) == (Upper+Lower)) )
         return Dg;
   else return (attribute & mt.attribute & ~Symmetric)
      | (attribute & OneRow) | (mt.attribute & OneCol);
}

MatrixType MatrixType::i() const
{ return attribute & ~(Band+LUDeco) ; }

MatrixType MatrixType::t() const
{
   int a = attribute & ~(Upper+Lower+OneRow+OneCol);
   if (attribute & Upper) a |= Lower;
   if (attribute & Lower) a |= Upper;
   if (attribute & OneRow) a |= OneCol;
   if (attribute & OneCol) a |= OneRow;
   return a;
}

MatrixType MatrixType::AddEqualEl() const
{ return attribute & (Valid + Symmetric + OneRow + OneCol); }

MatrixType MatrixType::MultRHS() const
{
   if (attribute & (Upper + Lower)) return attribute;
   else return attribute & ~Symmetric;
}

MatrixType MatrixType::sub() const
{ return attribute & (Valid + OneRow + OneCol); }


MatrixType MatrixType::ssub() const
{
   return *this;                  // won't work for selection matrix
}

MatrixType::operator char*() const
{
   switch (BestFit())
   {
   case USX:  return "UnSp ";
   case UTX:  return "UT   ";
   case LTX:  return "LT   ";
   case RtX:  return "Rect ";
   case SmX:  return "Sym  ";
   case DgX:  return "Diag ";
   case RVX:  return "RowV ";
   case CVX:  return "ColV ";
   case BMX:  return "Band ";
   case UBX:  return "UpBnd";
   case LBX:  return "LwBnd";
   case SBX:  return "SmBnd";
   case CtX:  return "Crout";
   case BCX:  return "BndLU";
   }
   return "?????";
}

GeneralMatrix* MatrixType::New() const
{
   GeneralMatrix* gm;
   switch (BestFit())
   {
   case USX:  Throw(CannotBuildException("UnSp"));
   case UTX:  gm = new UpperTriangularMatrix; break;
   case LTX:  gm = new LowerTriangularMatrix; break;
   case RtX:  gm = new Matrix; break;
   case SmX:  gm = new SymmetricMatrix; break;
   case DgX:  gm = new DiagonalMatrix; break;
   case RVX:  gm = new RowVector; break;
   case CVX:  gm = new ColumnVector; break;
   case BMX:  gm = new BandMatrix; break;
   case UBX:  gm = new UpperBandMatrix; break;
   case LBX:  gm = new LowerBandMatrix; break;
   case SBX:  gm = new SymmetricBandMatrix; break;
   case CtX:  Throw(CannotBuildException("Crout"));
   case BCX:  Throw(CannotBuildException("Band LU"));
   }
   MatrixErrorNoSpace(gm);
   gm->Protect(); return gm;
}

GeneralMatrix* MatrixType::New(int nr, int nc, BaseMatrix* bm) const
{
   Tracer tr("New");
   GeneralMatrix* gm;
   switch (BestFit())
   {
   case USX:  Throw(CannotBuildException("UnSp"));
   case UTX:  gm = new UpperTriangularMatrix(nr); break;
   case LTX:  gm = new LowerTriangularMatrix(nr); break;
   case RtX:  gm = new Matrix(nr, nc); break;
   case SmX:  gm = new SymmetricMatrix(nr); break;
   case DgX:  gm = new DiagonalMatrix(nr); break;
   case RVX:  if (nr!=1) Throw(VectorException());
	       gm = new RowVector(nc); break;
   case CVX:  if (nc!=1) Throw(VectorException());
	       gm = new ColumnVector(nr); break;
   case BMX:  {
                 MatrixBandWidth bw = bm->BandWidth();
                 gm = new BandMatrix(nr,bw.lower,bw.upper); break;
              }
   case UBX:  gm = new UpperBandMatrix(nr,bm->BandWidth().upper); break;
   case LBX:  gm = new LowerBandMatrix(nr,bm->BandWidth().lower); break;
   case SBX:  gm = new SymmetricBandMatrix(nr,bm->BandWidth().lower); break;
   case CtX:  Throw(CannotBuildException("Crout"));
   case BCX:  Throw(CannotBuildException("Band LU"));
   }
   MatrixErrorNoSpace(gm);
   gm->Protect(); return gm;
}