/[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 708 - (show annotations) (download) (as text)
Tue Jun 27 07:34:31 2006 UTC (17 years, 10 months ago) by johnpye
File MIME type: text/x-csrc
File size: 49173 byte(s)
Replaced some references to ascmalloc with ASC_NEW_ARRAY
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 "mtx.h"
40 #include "linsol.h"
41 #include <general/mathmacros.h>
42
43
44 #define LINSOL_DEBUG FALSE
45
46 struct rhs_list {
47 real64 *rhs; /* Vector of rhs values */
48 real64 *varvalue; /* Solution of the linear system */
49 boolean solved; /* ? has the rhs been solved */
50 boolean transpose; /* ? transpose the coef matrix */
51 struct rhs_list *next;
52 };
53
54 struct linsol_header {
55 int integrity;
56 int32 capacity; /* Capacity of arrays */
57 mtx_matrix_t coef; /* Coefficient matrix */
58 struct rhs_list *rl; /* List of rhs vectors */
59 int rlength; /* Length of rhs list */
60 real64 pivot_zero; /* Smallest exceptable pivot */
61 real64 ptol; /* Pivot selection tolerance */
62 mtx_matrix_t inv; /* Inverse matrix containing UL factors */
63 boolean inverted; /* ? has the matrix been inverted */
64 mtx_range_t rng; /* Pivot range */
65 mtx_region_t reg; /* Bounding region */
66 int32 rank; /* Rank of the matrix */
67 real64 smallest_pivot; /* Smallest pivot accepted */
68 };
69
70
71 #define OK ((int)439828676)
72 #define DESTROYED ((int)839276847)
73 static void check_system(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(struct rhs_list *rl)
99 /**
100 *** Destroys rhs list.
101 **/
102 {
103 while( rl != NULL ) {
104 struct rhs_list *p;
105 p = rl;
106 rl = rl->next;
107 if( p->varvalue != NULL )
108 ascfree( (POINTER)(p->varvalue) );
109 ascfree( (POINTER)p );
110 }
111 }
112
113 static struct rhs_list *find_rhs(struct rhs_list *rl, real64 *rhs)
114 /**
115 *** Searches for rhs in rhs list, returning it or NULL if not found.
116 **/
117 {
118 for( ; rl != NULL ; rl = rl->next )
119 if( rl->rhs == rhs )
120 break;
121 return(rl);
122 }
123
124
125 linsol_system_t linsol_create()
126 {
127 linsol_system_t sys;
128
129 sys = (linsol_system_t)ascmalloc( sizeof(struct linsol_header) );
130 sys->integrity = OK;
131 sys->capacity = 0;
132 sys->coef = NULL;
133 sys->rl = NULL;
134 sys->rlength = 0;
135 sys->pivot_zero = 1e-12; /* default value */
136 sys->ptol = 0.1; /* default value */
137 sys->inv = NULL;
138 sys->inverted = FALSE;
139 sys->rng.low = sys->rng.high = -1;
140 sys->reg.row.low = sys->reg.row.high = -1;
141 sys->reg.col.low = sys->reg.col.high = -1;
142 sys->rank = -1;
143 sys->smallest_pivot = 0.0;
144 return(sys);
145 }
146
147 void linsol_destroy(linsol_system_t sys)
148 {
149 check_system(sys);
150 if( sys->inv != NULL ) {
151 mtx_destroy(sys->inv);
152 }
153 destroy_rhs_list(sys->rl);
154 sys->integrity = DESTROYED;
155 ascfree( (POINTER)sys );
156 }
157
158 void linsol_set_matrix(linsol_system_t sys, mtx_matrix_t mtx)
159 {
160 check_system(sys);
161 sys->coef = mtx;
162 }
163
164 mtx_matrix_t linsol_get_matrix(linsol_system_t sys)
165 {
166 check_system(sys);
167 return( sys->coef );
168 }
169
170 mtx_matrix_t linsol_get_inverse(linsol_system_t sys)
171 {
172 check_system(sys);
173 return (sys->inv);
174 }
175
176 void linsol_add_rhs(linsol_system_t sys, real64 *rhs, boolean transpose)
177 {
178 struct rhs_list *rl;
179 int32 capacity;
180
181 check_system(sys);
182 rl = find_rhs(sys->rl,rhs);
183 if( rl != NULL ) {
184 return;
185 } else if( rhs != NULL ) {
186 rl = (struct rhs_list *)ascmalloc( sizeof(struct rhs_list) );
187 rl->rhs = rhs;
188 rl->varvalue = NULL;
189 rl->solved = FALSE;
190 rl->transpose = transpose;
191 rl->next = sys->rl;
192 sys->rl = rl;
193 ++(sys->rlength);
194 if (sys->coef) {
195 capacity = mtx_capacity(sys->coef);
196 if (capacity) {
197 rl->varvalue =(real64 *)
198 ascmalloc(capacity*sizeof(real64));
199 }
200 }
201 }
202 }
203
204 void linsol_remove_rhs(linsol_system_t sys, real64 *rhs)
205 {
206 struct rhs_list **q;
207
208 check_system(sys);
209 for( q = &(sys->rl) ; *q != NULL ; q = &((*q)->next) )
210 if( (*q)->rhs == rhs )
211 break;
212 if( *q != NULL ) {
213 struct rhs_list *p;
214 p = *q;
215 *q = p->next;
216 if( p->varvalue != NULL )
217 ascfree( (POINTER)(p->varvalue) );
218 ascfree( (POINTER)p );
219 --(sys->rlength);
220 } else if( rhs != NULL ) {
221 FPRINTF(stderr,"ERROR: (linsol) linsol_remove_rhs\n");
222 FPRINTF(stderr," Rhs does not exist.\n");
223 }
224 }
225
226 int linsol_number_of_rhs(linsol_system_t sys)
227 {
228 check_system(sys);
229 return( sys->rlength );
230 }
231
232 real64 *linsol_get_rhs(linsol_system_t sys, int n)
233 {
234 struct rhs_list *rl;
235 int count;
236
237 check_system(sys);
238
239 count = sys->rlength - 1 - n;
240 if( count < 0 ) return(NULL);
241 for( rl = sys->rl ; count > 0 && rl != NULL ; rl = rl->next )
242 --count;
243 return( rl == NULL ? NULL : rl->rhs );
244 }
245
246 void linsol_matrix_was_changed(linsol_system_t sys)
247 {
248 check_system(sys);
249 sys->inverted = FALSE;
250 }
251
252 void linsol_rhs_was_changed(linsol_system_t sys,real64 *rhs)
253 {
254 struct rhs_list *rl;
255
256 check_system(sys);
257 rl = find_rhs(sys->rl,rhs);
258 if( rl != NULL ) {
259 rl->solved = FALSE;
260 } else if( rhs != NULL ) {
261 FPRINTF(stderr,"ERROR: (linsol) linsol_rhs_was_modified\n");
262 FPRINTF(stderr," Rhs does not exist.\n");
263 }
264 }
265
266 void linsol_set_pivot_zero(linsol_system_t sys,real64 pivot_zero)
267 {
268 check_system(sys);
269 if( pivot_zero < 0.0 ) {
270 FPRINTF(stderr,"ERROR: (linsol) linsol_set_pivot_zero\n");
271 FPRINTF(stderr," Invalid pivot zero of %Lf\n",pivot_zero);
272 } else {
273 sys->pivot_zero = pivot_zero;
274 linsol_matrix_was_changed(sys);
275 }
276 }
277
278 real64 linsol_pivot_zero(linsol_system_t sys)
279 {
280 check_system(sys);
281 return( sys->pivot_zero );
282 }
283
284 void linsol_set_pivot_tolerance(linsol_system_t sys,real64 ptol)
285 {
286 check_system(sys);
287 if( ptol < 0.0 || ptol >= 1.0 ) {
288 FPRINTF(stderr,"ERROR: (linsol) linsol_set_pivot_tolerance\n");
289 FPRINTF(stderr," Invalid pivot tolerance of %Lf\n",ptol);
290 } else {
291 sys->ptol = ptol;
292 linsol_matrix_was_changed(sys);
293 }
294 }
295
296 real64 linsol_pivot_tolerance(linsol_system_t sys)
297 {
298 check_system(sys);
299 return( sys->ptol );
300 }
301
302
303 struct reorder_list { /* List of rows/columns and their counts. */
304 int32 ndx;
305 int32 count;
306 struct reorder_list *next;
307 };
308
309 struct reorder_vars {
310 mtx_matrix_t mtx;
311 mtx_region_t reg; /* Current active region */
312 int32 colhigh; /* Original value of reg.col.high */
313 struct reorder_list *tlist; /* Array of (enough) list elements */
314 int32 *rowcount; /* rowcount[reg.row.low .. reg.row.high] */
315 };
316
317 static void adjust_row_count(struct reorder_vars *vars,int32 removed_col)
318 /**
319 *** Adjusts the row counts to account for the (to be) removed column.
320 **/
321 {
322 mtx_coord_t nz;
323 real64 value;
324 nz.col = removed_col;
325 nz.row = mtx_FIRST;
326 while( value = mtx_next_in_col(vars->mtx,&nz,&(vars->reg.row)),
327 nz.row != mtx_LAST )
328 --(vars->rowcount[nz.row]);
329 }
330
331 static void assign_row_and_col(struct reorder_vars *vars,int32 row,int32 col)
332 /**
333 *** Assigns the given row to the given column, moving the row and column
334 *** to the beginning of the active region and removing them (readjusting
335 *** the region). The row counts are NOT adjusted. If col == mtx_NONE,
336 *** then no column is assigned and the row is moved to the end of the
337 *** active block instead. Otherwise, it is assumed that the column
338 *** is active.
339 **/
340 {
341 if( col == mtx_NONE ) {
342 mtx_swap_rows(vars->mtx,row,vars->reg.row.high);
343 vars->rowcount[row] = vars->rowcount[vars->reg.row.high];
344 --(vars->reg.row.high);
345 } else {
346 mtx_swap_rows(vars->mtx,row,vars->reg.row.low);
347 vars->rowcount[row] = vars->rowcount[vars->reg.row.low];
348 mtx_swap_cols(vars->mtx,col,vars->reg.col.low);
349 ++(vars->reg.row.low);
350 ++(vars->reg.col.low);
351 }
352 }
353
354 static void push_column_on_stack(struct reorder_vars *vars,int32 col)
355 /**
356 *** Pushes the given column onto the stack. It is assumed that the
357 *** column is active. Row counts are adjusted.
358 **/
359 {
360 adjust_row_count(vars,col);
361 mtx_swap_cols(vars->mtx,col,vars->reg.col.high);
362 --(vars->reg.col.high);
363 }
364
365 static int32 pop_column_from_stack(struct reorder_vars *vars)
366 /**
367 *** Pops the column on the "top" of the stack off of the stack and
368 *** returns the column index, where it now lies in the active region.
369 *** If the stack is empty, mtx_NONE is returned. Row counts are NOT
370 *** adjusted (this together with a subsequent assignment of this column
371 *** ==> no row count adjustment necessary).
372 **/
373 {
374 if( vars->reg.col.high < vars->colhigh )
375 return(++(vars->reg.col.high));
376 else
377 return( mtx_NONE );
378 }
379
380 static void assign_null_rows(struct reorder_vars *vars)
381 /**
382 *** Assigns empty rows, moving them to the assigned region. It is
383 *** assumed that row counts are correct. Columns are assigned off the
384 *** stack.
385 **/
386 {
387 int32 row;
388
389 for( row = vars->reg.row.low ; row <= vars->reg.row.high ; ++row )
390 if( vars->rowcount[row] == 0 )
391 assign_row_and_col(vars , row , pop_column_from_stack(vars));
392 }
393
394 static void forward_triangularize(struct reorder_vars *vars)
395 /**
396 *** Forward triangularizes the region, assigning singleton rows with their
397 *** one and only incident column until there are no more. The row counts
398 *** must be correct, and they are updated.
399 **/
400 {
401 boolean change;
402
403 do {
404 mtx_coord_t nz;
405 real64 value;
406 change = FALSE;
407 for( nz.row = vars->reg.row.low ;
408 nz.row <= vars->reg.row.high && vars->rowcount[nz.row] != 1;
409 ++nz.row ) ;
410 if( nz.row <= vars->reg.row.high ) {
411 /* found singleton row */
412 nz.col = mtx_FIRST;
413 value = mtx_next_in_row(vars->mtx,&nz,&(vars->reg.col));
414 adjust_row_count(vars,nz.col);
415 assign_row_and_col(vars,nz.row,nz.col);
416 change = TRUE;
417 }
418 } while( change );
419 }
420
421 static int32 select_row(struct reorder_vars *vars)
422 /**
423 *** Selects a row and returns its index. It is assumed that there is a
424 *** row. Row counts must be correct. vars->tlist will be used.
425 **/
426 {
427 int32 min_row_count;
428 int32 max_col_count;
429 int32 row;
430 int32 i, nties = 0; /* # elements currently defined in vars->tlist */
431
432 /* Set to something > any possible value */
433 min_row_count = vars->reg.col.high-vars->reg.col.low+2;
434 for( row = vars->reg.row.low ; row <= vars->reg.row.high ; ++row )
435 if( vars->rowcount[row] <= min_row_count ) {
436 if( vars->rowcount[row] < min_row_count ) {
437 min_row_count = vars->rowcount[row];
438 nties = 0;
439 }
440 vars->tlist[nties++].ndx = row;
441 }
442 /**
443 *** vars->tlist[0..nties-1] is a list of row numbers which tie for
444 *** minimum row count.
445 **/
446
447 max_col_count = -1; /* < any possible value */
448 i = 0;
449 while( i < nties ) {
450 int32 sum;
451 mtx_coord_t nz;
452 real64 value;
453
454 sum = 0;
455 nz.row = vars->tlist[i].ndx;
456 nz.col = mtx_FIRST;
457 while( value = mtx_next_in_row(vars->mtx,&nz,&(vars->reg.col)),
458 nz.col != mtx_LAST )
459 sum += mtx_nonzeros_in_col(vars->mtx,nz.col,&(vars->reg.row));
460 if( sum > max_col_count ) {
461 max_col_count = sum;
462 row = nz.row;
463 }
464 i++;
465 }
466 /**
467 *** Now row contains the row with the minimum row count, which has the
468 *** greatest total column count of incident columns among all rows with
469 *** the same (minimum) row count. Select it.
470 **/
471 return(row);
472 }
473
474 static void reorder(struct reorder_vars *vars)
475 /**
476 *** Reorders the assigned matrix vars->mtx within the specified bounding
477 *** block region vars->reg. The region is split into 6 subregions during
478 *** reordering: the rows are divided in two, and the columns divided in
479 *** three. Initially everything is in the active subregion. Ultimately,
480 *** everything will be assigned.
481 ***
482 *** <-- assigned -->|<-- active-->|<-- on stack -->|
483 *** ----+----------------+-------------+----------------+
484 *** a | | | |
485 *** s | | | |
486 *** s | | | |
487 *** i | (SQUARE) | | |
488 *** g | | | |
489 *** n | | | |
490 *** e | | | |
491 *** d | | | |
492 *** ----+----------------+-------------+----------------+
493 *** a | | | |
494 *** c | | ACTIVE | |
495 *** t | | REGION | |
496 *** i | | | |
497 *** v | | | |
498 *** e | | | |
499 *** ----+----------------+-------------+----------------+
500 ***
501 *** The algorithm is roughly as follows:
502 *** (1) select a row (by some criterion).
503 *** (2) push columns incident on that row onto the stack in decreasing
504 *** order of their length.
505 *** (3) pop first column off the stack and assign it to the selected
506 *** row.
507 *** (4) forward-triangularize (assign singleton rows with their one
508 *** and only incident column, until no longer possible).
509 ***
510 *** (1)-(4) should be repeated until the active subregion becomes empty.
511 ***
512 *** Everything above was written as though the entire matrix is
513 *** involved. In reality, only the relevant square region is involved.
514 **/
515 {
516 int32 row, size;
517 int32 *rowcount_array_origin;
518
519 size = vars->reg.row.high - vars->reg.row.low + 1;
520 size = MAX(size,vars->reg.col.high - vars->reg.col.low + 1);
521 vars->tlist = size > 0 ? (struct reorder_list *)
522 ascmalloc( size*sizeof(struct reorder_list) ) : NULL;
523 rowcount_array_origin = size > 0 ? (int32 *)
524 ascmalloc( size*sizeof(int32) ) : NULL;
525 vars->rowcount = rowcount_array_origin != NULL ?
526 rowcount_array_origin - vars->reg.row.low : NULL;
527
528 vars->colhigh = vars->reg.col.high;
529 /* Establish row counts */
530 for( row = vars->reg.row.low ; row <= vars->reg.row.high ; ++row )
531 vars->rowcount[row] =
532 mtx_nonzeros_in_row(vars->mtx,row,&(vars->reg.col));
533
534 while(TRUE) {
535 struct reorder_list *head;
536 int32 nelts; /* # elements "allocated" from vars->tlist */
537 mtx_coord_t nz;
538 real64 value;
539
540 forward_triangularize(vars);
541 assign_null_rows(vars);
542 if( vars->reg.row.low>vars->reg.row.high ||
543 vars->reg.col.low>vars->reg.col.high ) {
544 /* Active region is now empty, done */
545 if( vars->tlist != NULL )
546 ascfree( vars->tlist );
547 if( rowcount_array_origin != NULL )
548 ascfree( rowcount_array_origin );
549 return;
550 }
551
552 head = NULL;
553 nelts = 0;
554 nz.row = select_row(vars);
555 nz.col = mtx_FIRST;
556 while( value = mtx_next_in_row(vars->mtx,&nz,&(vars->reg.col)),
557 nz.col != mtx_LAST ) {
558 struct reorder_list **q,*p;
559
560 p = &(vars->tlist[nelts++]);
561 p->ndx = mtx_col_to_org(vars->mtx,nz.col);
562 p->count = mtx_nonzeros_in_col(vars->mtx,nz.col,&(vars->reg.row));
563 for( q = &head; *q && (*q)->count>p->count; q = &((*q)->next) );
564 p->next = *q;
565 *q = p;
566 }
567 /**
568 *** We now have a list of columns which intersect the selected row.
569 *** The list is in order of decreasing column count.
570 **/
571
572 /* Push incident columns on stack */
573 for( ; head != NULL ; head = head->next )
574 push_column_on_stack(vars,mtx_org_to_col(vars->mtx,head->ndx));
575
576 /* Pop column off stack and assign to selected row */
577 assign_row_and_col(vars , nz.row , pop_column_from_stack(vars));
578 }
579 /* Not reached. */
580 }
581
582 static boolean nonempty_row(mtx_matrix_t mtx,int32 row)
583 /**
584 *** ? row not empty in mtx.
585 **/
586 {
587 mtx_coord_t nz;
588 real64 value;
589 value = mtx_next_in_row(mtx, mtx_coord(&nz,row,mtx_FIRST), mtx_ALL_COLS);
590 return( nz.col != mtx_LAST );
591 }
592
593 static boolean nonempty_col(mtx_matrix_t mtx,int32 col)
594 /**
595 *** ? column not empty in mtx.
596 **/
597 {
598 mtx_coord_t nz;
599 real64 value;
600 value = mtx_next_in_col(mtx, mtx_coord(&nz,mtx_FIRST,col), mtx_ALL_ROWS);
601 return( nz.row != mtx_LAST );
602 }
603
604 static void determine_pivot_range(linsol_system_t sys)
605 /**
606 *** Calculates sys->rng from sys->coef.
607 **/
608 {
609 sys->reg.row.low = sys->reg.row.high = -1;
610 sys->reg.col.low = sys->reg.col.high = -1;
611
612 for( sys->rng.high=mtx_order(sys->coef) ; --(sys->rng.high) >= 0 ; ) {
613 if( nonempty_row(sys->coef,sys->rng.high) &&
614 sys->reg.row.high < 0 )
615 sys->reg.row.high = sys->rng.high;
616 if( nonempty_col(sys->coef,sys->rng.high) &&
617 sys->reg.col.high < 0 )
618 sys->reg.col.high = sys->rng.high;
619 if( nonempty_row(sys->coef,sys->rng.high) &&
620 nonempty_col(sys->coef,sys->rng.high) )
621 break;
622 }
623
624 for( sys->rng.low=0 ; sys->rng.low <= sys->rng.high ; ++(sys->rng.low) ) {
625 if( nonempty_row(sys->coef,sys->rng.low) &&
626 sys->reg.row.low < 0 )
627 sys->reg.row.low = sys->rng.low;
628 if( nonempty_col(sys->coef,sys->rng.low) &&
629 sys->reg.col.low < 0 )
630 sys->reg.col.low = sys->rng.low;
631 if( nonempty_row(sys->coef,sys->rng.low) &&
632 nonempty_col(sys->coef,sys->rng.low) )
633 break;
634 }
635 }
636
637 void linsol_reorder(linsol_system_t sys,mtx_region_t *region)
638 /**
639 *** The region to reorder is first isolated by truncating the region
640 *** provided to the largest square region confined to the matrix diagonal.
641 *** It is presumed it will contain no empty rows or columns and will
642 *** provide the basis of candidate pivots when inverting.
643 **/
644 {
645 struct reorder_vars vars;
646 check_system(sys);
647 if( region == mtx_ENTIRE_MATRIX )
648 determine_pivot_range(sys);
649 else {
650 sys->reg.row.low = region->row.low;
651 sys->reg.row.high = region->row.high;
652 sys->reg.col.low = region->col.low;
653 sys->reg.col.high = region->col.high;
654 sys->rng.low = MAX(region->row.low,region->col.low);
655 sys->rng.high = MIN(region->row.high,region->col.high);
656 }
657 vars.mtx = sys->coef;
658 vars.reg.row.low = vars.reg.col.low = sys->rng.low;
659 vars.reg.row.high = vars.reg.col.high = sys->rng.high;
660 reorder(&vars);
661 }
662
663
664 static void drag(mtx_matrix_t mtx,int32 n1,int32 n2)
665 /**
666 *** Drags row n1 to n2, moving row n1 to n2 and shifting all rows in
667 *** between back toward n1. Ditto for columns. This function works
668 *** regardless of whether n1 < n2 or n2 < n1.
669 **/
670 {
671 while( n1 < n2 ) {
672 mtx_swap_rows(mtx,n1,n1+1);
673 mtx_swap_cols(mtx,n1,n1+1);
674 ++n1;
675 }
676
677 while( n1 > n2 ) {
678 mtx_swap_rows(mtx,n1,n1-1);
679 mtx_swap_cols(mtx,n1,n1-1);
680 --n1;
681 }
682 }
683
684 static void eliminate_row(mtx_matrix_t mtx,mtx_range_t *rng,int32 row,real64 *tmp)
685 /**
686 *** Eliminates the given row to the left of the diagonal element, assuming
687 *** valid pivots for all of the diagonal elements above it (the elements
688 *** above those diagonal elements, if any exist, are assumed to be U
689 *** elements and therefore ignored). The temporary array is used by this
690 *** function to do its job. tmp[k], where rng->low <= k <= rng->high must
691 *** be defined (allocated) but need not be initialized.
692 **/
693 {
694 mtx_coord_t nz;
695 real64 value;
696 mtx_range_t left_to_eliminate,high_cols;
697
698 high_cols.low = row;
699 high_cols.high = rng->high;
700 left_to_eliminate.low = rng->low;
701 left_to_eliminate.high = row - 1;
702
703 /* Move non-zeros left of pivot from matrix to full array */
704 mem_zero_byte_cast(tmp+rng->low,0.0,(row-rng->low)*sizeof(real64));
705 nz.row = row;
706 nz.col = mtx_FIRST;
707 while( value = mtx_next_in_row(mtx,&nz,&left_to_eliminate),
708 nz.col != mtx_LAST )
709 tmp[nz.col] = value;
710 mtx_clear_row(mtx,row,&left_to_eliminate);
711
712 /* eliminates nzs from pivot, one by one, filling tmp with multipliers */
713 for( nz.row = row-1 ; nz.row >= rng->low ; --(nz.row) ) {
714 if( tmp[nz.row] == 0.0 )
715 continue; /* Nothing to do for this row */
716
717 nz.col = nz.row;
718 tmp[nz.row] /= mtx_value(mtx,&nz);
719 /* tmp[nz.row] now equals multiplier */
720
721 /* Perform "mtx_add_row" for full array part of the row */
722 left_to_eliminate.high = nz.row - 1;
723 nz.col = mtx_FIRST;
724 while( value = mtx_next_in_row(mtx,&nz,&left_to_eliminate),
725 nz.col != mtx_LAST )
726 tmp[nz.col] -= tmp[nz.row] * value;
727
728 /* Perform "mtx_add_row" on remaining part of row */
729 mtx_add_row(mtx,nz.row,row,-tmp[nz.row],&high_cols);
730 }
731
732 nz.row = row;
733 for( nz.col = rng->low ; nz.col < row ; ++(nz.col) )
734 if( tmp[nz.col] != 0.0 ) mtx_add_value(mtx,&nz,tmp[nz.col]);
735 }
736
737 #ifdef THIS_IS_AN_UNUSED_FUNCTION
738 static int32 top_of_spike(linsol_system_t sys,int32 col)
739 /**
740 *** Determines the top row (row of lowest index) in a possible spike
741 *** above the diagonal element in the given column. If there is no spike,
742 *** then (row = ) col is returned.
743 **/
744 {
745 mtx_range_t above_diagonal;
746 mtx_coord_t nz;
747 real64 value;
748 int32 top_row;
749
750 above_diagonal.low = sys->rng.low;
751 above_diagonal.high = col-1;
752 top_row = nz.col = col;
753 nz.row = mtx_FIRST;
754 while( value = mtx_next_in_col(sys->inv,&nz,&above_diagonal),
755 nz.row != mtx_LAST )
756 if( nz.row < top_row )
757 top_row = nz.row;
758 return( top_row );
759 }
760 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
761
762 static boolean col_is_a_spike(linsol_system_t sys,int32 col)
763 /**
764 *** Determines if the col is a spike, characterized by having any
765 *** nonzeros above the diagonal.
766 **/
767 {
768 mtx_range_t above_diagonal;
769 mtx_coord_t nz;
770 real64 value;
771 int32 top_row;
772
773 above_diagonal.low = sys->rng.low;
774 above_diagonal.high = col-1;
775 top_row = nz.col = col;
776 nz.row = mtx_FIRST;
777 while( value = mtx_next_in_col(sys->inv,&nz,&above_diagonal),
778 nz.row != mtx_LAST )
779 if( nz.row < col ) return TRUE;
780 return( FALSE );
781 }
782
783 static void invert(linsol_system_t sys)
784 /**
785 *** This is the heart of the linear equation solver. This function
786 *** factorizes the matrix into a lower (L) and upper (U) triangular
787 *** matrix. sys->smallest_pivot and sys->rank are calculated. The
788 *** RANKI method is utilized. At the end of elimination, the matrix A
789 *** is factored into A = U L, where U and L are stored as follows:
790 ***
791 *** <----- r ----> <- n-r ->
792 *** +--------------+---------+
793 *** | | |
794 *** | U | |
795 *** | | |
796 *** | L | | r
797 *** | | |
798 *** +--------------+---------+
799 *** | | |
800 *** | | 0 | n-r
801 *** | | |
802 *** +--------------+---------+
803 ***
804 *** The rows and columns have been permuted so that all of the pivoted
805 *** original rows and columns are found in the first r rows and columns
806 *** of the region. The last n-r rows and columns are unpivoted. U has
807 *** 1's on its diagonal, whereas L's diagonal is stored in the matrix.
808 ***
809 *** Everything above was written as though the entire matrix is
810 *** involved. In reality, only the relevant square region is involved.
811 **/
812 {
813 mtx_coord_t nz;
814 int32 last_row;
815 mtx_range_t pivot_candidates;
816 real64 *tmp,*tmp_array_origin;
817 int32 length;
818
819 length = sys->rng.high - sys->rng.low + 1;
820 tmp_array_origin = length > 0 ?
821 ASC_NEW_ARRAY(real64,length) : NULL;
822 tmp = tmp_array_origin != NULL ?
823 tmp_array_origin - sys->rng.low : NULL;
824
825 sys->smallest_pivot = MAXDOUBLE;
826 last_row = pivot_candidates.high = sys->rng.high;
827 for( nz.row = sys->rng.low ; nz.row <= last_row ; ) {
828 real64 pivot;
829
830 pivot_candidates.low = nz.col = nz.row;
831 pivot = mtx_value(sys->inv,&nz);
832 pivot = fabs(pivot);
833 if( pivot > sys->pivot_zero &&
834 pivot > sys->ptol * mtx_row_max(sys->inv,&nz,&pivot_candidates,NULL)
835 && !col_is_a_spike(sys,nz.row) ) {
836 /* Good pivot and not a spike: continue with next row */
837 if( pivot < sys->smallest_pivot )
838 sys->smallest_pivot = pivot;
839 ++(nz.row);
840 continue;
841 }
842
843 /**
844 *** Row is a spike row or will
845 *** be when a necessary column
846 *** exchange occurs.
847 **/
848 eliminate_row(sys->inv,&(sys->rng),nz.row,tmp);
849 if( (pivot=mtx_row_max(sys->inv,&nz,&pivot_candidates,NULL)) <=
850 sys->pivot_zero ) {
851 /* Dependent row, drag to the end */
852 drag(sys->inv,nz.row,last_row);
853 --last_row;
854 } else {
855 /* Independent row: nz contains best pivot */
856 mtx_swap_cols(sys->inv,nz.row,nz.col);
857 /* Move pivot to diagonal */
858 drag( sys->inv , nz.row , sys->rng.low );
859 if( pivot < sys->smallest_pivot )
860 sys->smallest_pivot = pivot;
861 ++(nz.row);
862 }
863
864 }
865 if( tmp_array_origin != NULL )
866 ascfree( (POINTER)tmp_array_origin );
867
868 sys->rank = last_row - sys->rng.low + 1;
869 }
870
871 static real64 *raise_capacity(real64 *vec,int32 oldcap,int32 newcap)
872 /**
873 *** Raises the capacity of the array and returns a new array.
874 *** It is assumed that oldcap < newcap. vec is destroyed or
875 *** returned as appropriate.
876 **/
877 {
878 real64 *newvec=NULL;
879 /*
880 real64 *newvec;
881 newvec = newcap > 0 ?
882 ASC_NEW_ARRAY(real64,newcap ) : NULL;
883 if( vec != NULL ) {
884 mem_move_cast(vec,newvec,oldcap*sizeof(real64));
885 ascfree( (POINTER)vec );
886 }
887 return(newvec);
888 */
889 if (newcap < oldcap)
890 return vec;
891 if (vec!=NULL) /* don't call realloc on null with newcap 0 or it frees */
892 newvec = (real64 *)ascrealloc(vec,(newcap * sizeof(real64)));
893 else
894 newvec=newcap > 0 ?
895 ASC_NEW_ARRAY(real64,newcap ) : NULL;
896 return newvec;
897 }
898
899 static void insure_capacity(linsol_system_t sys)
900 /**
901 *** Insures that the capacity of all of the solution vectors
902 *** for each rhs is large enough.
903 **/
904 {
905 int32 req_cap = mtx_capacity(sys->coef);
906
907 if( req_cap > sys->capacity ) {
908 struct rhs_list *rl;
909
910 for( rl = sys->rl ; rl != NULL ; rl = rl->next )
911 rl->varvalue = raise_capacity(rl->varvalue,sys->capacity,req_cap);
912 sys->capacity = req_cap;
913 }
914 }
915
916 static void forward_substitute(linsol_system_t sys,real64 *arr,boolean transpose)
917 /**
918 *** Forward substitute. It is assumed that the L (or U) part of
919 *** sys->inv is computed. This function converts c to x in place. The
920 *** values are stored in arr indexed by original row number (or original
921 *** column number).
922 ***
923 *** transpose = FALSE: transpose = TRUE:
924 *** T
925 *** L x = c U x = c
926 ***
927 *** The following formula holds:
928 *** 0<=k<r ==> x(k) = [c(k) - L(k,0..k-1) dot x(0..k-1)] / L(k,k)
929 *** or
930 *** 0<=k<r ==> x(k) = [c(k) - U(0..k-1,k) dot x(0..k-1)] / U(k,k)
931 **/
932 {
933 mtx_range_t dot_rng;
934 mtx_coord_t nz;
935 real64 value;
936
937 dot_rng.low = sys->rng.low;
938 if (transpose) { /* arr is indexed by original column number */
939 for( nz.col=dot_rng.low; nz.col < dot_rng.low+sys->rank; ++(nz.col) ) {
940 real64 sum;
941 register int32 org_col;
942
943 dot_rng.high = nz.col - 1;
944 sum = 0.0;
945 nz.row = mtx_FIRST;
946 while( value = mtx_next_in_col(sys->inv,&nz,&dot_rng),
947 nz.row != mtx_LAST )
948 sum += value*arr[mtx_col_to_org(sys->inv,nz.row)];
949
950 nz.row = nz.col;
951 org_col = mtx_col_to_org(sys->inv,nz.col);
952 arr[org_col] = (arr[org_col] - sum) / 1.0;
953 }
954
955 } else { /* arr is indexed by original row number */
956 for( nz.row=dot_rng.low; nz.row < dot_rng.low+sys->rank; ++(nz.row) ) {
957 real64 sum;
958 register int32 org_row;
959
960 dot_rng.high = nz.row - 1;
961 sum = 0.0;
962 nz.col = mtx_FIRST;
963 while( value = mtx_next_in_row(sys->inv,&nz,&dot_rng),
964 nz.col != mtx_LAST )
965 sum += value*arr[mtx_row_to_org(sys->inv,nz.col)];
966
967 nz.col = nz.row;
968 org_row = mtx_row_to_org(sys->inv,nz.row);
969 arr[org_row] = (arr[org_row] - sum) / mtx_value(sys->inv,&nz);
970 }
971 }
972 }
973
974 static void backward_substitute(linsol_system_t sys,real64 *arr,boolean transpose)
975 /**
976 *** Backward substitute. It is assumed that the U (or L) part of
977 *** sys->inv is computed. This function converts rhs to c in place. The
978 *** values are stored in arr indexed by original row number (or original
979 *** column number).
980 ***
981 *** transpose = FALSE: transpose = TRUE:
982 *** T
983 *** U c = rhs L c = rhs
984 ***
985 *** The following formula holds:
986 *** 0<=k<r ==> c(k) = [rhs(k) - U(k,k+1..r-1) dot c(k+1..r-1)] / U(k,k)
987 *** or
988 *** 0<=k<r ==> c(k) = [rhs(k) - L(k+1..r-1,k) dot c(k+1..r-1)] / L(k,k)
989 **/
990 {
991 mtx_range_t dot_rng;
992 mtx_coord_t nz;
993 real64 value;
994
995 dot_rng.high = sys->rng.low + sys->rank - 1;
996 if (transpose) { /* arr is indexed by original column number */
997 for( nz.col = dot_rng.high ; nz.col >= sys->rng.low ; --(nz.col) ) {
998 real64 sum;
999 register int32 org_col;
1000
1001 dot_rng.low = nz.col + 1;
1002 sum = 0.0;
1003 nz.row = mtx_FIRST;
1004 while( value = mtx_next_in_col(sys->inv,&nz,&dot_rng),
1005 nz.row != mtx_LAST )
1006 sum += value*arr[mtx_col_to_org(sys->inv,nz.row)];
1007
1008 nz.row = nz.col;
1009 org_col = mtx_col_to_org(sys->inv,nz.col);
1010 arr[org_col] = (arr[org_col] - sum) / mtx_value(sys->inv,&nz);
1011 }
1012
1013 } else { /* arr is indexed by original row number */
1014 for( nz.row = dot_rng.high ; nz.row >= sys->rng.low ; --(nz.row) ) {
1015 real64 sum;
1016 register int32 org_row;
1017
1018 dot_rng.low = nz.row + 1;
1019 sum = 0.0;
1020 nz.col = mtx_FIRST;
1021 while( value = mtx_next_in_row(sys->inv,&nz,&dot_rng),
1022 nz.col != mtx_LAST )
1023 sum += value*arr[mtx_row_to_org(sys->inv,nz.col)];
1024
1025 nz.col = nz.row;
1026 org_row = mtx_row_to_org(sys->inv,nz.row);
1027 arr[org_row] = (arr[org_row] - sum) / 1.0;
1028 }
1029 }
1030 }
1031
1032
1033 static void calc_dependent_rows(linsol_system_t sys)
1034 {
1035 mtx_coord_t nz;
1036 real64 value;
1037 mtx_range_t colrange;
1038 real64 *lc;
1039
1040 if( (sys->reg.row.low == sys->rng.low) &&
1041 ( sys->reg.row.high == sys->rng.low+sys->rank-1 ) )
1042 return;
1043 if (sys->rank==0) return;
1044
1045 lc = sys->capacity > 0 ? (real64 *)
1046 ascmalloc(sys->capacity*sizeof(real64)) : NULL;
1047 colrange.low = sys->rng.low;
1048 colrange.high = colrange.low + sys->rank - 1;
1049
1050 nz.row = sys->reg.row.low;
1051 for( ; nz.row <= sys->reg.row.high; nz.row++ ) {
1052 if( nz.row == sys->rng.low ) {
1053 nz.row = sys->rng.low+sys->rank-1;
1054 continue;
1055 }
1056 mem_zero_byte_cast(lc,0.0,(sys->capacity)*sizeof(real64));
1057 nz.col = mtx_FIRST;
1058 while( value = mtx_next_in_row(sys->inv,&nz,&colrange),
1059 nz.col != mtx_LAST )
1060 lc[mtx_col_to_org(sys->inv,nz.col)] = value;
1061 if( nz.row < sys->rng.low+sys->rank || nz.row > sys->rng.high )
1062 backward_substitute(sys,lc,TRUE);
1063 forward_substitute(sys,lc,TRUE);
1064 mtx_clear_row(sys->inv,nz.row,&colrange);
1065 for( nz.col=colrange.low; nz.col <= colrange.high; nz.col++ ) {
1066 real64 value;
1067 value = lc[mtx_col_to_org(sys->inv,nz.col)];
1068 if( value != 0.0 ) mtx_add_value(sys->inv,&nz,value);
1069 }
1070 }
1071
1072 if( lc ) ascfree(lc);
1073 }
1074
1075
1076 static void calc_dependent_cols(linsol_system_t sys)
1077 {
1078 mtx_coord_t nz;
1079 real64 value;
1080 mtx_range_t rowrange;
1081 real64 *lc;
1082
1083 if( (sys->reg.col.low == sys->rng.low) &&
1084 ( sys->reg.col.high == sys->rng.low+sys->rank-1 ) )
1085 return;
1086 if (sys->rank==0) return;
1087
1088 lc = sys->capacity > 0 ? (real64 *)
1089 ascmalloc(sys->capacity*sizeof(real64)) : NULL;
1090 rowrange.low = sys->rng.low;
1091 rowrange.high = rowrange.low + sys->rank - 1;
1092
1093 nz.col = sys->reg.col.low;
1094 for( ; nz.col <= sys->reg.col.high; nz.col++ ) {
1095 if( nz.col == sys->rng.low ) {
1096 nz.col = sys->rng.low+sys->rank-1;
1097 continue;
1098 }
1099 mem_zero_byte_cast(lc,0.0,sys->capacity*sizeof(real64));
1100 nz.row = mtx_FIRST;
1101 while( value = mtx_next_in_col(sys->inv,&nz,&rowrange),
1102 nz.row != mtx_LAST )
1103 lc[mtx_row_to_org(sys->inv,nz.row)] = value;
1104 if( nz.col < sys->rng.low+sys->rank || nz.col > sys->rng.high )
1105 backward_substitute(sys,lc,FALSE);
1106 forward_substitute(sys,lc,FALSE);
1107 mtx_clear_col(sys->inv,nz.col,&rowrange);
1108 for( nz.row=rowrange.low; nz.row <= rowrange.high; nz.row++ ) {
1109 real64 value;
1110 value = lc[mtx_row_to_org(sys->inv,nz.row)];
1111 if( value != 0.0 ) mtx_add_value(sys->inv,&nz,value);
1112 }
1113 }
1114
1115 if( lc ) ascfree(lc);
1116 }
1117
1118 #ifdef THIS_IS_AN_UNUSED_FUNCTION
1119 static void debug_out_inverse(FILE *fp,linsol_system_t sys)
1120 /**
1121 *** Outputs permutation and values of the nonzero elements in the
1122 *** inverse matrix.
1123 **/
1124 {
1125 mtx_coord_t nz;
1126 real64 value;
1127
1128 nz.row = sys->rng.low;
1129 for( ; nz.row <= sys->rng.high; ++(nz.row) ) {
1130 FPRINTF(fp," Row %d (org %d)\n",
1131 nz.row,
1132 mtx_row_to_org(sys->inv,nz.row));
1133 nz.col = mtx_FIRST;
1134 while( value = mtx_next_in_row(sys->inv,&nz,&(sys->rng)),
1135 nz.col != mtx_LAST )
1136 FPRINTF(fp," Col %d (org %d) has value %g\n",
1137 nz.col, mtx_col_to_org(sys->inv,nz.col), value);
1138 }
1139 }
1140 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
1141
1142 #ifdef THIS_IS_AN_UNUSED_FUNCTION
1143 /* from joe equivalent to ranki_jz*/
1144 static void invert_alt(linsol_system_t sys)
1145 /**
1146 *** This is the heart of the linear equation solver. This function
1147 *** factorizes the matrix into a lower (L) and upper (U) triangular
1148 *** matrix. sys->smallest_pivot and sys->rank are calculated. The
1149 *** RANKI method is utilized. At the end of elimination, the matrix A
1150 *** is factored into A = U L, where U and L are stored as follows:
1151 ***
1152 *** <----- r ----> <- n-r ->
1153 *** +--------------+---------+
1154 *** | | |
1155 *** | U | |
1156 *** | | |
1157 *** | L | | r
1158 *** | | |
1159 *** +--------------+---------+
1160 *** | | |
1161 *** | | 0 | n-r
1162 *** | | |
1163 *** +--------------+---------+
1164 ***
1165 *** The rows and columns have been permuted so that all of the pivoted
1166 *** original rows and columns are found in the first r rows and columns
1167 *** of the region. The last n-r rows and columns are unpivoted. U has
1168 *** 1's on its diagonal, whereas L's diagonal is stored in the matrix.
1169 ***
1170 *** Everything above was written as though the entire matrix is
1171 *** involved. In reality, only the relevant square region is involved.
1172 **/
1173 {
1174 real64 pivot;
1175 real64 biggest;
1176 mtx_coord_t nz, best;
1177 mtx_region_t candidates;
1178 real64 *tmp;
1179 int32 length;
1180
1181 length = sys->rng.high - sys->rng.low + 1;
1182 tmp = length > 0 ?
1183 ASC_NEW_ARRAY(real64,length) : NULL;
1184
1185 sys->smallest_pivot = MAXDOUBLE;
1186 candidates.row.high = sys->rng.high;
1187 candidates.col.high = sys->rng.high;
1188 for( nz.row = sys->rng.low ; nz.row <= candidates.row.high ; ) {
1189 nz.col = nz.row;
1190 pivot = mtx_value(sys->inv,&nz);
1191 pivot = fabs(pivot);
1192 candidates.row.low = nz.row;
1193 candidates.col.low = nz.row;
1194
1195 if( !col_is_a_spike(sys,nz.row) ) {
1196 best.col = nz.row;
1197 biggest = mtx_col_max(sys->inv,&best,&(candidates.row),NULL);
1198 if( biggest >= sys->pivot_zero ) {
1199 if( pivot < sys->pivot_zero || pivot < sys->ptol*biggest ) {
1200 mtx_swap_rows(sys->inv,nz.row,best.row);
1201 pivot = biggest;
1202 }
1203 if( pivot < sys->smallest_pivot ) sys->smallest_pivot = pivot;
1204 ++(nz.row);
1205 continue;
1206 }
1207 }
1208
1209 /**
1210 *** Row is a spike row or will
1211 *** be when a necessary column
1212 *** exchange occurs.
1213 **/
1214 eliminate_row(sys->inv,&(sys->rng),nz.row,tmp - sys->rng.low);
1215 pivot = mtx_row_max(sys->inv,&nz,&(candidates.col),NULL);
1216 if( pivot < sys->pivot_zero ) { /* pivot is an epsilon */
1217 /* Dependent row, drag nz to lower right */
1218 mtx_drag(sys->inv,nz.row,candidates.row.high);
1219 --(candidates.row.high);
1220 } else {
1221 /* Independent row, drag nz to upper left */
1222 mtx_swap_cols(sys->inv,nz.row,nz.col);
1223 mtx_drag(sys->inv,nz.row,sys->rng.low);
1224 if( pivot < sys->smallest_pivot )
1225 sys->smallest_pivot = pivot;
1226 ++(nz.row);
1227 }
1228 }
1229 if( tmp != NULL ) ascfree( (POINTER)tmp );
1230 sys->rank = candidates.row.high - sys->rng.low + 1;
1231 }
1232 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
1233
1234 #define KAA_DEBUG 0
1235 void linsol_invert(linsol_system_t sys,mtx_region_t *region)
1236 /**
1237 *** The region to invert is first isolated by truncating the region
1238 *** provided to the largest square region confined to the matrix diagonal.
1239 *** It is presumed it will contain no empty rows or columns and that it has
1240 *** been previously reordered using linsol_reorder.
1241 **/
1242 {
1243 struct rhs_list *rl;
1244 #if KAA_DEBUG
1245 double time;
1246 #endif
1247
1248 check_system(sys);
1249 if( sys->inverted )
1250 return;
1251
1252 #if KAA_DEBUG
1253 time = tm_cpu_time();
1254 #endif
1255 if( sys->inv != NULL )
1256 mtx_destroy(sys->inv);
1257 sys->inv = mtx_copy(sys->coef);
1258 sys->rank = -1;
1259 sys->smallest_pivot = 0.0;
1260 for( rl = sys->rl ; rl != NULL ; rl = rl->next )
1261 rl->solved = FALSE;
1262 insure_capacity(sys);
1263 if( region == mtx_ENTIRE_MATRIX )
1264 determine_pivot_range(sys);
1265 else {
1266 sys->reg = *region;
1267 /*
1268 sys->reg.row.low = region->row.low;
1269 sys->reg.row.high = region->row.high;
1270 sys->reg.col.low = region->col.low;
1271 sys->reg.col.high = region->col.high;
1272 */
1273 sys->rng.low = MAX(region->row.low,region->col.low);
1274 sys->rng.high = MIN(region->row.high,region->col.high);
1275 }
1276 invert(sys);
1277 calc_dependent_rows(sys);
1278 calc_dependent_cols(sys);
1279 sys->inverted = TRUE;
1280 #if KAA_DEBUG
1281 time = tm_cpu_time() - time;
1282 FPRINTF(stderr,"Time for Inversion = %f\n",time);
1283 #endif /* KAA_DEBUG */
1284 }
1285
1286 int32 linsol_rank(linsol_system_t sys)
1287 {
1288 check_system(sys);
1289 if( !sys->inverted ) {
1290 FPRINTF(stderr,"ERROR: (linsol) linsol_rank\n");
1291 FPRINTF(stderr," System not inverted yet.\n");
1292 }
1293 return(sys->rank);
1294 }
1295
1296 real64 linsol_smallest_pivot(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 }

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