From e39259fbbd5918405ca65963f3334f42bd30a790 Mon Sep 17 00:00:00 2001 From: GillesDuvert Date: Wed, 20 Dec 2023 16:16:01 +0100 Subject: [PATCH] use of Ptr instead of slow strutures. Finished 90% work. --- src/libinit_ac.cpp | 5 ++- src/sparse_matrix.cpp | 99 ++++++++++++++++++++++++++++++++----------- 2 files changed, 77 insertions(+), 27 deletions(-) diff --git a/src/libinit_ac.cpp b/src/libinit_ac.cpp index 8882e1b11..cad5c3b4c 100644 --- a/src/libinit_ac.cpp +++ b/src/libinit_ac.cpp @@ -77,8 +77,9 @@ void LibInit_ac() #endif //sparse matrices by GD - const string sprsinKey[]={"COLUMN", "DOUBLE", "THRESHOLD", KLISTEND}; - new DLibFunRetNew(lib::sprsin_fun,string("SPRSIN"),4,sprsinKey); + 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}; diff --git a/src/sparse_matrix.cpp b/src/sparse_matrix.cpp index 560b1d454..a9f734f3d 100644 --- a/src/sparse_matrix.cpp +++ b/src/sparse_matrix.cpp @@ -24,10 +24,24 @@ typedef Eigen::SparseMatrix SPMATRowMajDbl; typedef Eigen::SparseVector SPVecDbl; static const float f_thr=1E-7; static const double d_thr=1E-14; -namespace lib { +//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 + 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); @@ -46,14 +60,15 @@ namespace lib { return Mat; } } - BaseGDL* convertToGDL(SPMATRowMajDbl Mat) { - int ncols = Mat.cols(); - int nrows = Mat.rows(); + //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(); + 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]; @@ -63,8 +78,9 @@ namespace lib { } } 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(); @@ -91,11 +107,21 @@ namespace lib { } 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); @@ -107,13 +133,32 @@ namespace lib { 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(nRows,nCols); + 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 convertToStruct(Mat); - } + 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* size = e->GetParAs(3); + SPMATRowMajDbl *Mat= new SPMATRowMajDbl((*size)[0],(*size)[0]); //only square matrices this way. + Mat->reserve(nCols); + for (auto i=0; i< nCols; ++i) Mat->insert((*rows)[i],(*cols)[i])=(*vals)[i]; + return convertToPtr(Mat); + } else e->Throw("Incorrect number of arguments."); + return new DLongGDL(0); } @@ -124,30 +169,34 @@ namespace lib { 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); + SPMATRowMajDbl *Mat1=getFromPtr(e, 0); + SPMATRowMajDbl *Mat2=getFromPtr(e, 1); + SPMATRowMajDbl *Mat3=new SPMATRowMajDbl((*Mat1)*(*Mat2)); + Mat3->prune(thresh); + return convertToPtr(Mat3); } + BaseGDL* sprsax_fun(EnvT* e) { SizeT nParam = e->NParam(2); - SPMATRowMajDbl Mat1 = getFromStruct(e, 0); + 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 != Mat1.cols()) e->Throw("Argument " + e->GetString(1) + " does not have correct size."); + 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(1); - SPMATRowMajDbl Mat=getFromStruct(e, 0); +// SPMATRowMajDbl Mat=getFromIndexInMatVec(e, 0); + SPMATRowMajDbl *Mat=getFromPtr(e, 0); return convertToGDL(Mat); } }