/[ascend]/trunk/ascend4/solver/ascgauss.c
ViewVC logotype

Contents of /trunk/ascend4/solver/ascgauss.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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