newmat7.cpp

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

      {
         REPORT
         // what if New throws and exception
         Try { gmx = mtd.New(nr,nc,this); }
         CatchAll
         {
            if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
            ReThrow;
         }
         SPDS(gmx,gm1,gm2);
         if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
         gmx->ReleaseAndDelete();
      }
   }
   return gmx;
}



//*************************** norm functions ****************************/

Real BaseMatrix::norm1() const
{
   // maximum of sum of absolute values of a column
   REPORT
   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   int nc = gm->Ncols(); Real value = 0.0;
   MatrixCol mc(gm, LoadOnEntry);
   while (nc--)
      { Real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
   gm->tDelete(); return value;
}

Real BaseMatrix::norm_infinity() const
{
   // maximum of sum of absolute values of a row
   REPORT
   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   int nr = gm->Nrows(); Real value = 0.0;
   MatrixRow mr(gm, LoadOnEntry);
   while (nr--)
      { Real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
   gm->tDelete(); return value;
}

//********************** Concatenation and stacking *************************/

GeneralMatrix* ConcatenatedMatrix::Evaluate(MatrixType mtx)
{
   REPORT
   Tracer tr("Concatenate");
      gm2 = ((BaseMatrix*&)bm2)->Evaluate();
      gm1 = ((BaseMatrix*&)bm1)->Evaluate();
      Compare(gm1->type() | gm2->type(),mtx);
      int nr=gm1->Nrows(); int nc = gm1->Ncols() + gm2->Ncols();
      if (nr != gm2->Nrows())
         Throw(IncompatibleDimensionsException(*gm1, *gm2));
      GeneralMatrix* gmx = mtx.New(nr,nc,this);
      MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
      MatrixRow mr(gmx, StoreOnExit+DirectPart);
      while (nr--) { mr.ConCat(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
}

GeneralMatrix* StackedMatrix::Evaluate(MatrixType mtx)
{
   REPORT
   Tracer tr("Stack");
      gm2 = ((BaseMatrix*&)bm2)->Evaluate();
      gm1 = ((BaseMatrix*&)bm1)->Evaluate();
      Compare(gm1->type() & gm2->type(),mtx);
      int nc=gm1->Ncols();
      int nr1 = gm1->Nrows(); int nr2 = gm2->Nrows();
      if (nc != gm2->Ncols())
         Throw(IncompatibleDimensionsException(*gm1, *gm2));
      GeneralMatrix* gmx = mtx.New(nr1+nr2,nc,this);
      MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
      MatrixRow mr(gmx, StoreOnExit+DirectPart);
      while (nr1--) { mr.Copy(mr1); mr1.Next(); mr.Next(); }
      while (nr2--) { mr.Copy(mr2); mr2.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
}

// ************************* equality of matrices ******************** //

static bool RealEqual(Real* s1, Real* s2, int n)
{
   int i = n >> 2;
   while (i--)
   {
      if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
      if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
   }
   i = n & 3; while (i--) if (*s1++ != *s2++) return false;
   return true;
}

static bool intEqual(int* s1, int* s2, int n)
{
   int i = n >> 2;
   while (i--)
   {
      if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
      if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
   }
   i = n & 3; while (i--) if (*s1++ != *s2++) return false;
   return true;
}


bool operator==(const BaseMatrix& A, const BaseMatrix& B)
{
   Tracer tr("BaseMatrix ==");
   REPORT
   GeneralMatrix* gmA = ((BaseMatrix&)A).Evaluate();
   GeneralMatrix* gmB = ((BaseMatrix&)B).Evaluate();

   if (gmA == gmB)                            // same matrix
      { REPORT gmA->tDelete(); return true; }

   if ( gmA->Nrows() != gmB->Nrows() || gmA->Ncols() != gmB->Ncols() )
                                              // different dimensions
      { REPORT gmA->tDelete(); gmB->tDelete(); return false; }

   // check for CroutMatrix or BandLUMatrix
   MatrixType AType = gmA->type(); MatrixType BType = gmB->type();
   if (AType.CannotConvert() || BType.CannotConvert() )
   {
      REPORT
      bool bx = gmA->IsEqual(*gmB);
      gmA->tDelete(); gmB->tDelete();
      return bx;
   }

   // is matrix storage the same
   // will need to modify if further matrix structures are introduced
   if (AType == BType && gmA->bandwidth() == gmB->bandwidth())
   {                                          // compare store
      REPORT
      bool bx = RealEqual(gmA->Store(),gmB->Store(),gmA->Storage());
      gmA->tDelete(); gmB->tDelete();
      return bx;
   }

   // matrix storage different - just subtract
   REPORT  return is_zero(*gmA-*gmB);
}

bool operator==(const GeneralMatrix& A, const GeneralMatrix& B)
{
   Tracer tr("GeneralMatrix ==");
   // May or may not call tDeletes
   REPORT

   if (&A == &B)                              // same matrix
      { REPORT return true; }

   if ( A.Nrows() != B.Nrows() || A.Ncols() != B.Ncols() )
      { REPORT return false; }                // different dimensions

   // check for CroutMatrix or BandLUMatrix
   MatrixType AType = A.Type(); MatrixType BType = B.Type();
   if (AType.CannotConvert() || BType.CannotConvert() )
      { REPORT  return A.IsEqual(B); }

   // is matrix storage the same
   // will need to modify if further matrix structures are introduced
   if (AType == BType && A.bandwidth() == B.bandwidth())
      { REPORT return RealEqual(A.Store(),B.Store(),A.Storage()); }

   // matrix storage different - just subtract
   REPORT  return is_zero(A-B);
}

bool GeneralMatrix::is_zero() const
{
   REPORT
   Real* s=store; int i = storage >> 2;
   while (i--)
   {
      if (*s++) return false; if (*s++) return false;
      if (*s++) return false; if (*s++) return false;
   }
   i = storage & 3; while (i--) if (*s++) return false;
   return true;
}

bool is_zero(const BaseMatrix& A)
{
   Tracer tr("BaseMatrix::is_zero");
   REPORT
   GeneralMatrix* gm1 = 0; bool bx;
   Try { gm1=((BaseMatrix&)A).Evaluate(); bx = gm1->is_zero(); }
   CatchAll { if (gm1) gm1->tDelete(); ReThrow; }
   gm1->tDelete();
   return bx;
}

// IsEqual functions - insist matrices are of same type
// as well as equal values to be equal

bool GeneralMatrix::IsEqual(const GeneralMatrix& A) const
{
   Tracer tr("GeneralMatrix IsEqual");
   if (A.type() != type())                       // not same types
      { REPORT return false; }
   if (&A == this)                               // same matrix
      { REPORT  return true; }
   if (A.nrows_val != nrows_val || A.ncols_val != ncols_val)
                                                 // different dimensions
   { REPORT return false; }
   // is matrix storage the same - compare store
   REPORT
   return RealEqual(A.store,store,storage);
}

bool CroutMatrix::IsEqual(const GeneralMatrix& A) const
{
   Tracer tr("CroutMatrix IsEqual");
   if (A.type() != type())                       // not same types
      { REPORT return false; }
   if (&A == this)                               // same matrix
      { REPORT  return true; }
   if (A.nrows_val != nrows_val || A.ncols_val != ncols_val)
                                                 // different dimensions
   { REPORT return false; }
   // is matrix storage the same - compare store
   REPORT
   return RealEqual(A.store,store,storage)
      && intEqual(((CroutMatrix&)A).indx, indx, nrows_val);
}


bool BandLUMatrix::IsEqual(const GeneralMatrix& A) const
{
   Tracer tr("BandLUMatrix IsEqual");
   if (A.type() != type())                       // not same types
      { REPORT  return false; }
   if (&A == this)                               // same matrix
      { REPORT  return true; }
   if ( A.Nrows() != nrows_val || A.Ncols() != ncols_val
      || ((BandLUMatrix&)A).m1 != m1 || ((BandLUMatrix&)A).m2 != m2 )
                                                 // different dimensions
   { REPORT  return false; }

   // matrix storage the same - compare store
   REPORT
   return RealEqual(A.Store(),store,storage)
      && RealEqual(((BandLUMatrix&)A).store2,store2,storage2)
      && intEqual(((BandLUMatrix&)A).indx, indx, nrows_val);
}


// ************************* cross products ******************** //

inline void crossproduct_body(Real* a, Real* b, Real* c)
{
   c[0] = a[1] * b[2] - a[2] * b[1];
   c[1] = a[2] * b[0] - a[0] * b[2];
   c[2] = a[0] * b[1] - a[1] * b[0];
}

Matrix crossproduct(const Matrix& A, const Matrix& B)
{
   REPORT
   int ac = A.Ncols(); int ar = A.Nrows();
   int bc = B.Ncols(); int br = B.Nrows();
   Real* a = A.Store(); Real* b = B.Store();
   if (ac == 3)
   {
      if (bc != 3 || ar != 1 || br != 1)
         { Tracer et("crossproduct"); IncompatibleDimensionsException(A, B); }
      REPORT
      RowVector C(3);  Real* c = C.Store(); crossproduct_body(a, b, c);
      return (Matrix&)C;
   }
   else
   {
      if (ac != 1 || bc != 1 || ar != 3 || br != 3)
         { Tracer et("crossproduct"); IncompatibleDimensionsException(A, B); }
      REPORT
      ColumnVector C(3);  Real* c = C.Store(); crossproduct_body(a, b, c);
      return (Matrix&)C;
   }
}

ReturnMatrix crossproduct_rows(const Matrix& A, const Matrix& B)
{
   REPORT
   int n = A.Nrows();
   if (A.Ncols() != 3 || B.Ncols() != 3 || n != B.Nrows())
   {
      Tracer et("crossproduct_rows"); IncompatibleDimensionsException(A, B);
   }
   Matrix C(n, 3);
   Real* a = A.Store(); Real* b = B.Store(); Real* c = C.Store();
   if (n--)
   {
      for (;;)
      {
         crossproduct_body(a, b, c);
         if (!(n--)) break;
         a += 3; b += 3; c += 3;
      }
   }

   return C.ForReturn();
}

ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B)
{
   REPORT
   int n = A.Ncols();
   if (A.Nrows() != 3 || B.Nrows() != 3 || n != B.Ncols())
   {
      Tracer et("crossproduct_columns");
      IncompatibleDimensionsException(A, B);
   }
   Matrix C(3, n);
   Real* a = A.Store(); Real* b = B.Store(); Real* c = C.Store();
   Real* an = a+n; Real* an2 = an+n;
   Real* bn = b+n; Real* bn2 = bn+n;
   Real* cn = c+n; Real* cn2 = cn+n;

   int i = n; 
   while (i--)
   {
      *c++   = *an    * *bn2   - *an2   * *bn;
      *cn++  = *an2++ * *b     - *a     * *bn2++;
      *cn2++ = *a++   * *bn++  - *an++  * *b++;
   }

   return C.ForReturn();
}


#ifdef use_namespace
}
#endif

///@}