Temporary Cpp Code

#define USE_EXISTING_YS true
#define MIN_VALUE -100
#define MAX_VALUE +100

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string>

#include <chrono>
#include <random>

#include "boost\filesystem.hpp"
#include "boost\regex.hpp"
#include "boost\algorithm\string\predicate.hpp"

#include "mmio.c"

double* cholesky(double *A, int n)
{
    double *L = (double*)calloc(n * n, sizeof(double));
    if (L == NULL) {
        exit(EXIT_FAILURE);
    }

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < (i + 1); j++) {
            double s = 0;
            for (int k = 0; k < j; k++)
                s += L[i * n + k] * L[j * n + k];
            L[i * n + j] = (i == j) ?
                sqrt(A[i * n + i] - s) :
                (1.0 / L[j * n + j] * (A[i * n + j] - s));
        }
    }
    return L;
}

void show_matrix(double *A, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++)
            printf("%2.5f ", A[i * n + j]);
        printf("\n");
    }
}

void show_vector(double *v, int n) {
    for (int i = 0; i < n; i++)
        printf("%2.5f ", v[i]);

    printf("\n");
}

double* forwardSubstitution(double *L, double *b, int n) {

    double *y = (double*)calloc(n, sizeof(double));
    if (y == NULL)
        exit(EXIT_FAILURE);

    for (int i = 0; i < n; i++) {
        y[i] = b[i];

        for (int j = 0; j <= i - 1; j++)
            y[i] -= L[i + j * n] * y[j];

        y[i] /= L[i + i * n];
    }
    return y;
}

double* backwardSubstitution(double *U, double *b, int n, bool isTransposed = true) {

    double *x = (double*)calloc(n, sizeof(double));
    if (x == NULL)
        exit(EXIT_FAILURE);

    // Same codes but with i and j swapped for matrix U
    if (isTransposed) {
        for (int i = n - 1; i >= 0; i--) {
            x[i] = b[i];

            for (int j = i + 1; j < n; j++)
                x[i] -= U[i + j * n] * x[j];

            x[i] /= U[i + i * n];
        }
    }
    else {
        for (int i = n - 1; i >= 0; i--) {
            x[i] = b[i];

            for (int j = i + 1; j < n; j++)
                x[i] -= U[j + i * n] * x[j];

            x[i] /= U[i + i * n];
        }
    }
    return x;
}

int main()
{
    using namespace std;
    using namespace boost;
    using namespace boost::filesystem;
    using namespace boost::algorithm;
    using namespace std::chrono;

    time_point<high_resolution_clock> start, finish;

    // reusable time variable
    long long time;

    // Random number generator for Ys
    random_device rd; // Random at each execution
    default_random_engine dre(rd());

    // Get the files with extension .mtx (specify full or relative path)
    path current_dir("..\\SparseMatrixProject\\matrices\\small");
    regex pattern("(.*\\.mtx)");

    // Iterate and read matrices
    for (recursive_directory_iterator iter(current_dir), end; iter != end; ++iter)
    {
        string name = iter->path().filename().string();

        if (regex_match(name, pattern))
        {
            // Or if it ends with 'b.mtx'!
            if (ends_with(name, "b.mtx") || ends_with(name, "coord.mtx")) {
                continue;
            }

            // Ignore this matrix
            if (starts_with(name, "bibd")) {
                continue;
            }

            FILE *f;
            int ret_code;
            MM_typecode matcode;
            int M, N, nz;
            int i, *I, *J;
            double *val, *inlineMat;

            string filePath = iter->path().string();

            // Open file
            if ((f = fopen(filePath.c_str(), "r")) == NULL)
                exit(1);

            if (mm_read_banner(f, &matcode) != 0)
            {
                printf("Could not process Matrix Market banner.\n");
                exit(1);
            }

            /*  This is how one can screen matrix types if their application */
            /*  only supports a subset of the Matrix Market data types.      */
            /*if (mm_is_complex(matcode) && mm_is_matrix(matcode) &&
                mm_is_sparse(matcode))
            {
                printf("Sorry, this application does not support ");
                printf("Market Market type: [%s]\n", mm_typecode_to_str(matcode));
                exit(1);
            }*/

            /* find out size of sparse matrix .... */
            if ((ret_code = mm_read_mtx_crd_size(f, &M, &N, &nz)) != 0)
                exit(1);

            // Ys
            double *b;
            b = (double *)malloc(sizeof(double)*N);
            //x = (double *)calloc(N, sizeof(double)); // allocates memory and fills with zeroes.

            if (USE_EXISTING_YS)
            {
                string ysPath = filePath.substr(0, filePath.find(".mtx", 0)) + "_Ys.txt";

                // Open Ys file
                FILE *ysFile;
                if ((ysFile = fopen(ysPath.c_str(), "r")) == NULL)
                    exit(1);

                for (size_t i = 0; i < N; i++)
                    fscanf(ysFile, "%lg\n", &b[i]);

                fclose(ysFile);
            }
            else {
                uniform_real_distribution<double> di(MIN_VALUE, MAX_VALUE);
                for (size_t i = 0; i < N; i++)
                    b[i] = di(dre);
            }

            /* reseve memory for matrices */
            I = (int *)malloc(nz * sizeof(int));
            J = (int *)malloc(nz * sizeof(int));
            val = (double *)malloc(nz * sizeof(double));
            inlineMat = (double *)malloc(N * N * sizeof(double));

            /* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
            /*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
            /*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */
            for (i = 0; i < nz; i++)
            {
                fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
                I[i]--;  /* adjust from 1-based to 0-based */
                J[i]--;

                inlineMat[I[i] + J[i] * N] = val[i];
            }

            if (f != stdin) fclose(f);

            /************************/
            /* now write out matrix */
            /************************/
            mm_write_banner(stdout, matcode);
            mm_write_mtx_crd_size(stdout, M, N, nz);
            /*for (i = 0; i < nz; i++) {
                fprintf(stdout, "%d %d %20.19g\n", I[i] + 1, J[i] + 1, val[i]);
            }*/

            // [1] Solve with "Cholevsky Decomposition"
            start = high_resolution_clock::now();

            // TODO compute cholevsky

            finish = high_resolution_clock::now();
            time = duration_cast<milliseconds>(finish - start).count();
            fprintf(stdout, "Cholevsky: %d milliseconds. \n", time);

            // [2] Solve with "Conjugate Gradient"
            start = high_resolution_clock::now();

            // TODO compute cg

            finish = high_resolution_clock::now();
            time = duration_cast<milliseconds>(finish - start).count();
            fprintf(stdout, "CG: %d milliseconds.\n", time);

            free(I);
            free(J);
            free(val);

            free(b);

            printf("\n");
        }
    }
}