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