newmat2.cxx

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

//$$ newmat2.cxx      Matrix row and column operations

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

#define WANT_MATH

#include "include.h"

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

//#define REPORT { static ExeCounter ExeCount(__LINE__,2); ExeCount++; }

#define REPORT {}

//#define MONITOR(what,storage,store) \
//   { cout << what << " " << storage << " at " << (long)store << "\n"; }

#define MONITOR(what,store,storage) {}

/************************** Matrix Row/Col functions ************************/

void MatrixRowCol::Add(const MatrixRowCol& mrc)
{
   REPORT
   int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
   if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
   if (l<=0) return;
   Real* elx=store+f; Real* el=mrc.store+f;
   while (l--) *elx++ += *el++;
}

void MatrixRowCol::AddScaled(const MatrixRowCol& mrc, Real x)
{
   REPORT
   int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
   if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
   if (l<=0) return;
   Real* elx=store+f; Real* el=mrc.store+f;
   while (l--) *elx++ += *el++ * x;
}

void MatrixRowCol::Sub(const MatrixRowCol& mrc)
{
   REPORT
   int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
   if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
   if (l<=0) return;
   Real* elx=store+f; Real* el=mrc.store+f;
   while (l--) *elx++ -= *el++;
}

void MatrixRowCol::Inject(const MatrixRowCol& mrc)
// copy stored elements only
{
   REPORT
   int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
   if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
   if (l<=0) return;
   Real* elx = store+f; Real* ely = mrc.store+f;
   while (l--) *elx++ = *ely++;
}

