Skip to content

Commit

Permalink
SUPERLU+ARPACK integration
Browse files Browse the repository at this point in the history
  • Loading branch information
@lisitsyn
lisitsyn committed Sep 19, 2011
1 parent 8b3602c commit 0af9076
Showing 1 changed file with 130 additions and 3 deletions.
133 changes: 130 additions & 3 deletions src/shogun/mathematics/arpack.cpp
View file
Expand Up @@ -18,14 +18,18 @@
#include <shogun/io/SGIO.h>
#include <string.h>

#ifdef HAVE_SUPERLU
#include <superlu/slu_ddefs.h>
#endif

using namespace shogun;

namespace shogun
{

void arpack_dsaeupd_wrap(double* matrix, double* rhs_diag, int n, int nev, const char* which,
int mode, bool pos, double shift, double tolerance, double* eigenvalues,
double* eigenvectors, int& status)
int mode, bool pos, double shift, double tolerance,
double* eigenvalues, double* eigenvectors, int& status)
{
// loop vars
int i,j;
Expand Down Expand Up @@ -92,6 +96,27 @@ void arpack_dsaeupd_wrap(double* matrix, double* rhs_diag, int n, int nev, const

// ipiv for LUP factorization
int* ipiv = NULL;

#ifdef HAVE_SUPERLU
char equed[1];
void* work;
int lwork = 0;
SuperMatrix slu_A, slu_L, slu_U, slu_B, slu_X;
superlu_options_t options;
SuperLUStat_t stat;
mem_usage_t mem_usage;
int* perm_c = SG_MALLOC(int, n);
int* perm_r = SG_MALLOC(int, n);
int* etree = SG_MALLOC(int, n);
double* R = SG_MALLOC(double, n);
double* C = SG_MALLOC(double, n);
double ferr;
double berr;
double rcond;
double rpg;
int slu_info;
#endif

// shift-invert mode init
if (mode==3)
{
Expand All @@ -110,6 +135,57 @@ void arpack_dsaeupd_wrap(double* matrix, double* rhs_diag, int n, int nev, const
matrix[i*n+i] -= shift;
}

#ifdef HAVE_SUPERLU
SG_SDEBUG("SUPERLU: Constructing sparse matrix.\n");
int nnz = 0;
for (i=0; i<n*n; i++)
{
if (matrix[i]!=0.0)
nnz++;
}
double* val = doubleMalloc(nnz);
int* rowind = intMalloc(nnz);
int* colptr = intMalloc(n+1);
nnz = 0;
for (i=0; i<n; i++)
{
colptr[i] = nnz;
for (j=0; j<n; j++)
{
if (matrix[i*n+j]!=0.0)
{
val[nnz] = matrix[i*n+j];
rowind[nnz] = j;
nnz++;
}
}
}
colptr[i] = nnz;
dCreate_CompCol_Matrix(&slu_A,n,n,nnz,val,rowind,colptr,SLU_NC,SLU_D,SLU_GE);

set_default_options(&options);
options.Equil = YES;
options.ColPerm = MMD_AT_PLUS_A;
options.SymmetricMode = YES;
StatInit(&stat);

slu_info = 0;
double* tmp_rhs = SG_MALLOC(double, n);
for (i=0; i<n; i++) tmp_rhs[i] = 20.0;
dCreate_Dense_Matrix(&slu_B,n,1,tmp_rhs,n,SLU_DN,SLU_D,SLU_GE);
dCreate_Dense_Matrix(&slu_X,n,1,tmp_rhs,n,SLU_DN,SLU_D,SLU_GE);
slu_B.ncol = 0;
SG_SDEBUG("SUPERLU: Factorizing\n");
dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
&slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond, &ferr, &berr,
&mem_usage, &stat, &slu_info);
if (slu_info)
{
SG_SERROR("SUPERLU: Failed to factorize matrix, got %d code\n", slu_info);
}
SG_FREE(tmp_rhs);
options.Fact = FACTORED;
#else
// compute factorization according to pos value
if (pos)
{
Expand All @@ -124,6 +200,7 @@ void arpack_dsaeupd_wrap(double* matrix, double* rhs_diag, int n, int nev, const
ipiv = SG_CALLOC(int, n);
clapack_dgetrf(CblasColMajor,n,n,matrix,n,ipiv);
}
#endif
}
// main computation loop
SG_SDEBUG("ARPACK: Starting main computation DSAUPD loop.\n");
Expand All @@ -147,35 +224,68 @@ void arpack_dsaeupd_wrap(double* matrix, double* rhs_diag, int n, int nev, const
{
if (!rhs_diag)
{
#ifdef HAVE_SUPERLU
dCreate_Dense_Matrix(&slu_B,n,1,workd+ipntr[0]-1,n,SLU_DN,SLU_D,SLU_GE);
dCreate_Dense_Matrix(&slu_X,n,1,workd+ipntr[1]-1,n,SLU_DN,SLU_D,SLU_GE);
slu_info = 0;
dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
&slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond,
&ferr, &berr, &mem_usage, &stat, &slu_info);
if (slu_info)
SG_SERROR("SUPERLU: GOT %d\n", slu_info);
#else
for (i=0; i<n; i++)
(workd+ipntr[1]-1)[i] = (workd+ipntr[0]-1)[i];

if (pos)
if (pos)
clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
else
clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
#endif
}
else
{
if (ido==-1)
{
#ifdef HAVE_SUPERLU
dCreate_Dense_Matrix(&slu_B,n,1,workd+ipntr[0]-1,n,SLU_DN,SLU_D,SLU_GE);
dCreate_Dense_Matrix(&slu_X,n,1,workd+ipntr[1]-1,n,SLU_DN,SLU_D,SLU_GE);
slu_info = 0;
dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
&slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond,
&ferr, &berr, &mem_usage, &stat, &slu_info);
if (slu_info)
SG_SERROR("SUPERLU: GOT %d\n", slu_info);
#else
for (i=0; i<n; i++)
(workd+ipntr[1]-1)[i] = rhs_diag[i]*(workd+ipntr[0]-1)[i];

if (pos)
clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
else
clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
#endif
}
if (ido==1)
{
#ifdef HAVE_SUPERLU
dCreate_Dense_Matrix(&slu_B,n,1,workd+ipntr[2]-1,n,SLU_DN,SLU_D,SLU_GE);
dCreate_Dense_Matrix(&slu_X,n,1,workd+ipntr[1]-1,n,SLU_DN,SLU_D,SLU_GE);
slu_info = 0;
dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
&slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond,
&ferr, &berr, &mem_usage, &stat, &slu_info);
if (slu_info)
SG_SERROR("SUPERLU: GOT %d\n", slu_info);
#else
for (i=0; i<n; i++)
(workd+ipntr[1]-1)[i] = (workd+ipntr[2]-1)[i];

if (pos)
clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
else
clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
#endif
}
if (ido==2)
{
Expand All @@ -189,6 +299,23 @@ void arpack_dsaeupd_wrap(double* matrix, double* rhs_diag, int n, int nev, const

if (!pos && mode==3) SG_FREE(ipiv);

#ifdef HAVE_SUPERLU
SUPERLU_FREE(perm_r);
SUPERLU_FREE(perm_c);
SUPERLU_FREE(R);
SUPERLU_FREE(C);
SUPERLU_FREE(etree);
if (mode==3)
{
Destroy_CompCol_Matrix(&slu_A);
StatFree(&stat);
Destroy_SuperMatrix_Store(&slu_B);
Destroy_SuperMatrix_Store(&slu_X);
Destroy_SuperNode_Matrix(&slu_L);
Destroy_CompCol_Matrix(&slu_U);
}
#endif

// check if DSAUPD failed
if (info<0)
{
Expand Down

0 comments on commit 0af9076

Please sign in to comment.