Category: Numerical Analysis

#include <vector>
#include <random>

class DblLinearAlgebra
{
public:
    static void DblPrintMatrix(
        int m, int n, std::vector<std::vector<double>>& A);
    static void DblAddition(
        size_t m, size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblSubtraction(
        size_t m, size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblMultiply(
        size_t m, size_t n, size_t p,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblAnticommutator(
        size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblCommutator(
        size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static double DblDeterminant(
        int n, int row, int col,
        std::vector<std::vector<double>>& A);
    static bool DblGaussianElimination(
        int m, int n, std::vector<std::vector<double>>& A,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblGaussianFactor(
        int n, std::vector<std::vector<double>>& M,
        std::vector<size_t>& pivot);
    static bool DblGaussianSolution(
        int n, std::vector<std::vector<double>>& M,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblSubstitution(
        int n, std::vector<std::vector<double>>& M,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblInverse(
        int n, std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& A);
    static void DblCharPolyAndAdjoint(
        int n,
        std::vector<std::vector<double>>& C,
        std::vector<std::vector<double>>& I,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& adjoint,
        std::vector<double>& a);
    static void DblMatrixKernel(
        int m, int n,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& X,
        size_t& r);
    static void DblMatrixImage(
        int m, int n,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& N,
        std::vector<std::vector<double>>& X,
        int rank);
    static void DblGenerateNonSingular(
        double scale, double& determinant,
        int n, unsigned int seed,
        std::vector<std::vector<double>>& Mr);
};

#include "DblLinearAlgebra.h"
#include <iomanip>
#include <iostream>

void DblLinearAlgebra::DblPrintMatrix(
    int m, int n, std::vector<std::vector<double>>& A)
{
    for (int i = 0; i < m; i++)
    {
        for (int j = 0; j < n; j++)
        {
            std::cout << std::setprecision(6) << std::setw(9);

            if (fabs(A[i][j]) > 1.0e-12)
            {
                std::cout << A[i][j] << ' ';
            }

            else
            {
                std::cout << 0 << ' ';
            }
        }

        std::cout << std::endl;
    }
}

void DblLinearAlgebra::DblAddition(
    size_t m, size_t n,
    std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& B,
    std::vector<std::vector<double>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            C[i][j] = A[i][j] + B[i][j];
        }
    }
}

void DblLinearAlgebra::DblSubtraction(
    size_t m, size_t n,
    std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& B,
    std::vector<std::vector<double>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            C[i][j] = A[i][j] - B[i][j];
        }
    }
}

void DblLinearAlgebra::DblMultiply(
    size_t m, size_t n, size_t p,
    std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& B,
    std::vector<std::vector<double>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            double sum = { 0 };

            for (size_t k = 0; k < p; k++)
            {
                sum += A[i][k] * B[k][j];
            }

            C[i][j] = sum;
        }
    }
}

void DblLinearAlgebra::DblAnticommutator(
    size_t n,
    std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& B,
    std::vector<std::vector<double>>& C)
{
    std::vector<std::vector<double>> D(n,
        std::vector<double>(n));
    std::vector<std::vector<double>> E(n,
        std::vector<double>(n));

    DblMultiply(n, n, n, A, B, D);
    DblMultiply(n, n, n, B, A, E);
    DblAddition(n, n, D, E, C);
}

void DblLinearAlgebra::DblCommutator(
    size_t n,
    std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& B,
    std::vector<std::vector<double>>& C)
{
    std::vector<std::vector<double>> D(n,
        std::vector<double>(n));
    std::vector<std::vector<double>> E(n,
        std::vector<double>(n));

    DblMultiply(n, n, n, A, B, D);
    DblMultiply(n, n, n, B, A, E);
    DblSubtraction(n, n, D, E, C);
}

// https://www.geeksforgeeks.org/dsa/determinant-of-a-matrix/

double getDet(std::vector<std::vector<double>>& mat, int n) {

    // Base case: if the matrix is 1x1
    if (n == 1) {
        return mat[0][0];
    }

    // Base case for 2x2 matrix
    if (n == 2) {
        return mat[0][0] * mat[1][1] -
            mat[0][1] * mat[1][0];
    }

    // Recursive case for larger matrices
    double  res = 0;
    for (int col = 0; col < n; ++col) {

        // Create a submatrix by removing the first 
        // row and the current column
        std::vector<std::vector<double>> sub(n - 1,
            std::vector<double>(n - 1));
        for (int i = 1; i < n; ++i) {
            int subcol = 0;
            for (int j = 0; j < n; ++j) {

                // Skip the current column
                if (j == col) continue;

                // Fill the submatrix
                sub[i - 1LL][subcol++] = mat[i][j];
            }
        }

        // Cofactor expansion
        int sign = (col % 2 == 0) ? 1 : -1;
        res += sign * mat[0][col] * getDet(sub, n - 1);
    }

    return res;
}

double DblLinearAlgebra::DblDeterminant(
    int n, int row, int col,
    std::vector<std::vector<double>>& A)
{
    return getDet(A, A.size());
}

bool DblLinearAlgebra::DblGaussianElimination(
    int m, int n, std::vector<std::vector<double>>& M,
    std::vector<double>& b, std::vector<double>& x, 
    std::vector<size_t>& pivot)
{
    bool failure = false;
    std::vector<double> c(m);

    b.resize(n);
    x.resize(n);

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

    for (size_t j = 0; j < n; j++)
    {
        bool found = false;
        size_t i = j;

        while (i < n && !found)
        {
            if (M[i][j] != 0)
                found = true;
            else
                i++;
        }

        if (!found)
        {
            failure = true;
            break;
        }

        if (i > j)
        {
            for (size_t l = j; l < n; l++)
            {
                double t = M[i][l];
                M[i][l] = M[j][l];
                M[j][l] = t;
                t = b[i];
                b[i] = b[j];
                b[j] = t;
            }
        }

        double d = 1.0 / M[j][j];

        for (size_t k = j + 1; k < n; k++)
            c[k] = d * M[k][j];

        for (size_t k = j + 1; k < n; k++)
        {
            for (size_t l = j + 1; l < n; l++)
                M[k][l] = M[k][l] - c[k] * M[j][l];

            b[k] = b[k] - c[k] * b[j];
        }
    }

    for (long long i = (long long)n - 1; i >= 0; i--)
    {
        double sum = 0;

        for (size_t j = i + 1; j < n; j++)
            sum += M[i][j] * x[j];

        x[i] = (b[i] - sum) / M[i][i];
    }

    return failure;
}

bool DblLinearAlgebra::DblGaussianFactor(
    int n, std::vector<std::vector<double>>& M,
    std::vector<size_t>& pivot)
{
    // returns false if matrix is singular
    std::vector<double> d(n);
    double awikod, col_max, ratio, row_max, temp;
    int flag = 1;
    size_t i_star, itemp;

    for (size_t i = 0; i < n; i++)
    {
        pivot[i] = i;
        row_max = 0;
        
        for (size_t j = 0; j < n; j++)
            row_max = fmax(row_max, fabs(M[i][j]));
        
        if (row_max == 0)
        {
            flag = 0;
            row_max = 1.0;
        }

        d[i] = row_max;
    }

    if (n <= 1) return flag != 0;
    
    // factorization
    
    for (size_t k = 0; k < (size_t)(n - 1LL); k++)
    {
        // determine pivot row the row i_star

        col_max = fabs(M[k][k]) / d[k];
        i_star = k;

        for (size_t i = k + 1; i < n; i++)
        {
            awikod = fabs(M[i][k]) / d[i];

            if (awikod > col_max)
            {
                col_max = awikod;
                i_star = i;
            }
        }

        if (col_max == 0)
            flag = 0;
        
        else
        {
            if (i_star > k)
            {
                // make k the pivot row by
                // interchanging with i_star

                flag *= -1;
                itemp = pivot[i_star];
                pivot[i_star] = pivot[k];
                pivot[k] = itemp;
                temp = d[i_star];
                d[i_star] = d[k];
                d[k] = temp;

                for (size_t j = 0; j < n; j++)
                {
                    temp = M[i_star][j];
                    M[i_star][j] = M[k][j];
                    M[k][j] = temp;
                }
            }

            // eliminate x[k]

            for (size_t i = k + 1; i < n; i++)
            {
                M[i][k] /= M[k][k];
                ratio = M[i][k];

                for (size_t j = k + 1; j < n; j++)
                    M[i][j] -= ratio * M[k][j];
            }
        }

        if (M[n - 1LL][n - 1LL] == 0) flag = 0;
    }

    if (flag == 0)
        return false;

    return true;
}

bool DblLinearAlgebra::DblGaussianSolution(
    int n, std::vector<std::vector<double>>& M,
    std::vector<double>& b, std::vector<double>& x,
    std::vector<size_t>& pivot)
{
    if (!DblGaussianFactor(n, M, pivot))
        return false;

    return DblSubstitution(n, M, b, x, pivot);
}

bool DblLinearAlgebra::DblSubstitution(
    int n, std::vector<std::vector<double>>& M,
    std::vector<double>& b, std::vector<double>& x,
    std::vector<size_t>& pivot)
{
    double sum = 0.0;
    size_t n1 = n - 1LL;

    if (n == 1)
    {
        x[0] = b[0] / M[0][0];
        return true;
    }

    // forward substitution

    x[0] = b[pivot[0]];

    for (size_t i = 1; i < n; i++)
    {
        double sum = 0.0;

        for (size_t j = 0; j < i; j++)
            sum += M[i][j] * x[j];

        x[i] = b[pivot[i]] - sum;
    }

    // backward substitution

    x[n1] /= M[n1][n1];

    for (long long i = n - 2LL; i >= 0; i--)
    {
        double sum = 0.0;

        for (size_t j = i + 1; j < n; j++)
            sum += M[i][j] * x[j];

        x[i] = (x[i] - sum) / M[i][i];
    }

    return true;
}

bool DblLinearAlgebra::DblInverse(
    int n, std::vector<std::vector<double>>& M,
    std::vector<std::vector<double>>& Mi)
{
    std::vector<double> b(n);
    std::vector<double> x(n);
    std::vector<size_t> pivot(n);
    std::vector<std::vector<double>> Mc(n,
        std::vector<double>(n));

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            Mc[i][j] = M[i][j];
        }
    }

