newmat8.cpp

各种矩阵算法库。支持UpperTriangularMatrix,LowerTriangularMatrix, DiagonalMatrix, SymmetricMatrix, BandMatrix,U

   REPORT
   int k; Real m = GeneralMatrix::MaximumAbsoluteValue1(k); k--;
   i = k / Ncols(); j = k - i * Ncols(); i++; j++;
   return m;
}

Real Matrix::MinimumAbsoluteValue2(int& i, int& j) const
{
   REPORT
   int k; Real m = GeneralMatrix::MinimumAbsoluteValue1(k); k--;
   i = k / Ncols(); j = k - i * Ncols(); i++; j++;
   return m;
}

Real Matrix::Maximum2(int& i, int& j) const
{
   REPORT
   int k; Real m = GeneralMatrix::Maximum1(k); k--;
   i = k / Ncols(); j = k - i * Ncols(); i++; j++;
   return m;
}

Real Matrix::Minimum2(int& i, int& j) const
{
   REPORT
   int k; Real m = GeneralMatrix::Minimum1(k); k--;
   i = k / Ncols(); j = k - i * Ncols(); i++; j++;
   return m;
}

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

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

Real SymmetricMatrix::Sum() const
{
   REPORT
   Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
   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 BaseMatrix::SumSquare() const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->SumSquare(); return s;
}

Real BaseMatrix::NormFrobenius() const
   { REPORT  return sqrt(SumSquare()); }

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

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

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

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

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

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

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

Real BaseMatrix::MinimumAbsoluteValue2(int& i, int& j) const
{
   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   Real s = gm->MinimumAbsoluteValue2(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) 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 trace("Trace");
   int i = nrows; int d = i+1;
   if (i != ncols) 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; Real sum = 0.0; Real* s = store;
   while (i--) sum += *s++;
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

Real SymmetricMatrix::Trace() const
{
   REPORT
   int i = nrows; 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; 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; 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; int w = lower+upper+1;
   Real sum = 0.0; Real* s = store+lower;
   // 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; int w = lower+1;
   Real sum = 0.0; Real* s = store+lower;
   // while (i--) { sum += *s; s += w; }
   if (i) for (;;) { sum += *s; if (!(--i)) break; s += w; }
   ((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_value += log(x); }
   else if (x < 0.0) { log_value += log(-x); sign = -sign; }
   else sign = 0;
}

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

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

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

LogAndSign LowerTriangularMatrix::LogDeterminant() const
{
   REPORT
   int i = nrows; 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::LogDeterminant() const
{
   REPORT
   int i = nrows; LogAndSign sum; Real* s = store;
   while (i) { sum *= *s; s += i--; }
   ((GeneralMatrix&)*this).tDelete(); return sum;
}

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

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

LogAndSign CroutMatrix::LogDeterminant() const
{
   REPORT
   if (sing) return 0.0;
   int i = nrows; 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->LogDeterminant();
   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(); }
}


#ifdef use_namespace
}
#endif