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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
#-------------------------------------------------------------------------------
# SPEX/Python/SPEX/spex_connect.py: link SPEX to use in Python
#-------------------------------------------------------------------------------
# SPEX: (c) 2022-2024, Christopher Lourenco, Jinhao Chen,
# Lorena Mejia Domenzain, Erick Moreno-Centeno, and Timothy A. Davis.
# All Rights Reserved.
# SPDX-License-Identifier: GPL-2.0-or-later or LGPL-3.0-or-later
#------------------------------------------------------------------------------
import ctypes
import numpy as np
from numpy.ctypeslib import ndpointer
from .SPEX_error import *
def spex_connect( A, b, order, charOut, algorithm ):
## A is a scipy.sparse.csc_matrix (data must be float64) #technically it only needs to be numerical
## b is a numpy.array (data must be float64)
##--------------------------------------------------------------------------
## Load the library with the "C bridge code"
##--------------------------------------------------------------------------
lib = ctypes.CDLL('libspexpython.so.3')
c_backslash = lib.spex_python
##--------------------------------------------------------------------------
## Specify the parameter types and return type of the C function
##--------------------------------------------------------------------------
c_backslash.argtypes = [ctypes.POINTER(ctypes.c_void_p),
ndpointer(dtype=np.int64, ndim=1, flags=None),
ndpointer(dtype=np.int64, ndim=1, flags=None),
ndpointer(dtype=np.float64, ndim=1, flags=None),
ndpointer(dtype=np.float64, ndim=1, flags=None),
ctypes.c_int,
ctypes.c_int,
ctypes.c_int,
ctypes.c_int,
ctypes.c_int,
ctypes.c_bool]
c_backslash.restype = ctypes.c_int
n=A.shape[0] #number of columns/rows of A
x_v = (ctypes.c_void_p*n)()
##--------------------------------------------------------------------------
## Solve Ax=b using REF Sparse Cholesky Factorization
##--------------------------------------------------------------------------
ok=c_backslash(x_v,
A.indptr.astype(np.int64), #without the cast it would be int32 and it would not be compatible with the C method
A.indices.astype(np.int64),
A.data.astype(np.float64),
b,
n,
n,
A.nnz,
order,
algorithm,
charOut)
if ok!=0:
raise SPEX_error(determine_error(ok))
##--------------------------------------------------------------------------
## Cast solution into correct type (string or double)
##--------------------------------------------------------------------------
if charOut:
val = ctypes.cast(x_v, ctypes.POINTER(ctypes.c_char_p))
x=[]
for i in range(n):
x.append(val[i])
else:
#x = ctypes.cast(x_v, ctypes.POINTER(ctypes.c_double))
x=[]
for i in range(n):
val=ctypes.cast(x_v[i], ctypes.POINTER(ctypes.c_double))
x.append(val[0]) ##this can also be changed to be a numpy array instead of a list
return np.array(x)
|
