/[ascend]/branches/jacob/disused/solver/ascgauss.c
ViewVC logotype

Contents of /branches/jacob/disused/solver/ascgauss.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2933 - (show annotations) (download) (as text)
Mon Jun 1 20:38:42 2015 UTC (9 years, 1 month ago) by jacob
File MIME type: text/x-csrc
File size: 13645 byte(s)
Create a new branch from trunk r2932, with name 'jacob'
To be used by Jacob Shealy during GSOC2015 for development of ideal 
solution capabilities in FPROPS

1 /*
2 * ascgauss: sparse gauss elimination
3 * by Ken Tyner and Ben Allan
4 * Created: 4/14/97
5 * Version: $Revision: 1.6 $
6 * Version control file: $RCSfile: ascgauss.c,v $
7 * Date last modified: $Date: 2000/01/25 02:26:48 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the SLV solver.
11 * The SLV solver is free software; you can redistribute
12 * it and/or modify it under the terms of the GNU General Public License as
13 * published by the Free Software Foundation; either version 2 of the
14 * License, or (at your option) any later version.
15 *
16 * The SLV solver is distributed in hope that it will be
17 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with this program. If not, see <http://www.gnu.org/licenses/>.
23 */
24
25 /*
26 * This file contains gauss and modified gauss implementations.
27 * The first gauss is a partial column pivoted gauss L\U.
28 * Initially (4/97) it is not the fastest possible code, but is
29 * to be easily understood.
30 * It needs a header.
31 */
32 #define GAUSSALONE 0
33 #if GAUSSALONE
34 #include <math.h>
35 #include <ascend/general/platform.h>
36 #include <ascend/compiler/compiler.h>
37 #include <ascend/general/ascMalloc.h>
38 #include <ascend/general/mem.h>
39 #include <ascend/general/tm_time.h>
40 #include "mtx.h"
41 #include "linsolqr.h"
42 #define D_ZERO (real64)0.0
43 #define ZERO (int32)0
44 #define D_ONE (real64)1.0
45 #endif
46
47 #define GAUSS_DEBUG TRUE
48
49 /*
50 * See the comments attached to forward_substitute.
51 * This code otherwise is the same. It uses a 2 bodied
52 * matrix as well as making use of mtx_ALL_COLS and
53 * mtx_ALL_ROWS whereever possible.
54 * We make the following additional assumptions here:
55 * (if they do not hold, do NOT use this function)
56 * - sys->inverse (L) has no incidence that is not multipliers
57 * (and the diagonal of 1s is NOT in sys->inverse.)
58 * - sys->factors (U) has no incidence on the upper triangle,
59 * including the diagonal, or outside the factored region.
60 * relaxation: incidence anywhere allowed if value = 0.0
61 * since 0 doesn't contribute to a dot product
62 * and the only thing we do with triangles is dot them.
63 * - There may be singular rows and columns in the factorization,
64 * but any additions coming from these rows/columns during
65 * mtx_ALL_*O*S operations will not contribute to sums because the
66 * user zeroed the arr entries corresponding to these before
67 * calling this function.
68 */
69 static void forward_substitute_gauss2(linsolqr_system_t sys,
70 real64 *arr,
71 boolean transpose)
72 {
73 mtx_coord_t nz;
74 real64 sum, *pivlist;
75 mtx_matrix_t mtx;
76 int32 dotlim;
77 boolean nonzero_found=FALSE;
78
79 pivlist=sys->ludata->pivlist;
80 dotlim = sys->rng.low+sys->rank;
81 if (transpose) { /* arr is indexed by original column number */
82 mtx=sys->factors;
83 for( nz.col=sys->rng.low; nz.col < dotlim; ++(nz.col) ) {
84 register int32 org_col;
85
86 org_col = mtx_col_to_org(mtx,nz.col);
87 if (arr[org_col]!=D_ZERO) nonzero_found=TRUE;
88 if (nonzero_found) {
89 sum=mtx_col_dot_full_org_vec(mtx,nz.col,arr,mtx_ALL_ROWS,TRUE);
90 arr[org_col] = (arr[org_col] - sum) / pivlist[nz.col];
91 }
92 }
93 } else { /* arr is indexed by original row number */
94 mtx=sys->inverse;
95 for( nz.row=sys->rng.low; nz.row < dotlim; ++(nz.row) ) {
96 register int32 org_row;
97
98 org_row = mtx_row_to_org(mtx,nz.row);
99 if (arr[org_row]!=D_ZERO) nonzero_found=TRUE;
100 if (nonzero_found) {
101 sum = mtx_row_dot_full_org_vec(mtx,nz.row,arr,mtx_ALL_COLS,TRUE);
102 arr[org_row] -= sum;
103 }
104 }
105 }
106 }
107
108 /*
109 * See the comments attached to backward_substitute and
110 * forward_substitute2.
111 * When solving for the transpose, then we are actually
112 * running of the lower triangle, hence we use sys->factors.
113 * Otherwise we use sys->inverse which stores U.
114 */
115 static void backward_substitute_gauss2(linsolqr_system_t sys,
116 real64 *arr,
117 boolean transpose)
118 {
119 mtx_coord_t nz;
120 real64 sum, *pivlist;
121 mtx_matrix_t mtx;
122 int32 dotlim;
123 boolean nonzero_found=FALSE; /* once TRUE, substitution must be done
124 over remaining rows/cols */
125
126 dotlim=sys->rng.low;
127 pivlist=sys->ludata->pivlist;
128 if (transpose) { /* arr is indexed by original column number */
129 mtx = sys->inverse;
130 for( nz.col = sys->rng.low+sys->rank-1; nz.col >= dotlim ; --(nz.col) ) {
131 register int32 org_col;
132
133 org_col = mtx_col_to_org(mtx,nz.col);
134 if (arr[org_col] != D_ZERO) nonzero_found=TRUE;
135 if (nonzero_found) {
136 sum = mtx_col_dot_full_org_vec(mtx,nz.col,arr,mtx_ALL_ROWS,TRUE);
137 arr[org_col] -= sum;
138 }
139 }
140 } else { /* arr is indexed by original row number */
141 /* we are working from the bottom up */
142 mtx = sys->factors;
143 for( nz.row = sys->rng.low+sys->rank-1; nz.row >= dotlim ; --(nz.row) ) {
144 register int32 org_row;
145
146 org_row = mtx_row_to_org(mtx,nz.row);
147 if (arr[org_row]!=D_ZERO) nonzero_found=TRUE;
148 if (nonzero_found) {
149 sum= mtx_row_dot_full_org_vec(mtx,nz.row,arr,mtx_ALL_COLS,TRUE);
150 arr[org_row] = (arr[org_row] - sum) / pivlist[nz.row];
151 }
152 }
153 }
154 }
155
156 static int gauss2_solve(linsolqr_system_t sys, struct rhs_list *rl)
157 {
158 /* zero any unsolved for vars first so they don't contaminate
159 mtx_ALL_*O*S dot products.
160 */
161 zero_unpivoted_vars(sys,rl->varvalue,rl->transpose);
162 forward_substitute_gauss2(sys,rl->varvalue,rl->transpose);
163 backward_substitute_gauss2(sys,rl->varvalue,rl->transpose);
164 return 0;
165 }
166
167
168
169 /* this function eliminates a column */
170 static
171 int eliminate_column2(
172 mtx_matrix_t mtx, /* the matrix being factored (upper triangularized) */
173 mtx_matrix_t lmtx, /* the multipliers of the factorization, L */
174 int32 column, /* the column being eliminated (the pivot is (column,column))*/
175 int32 lastrow, /* rows from column+1 to lastrow will be eliminated in column*/
176 real64 pivot, /* the pivot value */
177 mtx_sparse_t *pivotrow /* the pivot row, less the pivot entry, orgindexed */
178 )
179 {
180 mtx_coord_t nz;
181 mtx_range_t elimrng;
182 real64 multiplier;
183
184 /* for rows in column+1 .. lastrow */
185 elimrng.low = column+1;
186 elimrng.high = lastrow;
187
188 nz.row = mtx_FIRST;
189 while ((void)mtx_next_in_col(mtx,&nz,&elimrng),
190 nz.row !=mtx_LAST) {
191 multiplier = mtx_value(mtx,&nz)/pivot;
192 mtx_add_row_sparse(mtx,nz.row,-multiplier,pivotrow);
193 mtx_fill_value(lmtx,&nz,multiplier);
194 }
195 mtx_clear_col(mtx,column,&elimrng);
196 return 0;
197 }
198
199 /* This factors a square region given in sys->rng.
200 * We might like to have a wider pivot region later, so does not assume
201 * that the region to the right of the factorization is empty.
202 * Assumes the region to the left of the factorization IS empty.
203 */
204 int gaussba2_factor_square(linsolqr_system_t sys)
205 {
206 mtx_coord_t nz;
207 int32 lastrow,lastcolumn,column;
208 mtx_range_t pivot_candidates;
209 real64 *tmp;
210 real64 pivot, *pivots;
211 int32 length, defect = 0;
212 mtx_matrix_t mtx, lmtx;
213 mtx_sparse_t *pivotrow;
214 int err = 0;
215
216 length = sys->rng.high - sys->rng.low + 1; /* ? */
217 sys->smallest_pivot = MAXDOUBLE;
218 lastrow = lastcolumn = pivot_candidates.high = sys->rng.high;
219 pivot_candidates.low = sys->rng.low;
220
221 tmp = sys->ludata->tmp;
222 pivots = sys->ludata->pivlist;
223 mtx = sys->factors;
224 lmtx = sys->inverse;
225 pivotrow = sys->ludata->tmp_sparse;
226
227 /* need pivotrow allocation */
228
229 for (column = sys->rng.low; column <= lastcolumn; ) {
230 /* pick pivot. */
231 nz.row = column;
232 pivot = mtx_get_pivot_col(mtx,&nz,&pivot_candidates,&(pivots[nz.row]),
233 sys->ptol,sys->pivot_zero);
234
235 if ( pivot < sys->pivot_zero ) { /* pivot is an epsilon */
236 /* If no pivot, drag to bottom and lastcolumn-- */
237 mtx_drag(mtx, nz.row, lastcolumn);
238 number_drag(pivots, nz.row, lastcolumn);
239 lastcolumn--;
240 pivot_candidates.high--; /* is this correct? */
241 defect++;
242 } else {
243 /* collect pivot row. pivot removed to reduce numerics cost*/
244 mtx_clear_coord(mtx,nz.row,nz.col);
245 pivotrow = mtx_org_row_sparse(mtx,nz.row,pivotrow,
246 mtx_ALL_COLS,mtx_IGNORE_ZEROES); /* check this */
247 /* eliminate column */
248 #if (GAUSS_DEBUG)
249 {
250 FILE *fp;
251 char fname[80];
252
253 sprintf(fname,"/tmp/lu/uga.%d",column");
254 fp = fopen(fname,"w+");
255 mtx_write_region_human_rows(fp,mtx,mtx_ENTIRE_MATRIX);
256 fclose(fp);
257
258 sprintf(fname,"/tmp/lu/lga.%d",column");
259 fp = fopen(fname,"w+");
260 mtx_write_region_human_rows(fp,lmtx,mtx_ENTIRE_MATRIX);
261 fclose(fp);
262
263 sprintf(fname,"/tmp/lu/aga.plot%d",column);
264 fp=fopen(fname,"w+");
265 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
266 fclose(fp);
267 }
268 #endif
269 err = eliminate_column2(mtx,lmtx,column,lastrow,pivot,pivotrow);
270 if (err != 0) {
271 break;
272 }
273 pivot_candidates.low++;
274 column++;
275 }
276 }
277 sys->rank = pivot_candidates.high - sys->rng.low + 1;
278 return err; /* err = 0 is good */
279 }
280
281
282 int gauss2_entry(linsolqr_system_t sys, mtx_region_t *region)
283 {
284 struct rhs_list *rl;
285 double time;
286
287 CHECK_SYSTEM(sys);
288 if (sys->factored) return 0;
289 switch(sys->fmethod) {
290 case gauss_ba2: /* add new methods here */
291 break;
292 default:
293 return 1;
294 }
295
296 if (ISNULL(sys->ludata)) return 1;
297 if (NOTNULL(sys->inverse)) mtx_destroy(sys->inverse);
298 sys->inverse = NULL;
299 if (NOTNULL(sys->factors)) mtx_destroy(sys->factors);
300 if (region == mtx_ENTIRE_MATRIX) determine_pivot_range(sys);
301 else square_region(sys,region);
302
303 sys->factors = mtx_copy_region(sys->coef,region);
304 sys->inverse = mtx_create_slave(sys->factors);
305 sys->rank = -1;
306 sys->smallest_pivot = MAXDOUBLE;
307 for (rl = sys->rl ; NOTNULL(rl) ; rl = rl->next)
308 rl->solved = FALSE;
309 insure_capacity(sys);
310 insure_lu_capacity(sys);
311
312 time = tm_cpu_time();
313 switch(sys->fmethod) {
314 case gauss_ba2:
315 gaussba2_factor_square(sys);
316 break;
317 default:
318 return 1;
319 }
320 sys->factored = TRUE;
321
322 #define KAA_DEBUG
323 #ifdef KAA_DEBUG
324 time = tm_cpu_time() - time;
325 FPRINTF(stderr,"Time for Inversion = %f\n",time);
326 FPRINTF(stderr,"Non-zeros in A = %d\n",
327 mtx_nonzeros_in_region(sys->coef,region));
328 FPRINTF(stderr,"Non-zeros in U+L = %d\n",
329 mtx_nonzeros_in_region(sys->factors,region) +
330 mtx_nonzeros_in_region(sys->inverse,0));
331 #endif /* KAA_DEBUG */
332 #undef KAA_DEBUG
333 return 0;
334 }
335
336
337
338
339
340
341
342
343
344 /* The following crap is from kirk. we should be ditching it and writing a
345 * real gauss suitable to ken's and other needs.
346 */
347 #define GAUSS_ELIMINATE_ROWS_IMPLEMENTED 0
348 #if GAUSS_ELIMINATE_ROWS_IMPLEMENTED
349 static void gauss_eliminate_rows(mtx_matrix_t mtx,
350 mtx_matrix_t lower_mtx,
351 int32 col,
352 real64 *value,
353 int32 *index,
354 real64 pivot)
355 {
356 mtx_coord_t nz;
357 real64 multiplier;
358 int32 i, ii;
359 int colcount, row;
360
361 /*
362 colcount = mtx_cur_col_to_csr(mtx,col,value,index);
363 */
364 mtx_clear_col(mtx,col,mtx_ALL_ROWS);
365
366 row = col;
367 nz.col = col;
368 for (ii=0; ii<colcount; ii++) {
369 i = index[ii];
370 if (i <= row) {
371 /* nothing to do - we would be back patching U here. */
372 continue;
373 }
374 else{
375 multiplier = value[ii]/pivot;
376 mtx_add_row(mtx,row,i,-multiplier,mtx_ALL_COLS);
377 nz.row = i;
378 mtx_fill_value(lower_mtx,&nz,multiplier);
379 }
380 }
381 }
382
383
384 /*
385 * At the moment this code assumes that we are given a
386 * square region. It also assumes that the submatrix A11 will
387 * be provided; in which case pivot selection should be restricted
388 * to A11. Later on we will modify the code to accept a rectangular
389 * region, we will then make use of the 'last_pivotable column',
390 * and sys->reg. At the moment only sys->rng and A11 are used.
391 */
392 static int kirk_gauss_factor(linsolqr_system_t sys,mtx_region_t *A11)
393 {
394 /*
395 * THIS HAS BUG -- row can not be updated.
396 */
397 mtx_matrix_t mtx, lower_mtx;
398 mtx_coord_t nz;
399 int32 last_col;
400 mtx_range_t pivot_candidates;
401 real64 pivot, abs_pivot, *pivots, max;
402 int32 row,col;
403 int32 *index = NULL; /* used by csr */
404 real64 *value = NULL; /* used by csr */
405 int length;
406
407 mtx = sys->factors;
408 lower_mtx = sys->inverse;
409 pivots = sys->ludata->pivlist;
410
411 last_col = sys->rng.high;
412 pivot_candidates.high = A11->col.high; /* initialize */
413
414 length = sys->rng.high - sys->rng.low + 1;
415 value = ASC_NEW_ARRAY(real64,length);
416 index = ASC_NEW_ARRAY(int32,length);
417
418 /*
419 * First find the maximum element in the col *regardless*. This
420 * will leave nz at the correct position in case we need to a
421 * column swap. We have to adjust the pivot range as we go along.
422 */
423 for (col = sys->rng.low; col < last_col; ) {
424 nz.row = nz.col = col;
425 pivot_candidates.low = col;
426 if (nz.col > A11->col.high) {
427 pivot_candidates.high = sys->rng.high;
428 }
429 pivot = mtx_value(mtx,&nz);
430
431 max = mtx_col_max(mtx, &nz, &pivot_candidates, &pivots[col]);
432 if (max <= sys->pivot_zero) {
433 FPRINTF(stderr,"Matrix appears numerically singular at %d %g %g\n",
434 row,pivot,max);
435 /* patch for the time being */
436 pivot = pivots[col] = 1.0;
437 col++;
438 continue;
439 }
440
441 abs_pivot = fabs(pivot);
442 if ((abs_pivot >= sys->ptol * max)) {
443 pivots[col] = pivot;
444 }
445 else {
446 mtx_swap_rows(mtx, row, nz.row);
447 pivot = pivots[col] = max;
448 }
449
450 gauss_eliminate_rows(mtx, lower_mtx, col, value, index, pivot);
451 col++;
452 }
453
454 if (value) ascfree((char *)value);
455 if (index) ascfree((char *)index);
456
457 sys->rank = last_col - sys->rng.low + 1;
458 return sys->rank;
459 }
460 #endif /* GAUSS_ELIMINATE_ROWS_IMPLEMENTED */

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22