    if (!DblGaussianFactor(n, Mc, pivot))
        return false;

    for (size_t i = 0; i < n; i++)
    {
        b[i] = 0;
    }

    for (size_t i = 0; i < n; i++)
    {
        b[i] = 1;

        if (!DblSubstitution(n, Mc, b, x, pivot))
            return false;

        b[i] = 0;

        for (size_t j = 0; j < n; j++)
            Mi[j][i] = x[j];
    }

    return true;
}

void DblLinearAlgebra::DblCharPolyAndAdjoint(
    int n,
    std::vector<std::vector<double>>& C,
    std::vector<std::vector<double>>& I,
    std::vector<std::vector<double>>& M,
    std::vector<std::vector<double>>& adjoint,
    std::vector<double>& a)
{
    C.resize(n, std::vector<double>(n));
    I.resize(n, std::vector<double>(n));

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            C[i][j] = I[i][j] = 0;
    }

    for (size_t i = 0; i < n; i++)
        C[i][i] = I[i][i] = 1;

    a[0] = 1;

    for (size_t i = 1; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
            {
                double sum = 0.0;

                for (size_t l = 0; l < n; l++)
                    sum += M[j][l] * C[l][k];

                C[j][k] = sum;
            }
        }

        double tr = 0.0;

        for (size_t j = 0; j < n; j++)
            tr += C[j][j];

        a[i] = -tr / i;

        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
                C[j][k] += a[i] * I[j][k];
        }
    }

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            double sum = 0.0;

            for (size_t k = 0; k < n; k++)
                sum += M[i][k] * C[k][j];

            C[i][j] = sum;
        }
    }

    double trace = 0.0;

    for (size_t i = 0; i < n; i++)
        trace += C[i][i];

    trace /= n;
    a[n - 1LL] = -trace;

    double factor = 1.0;

    if ((n - 1) % 2 != 0)
        factor = -1.0;

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            adjoint[i][j] = factor * C[i][j];
    }
}

