Project-2/src/main.cpp

/** @file main.cpp
 *  @brief Main program for Project 2
 *
 *  The program performs the Jacobi rotation method. 
 *  The size of the matrix, and number of transformations
 *  performed are written to file. 
 *  Eigenvector corresponding to the 3 smallest eigenvalues 
 *  for matrices of size 6x6 and 100x100 are written to file.
 *
 *  @author Cory Alexander Balaton (coryab)
 *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
 *  @bug No known bugs
 */
#include <cassert>
#include <cmath>
#include <ctime>
#include <iostream>
#include <omp.h>
#include <ostream>

#include "utils.hpp"
#include "matrix.hpp"
#include "jacobi.hpp"


void write_transformation_dense(int N)
{   
    std::ofstream ofile;

    ofile.open("../latex/output/transform_dense.csv");
    ofile << "N,T" << std::endl;

    // Increase size of matrix, start at 5 to avoid logic_error of N=4
    #pragma omp parallel for ordered schedule(static, 1)
    for (int i = 5; i <= N; i++) {
        arma::mat A = arma::mat(i, i).randn();
        A = arma::symmatu(A);
        arma::vec eigval;
        arma::mat eigvec;

        int iters;
        bool converged;
        
        jacobi_eigensolver(A, 10e-14, eigval, eigvec, 100000, iters, converged);

        // Write size, and number of iterations to file
        #pragma omp ordered
        {
            ofile << i << "," << iters << std::endl;
        }
    }
    ofile.close();
}

    
void write_transformation_tridiag(int N)
{   
    std::ofstream ofile;

    double h;
    double a, d;

    ofile.open("../latex/output/transform_tridiag.csv");
    // Write header line to file 
    ofile << "N,T" << std::endl;

    // Increase size of matrix, start at 5 to avoid logic_error of N=4
    #pragma omp parallel for ordered schedule(static, 1)
    for (int i = 5; i <= N; i++) {
        h = 1. / (double) (i+1);
        a = -1. / (h*h), d = 2. / (h*h);

        arma::mat A = create_symmetric_tridiagonal(i, a, d);
        arma::vec eigval;
        arma::mat eigvec;

        int iters;
        bool converged;
        
        jacobi_eigensolver(A, 10e-14, eigval, eigvec, 100000, iters, converged);

        // Write size, and number of iterations to file
        #pragma omp ordered
        {
            ofile << i << "," << iters << std::endl;
        }
    }
    ofile.close();
}


void write_eigenvec(int N)
{
    double h = 1. / (double) (N+1);
    double a = -1. / (h*h);
    double d = 2. / (h*h);
    double x = 0.;

    // Create tridiagonal matrix
    arma::mat A = create_symmetric_tridiagonal(N, a, d);
    arma::mat analytic = arma::mat(N, N);
    arma::vec eigval;
    arma::mat eigvec;

    int iters;
    bool converged;
    
    // Solve using Jacobi rotation method
    jacobi_eigensolver(A, 10e-14, eigval, eigvec, 100000, iters, converged);

    // Build analytic eigenvectors 
    arma::vec v, analytic_vec = arma::vec(N);
    for (int i=0; i < N; i++) {
        v = eigvec.col(i);

        for (int j=0; j < N; j++) {
            analytic_vec(j) = std::sin(((j+1.)*(i+1.)*M_PI) / (N+1.));
        }

        analytic_vec = arma::normalise(analytic_vec);

        // Flip the sign of the analytic vector if they are different
        if (analytic_vec(0)*v(0) < 0.) {
            analytic_vec *= -1;
        }
        analytic.col(i) = analytic_vec;
    }

    std::ofstream ofile;
    // Create file based on matrix size, and write header line to file 
    ofile.open("../latex/output/eigenvector_" + std::to_string(N) + ".csv");
    ofile << "x," 
        << "Vector 1,Vector 2,Vector 3," 
        << "Analytic 1,Analytic 2,Analytic 3" << std::endl;

    // Add boundary value for x=0 
    ofile << scientific_format(0., 16) << "," 
            << scientific_format(0., 16) << "," 
            << scientific_format(0., 16) << "," 
            << scientific_format(0., 16) << "," 
            << scientific_format(0., 16) << "," 
            << scientific_format(0., 16) << ","
            << scientific_format(0., 16) << std::endl;

    // Add x-value and element i of each eigenvector to same line
    for (int i = 0; i < N; i++) {
        x += h;
        ofile << scientific_format(x, 16)<< "," 
            << scientific_format(eigvec(i,0), 16) << "," 
            << scientific_format(eigvec(i,1), 16) << "," 
            << scientific_format(eigvec(i,2), 16) << ","
            << scientific_format(analytic(i,0), 16) << ","
            << scientific_format(analytic(i,1), 16) << ","
            << scientific_format(analytic(i,2), 16) << std::endl;
    }
    // Add boundary value for x=1
    ofile << scientific_format(1., 16) << "," 
            << scientific_format(0., 16) << "," 
            << scientific_format(0., 16) << "," 
            << scientific_format(0., 16) << ","
            << scientific_format(0., 16) << ","
            << scientific_format(0., 16) << ","
            << scientific_format(0., 16) << std::endl;
    ofile.close();
}


int main()
{
    write_transformation_tridiag(100);
    write_transformation_dense(100);
    write_eigenvec(10);
    write_eigenvec(100);
    return 0;
}