Skip to content

Commit

Permalink
use of Ptr instead of slow strutures.
Browse files Browse the repository at this point in the history 
Finished 90% work.
  • Loading branch information
@GillesDuvert
GillesDuvert committed Dec 20, 2023
1 parent 797d3d3 commit e39259f
Show file tree
Hide file tree
Showing 2 changed files with 77 additions and 27 deletions.
5 changes: 3 additions & 2 deletions src/libinit_ac.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down
99 changes: 74 additions & 25 deletions src/sparse_matrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,10 +24,24 @@ typedef Eigen::SparseMatrix<double, Eigen::RowMajor> SPMATRowMajDbl;
typedef Eigen::SparseVector<double> 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<SPMATRowMajDbl>::iterator MatVecIterator;
namespace lib {

SPMATRowMajDbl* getFromPtr(EnvT* e, DUInt i) {
BaseGDL* p0 = e->GetParDefined(i); // must be struct
DType varType = p0->Type();

Check warning on line 36 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L34-L36 

Added lines #L34 - L36 were not covered by tests
if (varType != GDL_PTR) e->Throw("Expression " + e->GetString(i) + "must be a PTR in this context.");
DPtrGDL* s = e->GetParAs<DPtrGDL>(i);
return *((SPMATRowMajDbl**)(s->DataAddr()));

Check warning on line 39 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L38-L39 

Added lines #L38 - L39 were not covered by tests
}

//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();

Check warning on line 45 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L43-L45 

Added lines #L43 - L45 were not covered by tests
if (varType != GDL_STRUCT) e->Throw("Expression " + e->GetString(i) + "must be a structure in this context.");
DStructGDL* s = e->GetParAs<DStructGDL>(i);
Expand All @@ -46,14 +60,15 @@ namespace lib {
return Mat;
}

Check warning on line 61 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L59-L61 

Added lines #L59 - L61 were not covered by tests
}
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();

Check warning on line 68 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L64-L68 

Added lines #L64 - L68 were not covered by tests
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];

Check warning on line 74 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L73-L74 

Added lines #L73 - L74 were not covered by tests
Expand All @@ -63,8 +78,9 @@ namespace lib {
}
}
return ret;

Check warning on line 80 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L80 

Added line #L80 was not covered by tests
}
}

//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();
Expand All @@ -91,11 +107,21 @@ namespace lib {
}
return s;

Check warning on line 108 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L108 

Added line #L108 was not covered by tests
}


DPtrGDL* convertToPtr(const SPMATRowMajDbl *Mat) {
DPtr pointer = (DPtr)Mat; //(DPtr)(MatVec.data());
return new DPtrGDL(pointer);

Check warning on line 113 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L111-L113 

Added lines #L111 - L113 were not covered by tests
}

BaseGDL* sprsin_fun(EnvT* e) {

Check warning on line 116 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L116 

Added line #L116 was not covered by tests
static bool warned=false;
SizeT nParam = e->NParam(); //1 or 4

Check warning on line 118 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L118 

Added line #L118 was not covered by tests
if (nParam != 1 && nParam != 4) e->Throw("Incorrect number of arguments.");
double thresh=d_thr;
static int COLUMN = e->KeywordIx("COLUMN");

Check warning on line 121 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L120-L121 

Added lines #L120 - L121 were not covered by tests
if (e->KeywordSet(COLUMN)) {
e->Throw("GDL's SPRSIN does not yet support the COLUMN keyword. Please report or use transposed arrays.");

Check warning on line 123 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L123 

Added line #L123 was not covered by tests
}
static int THRESHOLD=e->KeywordIx("THRESHOLD");

Check warning on line 125 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L125 

Added line #L125 was not covered by tests
if (e->KeywordSet(THRESHOLD)){
e->AssureDoubleScalarKW(THRESHOLD,thresh);

Check warning on line 127 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L127 

Added line #L127 was not covered by tests
Expand All @@ -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);

Check warning on line 135 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L134-L135 

Added lines #L134 - L135 were not covered by tests
if (nCols != nRows && !warned) {
Message("NOTE: use of SPRSIN with non-square matrices is a GDL extension.");
warned=true;

Check warning on line 138 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L137-L138 

Added lines #L137 - L138 were not covered by tests
}
DDoubleGDL* var = e->GetParAs<DDoubleGDL>(0);
SPMATRowMajDbl Mat(nRows,nCols);
SPMATRowMajDbl *Mat= new SPMATRowMajDbl(nRows,nCols);

Check warning on line 141 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L140-L141 

Added lines #L140 - L141 were not covered by tests
//import and prune wrt. threshold
for (auto i=0; i<nCols; ++i) for (auto j=0; j<nRows; ++j) if (fabs((*var)[j*nCols+i])>=thresh) Mat.insert(j,i)= (*var)[j*nCols+i];
Mat.makeCompressed();
return convertToStruct(Mat);
}
for (auto i=0; i<nCols; ++i) for (auto j=0; j<nRows; ++j) if (fabs((*var)[j*nCols+i])>=thresh) Mat->insert(j,i)= (*var)[j*nCols+i];
Mat->makeCompressed();
return convertToPtr(Mat);