void DblLinearAlgebra::DblMatrixKernel(
    int m, int n,
    std::vector<std::vector<double>>& M,
    std::vector<std::vector<double>>& X,
    size_t& r)
{
    double D = 0.0;
    std::vector <int> c(m);
    std::vector <int> d(n);

    r = 0;

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

    size_t j, k = 1;

Step2:
    
    for (j = 0; j < m; j++)
    {
        if (M[j][k] != 0 && c[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d[k] = 0;
        goto Step4;
    }

    D = -1.0 / M[j][k];

    M[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        M[j][s] = D * M[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = M[i][k];
                M[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            M[i][s] += D * M[j][s];
        }
    }

    c[j] = (int)k;
    d[k] = (int)j;

Step4:
    
    if (k < n - 1)
    {
        k++;
        goto Step2;
    }

    X.resize(n, std::vector<double>(n));

    if (r != 0)
    {
        for (k = 0; k < n; k++)
        {
            if (d[k] == 0)
            {
                for (size_t i = 0; i < n; i++)
                {
                    if (d[i] > 0)
                        X[k][i] = M[d[i]][k];
                    else if (i == k)
                        X[k][i] = 1;
                    else
                        X[k][i] = 0;
                }
            }
        }
    }
}

void DblLinearAlgebra::DblMatrixImage(
    int m, int n,
    std::vector<std::vector<double>>& M,
    std::vector<std::vector<double>>& N,
    std::vector<std::vector<double>>& X,
    int rank)
{
    double D = 0.0;
    size_t r = 0;
    std::vector<std::vector<double>> copyM(
        m, std::vector<double>(n));
    std::vector <int> c(m);
    std::vector <int> d(n);

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

    size_t j = 0, k = 1;

    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            N[i][j] = copyM[i][j] = M[i][j];
        }
    }

Step2:
    
    for (size_t j = 0; j < m; j++)
    {
        if (copyM[j][k] != 0 && c[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d[k] = 0;
        goto Step4;
    }

    D = -1.0 / copyM[j][k];

    copyM[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        copyM[j][s] = D * copyM[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = copyM[i][k];
                copyM[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            copyM[i][s] += D * copyM[j][s];
        }
    }

    c[j] = (int)k;
    d[k] = (int)j;

Step4:
    
    if (k < (size_t)(n - 1LL))
    {
        k++;
        goto Step2;
    }

    rank = (int)(n - r) ;

    for (j = 0; j < m; j++)
    {
        if (c[j] != 0)
        {
            for (size_t i = 0; i < m; i++)
            {
                N[i][c[j]] = M[i][c[j]];
            }
        }
    }
}

void DblLinearAlgebra::DblGenerateNonSingular(
    double scale, double& determinant,
    int n, unsigned int seed,
    std::vector<std::vector<double>>& Mr)
{
    bool failure = false;
    std::mt19937 rng(seed);
    std::uniform_real_distribution<double> dist(0.0, 1.0);

    while (true)
    {
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                Mr[i][j] = scale * dist(rng);
            }
        }

        determinant = DblDeterminant(n, 0, 0, Mr);
        failure = determinant == 0;

        if (!failure)
            return;
    }
}

#include <complex>
#include <vector>

class CmpLinearAlgebra
{
public:
    static void CmpPrintMatrix(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& Ac);
    static void CmpAddition(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpSubtraction(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpMultiply(
        size_t m, size_t n, size_t p,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpAnticommutator(
        size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpCommutator(
        size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static std::complex<double> CmpDeterminant(
        int n,
        std::vector<std::vector<std::complex<double>>>& Ac);
    static void CmpAdjoint(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& Ac,
        std::vector<std::vector<std::complex<double>>>& Ad);
    static bool CmpGaussianElimination(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpGaussianFactor(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<size_t>& pivot);
    static bool CmpGaussianSolution(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpSubstitution(
        int m, int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpInverse(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& Mi);
    static void CmpCharPolyAndAdjoint(
        int n,
        std::vector<std::vector<std::complex<double>>>& C,
        std::vector<std::vector<std::complex<double>>>& I,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& adjoint,
        std::vector<std::complex<double>>& a);
    static void CmpMatrixKernel(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& X,
        size_t& r);
    static void CmpMatrixImage(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& N,
        std::vector<std::vector<std::complex<double>>>& X,
        int rank);
    static void CmpGenerateNonSingular(
        double scale, std::complex<double>& determinant,
        int n, unsigned int seed,
        std::vector<std::vector<std::complex<double>>>& Mc);
};

#include "CmpLinearAlgebra.h"
#include <iomanip>
#include <iostream>
#include <random>

void CmpLinearAlgebra::CmpPrintMatrix(
    int m, int n,
    std::vector<std::vector<std::complex<double>>>& Ac)
{
    for (int i = 0; i < m; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (Ac[i][j]._Val[0] >= 0)
                std::cout << '+';
            else
                std::cout << '-';

            std::cout << std::setprecision(6) << std::setw(9);

            if (fabs(Ac[i][j]._Val[0]) > 1.0e-12)
            {
                std::cout << fabs(Ac[i][j]._Val[0]) << ' ';
            }

            else
            {
                std::cout << 0 << ' ';
            }

            if (Ac[i][j]._Val[1] >= 0)
                std::cout << '+';
            else
                std::cout << '-';

            std::cout << std::setprecision(6) << std::setw(9);

            if (fabs(Ac[i][j]._Val[1]) > 1.0e-12)
            {
                std::cout << fabs(Ac[i][j]._Val[1]) << "i\t";
            }

            else
            {
                std::cout << 0 << "i\t";
            }
        }

        std::cout << std::endl;
    }
}

void CmpLinearAlgebra::CmpAddition(
    size_t m, size_t n,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::vector<std::complex<double>>>& B,
    std::vector<std::vector<std::complex<double>>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            C[i][j] = A[i][j] + B[i][j];
        }
    }
}

void CmpLinearAlgebra::CmpSubtraction(
    size_t m, size_t n,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::vector<std::complex<double>>>& B,
    std::vector<std::vector<std::complex<double>>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            C[i][j] = A[i][j] - B[i][j];
        }
    }
}

void CmpLinearAlgebra::CmpMultiply(
    size_t m, size_t n, size_t p,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::vector<std::complex<double>>>& B,
    std::vector<std::vector<std::complex<double>>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            std::complex<double> sum = 0;

            for (size_t k = 0; k < p; k++)
            {
                sum += A[i][k] * B[k][j];
            }

            C[i][j] = sum;
        }
    }
}

void CmpLinearAlgebra::CmpAnticommutator(
    size_t n,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::vector<std::complex<double>>>& B,
    std::vector<std::vector<std::complex<double>>>& C)
{
    std::vector<std::vector<std::complex<double>>> D(n,
        std::vector<std::complex<double>>(n));
    std::vector<std::vector<std::complex<double>>> E(n,
        std::vector<std::complex<double>>(n));

    CmpMultiply(n, n, n, A, B, D);
    CmpMultiply(n, n, n, B, A, E);
    CmpAddition(n, n, D, E, C);
}

void CmpLinearAlgebra::CmpCommutator(
    size_t n,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::vector<std::complex<double>>>& B,
    std::vector<std::vector<std::complex<double>>>& C)
{
    std::vector<std::vector<std::complex<double>>> D(n,
        std::vector<std::complex<double>>(n));
    std::vector<std::vector<std::complex<double>>> E(n,
        std::vector<std::complex<double>>(n));

    CmpMultiply(n, n, n, A, B, D);
    CmpMultiply(n, n, n, B, A, E);
    CmpSubtraction(n, n, D, E, C);
}

void CmpLinearAlgebra::CmpAdjoint(
    size_t m, size_t n,
    std::vector<std::vector<std::complex<double>>>& Ac,
    std::vector<std::vector<std::complex<double>>>& Ad)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            Ad[j][i] = std::conj(Ac[i][j]);
        }
    }
}

