newmat8.cpp

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

/// \ingroup newmat
///@{

/// \file newmat8.cpp
/// LU transform, scalar functions of matrices.

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

#define WANT_MATH

#include "include.h"

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

#ifdef use_namespace
namespace NEWMAT {
#endif


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


/************************** LU transformation ****************************/

void CroutMatrix::ludcmp()
// LU decomposition from Golub & Van Loan, algorithm 3.4.1, (the "outer
// product" version).
// This replaces the code derived from Numerical Recipes in C in previous
// versions of newmat and being row oriented runs much faster with large
// matrices.
{
   REPORT
   Tracer tr( "Crout(ludcmp)" ); sing = false;
   Real* akk = store;                    // runs down diagonal

   Real big = fabs(*akk); int mu = 0; Real* ai = akk; int k;

   for (k = 1; k < nrows_val; k++)
   {
      ai += nrows_val; const Real trybig = fabs(*ai);
      if (big < trybig) { big = trybig; mu = k; }
   }


   if (nrows_val) for (k = 0;;)
   {
      /*
      int mu1;
      {
         Real big = fabs(*akk); mu1 = k; Real* ai = akk; int i;

         for (i = k+1; i < nrows_val; i++)
         {
            ai += nrows_val; const Real trybig = fabs(*ai);
            if (big < trybig) { big = trybig; mu1 = i; }
         }
      }
      if (mu1 != mu) cout << k << " " << mu << " " << mu1 << endl;
      */

      indx[k] = mu;

      if (mu != k)                       //row swap
      {
         Real* a1 = store + nrows_val * k;
         Real* a2 = store + nrows_val * mu; d = !d;
         int j = nrows_val;
         while (j--) { const Real temp = *a1; *a1++ = *a2; *a2++ = temp; }
      }

      Real diag = *akk; big = 0; mu = k + 1;
      if (diag != 0)
      {
         ai = akk; int i = nrows_val - k - 1;
         while (i--)
         {
            ai += nrows_val; Real* al = ai;
            Real mult = *al / diag; *al = mult;
            int l = nrows_val - k - 1; Real* aj = akk;
            // work out the next pivot as part of this loop
            // this saves a column operation
            if (l-- != 0)
            {
               *(++al) -= (mult * *(++aj));
               const Real trybig = fabs(*al);
               if (big < trybig) { big = trybig; mu = nrows_val - i - 1; }
               while (l--) *(++al) -= (mult * *(++aj));
            }
         }
      }
      else sing = true;
      if (++k == nrows_val) break;          // so next line won't overflow
      akk += nrows_val + 1;
   }
}

void CroutMatrix::lubksb(Real* B, int mini)
{
   REPORT
   // this has been adapted from Numerical Recipes in C. The code has been
   // substantially streamlined, so I do not think much of the original
   // copyright remains. However there is not much opportunity for
   // variation in the code, so it is still similar to the NR code.
   // I follow the NR code in skipping over initial zeros in the B vector.

   Tracer tr("Crout(lubksb)");
   if (sing) Throw(SingularException(*this));
   int i, j, ii = nrows_val;       // ii initialised : B might be all zeros


   // scan for first non-zero in B
   for (i = 0; i < nrows_val; i++)
   {
      int ip = indx[i]; Real temp = B[ip]; B[ip] = B[i]; B[i] = temp;
      if (temp != 0.0) { ii = i; break; }
   }

   Real* bi; Real* ai;
   i = ii + 1;

   if (i < nrows_val)
   {
      bi = B + ii; ai = store + ii + i * nrows_val;
      for (;;)
      {
         int ip = indx[i]; Real sum = B[ip]; B[ip] = B[i];
         Real* aij = ai; Real* bj = bi; j = i - ii;
         while (j--) sum -= *aij++ * *bj++;
         B[i] = sum;
         if (++i == nrows_val) break;
         ai += nrows_val;
      }
   }

   ai = store + nrows_val * nrows_val;

   for (i = nrows_val - 1; i >= mini; i--)
   {
      Real* bj = B+i; ai -= nrows_val; Real* ajx = ai+i;
      Real sum = *bj; Real diag = *ajx;
      j = nrows_val - i; while(--j) sum -= *(++ajx) * *(++bj);
      B[i] = sum / diag;
   }
}

/****************************** scalar functions ****************************/

inline Real square(Real x) { return x*x; }

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

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

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

// maxima and minima

// There are three sets of routines
// maximum_absolute_value, minimum_absolute_value, maximum, minimum
// ... these find just the maxima and minima
// maximum_absolute_value1, minimum_absolute_value1, maximum1, minimum1
// ... these find the maxima and minima and their locations in a
//     one dimensional object
// maximum_absolute_value2, minimum_absolute_value2, maximum2, minimum2
// ... these find the maxima and minima and their locations in a
//     two dimensional object

// If the matrix has no values throw an exception

// If we do not want the location find the maximum or minimum on the
// array stored by GeneralMatrix
// This won't work for BandMatrices. We call ClearCorner for
// maximum_absolute_value but for the others use the absolute_minimum_value2
// version and discard the location.

// For one dimensional objects, when we want the location of the
// maximum or minimum, work with the array stored by GeneralMatrix

// For two dimensional objects where we want the location of the maximum or
// minimum proceed as follows:

// For rectangular matrices use the array stored by GeneralMatrix and
// deduce the location from the location in the GeneralMatrix

// For other two dimensional matrices use the Matrix Row routine to find the
// maximum or minimum for each row.

static void NullMatrixError(const GeneralMatrix* gm)
{
   ((GeneralMatrix&)*gm).tDelete();
   Throw(ProgramException("Maximum or minimum of null matrix"));
}

Real GeneralMatrix::maximum_absolute_value() const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   Real maxval = 0.0; int l = storage; Real* s = store;
   while (l--) { Real a = fabs(*s++); if (maxval < a) maxval = a; }
   ((GeneralMatrix&)*this).tDelete(); return maxval;
}

