/[ascend]/trunk/ascend4/solver/mtx.h
ViewVC logotype

Annotation of /trunk/ascend4/solver/mtx.h

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

Properties

Name Value
svn:executable *

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