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
|
#-------------------------------------------------------------------------------
# SPEX/Python/spex_python_demo.py: demo of 3 backslash functions with different input
# matrices
#-------------------------------------------------------------------------------
# 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
#------------------------------------------------------------------------------
# SPEX is a package for solving sparse linear systems of equations
# with a roundoff-free integer-preserving method. The result is
# always exact, unless the matrix A is perfectly singular.
# Import SPEX
import SPEXpy as SPEX
from SPEXpy import Options
# Import scientific computing
import numpy as np
from numpy.random import default_rng
from scipy.sparse import csc_matrix
from scipy.sparse import random
from scipy import stats
##--------------------------------------------------------------------------
## Cholesky
##--------------------------------------------------------------------------
# Create A and B
row = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2])
col = np.array([0, 1, 2, 0, 1, 2, 0, 1, 2])
data = np.array([4, 12, -16, 12, 37, -43, -16, -43, 98],dtype=np.float64)
A=csc_matrix((data, (row, col)), shape=(3, 3))
b=np.ones(3,dtype=np.float64)
# Solve
x=SPEX.cholesky_backslash(A,b)
print(x)
##--------------------------------------------------------------------------
## Left LU
##--------------------------------------------------------------------------
# Generate a random sparse matrix A and populate b
n=10
rng = default_rng()
rvs = stats.poisson(25, loc=10).rvs
S = random(n, n, density=0.7, random_state=rng, data_rvs=rvs)
S2=S+np.eye(n)
A=csc_matrix(S2)
b=np.ones(n,dtype=np.float64)
# Solve
options=Options("string")
x=SPEX.lu_backslash(A,b,options)
print(x)
##--------------------------------------------------------------------------
## Backslash
##--------------------------------------------------------------------------
# Use the previous matrices
# Solve
x=SPEX.backslash(A,b)
print(x)
#SPEX.backslash always returns the output as float64
|