// https://www.geeksforgeeks.org/dsa/determinant-of-a-matrix/

std::complex<double> getDet(
    std::vector<std::vector<std::complex<double>>>& mat, int n) {

    // Base case: if the matrix is 1x1
    if (n == 1) {
        return mat[0][0];
    }

    // Base case for 2x2 matrix
    if (n == 2) {
        return mat[0][0] * mat[1][1] -
            mat[0][1] * mat[1][0];
    }

    // Recursive case for larger matrices
    std::complex<double> res = 0;
    for (int col = 0; col < n; ++col) {

        // Create a submatrix by removing the first 
        // row and the current column
        std::vector<std::vector<std::complex<double>>> sub(n - 1,
            std::vector<std::complex<double>>(n - 1));
        for (int i = 1; i < n; ++i) {
            int subcol = 0;
            for (int j = 0; j < n; ++j) {

                // Skip the current column
                if (j == col) continue;

                // Fill the submatrix
                sub[i - 1LL][subcol++] = mat[i][j];
            }
        }

        // Cofactor expansion
        int sign = (col % 2 == 0) ? 1 : -1;
        std::complex<double> csign(sign, 0.0);
        res = res + csign * mat[0][col] * getDet(sub, n - 1);
    }

    return res;
}

std::complex<double> CmpLinearAlgebra::CmpDeterminant(
    int n, std::vector<std::vector<std::complex<double>>>& A)
{
    return getDet(A, A.size());
}

