simple API to SuperLU

SuperLU_API.cpp

#include "slu_ddefs.h"
#include <vector>

struct Interal_SuperLU
{
	superlu_options_t options;

	SuperLUStat_t stat;

	std::vector<int> etree;
	std::vector<int> perm_c;
	std::vector<int> perm_r;

	SuperMatrix l;
	SuperMatrix u;
};

void dss_sym_destroy(Interal_SuperLU *sup);

Interal_SuperLU *dss_sym_create();

void dss_sym_reorder(Interal_SuperLU *sup, int dim, int nnz, int *csc_columns, int *csc_rows, double *csc_values);

/*
relax(input) int
To control degree of relaxing supernodes.If the number
of nodes(columns) in a subtree of the elimination tree is less
than relax, this subtree is considered as one supernode,
regardless of the row structures of those columns.

panel_size(input) int
A panel consists of at most panel_size consecutive columns.
*/
int dss_sym_factor(Interal_SuperLU *sup, int dim, int nnz, int *csc_columns, int *csc_rows, double *csc_values,
	int relax, int panel_size);

int dss_sym_refactor(Interal_SuperLU *sup, int dim, int nnz, int *csc_columns, int *csc_rows, double *csc_values,
	int relax, int panel_size);

int dss_sym_solve(Interal_SuperLU *sup, int dim, double *b, double *x);

void dss_sym_destroy(Interal_SuperLU *sup)
{
	if (sup->l.Store)
		Destroy_SuperNode_Matrix(&sup->l);
	if (sup->u.Store)
		Destroy_CompCol_Matrix(&sup->u);

	StatFree(&sup->stat);

	delete sup;
}

Interal_SuperLU *dss_sym_create()
{
	Interal_SuperLU *sup = new Interal_SuperLU();

	set_default_options(&sup->options);

	sup->options.Fact = DOFACT;

	sup->options.ColPerm = COLAMD;

	sup->options.Equil = NO;
	sup->options.DiagPivotThresh = 1;
	sup->options.ReplaceTinyPivot = YES;
	sup->options.PivotGrowth = NO;
	sup->options.ConditionNumber = NO;
	sup->options.SymmetricMode = YES;
	sup->options.PrintStat = NO;
	sup->options.SymPattern = YES;

	sup->l.Store = 0;
	sup->u.Store = 0;

	StatInit(&sup->stat);

	return sup;
}

void dss_sym_reorder(Interal_SuperLU *sup, int dim, int nnz, int *csc_columns, int *csc_rows, double *csc_values)
{
	SuperMatrix a;

	dCreate_CompCol_Matrix(&a, dim, dim, nnz, csc_values, csc_rows, csc_columns, SLU_NC, SLU_D, SLU_GE);

	sup->etree.resize(dim);
	sup->perm_c.resize(dim);
	sup->perm_r.resize(dim);

	get_perm_c(sup->options.ColPerm, &a, sup->perm_c.data());

	Destroy_SuperMatrix_Store(&a);
}

int dss_sym_factor(Interal_SuperLU *sup, int dim, int nnz, int *csc_columns, int *csc_rows, double *csc_values,
	int relax, int panel_size)
{
	int ret_code;
	SuperMatrix a;
	SuperMatrix ap; // a permuted

	// destroy previously created factor
	if (sup->l.Store)
 	{
		Destroy_SuperNode_Matrix(&sup->l);
	 	sup->l.Store = 0;
	}
	if (sup->u.Store)
	{
		Destroy_CompCol_Matrix(&sup->u);
	 	sup->u.Store = 0;
	}
	sup->options.Fact = DOFACT;

	dCreate_CompCol_Matrix(&a, dim, dim, nnz, csc_values, csc_rows, csc_columns, SLU_NC, SLU_D, SLU_GE);

	sp_preorder(&sup->options, &a, sup->perm_c.data(), sup->etree.data(), &ap);

	dgstrf(&sup->options, &ap, relax, panel_size, sup->etree.data(), NULL, 0,
		sup->perm_c.data(), sup->perm_r.data(), &sup->l, &sup->u, &sup->stat, &ret_code);

	Destroy_SuperMatrix_Store(&a);
	Destroy_CompCol_Permuted(&ap);

	return ret_code;
}

int dss_sym_refactor(Interal_SuperLU *sup, int dim, int nnz, int *csc_columns, int *csc_rows, double *csc_values,
	int relax, int panel_size)
{
	int ret_code;
	SuperMatrix a;
	SuperMatrix ap; // `a' permuted

	sup->options.Fact = SamePattern_SameRowPerm;

	dCreate_CompCol_Matrix(&a, dim, dim, nnz, csc_values, csc_rows, csc_columns, SLU_NC, SLU_D, SLU_GE);

	sp_preorder(&sup->options, &a, sup->perm_c.data(), sup->etree.data(), &ap);

	dgstrf(&sup->options, &ap, relax, panel_size, sup->etree.data(), NULL, 0,
		sup->perm_c.data(), sup->perm_r.data(), &sup->l, &sup->u, &sup->stat, &ret_code);

	Destroy_SuperMatrix_Store(&a);
	Destroy_CompCol_Permuted(&ap);

	return ret_code;
}

int dss_sym_solve(Interal_SuperLU *sup, int dim, double *b, double *x)
{
	int ret_code;
	SuperMatrix rhs;

	copy(b, b + dim, x);

	dCreate_Dense_Matrix(&rhs, dim, 1, x, dim, SLU_DN, SLU_D, SLU_GE);

	dgstrs(NOTRANS, &sup->l, &sup->u, sup->perm_c.data(), sup->perm_r.data(), &rhs, &sup->stat, &ret_code);

	Destroy_SuperMatrix_Store(&rhs);

	Destroy_SuperNode_Matrix(&sup->l);
	Destroy_CompCol_Matrix(&sup->u);

	return ret_code;
}