1 //------------------------------------------------------------------------------
2 // SLIP_LU/slip_back_sub: sparse REF backward substitution (x = U\x)
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 /* Purpose: This function performs sparse REF backward substitution, solving
12 * the system Ux = b. Note that prior to this, x is multiplied by
13 * the determinant of A.
14 *
15 * U is a sparse mpz matrix, and bx is a dense mpz matrix. The diagonal entry
16 * of U must appear as the last entry in each column.
17 *
18 * The input argument bx contains b on input, and it is overwritten on output
19 * by the solution x.
20 */
21
22 #include "slip_internal.h"
23
slip_back_sub(const SLIP_matrix * U,SLIP_matrix * bx)24 SLIP_info slip_back_sub // performs sparse REF backward substitution
25 (
26 const SLIP_matrix *U, // input upper triangular matrix
27 SLIP_matrix *bx // right hand side matrix
28 )
29 {
30
31 //--------------------------------------------------------------------------
32 // check inputs
33 //--------------------------------------------------------------------------
34
35 SLIP_info info ;
36 SLIP_REQUIRE (U, SLIP_CSC, SLIP_MPZ) ;
37 SLIP_REQUIRE (bx, SLIP_DENSE, SLIP_MPZ) ;
38
39 //--------------------------------------------------------------------------
40
41 int sgn;
42 mpz_t *Ux = U->x.mpz;
43 int64_t *Ui = U->i;
44 int64_t *Up = U->p;
45
46 for (int64_t k = 0; k < bx->n; k++)
47 {
48 // Start at bx[n]
49 for (int64_t j = U->n-1; j >= 0; j--)
50 {
51 // If bx[j] is zero skip this iteration
52 SLIP_CHECK( SLIP_mpz_sgn( &sgn, SLIP_2D( bx, j, k, mpz)));
53 if (sgn == 0) {continue;}
54
55 // Obtain bx[j]
56 SLIP_CHECK(SLIP_mpz_divexact( SLIP_2D(bx, j, k, mpz),
57 SLIP_2D(bx, j, k, mpz),
58 Ux[Up[j+1]-1]));
59 for (int64_t i = Up[j]; i < Up[j+1]-1; i++)
60 {
61 SLIP_CHECK(SLIP_mpz_sgn(&sgn, Ux[i]));
62 if (sgn == 0) {continue;}
63 // bx[i] = bx[i] - Ux[i]*bx[j]
64 SLIP_CHECK(SLIP_mpz_submul( SLIP_2D(bx, Ui[i], k, mpz),
65 Ux[i], SLIP_2D(bx, j, k, mpz)));
66 }
67 }
68 }
69
70 return (SLIP_OK) ;
71 }
72
73