bool CmpLinearAlgebra::CmpGaussianElimination(
    int m, int n,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::complex<double>>& b,
    std::vector<std::complex<double>>& x,
    std::vector<size_t>& pivot)
{
    bool failure = false;
    std::vector<std::complex<double>> c(m);

    b.resize(n);
    x.resize(n);

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

    for (size_t j = 0; j < n; j++)
    {
        bool found = false;
        size_t i = j;

        while (i < n && !found)
        {
            if (abs(A[i][j]) != 0)
                found = true;
            else
                i++;
        }

        if (!found)
        {
            failure = true;
            break;
        }

        if (i > j)
        {
            for (size_t l = j; l < n; l++)
            {
                std::complex<double> t = A[i][l];
                A[i][l] = A[j][l];
                A[j][l] = t;
                t = b[i];
                b[i] = b[j];
                b[j] = t;
            }
        }

        std::complex<double> d = 1.0 / A[j][j];

        for (size_t k = j + 1; k < n; k++)
            c[k] = d * A[k][j];

        for (size_t k = j + 1; k < n; k++)
        {
            for (size_t l = j + 1; l < n; l++)
                A[k][l] = A[k][l] - c[k] * A[j][l];

            b[k] = b[k] - c[k] * b[j];
        }
    }

    for (long long i = (long long)n - 1; i >= 0; i--)
    {
        std::complex<double> sum = 0;

        for (size_t j = i + 1; j < n; j++)
            sum += A[i][j] * x[j];

        x[i] = (b[i] - sum) / A[i][i];
    }

    return failure;
}

bool CmpLinearAlgebra::CmpSubstitution(
    int m, int n, std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::complex<double>>& b,
    std::vector<std::complex<double>>& x,
    std::vector<size_t>& pivot)
{
    std::complex<double> sum = 0;
    size_t n1 = n - 1LL;

    if (n == 1)
    {
        x[0] = b[0] / M[0][0];
        return true;
    }

    // forward substitution

    x[0] = b[pivot[0]];

    for (size_t i = 1; i < n; i++)
    {
        std::complex<double> sum = 0;

        for (size_t j = 0; j < i; j++)
            sum += M[i][j] * x[j];

        x[i] = b[pivot[i]] - sum;
    }

    // backward substitution

    x[n1] /= M[n1][n1];

    for (long long i = n - 2LL; i >= 0; i--)
    {
        std::complex<double> sum = 0;

        for (size_t j = i + 1; j < n; j++)
            sum += M[i][j] * x[j];

        x[i] = (x[i] - sum) / M[i][i];
    }

    return true;
}

static std::complex<double> complex_max(
    std::complex<double> a, std::complex<double> b) {
    return (std::abs(a) > std::abs(b)) ? a : b;
}

bool CmpLinearAlgebra::CmpGaussianFactor(
    int n, std::vector<std::vector<std::complex<double>>>& M,
    std::vector<size_t>& pivot)
{
    // returns false if matrix is singular
    std::vector<std::complex<double>> d(n);
    std::complex<double> awikod = 0, col_max = 0, ratio = 0, row_max = 0, temp = 0;
    int flag = 1;
    size_t i_star, itemp;

    for (size_t i = 0; i < n; i++)
    {
        pivot[i] = i;
        row_max = 0;

        for (size_t j = 0; j < n; j++)
            row_max = complex_max(row_max, abs(M[i][j]));
        
        if (abs(row_max) == 0)
        {
            flag = 0;
            row_max = 1;
        }

        d[i] = row_max;
    }
    if (n <= 1) return flag != 0;
    
    // factorization
    
    for (size_t k = 0; k < (size_t)n - 1LL; k++)
    {
        // determine pivot row the row i_star

        col_max = abs(M[k][k]) / d[k];
        i_star = k;

        for (size_t i = k + 1; i < n; i++)
        {
            awikod = abs(M[i][k]) / d[i];

            if (abs(awikod) > abs(col_max))
            {
                col_max = awikod;
                i_star = i;
            }
        }
        
        if (abs(col_max) == 0)
            flag = 0;
        
        else
        {
            if (i_star > k)
            {
                // make k the pivot row by
                // interchanging with i_star
                flag *= -1;
                itemp = pivot[i_star];
                pivot[i_star] = pivot[k];
                pivot[k] = itemp;
                temp = d[i_star];
                d[i_star] = d[k];
                d[k] = temp;

                for (size_t j = 0; j < n; j++)
                {
                    temp = M[i_star][j];
                    M[i_star][j] = M[k][j];
                    M[k][j] = temp;
                }
            }

            // eliminate x[k]
            
            for (size_t i = k + 1; i < n; i++)
            {
                M[i][k] /= M[k][k];
                ratio = M[i][k];
                
                for (size_t j = k + 1; j < n; j++)
                    M[i][j] -= ratio * M[k][j];
            }
        }

        if (abs(M[n - 1LL][n - 1LL]) == 0) flag = 0;
    }

    if (flag == 0)
        return false;

    return true;
}

