//$$newmatrm.cxx                         rectangular matrix operations

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

#include "include.hxx"

#include "newmat.hxx"
#include "newmatrm.hxx"


// operations on rectangular matrices


void RectMatrixRow::Reset (const Matrix& M, int row, int skip, int length)
{
   RectMatrixRowCol::Reset
      ( M.Store()+row*M.Ncols()+skip, length, 1, M.Ncols() );
}

void RectMatrixRow::Reset (const Matrix& M, int row)
{
   RectMatrixRowCol::Reset( M.Store()+row*M.Ncols(), M.Ncols(), 1, M.Ncols() );
}

void RectMatrixCol::Reset (const Matrix& M, int skip, int col, int length)
{
   RectMatrixRowCol::Reset
      ( M.Store()+col+skip*M.Ncols(), length, M.Ncols(), 1 );
}

void RectMatrixCol::Reset (const Matrix& M, int col)
   { RectMatrixRowCol::Reset( M.Store()+col, M.Nrows(), M.Ncols(), 1 ); }


real RectMatrixRowCol::SumSquare() const
{
   long_real sum = 0.0; int i = n; real* s = store; int d = spacing;
   while (i--) { sum += (long_real)*s * *s; s += d; }
   return (real)sum;
}

real RectMatrixRowCol::operator*(const RectMatrixRowCol& rmrc) const
{
   long_real sum = 0.0; int i = n;
   real* s = store; int d = spacing;
   real* s1 = rmrc.store; int d1 = rmrc.spacing;
   if (i!=rmrc.n) MatrixError("Row/column length error in *");
   while (i--) { sum += (long_real)*s * *s1; s += d; s1 += d1; }
   return (real)sum;
}

void RectMatrixRowCol::AddScaled(const RectMatrixRowCol& rmrc, real r)
{
   int i = n; real* s = store; int d = spacing;
   real* s1 = rmrc.store; int d1 = rmrc.spacing;
   if (i!=rmrc.n) MatrixError("Row/column length error in AddScaled");
   while (i--) { *s += *s1 * r; s += d; s1 += d1; }
}

void RectMatrixRowCol::Divide(const RectMatrixRowCol& rmrc, real r)
{
   int i = n; real* s = store; int d = spacing;
   real* s1 = rmrc.store; int d1 = rmrc.spacing;
   if (i!=rmrc.n) MatrixError("Row/column length error in Divide");
   while (i--) { *s = *s1 / r; s += d; s1 += d1; }
}

void RectMatrixRowCol::Divide(real r)
{
   int i = n; real* s = store; int d = spacing;
   while (i--) { *s /= r; s += d; }
}

void RectMatrixRowCol::Negate()
{
   int i = n; real* s = store; int d = spacing;
   while (i--) { *s = - *s; s += d; }
}

void RectMatrixRowCol::Zero()
{
   int i = n; real* s = store; int d = spacing;
   while (i--) { *s = 0.0; s += d; }
}

void ComplexScale(RectMatrixCol& U, RectMatrixCol& V, real x, real y)
{
   int n = U.n;
   if (n != V.n) MatrixError("Columns different length in ComplexScale");
   real* u = U.store; real* v = V.store; 
   int su = U.spacing; int sv = V.spacing;
   while (n--)
   {
      real z = *u * x - *v * y;  *v =  *u * y + *v * x;  *u = z;
      u += su;  v += sv;
   }
}

void Rotate(RectMatrixCol& U, RectMatrixCol& V, real tau, real s)
{
   //  (U, V) = (U, V) * (c, s)  where  tau = s/(1+c), c^2 + s^2 = 1
   int n = U.n;
   if (n != V.n) MatrixError("Columns different length in Rotate");
   real* u = U.store; real* v = V.store; 
   int su = U.spacing; int sv = V.spacing;
   while (n--)
   {
      real zu = *u; real zv = *v;
      *u -= s * (zv + zu * tau); *v += s * (zu - zv * tau);
      u += su;  v += sv;
   }
}


