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

Annotation of /trunk/ascend4/solver/linsol.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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