1 //------------------------------------------------------------------------------
2 // SLIP_LU/slip_get_pivot: find a pivot entry in a column
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 /* This function performs the pivoting for the SLIP LU factorization.
12 * The optional Order is:
13 *
14 * SLIP_SMALLEST = 0, Smallest pivot
15 * SLIP_DIAGONAL = 1, Diagonal pivoting
16 * SLIP_FIRST_NONZERO = 2, First nonzero per column chosen as pivot
17 * SLIP_TOL_SMALLEST = 3, Diagonal pivoting with tolerance for pivot. (Default)
18 * SLIP_TOL_LARGEST = 4, Diagonal pivoting with tolerance for largest pivot
19 * SLIP_LARGEST = 5 Largest pivot
20 *
21 * Options 2, 4 and 5 are not recommended and may lead to significant drops in
22 * performance.
23 *
24 * On output, the pivs, rhos, pinv, and row_perm arrays are all modified.
25 */
26
27 #define SLIP_FREE_ALL \
28 SLIP_MPQ_CLEAR (tol) ; \
29 SLIP_MPQ_CLEAR (ratio) ;
30
31 #include "slip_internal.h"
32
slip_get_pivot(int64_t * pivot,SLIP_matrix * x,int64_t * pivs,int64_t n,int64_t top,int64_t * xi,int64_t col,int64_t k,SLIP_matrix * rhos,int64_t * pinv,int64_t * row_perm,const SLIP_options * option)33 SLIP_info slip_get_pivot
34 (
35 int64_t *pivot, // found index of pivot entry
36 SLIP_matrix* x, // kth column of L and U
37 int64_t* pivs, // vector indicating which rows have been pivotal
38 int64_t n, // dimension of the problem
39 int64_t top, // nonzero pattern is located in xi[top..n-1]
40 int64_t* xi, // nonzero pattern of x
41 int64_t col, // current column of A (real kth column i.e., q[k])
42 int64_t k, // iteration of the algorithm
43 SLIP_matrix* rhos, // vector of pivots
44 int64_t* pinv, // row permutation
45 int64_t* row_perm, // opposite of pinv.
46 // if pinv[i] = j then row_perm[j] = i
47 const SLIP_options* option // command options
48 )
49 {
50
51 //--------------------------------------------------------------------------
52 // check inputs
53 //--------------------------------------------------------------------------
54
55 SLIP_info info ;
56 SLIP_REQUIRE(rhos, SLIP_DENSE, SLIP_MPZ);
57 SLIP_REQUIRE(x, SLIP_DENSE, SLIP_MPZ);
58
59 // inputs have been checked in the only caller SLIP_LU_factorize
60 // they are kept here for future reference
61 #if 0
62 if (!pivot || !pivs || !xi || !pinv || !row_perm )
63 {
64 return SLIP_INCORRECT_INPUT;
65 }
66 #endif
67 // pivoting method to use (see above description)
68 SLIP_pivot order = SLIP_OPTION_PIVOT(option);
69 // tolerance used if some tol-based pivoting is used
70 double tolerance = SLIP_OPTION_TOL(option);
71
72 //--------------------------------------------------------------------------
73 // allocate workspace
74 //--------------------------------------------------------------------------
75
76 int sgn, r;
77 mpq_t tol, ratio;
78 SLIP_MPQ_SET_NULL(tol);
79 SLIP_MPQ_SET_NULL(ratio);
80
81 //--------------------------------------------------------------------------
82 // Smallest pivot
83 //--------------------------------------------------------------------------
84
85 if (order == SLIP_SMALLEST)
86 {
87 SLIP_CHECK(slip_get_smallest_pivot(pivot, x, pivs, n, top, xi));
88 }
89
90 //--------------------------------------------------------------------------
91 // Diagonal
92 //--------------------------------------------------------------------------
93 else if (order == SLIP_DIAGONAL)
94 {
95 // Check if x[col] is eligible. take smallest pivot if not
96 SLIP_CHECK (SLIP_mpz_sgn(&sgn, x->x.mpz[col]));
97 if (sgn != 0 && pivs[col] < 0)
98 {
99 *pivot = col;
100 }
101 else
102 {
103 SLIP_CHECK (slip_get_smallest_pivot(pivot, x, pivs, n, top, xi));
104 }
105 }
106
107 //--------------------------------------------------------------------------
108 // First nonzero
109 //--------------------------------------------------------------------------
110 else if (order == SLIP_FIRST_NONZERO)
111 {
112 SLIP_CHECK (slip_get_nonzero_pivot(pivot, x, pivs, n, top, xi));
113 }
114
115 //--------------------------------------------------------------------------
116 // Tolerance with largest pivot
117 //--------------------------------------------------------------------------
118 else if (order == SLIP_TOL_LARGEST)
119 {
120 SLIP_CHECK (slip_get_largest_pivot(pivot, x, pivs, n, top, xi));
121
122 //----------------------------------------------------------------------
123 // Check x[col] vs largest potential pivot
124 //----------------------------------------------------------------------
125 SLIP_CHECK (SLIP_mpz_sgn(&sgn, x->x.mpz[col]));
126 if (sgn != 0 && pivs[col] < 0)
127 {
128 SLIP_CHECK(SLIP_mpq_init(tol));
129 SLIP_CHECK(SLIP_mpq_init(ratio));
130 // tol = user specified tolerance
131 SLIP_CHECK(SLIP_mpq_set_d(tol, tolerance));
132 // ratio = diagonal/largest
133 SLIP_CHECK(SLIP_mpq_set_num(ratio, x->x.mpz[col]));
134 SLIP_CHECK(SLIP_mpq_set_den(ratio, x->x.mpz[*pivot]));
135 // ratio = |ratio|
136 SLIP_CHECK(SLIP_mpq_abs(ratio, ratio));
137
138 // Is ratio >= tol?
139 SLIP_CHECK(SLIP_mpq_cmp(&r, ratio, tol));
140 if (r >= 0)
141 {
142 *pivot = col;
143 }
144 }
145 }
146
147 //--------------------------------------------------------------------------
148 // Use the largest potential pivot
149 //--------------------------------------------------------------------------
150 else if (order == SLIP_LARGEST)
151 {
152 SLIP_CHECK (slip_get_largest_pivot(pivot, x, pivs, n, top, xi));
153 }
154
155 //--------------------------------------------------------------------------
156 // Tolerance with smallest pivot (default option)
157 //--------------------------------------------------------------------------
158 else // if (order == SLIP_TOL_SMALLEST)
159 {
160 SLIP_CHECK (slip_get_smallest_pivot(pivot, x, pivs, n, top, xi)) ;
161
162 //----------------------------------------------------------------------
163 // Checking x[col] vs smallest pivot
164 //----------------------------------------------------------------------
165 SLIP_CHECK (SLIP_mpz_sgn(&sgn, x->x.mpz[col]));
166 if (sgn != 0 && pivs[col] < 0)
167 {
168
169 // Initialize tolerance and ratio
170 SLIP_CHECK(SLIP_mpq_init(tol));
171 SLIP_CHECK(SLIP_mpq_init(ratio));
172
173 // ratio = |smallest/diagonal|
174 SLIP_CHECK(SLIP_mpz_abs(SLIP_MPQ_NUM(ratio), x->x.mpz[*pivot]));
175 SLIP_CHECK(SLIP_mpz_abs(SLIP_MPQ_DEN(ratio), x->x.mpz[col]));
176
177 // Set user specified tolerance
178 SLIP_CHECK(SLIP_mpq_set_d(tol, tolerance));
179
180 // Is ratio >= tol?
181 SLIP_CHECK(SLIP_mpq_cmp(&r, ratio, tol));
182 if (r >= 0)
183 {
184 *pivot = col;
185 }
186 }
187 }
188
189 //--------------------------------------------------------------------------
190 // Reflect changes in row location & row_perm
191 //--------------------------------------------------------------------------
192 // Must move pivot into position k
193 int64_t intermed = pinv[*pivot];
194 int64_t intermed2 = row_perm[k];
195
196 //--------------------------------------------------------------------------
197 // Set row_perm[k] = pivot and row_perm[pinv[pivot]] = row_perm[k]
198 // Also, set pinv[pivot] = k and pinv[row_perm[k]] = pinv[pivot]
199 //--------------------------------------------------------------------------
200 row_perm[k] = *pivot;
201 row_perm[intermed] = intermed2;
202 pinv[*pivot] = k;
203 pinv[intermed2] = intermed;
204 // Row pivot is now pivotal
205 pivs[*pivot] = 1;
206
207 // Set the kth pivot.
208 size_t size;
209 // Get the size of x[pivot]
210 SLIP_CHECK(SLIP_mpz_sizeinbase(&size, x->x.mpz[*pivot], 2));
211 // GMP manual: Allocated size should be size+2
212 SLIP_CHECK(SLIP_mpz_init2(rhos->x.mpz[k], size+2));
213 // The kth pivot is x[pivot]
214 SLIP_CHECK (SLIP_mpz_set(rhos->x.mpz[k], x->x.mpz[*pivot]));
215
216 // Free memory
217 SLIP_FREE_ALL;
218 return SLIP_OK;
219 }
220
221