bandmat.cpp

C++矩阵算法库

//$$ bandmat.cpp                     Band matrix definitions

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

#define WANT_MATH                    // include.h will get math fns

//#define WANT_STREAM

#include "include.h"

#include "newmat.h"
#include "newmatrc.h"

#ifdef use_namespace
namespace NEWMAT {
#endif



#ifdef DO_REPORT
#define REPORT { static ExeCounter ExeCount(__LINE__,10); ++ExeCount; }
#else
#define REPORT {}
#endif

static inline int my_min(int x, int y) { return x < y ? x : y; }
static inline int my_max(int x, int y) { return x > y ? x : y; }


BandMatrix::BandMatrix(const BaseMatrix& M)
{
   REPORT // CheckConversion(M);
   // MatrixConversionCheck mcc;
   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::BM);
   GetMatrix(gmx); CornerClear();
}

void BandMatrix::SetParameters(const GeneralMatrix* gmx)
{
   REPORT
   MatrixBandWidth bw = gmx->BandWidth();
   lower = bw.lower; upper = bw.upper;
}

void BandMatrix::ReSize(int n, int lb, int ub)
{
   REPORT
   Tracer tr("BandMatrix::ReSize");
   if (lb<0 || ub<0) Throw(ProgramException("Undefined bandwidth"));
   lower = (lb<=n) ? lb : n-1; upper = (ub<=n) ? ub : n-1;
   GeneralMatrix::ReSize(n,n,n*(lower+1+upper)); CornerClear();
}

// SimpleAddOK shows when we can add etc 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
// For SP this will still be valid if we swap 1 and 2

short BandMatrix::SimpleAddOK(const GeneralMatrix* gm)
{
   const BandMatrix* bm = (const BandMatrix*)gm;
   if (bm->lower == lower && bm->upper == upper) { REPORT return 0; }
   else if (bm->lower >= lower && bm->upper >= upper) { REPORT return 1; }
   else if (bm->lower <= lower && bm->upper <= upper) { REPORT return 2; }
   else { REPORT return 3; }
}

short SymmetricBandMatrix::SimpleAddOK(const GeneralMatrix* gm)
{
   const SymmetricBandMatrix* bm = (const SymmetricBandMatrix*)gm;
   if (bm->lower == lower) { REPORT return 0; }
   else if (bm->lower > lower) { REPORT return 1; }
   else { REPORT return 2; }
}

void UpperBandMatrix::ReSize(int n, int lb, int ub)
{
   REPORT
   if (lb != 0)
   {
      Tracer tr("UpperBandMatrix::ReSize");
      Throw(ProgramException("UpperBandMatrix with non-zero lower band" ));
   }
   BandMatrix::ReSize(n, lb, ub);
}

void LowerBandMatrix::ReSize(int n, int lb, int ub)
{
   REPORT
   if (ub != 0)
   {
      Tracer tr("LowerBandMatrix::ReSize");
      Throw(ProgramException("LowerBandMatrix with non-zero upper band" ));
   }
   BandMatrix::ReSize(n, lb, ub);
}

void BandMatrix::ReSize(const GeneralMatrix& A)
{
   REPORT
   int n = A.Nrows();
   if (n != A.Ncols())
   {
      Tracer tr("BandMatrix::ReSize(GM)");
      Throw(NotSquareException(*this));
   }
   MatrixBandWidth mbw = A.BandWidth();
   ReSize(n, mbw.Lower(), mbw.Upper());
}

bool BandMatrix::SameStorageType(const GeneralMatrix& A) const
{
   if (Type() != A.Type()) { REPORT return false; }
   REPORT
   return BandWidth() == A.BandWidth();
}

void BandMatrix::ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B)
{
   REPORT
   Tracer tr("BandMatrix::ReSizeForAdd");
   MatrixBandWidth A_BW = A.BandWidth(); MatrixBandWidth B_BW = B.BandWidth();
   if ((A_BW.Lower() < 0) | (A_BW.Upper() < 0) | (B_BW.Lower() < 0)
      | (A_BW.Upper() < 0))
         Throw(ProgramException("Can't ReSize to BandMatrix" ));
   // already know A and B are square
   ReSize(A.Nrows(), my_max(A_BW.Lower(), B_BW.Lower()),
      my_max(A_BW.Upper(), B_BW.Upper()));
}

void BandMatrix::ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B)
{
   REPORT
   Tracer tr("BandMatrix::ReSizeForSP");
   MatrixBandWidth A_BW = A.BandWidth(); MatrixBandWidth B_BW = B.BandWidth();
   if ((A_BW.Lower() < 0) | (A_BW.Upper() < 0) | (B_BW.Lower() < 0)
      | (A_BW.Upper() < 0))
         Throw(ProgramException("Can't ReSize to BandMatrix" ));
   // already know A and B are square
   ReSize(A.Nrows(), my_min(A_BW.Lower(), B_BW.Lower()),
      my_min(A_BW.Upper(), B_BW.Upper()));
}


