/[ascend]/trunk/base/generic/solver/mtx.h
ViewVC logotype

Annotation of /trunk/base/generic/solver/mtx.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1014 - (hide annotations) (download) (as text)
Wed Jan 3 04:49:17 2007 UTC (15 years, 7 months ago) by johnpye
File MIME type: text/x-chdr
File size: 23799 byte(s)
Fixing @addtogroup comments
1 johnpye 993 /* ASCEND modelling environment
2     Copyright (C) 1990 Karl Michael Westerberg
3     Copyright (C) 1993 Joseph Zaher
4     Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
5     Copyright (C) 1995 Benjamin Andrew Allan, Kirk Andre' Abbott
6     Copyright (C) 2006 Carnegie Mellon University
7    
8     This program is free software; you can redistribute it and/or modify
9     it under the terms of the GNU General Public License as published by
10     the Free Software Foundation; either version 2, or (at your option)
11     any later version.
12    
13     This program is distributed in the hope that it will be useful,
14     but WITHOUT ANY WARRANTY; without even the implied warranty of
15     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16     GNU General Public License for more details.
17    
18     You should have received a copy of the GNU General Public License
19     along with this program; if not, write to the Free Software
20     Foundation, Inc., 59 Temple Place - Suite 330,
21     Boston, MA 02111-1307, USA.
22     *//** @file
23     * mtx: Ascend Sparse Matrix Package.
24 aw0a 1 *
25 johnpye 993 * Description: This module allows the user to create and manipulate
26     * matrices. The following is list of what constitutes
27     * a "matrix":
28 aw0a 1 *
29 johnpye 993 * - A square nxn (where n is the order of the matrix)
30     * rectangular region of numbers, indexed by
31     * (row,col) pairs where 0<=row<n and 0<=col<n.
32 aw0a 1 *
33 johnpye 993 * - Row and column permutations which keep track of
34     * where a given row/column "came from" originally.
35 aw0a 1 *
36 johnpye 993 * The coefficient matrix (i.e. the nxn region of numbers)
37     * can be divided into two classes, non-zeros and zeros.
38     * Roughly speaking, a given element is a non-zero so long
39     * as its value has the POTENTIAL of not equaling zero.
40     * Thus, if the value of an element (r,c) is not equal to
41     * zero, then (r,c) must be classified as a non-zero:
42     * however the converse need not hold, unless either
43     * mtx_del_zr_in_row or mtx_del_zr_in_col has been
44     * recently called with r or c respectively.
45 jds 54 *
46 johnpye 993 * The mtx_matrix_t includes a data structure for block
47     * information storage, a block being a matrix subregion.
48     * This feature supports partitioning solvers.
49 jds 54 *
50 johnpye 993 * There are only 3 fundamental operations _on_ matrices
51     * vector scaling (mtx_mult_*)
52     * vector addition (mtx_add_*)
53     * permutation (mtx_swap_*)
54     * Matrix elements are maintained as
55     * relatively unordered linked lists. Mtx_next_in_* is
56     * most generally useful operator for doing sparse math
57     * with this matrix package. There are several operations
58     * that take matrices/vectors and return vectors/scalars.
59     * Sparse matrix-matrix computations are best coded with
60     * vector primitives and knowledge of the specific matrix,
61     * so they are not provided here.
62     *
63     * It's amazing for a package that only addresses the
64     * fundamentals, this still has over 100 calls and looks
65     * like an extremely well documented kitchen sink.
66     *
67     * The user may grep on extern to get a semiformatted list
68     * of the operator names in this file. Please follow the
69     * formatting if you add functions to this module.
70     * i.e. grep extern mtx2*h cooks up a reduced header.
71     *
72     ** vocabulary for this header file.
73     ** int32 a 32 bit signed integer (int usually, long on some archs.)
74     ** real64 a 64 bit real (double on most systems.)
75     ** By using these typedefs, all the code remains intact
76     ** across different systems/compilers: only the typedefs
77     ** have to change if there is a difference in data size.
78     **
79     ** mtx, matrix: a peculiar sparse matrix structure which should not
80     ** be manipulated except by means of functions in this
81     ** header. Attempts to go around by another path are
82     ** 99.44% likely to fail, and 95% likely to fail
83     ** catastrophically.
84     **
85     ** cur_row,
86     ** cur_col, or
87     ** row,col: an int32 denoting the ith(jth) row(column) as the
88     ** matrix is currently permuted (range 0..order-1).
89     ** org_row, org_col:
90     ** an int32 denoting the ith(jth) row(col) as the matrix
91     ** is ordered in the unpermuted state (range 0..order-1).
92     **
93     ** vec: an array of real64, generally of size
94     ** mtx->order or larger. Indexing scheme is indeterminate.
95     ** org_vec: a vec indexed by the org_row or org_col index.
96     ** cur_vec: a vec indexed by the (current) row or col index.
97     **
98     ** s,txxx: source or target xxx, where xxx may be a vec, row, col
99     ** org_row, or org_col.
100     **
101     ** order: the working size of a matrix.
102     ** order N ==> cols/rows 0->N-1 exist.
103     ** capacity: the length, essentially, of the mtx edges.
104     ** This number is >= order, normally ==.
105     **
106     ** coord: a pair of current indices (row,col).
107     ** range: a pair of current indices (low,high).
108     ** operations which take a range will do nothing if
109     ** low > high.
110     ** region: a pair of ranges (rowrange,colrange)
111     ** Coords, ranges and regions are passed by pointer to mtx functions.
112     **
113     ** sparse: a compressed vector, with indeterminate indexing scheme.
114     ** a sparse has a capacity (maximum amount of data)
115     ** a length (the current amount of data)
116     ** an array of real64 (the data) and
117     ** an array of int32 (where the data correspond to.)
118     ** The array of int32 may be cur or org and this
119     ** is determined by the context in which the sparse is used.
120     **
121     ** Sparses are generally passed by pointer.
122     *//*
123     'mtx' by Karl Michael Westerberg, created 5/3/90
124     'mtx2' by Benjamin Andrew Allan
125     Last in CVS: $Revision: 1.11 $ $Date: 1998/07/05 20:46:27 $ $Author: ballan $
126    
127     Dates: KW 06/90 - original version
128     JZ 01/94 - added output assignment and partitioning
129     operations allowing them to be performed
130     individually.
131     JZ 08/94 - reduced the internal practice of repeatedly
132     finding the existence of an element by re-
133     writing mtx_next_in_row and mtx_next_in_col
134     to use a static element pointer which can
135     advance and be rewound. Dramatic cpu time
136     savings were obtained by having the values
137     of subsequent elements returned by these
138 jds 54 * functions to eliminate the need for calling
139     * mtx_value on the obtained coordinate. A new
140     * function mtx_add_value was added to replace
141     * mtx_set_value in applications where it is
142     * known by the user that no nonzero currently
143     * exists at the coordinate whose value is being
144     * set.
145     * BA 08/94 - Protoized header file to reduce coding errors.
146     * BA 11/94 - Added mtx_write_region as it is done elsewhere
147     * and seems a common lowlevel operation.
148     * Made error output file specifiable in __debug.
149     * Made functions explicitly void if needed.
150     * For the following, see headers for explanation.
151     * Added mtx_copy_complete and mtx_copy_wo_incidence.
152     * Added mtx_drag from JZ.
153     * Added mtx_write_region_plot from JZ for xgraphs.
154     * Added mtx_free_reused_mem.
155     * Added mtx_add_row_sparse and mtx_add_col_sparse.
156     * Added mtx_nonzeros_in_region.
157     * Added mtx_numbers_in_row.
158     * Added mtx_numbers_in_col.
159     * Added mtx_numbers_in_region.
160     * Annotated next_in operator descriptions.
161     * Renamed mtx_add_value to mtx_fill_value in
162     * anticipation of a set of fill commands and to
163     * keep the mtx_add namespace implying math only.
164     * BA 12/94 - Added a set of dense vector (double array)
165     * handling routines:
166     * mtx_org_row_vec \
167     * mtx_org_col_vec | stuff an array which is cur or
168     * mtx_cur_row_vec | org indexed from a row/col.
169     * mtx_cur_col_vec /
170     * mtx_set_org_row_vec \ UNIMPLEMENTED
171     * mtx_set_org_col_vec | set a mtx row/col from
172     * mtx_set_cur_row_vec | cur or org indexed array.
173     * mtx_set_cur_col_vec /
174     * mtx_fill_org_row_vec \
175     * mtx_fill_org_col_vec | set empty mtx row/col from
176     * mtx_fill_cur_row_vec | cur or org indexed array.
177     * mtx_fill_cur_col_vec /
178     * mtx_zr_org_vec_using_row \
179     * mtx_zr_org_vec_using_col | zero org/cur vec entries
180     * mtx_zr_cur_vec_using_row | matching row/col incid.
181     * mtx_zr_cur_vec_using_col /
182     * mtx_org_vec_add_row \
183     * mtx_org_vec_add_col | vec += mtx row/col is a
184     * mtx_cur_vec_add_row | cur/org indexed array.
185     * mtx_cur_vec_add_col /
186     * mtx_add_row_org_vec \ UNIMPLEMENTED
187     * mtx_add_col_org_vec | mtx row/col += vec. vec a
188     * mtx_add_row_cur_vec | cur/org indexed array.
189     * mtx_add_col_cur_vec /
190     * Added mtx_add_row/col_series for use in
191     * common elimination schemes.
192     * Added dot operators for dot products
193     * mtx_row_dot_full_org_vec \
194     * mtx_col_dot_full_org_vec | row/col dot vec indexed
195     * mtx_row_dot_full_cur_vec | by org or cur row or col
196     * mtx_col_dot_full_cur_vec /
197     *
198     * Declared a sparse vector type that can be mapped
199     * directly to compressed row or column storage
200     * called mtx_sparse_t.
201     *
202     * Added mtx_write_region_human.
203     * Added mtx_read_region which is inverse of
204     * mtx_write_region.
205     * Revised mtx_row_max, mtx_col_max to return
206     * signed value of element with max abs(value).
207     * Consolidated mtx_copy* functions into one with
208     * a set of macros for the various flavors.
209     *
210     * 6/95 ba Added mtx_org_permute for 'natural' orderings.
211     * Added mtx_reverse_diagonal which does that.
212     *
213     * 8/95 ba Added slave matrices and attendant functions.
214     * mtx_create_slave adds a slave.
215     * mtx_chattel_size counts slave memory.
216     *
217     * Added map back functions with a drop tolerance.
218     * mtx_dropfill_cur_row_vec \ cur indexed vector
219     * mtx_dropfill_cur_col_vec / filled in mtx if big.
220     * Added a coordinate list struct for passing lists
221     * of elements around.
222     *
223     * 10/95 ba Added sparse vector fetch routines and destroy.
224     * mtx_org_row_sparse \
225     * mtx_org_col_sparse | stuff a sparse which is cur or
226     * mtx_cur_row_sparse | org indexed from a row/col.
227     * mtx_cur_col_sparse /
228     * mtx_destroy_sparse > deallocate a sparse.
229     * 11/95 ba Added spec, but not code, for mtx_clear_rowlist
230     * in anticipation of nonlinear solvers recomputing
231     * the gradient of only the nonlinear rows.
232     * Added also:
233     * mtx_assemble, which does just that in much
234     * the finite element sense.
235     * mtx_fill_org_value which takes org coord.
236     * mtx_exception_recover which resets internals
237     * after floating point err has left things in
238     * an illdefined state. Also in this process
239     * unified several pieces of code that manage
240     * internal buffers.
241     *
242     * 12/95 ba mtx_steal_org_row_sparse \ clear incidence and
243     * mtx_steal_org_col_sparse | stuff a sparse which is
244     * mtx_steal_cur_row_sparse | cur or org indexed from
245     * mtx_steal_cur_col_sparse / a row/col.
246     * mtx_steal_org_row_vec \ clear incidence and
247     * mtx_steal_org_col_vec | stuff an array which is
248     * mtx_steal_cur_row_vec | cur or org indexed from
249     * mtx_steal_cur_col_vec / a row/col.
250     * mtx_get_pivot_col \ LU partial pivot selection
251     * mtx_get_pivot_row / with a sparsity tolerance.
252     * 1/96 ba mtx_fill_org_row_sparse \
253     * mtx_fill_org_col_sparse | set empty row/col from
254     * mtx_fill_cur_row_sparse | cur or org indexed sparse.
255     * mtx_fill_cur_col_sparse /
256     * mtx_write_region_human \ Revised write human
257     * mtx_write_region_human_rows| so we can have row
258     * mtx_write_region_human_cols| or col oriented print.
259     * mtx_write_region_human_f / Supports old usage.
260     *
261     * 5/96 ba Split mtx.[ch] into several files to make life
262     * much easier for maintaining. All the mtx_*.c
263     * files are intended to stand together as a
264     * sparse matrix package. The split is for header
265     * digestibility and for ease of maintenance.
266     * The decomposition is as follows:
267     *
268     * mtx.h.
269     *
270     * This file, which includes the rest.
271     * This describes general mtx concepts.
272     * There is no mtx.c.
273     *
274     * mtx_basic.c mtx_basic.h Most memory management.
275     * mtx_perms.c mtx_perms.h Permuation management.
276     * mtx_query.c mtx_query.h Most queries, readonly math.
277     * mtx_linal.c mtx_linal.h Unjustifiable pieces.
278     * mtx_internal_use_only.[ch] Nobody's business.
279     *
280     * 4/97 ba Added mtx_transpose, mtx_isa_transpose to perms.
281     * Added mtx_duplicate_region to make a slave copy.
282     *
283     * 7/98 ba Added mtx_write_region_matlab for harwell format.
284     *
285     *
286 johnpye 993 ** header conventions: (DON'T SKIP READING THIS SECTION!!) *yawn*
287 aw0a 1 ** -$- in place of *** implies as described in compilation flags below.
288     ** ! ! in place of *** implies a warning which should be read before
289     ** ! ! using the function under discussion in certain ways.
290     **
291     ** compilation flags:
292     ** MTX_DEBUG - if TRUE, check_matrix will be invoked all over the place
293     ** and -$- returns will be in effect.
294     ** check_sparse is also used where appropriate if TRUE.
295     ** - if FALSE, routines marked with -$- will assume that you
296     ** are sending in a good mtx_matrix_t. Check accordingly.
297     ** They will also assume any mtx_sparse_t * is good.
298     ** MTX_DEBUG is wired in file mtx.c, but shouldn't be.
299     ** MTX_DEBUG TRUE slows this code down incredibly and should not be
300     ** used when compiling production code.
301     **
302 johnpye 993 * right. let's get into it...
303     */
304 aw0a 1
305 johnpye 67 #ifndef ASC_MTX_H
306     #define ASC_MTX_H
307 jds 54
308 johnpye 1014 /** @defgroup linear
309 johnpye 1013 Ascend linear solver routines
310     @{
311     */
312    
313 johnpye 1012 #include <stdio.h>
314     #include <string.h>
315     #include <utilities/ascConfig.h>
316    
317 aw0a 1 /***********************************************************************\
318     public mtx data structures
319     \***********************************************************************/
320     typedef struct mtx_header *mtx_matrix_t;
321 ben.allan 33 /**< Handle to the matrix */
322 aw0a 1
323 jds 54 /* Used to index rows and columns of matrices */
324 aw0a 1
325 jds 54 /**
326     *** Refers to (row,col) element of a matrix
327     **/
328 aw0a 1 typedef struct mtx_coord_struct {
329     int32 row,col;
330     } mtx_coord_t;
331 jds 54
332     /**
333     *** provide a list type for giving functions a list of coordinates.
334 aw0a 1 **/
335     struct mtx_coord_list {
336     mtx_coord_t coord;
337     struct mtx_coord_list *next;
338     };
339    
340     typedef struct mtx_range_struct {
341     int32 low,high;
342     } mtx_range_t;
343 jds 54 /**<
344 aw0a 1 *** Range refered to is low..high, inclusive. If low>high, range is empty
345     **/
346    
347     typedef struct mtx_region_struct {
348     mtx_range_t row,col;
349     } mtx_region_t;
350 ben.allan 33 /**<
351 aw0a 1 *** Rectangular region of a matrix
352     **/
353    
354     typedef struct mtx_block_struct {
355     int32 nblocks;
356 ben.allan 33 mtx_region_t *block; /**< array of regions denoting blocks */
357 aw0a 1 } mtx_block_t;
358 ben.allan 33 /**<
359 aw0a 1 *** Block structure of a matrix.
360     **/
361    
362     typedef struct mtx_sparse_vector_struct {
363     real64 *data;
364     int32 *idata;
365     int32 cap;
366     int32 len;
367     } mtx_sparse_t;
368 ben.allan 33 /**<
369 aw0a 1 *** A structure for holding a compressed row or column vector.
370     *** The semantics imposed for this structure are very limited
371     *** so that one may be used (or reused) in many contexts.
372     *** An mtx_sparse_t doesn't know where its data came from, or by
373 jds 54 *** what scheme the data is indexed.<br><br>
374 aw0a 1 ***
375     *** The value of data[i] goes with mtx location somehow indexed idata[i].
376     ***
377     *** The values of idata[i] do not have to be in any particular
378     *** order, and the vector knows nothing about whether the indices
379     *** stored in idata are row indices, column indices, org_row indices
380     *** or org_column indices. This information is determined by the context
381 jds 54 *** of the mtx_sparse_t.<br><br>
382 aw0a 1 ***
383     *** cap is the size of data and idata.
384     *** len is the length of data and idata (starting from location 0)
385 jds 54 *** which contains valid data.<br><br>
386 aw0a 1 ***
387     *** A value of -1 for len indicates that pointers data and idata are invalid,
388 jds 54 *** as does a value of 0 for cap.<br><br>
389 aw0a 1 ***
390     *** NOTE: The mtx_sparse_t is not a pointer to a struct, it IS the struct.
391     **/
392    
393     typedef struct mtx_block_perm_structure *mtx_block_perm_t;
394 ben.allan 33 /**<
395 aw0a 1 *** The mtx_block_perm_structure contains a copy of the information
396     *** needed to apply a previously derived and saved ordering to the
397     *** matrix in question. Precisely what this information is is nobody's
398     *** business. We reserve the right to change at a moment's notice what
399 jds 54 *** the mtx_block_perm_structure actually is.<br><br>
400 aw0a 1 ***
401     *** The utility of the mtx_block_perm_structure is as follows:
402     *** Say in a sequence of nonlinear solves, one wishes to solve
403     *** many matrices of the same nonzero pattern, but different
404     *** nonzero values. One wants to reorder a specified block
405     *** of the matrix one time and reuse that initial ordering at each
406     *** factorization. Now suppose one also wants to use a linear
407     *** solution method that messes with the ordering during factorization.
408     *** One can do the following:
409 jds 54 *** Create the mtx in question, output assign and partition it.<br><br>
410     *** <pre>
411 aw0a 1 *** Call bp=mtx_create_block_perm(mtx).
412     *** This returns a copy of the information we need to store.
413     *** Working with individual blocks is then as follows:
414     *** 1- Reorder the block in question in the mtx in any way you like.
415     *** 2- Call mtx_update_block_perm(mtx,block_number,bp)
416     *** This updates the saved information for rows and columns
417     *** *** within *** the block specified by copying from the mtx.
418     *** All sorts of inexpensive sanity checks will be performed to
419     *** ensure that this is a wise thing to do.
420     *** 3- Factor, or otherwise do anything you want to the permutation
421     *** of the original mtx, so long as you do not mess with more than
422     *** one block. Changing the order (size) of the mtx is not permitted.
423     *** 4- Call mtx_restore_block_perm(mtx,block_number,bp)
424     *** This takes information stored in bp and resets the permutation
425 jds 54 *** of the mtx within block_number to match it.
426 aw0a 1 *** All sorts of inexpensive sanity checks will be performed to
427     *** ensure that this is a wise thing to do.
428     *** 5- Repeat steps 3 and 4 ad nauseum, or loop back through 1 and 2
429     *** if moving from block to block.
430     *** 6- When done with the matrix, or at least the process in question,
431 jds 54 *** don't forget to do an mtx_destroy_block_perm(bp).
432     *** </pre>
433 aw0a 1 *** Other variations on this scheme are possible. See the headers for
434     *** details of the individual functions.
435     **/
436    
437     #define mtx_ENTIRE_MATRIX NULL
438 ben.allan 33 /**<
439 jds 54 *** Refers to the entire matrix in any function where a region is called for.
440 aw0a 1 **/
441    
442     #define mtx_ALL_ROWS NULL
443     #define mtx_ALL_COLS NULL
444 ben.allan 33 /**<
445 aw0a 1 *** Refers to all rows/columns in any function where a row/column range
446 jds 54 *** is called for.
447 aw0a 1 **/
448    
449     #define mtx_CREATE_SPARSE NULL
450 ben.allan 33 /**<
451 aw0a 1 *** In calls that take an mtx_sparse_t as an argument AND return a
452     *** pointer to a sparse AND do not note otherwise in the header
453     *** mtx_CREATE_SPARSE can be passed the argument and the call will
454     *** allocate the sparse_t data structure it returns. In such cases,
455     *** it is the users responsibility to see that the sparse is
456     *** eventually destroyed.
457     **/
458    
459     #define mtx_IGNORE_ZEROES (1)
460     #define mtx_SOFT_ZEROES (0)
461 ben.allan 33 /**<
462 aw0a 1 *** Those functions which fill a sparse_t might optionally return
463     *** soft zeroes given mtx_SOFT_ZEROES.
464     **/
465    
466 jds 54 #define mtx_FIRST ((int32)(-1))
467     #define mtx_LAST ((int32)(-1))
468     #define mtx_NONE ((int32)(-1))
469     #define mtx_ALL_BLOCKS ((int32)(-1))
470 aw0a 1
471     #include "mtx_basic.h"
472 jds 54 /**
473 aw0a 1 *** All operations involving accounting/memory
474     *** management and CHANGES to matrix data structure.
475     *** This includes data file i/o and floating point
476     *** exception recovery handling.
477     **/
478    
479     #include "mtx_perms.h"
480     #include "mtx_reorder.h"
481     /***
482     *** Functions for permuting matrices, reordering, etc.
483     *** These probably ought to be combined. mtx_reorder is currently
484     *** thinly disguised linsol code that needs to be redone in the mtx idiom.
485     **/
486    
487     #include "mtx_query.h"
488 jds 54 /***
489 aw0a 1 *** READ_ONLY functions that don't make more sense
490     *** in mtx2_perms or mtx2_basic. Most vector data
491     *** and basic math operations fall in this file.
492     **/
493    
494     #include "mtx_linal.h"
495     /***
496     *** Functions which probably shouldn't exist at all,
497     *** but are simply too damned expensive to do without
498     *** knowing all the internals of the mtx. These tend
499     *** to be ultracritical path functions in linsol.
500     *** These functions generally have very, very
501     *** complex semantics and preconditions and assumptions.
502     *** You really DON'T want to add functions to this
503     *** file or use those that are in it.
504     **/
505    
506 johnpye 1013 /** @} */
507    
508 johnpye 67 #endif /* ASC_MTX_H */

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