1# CHOLMOD wrapper for scikits.sparse 2 3# Copyright (C) 2008-2017 The scikit-sparse developers: 4# 5# 2008 David Cournapeau <cournape@gmail.com> 6# 2009-2015 Nathaniel Smith <njs@pobox.com> 7# 2010 Dag Sverre Seljebotn <dagss@student.matnat.uio.no> 8# 2014 Leon Barrett <lbarrett@climate.com> 9# 2015 Yuri <yuri@tsoft.com> 10# 2016-2017 Antony Lee <anntzer.lee@gmail.com> 11# 2016 Alex Grigorievskiy <alex.grigorievskiy@gmail.com> 12# 2016-2018 Joscha Reimer <jor@informatik.uni-kiel.de> 13# All rights reserved. 14# 15# Redistribution and use in source and binary forms, with or without 16# modification, are permitted provided that the following conditions are 17# met: 18# 19# - Redistributions of source code must retain the above copyright 20# notice, this list of conditions and the following disclaimer. 21# - Redistributions in binary form must reproduce the above 22# copyright notice, this list of conditions and the following 23# disclaimer in the documentation and/or other materials 24# provided with the distribution. 25# 26# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 27# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 28# INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 29# MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 30# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS 31# BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 32# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 33# TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 34# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 35# ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR 36# TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF 37# THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 38# SUCH DAMAGE. 39 40#cython: binding = True 41#cython: language_level = 3 42 43cimport numpy as np 44 45import warnings 46import numpy as np 47from scipy import sparse 48 49np.import_array() 50 51cdef extern from "numpy/arrayobject.h": 52 void PyArray_ENABLEFLAGS(np.ndarray arr, int flags) 53 54cdef extern from "cholmod_backward_compatible.h": 55 cdef enum: 56 CHOLMOD_INT, CHOLMOD_INTLONG, CHOLMOD_LONG 57 CHOLMOD_PATTERN, CHOLMOD_REAL, CHOLMOD_COMPLEX 58 CHOLMOD_DOUBLE 59 CHOLMOD_AUTO, CHOLMOD_SIMPLICIAL, CHOLMOD_SUPERNODAL 60 CHOLMOD_OK, CHOLMOD_NOT_POSDEF 61 CHOLMOD_NOT_INSTALLED, CHOLMOD_OUT_OF_MEMORY, CHOLMOD_TOO_LARGE, CHOLMOD_INVALID, CHOLMOD_GPU_PROBLEM 62 CHOLMOD_A, CHOLMOD_LDLt, CHOLMOD_LD, CHOLMOD_DLt, CHOLMOD_L 63 CHOLMOD_Lt, CHOLMOD_D, CHOLMOD_P, CHOLMOD_Pt 64 CHOLMOD_NATURAL, CHOLMOD_GIVEN, CHOLMOD_AMD, CHOLMOD_METIS, CHOLMOD_NESDIS, CHOLMOD_COLAMD, CHOLMOD_POSTORDERED 65 66 ctypedef int SuiteSparse_long 67 68 ctypedef struct cholmod_method_struct: 69 int ordering 70 71 ctypedef struct cholmod_common: 72 int supernodal 73 int status 74 int itype, dtype 75 int print 76 int nmethods, postorder 77 cholmod_method_struct * method 78 void (*error_handler)(int status, const char * file, int line, const char * msg) 79 80 ctypedef struct cholmod_sparse: 81 size_t nrow, ncol, nzmax 82 void * p # column pointers 83 void * i # row indices 84 void * x 85 int stype # 0 = regular, >0 = upper triangular, <0 = lower triangular 86 int itype, dtype, xtype 87 int sorted 88 int packed 89 90 ctypedef struct cholmod_dense: 91 size_t nrow, ncol, nzmax, d 92 void * x 93 int dtype, xtype 94 95 ctypedef struct cholmod_factor: 96 size_t n 97 size_t minor 98 void * Perm 99 int itype, xtype 100 int is_ll, is_super, is_monotonic 101 size_t xsize, nzmax, nsuper 102 void * x 103 void * p 104 void * super_ "super" 105 void * pi 106 void * px 107 108 int cholmod_start(cholmod_common *) except * 109 int cholmod_l_start(cholmod_common *) except * 110 111 int cholmod_finish(cholmod_common *) except * 112 int cholmod_l_finish(cholmod_common *) except * 113 114 int cholmod_check_common(cholmod_common *) except * 115 int cholmod_l_check_common(cholmod_common *) except * 116 117 int cholmod_print_common(const char *, cholmod_common *) except * 118 int cholmod_l_print_common(const char *, cholmod_common *) except * 119 120 int cholmod_free_sparse(cholmod_sparse **, cholmod_common *) except * 121 int cholmod_l_free_sparse(cholmod_sparse **, cholmod_common *) except * 122 123 int cholmod_check_sparse(cholmod_sparse *, cholmod_common *) except * 124 int cholmod_l_check_sparse(cholmod_sparse *, cholmod_common *) except * 125 126 int cholmod_print_sparse(cholmod_sparse *, const char *, cholmod_common *) except * 127 int cholmod_l_print_sparse(cholmod_sparse *, const char *, cholmod_common *) except * 128 129 int cholmod_free_dense(cholmod_dense **, cholmod_common *) except * 130 int cholmod_l_free_dense(cholmod_dense **, cholmod_common *) except * 131 132 int cholmod_check_dense(cholmod_dense *, cholmod_common *) except * 133 int cholmod_l_check_dense(cholmod_dense *, cholmod_common *) except * 134 135 int cholmod_print_dense(cholmod_dense *, const char *, cholmod_common *) except * 136 int cholmod_l_print_dense(cholmod_dense *, const char *, cholmod_common *) except * 137 138 int cholmod_free_factor(cholmod_factor **, cholmod_common *) except * 139 int cholmod_l_free_factor(cholmod_factor **, cholmod_common *) except * 140 141 int cholmod_check_factor(cholmod_factor *, cholmod_common *) except * 142 int cholmod_l_check_factor(cholmod_factor *, cholmod_common *) except * 143 144 int cholmod_print_factor(cholmod_factor *, const char *, cholmod_common *) except * 145 int cholmod_l_print_factor(cholmod_factor *, const char *, cholmod_common *) except * 146 147 cholmod_factor * cholmod_copy_factor(cholmod_factor *, cholmod_common *) except? NULL 148 cholmod_factor * cholmod_l_copy_factor(cholmod_factor *, cholmod_common *) except? NULL 149 150 cholmod_factor * cholmod_analyze(cholmod_sparse *, cholmod_common *) except? NULL 151 cholmod_factor * cholmod_l_analyze(cholmod_sparse *, cholmod_common *) except? NULL 152 153 int cholmod_factorize_p(cholmod_sparse *, double beta[2], 154 int * fset, size_t fsize, 155 cholmod_factor *, 156 cholmod_common *) except * 157 int cholmod_l_factorize_p(cholmod_sparse *, double beta[2], 158 SuiteSparse_long * fset, size_t fsize, 159 cholmod_factor *, 160 cholmod_common *) except * 161 162 cholmod_sparse * cholmod_submatrix(cholmod_sparse *, 163 int * rset, int rsize, 164 int * cset, int csize, 165 int values, int sorted, 166 cholmod_common *) except? NULL 167 cholmod_sparse * cholmod_l_submatrix(cholmod_sparse *, 168 SuiteSparse_long * rset, SuiteSparse_long rsize, 169 SuiteSparse_long * cset, SuiteSparse_long csize, 170 int values, int sorted, 171 cholmod_common *) except? NULL 172 173 int cholmod_updown(int update, cholmod_sparse *, cholmod_factor *, 174 cholmod_common *) except * 175 int cholmod_l_updown(int update, cholmod_sparse *, cholmod_factor *, 176 cholmod_common *) except * 177 178 cholmod_dense * cholmod_solve(int, cholmod_factor *, 179 cholmod_dense *, cholmod_common *) except? NULL 180 cholmod_dense * cholmod_l_solve(int, cholmod_factor *, 181 cholmod_dense *, cholmod_common *) except? NULL 182 183 cholmod_sparse * cholmod_spsolve(int, cholmod_factor *, 184 cholmod_sparse *, cholmod_common *) except? NULL 185 cholmod_sparse * cholmod_l_spsolve(int, cholmod_factor *, 186 cholmod_sparse *, cholmod_common *) except? NULL 187 188 int cholmod_change_factor(int to_xtype, int to_ll, int to_super, 189 int to_packed, int to_monotonic, 190 cholmod_factor *, cholmod_common *) except * 191 int cholmod_l_change_factor(int to_xtype, int to_ll, int to_super, 192 int to_packed, int to_monotonic, 193 cholmod_factor *, cholmod_common *) except * 194 195 cholmod_sparse * cholmod_factor_to_sparse(cholmod_factor *, 196 cholmod_common *) except? NULL 197 cholmod_sparse * cholmod_l_factor_to_sparse(cholmod_factor *, 198 cholmod_common *) except? NULL 199 200cdef class Common 201cdef class Factor 202 203class CholmodError(Exception): 204 pass 205 206class CholmodNotPositiveDefiniteError(CholmodError): 207 def __init__(self, message, column=None, factor=None): 208 super().__init__(message) 209 self.column = column 210 self.factor = factor 211 212class CholmodNotInstalledError(CholmodError): 213 pass 214 215class CholmodOutOfMemoryError(CholmodError): 216 pass 217 218class CholmodTooLargeError(CholmodError): 219 pass 220 221class CholmodInvalidError(CholmodError): 222 pass 223 224class CholmodGpuProblemError(CholmodError): 225 pass 226 227class CholmodWarning(UserWarning): 228 pass 229 230class CholmodTypeConversionWarning( 231 CholmodWarning, sparse.SparseEfficiencyWarning): 232 pass 233 234cdef int _integer_typenum = np.NPY_INT32 235cdef object _integer_py_dtype = np.dtype(np.int32) 236assert sizeof(int) == _integer_py_dtype.itemsize == 4 237 238cdef int _long_typenum = np.NPY_INT64 239cdef object _long_py_dtype = np.dtype(np.int64) 240assert sizeof(SuiteSparse_long) == _long_py_dtype.itemsize == 8 241 242cdef int _real_typenum = np.NPY_FLOAT64 243cdef object _real_py_dtype = np.dtype(np.float64) 244assert sizeof(double) == _real_py_dtype.itemsize == 8 245 246cdef int _complex_typenum = np.NPY_COMPLEX128 247cdef object _complex_py_dtype = np.dtype(np.complex128) 248assert 2 * sizeof(double) == _complex_py_dtype.itemsize == 16 249 250cdef _require_1d_integer(a, dtype): 251 dtype = np.dtype(dtype) 252 if a.dtype.itemsize != dtype.itemsize: 253 warnings.warn("array contains %s bit integers; " 254 "but %s bit integers are needed; " 255 "slowing down due to converting" 256 % (a.dtype.itemsize * 8, 257 dtype.itemsize * 8), 258 CholmodTypeConversionWarning) 259 a = np.ascontiguousarray(a, dtype=dtype) 260 assert a.dtype == dtype 261 assert a.ndim == 1 262 return a 263 264########## 265# Cholmod -> Python conversion: 266########## 267 268cdef int _np_typenum_for_data(int xtype): 269 if xtype == CHOLMOD_COMPLEX: 270 return _complex_typenum 271 elif xtype == CHOLMOD_REAL: 272 return _real_typenum 273 else: 274 raise CholmodError("cholmod->numpy type conversion failed") 275 276cdef type _np_dtype_for_data(int xtype): 277 return np.PyArray_TypeObjectFromType(_np_typenum_for_data(xtype)) 278 279cdef int _np_typenum_for_index(int itype): 280 if itype == CHOLMOD_INT: 281 return _integer_typenum 282 elif itype == CHOLMOD_LONG: 283 return _long_typenum 284 else: 285 raise CholmodError("cholmod->numpy type conversion failed") 286 287cdef type _np_dtype_for_index(int itype): 288 return np.PyArray_TypeObjectFromType(_np_typenum_for_index(itype)) 289 290cdef class _CholmodSparseDestructor: 291 """This is a destructor for NumPy arrays based on sparse data of Cholmod. 292 Use this only once for each Cholmod sparse array. Otherwise memory will be 293 freed multiple times.""" 294 cdef cholmod_sparse* _sparse 295 cdef Common _common 296 297 cdef init(self, cholmod_sparse* m, Common common): 298 assert m is not NULL 299 assert common is not None 300 self._sparse = m 301 self._common = common 302 303 def __dealloc__(self): 304 if self._common._use_long: 305 cholmod_c_free_sparse = cholmod_l_free_sparse 306 else: 307 cholmod_c_free_sparse = cholmod_free_sparse 308 cholmod_c_free_sparse(&self._sparse, &self._common._common) 309 310cdef _cholmod_sparse_to_scipy_sparse(cholmod_sparse * m, Common common): 311 """Build a scipy.sparse.csc_matrix that's a view onto m, with a 'base' with 312 appropriate destructor. 'm' must have been allocated by cholmod.""" 313 314 # This is a little tricky: We build 3 arrays, views on each part of the 315 # cholmod_dense object. They all have the same _CholmodSparseDestructor 316 # object as base. So none of them will be deallocated until they have all 317 # become unused. Then those are built into a csc_matrix. 318 319 # init destructor for cholmod data 320 cdef _CholmodSparseDestructor base = _CholmodSparseDestructor() 321 base.init(m, common) 322 323 # convert to NumPy arrays 324 assert m.itype == common._common.itype 325 cdef np.ndarray indptr = np.PyArray_SimpleNewFromData( 326 1, [m.ncol + 1], _np_typenum_for_index(m.itype), m.p) 327 cdef np.ndarray indices = np.PyArray_SimpleNewFromData( 328 1, [m.nzmax], _np_typenum_for_index(m.itype), m.i) 329 cdef np.ndarray data = np.PyArray_SimpleNewFromData( 330 1, [m.nzmax], _np_typenum_for_data(m.xtype), m.x) 331 332 # set destructor and check if writeable 333 for array in (indptr, indices, data): 334 np.set_array_base(array, base) 335 assert np.PyArray_ISWRITEABLE(indptr) 336 337 # return sparse matrix 338 return sparse.csc_matrix((data, indices, indptr), shape=(m.nrow, m.ncol)) 339 340cdef class _CholmodDenseDestructor: 341 """This is a destructor for NumPy arrays based on dense data of Cholmod. 342 Use this only once for each Cholmod dense array. Otherwise memory will be 343 freed multiple times.""" 344 cdef cholmod_dense* _dense 345 cdef Common _common 346 347 cdef init(self, cholmod_dense* m, Common common): 348 assert m is not NULL 349 assert common is not None 350 self._dense = m 351 self._common = common 352 353 def __dealloc__(self): 354 if self._common._use_long: 355 cholmod_c_free_dense = cholmod_l_free_dense 356 else: 357 cholmod_c_free_dense = cholmod_free_dense 358 cholmod_c_free_dense(&self._dense, &self._common._common) 359 360cdef _cholmod_dense_to_numpy_array(cholmod_dense* m, Common common): 361 """Converts Cholmod dense array to NumPy array. 362 Use this only once for each Cholmod dense array. Otherwise memory will be 363 freed multiple times.""" 364 # init destructor for cholmod data 365 cdef _CholmodDenseDestructor base = _CholmodDenseDestructor() 366 base.init(m, common) 367 # convert cholmod array to numpy array 368 cdef np.ndarray array = np.PyArray_SimpleNewFromData( 369 1, [m.ncol * m.nrow], _np_typenum_for_data(m.xtype), m.x) 370 # set destructor that frees the memory as base 371 np.set_array_base(array, base) 372 # reshape array 373 array = array.reshape((m.ncol, m.nrow)).T 374 # check if array is writeable 375 assert np.PyArray_ISWRITEABLE(array) 376 # return array 377 return array 378 379cdef void _error_handler( 380 int status, const char * file, int line, const char * msg) except * with gil: 381 if status == CHOLMOD_OK: 382 return 383 full_msg = "{}:{:d}: {} (code {:d})".format(file.decode(), line, msg.decode(), status) 384 # known errors 385 if status == CHOLMOD_NOT_POSDEF: 386 raise CholmodNotPositiveDefiniteError(full_msg) 387 elif status == CHOLMOD_NOT_INSTALLED: 388 raise CholmodNotInstalledError(full_msg) 389 elif status == CHOLMOD_OUT_OF_MEMORY: 390 raise CholmodOutOfMemoryError(full_msg) 391 elif status == CHOLMOD_TOO_LARGE: 392 raise CholmodTooLargeError(full_msg) 393 elif status == CHOLMOD_INVALID: 394 raise CholmodInvalidError(full_msg) 395 elif status == CHOLMOD_GPU_PROBLEM: 396 raise CholmodGpuProblemError(full_msg) 397 # unknown errors 398 if status < 0: 399 raise CholmodError(full_msg) 400 # warnings 401 else: 402 warnings.warn(full_msg, CholmodWarning) 403 404def _check_for_csc(m): 405 if not sparse.isspmatrix_csc(m): 406 warnings.warn("converting matrix of class %s to CSC format" 407 % (m.__class__.__name__,), 408 CholmodTypeConversionWarning) 409 m = m.tocsc() 410 assert sparse.isspmatrix_csc(m) 411 return m 412 413cdef class Common: 414 cdef cholmod_common _common 415 cdef int _complex 416 cdef int _xtype 417 cdef int _use_long 418 419 def __cinit__(self, _complex, _use_long): 420 self._complex = _complex 421 if self._complex: 422 self._xtype = CHOLMOD_COMPLEX 423 else: 424 self._xtype = CHOLMOD_REAL 425 self._use_long = _use_long 426 if self._use_long: 427 cholmod_c_start = cholmod_l_start 428 else: 429 cholmod_c_start = cholmod_start 430 cholmod_c_start(&self._common) 431 assert (_use_long == 0 and self._common.itype == CHOLMOD_INT) \ 432 or (_use_long == 1 and self._common.itype == CHOLMOD_LONG) 433 self._common.print = 0 434 self._common.error_handler = ( 435 <void (*)(int, const char *, int, const char *)>_error_handler) 436 437 def __dealloc__(self): 438 if self._use_long: 439 cholmod_c_finish = cholmod_l_finish 440 else: 441 cholmod_c_finish = cholmod_finish 442 cholmod_c_finish(&self._common) 443 444 # Debugging: 445 def _print(self): 446 if self._use_long: 447 cholmod_c_check_common = cholmod_l_check_common 448 cholmod_c_print_common = cholmod_l_print_common 449 else: 450 cholmod_c_check_common = cholmod_check_common 451 cholmod_c_print_common = cholmod_print_common 452 453 print(cholmod_c_check_common(&self._common)) 454 name = repr(self).encode() 455 return cholmod_c_print_common(name, &self._common) 456 457 def _print_sparse(self, name, symmetric, matrix): 458 if self._use_long: 459 cholmod_c_check_sparse = cholmod_l_check_sparse 460 cholmod_c_print_sparse = cholmod_l_print_sparse 461 else: 462 cholmod_c_check_sparse = cholmod_check_sparse 463 cholmod_c_print_sparse = cholmod_print_sparse 464 cdef cholmod_sparse m 465 cdef object ref = self._init_view_sparse(&m, matrix, symmetric) 466 print(cholmod_c_check_sparse(&m, &self._common)) 467 return cholmod_c_print_sparse(&m, name, &self._common) 468 469 def _print_dense(self, name, matrix): 470 if self._use_long: 471 cholmod_c_check_dense = cholmod_l_check_dense 472 cholmod_c_print_dense = cholmod_l_print_dense 473 else: 474 cholmod_c_check_dense = cholmod_check_dense 475 cholmod_c_print_dense = cholmod_print_dense 476 cdef cholmod_dense m 477 cdef object ref = self._init_view_dense(&m, matrix) 478 print(cholmod_c_check_dense(&m, &self._common)) 479 return cholmod_c_print_dense(&m, name, &self._common) 480 481 ########## 482 # Python -> Cholmod conversion: 483 ########## 484 cdef np.ndarray _cast(self, np.ndarray arr): 485 if not issubclass(arr.dtype.type, np.number): 486 raise CholmodError("non-numeric dtype %s" % (arr.dtype,)) 487 if self._complex: 488 # All numeric types can be upcast to complex: 489 return np.asfortranarray(arr, dtype=_complex_py_dtype) 490 else: 491 # Refuse to downcast complex types to real: 492 if issubclass(arr.dtype.type, np.complexfloating): 493 raise CholmodError("inconsistent use of complex array") 494 else: 495 return np.asfortranarray(arr, dtype=_real_py_dtype) 496 497 # Some memory allocated for the init'd sparse matrix is refcounted by the 498 # returned object; do not let it be garbage collected as long as you want 499 # to use the matrix. 500 cdef object _init_view_sparse(self, cholmod_sparse *out, m, symmetric): 501 if symmetric and m.shape[0] != m.shape[1]: 502 raise CholmodError("supposedly symmetric matrix is not square!") 503 m = _check_for_csc(m) 504 m.sort_indices() 505 cdef np.ndarray indptr = _require_1d_integer(m.indptr, _np_dtype_for_index(self._common.itype)) 506 cdef np.ndarray indices = _require_1d_integer(m.indices, _np_dtype_for_index(self._common.itype)) 507 cdef np.ndarray data = self._cast(m.data) 508 out.nrow, out.ncol = m.shape 509 out.nzmax = m.nnz 510 out.p = indptr.data 511 out.i = indices.data 512 out.x = data.data 513 if symmetric: 514 out.stype = -1 515 else: 516 out.stype = 0 517 out.itype = self._common.itype 518 out.dtype = CHOLMOD_DOUBLE 519 out.xtype = self._xtype 520 out.sorted = 1 521 out.packed = 1 522 return m, indptr, indices, data 523 524 # Some memory allocated for the init'd dense matrix is refcounted by the 525 # returned object; do not let it be garbage collected as long as you want 526 # to use the matrix. 527 cdef object _init_view_dense(self, cholmod_dense *out, np.ndarray m): 528 if m.ndim != 2: 529 raise CholmodError("array has %s dimensions (expected 2)" % m.ndim) 530 m = self._cast(m) 531 # The leading dimension is equal to m.shape[0] because `_cast` ensures 532 # that m is Fortran-contiguous. It is not necessarily equal to 533 # `m.strides[1] // m.itemsize` when relaxed stride checking is on. 534 out.nrow = out.d = m.shape[0] 535 out.ncol = m.shape[1] 536 out.nzmax = m.size 537 out.x = m.data 538 out.dtype = CHOLMOD_DOUBLE 539 out.xtype = self._xtype 540 return m 541 542cdef object factor_secret_handshake = object() 543 544cdef class Factor: 545 """This class represents a Cholesky decomposition with a particular 546 fill-reducing permutation. It cannot be instantiated directly; see 547 :func:`analyze` and :func:`cholesky`, both of which return objects of type 548 Factor. 549 """ 550 551 cdef readonly Common _common 552 cdef cholmod_factor * _factor 553 554 def __init__(self, handshake): 555 if handshake is not factor_secret_handshake: 556 raise CholmodError("Factor may not be constructed directly; use analyze()") 557 558 def __dealloc__(self): 559 if self._common._use_long: 560 cholmod_c_free_factor = cholmod_l_free_factor 561 else: 562 cholmod_c_free_factor = cholmod_free_factor 563 cholmod_c_free_factor(&self._factor, &self._common._common) 564 565 def _print(self): 566 if self._common._use_long: 567 cholmod_c_check_factor = cholmod_l_check_factor 568 cholmod_c_print_factor = cholmod_l_print_factor 569 else: 570 cholmod_c_check_factor = cholmod_check_factor 571 cholmod_c_print_factor = cholmod_print_factor 572 print(cholmod_c_check_factor(self._factor, &self._common._common)) 573 name = repr(self).encode() 574 return cholmod_c_print_factor(self._factor, name, &self._common._common) 575 576 def cholesky_inplace(self, A, beta=0): 577 """Updates this Factor so that it represents the Cholesky 578 decomposition of :math:`A + \\beta I`, rather than whatever it 579 contained before. 580 581 :math:`A` must have the same pattern of non-zeros as the matrix used 582 to create this factor originally.""" 583 return self._cholesky_inplace(A, True, beta=beta) 584 585 def cholesky_AAt_inplace(self, A, beta=0): 586 """The same as :meth:`cholesky_inplace`, except it factors :math:`AA' 587 + \\beta I` instead of :math:`A + \\beta I`.""" 588 return self._cholesky_inplace(A, False, beta=beta) 589 590 def _cholesky_inplace(self, A, symmetric, beta=0, **kwargs): 591 cdef cholmod_sparse c_A 592 cdef object ref = self._common._init_view_sparse(&c_A, A, symmetric) 593 try: 594 if self._common._use_long: 595 cholmod_l_factorize_p(&c_A, [beta, 0], NULL, 0, 596 self._factor, &self._common._common) 597 else: 598 cholmod_factorize_p(&c_A, [beta, 0], NULL, 0, 599 self._factor, &self._common._common) 600 except CholmodNotPositiveDefiniteError as e: 601 e.factor = self 602 e.column = self._factor.minor 603 raise 604 assert self._common._common.status == CHOLMOD_OK 605 606 def copy(self): 607 """Copies the current :class:`Factor`. 608 609 :returns: A new :class:`Factor` object.""" 610 if self._common._use_long: 611 cholmod_c_copy_factor = cholmod_l_copy_factor 612 else: 613 cholmod_c_copy_factor = cholmod_copy_factor 614 615 cdef cholmod_factor * c_clone = cholmod_c_copy_factor(self._factor, 616 &self._common._common) 617 assert c_clone 618 cdef Factor clone = Factor(factor_secret_handshake) 619 clone._common = self._common 620 clone._factor = c_clone 621 assert self._factor.itype == clone._factor.itype == self._common._common.itype 622 assert self._factor.xtype == clone._factor.xtype 623 return clone 624 625 def cholesky(self, A, beta=0): 626 """The same as :meth:`cholesky_inplace` except that it first creates 627 a copy of the current :class:`Factor` and modifes the copy. 628 629 :returns: The new :class:`Factor` object.""" 630 clone = self.copy() 631 clone.cholesky_inplace(A, beta=beta) 632 return clone 633 634 def cholesky_AAt(self, A, beta=0): 635 """The same as :meth:`cholesky_AAt_inplace` except that it first 636 creates a copy of the current :class:`Factor` and modifes the copy. 637 638 :returns: The new :class:`Factor` object.""" 639 clone = self.copy() 640 clone.cholesky_AAt_inplace(A, beta=beta) 641 return clone 642 643 def update_inplace(self, C, bint subtract=False): 644 """Incremental building of :math:`AA'` decompositions. 645 646 Updates this factor so that instead of representing the decomposition 647 of :math:`A` (:math:`AA'`), it computes the decomposition of 648 :math:`A + CC'` (:math:`AA' + CC'`) for ``subtract=False`` which is the 649 default, or :math:`A - CC'` (:math:`AA' - CC'`) for 650 ``subtract=True``. This method does not require that the 651 :class:`Factor` was created with :func:`cholesky_AAt`, though that 652 is the common case. 653 654 The usual use for this is to factor AA' when A has a large number of 655 columns, or those columns become available incrementally. Instead of 656 loading all of A into memory, one can load in 'strips' of columns and 657 pass them to this method one at a time. 658 659 Note that no fill-reduction analysis is done; whatever permutation was 660 chosen by the initial call to :func:`analyze` will be used regardless 661 of the pattern of non-zeros in C.""" 662 # permute C 663 if self._common._use_long: 664 cholmod_c_updown = cholmod_l_updown 665 cholmod_c_free_sparse = cholmod_l_free_sparse 666 else: 667 cholmod_c_updown = cholmod_updown 668 cholmod_c_free_sparse = cholmod_free_sparse 669 670 cdef cholmod_sparse c_C 671 cdef object ref = self._common._init_view_sparse(&c_C, C, False) 672 cdef cholmod_sparse * C_perm 673 674 if self._common._use_long: 675 C_perm = cholmod_l_submatrix( 676 &c_C, <SuiteSparse_long *> self._factor.Perm, self._factor.n, NULL, -1, 677 True, True, &self._common._common) 678 else: 679 C_perm = cholmod_submatrix( 680 &c_C, <int *> self._factor.Perm, self._factor.n, NULL, -1, 681 True, True, &self._common._common) 682 assert C_perm 683 try: 684 cholmod_c_updown(not subtract, C_perm, self._factor, 685 &self._common._common) 686 finally: 687 cholmod_c_free_sparse(&C_perm, &self._common._common) 688 689 # Everything below here will fail for matrices that were only analyzed, 690 # not factorized. 691 def P(self): 692 """Returns the fill-reducing permutation P, as a vector of indices. 693 694 The decomposition :math:`LL'` or :math:`LDL'` is of:: 695 696 A[P[:, np.newaxis], P[np.newaxis, :]] 697 698 (or similar for AA').""" 699 if self._factor.Perm is NULL: 700 raise CholmodError("you must analyze a matrix first") 701 assert self._factor.itype == self._common._common.itype 702 703 cdef np.ndarray out = np.PyArray_SimpleNewFromData( 704 1, [self._factor.n], _np_typenum_for_index(self._factor.itype), self._factor.Perm) 705 np.set_array_base(out, self) 706 707 return out 708 709 def _ensure_L_or_LD_inplace(self, want_L): 710 # In CHOLMOD, supernodal factorizations are always LL'. If we request 711 # to change to a supernodal LDL' factorization, cholmod_change_factor 712 # will silently do nothing! So we can only stay supernodal when LL' is 713 # requested: 714 if self._common._use_long: 715 cholmod_c_change_factor = cholmod_l_change_factor 716 else: 717 cholmod_c_change_factor = cholmod_change_factor 718 want_super = self._factor.is_super and want_L 719 cholmod_c_change_factor(self._factor.xtype, 720 want_L, # to_ll 721 want_super, 722 True, # to_packed 723 self._factor.is_monotonic, 724 self._factor, 725 &self._common._common) 726 assert bool(self._factor.is_ll) == want_L 727 728 def _L_or_LD(self, want_L): 729 if self._common._use_long: 730 cholmod_c_factor_to_sparse = cholmod_l_factor_to_sparse 731 else: 732 cholmod_c_factor_to_sparse = cholmod_factor_to_sparse 733 cdef Factor f = self.copy() 734 cdef cholmod_sparse * l 735 f._ensure_L_or_LD_inplace(want_L) 736 l = cholmod_c_factor_to_sparse(f._factor, 737 &f._common._common) 738 assert l 739 return _cholmod_sparse_to_scipy_sparse(l, self._common) 740 741 def D(self): 742 """Converts this factorization to the style 743 744 .. math:: LDL' = PAP' 745 746 or 747 748 .. math:: LDL' = PAA'P' 749 750 and then returns the diagonal matrix D *as a 1d vector*. 751 752 .. note:: This method uses an efficient implementation that extracts 753 the diagonal D directly from CHOLMOD's internal 754 representation. It never makes a copy of the factor matrices, or 755 actually converts a full `LL'` factorization into an `LDL'` 756 factorization just to extract `D`. 757 758 """ 759 760 if self._factor.xtype == CHOLMOD_PATTERN: 761 raise CholmodError("cannot extract diagonal from a symbolic " 762 "factor; call a cholesky*() method first.") 763 764 cdef np.ndarray x = np.PyArray_SimpleNewFromData( 765 1, [self._factor.xsize if self._factor.is_super else self._factor.nzmax], 766 _np_typenum_for_data(self._factor.xtype), self._factor.x) 767 768 cdef size_t i 769 cdef np.npy_intp n 770 cdef SuiteSparse_long * l_super_ = <SuiteSparse_long *> self._factor.super_ 771 cdef SuiteSparse_long * l_pi = <SuiteSparse_long *> self._factor.pi 772 cdef SuiteSparse_long * l_px = <SuiteSparse_long *> self._factor.px 773 cdef int * i_super_ = <int *> self._factor.super_ 774 cdef int * i_pi = <int *> self._factor.pi 775 cdef int * i_px = <int *> self._factor.px 776 777 if self._factor.is_super: 778 # This is a supernodal factorization, which is stored as a bunch 779 # of dense, lower-triangular, column-major arrays packed into the 780 # x vector. This is not documented in the CHOLMOD user-guide, or 781 # anywhere else as far as I can tell; I got the details from 782 # CVXOPT's C/cholmod.c. 783 d = np.empty(self._factor.n, dtype=_np_dtype_for_data(self._factor.xtype)) 784 filled = 0 785 for i in xrange(self._factor.nsuper): 786 if self._common._use_long: 787 ncols = l_super_[i + 1] - l_super_[i] 788 nrows = l_pi[i + 1] - l_pi[i] 789 px_i = l_px[i] 790 else: 791 ncols = i_super_[i + 1] - i_super_[i] 792 nrows = i_pi[i + 1] - i_pi[i] 793 px_i = i_px[i] 794 d[filled:filled + ncols] = x[px_i 795 :px_i + nrows * ncols 796 :nrows + 1] 797 filled += ncols 798 else: 799 # This is a simplicial factorization, which is simply stored as a 800 # sparse CSC matrix in x, p, i. We want the diagonal, which is 801 # just the first entry in each column; p gives the offsets in x to 802 # the beginning of each column. 803 804 # The ->p array actually has n+1 entries, but only the first n 805 # entries actually point to real columns (the last entry is a 806 # sentinel), so we just create a view onto those: 807 assert self._factor.itype == self._common._common.itype 808 p = np.PyArray_SimpleNewFromData( 809 1, [self._factor.n], _np_typenum_for_index(self._factor.itype), self._factor.p) 810 811 d = x[p] 812 if self._factor.is_ll: 813 return d ** 2 814 else: 815 return d 816 817 def L(self): 818 """If necessary, converts this factorization to the style 819 820 .. math:: LL' = PAP' 821 822 or 823 824 .. math:: LL' = PAA'P' 825 826 and then returns the sparse lower-triangular matrix L. 827 828 .. warning:: The L matrix returned by this method and the one returned 829 by :meth:`L_D` are different! 830 """ 831 return self._L_or_LD(True) 832 833 def LD(self): 834 """If necessary, converts this factorization to the style 835 836 .. math:: LDL' = PAP' 837 838 or 839 840 .. math:: LDL' = PAA'P' 841 842 and then returns a sparse lower-triangular matrix "LD", which contains 843 the D matrix on its diagonal, plus the below-diagonal part of L (the 844 actual diagonal of L is all-ones). 845 846 See :meth:`L_D` for a more convenient interface.""" 847 return self._L_or_LD(False) 848 849 def L_D(self): 850 """If necessary, converts this factorization to the style 851 852 .. math:: LDL' = PAP' 853 854 or 855 856 .. math:: LDL' = PAA'P' 857 858 and then returns the pair (L, D) where L is a sparse lower-triangular 859 matrix and D is a sparse diagonal matrix. 860 861 .. warning:: The L matrix returned by this method and the one returned 862 by :meth:`L` are different! 863 """ 864 ld = self.LD() 865 l = sparse.tril(ld, -1) + sparse.eye(*ld.shape) 866 d = sparse.dia_matrix((ld.diagonal(), [0]), shape=ld.shape) 867 return (l, d) 868 869 def solve_A(self, b): 870 """ Solves a linear system. 871 872 :param b: right-hand-side 873 874 :returns: math:`x`, where :math:`Ax = b` (or :math:`AA'x = b`, if 875 you used :func:`cholesky_AAt`). 876 877 :meth:`__call__` is an alias for this function, i.e., you can simply 878 call the :class:`Factor` object like a function to solve :math:`Ax = 879 b`.""" 880 return self._solve(b, CHOLMOD_A) 881 882 def __call__(self, b): 883 "Alias for :meth:`solve_A`." 884 return self.solve_A(b) 885 886 def solve_LDLt(self, b): 887 """ Solves a linear system. 888 889 :param b: right-hand-side 890 891 :returns: math:`x`, where :math:`LDL'x = b`. 892 893 (This is different from :meth:`solve_A` because it does not correct 894 for the fill-reducing permutation.)""" 895 return self._solve(b, CHOLMOD_LDLt) 896 897 def solve_LD(self, b): 898 """ Solves a linear system. 899 900 :param b: right-hand-side 901 902 :returns: math:`x`, where :math:`LDx = b`.""" 903 self._ensure_L_or_LD_inplace(False) 904 return self._solve(b, CHOLMOD_LD) 905 906 def solve_DLt(self, b): 907 """ Solves a linear system. 908 909 :param b: right-hand-side 910 911 :returns: math:`x`, where :math:`DL'x = b`.""" 912 913 self._ensure_L_or_LD_inplace(False) 914 return self._solve(b, CHOLMOD_DLt) 915 916 def solve_L(self, b, use_LDLt_decomposition=True): 917 """ Solves a linear system. 918 919 :param b: right-hand-side 920 921 :param use_LDLt_decomposition: If True, use the `L` of the `LDL'` 922 decomposition. If False, use the `L` of the `LL'` decomposition. 923 924 :returns: math:`x`, where :math:`Lx = b`.""" 925 self._ensure_L_or_LD_inplace(not use_LDLt_decomposition) 926 return self._solve(b, CHOLMOD_L) 927 928 def solve_Lt(self, b, use_LDLt_decomposition=True): 929 """ Solves a linear system. 930 931 :param b: right-hand-side 932 933 :param use_LDLt_decomposition: If True, use the `L` of the `LDL'` 934 decomposition. If False, use the `L` of the `LL'` decomposition. 935 936 :returns: math:`x`, where :math:`L'x = b`.""" 937 self._ensure_L_or_LD_inplace(not use_LDLt_decomposition) 938 return self._solve(b, CHOLMOD_Lt) 939 940 def solve_D(self, b): 941 "Returns :math:`x`, where :math:`Dx = b`." 942 return self._solve(b, CHOLMOD_D) 943 944 # CHOLMOD API is quite confusing here -- unlike all the other solve 945 # magic constants, CHOLMOD_P and CHOLMOD_Pt actually apply the matrix to b 946 # rather than performing a matrix solve. Basically their names are 947 # backwards... therefore let's call the functions `apply_P` and `apply_Pt`. 948 def apply_P(self, b): 949 "Returns :math:`x`, where :math:`x = Pb`." 950 return self._solve(b, CHOLMOD_P) 951 952 def apply_Pt(self, b): 953 "Returns :math:`x`, where :math:`x = P'b`." 954 return self._solve(b, CHOLMOD_Pt) 955 956 def _solve(self, b, system): 957 if sparse.issparse(b): 958 return self._solve_sparse(b, system) 959 else: 960 return self._solve_dense(b, system) 961 962 def _solve_sparse(self, b, system): 963 if self._common._use_long: 964 cholmod_c_spsolve = cholmod_l_spsolve 965 else: 966 cholmod_c_spsolve = cholmod_spsolve 967 cdef cholmod_sparse c_b 968 cdef object ref = self._common._init_view_sparse(&c_b, b, False) 969 cdef cholmod_sparse *out = cholmod_c_spsolve( 970 system, self._factor, &c_b, &self._common._common) 971 return _cholmod_sparse_to_scipy_sparse(out, self._common) 972 973 def _solve_dense(self, b, system): 974 if self._common._use_long: 975 cholmod_c_solve = cholmod_l_solve 976 else: 977 cholmod_c_solve = cholmod_solve 978 979 b = np.asarray(b) 980 ndim = b.ndim 981 if b.ndim == 1: 982 b = b[:, np.newaxis] 983 cdef cholmod_dense c_b 984 cdef object ref = self._common._init_view_dense(&c_b, b) 985 cdef cholmod_dense *out = cholmod_c_solve( 986 system, self._factor, &c_b, &self._common._common) 987 py_out = _cholmod_dense_to_numpy_array(out, self._common) 988 if ndim == 1: 989 py_out = py_out[:, 0] 990 return py_out 991 992 def slogdet(self): 993 """Computes the log-determinant of the matrix A, with the same API as 994 :meth:`numpy.linalg.slogdet`. 995 996 This returns a tuple `(sign, logdet)`, where `sign` is always the 997 number 1.0 (because the determinant of a positive-definite matrix is 998 always a positive real number), and `logdet` is the (natural) 999 logarithm of the determinant of the matrix A. 1000 1001 .. versionadded:: 0.2 1002 """ 1003 return (1.0, self.logdet()) 1004 1005 def logdet(self): 1006 """Computes the (natural) log of the determinant of the matrix A. 1007 1008 If `f` is a factor, then `f.logdet()` is equivalent to 1009 `np.sum(np.log(f.D()))`. 1010 1011 .. versionadded:: 0.2 1012 """ 1013 return np.sum(np.log(self.D())) 1014 1015 def det(self): 1016 """Computes the determinant of the matrix A. 1017 1018 Consider using :meth:`logdet` instead, for improved numerical 1019 stability. (In particular, determinants are often prone to problems 1020 with underflow or overflow.) 1021 1022 .. versionadded:: 0.2 1023 """ 1024 return np.exp(self.logdet()) 1025 1026 def inv(self): 1027 """Returns the inverse of the matrix A, as a sparse (CSC) matrix. 1028 1029 .. warning:: For most purposes, it is better to use :meth:`solve` 1030 instead of computing the inverse explicitly. That is, the 1031 following two pieces of code produce identical results:: 1032 1033 x = f.solve(b) 1034 x = f.inv() * b # DON'T DO THIS! 1035 1036 But the first line is both faster and produces more accurate 1037 results. 1038 1039 Sometimes, though, you really do need the inverse explicitly (e.g., 1040 for calculating standard errors in least squares regression), so if 1041 that's your situation, here you go. 1042 1043 .. versionadded:: 0.2 1044 """ 1045 return self(sparse.eye(self._factor.n, self._factor.n, 1046 dtype=_np_dtype_for_data(self._factor.xtype), 1047 format="csc")) 1048 1049def analyze(A, mode="auto", ordering_method="default", use_long=None): 1050 """Computes the optimal fill-reducing permutation for the symmetric matrix 1051 A, but does *not* factor it (i.e., it performs a "symbolic Cholesky 1052 decomposition"). This function ignores the actual contents of the matrix 1053 A. All it cares about are (1) which entries are non-zero, and (2) whether 1054 A has real or complex type. 1055 1056 :param A: The matrix to be analyzed. 1057 1058 :param mode: Specifies which algorithm should be used to (eventually) 1059 compute the Cholesky decomposition -- one of "simplicial", "supernodal", 1060 or "auto". See the CHOLMOD documentation for details on how "auto" chooses 1061 the algorithm to be used. 1062 1063 :param ordering_method: Specifies which ordering algorithm should be used to 1064 (eventually) order the matrix A -- one of "natural", "amd", "metis", 1065 "nesdis", "colamd", "default" and "best". "natural" means no permutation. 1066 See the CHOLMOD documentation for more details. 1067 1068 :param use_long: Specifies if the long type (64 bit) or the int type 1069 (32 bit) should be used for the indices of the sparse matrices. If 1070 use_long is None try to estimate if long type is needed. 1071 1072 :returns: A :class:`Factor` object representing the analysis. Many 1073 operations on this object will fail, because it does not yet hold a full 1074 decomposition. Use :meth:`Factor.cholesky_inplace` (or similar) to 1075 actually factor a matrix. 1076 """ 1077 return _analyze(A, True, mode=mode, ordering_method=ordering_method, use_long=use_long) 1078 1079def analyze_AAt(A, mode="auto", ordering_method="default", use_long=None): 1080 """Computes the optimal fill-reducing permutation for the symmetric matrix 1081 :math:`AA'`, but does *not* factor it (i.e., it performs a "symbolic 1082 Cholesky decomposition"). This function ignores the actual contents of the 1083 matrix A. All it cares about are (1) which entries are non-zero, and (2) 1084 whether A has real or complex type. 1085 1086 :param A: The matrix to be analyzed. 1087 1088 :param mode: Specifies which algorithm should be used to (eventually) 1089 compute the Cholesky decomposition -- one of "simplicial", "supernodal", 1090 or "auto". See the CHOLMOD documentation for details on how "auto" chooses 1091 the algorithm to be used. 1092 1093 :param ordering_method: Specifies which ordering algorithm should be used to 1094 (eventually) order the matrix A -- one of "natural", "amd", "metis", 1095 "nesdis", "colamd", "default" and "best". "natural" means no permutation. 1096 See the CHOLMOD documentation for more details. 1097 1098 :param use_long: Specifies if the long type (64 bit) or the int type 1099 (32 bit) should be used for the indices of the sparse matrices. If 1100 use_long is None try to estimate if long type is needed. 1101 1102 :returns: A :class:`Factor` object representing the analysis. Many 1103 operations on this object will fail, because it does not yet hold a full 1104 decomposition. Use :meth:`Factor.cholesky_AAt_inplace` (or similar) to 1105 actually factor a matrix. 1106 """ 1107 return _analyze(A, False, mode=mode, ordering_method=ordering_method, use_long=use_long) 1108 1109_modes = { 1110 "simplicial": CHOLMOD_SIMPLICIAL, 1111 "supernodal": CHOLMOD_SUPERNODAL, 1112 "auto": CHOLMOD_AUTO, 1113 } 1114_ordering_methods = { 1115 "natural": CHOLMOD_NATURAL, 1116 "amd": CHOLMOD_AMD, 1117 "metis": CHOLMOD_METIS, 1118 "nesdis": CHOLMOD_NESDIS, 1119 "colamd": CHOLMOD_COLAMD, 1120 "default": None, 1121 "best": None, 1122} 1123 1124def _analyze(A, symmetric, mode, ordering_method="default", use_long=None): 1125 A = _check_for_csc(A) 1126 if use_long is None: 1127 assert A.indices.dtype == A.indptr.dtype 1128 use_long = A.indices.dtype == np.int64 1129 elif not use_long: 1130 INT_MAX = np.iinfo(np.int32).max 1131 if (A.nnz > INT_MAX or np.any(np.array(A.shape) > INT_MAX)): 1132 warnings.warn('Problem to large for int, switching to long.') 1133 use_long = True 1134 cdef Common common = Common(issubclass(A.dtype.type, np.complexfloating), use_long) 1135 cdef cholmod_sparse c_A 1136 cdef object ref = common._init_view_sparse(&c_A, A, symmetric) 1137 if mode in _modes: 1138 common._common.supernodal = _modes[mode] 1139 else: 1140 raise CholmodError("Unknown mode '%s', must be one of %s" % 1141 (mode, ", ".join(_modes.keys()))) 1142 1143 if ordering_method in _ordering_methods: 1144 if ordering_method == "default": 1145 common._common.nmethods = 0 1146 elif ordering_method == "best": 1147 common._common.nmethods = 9 1148 else: 1149 common._common.nmethods = 1 1150 common._common.method [0].ordering = _ordering_methods[ordering_method] 1151 common._common.postorder = ordering_method != "natural" 1152 elif ordering_method is not None: 1153 raise CholmodError, ("Unknown ordering method '%s', must be one of %s" 1154 % (ordering_method, ", ".join(_ordering_methods.keys()))) 1155 1156 if common._use_long: 1157 cholmod_c_analyze = cholmod_l_analyze 1158 else: 1159 cholmod_c_analyze = cholmod_analyze 1160 cdef cholmod_factor *c_f = cholmod_c_analyze(&c_A, &common._common) 1161 if c_f is NULL: 1162 raise CholmodError("Error in cholmod_analyze") 1163 1164 cdef Factor f = Factor(factor_secret_handshake) 1165 f._common = common 1166 f._factor = c_f 1167 return f 1168 1169def cholesky(A, beta=0, mode="auto", ordering_method="default", use_long=None): 1170 """Computes the fill-reducing Cholesky decomposition of 1171 1172 .. math:: A + \\beta I 1173 1174 where ``A`` is a sparse, symmetric, positive-definite matrix, preferably 1175 in CSC format, and ``beta`` is any real scalar (usually 0 or 1). (And 1176 :math:`I` denotes the identity matrix.) 1177 1178 Only the lower triangular part of ``A`` is used. 1179 1180 ``mode`` is passed to :func:`analyze`. 1181 1182 ``ordering_method`` is passed to :func:`analyze`. 1183 1184 ``use_long`` is passed to :func:`analyze`. 1185 1186 :returns: A :class:`Factor` object represented the decomposition. 1187 """ 1188 return _cholesky(A, True, beta=beta, mode=mode, ordering_method=ordering_method, use_long=use_long) 1189 1190def cholesky_AAt(A, beta=0, mode="auto", ordering_method="default", use_long=None): 1191 """Computes the fill-reducing Cholesky decomposition of 1192 1193 .. math:: AA' + \\beta I 1194 1195 where ``A`` is a sparse matrix, preferably in CSC format, and ``beta`` is 1196 any real scalar (usually 0 or 1). (And :math:`I` denotes the identity 1197 matrix.) 1198 1199 Note that if you are solving a conventional least-squares problem, you 1200 will need to transpose your matrix before calling this function, and 1201 therefore it will be somewhat more efficient to construct your matrix in 1202 CSR format (so that its transpose will be in CSC format). 1203 1204 ``mode`` is passed to :func:`analyze_AAt`. 1205 1206 ``ordering_method`` is passed to :func:`analyze_AAt`. 1207 1208 ``use_long`` is passed to :func:`analyze_AAt`. 1209 1210 :returns: A :class:`Factor` object represented the decomposition. 1211 """ 1212 return _cholesky(A, False, beta=beta, mode=mode, ordering_method=ordering_method, use_long=use_long) 1213 1214def _cholesky(A, symmetric, beta, mode, ordering_method="default", use_long=None): 1215 f = _analyze(A, symmetric, mode=mode, ordering_method=ordering_method, use_long=use_long) 1216 f._cholesky_inplace(A, symmetric, beta=beta) 1217 return f 1218 1219__all__ = ["analyze", "analyze_AAt", "cholesky", "cholesky_AAt"] 1220