void BandMatrix::operator=(const BaseMatrix& X)
{
   REPORT // CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::BM); CornerClear();
}

void BandMatrix::CornerClear() const
{
   // set unused parts of BandMatrix to zero
   REPORT
   int i = lower; Real* s = store; int bw = lower + 1 + upper;
   while (i)
      { int j = i--; Real* sj = s; s += bw; while (j--) *sj++ = 0.0; }
   i = upper; s = store + storage;
   while (i)
      { int j = i--; Real* sj = s; s -= bw; while (j--) *(--sj) = 0.0; }
}

MatrixBandWidth MatrixBandWidth::operator+(const MatrixBandWidth& bw) const
{
   REPORT
   int l = bw.lower; int u = bw.upper;
   l = (lower < 0 || l < 0) ? -1 : (lower > l) ? lower : l;
   u = (upper < 0 || u < 0) ? -1 : (upper > u) ? upper : u;
   return MatrixBandWidth(l,u);
}

MatrixBandWidth MatrixBandWidth::operator*(const MatrixBandWidth& bw) const
{
   REPORT
   int l = bw.lower; int u = bw.upper;
   l = (lower < 0 || l < 0) ? -1 : lower+l;
   u = (upper < 0 || u < 0) ? -1 : upper+u;
   return MatrixBandWidth(l,u);
}

MatrixBandWidth MatrixBandWidth::minimum(const MatrixBandWidth& bw) const
{
   REPORT
   int l = bw.lower; int u = bw.upper;
   if ((lower >= 0) && ( (l < 0) || (l > lower) )) l = lower;
   if ((upper >= 0) && ( (u < 0) || (u > upper) )) u = upper;
   return MatrixBandWidth(l,u);
}

UpperBandMatrix::UpperBandMatrix(const BaseMatrix& M)
{
   REPORT // CheckConversion(M);
   // MatrixConversionCheck mcc;
   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UB);
   GetMatrix(gmx); CornerClear();
}

void UpperBandMatrix::operator=(const BaseMatrix& X)
{
   REPORT // CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::UB); CornerClear();
}

LowerBandMatrix::LowerBandMatrix(const BaseMatrix& M)
{
   REPORT // CheckConversion(M);
   // MatrixConversionCheck mcc;
   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LB);
   GetMatrix(gmx); CornerClear();
}

void LowerBandMatrix::operator=(const BaseMatrix& X)
{
   REPORT // CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::LB); CornerClear();
}

BandLUMatrix::BandLUMatrix(const BaseMatrix& m)
{
   REPORT
   Tracer tr("BandLUMatrix");
   storage2 = 0; store2 = 0;  // in event of exception during build
   GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate(MatrixType::BM);
   m1 = ((BandMatrix*)gm)->lower; m2 = ((BandMatrix*)gm)->upper;
   GetMatrix(gm);
   if (nrows!=ncols) Throw(NotSquareException(*this));
   d = true; sing = false;
   indx = new int [nrows]; MatrixErrorNoSpace(indx);
   MONITOR_INT_NEW("Index (BndLUMat)",nrows,indx)
   storage2 = nrows * m1;
   store2 = new Real [storage2]; MatrixErrorNoSpace(store2);
   MONITOR_REAL_NEW("Make (BandLUMat)",storage2,store2)
   ludcmp();
}

BandLUMatrix::~BandLUMatrix()
{
   REPORT
   MONITOR_INT_DELETE("Index (BndLUMat)",nrows,indx)
   MONITOR_REAL_DELETE("Delete (BndLUMt)",storage2,store2)
   delete [] indx; delete [] store2;
}

MatrixType BandLUMatrix::Type() const { REPORT return MatrixType::BC; }


LogAndSign BandLUMatrix::LogDeterminant() const
{
   REPORT
   if (sing) return 0.0;
   Real* a = store; int w = m1+1+m2; LogAndSign sum; int i = nrows;
   // while (i--) { sum *= *a; a += w; }
   if (i) for (;;) { sum *= *a; if (!(--i)) break; a += w; }
   if (!d) sum.ChangeSign(); return sum;
}

GeneralMatrix* BandMatrix::MakeSolver()
{
   REPORT
   GeneralMatrix* gm = new BandLUMatrix(*this);
   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
}


void BandLUMatrix::ludcmp()
{
   REPORT
   Real* a = store2; int i = storage2;
   // clear store2 - so unused locations are always zero -
   // required by operator==
   while (i--) *a++ = 0.0;
   a = store;
   i = m1; int j = m2; int k; int n = nrows; int w = m1 + 1 + m2;
   while (i)
   {
      Real* ai = a + i;
      k = ++j; while (k--) *a++ = *ai++;
      k = i--; while (k--) *a++ = 0.0;