From 797d3d34b7636a49b39f1f05da080fee990ad9b4 Mon Sep 17 00:00:00 2001 From: GillesDuvert Date: Tue, 19 Dec 2023 23:35:08 +0100 Subject: [PATCH] Sparse Matrix Support, minimal double version. --- src/CMakeLists.txt | 1 + src/libinit_ac.cpp | 10 +++ src/sparse_matrix.cpp | 153 ++++++++++++++++++++++++++++++++++++++++++ src/sparse_matrix.hpp | 32 +++++++++ 4 files changed, 196 insertions(+) create mode 100644 src/sparse_matrix.cpp create mode 100644 src/sparse_matrix.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 968ead9d9..fba492a6a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -155,6 +155,7 @@ saverestore.cpp semshm.cpp sigfpehandler.cpp sorting.cpp +sparse_matrix.cpp str.cpp terminfo.cpp tiff.cxx diff --git a/src/libinit_ac.cpp b/src/libinit_ac.cpp index 416c86544..8882e1b11 100644 --- a/src/libinit_ac.cpp +++ b/src/libinit_ac.cpp @@ -24,6 +24,7 @@ #include "math_fun_ac.hpp" #include "gsl_matrix.hpp" +#include "sparse_matrix.hpp" #include "gsl_fun.hpp" #include "smooth.hpp" @@ -74,6 +75,15 @@ void LibInit_ac() new DLibFunRetNew(lib::fx_root_fun,string("FX_ROOT"),2,fx_rootKey); #endif + + //sparse matrices by GD + const string sprsinKey[]={"COLUMN", "DOUBLE", "THRESHOLD", KLISTEND}; + new DLibFunRetNew(lib::sprsin_fun,string("SPRSIN"),4,sprsinKey); + const string sprsabKey[]={"DOUBLE", "THRESHOLD",KLISTEND}; + new DLibFunRetNew(lib::sprsab_fun,string("SPRSAB"),2,sprsabKey); + const string sprsaxKey[]={"DOUBLE",KLISTEND}; + new DLibFunRetNew(lib::sprsax_fun,string("SPRSAX"),2,sprsaxKey); + new DLibFunRetNew(lib::fulstr_fun,string("FULSTR"),1); const string spl1Key[]={"YP0","YPN_1","YP1","DOUBLE","HELP",KLISTEND}; //YP1 is old value for YP0 new DLibFunRetNew(lib::spl_init_fun,string("SPL_INIT"),2,spl1Key); diff --git a/src/sparse_matrix.cpp b/src/sparse_matrix.cpp new file mode 100644 index 000000000..560b1d454 --- /dev/null +++ b/src/sparse_matrix.cpp @@ -0,0 +1,153 @@ +/*************************************************************************** + sparse_matrix.cpp - GDL sparse matrix functions + ------------------- + begin : Dec 9 2023 + copyright : (C) 2023 by Gilles Duvert + email : surname dot name at free dot fr + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include +#include +#include +#include "sparse_matrix.hpp" + +typedef Eigen::SparseMatrix SPMATRowMajDbl; +typedef Eigen::SparseVector SPVecDbl; +static const float f_thr=1E-7; +static const double d_thr=1E-14; +namespace lib { + + SPMATRowMajDbl getFromStruct(EnvT* e, DUInt i) { + BaseGDL* p0 = e->GetParDefined(i); // must be struct + DType varType = p0->Type(); + if (varType != GDL_STRUCT) e->Throw("Expression " + e->GetString(i) + "must be a structure in this context."); + DStructGDL* s = e->GetParAs(i); + DLongGDL* N = static_cast (s->GetTag(0)); + int nCols = (*N)[0]; + int nRows = (*N)[1]; + int nnz = (*static_cast (s->GetTag(1)))[0]; + if (nnz > 0 ) { //protect against 0 + double* values = static_cast (s->GetTag(2)->DataAddr()); + int* innerIndicesPtr = static_cast (s->GetTag(3)->DataAddr()); + int* outerIndexPtr = static_cast (s->GetTag(4)->DataAddr()); + Eigen::Map Mat(nRows, nCols, nnz, outerIndexPtr, innerIndicesPtr, values); + return Mat; + } else { + SPMATRowMajDbl Mat(nRows, nCols); + return Mat; + } + } + BaseGDL* convertToGDL(SPMATRowMajDbl Mat) { + int ncols = Mat.cols(); + int nrows = Mat.rows(); + DDoubleGDL* ret = new DDoubleGDL(dimension(ncols, nrows), BaseGDL::ZERO); + const int outs=Mat.outerSize(); + const int* outerStartIndexPtr = Mat.outerIndexPtr(); + const int* innerIndicesPtr = Mat.innerIndexPtr(); + const double* values = Mat.valuePtr(); + for (auto iRow = 0; iRow < outs; ++iRow) { + int jValmin = outerStartIndexPtr[iRow]; + int jValmax = outerStartIndexPtr[iRow + 1]; + for (int kk = jValmin; kk < jValmax; ++kk) { //outerstart + int iCol = innerIndicesPtr[kk]; //row + (*ret)[iRow * ncols + iCol] = values[kk]; + } + } + return ret; + } + + DStructGDL* convertToStruct(const SPMATRowMajDbl Mat) { + int nCols = Mat.cols(); + int nRows = Mat.rows(); + DStructDesc* sd = new DStructDesc("$truct"); + DStructGDL* s = new DStructGDL(sd, dimension(1)); + DLongGDL* Dim = new DLongGDL(dimension(2), BaseGDL::NOZERO); + (*Dim)[0] = nCols; + (*Dim)[1] = nRows; + s->NewTag("DIM", Dim); + DLong nnz = Mat.nonZeros(); + s->NewTag("NNZ", new DLongGDL(nnz)); + //protect against 0 + if (nnz > 0) { + DDoubleGDL* Values = new DDoubleGDL(dimension(nnz), BaseGDL::NOZERO); + for (auto i = 0; i < nnz; ++i) (*Values)[i] = Mat.valuePtr()[i]; + s->NewTag("VALUES", Values); + DLongGDL* InnerIndices = new DLongGDL(dimension(nnz), BaseGDL::NOZERO); + for (auto i = 0; i < nnz; ++i) (*InnerIndices)[i] = Mat.innerIndexPtr()[i]; + s->NewTag("INNER_INDICES", InnerIndices); + int outs=Mat.outerSize(); + DLongGDL* OuterStarts = new DLongGDL(dimension(outs+1), BaseGDL::NOZERO); + for (auto i = 0; i < outs+1; ++i) (*OuterStarts)[i] = Mat.outerIndexPtr()[i]; + s->NewTag("OUTER_STARTS", OuterStarts); + } + return s; + } + + BaseGDL* sprsin_fun(EnvT* e) { + SizeT nParam = e->NParam(); //1 or 4 + if (nParam != 1 && nParam != 4) e->Throw("Incorrect number of arguments."); + double thresh=d_thr; + static int THRESHOLD=e->KeywordIx("THRESHOLD"); + if (e->KeywordSet(THRESHOLD)){ + e->AssureDoubleScalarKW(THRESHOLD,thresh); + } + if (nParam == 1) { + BaseGDL* p0 = e->GetParDefined(0); // matrix + DType varType = p0->Type(); + if (p0->Dim().Rank() != 2) e->Throw("Argument " + e->GetString(0) + " must be a square matrix."); + if (varType == GDL_STRING) e->Throw("Argument " + e->GetString(0) + " must not be of STRING type."); + int nCols = p0->Dim(0); + int nRows = p0->Dim(1); + DDoubleGDL* var = e->GetParAs(0); + SPMATRowMajDbl Mat(nRows,nCols); + //import and prune wrt. threshold + for (auto i=0; i=thresh) Mat.insert(j,i)= (*var)[j*nCols+i]; + Mat.makeCompressed(); + return convertToStruct(Mat); + } + return new DLongGDL(0); + } + + BaseGDL* sprsab_fun(EnvT* e){ + SizeT nParam = e->NParam(2); + double thresh=d_thr; + static int THRESHOLD=e->KeywordIx("THRESHOLD"); + if (e->KeywordSet(THRESHOLD)){ + e->AssureDoubleScalarKW(THRESHOLD,thresh); + } + SPMATRowMajDbl Mat1=getFromStruct(e, 0); + SPMATRowMajDbl Mat2=getFromStruct(e, 1); + SPMATRowMajDbl Mat3=(Mat1*Mat2).pruned(thresh); + return convertToStruct(Mat3); + } + BaseGDL* sprsax_fun(EnvT* e) { + SizeT nParam = e->NParam(2); + SPMATRowMajDbl Mat1 = getFromStruct(e, 0); + BaseGDL* p1 = e->GetParDefined(1); // vector + DType varType = p1->Type(); + if (p1->Dim().Rank() != 1) e->Throw("Argument " + e->GetString(1) + " must be a vector."); + if (varType == GDL_STRING) e->Throw("Argument " + e->GetString(1) + " must not be of STRING type."); + int m = p1->Dim(0); + if (m != Mat1.cols()) e->Throw("Argument " + e->GetString(1) + " does not have correct size."); + DDoubleGDL* var = e->GetParAs(1); + SPVecDbl Vec(m); + //import + for (auto i=0; iNParam(1); + SPMATRowMajDbl Mat=getFromStruct(e, 0); + return convertToGDL(Mat); + } +} diff --git a/src/sparse_matrix.hpp b/src/sparse_matrix.hpp new file mode 100644 index 000000000..7c1c59779 --- /dev/null +++ b/src/sparse_matrix.hpp @@ -0,0 +1,32 @@ +/*************************************************************************** + sparse_matrix.hpp - GDL sparse matrix functions + ------------------- + begin : Dec 9 2023 + copyright : (C) 2023 by Gilles Duvert + email : surname dot name at free dot fr + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +#ifndef SPARSE_MATRIX_HPP_ +#define SPARSE_MATRIX_HPP_ + +#include "includefirst.hpp" +#include "datatypes.hpp" +#include "envt.hpp" + +namespace lib { + + BaseGDL* sprsin_fun(EnvT* e); + BaseGDL* sprsax_fun(EnvT* e); + BaseGDL* sprsab_fun(EnvT* e); + BaseGDL* fulstr_fun(EnvT* e); + +} +#endif