solver_superlu.cpp 2.75 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
3
4
5
6
7
8
9
10
11
12
#define _CRT_SECURE_NO_WARNINGS
#include "inmost_solver.h"
#if defined(USE_SOLVER)
#if defined(USE_SOLVER_SUPERLU)
#include "superlu/slu_ddefs.h"
#include "solver_superlu.h"

struct SuperLU_data
{
	SuperMatrix A, L, U;
	int * perm_r;
	int * perm_c;
13
	int * remap;
Kirill Terekhov's avatar
Kirill Terekhov committed
14
15
16
17
18
19
20
21
22
23
24
	int size, info;
	superlu_options_t options;
	SuperLUStat_t stat;
};

void SolverInitializeSuperLU(int * argc,char *** argv, const char * file_options)
{
	//read options from file and arguments
}


25
void MatrixFillSuperLU(void * data, int size, int nnz, int * col_ranges, int * col_positions, double * col_values, int * remap)
Kirill Terekhov's avatar
Kirill Terekhov committed
26
27
28
{
	SuperLU_data * m = static_cast<SuperLU_data *>(data);
	dCreate_CompCol_Matrix(&m->A,size,size,nnz,col_values,col_positions,col_ranges,SLU_NR,SLU_D,SLU_GE);
29
	//dPrint_CompCol_Matrix("A", &m->A);
Kirill Terekhov's avatar
Kirill Terekhov committed
30
	m->size = size;
31
	m->remap = remap;
Kirill Terekhov's avatar
Kirill Terekhov committed
32
33
34
35
	m->perm_c = new int [size];
	m->perm_r = new int [size];
}

36
37
38
39
40
41
42
43
44
int * MatrixRemapArraySuperLU(void * data)
{
	return static_cast<SuperLU_data *>(data)->remap;
}

int MatrixSizeSuperLU(void * data)
{
	return static_cast<SuperLU_data *>(data)->size;
}
Kirill Terekhov's avatar
Kirill Terekhov committed
45
46
47
48
49
50
51
52
53
54

void SolverDestroyDataSuperLU(void ** data)
{
	SuperLU_data * m = static_cast<SuperLU_data *>(*data);
	Destroy_CompCol_Matrix(&m->A);
	Destroy_SuperNode_Matrix(&m->L);
	Destroy_CompCol_Matrix(&m->U);
	StatFree(&m->stat);
	if( m->perm_c != NULL ) delete [] m->perm_c;
	if( m->perm_r != NULL ) delete [] m->perm_r;
55
	if( m->remap != NULL ) delete [] m->remap; //allocated outside
Kirill Terekhov's avatar
Kirill Terekhov committed
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
	delete m;
	*data = NULL;
}

void SolverInitDataSuperLU(void ** data)
{
	if( data == NULL ) throw INMOST::DataCorruptedInSolver;
	if( *data != NULL ) 
		SolverDestroyDataSuperLU(data);
	*data = static_cast<void *>(new SuperLU_data);
	SuperLU_data * m = static_cast<SuperLU_data *>(*data);
	set_default_options(&(m->options));
	StatInit(&(m->stat));
	m->perm_c = NULL;
	m->perm_r = NULL;
	m->size = 0;
}


75
bool SolverSolveSuperLU(void * data, void * rhs_in_sol_out_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
76
77
78
79
{
	SuperLU_data * m = static_cast<SuperLU_data *>(data);
	assert(m->size != 0);
	SuperMatrix B;
80
	dCreate_Dense_Matrix(&B,m->size,1,(double *)rhs_in_sol_out_data,m->size,SLU_DN,SLU_D,SLU_GE);
Kirill Terekhov's avatar
Kirill Terekhov committed
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
	dgssv(&m->options,&m->A,m->perm_c,m->perm_r,&m->L,&m->U,&B,&m->stat,&m->info);
	Destroy_SuperMatrix_Store(&B);
	return m->info == 0;
}

const char * SolverConvergedReasonSuperLU(void * data)
{
	SuperLU_data * m = static_cast<SuperLU_data *>(data);
	static char reason_str[4096];
	if( m->info <= m->size )
		sprintf(reason_str,"diagonal entry of U-factor is exactly singular at %d/%d",m->info,m->size);
	else if( m->info > m->size )
		sprintf(reason_str,"memory allocation failed after %d bytes were allocated",m->info-m->size);
	else if( m->info == 0 )
		strcpy(reason_str,"factorization was successfull");
	else
		sprintf(reason_str,"unknown exit code %d",m->info);
	return reason_str;
}

#endif
#endif