#include <iostream>
#include <iomanip>

#include "TriMatrix.h"

using namespace std;

TriMatrix::TriMatrix(int N) {
  
}

TriMatrix::TriMatrix(const TriMatrix &T) {
  
}

TriMatrix::~TriMatrix() {
  
}

TriMatrix& TriMatrix::operator= (const TriMatrix &T) {
  
}

double& TriMatrix::operator()(int i, int j) {
  
}

TriMatrix TriMatrix::operator+(const TriMatrix &T) const {
  
}

TriMatrix TriMatrix::operator*(const double a) const {
  
}

TriMatrix TriMatrix::operator-(const TriMatrix &T) const {
  
}

Field TriMatrix::operator*(Field& F) const {
  
}

/** 
 * This "/" operator is really going to be a tridiagonal sovle as we are not
 * going to really invert the matrix but solve it efficiently using the thomas
 * algorithm/Tridiagonal solver 
 * 
 * @param x The RHS of the linear system to solve.
 */
Field TriMatrix::operator/(const Field& x) const {
  return triSolve(x);
}

Field TriMatrix::triSolve(const Field &x) const {
  Field y(N);
  
  int     j;
  double  beta;
  double *gamma = new double[N];
  
  if (mDiag[0] == 0.0) {
    cout << "Error 1 in Tridag()" << endl;
  }
  
  beta = mDiag[0];
  y[0] = x.data[0] / beta;
  for (j = 1; j < N; j++) {
    gamma[j] = mUpper[j - 1] / beta;
    beta = mDiag[j] - mLower[j - 1] * gamma[j];
    if (beta == 0.0) {
      cout << "Error 2 in Tridag()" << endl;
    }
    y[j] = (x.data[j] - mLower[j - 1] * y[j - 1]) / beta;
  }
  
  for (j = N-2; j >= 0; j--) {
    y[j] -= gamma[j + 1] * y[j + 1];
  }
  
  delete[] gamma;
  return y;
}

void TriMatrix::print() {
  int i, j;
  
  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      if (abs(i-j) <= 1) 
        cout << setw(5) << setprecision(2) << (*this)(i,j) << " ";
      else
        cout << setw(5) << setprecision(2) << 0.0 << " ";
    }
    cout << endl;
  }
}

TriMatrix TriMatrix::id(int N) {
  int       i;
  TriMatrix T(N);
  
  for (i = 0; i < N-1; i++) {
    T.mDiag[i]  = 1.0;
    T.mLower[i] = 0.0;
    T.mUpper[i] = 0.0;
  }
  T.mDiag[N-1]  = 1.0;
  
  return T;
}