/** 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;
}