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..604c3dd90 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,20 @@ void LibInit_ac() new DLibFunRetNew(lib::fx_root_fun,string("FX_ROOT"),2,fx_rootKey); #endif + + //sparse matrices by GD + const string sprsinKey[]={"THRESHOLD", "COLUMN", KLISTEND}; + const string sprsinIgnoreKey[]={"DOUBLE", KLISTEND}; + new DLibFunRetNew(lib::sprsin_fun,string("SPRSIN"),4,sprsinKey,sprsinIgnoreKey); + 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); + new DLibFunRetNew(lib::sprstp_fun,string("SPRSTP"),1); + const string linbcgKey[]={"ITER", "DOUBLE", KLISTEND}; + const string linbcgIgnoreKey[]={"ITOL", "TOL", "ITMAX", KLISTEND}; + new DLibFunRetNew(lib::linbcg_fun,string("LINBCG"),3,linbcgKey,linbcgIgnoreKey); 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..d32022f2c --- /dev/null +++ b/src/sparse_matrix.cpp @@ -0,0 +1,246 @@ +/*************************************************************************** + 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; + +//To save time, the 'external' (i.e., as a GDL variable) representation of a Eigen::SPMatrix is just a DPtr to the Eigen::SPMatrix itself +//I have not checked if the reference counting of Ptrs is sufficient to destroy the SPMatrix when the Ptr is destroyed. +//Otherwise this code may lead to leaks. +typedef std::vector::iterator MatVecIterator; +namespace lib { + + SPMATRowMajDbl* getFromPtr(EnvT* e, DUInt i) { + BaseGDL* p0 = e->GetParDefined(i); // must be struct + DType varType = p0->Type(); + if (varType != GDL_PTR) e->Throw("Expression " + e->GetString(i) + "must be a PTR in this context."); + DPtrGDL* s = e->GetParAs(i); + return *((SPMATRowMajDbl**)(s->DataAddr())); + } + + //Version where the variable exchanged is a special structure. useful for save/restore, but slows down process in all other cases + 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; + } + } + //takes a SparseMatrix and returns the full GDL variable. + 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; + } + + //Version where the variable exchanged is a special structure. useful for save/restore, but slows down process in all other cases + 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; + } + + DPtrGDL* convertToPtr(const SPMATRowMajDbl *Mat) { + DPtr pointer = (DPtr)Mat; //(DPtr)(MatVec.data()); + return new DPtrGDL(pointer); + } + + BaseGDL* sprsin_fun(EnvT* e) { + static bool warned=false; + SizeT nParam = e->NParam(); //1 or 4 + if (nParam != 1 && nParam != 4) e->Throw("Incorrect number of arguments."); + double thresh=d_thr; + static int COLUMN = e->KeywordIx("COLUMN"); + if (e->KeywordSet(COLUMN)) { + e->Throw("GDL's SPRSIN does not yet support the COLUMN keyword. Please report or use transposed arrays."); + } + 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); + if (nCols != nRows && !warned) { + Message("NOTE: use of SPRSIN with non-square matrices is a GDL extension."); + warned=true; + } + DDoubleGDL* var = e->GetParAs(0); + SPMATRowMajDbl *Mat= new SPMATRowMajDbl(nRows,nCols); + //import and prune wrt. threshold + for (auto i=0; i=thresh) Mat->insert(j,i)= (*var)[j*nCols+i]; + Mat->makeCompressed(); + return convertToPtr(Mat); + } else if (nParam == 4) { + DLongGDL* cols = e->GetParAs(0); + int nCols=cols->N_Elements(); + DLongGDL* rows = e->GetParAs(1); + int nRows=rows->N_Elements(); + if (nCols != nRows) e->Throw("Vector "+e->GetString(1) + " must have "+ i2s(nCols) + " elements."); + DDoubleGDL* vals = e->GetParAs(2); + int nVals=vals->N_Elements(); + if (nVals != nRows) e->Throw("Vector "+e->GetString(2) + " must have "+ i2s(nCols) + " elements."); + DLongGDL* sizegdl = e->GetParAs(3); + DLong size=(*sizegdl)[0]; + SPMATRowMajDbl *Mat= new SPMATRowMajDbl(size,size); //only square matrices this way. + Mat->reserve(nCols); + for (auto i=0; i< nCols; ++i) { + DLong irow=(*rows)[i]; + if (irow >= size ) e->Throw(" Out of range subscript encountered: "+e->GetString(0)+" ."); + DLong icol=(*cols)[i]; + if (icol >= size ) e->Throw(" Out of range subscript encountered: "+e->GetString(1)+" ."); + Mat->coeffRef(irow,icol)+=(*vals)[i]; //protect against doublons + if (Mat->coeffRef(irow,icol)!=(*vals)[i]) e->Throw("Duplicate subscript encountered in Columns and Rows arrays at element: "+i2s(i)); + } + return convertToPtr(Mat); + } else e->Throw("Incorrect number of arguments."); + + 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=getFromPtr(e, 0); + SPMATRowMajDbl *Mat2=getFromPtr(e, 1); + SPMATRowMajDbl *Mat3=new SPMATRowMajDbl((*Mat1)*(*Mat2)); + Mat3->prune(thresh); + return convertToPtr(Mat3); + } + + BaseGDL* sprstp_fun(EnvT* e){ + SizeT nParam = e->NParam(1); + SPMATRowMajDbl *Mat=getFromPtr(e, 0); + SPMATRowMajDbl *res=new SPMATRowMajDbl((*Mat).transpose()); + return convertToPtr(res); + } + BaseGDL* sprsax_fun(EnvT* e) { + SizeT nParam = e->NParam(2); + SPMATRowMajDbl* Mat = getFromPtr(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 != Mat->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(3); + SPMATRowMajDbl* A = getFromPtr(e, 0); + BaseGDL* p1 = e->GetParDefined(1); // Right Hand B + 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 != A->cols()) e->Throw("Argument " + e->GetString(1) + " does not have correct size."); + DDoubleGDL* Bgdl = e->GetParAs(1); + Eigen::Map B(&(*Bgdl)[0], m); + BaseGDL* p2 = e->GetParDefined(2); // Initial Guess + varType = p2->Type(); + if (p2->Dim().Rank() != 1) e->Throw("Argument " + e->GetString(2) + " must be a vector."); + if (varType == GDL_STRING) e->Throw("Argument " + e->GetString(2) + " must not be of STRING type."); + int n = p2->Dim(0); + if (n != m) e->Throw("Argument " + e->GetString(2) + " does not have correct size."); + DDoubleGDL* Xgdl = e->GetParAs(2); + Eigen::Map X(&(*Xgdl)[0], n); + //solve ax=b + Eigen::SparseLU> sparseLU; + sparseLU.compute(*A); + if(sparseLU.info()!=Eigen::Success) e->Throw("Matrix decomposition failed."); + X=sparseLU.solve(B); + if(sparseLU.info()!=Eigen::Success) e->Throw("No solution found."); + DDoubleGDL* resD =new DDoubleGDL(n, BaseGDL::NOZERO); + Eigen::Map(&(*resD)[0], n) = X; + return resD; + } + + BaseGDL* fulstr_fun(EnvT* e){ + SizeT nParam = e->NParam(1); +// SPMATRowMajDbl Mat=getFromIndexInMatVec(e, 0); + SPMATRowMajDbl *Mat=getFromPtr(e, 0); + return convertToGDL(Mat); + } +} diff --git a/src/sparse_matrix.hpp b/src/sparse_matrix.hpp new file mode 100644 index 000000000..cb648463d --- /dev/null +++ b/src/sparse_matrix.hpp @@ -0,0 +1,33 @@ +/*************************************************************************** + 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* sprstp_fun(EnvT* e); + BaseGDL* fulstr_fun(EnvT* e); + BaseGDL* linbcg_fun(EnvT* e); +} +#endif diff --git a/testsuite/LIST b/testsuite/LIST index f0f309539..3b6b821d0 100644 --- a/testsuite/LIST +++ b/testsuite/LIST @@ -167,6 +167,7 @@ test_scope_varname.pro test_simplex.pro test_size.pro test_sort.pro +test_sparse.pro test_spher_harm.pro test_spl.pro test_standardize.pro diff --git a/testsuite/test_sparse.pro b/testsuite/test_sparse.pro new file mode 100644 index 000000000..b76e81186 --- /dev/null +++ b/testsuite/test_sparse.pro @@ -0,0 +1,20 @@ +; no a test, really. Just a suggestion for test, plu sthe fact +; that it exercises the coverage. +pro test_sparse + + a = [[ 2.0, 1.0, 1.0],[ 4.0, -6.0, 0.0],[-2.0, 7.0, 2.0]] + z=sprsin(a, thresh=0.5) + zz=sprstp(z) + q=fulstr(z) + if abs(total(a-q)) gt 1E-6 then exit,status=1 + res=sprsab(z,zz) + result=fulstr(res) + if total(result) ne 29 then exit,status=1 + if total(sprsax(z,[1,1,1])) ne 9 then exit,status=1 + aludc=a + LUDC, aludc, index + b = [3,-8.0,10] + x = LUSOL(aludc, index, b) + r= LINBCG(SPRSIN(a), B, X) + if abs(total(r-x)) gt 1E-6 then exit,status=1 +end