newmat8.cpp

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

{
   REPORT
   Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows_val;
   for (int i = 0; i<nr; i++)
   {
      int j = i;
      while (j--) sum2 += *s++;
      sum1 += *s++;
   }
   ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
}

Real IdentityMatrix::sum_square() const
{
   Real sum = *store * *store * nrows_val;
   ((GeneralMatrix&)*this).tDelete(); return sum;
}


Real BaseMatrix::sum_square() const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->sum_square(); return s;
}

Real BaseMatrix::norm_Frobenius() const
   { REPORT  return sqrt(sum_square()); }

Real BaseMatrix::sum_absolute_value() const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->sum_absolute_value(); return s;
}

Real BaseMatrix::sum() const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->sum(); return s;
}

Real BaseMatrix::maximum_absolute_value() const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->maximum_absolute_value(); return s;
}

Real BaseMatrix::maximum_absolute_value1(int& i) const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->maximum_absolute_value1(i); return s;
}

Real BaseMatrix::maximum_absolute_value2(int& i, int& j) const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->maximum_absolute_value2(i, j); return s;
}

Real BaseMatrix::minimum_absolute_value() const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->minimum_absolute_value(); return s;
}

Real BaseMatrix::minimum_absolute_value1(int& i) const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->minimum_absolute_value1(i); return s;
}

Real BaseMatrix::minimum_absolute_value2(int& i, int& j) const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->minimum_absolute_value2(i, j); return s;
}

Real BaseMatrix::maximum() const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->maximum(); return s;
}

Real BaseMatrix::maximum1(int& i) const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->maximum1(i); return s;
}

Real BaseMatrix::maximum2(int& i, int& j) const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->maximum2(i, j); return s;
}

Real BaseMatrix::minimum() const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->minimum(); return s;
}

Real BaseMatrix::minimum1(int& i) const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->minimum1(i); return s;
}

Real BaseMatrix::minimum2(int& i, int& j) const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->minimum2(i, j); return s;
}

Real dotproduct(const Matrix& A, const Matrix& B)
{
   REPORT
   int n = A.storage;
   if (n != B.storage)
   {
      Tracer tr("dotproduct");
      Throw(IncompatibleDimensionsException(A,B));
   }
   Real sum = 0.0; Real* a = A.store; Real* b = B.store;
   while (n--) sum += *a++ * *b++;
   return sum;
}

