To get SVD with GSL -- GNU Scientific Library

gsl_svd.cpp

#include <stdio.h>
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <sys/time.h>
#include <gsl/gsl_vector.h>
#include "gsl_SVD.h"


using namespace std;
using std::cout;
using std::endl;
using std::fixed;
using std::setprecision;


void linearSolve_SVD(const gsl_matrix *A, const gsl_vector *b, gsl_vector *x) {
    int rows = A -> size1;
    int cols = A -> size2;

    gsl_vector * work = gsl_vector_alloc(cols);
    gsl_vector * S = gsl_vector_alloc(cols);
    gsl_matrix * U = gsl_matrix_alloc(rows, cols);
    gsl_matrix * V = gsl_matrix_alloc(cols, cols);

    gsl_matrix_memcpy(U, A);  // 为了不破坏 A 中原始的数据，这里全都拷贝到 U 中

    gsl_linalg_SV_decomp(U, V, S, work);
    gsl_linalg_SV_solve(U, V, S, b, x);

    gsl_vector_free(work);
    gsl_vector_free(S);
    gsl_matrix_free(V);
    gsl_matrix_free(U);
}


int GslSVD::trimVectorS(double abseps) {
    int count = 0;
    for(int i = 0; i < S -> size; i++) {
        if(fabs(gsl_vector_get(S, i)) < abseps) {
            count ++;
            gsl_vector_set(S, i, 0);
        }
    }

    return count;
}


gsl_vector *GslSVD::getVectorS() {
    if(S == NULL) return NULL;
    gsl_vector *s = gsl_vector_alloc(S -> size);
    gsl_vector_memcpy(s, S);

    return s;
}


gsl_matrix *GslSVD::getMatrixU() {
    if(U == NULL) return NULL;
    gsl_matrix *u = gsl_matrix_alloc(U -> size1, U -> size2);
    gsl_matrix_memcpy(u, U);

    return u;
}


gsl_matrix *GslSVD::getMatrixV() {
    if(V == NULL) return NULL;
    gsl_matrix *v = gsl_matrix_alloc(V -> size1, V -> size2);
    gsl_matrix_memcpy(v, V);

    return v;
}


GslSVD::GslSVD() {
    S = NULL;
    U = NULL;
    V = NULL;
}


void GslSVD::alloc_suv(int rows, int cols) {
    if(S != NULL) {
        gsl_vector_free(S);
        gsl_matrix_free(U);
        gsl_matrix_free(V);
    }

    S = gsl_vector_alloc(cols);
    U = gsl_matrix_alloc(rows, cols);
    V = gsl_matrix_alloc(cols, cols);
}


int GslSVD::SV_decomp(const gsl_matrix *A) {
    int rows = A -> size1;
    int cols = A -> size2;

    gsl_vector *work = gsl_vector_alloc(cols);

    alloc_suv(rows, cols);
    gsl_matrix_memcpy(U, A);  // 为了不破坏 A 中原始的数据，这里全都拷贝到 U 中
    int ret = gsl_linalg_SV_decomp(U, V, S, work);

    gsl_vector_free(work);

    return ret;
}


int GslSVD::SV_decomp_mod(const gsl_matrix *A) {
    int rows = A -> size1;
    int cols = A -> size2;

    gsl_vector *work = gsl_vector_alloc(cols);
    gsl_matrix *X = gsl_matrix_alloc(cols, cols);

    alloc_suv(rows, cols);
    gsl_matrix_memcpy(U, A);  // 为了不破坏 A 中原始的数据，这里全都拷贝到 U 中
    int ret = gsl_linalg_SV_decomp_mod(U, X, V, S, work);

    gsl_matrix_free(X);
    gsl_vector_free(work);

    return ret;
}


int GslSVD::SV_decomp_jacobi(gsl_matrix *A) {
    int rows = A -> size1;
    int cols = A -> size2;
    alloc_suv(rows, cols);
    gsl_matrix_memcpy(U, A);  // 为了不破坏 A 中原始的数据，这里全都拷贝到 U 中
    int ret = gsl_linalg_SV_decomp_jacobi(U, V, S);

    return ret;
}


int GslSVD::SV_solve(const gsl_vector *b, gsl_vector *x) {
    if(U != NULL) {
        return gsl_linalg_SV_solve(U, V, S, b, x);
    }

    return -1;
}


GslSVD::~GslSVD() {
    if(S != NULL) {
        gsl_vector_free(S);
        gsl_matrix_free(V);
        gsl_matrix_free(U);
    }
}


// function declarations
// gsl_block read_csv(const string& fileName) {
    // FILE* inputfile;
    // fopen(fileName.c_str(), "r");  /* open file for reading */
    // if (inputfile == NULL)
    //    std::cout << "file does not exist n";

    // fseek(inputfile, 0L, SEEK_END);  // go until the end
    // int file_size = ftell(inputfile);  // in order to tell the size of the file
    // gsl_block* b = gsl_block_alloc(file_size);

    // if(inputfile)
    //    gsl_block_fscanf(inputfile, b);

    // fclose(inputfile);

    // return *b;
// }
// end function declarations


int main(void) {
    struct timeval start, end;

    // start timer
    gettimeofday(&start, NULL);

    // READ THE CSV INTO A BLOCK
    string filename = "./LBPhistogram_Binary_20220609.csv";
    // gsl_block data;
    // data = read_csv(filename);
    // use gsl_vector_alloc_from_block() but how?

    int rows = 5000;
    int cols = 131072;

    // Read csv into matrix
    gsl_matrix *data = gsl_matrix_alloc(rows, cols);  // specify input data dimensions

    FILE *f = fopen(filename.c_str(), "rb");
    gsl_matrix_fread(f, data);
    fclose(f);

    // gsl_matrix_view A = gsl_matrix_alloc_from_block(data, rows, cols);
    // gsl_block_free(gsl_block * data);  // 通过指针，内存释放
    // gsl_matrix_view A = data;

    /*
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            double dij = gsl_matrix_get(data, i, j);
            cout << dij << endl;
        }
    }
    */
    // printf("rows=%ld, cols=%ld\n", data->size1, data->size2);

    // gsl_matrix_free(data);

    // double a_data[] = {1, 2, 2, 4};
    // gsl_matrix_view A = gsl_matrix_view_array(a_data, 2, 2);
    GslSVD svd;
    // svd.SV_decomp(&A.matrix);
    svd.SV_decomp(data);

    gsl_matrix_free(data);

    puts("S = \n");
    gsl_vector_fprintf(stdout, svd.getVectorS(), "%f");

    // puts("\nU = ");
    // gsl_matrix_fprintf(stdout, svd.getMatrixU(), "%f");

    // puts("\nV = ");
    // gsl_matrix_fprintf(stdout, svd.getMatrixV(), "%f");

    // double b_data[] = {3, 6};
    // gsl_vector_view b = gsl_vector_view_array(b_data, 2);
    // gsl_vector *x = gsl_vector_alloc(2);
    // svd.SV_solve(&b.vector, x);

    // puts("\nx = ");
    // gsl_vector_fprintf(stdout, x, "%f");

    // stop timer.
    gettimeofday(&end, NULL);

    // Calculating total time taken by the program.
    double time_taken;
    time_taken = (end.tv_sec - start.tv_sec) * 1e6;
    time_taken = (time_taken + (end.tv_usec - start.tv_usec)) * 1e-6;
    cout << "\nTime taken by program is: " << fixed << setprecision(6) << time_taken << " seconds." << endl;

    return 0;
}