Coverage for pyxfem/core/sparse_tools.py: 55%
109 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4Tools for sparse matrix, not present or ineficient for my usage in scipy
6Created on Wed Feb 28 09:49:32 2024
7@author: nchevaug
8"""
9from typing import cast
10import numpy as np
11import scipy as sp
12import scipy.sparse.linalg as spla
13from scipy.linalg import cho_solve, cho_factor, lu_factor, lu_solve
14from scipy.linalg import lapack
15from pyxfem.core.config import HAVE_UMFPACK
16if HAVE_UMFPACK:
17 from scikits import umfpack
18#pylint: disable=invalid-name
20#USE_UMFPACK = True
21spla.use_solver(useUmfpack = HAVE_UMFPACK)
22def factorized(A, assume_a="gen", overwrite_a=False):
23 """Returns a function solve(b) using pre-factored dense LAPACK routines.
24 or scipy.sparse.factorized
25 assume_a : 'gen' (LU), 'pos' (Cholesky), or 'sym' (LDL^T)
26 """
27 if sp.sparse.issparse(A):
28 if HAVE_UMFPACK:
29 A = sp.sparse.csc_matrix(A, dtype = float)
30 # fall back since umfpack.splu is not fully compatible with csc_array
31 lu = umfpack.splu(A)
32 def solve(b):
33 return lu.solve(b)
34 return solve
35 return spla.factorized(A)
37 if assume_a in ["pos", "positive"]:
38 # Cholesky LL^T
39 c, low = cho_factor(A, lower=True, overwrite_a=overwrite_a)
40 return lambda b: cho_solve((c, low), b)
41 elif assume_a in ["sym", "symmetric"]:
42 # Pre-factors symmetric matrix A using LAPACK dsytrf
43 # and returns a fast solve(b) function.
45 # dsytrf requires Fortran-contiguous array, and float64 (or float32 for ssytrf)
46 A_f = np.asfortranarray(A, dtype=np.float64)
47 #uplo = 1 if lower else 0
48 uplo =1
49 # 1. Factorize once (O(N^3))
50 (sytrf,) = lapack.get_lapack_funcs(('sytrf',), (A,))
51 ldl_packed, ipiv, info = sytrf(A_f, lower=uplo, overwrite_a = overwrite_a)
52 if info != 0:
53 raise np.linalg.LinAlgError(f"dsytrf failed with info={info} (singular matrix)")
54 # 2. Return direct solve closure (O(N^2))
55 def ldl_solve(b):
56 b_f = np.asfortranarray(b, dtype=np.float64)
57 (sytrs,) = lapack.get_lapack_funcs(('sytrs',), (A,))
58 x, solve_info = sytrs(ldl_packed, ipiv, b_f,
59 lower=uplo, overwrite_a = overwrite_a)
60 if solve_info != 0:
61 raise np.linalg.LinAlgError(f"dsytrs failed with info={solve_info}")
62 return x
63 return ldl_solve
64 else:
65 # General LU
66 lu, piv = lu_factor(A)
67 return lambda b: lu_solve((lu, piv), b)
70def regular_block_diag_to_csr(block_diag):
71 ''' from an np.array of shape (nblock, sblock, sblock) representing a block diagonal
72 matrix A of shape (nblock*sblock) with block size sblock
73 constuct the csr_array representing A.
74 '''
75 fun_name = regular_block_diag_to_csr.__name__
76 if len(block_diag.shape) != 3:
77 raise ValueError('Error in '+fun_name+' input block_diag must be 3d.')
78 nblock = block_diag.shape[0] #number of diagonal block
79 xblock = block_diag.shape[1] #size of the block
80 yblock = block_diag.shape[2] #size of the block
81 lin, col, shape = block_diag_compute_index(nblock, xblock, yblock)
82 return sp.sparse.csr_array((block_diag.flatten(), (lin, col)), shape)
84def regular_block_diag_to_csc(block_diag):
85 ''' from an np.array of shape (nblock, sblock, sblock) representing a block diagonal
86 matrix A of shape (nblock*sblock) with block size sblock
87 constuct the csr_array representing A.
88 '''
89 fun_name = regular_block_diag_to_csc.__name__
90 if len(block_diag.shape) != 3:
91 raise ValueError('Error in '+fun_name+' input block_diag must be 3d.')
92 nblock = block_diag.shape[0] #number of diagonal block
93 xblock = block_diag.shape[1] #size of the block
94 yblock = block_diag.shape[2] #size of the block
95 lin, col, shape = block_diag_compute_index(nblock, xblock, yblock)
96 return sp.sparse.csc_array((block_diag.flatten(), (lin, col)), shape)
98def block_diag_compute_index(nblock, xblock, yblock):
99 """ compute index of block diag matrix in order to pass to csc_array """
100 col = np.repeat(np.arange(yblock*nblock).reshape((-1, yblock)), xblock, axis=0).flatten()
101 lin = np.arange(xblock*nblock).repeat(yblock)
102 shape = (xblock*nblock, yblock*nblock)
103 return lin, col, shape
105def vblock(x0, x1):
106 return [x0, x1]
108def split_vblock(x):
109 return x[0], x[1]
111def split_mblock(A):
112 A00 = A[0][0]
113 A01 = A[0][1]
114 A10 = A[1][0]
115 A11 = A[1][1]
116 return A00, A01, A10, A11
118def mblock(A00,A01, A10, A11):
119 return [[A00, A01], [A10, A11]]
121def block_norm(x):
122 return np.linalg.norm(np.concatenate([x[0], x[1]]))
124def block_isub(x,y):
125 x[0] -= y[0]
126 x[1] -= y[1]
128def block_linsolve(A, b, solveA00 = None):
129 A00, A01, A10, A11 = split_mblock(A)
130 b0, b1 = split_vblock(b)
131 if solveA00 is None:
132 solveA00 = factorized(A00)
133 x0, x1 = linsolve_schurr(solveA00, A01, A10, A11, b0, b1)
134 return vblock(x0, x1)
136def linsolve_schurr(solveA11, A12, A21, A22, b1, b2):
137 ''' Solve block linear system using shcurr complement '''
138 A12dense = A12 if isinstance(A12, np.ndarray) else A12.todense()
139 K = A22 - A21@solveA11(A12dense)
140 y = b2 - A21@solveA11(b1)
141 x2 = sp.linalg.solve(K, y, lower=False, overwrite_a=True,
142 overwrite_b=True, check_finite=False,
143 assume_a=None, transposed=False)
144 x1 = solveA11(b1 - A12@x2)
145 return x1, x2
147def rcond_estimate(A, luA = None):
148 ''' return an estimate of the reciproque condition number
149 need A and an lu factorization of A.
150 Warning, if luA is not given, it will be computed (and destroyed ..)
151 inside the function
152 '''
153 if luA is None :
154 luA = spla.splu(A)
155 norm_A = spla.norm(A, 1)
157 class LUInverseOperator(spla.LinearOperator):
158 """Opérateur linéaire représentant l'inverse d'une factorisation LU."""
159 def __init__(self, luA, shape, dtype):
160 super().__init__(dtype=dtype, shape=shape)
161 self.luA = luA
163 def _matvec(self, x):
164 return self.luA.solve(x)
166 def _rmatvec(self, x):
167 return self.luA.solve(x, trans='T')
170 A_inv_op = LUInverseOperator(luA, shape=A.shape, dtype=A.dtype)
171 norm_A_inv = cast(float,spla.onenormest(A_inv_op))
172 cond = norm_A * norm_A_inv
173 return 1./cond