bool CmpLinearAlgebra::CmpGaussianSolution(
    int n, std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::complex<double>>& b,
    std::vector<std::complex<double>>& x,
    std::vector<size_t>& pivot)
{
    if (!CmpGaussianFactor(n, M, pivot))
        return false;

    return CmpSubstitution(n, n, M, b, x, pivot);
}

bool CmpLinearAlgebra::CmpInverse(
    int n, std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::vector<std::complex<double>>>& Mi)
{
    std::vector<std::complex<double>> b(n);
    std::vector<std::complex<double>> x(n);
    std::vector<size_t> pivot(n);
    std::vector<std::vector<std::complex<double>>> Mc(n,
        std::vector<std::complex<double>>(n));

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            Mc[i][j] = M[i][j];
        }
    }

    if (!CmpGaussianFactor(n, Mc, pivot))
        return false;

    for (size_t i = 0; i < n; i++)
    {
        b[i] = 0;
    }

    for (size_t i = 0; i < n; i++)
    {
        b[i] = 1;

        if (!CmpSubstitution(n, n, Mc, b, x, pivot))
            return false;

        b[i] = 0;

        for (size_t j = 0; j < n; j++)
            Mi[j][i] = x[j];
    }

    return true;
}

void CmpLinearAlgebra::CmpCharPolyAndAdjoint(
    int n,
    std::vector<std::vector<std::complex<double>>>& C,
    std::vector<std::vector<std::complex<double>>>& I,
    std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::vector<std::complex<double>>>& adjoint,
    std::vector<std::complex<double>>& a)
{
    C.resize(n, std::vector<std::complex<double>>(n));
    I.resize(n, std::vector<std::complex<double>>(n));

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            C[i][j] = I[i][j] = 0;
    }

    for (size_t i = 0; i < n; i++)
        C[i][i] = I[i][i] = 1;

    a[0] = 1;

    for (size_t i = 1; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
            {
                std::complex<double> sum = 0.0;

                for (size_t l = 0; l < n; l++)
                    sum += M[j][l] * C[l][k];

                C[j][k] = sum;
            }
        }

        std::complex<double> tr = 0.0;

        for (size_t j = 0; j < n; j++)
            tr += C[j][j];

        std::complex<double> ci = 0;
        ci._Val[0] = (double)i;

        a[i] = -tr / ci;

        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
                C[j][k] += a[i] * I[j][k];
        }
    }

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            std::complex<double> sum = 0.0;

            for (size_t k = 0; k < n; k++)
                sum += M[i][k] * C[k][j];

            C[i][j] = sum;
        }
    }

    std::complex<double> trace = 0.0;

    for (size_t i = 0; i < n; i++)
        trace += C[i][i];

    trace /= n;
    a[n - 1LL] = -trace;

    std::complex<double> factor = 1.0;

    if ((n - 1) % 2 != 0)
        factor = -1.0;

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            adjoint[i][j] = factor * C[i][j];
    }
}