Real Matrix::trace() const
{
   REPORT
   Tracer tr("trace");
   int i = nrows_val; int d = i+1;
   if (i != ncols_val) Throw(NotSquareException(*this));
   Real sum = 0.0; Real* s = store;
//   while (i--) { sum += *s; s += d; }
   if (i) for (;;) { sum += *s; if (!(--i)) break; s += d; }
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

Real DiagonalMatrix::trace() const
{
   REPORT
   int i = nrows_val; Real sum = 0.0; Real* s = store;
   while (i--) sum += *s++;
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

Real SymmetricMatrix::trace() const
{
   REPORT
   int i = nrows_val; Real sum = 0.0; Real* s = store; int j = 2;
   // while (i--) { sum += *s; s += j++; }
   if (i) for (;;) { sum += *s; if (!(--i)) break; s += j++; }
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

Real LowerTriangularMatrix::trace() const
{
   REPORT
   int i = nrows_val; Real sum = 0.0; Real* s = store; int j = 2;
   // while (i--) { sum += *s; s += j++; }
   if (i) for (;;) { sum += *s; if (!(--i)) break; s += j++; }
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

Real UpperTriangularMatrix::trace() const
{
   REPORT
   int i = nrows_val; Real sum = 0.0; Real* s = store;
   while (i) { sum += *s; s += i--; }             // won t cause a problem
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

Real BandMatrix::trace() const
{
   REPORT
   int i = nrows_val; int w = lower_val+upper_val+1;
   Real sum = 0.0; Real* s = store+lower_val;
   // while (i--) { sum += *s; s += w; }
   if (i) for (;;) { sum += *s; if (!(--i)) break; s += w; }
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

Real SymmetricBandMatrix::trace() const
{
   REPORT
   int i = nrows_val; int w = lower_val+1;
   Real sum = 0.0; Real* s = store+lower_val;
   // while (i--) { sum += *s; s += w; }
   if (i) for (;;) { sum += *s; if (!(--i)) break; s += w; }
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

Real IdentityMatrix::trace() const
{
   Real sum = *store * nrows_val;
   ((GeneralMatrix&)*this).tDelete(); return sum;
}


Real BaseMatrix::trace() const
{
   REPORT
   MatrixType Diag = MatrixType::Dg; Diag.SetDataLossOK();
   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(Diag);
   Real sum = gm->trace(); return sum;
}

void LogAndSign::operator*=(Real x)
{
   if (x > 0.0) { log_val += log(x); }
   else if (x < 0.0) { log_val += log(-x); sign_val = -sign_val; }
   else sign_val = 0;
}

void LogAndSign::pow_eq(int k)
{
   if (sign_val)
   {
      log_val *= k;
      if ( (k & 1) == 0 ) sign_val = 1;
   }
}

Real LogAndSign::value() const
{
   Tracer et("LogAndSign::value");
   if (log_val >= FloatingPointPrecision::LnMaximum())
      Throw(OverflowException("Overflow in exponential"));
   return sign_val * exp(log_val);
}

LogAndSign::LogAndSign(Real f)
{
   if (f == 0.0) { log_val = 0.0; sign_val = 0; return; }
   else if (f < 0.0) { sign_val = -1; f = -f; }
   else sign_val = 1;
   log_val = log(f);
}

LogAndSign DiagonalMatrix::log_determinant() const
{
   REPORT
   int i = nrows_val; LogAndSign sum; Real* s = store;
   while (i--) sum *= *s++;
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

LogAndSign LowerTriangularMatrix::log_determinant() const
{
   REPORT
   int i = nrows_val; LogAndSign sum; Real* s = store; int j = 2;
   // while (i--) { sum *= *s; s += j++; }
   if (i) for(;;) { sum *= *s; if (!(--i)) break; s += j++; }
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

LogAndSign UpperTriangularMatrix::log_determinant() const
{
   REPORT
   int i = nrows_val; LogAndSign sum; Real* s = store;
   while (i) { sum *= *s; s += i--; }
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

LogAndSign IdentityMatrix::log_determinant() const
{
   REPORT
   int i = nrows_val; LogAndSign sum;
   if (i > 0) { sum = *store; sum.PowEq(i); }
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

LogAndSign BaseMatrix::log_determinant() const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   LogAndSign sum = gm->log_determinant(); return sum;
}

LogAndSign GeneralMatrix::log_determinant() const
{
   REPORT
   Tracer tr("log_determinant");
   if (nrows_val != ncols_val) Throw(NotSquareException(*this));
   CroutMatrix C(*this); return C.log_determinant();
}

LogAndSign CroutMatrix::log_determinant() const
{
   REPORT
   if (sing) return 0.0;
   int i = nrows_val; int dd = i+1; LogAndSign sum; Real* s = store;
   if (i) for(;;)
   {
      sum *= *s;
      if (!(--i)) break;
      s += dd;
   }
   if (!d) sum.ChangeSign(); return sum;

}

Real BaseMatrix::determinant() const
{
   REPORT
   Tracer tr("determinant");
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   LogAndSign ld = gm->log_determinant();
   return ld.Value();
}

LinearEquationSolver::LinearEquationSolver(const BaseMatrix& bm)
{
   gm = ( ((BaseMatrix&)bm).Evaluate() )->MakeSolver();
   if (gm==&bm) { REPORT  gm = gm->Image(); }
   // want a copy if  *gm is actually bm
   else { REPORT  gm->Protect(); }
}

ReturnMatrix BaseMatrix::sum_square_rows() const
{
   REPORT
   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   int nr = gm->nrows();
   ColumnVector ssq(nr);
   if (gm->size() == 0) { REPORT ssq = 0.0; }
   else
   {
      MatrixRow mr(gm, LoadOnEntry);
      for (int i = 1; i <= nr; ++i)
      {
         Real sum = 0.0;
         int s = mr.Storage();
         Real* in = mr.Data();
         while (s--) sum += square(*in++);
         ssq(i) = sum;   
         mr.Next();
      }
   }
   gm->tDelete();
   ssq.release(); return ssq.for_return();
}

ReturnMatrix BaseMatrix::sum_square_columns() const
{
   REPORT
   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   int nr = gm->nrows(); int nc = gm->ncols();
   RowVector ssq(nc); ssq = 0.0;
   if (gm->size() != 0)
   {
      MatrixRow mr(gm, LoadOnEntry);
      for (int i = 1; i <= nr; ++i)
      {
         int s = mr.Storage();
         Real* in = mr.Data(); Real* out = ssq.data() + mr.Skip();
         while (s--) *out++ += square(*in++);
         mr.Next();
      }
   }
   gm->tDelete();
   ssq.release(); return ssq.for_return();
}

ReturnMatrix BaseMatrix::sum_rows() const
{
   REPORT
   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   int nr = gm->nrows();
   ColumnVector sum_vec(nr);
   if (gm->size() == 0) { REPORT sum_vec = 0.0; }
   else
   {
      MatrixRow mr(gm, LoadOnEntry);
      for (int i = 1; i <= nr; ++i)
      {
         Real sum = 0.0;
         int s = mr.Storage();
         Real* in = mr.Data();
         while (s--) sum += *in++;
         sum_vec(i) = sum;   
         mr.Next();
      }
   }
   gm->tDelete();
   sum_vec.release(); return sum_vec.for_return();
}

ReturnMatrix BaseMatrix::sum_columns() const
{
   REPORT
   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   int nr = gm->nrows(); int nc = gm->ncols();
   RowVector sum_vec(nc); sum_vec = 0.0;
   if (gm->size() != 0)
   {
      MatrixRow mr(gm, LoadOnEntry);
      for (int i = 1; i <= nr; ++i)
      {
         int s = mr.Storage();
         Real* in = mr.Data(); Real* out = sum_vec.data() + mr.Skip();
         while (s--) *out++ += *in++;
         mr.Next();
      }
   }
   gm->tDelete();
   sum_vec.release(); return sum_vec.for_return();
}


#ifdef use_namespace
}
#endif


///}