1 #ifndef _LU_INTERNAL_H 2 #define _LU_INTERNAL_H 3 4 #include "lu_def.h" 5 6 /* -------------------------------------------------------------------------- */ 7 /* struct lu */ 8 /* -------------------------------------------------------------------------- */ 9 10 /* 11 This data structure provides access to istore, xstore. 12 13 lu_* routines do not access istore, xstore directly. Instead, they operate 14 on a struct lu object. Scalar quantities stored in istore, xstore are copied 15 to a struct lu object by lu_load() and copied back by lu_save(). Subarrays 16 of istore, xstore and the user arrays Li, Lx, Ui, Ux, Wi, Wx are aliased by 17 pointers in struct lu. 18 */ 19 20 struct lu 21 { 22 /* user parameters, not modified */ 23 lu_int Lmem; 24 lu_int Umem; 25 lu_int Wmem; 26 double droptol; 27 double abstol; 28 double reltol; 29 lu_int nzbias; 30 lu_int maxsearch; 31 lu_int pad; 32 double stretch; 33 double compress_thres; 34 double sparse_thres; 35 lu_int search_rows; 36 37 /* user readable */ 38 lu_int m; 39 lu_int addmemL; 40 lu_int addmemU; 41 lu_int addmemW; 42 43 lu_int nupdate; 44 lu_int nforrest; 45 lu_int nfactorize; 46 lu_int nupdate_total; 47 lu_int nforrest_total; 48 lu_int nsymperm_total; 49 lu_int Lnz; /* nz in L excluding diagonal */ 50 lu_int Unz; /* nz in U excluding diagonal */ 51 lu_int Rnz; /* nz in update etas excluding diagonal */ 52 double min_pivot; 53 double max_pivot; 54 double max_eta; 55 double update_cost_numer; 56 double update_cost_denom; 57 double time_factorize; 58 double time_solve; 59 double time_update; 60 double time_factorize_total; 61 double time_solve_total; 62 double time_update_total; 63 lu_int Lflops; 64 lu_int Uflops; 65 lu_int Rflops; 66 double condestL; 67 double condestU; 68 double normL; 69 double normU; 70 double normestLinv; 71 double normestUinv; 72 double onenorm; /* 1-norm and inf-norm of matrix after fresh */ 73 double infnorm; /* factorization with dependent cols replaced */ 74 double residual_test; /* computed by lu_residual_test() */ 75 76 lu_int matrix_nz; /* nz in basis matrix when factorized */ 77 lu_int rank; /* rank of basis matrix when factorized */ 78 lu_int bump_size; 79 lu_int bump_nz; 80 lu_int nsearch_pivot; /* # rows/cols searched for pivot */ 81 lu_int nexpand; /* # rows/cols expanded in factorize */ 82 lu_int ngarbage; /* # garbage collections in factorize */ 83 lu_int factor_flops; /* # flops in factorize */ 84 double time_singletons; 85 double time_search_pivot; 86 double time_elim_pivot; 87 88 double pivot_error; /* error estimate for pivot in last update */ 89 90 /* private */ 91 lu_int task; /* the part of factorization in progress */ 92 lu_int pivot_row; /* chosen pivot row */ 93 lu_int pivot_col; /* chosen pivot column */ 94 lu_int ftran_for_update; /* >= 0 if FTRAN prepared for update */ 95 lu_int btran_for_update; /* >= 0 if BTRAN prepared for update */ 96 lu_int marker; /* see @marked, below */ 97 lu_int pivotlen; /* length of @pivotcol, @pivotrow; <= 2*m */ 98 lu_int rankdef; /* # columns removed from active submatrix 99 because maximum was 0 or < abstol */ 100 lu_int min_colnz; /* colcount lists 1..min_colnz-1 are empty */ 101 lu_int min_rownz; /* rowcount lists 1..min_rownz-1 are empty */ 102 103 /* aliases to user arrays */ 104 lu_int *Lindex, *Uindex, *Windex; 105 double *Lvalue, *Uvalue, *Wvalue; 106 107 /* 108 * pointers into istore 109 * 110 * When two declaration lists are on one line, then the arrays from the 111 * second list share memory with the array from the first list. The arrays 112 * from the first lists are used during factorization, the arrays from the 113 * second lists are used during solves/updates. 114 */ 115 116 /* documented in lu_singletons.c, lu_setup_bump.c, lu_build_factors.c */ 117 lu_int *colcount_flink; lu_int *pivotcol; 118 lu_int *colcount_blink; lu_int *pivotrow; 119 lu_int *rowcount_flink; lu_int *Rbegin, *eta_row; 120 lu_int *rowcount_blink; lu_int *iwork1; 121 lu_int *Wbegin; lu_int *Lbegin; /* + Wbegin reused */ 122 lu_int *Wend; lu_int *Ltbegin; /* + Wend reused */ 123 lu_int *Wflink; lu_int *Ltbegin_p; /* + Wflink reused */ 124 lu_int *Wblink; lu_int *p; /* + Wblink reused */ 125 lu_int *pinv; lu_int *pmap; 126 lu_int *qinv; lu_int *qmap; 127 lu_int *Lbegin_p; /* Lbegin_p reused */ 128 lu_int *Ubegin; /* Ubegin reused */ 129 130 lu_int *iwork0; lu_int *marked; 131 /* iwork0: size m workspace, zeroed */ 132 /* marked: size m workspace, 0 <= marked[i] <= @marker */ 133 134 /* pointers into xstore */ 135 double *work0; /* size m workspace, zeroed */ 136 double *work1; /* size m workspace, uninitialized */ 137 double *col_pivot; /* pivot elements by column index */ 138 double *row_pivot; /* pivot elements by row index */ 139 }; 140 141 142 /* -------------------------------------------------------------------------- */ 143 /* Internal function prototypes */ 144 /* -------------------------------------------------------------------------- */ 145 146 lu_int lu_load( 147 struct lu *this, lu_int *istore, double *xstore, lu_int *Li, double *Lx, 148 lu_int *Ui, double *Ux, lu_int *Wi, double *Wx); 149 150 lu_int lu_save( 151 const struct lu *this, lu_int *istore, double *xstore, lu_int status); 152 153 void lu_reset(struct lu *this); 154 155 void lu_initialize(lu_int m, lu_int *istore, double *xstore); 156 157 lu_int lu_factorize_bump (struct lu *this); 158 159 lu_int lu_build_factors(struct lu *this); 160 161 void lu_garbage_perm(struct lu *this); 162 163 lu_int lu_markowitz(struct lu *this); 164 165 lu_int lu_pivot(struct lu *this); 166 167 lu_int lu_setup_bump( 168 struct lu *this, const lu_int *Bbegin, const lu_int *Bend, const lu_int *Bi, 169 const double *Bx); 170 171 lu_int lu_singletons( 172 struct lu *this, const lu_int *Bbegin, const lu_int *Bend, const lu_int *Bi, 173 const double *Bx); 174 175 void lu_solve_dense( 176 struct lu *this, const double *rhs, double *lhs, char trans); 177 178 lu_int lu_solve_for_update( 179 struct lu *this, const lu_int nrhs, const lu_int *irhs, const double *xrhs, 180 lu_int *nlhs, lu_int *ilhs, double *xlhs, char trans); 181 182 void lu_solve_sparse( 183 struct lu *this, const lu_int nrhs, const lu_int *irhs, const double *xrhs, 184 lu_int *nlhs, lu_int *ilhs, double *xlhs, char trans); 185 186 lu_int lu_dfs( 187 lu_int i, const lu_int *begin, const lu_int *end, const lu_int *index, 188 lu_int top, lu_int *xi, lu_int *pstack, lu_int *marked, const lu_int M); 189 190 lu_int lu_solve_symbolic( 191 const lu_int m, const lu_int *begin, const lu_int *end, const lu_int *index, 192 const lu_int nrhs, const lu_int *irhs, lu_int *ilhs, lu_int *pstack, 193 lu_int *marked, const lu_int M); 194 195 lu_int lu_solve_triangular( 196 const lu_int nz_symb, const lu_int *pattern_symb, const lu_int *begin, 197 const lu_int *end, const lu_int *index, const double *value, 198 const double *pivot, const double droptol, double *lhs, lu_int *pattern, 199 lu_int *flops); 200 201 lu_int lu_update(struct lu *this, double xtbl); 202 203 double lu_condest( 204 lu_int m, const lu_int *Ubegin, const lu_int *Ui, const double *Ux, 205 const double *pivot, const lu_int *perm, int upper, double *work, 206 double *norm, double *norminv); 207 208 double lu_normest( 209 lu_int m, const lu_int *Ubegin, const lu_int *Ui, const double *Ux, 210 const double *pivot, const lu_int *perm, int upper, double *work); 211 212 void lu_matrix_norm( 213 struct lu *this, const lu_int *Bbegin, const lu_int *Bend, const lu_int *Bi, 214 const double *Bx); 215 216 void lu_residual_test( 217 struct lu *this, const lu_int *Bbegin, const lu_int *Bend, const lu_int *Bi, 218 const double *Bx); 219 220 #endif 221