Skip to content

Commit

Permalink
Accommodate the following changes in R package Matrix
Browse files Browse the repository at this point in the history
  • Loading branch information
eromero-vlc committed Sep 28, 2023
1 parent 4a24d5e commit 541d59e
Showing 1 changed file with 23 additions and 12 deletions.
35 changes: 23 additions & 12 deletions R/src/primmeR.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
#include "c_interface/primme.h"
#include "PRIMME_types.h"
#include <Rconfig.h>
#include <Rinternals.h>
#include <R_ext/BLAS.h> // for BLAS and F77_NAME

#include "Matrix.h"
Expand Down Expand Up @@ -204,6 +205,16 @@ void copyMatrix<double, ComplexMatrix>(ComplexMatrix mat, double *x, PRIMME_INT
stop("Unsupported to return complex values when using dprimme/dprimme_svds");
}

int is_ge(SEXP x) {
static const char *valid[] = {"ngeMatrix", "lgeMatrix", "dgeMatrix", ""};
return R_check_class_etc(x, valid) >= 0;
}
int is_Csparse(SEXP x) {
static const char *valid[] = {"ngCMatrix", "lgCMatrix", "dgCMatrix",
"nsCMatrix", "lsCMatrix", "dsCMatrix", "ntCMatrix", "ltCMatrix",
"dtCMatrix", ""};
return R_check_class_etc(x, valid) >= 0;
}

template <typename T>
void copyMatrix_SEXP(SEXP mat, T *x, PRIMME_INT m, int n, PRIMME_INT ld,
Expand All @@ -215,7 +226,7 @@ void copyMatrix_SEXP(SEXP mat, T *x, PRIMME_INT m, int n, PRIMME_INT ld,
} else if (is<ComplexMatrix>(mat)) {
copyMatrix(as<ComplexMatrix>(mat), x, m, n, ld, checkDimensions);
return;
} else if (!Matrix_isclass_ge_dense(mat)) {
} else if (!is_ge(mat)) {
stop("Vector/matrix type not supported");
}

Expand Down Expand Up @@ -559,7 +570,7 @@ void matrixMatvecEigs_CHM_SP(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy,
chy.z = NULL;
chy.xtype = (sizeof(T) == sizeof(double) ? CHOLMOD_REAL : CHOLMOD_COMPLEX);
chy.dtype = CHOLMOD_DOUBLE;
const double ONEf[] = {1.0, 0.0}, ZEROf[] = {0.0, 0.0};
double ONEf[] = {1.0, 0.0}, ZEROf[] = {0.0, 0.0};
CHM_CM chol_c = (CHM_CM)((void**)F::get(primme))[1];

M_cholmod_sdmult(chm, 0, ONEf, ZEROf, (const_CHM_DN)&chx, &chy, chol_c);
Expand Down Expand Up @@ -652,10 +663,10 @@ static List xprimme(Matrix<S> ortho, Matrix<S> init, SEXP A, SEXP B,
} else if (is<ComplexMatrix>(A)) {
primme->matrix = COMPLEX(A);
primme->matrixMatvec = matrixMatvecEigs_Matrix<TS, getMatrixField>;
} else if (Matrix_isclass_ge_dense(A)) {
} else if (is_ge(A)) {
primme->matrix = AS_CHM_DN(A);
primme->matrixMatvec = matrixMatvecEigs_CHM_DN<TS, getMatrixField>;
} else if (Matrix_isclass_Csparse(A)) {
} else if (is_Csparse(A)) {
Aaux[0] = AS_CHM_SP(A);
Aaux[1] = &Achol_c;
M_R_cholmod_start(&Achol_c);
Expand Down Expand Up @@ -683,10 +694,10 @@ static List xprimme(Matrix<S> ortho, Matrix<S> init, SEXP A, SEXP B,
} else if (is<ComplexMatrix>(B)) {
primme->massMatrix = COMPLEX(B);
primme->massMatrixMatvec = matrixMatvecEigs_Matrix<TS, getMassMatrixField>;
} else if (Matrix_isclass_ge_dense(B)) {
} else if (is_ge(B)) {
primme->massMatrix = AS_CHM_DN(B);
primme->massMatrixMatvec = matrixMatvecEigs_CHM_DN<TS, getMassMatrixField>;
} else if (Matrix_isclass_Csparse(B)) {
} else if (is_Csparse(B)) {
Baux[0] = AS_CHM_SP(B);
Baux[1] = &Bchol_c;
M_R_cholmod_start(&Bchol_c);
Expand Down Expand Up @@ -723,13 +734,13 @@ static List xprimme(Matrix<S> ortho, Matrix<S> init, SEXP A, SEXP B,
if (Ac) delete Ac;
if (An) delete An;
if (Af) delete Af;
if (Matrix_isclass_Csparse(A)) {
if (is_Csparse(A)) {
M_cholmod_finish(&Achol_c);
}
if (Bc) delete Bc;
if (Bn) delete Bn;
if (Bf) delete Bf;
if (Matrix_isclass_Csparse(B)) {
if (is_Csparse(B)) {
M_cholmod_finish(&Bchol_c);
}
if (fprec) delete fprec;
Expand Down Expand Up @@ -1123,7 +1134,7 @@ static void matrixMatvecSvds_CHM_SP(void *x, PRIMME_INT *ldx, void *y, PRIMME_IN
chy.z = NULL;
chy.xtype = (sizeof(T) == sizeof(double) ? CHOLMOD_REAL : CHOLMOD_COMPLEX);
chy.dtype = CHOLMOD_DOUBLE;
const double ONEf[] = {1.0, 0.0}, ZEROf[] = {0.0, 0.0};
double ONEf[] = {1.0, 0.0}, ZEROf[] = {0.0, 0.0};
CHM_CM chol_c = (CHM_CM)((void**)primme_svds->matrix)[1];

M_cholmod_sdmult(chm, *transpose?1:0, ONEf, ZEROf, (const_CHM_DN)&chx, &chy,
Expand Down Expand Up @@ -1194,10 +1205,10 @@ static List xprimme_svds(Matrix<S> orthol, Matrix<S> orthor, Matrix<S> initl,
if (is<Matrix<S> >(A)) {
primme_svds->matrix = Am = new Matrix<S>(A);
primme_svds->matrixMatvec = matrixMatvecSvds_Matrix<T, S, TS>;
} else if (Matrix_isclass_ge_dense(A)) {
} else if (is_ge(A)) {
primme_svds->matrix = AS_CHM_DN(A);
primme_svds->matrixMatvec = matrixMatvecSvds_CHM_DN<TS>;
} else if (Matrix_isclass_Csparse(A)) {
} else if (is_Csparse(A)) {
aux[0] = AS_CHM_SP(A);
aux[1] = &chol_c;
M_R_cholmod_start(&chol_c);
Expand Down Expand Up @@ -1225,7 +1236,7 @@ static List xprimme_svds(Matrix<S> orthol, Matrix<S> orthor, Matrix<S> initl,
// Destroy auxiliary memory
if (Am) delete Am;
if (Af) delete Af;
if (Matrix_isclass_Csparse(A)) {
if (is_Csparse(A)) {
M_cholmod_finish(&chol_c);
}
if (fprec) delete fprec;
Expand Down

0 comments on commit 541d59e

Please sign in to comment.