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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-chdr
File size: 24989 byte(s)
Setting up web subdirectory in repository
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