/[ascend]/trunk/base/generic/solver/linsol.c
ViewVC logotype

Contents of /trunk/base/generic/solver/linsol.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 101 - (show annotations) (download) (as text)
Sat Dec 10 04:22:07 2005 UTC (15 years, 4 months ago) by jds
File MIME type: text/x-csrc
File size: 49226 byte(s)
A little more progress killing compiler warnings.
1 /*
2 * linsol: Ascend Linear Solver
3 * by Karl Michael Westerberg
4 * Created: 2/6/90
5 * Version: $Revision: 1.8 $
6 * Version control file: $RCSfile: linsol.c,v $
7 * Date last modified: $Date: 2000/01/25 02:26:54 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the SLV solver.
11 *
12 * Copyright (C) 1990 Karl Michael Westerberg
13 * Copyright (C) 1993 Joseph Zaher
14 * Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
15 *
16 * The SLV solver is free software; you can redistribute
17 * it and/or modify it under the terms of the GNU General Public License as
18 * published by the Free Software Foundation; either version 2 of the
19 * License, or (at your option) any later version.
20 *
21 * The SLV solver is distributed in hope that it will be
22 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
24 * General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License along with
27 * the program; if not, write to the Free Software Foundation, Inc., 675
28 * Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
29 * COPYING is found in ../compiler.
30 */
31
32 #include <math.h>
33 #include "utilities/ascConfig.h"
34 #include "compiler/compiler.h"
35 #include "utilities/ascMalloc.h"
36 #include "utilities/mem.h"
37 #include "utilities/set.h"
38 #include "general/tm_time.h"
39 #include "solver/mtx.h"
40 #include "solver/linsol.h"
41
42
43 #define LINSOL_DEBUG FALSE
44
45 struct rhs_list {
46 real64 *rhs; /* Vector of rhs values */
47 real64 *varvalue; /* Solution of the linear system */
48 boolean solved; /* ? has the rhs been solved */
49 boolean transpose; /* ? transpose the coef matrix */
50 struct rhs_list *next;
51 };
52
53 struct linsol_header {
54 int integrity;
55 int32 capacity; /* Capacity of arrays */
56 mtx_matrix_t coef; /* Coefficient matrix */
57 struct rhs_list *rl; /* List of rhs vectors */
58 int rlength; /* Length of rhs list */
59 real64 pivot_zero; /* Smallest exceptable pivot */
60 real64 ptol; /* Pivot selection tolerance */
61 mtx_matrix_t inv; /* Inverse matrix containing UL factors */
62 boolean inverted; /* ? has the matrix been inverted */
63 mtx_range_t rng; /* Pivot range */
64 mtx_region_t reg; /* Bounding region */
65 int32 rank; /* Rank of the matrix */
66 real64 smallest_pivot; /* Smallest pivot accepted */
67 };
68
69
70 #define OK ((int)439828676)
71 #define DESTROYED ((int)839276847)
72 static void check_system(linsol_system_t sys)
73 /**
74 *** Checks the system handle
75 **/
76 {
77 if( sys == NULL ) {
78 FPRINTF(stderr,"ERROR: (linsol) check_system\n");
79 FPRINTF(stderr," NULL system handle.\n");
80 return;
81 }
82
83 switch( sys->integrity ) {
84 case OK:
85 break;
86 case DESTROYED:
87 FPRINTF(stderr,"ERROR: (linsol) check_system\n");
88 FPRINTF(stderr," System was recently destroyed.\n");
89 break;
90 default:
91 FPRINTF(stderr,"ERROR: (linsol) check_system\n");
92 FPRINTF(stderr," System was reused or never created.\n");
93 break;
94 }
95 }
96
97 static void destroy_rhs_list(struct rhs_list *rl)
98 /**
99 *** Destroys rhs list.
100 **/
101 {
102 while( rl != NULL ) {
103 struct rhs_list *p;
104 p = rl;
105 rl = rl->next;
106 if( p->varvalue != NULL )
107 ascfree( (POINTER)(p->varvalue) );
108 ascfree( (POINTER)p );
109 }
110 }
111
112 static struct rhs_list *find_rhs(struct rhs_list *rl, real64 *rhs)
113 /**
114 *** Searches for rhs in rhs list, returning it or NULL if not found.
115 **/
116 {
117 for( ; rl != NULL ; rl = rl->next )
118 if( rl->rhs == rhs )
119 break;
120 return(rl);
121 }
122
123
124 linsol_system_t linsol_create()
125 {
126 linsol_system_t sys;
127
128 sys = (linsol_system_t)ascmalloc( sizeof(struct linsol_header) );
129 sys->integrity = OK;
130 sys->capacity = 0;
131 sys->coef = NULL;
132 sys->rl = NULL;
133 sys->rlength = 0;
134 sys->pivot_zero = 1e-12; /* default value */
135 sys->ptol = 0.1; /* default value */
136 sys->inv = NULL;
137 sys->inverted = FALSE;
138 sys->rng.low = sys->rng.high = -1;
139 sys->reg.row.low = sys->reg.row.high = -1;
140 sys->reg.col.low = sys->reg.col.high = -1;
141 sys->rank = -1;
142 sys->smallest_pivot = 0.0;
143 return(sys);
144 }
145
146 void linsol_destroy(linsol_system_t sys)
147 {
148 check_system(sys);
149 if( sys->inv != NULL ) {
150 mtx_destroy(sys->inv);
151 }
152 destroy_rhs_list(sys->rl);
153 sys->integrity = DESTROYED;
154 ascfree( (POINTER)sys );
155 }
156
157 void linsol_set_matrix(linsol_system_t sys, mtx_matrix_t mtx)
158 {
159 check_system(sys);
160 sys->coef = mtx;
161 }
162
163 mtx_matrix_t linsol_get_matrix(linsol_system_t sys)
164 {
165 check_system(sys);
166 return( sys->coef );
167 }
168
169 mtx_matrix_t linsol_get_inverse(linsol_system_t sys)
170 {
171 check_system(sys);
172 return (sys->inv);
173 }
174
175 void linsol_add_rhs(linsol_system_t sys, real64 *rhs, boolean transpose)
176 {
177 struct rhs_list *rl;
178 int32 capacity;
179
180 check_system(sys);
181 rl = find_rhs(sys->rl,rhs);
182 if( rl != NULL ) {
183 return;
184 } else if( rhs != NULL ) {
185 rl = (struct rhs_list *)ascmalloc( sizeof(struct rhs_list) );
186 rl->rhs = rhs;
187 rl->varvalue = NULL;
188 rl->solved = FALSE;
189 rl->transpose = transpose;
190 rl->next = sys->rl;
191 sys->rl = rl;
192 ++(sys->rlength);
193 if (sys->coef) {
194 capacity = mtx_capacity(sys->coef);
195 if (capacity) {
196 rl->varvalue =(real64 *)
197 ascmalloc(capacity*sizeof(real64));
198 }
199 }
200 }
201 }
202
203 void linsol_remove_rhs(linsol_system_t sys, real64 *rhs)
204 {
205 struct rhs_list **q;
206
207 check_system(sys);
208 for( q = &(sys->rl) ; *q != NULL ; q = &((*q)->next) )
209 if( (*q)->rhs == rhs )
210 break;
211 if( *q != NULL ) {
212 struct rhs_list *p;
213 p = *q;
214 *q = p->next;
215 if( p->varvalue != NULL )
216 ascfree( (POINTER)(p->varvalue) );
217 ascfree( (POINTER)p );
218 --(sys->rlength);
219 } else if( rhs != NULL ) {
220 FPRINTF(stderr,"ERROR: (linsol) linsol_remove_rhs\n");
221 FPRINTF(stderr," Rhs does not exist.\n");
222 }
223 }
224
225 int linsol_number_of_rhs(linsol_system_t sys)
226 {
227 check_system(sys);
228 return( sys->rlength );
229 }
230
231 real64 *linsol_get_rhs(linsol_system_t sys, int n)
232 {
233 struct rhs_list *rl;
234 int count;
235
236 check_system(sys);
237
238 count = sys->rlength - 1 - n;
239 if( count < 0 ) return(NULL);
240 for( rl = sys->rl ; count > 0 && rl != NULL ; rl = rl->next )
241 --count;
242 return( rl == NULL ? NULL : rl->rhs );
243 }
244
245 void linsol_matrix_was_changed(linsol_system_t sys)
246 {
247 check_system(sys);
248 sys->inverted = FALSE;
249 }
250
251 void linsol_rhs_was_changed(linsol_system_t sys,real64 *rhs)
252 {
253 struct rhs_list *rl;
254
255 check_system(sys);
256 rl = find_rhs(sys->rl,rhs);
257 if( rl != NULL ) {
258 rl->solved = FALSE;
259 } else if( rhs != NULL ) {
260 FPRINTF(stderr,"ERROR: (linsol) linsol_rhs_was_modified\n");
261 FPRINTF(stderr," Rhs does not exist.\n");
262 }
263 }
264
265 void linsol_set_pivot_zero(linsol_system_t sys,real64 pivot_zero)
266 {
267 check_system(sys);
268 if( pivot_zero < 0.0 ) {
269 FPRINTF(stderr,"ERROR: (linsol) linsol_set_pivot_zero\n");
270 FPRINTF(stderr," Invalid pivot zero of %Lf\n",pivot_zero);
271 } else {
272 sys->pivot_zero = pivot_zero;
273 linsol_matrix_was_changed(sys);
274 }
275 }
276
277 real64 linsol_pivot_zero(linsol_system_t sys)
278 {
279 check_system(sys);
280 return( sys->pivot_zero );
281 }
282
283 void linsol_set_pivot_tolerance(linsol_system_t sys,real64 ptol)
284 {
285 check_system(sys);
286 if( ptol < 0.0 || ptol >= 1.0 ) {
287 FPRINTF(stderr,"ERROR: (linsol) linsol_set_pivot_tolerance\n");
288 FPRINTF(stderr," Invalid pivot tolerance of %Lf\n",ptol);
289 } else {
290 sys->ptol = ptol;
291 linsol_matrix_was_changed(sys);
292 }
293 }
294
295 real64 linsol_pivot_tolerance(linsol_system_t sys)
296 {
297 check_system(sys);
298 return( sys->ptol );
299 }
300
301
302 struct reorder_list { /* List of rows/columns and their counts. */
303 int32 ndx;
304 int32 count;
305 struct reorder_list *next;
306 };
307
308 struct reorder_vars {
309 mtx_matrix_t mtx;
310 mtx_region_t reg; /* Current active region */
311 int32 colhigh; /* Original value of reg.col.high */
312 struct reorder_list *tlist; /* Array of (enough) list elements */
313 int32 *rowcount; /* rowcount[reg.row.low .. reg.row.high] */
314 };
315
316 static void adjust_row_count(struct reorder_vars *vars,int32 removed_col)
317 /**
318 *** Adjusts the row counts to account for the (to be) removed column.
319 **/
320 {
321 mtx_coord_t nz;
322 real64 value;
323 nz.col = removed_col;
324 nz.row = mtx_FIRST;
325 while( value = mtx_next_in_col(vars->mtx,&nz,&(vars->reg.row)),
326 nz.row != mtx_LAST )
327 --(vars->rowcount[nz.row]);
328 }
329
330 static void assign_row_and_col(struct reorder_vars *vars,int32 row,int32 col)
331 /**
332 *** Assigns the given row to the given column, moving the row and column
333 *** to the beginning of the active region and removing them (readjusting
334 *** the region). The row counts are NOT adjusted. If col == mtx_NONE,
335 *** then no column is assigned and the row is moved to the end of the
336 *** active block instead. Otherwise, it is assumed that the column
337 *** is active.
338 **/
339 {
340 if( col == mtx_NONE ) {
341 mtx_swap_rows(vars->mtx,row,vars->reg.row.high);
342 vars->rowcount[row] = vars->rowcount[vars->reg.row.high];
343 --(vars->reg.row.high);
344 } else {
345 mtx_swap_rows(vars->mtx,row,vars->reg.row.low);
346 vars->rowcount[row] = vars->rowcount[vars->reg.row.low];
347 mtx_swap_cols(vars->mtx,col,vars->reg.col.low);
348 ++(vars->reg.row.low);
349 ++(vars->reg.col.low);
350 }
351 }
352
353 static void push_column_on_stack(struct reorder_vars *vars,int32 col)
354 /**
355 *** Pushes the given column onto the stack. It is assumed that the
356 *** column is active. Row counts are adjusted.
357 **/
358 {
359 adjust_row_count(vars,col);
360 mtx_swap_cols(vars->mtx,col,vars->reg.col.high);
361 --(vars->reg.col.high);
362 }
363
364 static int32 pop_column_from_stack(struct reorder_vars *vars)
365 /**
366 *** Pops the column on the "top" of the stack off of the stack and
367 *** returns the column index, where it now lies in the active region.
368 *** If the stack is empty, mtx_NONE is returned. Row counts are NOT
369 *** adjusted (this together with a subsequent assignment of this column
370 *** ==> no row count adjustment necessary).
371 **/
372 {
373 if( vars->reg.col.high < vars->colhigh )
374 return(++(vars->reg.col.high));
375 else
376 return( mtx_NONE );
377 }
378
379 static void assign_null_rows(struct reorder_vars *vars)
380 /**
381 *** Assigns empty rows, moving them to the assigned region. It is
382 *** assumed that row counts are correct. Columns are assigned off the
383 *** stack.
384 **/
385 {
386 int32 row;
387
388 for( row = vars->reg.row.low ; row <= vars->reg.row.high ; ++row )
389 if( vars->rowcount[row] == 0 )
390 assign_row_and_col(vars , row , pop_column_from_stack(vars));
391 }
392
393 static void forward_triangularize(struct reorder_vars *vars)
394 /**
395 *** Forward triangularizes the region, assigning singleton rows with their
396 *** one and only incident column until there are no more. The row counts
397 *** must be correct, and they are updated.
398 **/
399 {
400 boolean change;
401
402 do {
403 mtx_coord_t nz;
404 real64 value;
405 change = FALSE;
406 for( nz.row = vars->reg.row.low ;
407 nz.row <= vars->reg.row.high && vars->rowcount[nz.row] != 1;
408 ++nz.row ) ;
409 if( nz.row <= vars->reg.row.high ) {
410 /* found singleton row */
411 nz.col = mtx_FIRST;
412 value = mtx_next_in_row(vars->mtx,&nz,&(vars->reg.col));
413 adjust_row_count(vars,nz.col);
414 assign_row_and_col(vars,nz.row,nz.col);
415 change = TRUE;
416 }
417 } while( change );
418 }
419
420 static int32 select_row(struct reorder_vars *vars)
421 /**
422 *** Selects a row and returns its index. It is assumed that there is a
423 *** row. Row counts must be correct. vars->tlist will be used.
424 **/
425 {
426 int32 min_row_count;
427 int32 max_col_count;
428 int32 row;
429 int32 i, nties = 0; /* # elements currently defined in vars->tlist */
430
431 /* Set to something > any possible value */
432 min_row_count = vars->reg.col.high-vars->reg.col.low+2;
433 for( row = vars->reg.row.low ; row <= vars->reg.row.high ; ++row )
434 if( vars->rowcount[row] <= min_row_count ) {
435 if( vars->rowcount[row] < min_row_count ) {
436 min_row_count = vars->rowcount[row];
437 nties = 0;
438 }
439 vars->tlist[nties++].ndx = row;
440 }
441 /**
442 *** vars->tlist[0..nties-1] is a list of row numbers which tie for
443 *** minimum row count.
444 **/
445
446 max_col_count = -1; /* < any possible value */
447 i = 0;
448 while( i < nties ) {
449 int32 sum;
450 mtx_coord_t nz;
451 real64 value;
452
453 sum = 0;
454 nz.row = vars->tlist[i].ndx;
455 nz.col = mtx_FIRST;
456 while( value = mtx_next_in_row(vars->mtx,&nz,&(vars->reg.col)),
457 nz.col != mtx_LAST )
458 sum += mtx_nonzeros_in_col(vars->mtx,nz.col,&(vars->reg.row));
459 if( sum > max_col_count ) {
460 max_col_count = sum;
461 row = nz.row;
462 }
463 i++;
464 }
465 /**
466 *** Now row contains the row with the minimum row count, which has the
467 *** greatest total column count of incident columns among all rows with
468 *** the same (minimum) row count. Select it.
469 **/
470 return(row);
471 }
472
473 static void reorder(struct reorder_vars *vars)
474 /**
475 *** Reorders the assigned matrix vars->mtx within the specified bounding
476 *** block region vars->reg. The region is split into 6 subregions during
477 *** reordering: the rows are divided in two, and the columns divided in
478 *** three. Initially everything is in the active subregion. Ultimately,
479 *** everything will be assigned.
480 ***
481 *** <-- assigned -->|<-- active-->|<-- on stack -->|
482 *** ----+----------------+-------------+----------------+
483 *** a | | | |
484 *** s | | | |
485 *** s | | | |
486 *** i | (SQUARE) | | |
487 *** g | | | |
488 *** n | | | |
489 *** e | | | |
490 *** d | | | |
491 *** ----+----------------+-------------+----------------+
492 *** a | | | |
493 *** c | | ACTIVE | |
494 *** t | | REGION | |
495 *** i | | | |
496 *** v | | | |
497 *** e | | | |
498 *** ----+----------------+-------------+----------------+
499 ***
500 *** The algorithm is roughly as follows:
501 *** (1) select a row (by some criterion).
502 *** (2) push columns incident on that row onto the stack in decreasing
503 *** order of their length.
504 *** (3) pop first column off the stack and assign it to the selected
505 *** row.
506 *** (4) forward-triangularize (assign singleton rows with their one
507 *** and only incident column, until no longer possible).
508 ***
509 *** (1)-(4) should be repeated until the active subregion becomes empty.
510 ***
511 *** Everything above was written as though the entire matrix is
512 *** involved. In reality, only the relevant square region is involved.
513 **/
514 {
515 int32 row, size;
516 int32 *rowcount_array_origin;
517
518 size = vars->reg.row.high - vars->reg.row.low + 1;
519 size = MAX(size,vars->reg.col.high - vars->reg.col.low + 1);
520 vars->tlist = size > 0 ? (struct reorder_list *)
521 ascmalloc( size*sizeof(struct reorder_list) ) : NULL;
522 rowcount_array_origin = size > 0 ? (int32 *)
523 ascmalloc( size*sizeof(int32) ) : NULL;
524 vars->rowcount = rowcount_array_origin != NULL ?
525 rowcount_array_origin - vars->reg.row.low : NULL;
526
527 vars->colhigh = vars->reg.col.high;
528 /* Establish row counts */
529 for( row = vars->reg.row.low ; row <= vars->reg.row.high ; ++row )
530 vars->rowcount[row] =
531 mtx_nonzeros_in_row(vars->mtx,row,&(vars->reg.col));
532
533 while(TRUE) {
534 struct reorder_list *head;
535 int32 nelts; /* # elements "allocated" from vars->tlist */
536 mtx_coord_t nz;
537 real64 value;
538
539 forward_triangularize(vars);
540 assign_null_rows(vars);
541 if( vars->reg.row.low>vars->reg.row.high ||
542 vars->reg.col.low>vars->reg.col.high ) {
543 /* Active region is now empty, done */
544 if( vars->tlist != NULL )
545 ascfree( vars->tlist );
546 if( rowcount_array_origin != NULL )
547 ascfree( rowcount_array_origin );
548 return;
549 }
550
551 head = NULL;
552 nelts = 0;
553 nz.row = select_row(vars);
554 nz.col = mtx_FIRST;
555 while( value = mtx_next_in_row(vars->mtx,&nz,&(vars->reg.col)),
556 nz.col != mtx_LAST ) {
557 struct reorder_list **q,*p;
558
559 p = &(vars->tlist[nelts++]);
560 p->ndx = mtx_col_to_org(vars->mtx,nz.col);
561 p->count = mtx_nonzeros_in_col(vars->mtx,nz.col,&(vars->reg.row));
562 for( q = &head; *q && (*q)->count>p->count; q = &((*q)->next) );
563 p->next = *q;
564 *q = p;
565 }
566 /**
567 *** We now have a list of columns which intersect the selected row.
568 *** The list is in order of decreasing column count.
569 **/
570
571 /* Push incident columns on stack */
572 for( ; head != NULL ; head = head->next )
573 push_column_on_stack(vars,mtx_org_to_col(vars->mtx,head->ndx));
574
575 /* Pop column off stack and assign to selected row */
576 assign_row_and_col(vars , nz.row , pop_column_from_stack(vars));
577 }
578 /* Not reached. */
579 }
580
581 static boolean nonempty_row(mtx_matrix_t mtx,int32 row)
582 /**
583 *** ? row not empty in mtx.
584 **/
585 {
586 mtx_coord_t nz;
587 real64 value;
588 value = mtx_next_in_row(mtx, mtx_coord(&nz,row,mtx_FIRST), mtx_ALL_COLS);
589 return( nz.col != mtx_LAST );
590 }
591
592 static boolean nonempty_col(mtx_matrix_t mtx,int32 col)
593 /**
594 *** ? column not empty in mtx.
595 **/
596 {
597 mtx_coord_t nz;
598 real64 value;
599 value = mtx_next_in_col(mtx, mtx_coord(&nz,mtx_FIRST,col), mtx_ALL_ROWS);
600 return( nz.row != mtx_LAST );
601 }
602
603 static void determine_pivot_range(linsol_system_t sys)
604 /**
605 *** Calculates sys->rng from sys->coef.
606 **/
607 {
608 sys->reg.row.low = sys->reg.row.high = -1;
609 sys->reg.col.low = sys->reg.col.high = -1;
610
611 for( sys->rng.high=mtx_order(sys->coef) ; --(sys->rng.high) >= 0 ; ) {
612 if( nonempty_row(sys->coef,sys->rng.high) &&
613 sys->reg.row.high < 0 )
614 sys->reg.row.high = sys->rng.high;
615 if( nonempty_col(sys->coef,sys->rng.high) &&
616 sys->reg.col.high < 0 )
617 sys->reg.col.high = sys->rng.high;
618 if( nonempty_row(sys->coef,sys->rng.high) &&
619 nonempty_col(sys->coef,sys->rng.high) )
620 break;
621 }
622
623 for( sys->rng.low=0 ; sys->rng.low <= sys->rng.high ; ++(sys->rng.low) ) {
624 if( nonempty_row(sys->coef,sys->rng.low) &&
625 sys->reg.row.low < 0 )
626 sys->reg.row.low = sys->rng.low;
627 if( nonempty_col(sys->coef,sys->rng.low) &&
628 sys->reg.col.low < 0 )
629 sys->reg.col.low = sys->rng.low;
630 if( nonempty_row(sys->coef,sys->rng.low) &&
631 nonempty_col(sys->coef,sys->rng.low) )
632 break;
633 }
634 }
635
636 void linsol_reorder(linsol_system_t sys,mtx_region_t *region)
637 /**
638 *** The region to reorder is first isolated by truncating the region
639 *** provided to the largest square region confined to the matrix diagonal.
640 *** It is presumed it will contain no empty rows or columns and will
641 *** provide the basis of candidate pivots when inverting.
642 **/
643 {
644 struct reorder_vars vars;
645 check_system(sys);
646 if( region == mtx_ENTIRE_MATRIX )
647 determine_pivot_range(sys);
648 else {
649 sys->reg.row.low = region->row.low;
650 sys->reg.row.high = region->row.high;
651 sys->reg.col.low = region->col.low;
652 sys->reg.col.high = region->col.high;
653 sys->rng.low = MAX(region->row.low,region->col.low);
654 sys->rng.high = MIN(region->row.high,region->col.high);
655 }
656 vars.mtx = sys->coef;
657 vars.reg.row.low = vars.reg.col.low = sys->rng.low;
658 vars.reg.row.high = vars.reg.col.high = sys->rng.high;
659 reorder(&vars);
660 }
661
662
663 static void drag(mtx_matrix_t mtx,int32 n1,int32 n2)
664 /**
665 *** Drags row n1 to n2, moving row n1 to n2 and shifting all rows in
666 *** between back toward n1. Ditto for columns. This function works
667 *** regardless of whether n1 < n2 or n2 < n1.
668 **/
669 {
670 while( n1 < n2 ) {
671 mtx_swap_rows(mtx,n1,n1+1);
672 mtx_swap_cols(mtx,n1,n1+1);
673 ++n1;
674 }
675
676 while( n1 > n2 ) {
677 mtx_swap_rows(mtx,n1,n1-1);
678 mtx_swap_cols(mtx,n1,n1-1);
679 --n1;
680 }
681 }
682
683 static void eliminate_row(mtx_matrix_t mtx,mtx_range_t *rng,int32 row,real64 *tmp)
684 /**
685 *** Eliminates the given row to the left of the diagonal element, assuming
686 *** valid pivots for all of the diagonal elements above it (the elements
687 *** above those diagonal elements, if any exist, are assumed to be U
688 *** elements and therefore ignored). The temporary array is used by this
689 *** function to do its job. tmp[k], where rng->low <= k <= rng->high must
690 *** be defined (allocated) but need not be initialized.
691 **/
692 {
693 mtx_coord_t nz;
694 real64 value;
695 mtx_range_t left_to_eliminate,high_cols;
696
697 high_cols.low = row;
698 high_cols.high = rng->high;
699 left_to_eliminate.low = rng->low;
700 left_to_eliminate.high = row - 1;
701
702 /* Move non-zeros left of pivot from matrix to full array */
703 mem_zero_byte_cast(tmp+rng->low,0.0,(row-rng->low)*sizeof(real64));
704 nz.row = row;
705 nz.col = mtx_FIRST;
706 while( value = mtx_next_in_row(mtx,&nz,&left_to_eliminate),
707 nz.col != mtx_LAST )
708 tmp[nz.col] = value;
709 mtx_clear_row(mtx,row,&left_to_eliminate);
710
711 /* eliminates nzs from pivot, one by one, filling tmp with multipliers */
712 for( nz.row = row-1 ; nz.row >= rng->low ; --(nz.row) ) {
713 if( tmp[nz.row] == 0.0 )
714 continue; /* Nothing to do for this row */
715
716 nz.col = nz.row;
717 tmp[nz.row] /= mtx_value(mtx,&nz);
718 /* tmp[nz.row] now equals multiplier */
719
720 /* Perform "mtx_add_row" for full array part of the row */
721 left_to_eliminate.high = nz.row - 1;
722 nz.col = mtx_FIRST;
723 while( value = mtx_next_in_row(mtx,&nz,&left_to_eliminate),
724 nz.col != mtx_LAST )
725 tmp[nz.col] -= tmp[nz.row] * value;
726
727 /* Perform "mtx_add_row" on remaining part of row */
728 mtx_add_row(mtx,nz.row,row,-tmp[nz.row],&high_cols);
729 }
730
731 nz.row = row;
732 for( nz.col = rng->low ; nz.col < row ; ++(nz.col) )
733 if( tmp[nz.col] != 0.0 ) mtx_add_value(mtx,&nz,tmp[nz.col]);
734 }
735
736 #ifdef THIS_IS_AN_UNUSED_FUNCTION
737 static int32 top_of_spike(linsol_system_t sys,int32 col)
738 /**
739 *** Determines the top row (row of lowest index) in a possible spike
740 *** above the diagonal element in the given column. If there is no spike,
741 *** then (row = ) col is returned.
742 **/
743 {
744 mtx_range_t above_diagonal;
745 mtx_coord_t nz;
746 real64 value;
747 int32 top_row;
748
749 above_diagonal.low = sys->rng.low;
750 above_diagonal.high = col-1;
751 top_row = nz.col = col;
752 nz.row = mtx_FIRST;
753 while( value = mtx_next_in_col(sys->inv,&nz,&above_diagonal),
754 nz.row != mtx_LAST )
755 if( nz.row < top_row )
756 top_row = nz.row;
757 return( top_row );
758 }
759 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
760
761 static boolean col_is_a_spike(linsol_system_t sys,int32 col)
762 /**
763 *** Determines if the col is a spike, characterized by having any
764 *** nonzeros above the diagonal.
765 **/
766 {
767 mtx_range_t above_diagonal;
768 mtx_coord_t nz;
769 real64 value;
770 int32 top_row;
771
772 above_diagonal.low = sys->rng.low;
773 above_diagonal.high = col-1;
774 top_row = nz.col = col;
775 nz.row = mtx_FIRST;
776 while( value = mtx_next_in_col(sys->inv,&nz,&above_diagonal),
777 nz.row != mtx_LAST )
778 if( nz.row < col ) return TRUE;
779 return( FALSE );
780 }
781
782 static void invert(linsol_system_t sys)
783 /**
784 *** This is the heart of the linear equation solver. This function
785 *** factorizes the matrix into a lower (L) and upper (U) triangular
786 *** matrix. sys->smallest_pivot and sys->rank are calculated. The
787 *** RANKI method is utilized. At the end of elimination, the matrix A
788 *** is factored into A = U L, where U and L are stored as follows:
789 ***
790 *** <----- r ----> <- n-r ->
791 *** +--------------+---------+
792 *** | | |
793 *** | U | |
794 *** | | |
795 *** | L | | r
796 *** | | |
797 *** +--------------+---------+
798 *** | | |
799 *** | | 0 | n-r
800 *** | | |
801 *** +--------------+---------+
802 ***
803 *** The rows and columns have been permuted so that all of the pivoted
804 *** original rows and columns are found in the first r rows and columns
805 *** of the region. The last n-r rows and columns are unpivoted. U has
806 *** 1's on its diagonal, whereas L's diagonal is stored in the matrix.
807 ***
808 *** Everything above was written as though the entire matrix is
809 *** involved. In reality, only the relevant square region is involved.
810 **/
811 {
812 mtx_coord_t nz;
813 int32 last_row;
814 mtx_range_t pivot_candidates;
815 real64 *tmp,*tmp_array_origin;
816 int32 length;
817
818 length = sys->rng.high - sys->rng.low + 1;
819 tmp_array_origin = length > 0 ?
820 (real64 *)ascmalloc( length*sizeof(real64) ) : NULL;
821 tmp = tmp_array_origin != NULL ?
822 tmp_array_origin - sys->rng.low : NULL;
823
824 sys->smallest_pivot = MAXDOUBLE;
825 last_row = pivot_candidates.high = sys->rng.high;
826 for( nz.row = sys->rng.low ; nz.row <= last_row ; ) {
827 real64 pivot;
828
829 pivot_candidates.low = nz.col = nz.row;
830 pivot = mtx_value(sys->inv,&nz);
831 pivot = fabs(pivot);
832 if( pivot > sys->pivot_zero &&
833 pivot > sys->ptol * mtx_row_max(sys->inv,&nz,&pivot_candidates,NULL)
834 && !col_is_a_spike(sys,nz.row) ) {
835 /* Good pivot and not a spike: continue with next row */
836 if( pivot < sys->smallest_pivot )
837 sys->smallest_pivot = pivot;
838 ++(nz.row);
839 continue;
840 }
841
842 /**
843 *** Row is a spike row or will
844 *** be when a necessary column
845 *** exchange occurs.
846 **/
847 eliminate_row(sys->inv,&(sys->rng),nz.row,tmp);
848 if( (pivot=mtx_row_max(sys->inv,&nz,&pivot_candidates,NULL)) <=
849 sys->pivot_zero ) {
850 /* Dependent row, drag to the end */
851 drag(sys->inv,nz.row,last_row);
852 --last_row;
853 } else {
854 /* Independent row: nz contains best pivot */
855 mtx_swap_cols(sys->inv,nz.row,nz.col);
856 /* Move pivot to diagonal */
857 drag( sys->inv , nz.row , sys->rng.low );
858 if( pivot < sys->smallest_pivot )
859 sys->smallest_pivot = pivot;
860 ++(nz.row);
861 }
862
863 }
864 if( tmp_array_origin != NULL )
865 ascfree( (POINTER)tmp_array_origin );
866
867 sys->rank = last_row - sys->rng.low + 1;
868 }
869
870 static real64 *raise_capacity(real64 *vec,int32 oldcap,int32 newcap)
871 /**
872 *** Raises the capacity of the array and returns a new array.
873 *** It is assumed that oldcap < newcap. vec is destroyed or
874 *** returned as appropriate.
875 **/
876 {
877 real64 *newvec=NULL;
878 /*
879 real64 *newvec;
880 newvec = newcap > 0 ?
881 (real64 *)ascmalloc( newcap * sizeof(real64) ) : NULL;
882 if( vec != NULL ) {
883 mem_move_cast(vec,newvec,oldcap*sizeof(real64));
884 ascfree( (POINTER)vec );
885 }
886 return(newvec);
887 */
888 if (newcap < oldcap)
889 return vec;
890 if (vec!=NULL) /* don't call realloc on null with newcap 0 or it frees */
891 newvec = (real64 *)ascrealloc(vec,(newcap * sizeof(real64)));
892 else
893 newvec=newcap > 0 ?
894 (real64 *)ascmalloc( newcap * sizeof(real64) ) : NULL;
895 return newvec;
896 }
897
898 static void insure_capacity(linsol_system_t sys)
899 /**
900 *** Insures that the capacity of all of the solution vectors
901 *** for each rhs is large enough.
902 **/
903 {
904 int32 req_cap = mtx_capacity(sys->coef);
905
906 if( req_cap > sys->capacity ) {
907 struct rhs_list *rl;
908
909 for( rl = sys->rl ; rl != NULL ; rl = rl->next )
910 rl->varvalue = raise_capacity(rl->varvalue,sys->capacity,req_cap);
911 sys->capacity = req_cap;
912 }
913 }
914
915 static void forward_substitute(linsol_system_t sys,real64 *arr,boolean transpose)
916 /**
917 *** Forward substitute. It is assumed that the L (or U) part of
918 *** sys->inv is computed. This function converts c to x in place. The
919 *** values are stored in arr indexed by original row number (or original
920 *** column number).
921 ***
922 *** transpose = FALSE: transpose = TRUE:
923 *** T
924 *** L x = c U x = c
925 ***
926 *** The following formula holds:
927 *** 0<=k<r ==> x(k) = [c(k) - L(k,0..k-1) dot x(0..k-1)] / L(k,k)
928 *** or
929 *** 0<=k<r ==> x(k) = [c(k) - U(0..k-1,k) dot x(0..k-1)] / U(k,k)
930 **/
931 {
932 mtx_range_t dot_rng;
933 mtx_coord_t nz;
934 real64 value;
935
936 dot_rng.low = sys->rng.low;
937 if (transpose) { /* arr is indexed by original column number */
938 for( nz.col=dot_rng.low; nz.col < dot_rng.low+sys->rank; ++(nz.col) ) {
939 real64 sum;
940 register int32 org_col;
941
942 dot_rng.high = nz.col - 1;
943 sum = 0.0;
944 nz.row = mtx_FIRST;
945 while( value = mtx_next_in_col(sys->inv,&nz,&dot_rng),
946 nz.row != mtx_LAST )
947 sum += value*arr[mtx_col_to_org(sys->inv,nz.row)];
948
949 nz.row = nz.col;
950 org_col = mtx_col_to_org(sys->inv,nz.col);
951 arr[org_col] = (arr[org_col] - sum) / 1.0;
952 }
953
954 } else { /* arr is indexed by original row number */
955 for( nz.row=dot_rng.low; nz.row < dot_rng.low+sys->rank; ++(nz.row) ) {
956 real64 sum;
957 register int32 org_row;
958
959 dot_rng.high = nz.row - 1;
960 sum = 0.0;
961 nz.col = mtx_FIRST;
962 while( value = mtx_next_in_row(sys->inv,&nz,&dot_rng),
963 nz.col != mtx_LAST )
964 sum += value*arr[mtx_row_to_org(sys->inv,nz.col)];
965
966 nz.col = nz.row;
967 org_row = mtx_row_to_org(sys->inv,nz.row);
968 arr[org_row] = (arr[org_row] - sum) / mtx_value(sys->inv,&nz);
969 }
970 }
971 }
972
973 static void backward_substitute(linsol_system_t sys,real64 *arr,boolean transpose)
974 /**
975 *** Backward substitute. It is assumed that the U (or L) part of
976 *** sys->inv is computed. This function converts rhs to c in place. The
977 *** values are stored in arr indexed by original row number (or original
978 *** column number).
979 ***
980 *** transpose = FALSE: transpose = TRUE:
981 *** T
982 *** U c = rhs L c = rhs
983 ***
984 *** The following formula holds:
985 *** 0<=k<r ==> c(k) = [rhs(k) - U(k,k+1..r-1) dot c(k+1..r-1)] / U(k,k)
986 *** or
987 *** 0<=k<r ==> c(k) = [rhs(k) - L(k+1..r-1,k) dot c(k+1..r-1)] / L(k,k)
988 **/
989 {
990 mtx_range_t dot_rng;
991 mtx_coord_t nz;
992 real64 value;
993
994 dot_rng.high = sys->rng.low + sys->rank - 1;
995 if (transpose) { /* arr is indexed by original column number */
996 for( nz.col = dot_rng.high ; nz.col >= sys->rng.low ; --(nz.col) ) {
997 real64 sum;
998 register int32 org_col;
999
1000 dot_rng.low = nz.col + 1;
1001 sum = 0.0;
1002 nz.row = mtx_FIRST;
1003 while( value = mtx_next_in_col(sys->inv,&nz,&dot_rng),
1004 nz.row != mtx_LAST )
1005 sum += value*arr[mtx_col_to_org(sys->inv,nz.row)];
1006
1007 nz.row = nz.col;
1008 org_col = mtx_col_to_org(sys->inv,nz.col);
1009 arr[org_col] = (arr[org_col] - sum) / mtx_value(sys->inv,&nz);
1010 }
1011
1012 } else { /* arr is indexed by original row number */
1013 for( nz.row = dot_rng.high ; nz.row >= sys->rng.low ; --(nz.row) ) {
1014 real64 sum;
1015 register int32 org_row;
1016
1017 dot_rng.low = nz.row + 1;
1018 sum = 0.0;
1019 nz.col = mtx_FIRST;
1020 while( value = mtx_next_in_row(sys->inv,&nz,&dot_rng),
1021 nz.col != mtx_LAST )
1022 sum += value*arr[mtx_row_to_org(sys->inv,nz.col)];
1023
1024 nz.col = nz.row;
1025 org_row = mtx_row_to_org(sys->inv,nz.row);
1026 arr[org_row] = (arr[org_row] - sum) / 1.0;
1027 }
1028 }
1029 }
1030
1031
1032 static void calc_dependent_rows(linsol_system_t sys)
1033 {
1034 mtx_coord_t nz;
1035 real64 value;
1036 mtx_range_t colrange;
1037 real64 *lc;
1038
1039 if( (sys->reg.row.low == sys->rng.low) &&
1040 ( sys->reg.row.high == sys->rng.low+sys->rank-1 ) )
1041 return;
1042 if (sys->rank==0) return;
1043
1044 lc = sys->capacity > 0 ? (real64 *)
1045 ascmalloc(sys->capacity*sizeof(real64)) : NULL;
1046 colrange.low = sys->rng.low;
1047 colrange.high = colrange.low + sys->rank - 1;
1048
1049 nz.row = sys->reg.row.low;
1050 for( ; nz.row <= sys->reg.row.high; nz.row++ ) {
1051 if( nz.row == sys->rng.low ) {
1052 nz.row = sys->rng.low+sys->rank-1;
1053 continue;
1054 }
1055 mem_zero_byte_cast(lc,0.0,(sys->capacity)*sizeof(real64));
1056 nz.col = mtx_FIRST;
1057 while( value = mtx_next_in_row(sys->inv,&nz,&colrange),
1058 nz.col != mtx_LAST )
1059 lc[mtx_col_to_org(sys->inv,nz.col)] = value;
1060 if( nz.row < sys->rng.low+sys->rank || nz.row > sys->rng.high )
1061 backward_substitute(sys,lc,TRUE);
1062 forward_substitute(sys,lc,TRUE);
1063 mtx_clear_row(sys->inv,nz.row,&colrange);
1064 for( nz.col=colrange.low; nz.col <= colrange.high; nz.col++ ) {
1065 real64 value;
1066 value = lc[mtx_col_to_org(sys->inv,nz.col)];
1067 if( value != 0.0 ) mtx_add_value(sys->inv,&nz,value);
1068 }
1069 }
1070
1071 if( lc ) ascfree(lc);
1072 }
1073
1074
1075 static void calc_dependent_cols(linsol_system_t sys)
1076 {
1077 mtx_coord_t nz;
1078 real64 value;
1079 mtx_range_t rowrange;
1080 real64 *lc;
1081
1082 if( (sys->reg.col.low == sys->rng.low) &&
1083 ( sys->reg.col.high == sys->rng.low+sys->rank-1 ) )
1084 return;
1085 if (sys->rank==0) return;
1086
1087 lc = sys->capacity > 0 ? (real64 *)
1088 ascmalloc(sys->capacity*sizeof(real64)) : NULL;
1089 rowrange.low = sys->rng.low;
1090 rowrange.high = rowrange.low + sys->rank - 1;
1091
1092 nz.col = sys->reg.col.low;
1093 for( ; nz.col <= sys->reg.col.high; nz.col++ ) {
1094 if( nz.col == sys->rng.low ) {
1095 nz.col = sys->rng.low+sys->rank-1;
1096 continue;
1097 }
1098 mem_zero_byte_cast(lc,0.0,sys->capacity*sizeof(real64));
1099 nz.row = mtx_FIRST;
1100 while( value = mtx_next_in_col(sys->inv,&nz,&rowrange),
1101 nz.row != mtx_LAST )
1102 lc[mtx_row_to_org(sys->inv,nz.row)] = value;
1103 if( nz.col < sys->rng.low+sys->rank || nz.col > sys->rng.high )
1104 backward_substitute(sys,lc,FALSE);
1105 forward_substitute(sys,lc,FALSE);
1106 mtx_clear_col(sys->inv,nz.col,&rowrange);
1107 for( nz.row=rowrange.low; nz.row <= rowrange.high; nz.row++ ) {
1108 real64 value;
1109 value = lc[mtx_row_to_org(sys->inv,nz.row)];
1110 if( value != 0.0 ) mtx_add_value(sys->inv,&nz,value);
1111 }
1112 }
1113
1114 if( lc ) ascfree(lc);
1115 }
1116
1117 #ifdef THIS_IS_AN_UNUSED_FUNCTION
1118 static void debug_out_inverse(FILE *fp,linsol_system_t sys)
1119 /**
1120 *** Outputs permutation and values of the nonzero elements in the
1121 *** inverse matrix.
1122 **/
1123 {
1124 mtx_coord_t nz;
1125 real64 value;
1126
1127 nz.row = sys->rng.low;
1128 for( ; nz.row <= sys->rng.high; ++(nz.row) ) {
1129 FPRINTF(fp," Row %d (org %d)\n",
1130 nz.row,
1131 mtx_row_to_org(sys->inv,nz.row));
1132 nz.col = mtx_FIRST;
1133 while( value = mtx_next_in_row(sys->inv,&nz,&(sys->rng)),
1134 nz.col != mtx_LAST )
1135 FPRINTF(fp," Col %d (org %d) has value %g\n",
1136 nz.col, mtx_col_to_org(sys->inv,nz.col), value);
1137 }
1138 }
1139 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
1140
1141 #ifdef THIS_IS_AN_UNUSED_FUNCTION
1142 /* from joe equivalent to ranki_jz*/
1143 static void invert_alt(linsol_system_t sys)
1144 /**
1145 *** This is the heart of the linear equation solver. This function
1146 *** factorizes the matrix into a lower (L) and upper (U) triangular
1147 *** matrix. sys->smallest_pivot and sys->rank are calculated. The
1148 *** RANKI method is utilized. At the end of elimination, the matrix A
1149 *** is factored into A = U L, where U and L are stored as follows:
1150 ***
1151 *** <----- r ----> <- n-r ->
1152 *** +--------------+---------+
1153 *** | | |
1154 *** | U | |
1155 *** | | |
1156 *** | L | | r
1157 *** | | |
1158 *** +--------------+---------+
1159 *** | | |
1160 *** | | 0 | n-r
1161 *** | | |
1162 *** +--------------+---------+
1163 ***
1164 *** The rows and columns have been permuted so that all of the pivoted
1165 *** original rows and columns are found in the first r rows and columns
1166 *** of the region. The last n-r rows and columns are unpivoted. U has
1167 *** 1's on its diagonal, whereas L's diagonal is stored in the matrix.
1168 ***
1169 *** Everything above was written as though the entire matrix is
1170 *** involved. In reality, only the relevant square region is involved.
1171 **/
1172 {
1173 real64 pivot;
1174 real64 biggest;
1175 mtx_coord_t nz, best;
1176 mtx_region_t candidates;
1177 real64 *tmp;
1178 int32 length;
1179
1180 length = sys->rng.high - sys->rng.low + 1;
1181 tmp = length > 0 ?
1182 (real64 *)ascmalloc( length*sizeof(real64) ) : NULL;
1183
1184 sys->smallest_pivot = MAXDOUBLE;
1185 candidates.row.high = sys->rng.high;
1186 candidates.col.high = sys->rng.high;
1187 for( nz.row = sys->rng.low ; nz.row <= candidates.row.high ; ) {
1188 nz.col = nz.row;
1189 pivot = mtx_value(sys->inv,&nz);
1190 pivot = fabs(pivot);
1191 candidates.row.low = nz.row;
1192 candidates.col.low = nz.row;
1193
1194 if( !col_is_a_spike(sys,nz.row) ) {
1195 best.col = nz.row;
1196 biggest = mtx_col_max(sys->inv,&best,&(candidates.row),NULL);
1197 if( biggest >= sys->pivot_zero ) {
1198 if( pivot < sys->pivot_zero || pivot < sys->ptol*biggest ) {
1199 mtx_swap_rows(sys->inv,nz.row,best.row);
1200 pivot = biggest;
1201 }
1202 if( pivot < sys->smallest_pivot ) sys->smallest_pivot = pivot;
1203 ++(nz.row);
1204 continue;
1205 }
1206 }
1207
1208 /**
1209 *** Row is a spike row or will
1210 *** be when a necessary column
1211 *** exchange occurs.
1212 **/
1213 eliminate_row(sys->inv,&(sys->rng),nz.row,tmp - sys->rng.low);
1214 pivot = mtx_row_max(sys->inv,&nz,&(candidates.col),NULL);
1215 if( pivot < sys->pivot_zero ) { /* pivot is an epsilon */
1216 /* Dependent row, drag nz to lower right */
1217 mtx_drag(sys->inv,nz.row,candidates.row.high);
1218 --(candidates.row.high);
1219 } else {
1220 /* Independent row, drag nz to upper left */
1221 mtx_swap_cols(sys->inv,nz.row,nz.col);
1222 mtx_drag(sys->inv,nz.row,sys->rng.low);
1223 if( pivot < sys->smallest_pivot )
1224 sys->smallest_pivot = pivot;
1225 ++(nz.row);
1226 }
1227 }
1228 if( tmp != NULL ) ascfree( (POINTER)tmp );
1229 sys->rank = candidates.row.high - sys->rng.low + 1;
1230 }
1231 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
1232
1233 #define KAA_DEBUG 0
1234 void linsol_invert(linsol_system_t sys,mtx_region_t *region)
1235 /**
1236 *** The region to invert is first isolated by truncating the region
1237 *** provided to the largest square region confined to the matrix diagonal.
1238 *** It is presumed it will contain no empty rows or columns and that it has
1239 *** been previously reordered using linsol_reorder.
1240 **/
1241 {
1242 struct rhs_list *rl;
1243 #if KAA_DEBUG
1244 double time;
1245 #endif
1246
1247 check_system(sys);
1248 if( sys->inverted )
1249 return;
1250
1251 #if KAA_DEBUG
1252 time = tm_cpu_time();
1253 #endif
1254 if( sys->inv != NULL )
1255 mtx_destroy(sys->inv);
1256 sys->inv = mtx_copy(sys->coef);
1257 sys->rank = -1;
1258 sys->smallest_pivot = 0.0;
1259 for( rl = sys->rl ; rl != NULL ; rl = rl->next )
1260 rl->solved = FALSE;
1261 insure_capacity(sys);
1262 if( region == mtx_ENTIRE_MATRIX )
1263 determine_pivot_range(sys);
1264 else {
1265 sys->reg = *region;
1266 /*
1267 sys->reg.row.low = region->row.low;
1268 sys->reg.row.high = region->row.high;
1269 sys->reg.col.low = region->col.low;
1270 sys->reg.col.high = region->col.high;
1271 */
1272 sys->rng.low = MAX(region->row.low,region->col.low);
1273 sys->rng.high = MIN(region->row.high,region->col.high);
1274 }
1275 invert(sys);
1276 calc_dependent_rows(sys);
1277 calc_dependent_cols(sys);
1278 sys->inverted = TRUE;
1279 #if KAA_DEBUG
1280 time = tm_cpu_time() - time;
1281 FPRINTF(stderr,"Time for Inversion = %f\n",time);
1282 #endif /* KAA_DEBUG */
1283 }
1284
1285 int32 linsol_rank(linsol_system_t sys)
1286 {
1287 check_system(sys);
1288 if( !sys->inverted ) {
1289 FPRINTF(stderr,"ERROR: (linsol) linsol_rank\n");
1290 FPRINTF(stderr," System not inverted yet.\n");
1291 }
1292 return(sys->rank);
1293 }
1294
1295 real64 linsol_smallest_pivot(sys)
1296 linsol_system_t sys;
1297 {
1298 check_system(sys);
1299 #if LINSOL_DEBUG
1300 if( !sys->inverted ) {
1301 FPRINTF(stderr,"ERROR: (linsol) linsol_smallest_pivot\n");
1302 FPRINTF(stderr," System not inverted yet.\n");
1303 }
1304 #endif
1305 return(sys->smallest_pivot);
1306 }
1307
1308 int linsol_get_pivot_sets(linsol_system_t sys,unsigned *org_rowpivots,unsigned *org_colpivots)
1309 {
1310 int32 ndx;
1311
1312 check_system(sys);
1313 if( !sys->inverted ) {
1314 #if LINSOL_DEBUG
1315 FPRINTF(stderr,"ERROR: (linsol) linsol_get_pivot_sets\n");
1316 FPRINTF(stderr," System not inverted yet.\n");
1317 #endif
1318 return 0;
1319 }
1320 for( ndx = sys->rng.low ; ndx < sys->rng.low + sys->rank ; ++ndx ) {
1321 set_change_member(org_rowpivots,mtx_row_to_org(sys->inv,ndx),TRUE);
1322 set_change_member(org_colpivots,mtx_col_to_org(sys->inv,ndx),TRUE);
1323 }
1324 return 1;
1325 }
1326
1327 #define org_row_to_org_col(sys,org_row) \
1328 mtx_col_to_org((sys)->inv,mtx_org_to_row((sys)->inv,org_row))
1329 #define org_col_to_org_row(sys,org_col) \
1330 mtx_row_to_org((sys)->inv,mtx_org_to_col((sys)->inv,org_col))
1331
1332 int32 linsol_org_row_to_org_col(linsol_system_t sys,int32 org_row)
1333 {
1334 check_system(sys);
1335 if( !sys->inverted ) {
1336 FPRINTF(stderr,"ERROR: (linsol) linsol_org_row_to_org_col\n");
1337 FPRINTF(stderr," System not inverted yet.\n");
1338 }
1339 return( org_row_to_org_col(sys,org_row) );
1340 }
1341
1342 int32 linsol_org_col_to_org_row(linsol_system_t sys,int32 org_col)
1343 {
1344 check_system(sys);
1345 if( !sys->inverted ) {
1346 FPRINTF(stderr,"ERROR: (linsol) linsol_org_col_to_org_row\n");
1347 FPRINTF(stderr," System not inverted yet.\n");
1348 }
1349 return( org_col_to_org_row(sys,org_col) );
1350 }
1351
1352 real64 linsol_org_row_dependency(linsol_system_t sys,int32 dep,int32 ind)
1353 {
1354 mtx_coord_t nz;
1355
1356 check_system(sys);
1357 if( !sys->inverted ) {
1358 #if LINSOL_DEBUG
1359 FPRINTF(stderr,"ERROR: (linsol) linsol_org_row_dependency\n");
1360 FPRINTF(stderr," System not inverted yet.\n");
1361 FPRINTF(stderr," Returning 0.0.\n");
1362 #endif
1363 return(0.0);
1364 }
1365
1366 nz.col = mtx_org_to_row(sys->inv,ind);
1367 if( (sys->rng.low > nz.col) || (nz.col > sys->rng.low+sys->rank-1) ) {
1368 FPRINTF(stderr,"ERROR: (linsol) linsol_org_row_dependency\n");
1369 FPRINTF(stderr," Original row %d is not independent.\n", ind);
1370 FPRINTF(stderr," Returning 0.0.\n");
1371 return(0.0);
1372 }
1373
1374 nz.row = mtx_org_to_row(sys->inv,dep);
1375 if( (nz.row <= sys->rng.low+sys->rank-1) && (nz.row >= sys->rng.low) )
1376 return( ind == dep ? 1.0 : 0.0 );
1377
1378 return(mtx_value(sys->inv,&nz));
1379 }
1380
1381 real64 linsol_org_col_dependency(linsol_system_t sys,int32 dep,int32 ind)
1382 {
1383 mtx_coord_t nz;
1384
1385 check_system(sys);
1386 if( !sys->inverted ) {
1387 #if LINSOL_DEBUG
1388 FPRINTF(stderr,"ERROR: (linsol) linsol_org_col_dependency\n");
1389 FPRINTF(stderr," System not inverted yet.\n");
1390 FPRINTF(stderr," Returning 0.0.\n");
1391 return(0.0);
1392 #endif
1393 }
1394
1395 nz.row = mtx_org_to_col(sys->inv,ind);
1396 if( (sys->rng.low > nz.row) || (nz.row > sys->rng.low+sys->rank-1) ) {
1397 FPRINTF(stderr,"ERROR: (linsol) linsol_org_col_dependency\n");
1398 FPRINTF(stderr," Original col %d is not independent.\n",ind);
1399 FPRINTF(stderr," Returning 0.0.\n");
1400 return(0.0);
1401 }
1402
1403 nz.col = mtx_org_to_col(sys->inv,dep);
1404 if( (nz.col <= sys->rng.low+sys->rank-1) && (nz.col >= sys->rng.low) )
1405 return( ind == dep ? 1.0 : 0.0 );
1406
1407 return(mtx_value(sys->inv,&nz));
1408 }
1409
1410
1411 static void zero_unpivoted_vars(linsol_system_t sys,real64 *varvalues,boolean transpose)
1412 /**
1413 *** Sets the values of unpivoted variables to zero.
1414 **/
1415 {
1416 int32 ndx,order;
1417
1418 order = mtx_order(sys->inv);
1419 for( ndx = 0 ; ndx < sys->rng.low ; ++ndx )
1420 if (transpose)
1421 varvalues[mtx_col_to_org(sys->inv,ndx)] = 0.0;
1422 else
1423 varvalues[mtx_row_to_org(sys->inv,ndx)] = 0.0;
1424
1425 for( ndx = sys->rng.low + sys->rank ; ndx < order ; ++ndx )
1426 if (transpose)
1427 varvalues[mtx_col_to_org(sys->inv,ndx)] = 0.0;
1428 else
1429 varvalues[mtx_row_to_org(sys->inv,ndx)] = 0.0;
1430 }
1431
1432 void linsol_solve(linsol_system_t sys,real64 *rhs)
1433 /**
1434 *** Assuming the bounding block region of the matrix has been previously
1435 *** inverted, the specified rhs can then be applied.
1436 ***
1437 *** A x = U L x = rhs. Define c := L x, so that U c = rhs and L x = c.
1438 ***
1439 *** or
1440 *** T T T T T T
1441 *** A x = L U x = rhs. Define c := U x, so that L c = rhs and U x = c.
1442 ***
1443 *** The variables associated with any of the unpivoted original rows
1444 *** and columns are assigned the value of zero by convention.
1445 **/
1446 {
1447 struct rhs_list *rl;
1448
1449 check_system(sys);
1450 if( !sys->inverted ) {
1451 FPRINTF(stderr,"ERROR: (linsol) linsol_solve\n");
1452 FPRINTF(stderr," System not inverted yet.\n");
1453 }
1454 rl = find_rhs(sys->rl,rhs);
1455 if( rl != NULL ) {
1456 if( rl->solved )
1457 return;
1458 /* by my reading, these are disjoint
1459 mem_move_cast(rl->rhs,rl->varvalue,sys->capacity*sizeof(real64));
1460 */
1461 mem_copy_cast(rl->rhs,rl->varvalue,sys->capacity*sizeof(real64));
1462
1463 backward_substitute(sys,rl->varvalue,rl->transpose);
1464 forward_substitute(sys,rl->varvalue,rl->transpose);
1465
1466 zero_unpivoted_vars(sys,rl->varvalue,rl->transpose);
1467 rl->solved = TRUE;
1468 } else if( rhs != NULL ) {
1469 FPRINTF(stderr,"ERROR: (linsol) linsol_solve\n");
1470 FPRINTF(stderr," Rhs does not exist.\n");
1471 }
1472 }
1473
1474 real64 linsol_var_value(linsol_system_t sys, real64 *rhs, int32 ndx)
1475 {
1476 struct rhs_list *rl;
1477
1478 check_system(sys);
1479 if( !sys->inverted ) {
1480 FPRINTF(stderr,"ERROR: (linsol) linsol_var_value\n");
1481 FPRINTF(stderr," System not inverted yet.\n");
1482 }
1483 rl = find_rhs(sys->rl,rhs);
1484 if( rl != NULL ) {
1485 if( !(rl->solved) ) {
1486 FPRINTF(stderr,"ERROR: (linsol) linsol_var_value\n");
1487 FPRINTF(stderr," Rhs not solved yet.\n");
1488 }
1489 if( rl->transpose ) {
1490 /* ndx is an original row index */
1491 return( rl->varvalue[org_row_to_org_col(sys,ndx)] );
1492 } else {
1493 /* ndx is an original column index */
1494 return( rl->varvalue[org_col_to_org_row(sys,ndx)] );
1495 }
1496 } else if( rhs != NULL ) {
1497 FPRINTF(stderr,"ERROR: (linsol) linsol_var_value\n");
1498 FPRINTF(stderr," Rhs does not exist.\n");
1499 }
1500 return 0.0; /* Function had no return value. Added this line. JDS */
1501 }
1502
1503 boolean linsol_copy_solution(linsol_system_t sys,real64 *rhs,real64 *vector)
1504 {
1505 struct rhs_list *rl;
1506 real64 *varvalue;
1507 int ndx,size;
1508
1509 check_system(sys);
1510 if( !sys->inverted ) {
1511 FPRINTF(stderr,"ERROR: (linsol) linsol_var_value\n");
1512 FPRINTF(stderr," System not inverted yet.\n");
1513 return TRUE;
1514 }
1515 rl = find_rhs(sys->rl,rhs);
1516 if( rl != NULL ) {
1517 if( !(rl->solved) ) {
1518 FPRINTF(stderr,"ERROR: (linsol) linsol_var_value\n");
1519 FPRINTF(stderr," Rhs not solved yet.\n");
1520 return TRUE;
1521 }
1522 size = sys->capacity;
1523 varvalue = rl->varvalue;
1524 if( rl->transpose ) { /* ndx is an original row index */
1525 for (ndx = 0;ndx < size;ndx++) {
1526 vector[ndx] = varvalue[org_row_to_org_col(sys,ndx)];
1527 }
1528 }
1529 else{ /* ndx is an original column index */
1530 for (ndx = 0;ndx < size;ndx++) {
1531 vector[ndx] = varvalue[org_col_to_org_row(sys,ndx)];
1532 }
1533 }
1534 } else if( rhs != NULL ) {
1535 FPRINTF(stderr,"ERROR: (linsol) linsol_copy_solution\n");
1536 FPRINTF(stderr," Rhs does not exist.\n");
1537 return TRUE;
1538 }
1539 return FALSE;
1540 }
1541
1542 real64 linsol_eqn_residual(linsol_system_t sys, real64 *rhs, int32 ndx)
1543 {
1544 struct rhs_list *rl;
1545 mtx_coord_t nz;
1546 real64 value;
1547 real64 lhs;
1548
1549 check_system(sys);
1550 if( !sys->inverted ) {
1551 FPRINTF(stderr,"ERROR: (linsol) linsol_eqn_residual\n");
1552 FPRINTF(stderr," System not inverted yet.\n");
1553 }
1554 rl = find_rhs(sys->rl,rhs);
1555 if( rl != NULL ) {
1556 if( !(rl->solved) ) {
1557 FPRINTF(stderr,"ERROR: (linsol) linsol_eqn_residual\n");
1558 FPRINTF(stderr," Rhs not solved yet.\n");
1559 }
1560 if (rl->transpose) {
1561 /* ndx is an original column index */
1562 lhs = 0.0;
1563 nz.col = mtx_org_to_col(sys->coef,ndx);
1564 nz.row = mtx_FIRST;
1565 while( value = mtx_next_in_col(sys->coef,&nz,&(sys->reg.row)),
1566 nz.row != mtx_LAST )
1567 lhs += value * rl->varvalue
1568 [org_row_to_org_col(sys,mtx_row_to_org(sys->coef,nz.row))];
1569 return( lhs - rl->rhs[ndx] );
1570 } else {
1571 /* ndx is an original row index */
1572 lhs = 0.0;
1573 nz.row = mtx_org_to_row(sys->coef,ndx);
1574 nz.col = mtx_FIRST;
1575 while( value = mtx_next_in_row(sys->coef,&nz,&(sys->reg.col)),
1576 nz.col != mtx_LAST )
1577 lhs += value * rl->varvalue
1578 [org_col_to_org_row(sys,mtx_col_to_org(sys->coef,nz.col))];
1579 return( lhs - rl->rhs[ndx] );
1580 }
1581 } else if( rhs != NULL ) {
1582 FPRINTF(stderr,"ERROR: (linsol) linsol_eqn_residual\n");
1583 FPRINTF(stderr," Rhs does not exist.\n");
1584 }
1585 return 0.0; /* Function had no return value. Added this line. JDS */
1586 }

Properties

Name Value
svn:executable *

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