#include "libscl.h"

using namespace scl;
using namespace std;

// Gram-Schmidt
// http://web.mit.edu/18.06/www/Essays/gramschmidtmat.pdf

// On return A = Q*R and Q has orthonormal columns

void gs(const realmat& A, realmat& Q, realmat& R)
{
  INTEGER rows = A.nrow();
  INTEGER cols = A.ncol();

  if (cols > rows) error("Error, gs, cols > rows");

  realmat v(rows,1);
  Q.resize(rows,cols,0.0);
  R.resize(cols,cols,0.0);

  for (INTEGER j=1; j<=cols; ++j) {
    for (INTEGER k=1; k<=rows; ++k) v[k] = A(k,j);
    for (INTEGER i=1; i<=j-1; ++i) {
      REAL sum = 0.0;
      for (INTEGER k=1; k<=rows; ++k) sum += Q(k,i)*A(k,j);
      for (INTEGER k=1; k<=rows; ++k) v[k] -= sum*Q(k,i);
      R(i,j) = sum;
    } 
    REAL sum = 0.0;
    for (INTEGER k=1; k<=rows; ++k) sum += pow(v[k],2);
    REAL norm = sqrt(sum);
    for (INTEGER k=1; k<=rows; ++k) Q(k,j) = v[k]/norm;
    R(j,j) = norm;
  }
}

// Random othogonal matrix

void rom(INTEGER rows, INTEGER cols, realmat& Q, INT_32BIT& seed)
{
  realmat A(rows,cols);
  INTEGER l = rows*cols;
  for (INTEGER i=1; i<=l; ++i) A[i] = unsk(seed);
  realmat R;
  gs(A,Q,R);
}

