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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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