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

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

Properties

Name Value
svn:executable *

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