Check warning on line 145 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L144-L145 

Added lines #L144 - L145 were not covered by tests
} else if (nParam == 4) {
DLongGDL* cols = e->GetParAs<DLongGDL>(0);
int nCols=cols->N_Elements();
DLongGDL* rows = e->GetParAs<DLongGDL>(1);
int nRows=rows->N_Elements();

Check warning on line 150 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L147-L150 

Added lines #L147 - L150 were not covered by tests
if (nCols != nRows) e->Throw("Vector "+e->GetString(1) + " must have "+ i2s(nCols) + " elements.");
DDoubleGDL* vals = e->GetParAs<DDoubleGDL>(2);
int nVals=vals->N_Elements();

Check warning on line 153 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L152-L153 

Added lines #L152 - L153 were not covered by tests
if (nVals != nRows) e->Throw("Vector "+e->GetString(2) + " must have "+ i2s(nCols) + " elements.");
DLongGDL* size = e->GetParAs<DLongGDL>(3);
SPMATRowMajDbl *Mat= new SPMATRowMajDbl((*size)[0],(*size)[0]); //only square matrices this way.
Mat->reserve(nCols);

Check warning on line 157 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L155-L157 

Added lines #L155 - L157 were not covered by tests
for (auto i=0; i< nCols; ++i) Mat->insert((*rows)[i],(*cols)[i])=(*vals)[i];
return convertToPtr(Mat);

Check warning on line 159 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L159 

Added line #L159 was not covered by tests
} else e->Throw("Incorrect number of arguments.");

return new DLongGDL(0);
}

Expand All @@ -124,30 +169,34 @@ namespace lib {
if (e->KeywordSet(THRESHOLD)){
e->AssureDoubleScalarKW(THRESHOLD,thresh);

Check warning on line 170 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L170 

Added line #L170 was not covered by tests
}
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));

Check warning on line 174 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L172-L174 

Added lines #L172 - L174 were not covered by tests
Mat3->prune(thresh);
return convertToPtr(Mat3);

Check warning on line 176 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L176 

Added line #L176 was not covered by tests
}

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();

Check warning on line 183 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L179-L183 

Added lines #L179 - L183 were not covered by tests
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);

Check warning on line 186 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L186 

Added line #L186 was not covered by tests
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<DDoubleGDL>(1);

Check warning on line 188 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L188 

Added line #L188 was not covered by tests
SPVecDbl Vec(m);
//import
for (auto i=0; i<m; ++i) Vec.insert(i)= (*var)[i];
SPMATRowMajDbl Mat3 = Mat1*Vec;
return convertToGDL(Mat3.transpose());
SPMATRowMajDbl* Mat3 = new SPMATRowMajDbl(((*Mat)*Vec).transpose());
return convertToGDL(Mat3);

Check warning on line 193 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L192-L193 

Added lines #L192 - L193 were not covered by tests
}

BaseGDL* fulstr_fun(EnvT* e){
SizeT nParam = e->NParam(1);

Check warning on line 197 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L196-L197 

Added lines #L196 - L197 were not covered by tests
SPMATRowMajDbl Mat=getFromStruct(e, 0);
// SPMATRowMajDbl Mat=getFromIndexInMatVec(e, 0);
SPMATRowMajDbl *Mat=getFromPtr(e, 0);
return convertToGDL(Mat);

Check warning on line 200 in src/sparse_matrix.cpp

View check run for this annotation

Codecov / codecov/patch

src/sparse_matrix.cpp#L199-L200 

Added lines #L199 - L200 were not covered by tests
}
}

0 comments on commit e39259f

Please sign in to comment.