void CmpLinearAlgebra::CmpMatrixKernel(
    int m, int n,
    std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::vector<std::complex<double>>>& X,
    size_t& r)
{
    std::complex<double> D = 0;
    std::vector<int> c(m);
    std::vector<int> d(n);

    r = 0;

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

    size_t j = 0, k = 1;

Step2:

    for (j = 0; j < m; j++)
    {
        if (abs(M[j][k]) != 0 && c[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d[k] = 0;
        goto Step4;
    }

    D = -1.0 / M[j][k];

    M[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        M[j][s] = D * M[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = M[i][k];
                M[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            M[i][s] += D * M[j][s];
        }
    }

    c[j] = (int)k;
    d[k] = (int)j;

Step4:

    if (k < n - 1)
    {
        k++;
        goto Step2;
    }

    X.resize(n, std::vector<std::complex<double>>(n));

    if (r != 0)
    {
        for (k = 0; k < n; k++)
        {
            if (d[k] == 0)
            {
                for (size_t i = 0; i < n; i++)
                {
                    if (d[i] > 0)
                        X[k][i] = M[d[i]][k];
                    else if (i == k)
                        X[k][i] = 1;
                    else
                        X[k][i] = 0;
                }
            }
        }
    }
}

void CmpLinearAlgebra::CmpMatrixImage(
    int m, int n,
    std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::vector<std::complex<double>>>& N,
    std::vector<std::vector<std::complex<double>>>& X,
    int rank)
{
    std::complex<double> D = 0.0;
    size_t r = 0;
    std::vector<std::vector<std::complex<double>>> copyM(
        m, std::vector<std::complex<double>>(n));
    std::vector<int> c(m);
    std::vector<int> d(n);

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

    size_t j = 0, k = 1;

    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            N[i][j] = copyM[i][j] = M[i][j];
        }
    }

Step2:

    for (size_t j = 0; j < m; j++)
    {
        if (abs(copyM[j][k]) != 0 && c[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d[k] = 0;
        goto Step4;
    }

    D = -1.0 / copyM[j][k];

    copyM[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        copyM[j][s] = D * copyM[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = copyM[i][k];
                copyM[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            copyM[i][s] += D * copyM[j][s];
        }
    }

    c[j] = (int)k;
    d[k] = (int)j;

Step4:
    
    if (k < (size_t)(n - 1LL))
    {
        k++;
        goto Step2;
    }

    rank = (int)(n - r);

    for (j = 0; j < m; j++)
    {
        if (c[j] != 0)
        {
            for (size_t i = 0; i < m; i++)
            {
                N[i][c[j]] = M[i][c[j]];
            }
        }
    }
}

void CmpLinearAlgebra::CmpGenerateNonSingular(
    double scale, std::complex<double>& cDeterminant,
    int n, unsigned int seed,
    std::vector<std::vector<std::complex<double>>>& Mc)
{
    bool failure = false;
    std::mt19937 rng(seed);
    std::uniform_real_distribution<double> dist(0.0, 1.0);

    while (true)
    {
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                Mc[i][j]._Val[0] = scale * dist(rng);
                Mc[i][j]._Val[1] = scale * dist(rng);
            }
        }

       cDeterminant = CmpDeterminant(n, Mc);

       if (cDeterminant._Val[0] != 0 || cDeterminant._Val[1] != 0)
           break;
    }
}

// Exercises from "Modern Quantum Chemistry an Introduction to Advanced
// Electronic Structure Theory" by Attila Szabo and Neil S. Ostlund
// https://chemistlibrary.wordpress.com/wp-content/uploads/2015/02/modern-quantum-chemistry.pdf
// Program (c) Tuesday, August 12, 2025 by James Pate Williams, Jr.
// Program (c) Tuesday, August 19, 2025 by James Pate Williams, Jr.

#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>
#include "DblLinearAlgebra.h"
#include "CmpLinearAlgebra.h"

int main()
{
	// static data matrices
	double AArcb[3][3] = { { 2, 3, -1 }, { 4, 4, -3 }, { -2, 3, -1 } };
	double AArso[3][3] = { { 1, 1, 0 }, { 1, 2, 2 }, { 0, 2, -1 } };
	double BBrso[3][3] = { { 1, -1, 1 }, { -1, 0, 0 }, { 1, 0, 1} };
	double BBr[3][3] = { { 1, -1, 1 }, { -1 , 0, 0 }, { 1, 0, 1 } };
	double AAcr[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
	double AAci[3][3] = { { 1, 1, 2 }, { 3, 0, 1 }, { 0, 2, 4 } };
	double BBcr[3][3] = { { 1, 0, 1 }, { 1 , 1, 0 }, { 0, 1, 1 } };
	double BBci[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
	// some array dimensions
	int m = 3, n = 3, p = 3;
	// a couple of 3x1 vectors
	std::vector<double> br(3);
	std::vector<size_t> pivot(3);
	// 3x3 real matrices
	std::vector<std::vector<double>> Arcb(3, std::vector<double>(3));
	std::vector<std::vector<double>> Arso(3, std::vector<double>(3));
	std::vector<std::vector<double>> Brso(3, std::vector<double>(3));
	std::vector<std::vector<double>> Br(3, std::vector<double>(3));
	std::vector<std::vector<double>> Cr(3, std::vector<double>(3));
	std::vector<std::vector<double>> Ai(3, std::vector<double>(3));
	std::vector<std::vector<double>> Ari(3, std::vector<double>(3));
	// a 4x4 real matrix
	std::vector<std::vector<double>> Mr(4, std::vector<double>(4));
	// 3x3 complex matrices
	std::vector<std::vector<std::complex<double>>> Ac(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Bc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Cc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Dc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Ec(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Fc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Gc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Mc(4,
		std::vector<std::complex<double>>(4));
	// copy static real matrices to dynamic matrices
	for (int i = 0; i < m; i++)
	{
		for (int j = 0; j < p; j++)
		{
			Arcb[i][j] = AArcb[i][j];
			Arso[i][j] = AArso[i][j];
			Brso[i][j] = BBrso[i][j];
			Ac[i][j]._Val[0] = AAcr[i][j];
			Ac[i][j]._Val[1] = AAci[i][j];
		}
	}
	// copy static complex matrices to dynamic matrices
	for (int i = 0; i < p; i++)
	{
		for (int j = 0; j < n; j++)
		{
			Br[i][j] = BBr[i][j];
			Bc[i][j]._Val[0] = BBcr[i][j];
			Bc[i][j]._Val[1] = BBci[i][j];
		}
	}
	// See "Elementary Numerical Analysis an
	// Algorithmic Approach" (c) 1980 by S. D. Conte
	// and Carl de Boor
	DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Br, Cr);
	std::cout << "Ar * Br = Cr Conte & de Boor" << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	// complex matrix multiplication
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
	std::cout << "Ac * Bc = Cc" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
	std::cout << std::endl;
	// Exercise 1.2 from Szabo and Ostlund
	std::cout << "Exercise 1.2 page 5 Commutator" << std::endl;
	DblLinearAlgebra::DblCommutator(3, Arso, Brso, Cr);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	std::cout << "Exercise 1.2 page 5 Anticommutator" << std::endl;
	DblLinearAlgebra::DblAnticommutator(3, Arso, Brso, Cr);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpAdjoint(3, 3, Cc, Dc);
	std::cout << "Exercise 1.3 page 6 Cc adjoint" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Dc);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpAdjoint(3, 3, Ac, Ec);
	CmpLinearAlgebra::CmpAdjoint(3, 3, Bc, Fc);
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Fc, Ec, Gc);
	std::cout << "Exercise 1.3 page 6 Bc adjoint * Ac adjoint"
		<< std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Gc);
	std::cout << std::endl;
	std::cout << "Ar matrix" << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Arcb);
	bool inv = DblLinearAlgebra::DblInverse(n, Arcb, Ai);
	std::cout << std::endl;
	std::cout << "Ar Conte & de Boor inverse flag = "
		<< inv << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Ai);
	std::cout << std::endl;
	std::cout << "Ar * Ar inverse" << std::endl;
	DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Ai, Ari);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Ari);
	std::cout << std::endl;
	double rDeterminant = 0;
	DblLinearAlgebra::DblGenerateNonSingular(
		2.0, rDeterminant, 4, 1, Mr);
	std::cout << "rDeterminant = ";
	std::cout << rDeterminant << std::endl;
	std::cout << std::endl;
	std::cout << "Ac" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Ac);
	std::cout << std::endl;
	inv = CmpLinearAlgebra::CmpInverse(3, Ac, Bc);
	std::cout << "Ac inverse flag = " << inv << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Bc);
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
	std::cout << std::endl;
	std::cout << "Ac * Ac inverse" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
	std::complex<double> cDeterminant = 0;
	CmpLinearAlgebra::CmpGenerateNonSingular(
		2.0, cDeterminant, 4, 1, Mc);
	std::cout << std::endl;
	std::cout << "complex determinant = ";
	std::cout << cDeterminant << std::endl;
	double rDeterminantA = 0;
	std::vector<std::vector<double>> A44r(4,
		std::vector<double>(4));
	DblLinearAlgebra::DblGenerateNonSingular(
		2.0, rDeterminantA, 4, 2, A44r);
	double rDeterminantB = 0;
	std::vector<std::vector<double>> B44r(4,
		std::vector<double>(4));
	DblLinearAlgebra::DblGenerateNonSingular(
		2.0, rDeterminantB, 4, 3, B44r);
	std::cout << std::endl;
	std::vector<std::vector<double>> C44r(4,
		std::vector<double>(4));
	DblLinearAlgebra::DblMultiply(4, 4, 4, A44r, B44r, C44r);
	std::cout << "|A| = " << rDeterminantA << std::endl;
	std::cout << "|B| = " << rDeterminantB << std::endl;
	bool failure = false;
	double rDeterminantC =
		DblLinearAlgebra::DblDeterminant(4, 0, 0, C44r);
	std::cout << "|AB| = " << rDeterminantC << std::endl;
	std::cout << "|A||B| = " << rDeterminantA *
		rDeterminantB << std::endl;
	// Exercise 1.6 with 4x4 complex determinants
	std::vector<std::vector<std::complex<double>>> A44c(4,
		std::vector<std::complex<double>>(4));
	std::cout << std::endl;
	std::complex<double> cDeterminantA = 0;
	CmpLinearAlgebra::CmpGenerateNonSingular(
		2.0, cDeterminantA, 4, 2, A44c);
	std::vector<std::vector<std::complex<double>>> B44c(4,
		std::vector<std::complex<double>>(4));
	std::complex<double> cDeterminantB = 0;
	CmpLinearAlgebra::CmpGenerateNonSingular(
		2.0, cDeterminantB, 4, 3, B44c);
	std::vector<std::vector<std::complex<double>>> C44c(4,
		std::vector<std::complex<double>>(4));
	CmpLinearAlgebra::CmpMultiply(4, 4, 4, A44c, B44c, C44c);
	std::cout << "|A| = " << cDeterminantA << std::endl;
	std::cout << "|B| = " << cDeterminantB << std::endl;
	failure = false;
	std::complex<double> cDeterminantC =
		CmpLinearAlgebra::CmpDeterminant(4, C44c);
	std::cout << "|AB| = " << cDeterminantC << std::endl;
	std::cout << "|A||B| = " << cDeterminantA *
		cDeterminantB << std::endl;
	std::cout << "\nEnter any key to halt: ";
	char line[128] = "";
	std::cin.getline(line, 128);
}