Real GeneralMatrix::maximum_absolute_value1(int& i) const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   Real maxval = 0.0; int l = storage; Real* s = store; int li = storage;
   while (l--)
      { Real a = fabs(*s++); if (maxval <= a) { maxval = a; li = l; }  }
   i = storage - li;
   ((GeneralMatrix&)*this).tDelete(); return maxval;
}

Real GeneralMatrix::minimum_absolute_value() const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   int l = storage - 1; Real* s = store; Real minval = fabs(*s++);
   while (l--) { Real a = fabs(*s++); if (minval > a) minval = a; }
   ((GeneralMatrix&)*this).tDelete(); return minval;
}

Real GeneralMatrix::minimum_absolute_value1(int& i) const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   int l = storage - 1; Real* s = store; Real minval = fabs(*s++); int li = l;
   while (l--)
      { Real a = fabs(*s++); if (minval >= a) { minval = a; li = l; }  }
   i = storage - li;
   ((GeneralMatrix&)*this).tDelete(); return minval;
}

Real GeneralMatrix::maximum() const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   int l = storage - 1; Real* s = store; Real maxval = *s++;
   while (l--) { Real a = *s++; if (maxval < a) maxval = a; }
   ((GeneralMatrix&)*this).tDelete(); return maxval;
}

Real GeneralMatrix::maximum1(int& i) const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   int l = storage - 1; Real* s = store; Real maxval = *s++; int li = l;
   while (l--) { Real a = *s++; if (maxval <= a) { maxval = a; li = l; } }
   i = storage - li;
   ((GeneralMatrix&)*this).tDelete(); return maxval;
}

Real GeneralMatrix::minimum() const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   int l = storage - 1; Real* s = store; Real minval = *s++;
   while (l--) { Real a = *s++; if (minval > a) minval = a; }
   ((GeneralMatrix&)*this).tDelete(); return minval;
}

Real GeneralMatrix::minimum1(int& i) const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   int l = storage - 1; Real* s = store; Real minval = *s++; int li = l;
   while (l--) { Real a = *s++; if (minval >= a) { minval = a; li = l; } }
   i = storage - li;
   ((GeneralMatrix&)*this).tDelete(); return minval;
}

Real GeneralMatrix::maximum_absolute_value2(int& i, int& j) const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   Real maxval = 0.0; int nr = Nrows();
   MatrixRow mr((GeneralMatrix*)this, LoadOnEntry+DirectPart);
   for (int r = 1; r <= nr; r++)
   {
      int c; maxval = mr.MaximumAbsoluteValue1(maxval, c);
      if (c > 0) { i = r; j = c; }
      mr.Next();
   }
   ((GeneralMatrix&)*this).tDelete(); return maxval;
}

Real GeneralMatrix::minimum_absolute_value2(int& i, int& j) const
{
   REPORT
   if (storage == 0)  NullMatrixError(this);
   Real minval = FloatingPointPrecision::Maximum(); int nr = Nrows();
   MatrixRow mr((GeneralMatrix*)this, LoadOnEntry+DirectPart);
   for (int r = 1; r <= nr; r++)
   {
      int c; minval = mr.MinimumAbsoluteValue1(minval, c);
      if (c > 0) { i = r; j = c; }
      mr.Next();
   }
   ((GeneralMatrix&)*this).tDelete(); return minval;
}

Real GeneralMatrix::maximum2(int& i, int& j) const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   Real maxval = -FloatingPointPrecision::Maximum(); int nr = Nrows();
   MatrixRow mr((GeneralMatrix*)this, LoadOnEntry+DirectPart);
   for (int r = 1; r <= nr; r++)
   {
      int c; maxval = mr.Maximum1(maxval, c);
      if (c > 0) { i = r; j = c; }
      mr.Next();
   }
   ((GeneralMatrix&)*this).tDelete(); return maxval;
}

Real GeneralMatrix::minimum2(int& i, int& j) const
{
   REPORT
   if (storage == 0) NullMatrixError(this);
   Real minval = FloatingPointPrecision::Maximum(); int nr = Nrows();
   MatrixRow mr((GeneralMatrix*)this, LoadOnEntry+DirectPart);
   for (int r = 1; r <= nr; r++)
   {
      int c; minval = mr.Minimum1(minval, c);
      if (c > 0) { i = r; j = c; }
      mr.Next();
   }
   ((GeneralMatrix&)*this).tDelete(); return minval;
}

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

Real Matrix::minimum_absolute_value2(int& i, int& j) const
{
   REPORT
   int k; Real m = GeneralMatrix::minimum_absolute_value1(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::sum_square() const
{
   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 += square(*s++);
      sum1 += square(*s++);
   }
   ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
}

Real SymmetricMatrix::sum_absolute_value() const
{
   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 += fabs(*s++);
      sum1 += fabs(*s++);
   }
   ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
}

Real IdentityMatrix::sum_absolute_value() const
   { REPORT  return fabs(trace()); }    // no need to do tDelete?

Real SymmetricMatrix::sum() const