//$$ hholder.cxx                       Householder triangularisation

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

#define WANT_MATH

#include "include.hxx"

#include "newmatap.hxx"


/********************** householder triangularisation *********************/

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

void HHDecompose(Matrix& X, LowerTriangularMatrix& L)
{
   int n = X.Ncols(); int s = X.Nrows(); L.ReDimension(s);
   real* xi = X.Store(); int k;
   for (int i=0; i<s; i++)
   {
      real sum = 0.0;
      real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
      sum = sqrt(sum);
      L.element(i,i) = sum;
      real* xj0=xi0; k=n; while(k--) { *xj0++ /= sum; }
      for (int j=i+1; j<s; j++)
      {
         sum=0.0;
	 xi=xi0; real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
	 xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
	 L.element(j,i) = sum;
      }
   }
}

void HHDecompose(const Matrix& X, Matrix& Y, Matrix& M)
{
   int n = X.Ncols(); int s = X.Nrows(); int t = Y.Nrows();
   if (Y.Ncols() != n) MatrixError("Incompatible dimensions in HHDecompose");
   M.ReDimension(t,s);
   real* xi = X.Store(); int k;
   for (int i=0; i<s; i++)
   {
      real* xj0 = Y.Store(); real* xi0 = xi;
      for (int j=0; j<t; j++)
      {
         real sum=0.0;
         xi=xi0; real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
	 xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
	 M.element(j,i) = sum;
      }
   }
}

