1 //------------------------------------------------------------------------------ 2 // SLIP_LU/Include/SLIP_LU.h: user #include file for SLIP_LU. 3 //------------------------------------------------------------------------------ 4 5 // SLIP_LU: (c) 2019-2020, Chris Lourenco, Jinhao Chen, Erick Moreno-Centeno, 6 // Timothy A. Davis, Texas A&M University. All Rights Reserved. See 7 // SLIP_LU/License for the license. 8 9 //------------------------------------------------------------------------------ 10 11 #ifndef SLIP_LU_H 12 #define SLIP_LU_H 13 14 // This software package exactly solves a sparse system of linear equations 15 // using the SLIP LU factorization. This code accompanies the paper (submitted 16 // to ACM Transactions on Mathematical Software): 17 18 // "Algorithm 1xxx: SLIP LU: A Sparse Left-Looking Integer-Preserving LU 19 // Factorization for Exactly Solving Sparse Linear Systems", 20 // C. Lourenco, J. Chen, E. Moreno-Centeno, T. Davis, under submission, 21 // ACM Trans. Mathematical Software. 22 23 // The theory associated with this software can be found in the paper 24 // (published in SIAM journal on matrix analysis and applications): 25 26 // "Exact Solution of Sparse Linear Systems via Left-Looking 27 // Roundoff-Error-Free LU Factorization in Time Proportional to 28 // Arithmetic Work", C. Lourenco, A. R. Escobedo, E. Moreno-Centeno, 29 // T. Davis, SIAM J. Matrix Analysis and Applications. pp 609-638, 30 // vol 40, no 2, 2019. 31 32 // If you use this code, you must first download and install the GMP and 33 // MPFR libraries. GMP and MPFR can be found at: 34 // https://gmplib.org/ 35 // http://www.mpfr.org/ 36 37 // If you use SLIP LU for a publication, we request that you please cite 38 // the above two papers. 39 40 //------------------------------------------------------------------------------ 41 //------------------------------------------------------------------------------ 42 //-------------------------Authors---------------------------------------------- 43 //------------------------------------------------------------------------------ 44 //------------------------------------------------------------------------------ 45 46 // Christopher Lourenco, Jinhao Chen, Erick Moreno-Centeno, and Timothy Davis 47 // 48 49 //------------------------------------------------------------------------------ 50 //------------------------------------------------------------------------------ 51 //-------------------------Contact Information---------------------------------- 52 //------------------------------------------------------------------------------ 53 //------------------------------------------------------------------------------ 54 55 // Please contact Chris Lourenco (chrisjlourenco@gmail.com) 56 // or Tim Davis (timdavis@aldenmath.com, DrTimothyAldenDavis@gmail.com, 57 // davis@tamu.edu) 58 59 //------------------------------------------------------------------------------ 60 //------------------------------------------------------------------------------ 61 //-------------------------Copyright-------------------------------------------- 62 //------------------------------------------------------------------------------ 63 //------------------------------------------------------------------------------ 64 65 // SLIP LU is free software; you can redistribute it and/or modify 66 // it under the terms of either: 67 // 68 // * the GNU Lesser General Public License as published by the 69 // Free Software Foundation; either version 3 of the License, 70 // or (at your option) any later version. 71 // 72 // or 73 // 74 // * the GNU General Public License as published by the Free Software 75 // Foundation; either version 2 of the License, or (at your option) any 76 // later version. 77 // 78 // or both in parallel, as here. 79 // 80 // See license.txt for license info. 81 // 82 // This software is copyright by Christopher Lourenco, Jinhao Chen, Erick 83 // Moreno-Centeno and Timothy A. Davis. All Rights Reserved. 84 // 85 86 //------------------------------------------------------------------------------ 87 //------------------------------------------------------------------------------ 88 //---------------------------DISCLAIMER----------------------------------------- 89 //------------------------------------------------------------------------------ 90 //------------------------------------------------------------------------------ 91 92 // SLIP LU is distributed in the hope that it will be useful, but 93 // WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 94 // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 95 // for more details. 96 97 //------------------------------------------------------------------------------ 98 //------------------------------------------------------------------------------ 99 //--------------------------Summary--------------------------------------------- 100 //------------------------------------------------------------------------------ 101 //------------------------------------------------------------------------------ 102 103 // This software package solves the linear system Ax = b exactly. The input 104 // matrix and right hand side vectors are stored as either integers, double 105 // precision numbers, multiple precision floating points (through the mpfr 106 // library) or as rational numbers (as a collection of numerators and 107 // denominators using the GMP mpq_t data structure). Appropriate routines 108 // within the code transform the input into an integral matrix in compressed 109 // column form. 110 111 // This package computes the factorization PAQ = LDU. Note that we store the 112 // "functional" form of the factorization by only storing L and U. The user 113 // is given some freedom to select the permutation matrices P and Q. The 114 // recommended default settings select Q using the COLAMD column ordering 115 // and select P via a partial pivoting scheme in which the diagonal entry 116 // in column k is selected if it is the same magnitude as the smallest 117 // entry, otherwise the smallest entry is selected as the kth pivot. 118 // Alternative strategies allowed to select Q include the AMD column 119 // ordering or no column permutation (Q=I). For pivots, there are a variety 120 // of potential schemes including traditional partial pivoting, diagonal 121 // pivoting, tolerance pivoting etc. This package does not allow pivoting 122 // based on sparsity criterion. 123 124 // The factors L and U are computed via integer preserving operations via 125 // integer-preserving Gaussian elimination. The key part of this algorithm 126 // is a Roundoff Error Free (REF) sparse triangular solve function which 127 // exploits sparsity to reduce the number of operations that must be 128 // performed. 129 130 // Once L and U are computed, a simplified version of the triangular solve 131 // is performed which assumes the vector b is dense. The final solution 132 // vector x is gauranteed to be exact. This vector can be output in one of 133 // three ways: 1) full precision rational arithmetic (as a sequence of 134 // numerators and denominators) using the GMP mpq_t data type, 2) double 135 // precision while not exact will produce a solution accurate to machine 136 // roundoff unless the size of the associated solution exceeds double 137 // precision (i.e., the solution is 10^500 or something), 3) variable 138 // precision floating point using the GMP mpfr_t data type. The associated 139 // precision is user defined. 140 141 142 //------------------------------------------------------------------------------ 143 //------------------------------------------------------------------------------ 144 //---------------------Include files required by SLIP LU------------------------ 145 //------------------------------------------------------------------------------ 146 //------------------------------------------------------------------------------ 147 148 #include <stdio.h> 149 #include <stdbool.h> 150 #include <stdint.h> 151 #include <stdlib.h> 152 #include <string.h> 153 #include <gmp.h> 154 #include <mpfr.h> 155 #include "SuiteSparse_config.h" 156 157 //------------------------------------------------------------------------------ 158 // Version 159 //------------------------------------------------------------------------------ 160 161 // Current version of the code 162 #define SLIP_LU_VERSION "1.0.2" 163 #define SLIP_LU_VERSION_MAJOR 1 164 #define SLIP_LU_VERSION_MINOR 0 165 #define SLIP_LU_VERSION_SUB 2 166 167 //------------------------------------------------------------------------------ 168 // Error codes 169 //------------------------------------------------------------------------------ 170 171 // Most SLIP_LU functions return a code that indicates if it was successful 172 // or not. Otherwise the code returns a pointer to the object that was created 173 // or it returns void (in the case that an object was deleted) 174 175 typedef enum 176 { 177 SLIP_OK = 0, // all is well 178 SLIP_OUT_OF_MEMORY = -1, // out of memory 179 SLIP_SINGULAR = -2, // the input matrix A is singular 180 SLIP_INCORRECT_INPUT = -3, // one or more input arguments are incorrect 181 SLIP_INCORRECT = -4, // The solution is incorrect 182 SLIP_PANIC = -5 // SLIP_LU used without proper initialization 183 } 184 SLIP_info ; 185 186 //------------------------------------------------------------------------------ 187 // Pivot scheme codes 188 //------------------------------------------------------------------------------ 189 190 // A code in SLIP_options to tell SLIP LU what type of pivoting to use. 191 192 typedef enum 193 { 194 SLIP_SMALLEST = 0, // Smallest pivot 195 SLIP_DIAGONAL = 1, // Diagonal pivoting 196 SLIP_FIRST_NONZERO = 2, // First nonzero per column chosen as pivot 197 SLIP_TOL_SMALLEST = 3, // Diagonal pivoting with tol for smallest pivot. 198 // (Default) 199 SLIP_TOL_LARGEST = 4, // Diagonal pivoting with tol. for largest pivot 200 SLIP_LARGEST = 5 // Largest pivot 201 } 202 SLIP_pivot ; 203 204 //------------------------------------------------------------------------------ 205 // Column ordering scheme codes 206 //------------------------------------------------------------------------------ 207 208 // A code in SLIP_options to tell SLIP LU what column ordering to use. 209 210 typedef enum 211 { 212 SLIP_NO_ORDERING = 0, // None: A is factorized as-is 213 SLIP_COLAMD = 1, // COLAMD: Default 214 SLIP_AMD = 2 // AMD 215 } 216 SLIP_col_order ; 217 218 //------------------------------------------------------------------------------ 219 //------------------------------------------------------------------------------ 220 //-------------------------Data Structures-------------------------------------- 221 //------------------------------------------------------------------------------ 222 //------------------------------------------------------------------------------ 223 224 // This struct serves as a global struct to define all user-selectable options. 225 226 typedef struct SLIP_options 227 { 228 SLIP_pivot pivot ; // row pivoting scheme used. 229 SLIP_col_order order ; // column ordering scheme used 230 double tol ; // tolerance for the row-pivotin methods 231 // SLIP_TOL_SMALLEST and SLIP_TOL_LARGEST 232 int print_level ; // 0: print nothing, 1: just errors, 233 // 2: terse (basic stats from COLAMD/AMD and 234 // SLIP LU), 3: all, with matrices and results 235 int32_t prec ; // Precision used to output file if MPFR is chosen 236 mpfr_rnd_t round ; // Type of MPFR rounding used 237 bool check ; // Set true if the solution to the system should be 238 // checked. Intended for debugging only; SLIP_LU is 239 // guaranteed to return the exact solution. 240 } SLIP_options ; 241 242 // Purpose: Create and return SLIP_options object with default parameters 243 // upon successful allocation, which are defined in SLIP_LU_internal.h 244 // To free it, simply use SLIP_FREE (option). 245 SLIP_options* SLIP_create_default_options (void) ; 246 247 //------------------------------------------------------------------------------ 248 // SLIP_matrix: a sparse CSC, sparse triplet, or dense matrix 249 //------------------------------------------------------------------------------ 250 251 // SLIP LU uses a single matrix data type, SLIP_matrix, which can be held in 252 // one of three kinds of formats: sparse CSC (compressed sparse column), 253 // sparse triplet, and dense: 254 255 typedef enum 256 { 257 SLIP_CSC = 0, // matrix is in compressed sparse column format 258 SLIP_TRIPLET = 1, // matrix is in sparse triplet format 259 SLIP_DENSE = 2 // matrix is in dense format 260 } 261 SLIP_kind ; 262 263 // Each of the three formats can have values of 5 different data types: mpz_t, 264 // mpq_t, mpfr_t, int64_t, and double: 265 266 typedef enum 267 { 268 SLIP_MPZ = 0, // matrix of mpz_t integers 269 SLIP_MPQ = 1, // matrix of mpq_t rational numbers 270 SLIP_MPFR = 2, // matrix of mpfr_t 271 SLIP_INT64 = 3, // matrix of int64_t integers 272 SLIP_FP64 = 4 // matrix of doubles 273 } 274 SLIP_type ; 275 276 // This gives a total of 15 different matrix types. Not all functions accept 277 // all 15 matrices types, however. 278 279 // Suppose A is an m-by-n matrix with nz <= nzmax entries. 280 // The p, i, j, and x components are defined as: 281 282 // (0) SLIP_CSC: A sparse matrix in CSC (compressed sparse column) format. 283 // A->p is an int64_t array of size n+1, A->i is an int64_t array of size 284 // nzmax (with nz <= nzmax), and A->x.type is an array of size nzmax of 285 // matrix entries ('type' is one of mpz, mpq, mpfr, int64, or fp64). The 286 // row indices of column j appear in A->i [A->p [j] ... A->p [j+1]-1], and 287 // the values appear in the same locations in A->x.type. The A->j array 288 // is NULL. A->nz is ignored; nz is A->p [A->n]. 289 290 // (1) SLIP_TRIPLET: A sparse matrix in triplet format. A->i and A->j are 291 // both int64_t arrays of size nzmax, and A->x.type is an array of values 292 // of the same size. The kth tuple has row index A->i [k], column index 293 // A->j [k], and value A->x.type [k], with 0 <= k < A->nz. The A->p array 294 // is NULL. 295 296 // (2) SLIP_DENSE: A dense matrix. The integer arrays A->p, A->i, and A->j 297 // are all NULL. A->x.type is a pointer to an array of size m*n, stored 298 // in column-oriented format. The value of A(i,j) is A->x.type [p] 299 // with p = i + j*A->m. A->nz is ignored; nz is A->m * A->n. 300 301 // The SLIP_matrix may contain 'shallow' components, A->p, A->i, A->j, and 302 // A->x. For example, if A->p_shallow is true, then a non-NULL A->p is a 303 // pointer to a read-only array, and the A->p array is not freed by 304 // SLIP_matrix_free. If A->p is NULL (for a triplet or dense matrix), then 305 // A->p_shallow has no effect. 306 307 typedef struct 308 { 309 int64_t m ; // number of rows 310 int64_t n ; // number of columns 311 int64_t nzmax ; // size of A->i, A->j, and A->x 312 int64_t nz ; // # nonzeros in a triplet matrix . 313 // Ignored for CSC and dense matrices. 314 SLIP_kind kind ; // CSC, triplet, or dense 315 SLIP_type type ; // mpz, mpq, mpfr, int64, or fp64 (double) 316 317 int64_t *p ; // if CSC: column pointers, an array size is n+1. 318 // if triplet or dense: A->p is NULL. 319 bool p_shallow ; // if true, A->p is shallow. 320 321 int64_t *i ; // if CSC or triplet: row indices, of size nzmax. 322 // if dense: A->i is NULL. 323 bool i_shallow ; // if true, A->i is shallow. 324 325 int64_t *j ; // if triplet: column indices, of size nzmax. 326 // if CSC or dense: A->j is NULL. 327 bool j_shallow ; // if true, A->j is shallow. 328 329 union // A->x.type has size nzmax. 330 { 331 mpz_t *mpz ; // A->x.mpz 332 mpq_t *mpq ; // A->x.mpq 333 mpfr_t *mpfr ; // A->x.mpfr 334 int64_t *int64 ; // A->x.int64 335 double *fp64 ; // A->x.fp64 336 } x ; 337 bool x_shallow ; // if true, A->x.type is shallow. 338 339 mpq_t scale ; // scale factor for mpz matrices (never shallow) 340 // For all matrices who's type is not mpz, 341 // mpz_scale = 1. 342 343 } SLIP_matrix ; 344 345 //------------------------------------------------------------------------------ 346 // SLIP_matrix_allocate: allocate an m-by-n SLIP_matrix 347 //------------------------------------------------------------------------------ 348 349 // if shallow is false: All components (p,i,j,x) are allocated and set to zero, 350 // and then shallow flags are all false. 351 352 // if shallow is true: All components (p,i,j,x) are NULL, and their shallow 353 // flags are all true. The user can then set A->p, 354 // A->i, A->j, and/or A->x accordingly, from their own 355 // arrays. 356 357 SLIP_info SLIP_matrix_allocate 358 ( 359 SLIP_matrix **A_handle, // matrix to allocate 360 SLIP_kind kind, // CSC, triplet, or dense 361 SLIP_type type, // mpz, mpq, mpfr, int64, or double 362 int64_t m, // # of rows 363 int64_t n, // # of columns 364 int64_t nzmax, // max # of entries 365 bool shallow, // if true, matrix is shallow. A->p, A->i, A->j, 366 // A->x are all returned as NULL and must be set 367 // by the caller. All A->*_shallow are returned 368 // as true. 369 bool init, // If true, and the data types are mpz, mpq, or 370 // mpfr, the entries are initialized (using the 371 // appropriate SLIP_mp*_init function). If false, 372 // the mpz, mpq, and mpfr arrays are allocated but 373 // not initialized. 374 const SLIP_options *option 375 ) ; 376 377 //------------------------------------------------------------------------------ 378 // SLIP_matrix_free: free a SLIP_matrix 379 //------------------------------------------------------------------------------ 380 381 SLIP_info SLIP_matrix_free 382 ( 383 SLIP_matrix **A_handle, // matrix to free 384 const SLIP_options *option 385 ) ; 386 387 //------------------------------------------------------------------------------ 388 // SLIP_matrix_nnz: # of entries in a matrix 389 //------------------------------------------------------------------------------ 390 391 int64_t SLIP_matrix_nnz // return # of entries in A, or -1 on error 392 ( 393 const SLIP_matrix *A, // matrix to query 394 const SLIP_options *option 395 ) ; 396 397 //------------------------------------------------------------------------------ 398 // SLIP_matrix_copy: makes a copy of a matrix 399 //------------------------------------------------------------------------------ 400 401 // SLIP_matrix_copy: make a copy of a SLIP_matrix, into another kind and type. 402 403 SLIP_info SLIP_matrix_copy 404 ( 405 SLIP_matrix **C_handle, // matrix to create (never shallow) 406 // inputs, not modified: 407 SLIP_kind C_kind, // C->kind: CSC, triplet, or dense 408 SLIP_type C_type, // C->type: mpz_t, mpq_t, mpfr_t, int64_t, or double 409 SLIP_matrix *A, // matrix to make a copy of (may be shallow) 410 const SLIP_options *option 411 ) ; 412 413 //------------------------------------------------------------------------------ 414 // SLIP_matrix macros 415 //------------------------------------------------------------------------------ 416 417 // These macros simplify the access to entries in a SLIP_matrix. 418 // The type parameter is one of: mpq, mpz, mpfr, int64, or fp64. 419 420 // To access the kth entry in a SLIP_matrix using 1D linear addressing, 421 // in any matrix kind (CSC, triplet, or dense), in any type: 422 #define SLIP_1D(A,k,type) ((A)->x.type [k]) 423 424 // To access the (i,j)th entry in a 2D SLIP_matrix, in any type: 425 #define SLIP_2D(A,i,j,type) SLIP_1D (A, (i)+(j)*((A)->m), type) 426 427 //------------------------------------------------------------------------------ 428 // SLIP_LU_analysis: symbolic pre-analysis 429 //------------------------------------------------------------------------------ 430 431 // This struct stores the column permutation for LU and the estimate of the 432 // number of nonzeros in L and U. 433 434 typedef struct 435 { 436 int64_t *q ; // Column permutation for LU factorization, representing 437 // the permutation matrix Q. The matrix A*Q is factorized. 438 // If the kth column of L, U, and A*Q is column j of the 439 // unpermuted matrix A, then j = S->q [k]. 440 int64_t lnz ; // Approximate number of nonzeros in L. 441 int64_t unz ; // Approximate number of nonzeros in U. 442 // lnz and unz are used to allocate the initial space for 443 // L and U; the space is reallocated as needed. 444 } SLIP_LU_analysis ; 445 446 // The symbolic analysis object is created by SLIP_LU_analyze. 447 448 // SLIP_LU_analysis_free frees the SLIP_LU_analysis object. 449 SLIP_info SLIP_LU_analysis_free 450 ( 451 SLIP_LU_analysis **S, // Structure to be deleted 452 const SLIP_options *option 453 ) ; 454 455 //------------------------------------------------------------------------------ 456 // Memory management 457 //------------------------------------------------------------------------------ 458 459 // SLIP_LU relies on the SuiteSparse memory management functions, 460 // SuiteSparse_malloc, SuiteSparse_calloc, SuiteSparse_realloc, and 461 // SuiteSparse_free. 462 463 // Allocate and initialize memory space for SLIP_LU. 464 void *SLIP_calloc 465 ( 466 size_t nitems, // number of items to allocate 467 size_t size // size of each item 468 ) ; 469 470 // Allocate memory space for SLIP_LU. 471 void *SLIP_malloc 472 ( 473 size_t size // size of memory space to allocate 474 ) ; 475 476 // Free the memory allocated by SLIP_calloc, SLIP_malloc, or SLIP_realloc. 477 void SLIP_free 478 ( 479 void *p // pointer to memory space to free 480 ) ; 481 482 // Free a pointer and set it to NULL. 483 #define SLIP_FREE(p) \ 484 { \ 485 SLIP_free (p) ; \ 486 (p) = NULL ; \ 487 } 488 489 // SLIP_realloc is a wrapper for realloc. If p is non-NULL on input, it points 490 // to a previously allocated object of size old_size * size_of_item. The 491 // object is reallocated to be of size new_size * size_of_item. If p is NULL 492 // on input, then a new object of that size is allocated. On success, a 493 // pointer to the new object is returned. If the reallocation fails, p is not 494 // modified, and a flag is returned to indicate that the reallocation failed. 495 // If the size decreases or remains the same, then the method always succeeds 496 // (ok is returned as true). 497 498 // Typical usage: the following code fragment allocates an array of 10 int's, 499 // and then increases the size of the array to 20 int's. If the SLIP_malloc 500 // succeeds but the SLIP_realloc fails, then the array remains unmodified, 501 // of size 10. 502 // 503 // int *p ; 504 // p = SLIP_malloc (10 * sizeof (int)) ; 505 // if (p == NULL) { error here ... } 506 // printf ("p points to an array of size 10 * sizeof (int)\n") ; 507 // bool ok ; 508 // p = SLIP_realloc (20, 10, sizeof (int), p, &ok) ; 509 // if (ok) printf ("p has size 20 * sizeof (int)\n") ; 510 // else printf ("realloc failed; p still has size 10 * sizeof (int)\n") ; 511 // SLIP_free (p) ; 512 513 void *SLIP_realloc // pointer to reallocated block, or original block 514 // if the realloc failed 515 ( 516 int64_t nitems_new, // new number of items in the object 517 int64_t nitems_old, // old number of items in the object 518 size_t size_of_item, // sizeof each item 519 void *p, // old object to reallocate 520 bool *ok // true if success, false on failure 521 ) ; 522 523 //------------------------------------------------------------------------------ 524 // SLIP LU memory environment routines 525 //------------------------------------------------------------------------------ 526 527 // SLIP_initialize: initializes the working evironment for SLIP LU library. 528 // It must be called prior to calling any other SLIP_* function. 529 SLIP_info SLIP_initialize (void) ; 530 531 // SLIP_initialize_expert is the same as SLIP_initialize, except that it allows 532 // for a redefinition of custom memory functions that are used for SLIP_LU and 533 // GMP. The four inputs to this function are pointers to four functions with 534 // the same signatures as the ANSI C malloc, calloc, realloc, and free. 535 SLIP_info SLIP_initialize_expert 536 ( 537 void* (*MyMalloc) (size_t), // user-defined malloc 538 void* (*MyCalloc) (size_t, size_t), // user-defined calloc 539 void* (*MyRealloc) (void *, size_t), // user-defined realloc 540 void (*MyFree) (void *) // user-defined free 541 ) ; 542 543 // SLIP_finalize: This function finalizes the working evironment for SLIP LU 544 // library, and frees any internal workspace created by SLIP_LU. It must be 545 // called as the last SLIP_* function called. 546 SLIP_info SLIP_finalize (void) ; 547 548 //------------------------------------------------------------------------------ 549 // Primary factorization & solve routines 550 //------------------------------------------------------------------------------ 551 552 // SLIP_backslash solves the linear system Ax = b. This is the simplest way to 553 // use the SLIP LU package. This function encompasses both factorization and 554 // solve and returns the solution vector in the user desired type. It can be 555 // thought of as an exact version of MATLAB sparse backslash. 556 SLIP_info SLIP_backslash 557 ( 558 // Output 559 SLIP_matrix **X_handle, // Final solution vector 560 // Input 561 SLIP_type type, // Type of output desired: 562 // Must be SLIP_MPQ, SLIP_MPFR, 563 // or SLIP_FP64 564 const SLIP_matrix *A, // Input matrix 565 const SLIP_matrix *b, // Right hand side vector(s) 566 const SLIP_options* option 567 ) ; 568 569 // SLIP_LU_analyze performs the symbolic ordering and analysis for SLIP LU. 570 // Currently, there are three options: no ordering, COLAMD, and AMD. 571 SLIP_info SLIP_LU_analyze 572 ( 573 SLIP_LU_analysis **S, // symbolic analysis (column permutation and nnz L,U) 574 const SLIP_matrix *A, // Input matrix 575 const SLIP_options *option // Control parameters 576 ) ; 577 578 // SLIP_LU_factorize performs the SLIP LU factorization. This factorization is 579 // done via n iterations of the sparse REF triangular solve function. The 580 // overall factorization is PAQ = LDU. The determinant can be obtained as 581 // rhos->x.mpz[n-1]. 582 // 583 // L: undefined on input, created on output 584 // U: undefined on input, created on output 585 // rhos: undefined on input, created on output 586 // pinv: undefined on input, created on output 587 // 588 // A: input only, not modified 589 // S: input only, not modified 590 // option: input only, not modified 591 SLIP_info SLIP_LU_factorize 592 ( 593 // output: 594 SLIP_matrix **L_handle, // lower triangular matrix 595 SLIP_matrix **U_handle, // upper triangular matrix 596 SLIP_matrix **rhos_handle, // sequence of pivots 597 int64_t **pinv_handle, // inverse row permutation 598 // input: 599 const SLIP_matrix *A, // matrix to be factored 600 const SLIP_LU_analysis *S, // column permutation and estimates 601 // of nnz in L and U 602 const SLIP_options* option 603 ) ; 604 605 // SLIP_LU_solve solves the linear system LD^(-1)U x = b. 606 SLIP_info SLIP_LU_solve // solves the linear system LD^(-1)U x = b 607 ( 608 // Output 609 SLIP_matrix **X_handle, // rational solution to the system 610 // input: 611 const SLIP_matrix *b, // right hand side vector 612 const SLIP_matrix *A, // Input matrix 613 const SLIP_matrix *L, // lower triangular matrix 614 const SLIP_matrix *U, // upper triangular matrix 615 const SLIP_matrix *rhos, // sequence of pivots 616 const SLIP_LU_analysis *S, // symbolic analysis struct 617 const int64_t *pinv, // inverse row permutation 618 const SLIP_options* option 619 ) ; 620 621 // SLIP_matrix_check: check and print a SLIP_sparse matrix 622 SLIP_info SLIP_matrix_check // returns a SLIP_LU status code 623 ( 624 const SLIP_matrix *A, // matrix to check 625 const SLIP_options* option // defines the print level 626 ) ; 627 628 //------------------------------------------------------------------------------ 629 //---------------------------SLIP GMP/MPFR Functions---------------------------- 630 //------------------------------------------------------------------------------ 631 632 // The following functions are the SLIP LU interface to the GMP/MPFR libary. 633 // Each corresponding GMP/MPFR function is given a wrapper to ensure that no 634 // memory leaks or crashes occur. All covered GMP functions can be found in 635 // SLIP_gmp.c 636 637 // The GMP library does not handle out-of-memory failures. However, it does 638 // provide a mechanism for passing function pointers that replace GMP's use of 639 // malloc, realloc, and free. This mechanism is used to provide a try/catch 640 // mechanism for memory allocation errors, using setjmp and longjmp. 641 642 // When a GMP function is called, this wrapper keeps track of a list of objects 643 // allocated by that function. The list is started fresh each time a GMP 644 // function is called. If any allocation fails, the NULL pointer is not 645 // returned to GMP. Instead, all allocated blocks in the list are freed, 646 // and slip_gmp_allocate returns directly to wrapper. 647 648 SLIP_info SLIP_mpfr_asprintf (char **str, const char *format, ... ) ; 649 650 SLIP_info SLIP_gmp_fscanf (FILE *fp, const char *format, ... ) ; 651 652 SLIP_info SLIP_mpz_init (mpz_t x) ; 653 654 SLIP_info SLIP_mpz_init2(mpz_t x, const size_t size) ; 655 656 SLIP_info SLIP_mpz_set (mpz_t x, const mpz_t y) ; 657 658 SLIP_info SLIP_mpz_set_ui (mpz_t x, const uint64_t y) ; 659 660 SLIP_info SLIP_mpz_set_si (mpz_t x, const int64_t y) ; 661 662 SLIP_info SLIP_mpz_get_d (double *x, const mpz_t y) ; 663 664 SLIP_info SLIP_mpz_get_si (int64_t *x, const mpz_t y) ; 665 666 SLIP_info SLIP_mpz_set_q (mpz_t x, const mpq_t y) ; 667 668 SLIP_info SLIP_mpz_mul (mpz_t a, const mpz_t b, const mpz_t c) ; 669 670 SLIP_info SLIP_mpz_submul (mpz_t x, const mpz_t y, const mpz_t z) ; 671 672 SLIP_info SLIP_mpz_divexact (mpz_t x, const mpz_t y, const mpz_t z) ; 673 674 SLIP_info SLIP_mpz_gcd (mpz_t x, const mpz_t y, const mpz_t z) ; 675 676 SLIP_info SLIP_mpz_lcm (mpz_t lcm, const mpz_t x, const mpz_t y) ; 677 678 SLIP_info SLIP_mpz_abs (mpz_t x, const mpz_t y) ; 679 680 SLIP_info SLIP_mpz_cmp (int *r, const mpz_t x, const mpz_t y) ; 681 682 SLIP_info SLIP_mpz_cmpabs (int *r, const mpz_t x, const mpz_t y) ; 683 684 SLIP_info SLIP_mpz_cmp_ui (int *r, const mpz_t x, const uint64_t y) ; 685 686 SLIP_info SLIP_mpz_sgn (int *sgn, const mpz_t x) ; 687 688 SLIP_info SLIP_mpz_sizeinbase (size_t *size, const mpz_t x, int64_t base) ; 689 690 SLIP_info SLIP_mpq_init (mpq_t x) ; 691 692 SLIP_info SLIP_mpq_set (mpq_t x, const mpq_t y) ; 693 694 SLIP_info SLIP_mpq_set_z (mpq_t x, const mpz_t y) ; 695 696 SLIP_info SLIP_mpq_set_d (mpq_t x, const double y) ; 697 698 SLIP_info SLIP_mpq_set_ui (mpq_t x, const uint64_t y, const uint64_t z) ; 699 700 SLIP_info SLIP_mpq_set_si (mpq_t x, const int64_t y, const uint64_t z) ; 701 702 SLIP_info SLIP_mpq_set_num (mpq_t x, const mpz_t y) ; 703 704 SLIP_info SLIP_mpq_set_den (mpq_t x, const mpz_t y) ; 705 706 SLIP_info SLIP_mpq_get_den (mpz_t x, const mpq_t y) ; 707 708 SLIP_info SLIP_mpq_get_d (double *x, const mpq_t y) ; 709 710 SLIP_info SLIP_mpq_abs (mpq_t x, const mpq_t y) ; 711 712 SLIP_info SLIP_mpq_add (mpq_t x, const mpq_t y, const mpq_t z) ; 713 714 SLIP_info SLIP_mpq_mul (mpq_t x, const mpq_t y, const mpq_t z) ; 715 716 SLIP_info SLIP_mpq_div (mpq_t x, const mpq_t y, const mpq_t z) ; 717 718 SLIP_info SLIP_mpq_cmp (int *r, const mpq_t x, const mpq_t y) ; 719 720 SLIP_info SLIP_mpq_cmp_ui (int *r, const mpq_t x, 721 const uint64_t num, const uint64_t den) ; 722 723 SLIP_info SLIP_mpq_sgn (int *sgn, const mpq_t x) ; 724 725 SLIP_info SLIP_mpq_equal (int *r, const mpq_t x, const mpq_t y) ; 726 727 SLIP_info SLIP_mpfr_init2(mpfr_t x, const uint64_t size) ; 728 729 SLIP_info SLIP_mpfr_set (mpfr_t x, const mpfr_t y, const mpfr_rnd_t rnd) ; 730 731 SLIP_info SLIP_mpfr_set_d (mpfr_t x, const double y, const mpfr_rnd_t rnd) ; 732 733 SLIP_info SLIP_mpfr_set_si (mpfr_t x, int64_t y, const mpfr_rnd_t rnd) ; 734 735 SLIP_info SLIP_mpfr_set_q (mpfr_t x, const mpq_t y, const mpfr_rnd_t rnd) ; 736 737 SLIP_info SLIP_mpfr_set_z (mpfr_t x, const mpz_t y, const mpfr_rnd_t rnd) ; 738 739 SLIP_info SLIP_mpfr_get_z (mpz_t x, const mpfr_t y, const mpfr_rnd_t rnd) ; 740 741 SLIP_info SLIP_mpfr_get_q (mpq_t x, const mpfr_t y, const mpfr_rnd_t rnd) ; 742 743 SLIP_info SLIP_mpfr_get_d (double *x, const mpfr_t y, const mpfr_rnd_t rnd) ; 744 745 SLIP_info SLIP_mpfr_get_si (int64_t *x, const mpfr_t y, const mpfr_rnd_t rnd) ; 746 747 SLIP_info SLIP_mpfr_mul (mpfr_t x, const mpfr_t y, const mpfr_t z, 748 const mpfr_rnd_t rnd) ; 749 750 SLIP_info SLIP_mpfr_mul_d (mpfr_t x, const mpfr_t y, const double z, 751 const mpfr_rnd_t rnd) ; 752 753 SLIP_info SLIP_mpfr_div_d (mpfr_t x, const mpfr_t y, const double z, 754 const mpfr_rnd_t rnd) ; 755 756 SLIP_info SLIP_mpfr_ui_pow_ui (mpfr_t x, const uint64_t y, const uint64_t z, 757 const mpfr_rnd_t rnd) ; 758 759 SLIP_info SLIP_mpfr_sgn (int *sgn, const mpfr_t x) ; 760 761 SLIP_info SLIP_mpfr_free_cache (void) ; 762 763 SLIP_info SLIP_mpfr_free_str (char *str) ; 764 765 #if 0 766 // These functions are currently unused, but kept here for future reference. 767 SLIP_info SLIP_gmp_asprintf (char **str, const char *format, ... ) ; 768 SLIP_info SLIP_gmp_printf (const char *format, ... ) ; 769 SLIP_info SLIP_mpfr_printf ( const char *format, ... ) ; 770 SLIP_info SLIP_gmp_fprintf (FILE *fp, const char *format, ... ) ; 771 SLIP_info SLIP_mpfr_fprintf (FILE *fp, const char *format, ... ) ; 772 SLIP_info SLIP_mpz_set_d (mpz_t x, const double y) ; 773 SLIP_info SLIP_mpz_add (mpz_t a, const mpz_t b, const mpz_t c) ; 774 SLIP_info SLIP_mpz_addmul (mpz_t x, const mpz_t y, const mpz_t z) ; 775 SLIP_info SLIP_mpfr_log2(mpfr_t x, const mpfr_t y, const mpfr_rnd_t rnd) ; 776 #endif 777 778 #endif 779 780