1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
|
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mkl_dss.h"
#include "mkl_types.h"
#include "mkl_spblas.h"
#define KIVEC(A) int A##n, const int*A##p
#define KDVEC(A) int A##n, const double*A##p
#define DVEC(A) int A##n, double*A##p
#define OK return 0;
void check_error(int error)
{
if(error != MKL_DSS_SUCCESS) {
printf ("Solver returned error code %d\n", error);
exit (1);
}
}
int dss(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
MKL_INT nRows = rowsn-1, nCols = rn, nNonZeros = valsn, nRhs = 1;
MKL_INT *rowIndex = (MKL_INT*) rowsp;
MKL_INT *columns = (MKL_INT*) colsp;
double *values = (double*) valsp;
_DOUBLE_PRECISION_t *rhs = (_DOUBLE_PRECISION_t*) xp;
// _DOUBLE_PRECISION_t *obtrhs = (_DOUBLE_PRECISION_t*) malloc((nCols)*sizeof(_DOUBLE_PRECISION_t));
_DOUBLE_PRECISION_t *solValues = (_DOUBLE_PRECISION_t*) rp;
_MKL_DSS_HANDLE_t handle;
_INTEGER_t error;
// _CHARACTER_t *uplo;
MKL_INT opt;
opt = MKL_DSS_DEFAULTS;
error = dss_create(handle, opt);
check_error(error);
opt = MKL_DSS_NON_SYMMETRIC;
error = dss_define_structure(handle, opt, rowIndex, nRows, nCols, columns, nNonZeros);
check_error(error);
opt = MKL_DSS_DEFAULTS;
error = dss_reorder(handle, opt, 0);
check_error(error);
opt = MKL_DSS_INDEFINITE;
error = dss_factor_real(handle, opt, values);
check_error(error);
int j;
for (j = 0; j < nCols; j++) {
solValues[j] = 0.0;
}
// Solve system
opt = MKL_DSS_REFINEMENT_ON;
error = dss_solve_real(handle, opt, rhs, nRhs, solValues);
check_error(error);
OK
}
|
