1 //------------------------------------------------------------------------------
2 // SLIP_LU/slip_forward_sub: sparse forward substitution (x = (LD)\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 roundoff-error-free (REF) forward
12 * substitution This is essentially the same as the sparse REF triangular solve
13 * applied to each column of the right hand side vectors. Like the normal one,
14 * this function expects that the matrix x is dense. As a result,the nonzero
15 * pattern is not computed and each nonzero in x is iterated across. The
16 * system to solve is L*D*x_output = x_input, overwriting the right-hand-side
17 * with the solution.
18 *
19 * On output, the SLIP_matrix* x structure is modified.
20 */
21
22 #define SLIP_FREE_ALL \
23 SLIP_matrix_free(&h, NULL) ;
24
25 #include "slip_internal.h"
26
slip_forward_sub(const SLIP_matrix * L,SLIP_matrix * x,const SLIP_matrix * rhos)27 SLIP_info slip_forward_sub
28 (
29 const SLIP_matrix *L, // lower triangular matrix
30 SLIP_matrix *x, // right hand side matrix of size n*numRHS
31 const SLIP_matrix *rhos // sequence of pivots used in factorization
32 )
33 {
34
35 //--------------------------------------------------------------------------
36 // check inputs
37 //--------------------------------------------------------------------------
38
39 SLIP_info info ;
40 SLIP_REQUIRE(L, SLIP_CSC, SLIP_MPZ);
41 SLIP_REQUIRE(x, SLIP_DENSE, SLIP_MPZ);
42 SLIP_REQUIRE(rhos, SLIP_DENSE, SLIP_MPZ);
43
44 //--------------------------------------------------------------------------
45
46 int64_t i, hx, k, j, jnew;
47 int sgn ;
48
49 // Build the history matrix
50 SLIP_matrix *h;
51 SLIP_CHECK (SLIP_matrix_allocate(&h, SLIP_DENSE, SLIP_INT64, x->m, x->n,
52 x->nzmax, false, true, NULL));
53
54 // initialize entries of history matrix to be -1
55 for (i = 0; i < x->nzmax; i++)
56 {
57 h->x.int64[i] = -1;
58 }
59
60
61 //--------------------------------------------------------------------------
62 // Iterate across each RHS vector
63 //--------------------------------------------------------------------------
64
65 for (k = 0; k < x->n; k++)
66 {
67
68 //----------------------------------------------------------------------
69 // Iterate accross all nonzeros in x. Assume x is dense
70 //----------------------------------------------------------------------
71
72 for (i = 0; i < x->m; i++)
73 {
74 hx = SLIP_2D(h, i, k, int64);
75 // If x[i][k] = 0, can skip operations and continue to next i
76 SLIP_CHECK(SLIP_mpz_sgn(&sgn, SLIP_2D(x, i, k, mpz)));
77 if (sgn == 0) {continue;}
78
79 //------------------------------------------------------------------
80 // History Update
81 //------------------------------------------------------------------
82
83 if (hx < i-1)
84 {
85 // x[i] = x[i] * rhos[i-1]
86 SLIP_CHECK(SLIP_mpz_mul( SLIP_2D(x, i, k, mpz),
87 SLIP_2D(x, i, k, mpz),
88 SLIP_1D(rhos, i-1, mpz)));
89 // x[i] = x[i] / rhos[hx]
90 if (hx > -1)
91 {
92 SLIP_CHECK(SLIP_mpz_divexact( SLIP_2D(x, i, k, mpz),
93 SLIP_2D(x, i, k, mpz),
94 SLIP_1D(rhos, hx, mpz)));
95 }
96 }
97
98 //------------------------------------------------------------------
99 // IPGE updates
100 //------------------------------------------------------------------
101
102 // Access the Lji
103 for (j = L->p[i]; j < L->p[i+1]; j++)
104 {
105 // Location of Lji
106 jnew = L->i[j];
107
108 // skip if Lx[j] is zero
109 SLIP_CHECK(SLIP_mpz_sgn(&sgn, L->x.mpz[j]));
110 if (sgn == 0) {continue;}
111
112 // j > i
113 if (jnew > i)
114 {
115 // check if x[jnew] is zero
116 SLIP_CHECK(SLIP_mpz_sgn(&sgn, SLIP_2D(x, jnew, k, mpz)));
117 if (sgn == 0)
118 {
119 // x[j] = x[j] - lji xi
120 SLIP_CHECK(SLIP_mpz_submul(SLIP_2D(x, jnew, k, mpz),
121 SLIP_1D(L, j, mpz),
122 SLIP_2D(x, i, k, mpz)));
123 // x[j] = x[j] / rhos[i-1]
124 if (i > 0)
125 {
126 SLIP_CHECK(
127 SLIP_mpz_divexact(SLIP_2D(x, jnew, k, mpz),
128 SLIP_2D(x, jnew, k, mpz),
129 SLIP_1D(rhos, i-1, mpz)));
130 }
131 }
132 else
133 {
134 hx = SLIP_2D(h, jnew, k, int64);
135 // History update if necessary
136 if (hx < i-1)
137 {
138 // x[j] = x[j] * rhos[i-1]
139 SLIP_CHECK(SLIP_mpz_mul(SLIP_2D(x, jnew, k, mpz),
140 SLIP_2D(x, jnew, k, mpz),
141 SLIP_1D(rhos, i-1, mpz)));
142 // x[j] = x[j] / rhos[hx]
143 if (hx > -1)
144 {
145 SLIP_CHECK(
146 SLIP_mpz_divexact(SLIP_2D(x, jnew, k, mpz),
147 SLIP_2D(x, jnew, k, mpz),
148 SLIP_1D(rhos, hx, mpz)));
149 }
150 }
151 // x[j] = x[j] * rhos[i]
152 SLIP_CHECK(SLIP_mpz_mul(SLIP_2D(x, jnew, k, mpz),
153 SLIP_2D(x, jnew, k, mpz),
154 SLIP_1D(rhos, i, mpz)));
155 // x[j] = x[j] - lmi xi
156 SLIP_CHECK(SLIP_mpz_submul(SLIP_2D(x, jnew, k, mpz),
157 SLIP_1D(L, j, mpz),
158 SLIP_2D(x, i, k, mpz)));
159 // x[j] = x[j] / rhos[i-1]
160 if (i > 0)
161 {
162 SLIP_CHECK(
163 SLIP_mpz_divexact(SLIP_2D(x, jnew, k, mpz),
164 SLIP_2D(x, jnew, k, mpz),
165 SLIP_1D(rhos, i-1, mpz)));
166 }
167 }
168 //h[jnew][k] = i;
169 SLIP_2D(h, jnew, k, int64) = i;
170 }
171 }
172 }
173 }
174
175 //--------------------------------------------------------------------------
176 // Free h memory
177 //--------------------------------------------------------------------------
178
179 SLIP_FREE_ALL;
180 return SLIP_OK;
181 }
182
183