/** @file matrix.cpp
 *  @brief Function prototypes for creating tridiagonal matrices.
 *
 *
 *
 *  @author Cory Alexander Balaton (coryab)
 *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
 *  @bug No known bugs
 */
#include "matrix.hpp"

arma::mat create_tridiagonal(
    const arma::vec& a, 
    const arma::vec& d, 
    const arma::vec& e)
{
    int n = d.n_elem;
    arma::mat A = arma::mat(n, n).fill(0.);

    A(0, 0) = d(0);
    A(0, 1) = e(0);

    #pragma omp parallel for if(n > 50)
    for (int i = 1; i < n-1; i++) {
        A(i, i-1) = a(i-1);
        A(i, i) = d(i);
        A(i, i+1) = e(i);
    }

    A(n-1, n-2) = a(n-2);
    A(n-1, n-1) = d(n-1);

    return A;
}

arma::mat create_tridiagonal(int n, double a, double d, double e)
{
    arma::vec a_vec = arma::vec(n-1, arma::fill::ones) * a;
    arma::vec d_vec = arma::vec(n, arma::fill::ones) * d;
    arma::vec e_vec = arma::vec(n-1, arma::fill::ones) * e;

    return create_tridiagonal(a_vec, d_vec, e_vec);
}

arma::mat create_symmetric_tridiagonal(int n, double a, double d)
{
    return create_tridiagonal(n, a, d, a);
}

double max_offdiag_symmetric(arma::mat& A, int& k, int& l)
{
    int m = A.n_rows;
    int n = A.n_cols;
    double item, res = 0.;

    for (int i=0; i < m; i++) {
        for (int j=i+1; j < n; j++) {
            if ((item = std::abs(A(i,j))) > res) {
                res = item;
                k = i;
                l = j;
            }
        }
    }

    return res;
}