//$$ newmat8.cxx         Advanced LU transform, scalar functions

// Copyright (C) 1991: R B Davies and DSIR

#define WANT_MATH

#include "include.hxx"

#include "newmatap.hxx"

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

#define REPORT {}


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

void CroutMatrix::ludcmp()
// LU decomposition - from numerical recipes in C
{
   REPORT
   int i,j;

   real* vv=new real [nrows]; MatrixErrorNoSpace(vv);
   real* a;

   a=store;
   for (i=0;i<nrows;i++)
   {
      real big=0.0;
      j=nrows; while (j--) { real temp=fabs(*a++); if (temp > big) big=temp; }
      if (big == 0.0) MatrixError("Singular matrix");
      vv[i]=1.0/big;
   }

   real* aj=store;
   for (j=0;j<nrows;j++)
   {
      real* ai=store;
      for (i=0;i<j;i++)
      {
         real sum=*(ai+j); real* aix=ai; real* ajx=aj;
         int k=i; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
         *ajx = sum; ai += nrows;
      }

      real big=0.0; int imax;
      for (i=j;i<nrows;i++)
      {
         real sum=*(ai+j); real* aix=ai; real* ajx=aj;
         int k=j; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
         *aix = sum; ai += nrows;
         real dum=vv[i]*fabs(sum); if (dum >= big) { big=dum; imax=i; }
      }

      if (j != imax)
      {
         real* amax=store+imax*nrows; real* ajx=store+j*nrows;
         int k=nrows; while(k--) { real dum=*amax; *amax++=*ajx; *ajx++=dum; }
         d=!d; vv[imax]=vv[j];
      }

      indx[j]=imax; ai=aj+j*nrows;
      if (*ai == 0.0) MatrixError("Singular matrix");
      real dum=1.0/(*ai);
      i=nrows-j; while (--i) { ai += nrows; *ai *= dum; }

      aj++;
   }
   delete [nrows] vv;
}

void CroutMatrix::lubksb(real* B, int mini)
{
   REPORT
   int i,j; int ii=-1; real* ai=store;

   for (i=0;i<nrows;i++)
   {
      int ip=indx[i]; real sum=B[ip]; B[ip]=B[i];
      if (ii>=0)
      {
         real* aij=ai+ii; real* bj=B+ii; j=i-ii;
         while (j--) { sum -= *aij++ * *bj; bj++; }
      }
      else if (sum) ii=i;
      B[i]=sum; ai += nrows;
   }

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


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

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

real GeneralMatrix::SumSquare()
{
   REPORT
   real sum = 0.0; int i = storage; real* s = store;
   while (i--) sum += square(*s++);
   tDelete(); return sum;
}

real GeneralMatrix::SumAbsoluteValue()
{
   REPORT
   real sum = 0.0; int i = storage; real* s = store;
   while (i--) sum += fabs(*s++);
   tDelete(); return sum;
}

real GeneralMatrix::MaximumAbsoluteValue()
{
   REPORT
   real maxval = 0.0; int i = storage; real* s = store;
   while (i--) { real a = fabs(*s++); if (maxval < a) maxval = a; }
   tDelete(); return maxval;
}

real SymmetricMatrix::SumSquare()
{
   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++);
   }
   tDelete(); return sum1 + 2.0 * sum2;
}

real SymmetricMatrix::SumAbsoluteValue()
{
   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++);
   }
   tDelete(); return sum1 + 2.0 * sum2;
}

real BaseMatrix::SumSquare()
{
   REPORT GeneralMatrix* gm = Evaluate();
   real s = gm->SumSquare(); return s;
}

real BaseMatrix::SumAbsoluteValue()
{
   REPORT GeneralMatrix* gm = Evaluate();
   real s = gm->SumAbsoluteValue(); return s;
}

real BaseMatrix::MaximumAbsoluteValue()
{
   REPORT GeneralMatrix* gm = Evaluate();
   real s = gm->MaximumAbsoluteValue(); return s;
}

real Matrix::Trace()
{
   REPORT
   int i = nrows; int d = i+1;
   if (i != ncols) MatrixError("trace of non-square matrix");
   real sum = 0.0; real* s = store;
   while (i--) { sum += *s; s += d; }
   tDelete(); return sum;
}

real DiagonalMatrix::Trace()
{
   REPORT
   int i = nrows; real sum = 0.0; real* s = store;
   while (i--) sum += *s++;
   tDelete(); return sum;
}

real SymmetricMatrix::Trace()
{
   REPORT
   int i = nrows; real sum = 0.0; real* s = store; int j = 2;
   while (i--) { sum += *s; s += j++; }
   tDelete(); return sum;
}

real LowerTriangularMatrix::Trace()
{
   REPORT
   int i = nrows; real sum = 0.0; real* s = store; int j = 2;
   while (i--) { sum += *s; s += j++; }
   tDelete(); return sum;
}

real UpperTriangularMatrix::Trace()
{
   REPORT
   int i = nrows; real sum = 0.0; real* s = store;
   while (i) { sum += *s; s += i--; }
   tDelete(); return sum;
}

real BaseMatrix::Trace()
{
   REPORT
   GeneralMatrix* gm = Evaluate(MatrixType::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() { return sign * exp(log_value); }

LogAndSign DiagonalMatrix::LogDeterminant()
{
   REPORT
   int i = nrows; LogAndSign sum; real* s = store;
   while (i--) sum *= *s++;
   tDelete(); return sum;
}

LogAndSign LowerTriangularMatrix::LogDeterminant()
{
   REPORT
   int i = nrows; LogAndSign sum; real* s = store; int j = 2;
   while (i--) { sum *= *s; s += j++; }
   tDelete(); return sum;
}

LogAndSign UpperTriangularMatrix::LogDeterminant()
{
   REPORT
   int i = nrows; LogAndSign sum; real* s = store;
   while (i) { sum *= *s; s += i--; }
   tDelete(); return sum;
}

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

LogAndSign GeneralMatrix::LogDeterminant()
{
   REPORT
   if (nrows != ncols) MatrixError("determinant of non-square matrix"); 
   CroutMatrix C(*this); return C.LogDeterminant();
}

LogAndSign CroutMatrix::LogDeterminant()
{
   REPORT
   int i = nrows; int dd = i+1; LogAndSign sum; real* s = store;
   while (i--) { sum *= *s; s += dd; }
   if (!d) sum.ChangeSign(); return sum;

}