Real DotProd(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
{
   REPORT                                         // not accessed
   int f = mrc1.skip; int f2 = mrc2.skip;
   int l = f + mrc1.storage; int l2 = f2 + mrc2.storage;
   if (f < f2) f = f2; if (l > l2) l = l2; l -= f;
   if (l<=0) return 0.0;
   
   Real* el1=mrc1.store+f; Real* el2=mrc2.store+f;
   Real sum = 0.0;
   while (l--) sum += *el1++ * *el2++;
   return sum;
}

void MatrixRowCol::Add(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
{
   int f = skip; int l = skip + storage;
   int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
   if (f1<f) f1=f; if (l1>l) l1=l;
   int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
   if (f2<f) f2=f; if (l2>l) l2=l;
   Real* el = store + f;
   Real* el1 = mrc1.store+f1; Real* el2 = mrc2.store+f2;
   if (f1<f2)
   {
      register int i = f1-f; while (i--) *el++ = 0.0;
      if (l1<=f2)                              // disjoint
      {
	 REPORT                                // not accessed
         i = l1-f1;     while (i--) *el++ = *el1++;
         i = f2-l1;     while (i--) *el++ = 0.0;
         i = l2-f2;     while (i--) *el++ = *el2++;
         i = l-l2;      while (i--) *el++ = 0.0;
      }
      else
      {
         i = f2-f1;    while (i--) *el++ = *el1++;
         if (l1<=l2)
         {
            REPORT
            i = l1-f2; while (i--) *el++ = *el1++ + *el2++;
            i = l2-l1; while (i--) *el++ = *el2++;
            i = l-l2;  while (i--) *el++ = 0.0;
         }
         else
         {
            REPORT
            i = l2-f2; while (i--) *el++ = *el1++ + *el2++;
            i = l1-l2; while (i--) *el++ = *el1++;
            i = l-l1;  while (i--) *el++ = 0.0;
         }
      }
   }
   else
   {
      register int i = f2-f; while (i--) *el++ = 0.0;
      if (l2<=f1)                              // disjoint
      {
	 REPORT                                // not accessed
         i = l2-f2;     while (i--) *el++ = *el2++;
         i = f1-l2;     while (i--) *el++ = 0.0;
	 i = l1-f1;     while (i--) *el++ = *el1++;
	 i = l-l1;      while (i--) *el++ = 0.0;
      }
      else
      {
         i = f1-f2;    while (i--) *el++ = *el2++;
         if (l2<=l1)
         {
	    REPORT
            i = l2-f1; while (i--) *el++ = *el1++ + *el2++;
            i = l1-l2; while (i--) *el++ = *el1++;
            i = l-l1;  while (i--) *el++ = 0.0;
         }
         else
         {
	    REPORT
            i = l1-f1; while (i--) *el++ = *el1++ + *el2++;
            i = l2-l1; while (i--) *el++ = *el2++;
            i = l-l2;  while (i--) *el++ = 0.0;
         }
      }
   }
}

void MatrixRowCol::Sub(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
{
   int f = skip; int l = skip + storage;
   int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
   if (f1<f) f1=f; if (l1>l) l1=l;
   int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
   if (f2<f) f2=f; if (l2>l) l2=l;
   Real* el = store + f;
   Real* el1 = mrc1.store+f1; Real* el2 = mrc2.store+f2;
   if (f1<f2)
   {
      register int i = f1-f; while (i--) *el++ = 0.0;
      if (l1<=f2)                              // disjoint
      {
	 REPORT                                // not accessed
         i = l1-f1;     while (i--) *el++ = *el1++;
         i = f2-l1;     while (i--) *el++ = 0.0;
         i = l2-f2;     while (i--) *el++ = - *el2++;
         i = l-l2;      while (i--) *el++ = 0.0;
      }
      else
      {
         i = f2-f1;    while (i--) *el++ = *el1++;
         if (l1<=l2)
         {
	    REPORT
            i = l1-f2; while (i--) *el++ = *el1++ - *el2++;
            i = l2-l1; while (i--) *el++ = - *el2++;
            i = l-l2;  while (i--) *el++ = 0.0;
         }
         else
         {
	    REPORT
            i = l2-f2; while (i--) *el++ = *el1++ - *el2++;
            i = l1-l2; while (i--) *el++ = *el1++;
            i = l-l1;  while (i--) *el++ = 0.0;
         }
      }
   }
   else
   {
      register int i = f2-f; while (i--) *el++ = 0.0;
      if (l2<=f1)                              // disjoint
      {
	 REPORT                                // not accessed
         i = l2-f2;     while (i--) *el++ = - *el2++;
         i = f1-l2;     while (i--) *el++ = 0.0;
         i = l1-f1;     while (i--) *el++ = *el1++;
         i = l-l1;      while (i--) *el++ = 0.0;
      }
      else
      {
         i = f1-f2;    while (i--) *el++ = - *el2++;
         if (l2<=l1)
         {
	    REPORT
            i = l2-f1; while (i--) *el++ = *el1++ - *el2++;
            i = l1-l2; while (i--) *el++ = *el1++;
            i = l-l1;  while (i--) *el++ = 0.0;
         }
         else
         {
            REPORT
            i = l1-f1; while (i--) *el++ = *el1++ - *el2++;
            i = l2-l1; while (i--) *el++ = - *el2++;
            i = l-l2;  while (i--) *el++ = 0.0;
         }
      }
   }
}


void MatrixRowCol::Add(const MatrixRowCol& mrc1, Real x)
{
   REPORT
   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
   if (f < skip) { f = skip; if (l < f) l = f; }
   if (l > lx) { l = lx; if (f > lx) f = lx; }

   Real* elx = store+skip; Real* ely = mrc1.store+f;

   int l1 = f-skip;  while (l1--) *elx++ = x;
       l1 = l-f;     while (l1--) *elx++ = *ely++ + x;
       lx -= l;      while (lx--) *elx++ = x;
}

void MatrixRowCol::RevSub(const MatrixRowCol& mrc1)
{
   REPORT
   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
   if (f < skip) { f = skip; if (l < f) l = f; }
   if (l > lx) { l = lx; if (f > lx) f = lx; }

   Real* elx = store+skip; Real* ely = mrc1.store+f;

   int l1 = f-skip;  while (l1--) { *elx = - *elx; elx++; }
       l1 = l-f;     while (l1--) { *elx = *ely++ - *elx; elx++; }
       lx -= l;      while (lx--) { *elx = - *elx; elx++; }
}

void MatrixRowCol::Copy(const MatrixRowCol& mrc1)
{
   REPORT
   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
   if (f < skip) { f = skip; if (l < f) l = f; }
   if (l > lx) { l = lx; if (f > lx) f = lx; }

   Real* elx = store+skip; Real* ely = mrc1.store+f;

   int l1 = f-skip;  while (l1--) *elx++ = 0.0;
       l1 = l-f;     while (l1--) *elx++ = *ely++;
       lx -= l;      while (lx--) *elx++ = 0.0;
}

void MatrixRowCol::Negate(const MatrixRowCol& mrc1)
{
   REPORT
   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
   if (f < skip) { f = skip; if (l < f) l = f; }
   if (l > lx) { l = lx; if (f > lx) f = lx; }

   Real* elx = store+skip; Real* ely = mrc1.store+f;

   int l1 = f-skip;  while (l1--) *elx++ = 0.0;
       l1 = l-f;     while (l1--) *elx++ = - *ely++;
       lx -= l;      while (lx--) *elx++ = 0.0;
}

void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, Real s)
{
   REPORT
   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
   if (f < skip) { f = skip; if (l < f) l = f; }
   if (l > lx) { l = lx; if (f > lx) f = lx; }

   Real* elx = store+skip; Real* ely = mrc1.store+f;

   int l1 = f-skip;  while (l1--) *elx++ = 0.0;
       l1 = l-f;     while (l1--) *elx++ = *ely++ * s;
       lx -= l;      while (lx--) *elx++ = 0.0;
}

void DiagonalMatrix::Solver(MatrixRowCol& mrc, const MatrixRowCol& mrc1)
{
   REPORT
   int f = mrc1.skip; int f0 = mrc.skip;
   int l = f + mrc1.storage; int lx = f0 + mrc.storage;
   if (f < f0) { f = f0; if (l < f) l = f; }
   if (l > lx) { l = lx; if (f > lx) f = lx; }

   Real* elx = mrc.store+f0; Real* eld = store+f;

   int l1 = f-f0;    while (l1--) *elx++ = 0.0;
       l1 = l-f;     while (l1--) *elx++ /= *eld++;
       lx -= l;      while (lx--) *elx++ = 0.0;
   // Solver makes sure input and output point to same memory
}

void MatrixRowCol::Copy(const Real*& r)
{
   REPORT
   Real* elx = store+skip; const Real* ely = r+skip; r += length;
   int l = storage; while (l--) *elx++ = *ely++;
}

void MatrixRowCol::Copy(Real r)
{
   REPORT
   Real* elx = store+skip; int l = storage; while (l--) *elx++ = r;
}

Real MatrixRowCol::SumAbsoluteValue()
{
   REPORT
   Real sum = 0.0; Real* elx = store+skip; int l = storage;
   while (l--) sum += fabs(*elx++);
   return sum;
}

void MatrixRowCol::SubRowCol(MatrixRowCol& mrc, int skip1, int l1) const
{
   mrc.length = l1;  mrc.store = store + skip1;  int d = skip - skip1;
   mrc.skip = (d < 0) ? 0 : d;  d = skip + storage - skip1;
   d = ((l1 < d) ? l1 : d) - mrc.skip;  mrc.storage = (d < 0) ? 0 : d;
   mrc.cw = 0;
}

MatrixRowCol::~MatrixRowCol()
{
   if (+(cw*IsACopy) && !(cw*StoreHere))
   {
      Real* f = store+skip;
      MONITOR_REAL_DELETE("Free    (RowCol)",storage,f) 
#ifdef Version21
      delete [] f;
#else
      delete [storage] f;
#endif
   }
}

MatrixRow::~MatrixRow() { if (+(cw*StoreOnExit)) gm->RestoreRow(*this); }

MatrixCol::~MatrixCol() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }