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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1012 - (show annotations) (download) (as text)
Wed Jan 3 03:13:20 2007 UTC (13 years, 6 months ago) by johnpye
File MIME type: text/x-chdr
File size: 23728 byte(s)
Added idalinear.[ch], which will expose the linsolqr solver as a linear solver for IDA.
Modified ida.c to support this new linear solver under the name 'ASCEND'.
1 /* 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 *
25 * Description: This module allows the user to create and manipulate
26 * matrices. The following is list of what constitutes
27 * a "matrix":
28 *
29 * - 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 *
33 * - Row and column permutations which keep track of
34 * where a given row/column "came from" originally.
35 *
36 * 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 *
46 * 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 *
50 * 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 * 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 ** header conventions: (DON'T SKIP READING THIS SECTION!!) *yawn*
287 ** -$- 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 * right. let's get into it...
303 */
304
305 #ifndef ASC_MTX_H
306 #define ASC_MTX_H
307
308 #include <stdio.h>
309 #include <string.h>
310 #include <utilities/ascConfig.h>
311
312 /***********************************************************************\
313 public mtx data structures
314 \***********************************************************************/
315 typedef struct mtx_header *mtx_matrix_t;
316 /**< Handle to the matrix */
317
318 /* Used to index rows and columns of matrices */
319
320 /**
321 *** Refers to (row,col) element of a matrix
322 **/
323 typedef struct mtx_coord_struct {
324 int32 row,col;
325 } mtx_coord_t;
326
327 /**
328 *** provide a list type for giving functions a list of coordinates.
329 **/
330 struct mtx_coord_list {
331 mtx_coord_t coord;
332 struct mtx_coord_list *next;
333 };
334
335 typedef struct mtx_range_struct {
336 int32 low,high;
337 } mtx_range_t;
338 /**<
339 *** Range refered to is low..high, inclusive. If low>high, range is empty
340 **/
341
342 typedef struct mtx_region_struct {
343 mtx_range_t row,col;
344 } mtx_region_t;
345 /**<
346 *** Rectangular region of a matrix
347 **/
348
349 typedef struct mtx_block_struct {
350 int32 nblocks;
351 mtx_region_t *block; /**< array of regions denoting blocks */
352 } mtx_block_t;
353 /**<
354 *** Block structure of a matrix.
355 **/
356
357 typedef struct mtx_sparse_vector_struct {
358 real64 *data;
359 int32 *idata;
360 int32 cap;
361 int32 len;
362 } mtx_sparse_t;
363 /**<
364 *** A structure for holding a compressed row or column vector.
365 *** The semantics imposed for this structure are very limited
366 *** so that one may be used (or reused) in many contexts.
367 *** An mtx_sparse_t doesn't know where its data came from, or by
368 *** what scheme the data is indexed.<br><br>
369 ***
370 *** The value of data[i] goes with mtx location somehow indexed idata[i].
371 ***
372 *** The values of idata[i] do not have to be in any particular
373 *** order, and the vector knows nothing about whether the indices
374 *** stored in idata are row indices, column indices, org_row indices
375 *** or org_column indices. This information is determined by the context
376 *** of the mtx_sparse_t.<br><br>
377 ***
378 *** cap is the size of data and idata.
379 *** len is the length of data and idata (starting from location 0)
380 *** which contains valid data.<br><br>
381 ***
382 *** A value of -1 for len indicates that pointers data and idata are invalid,
383 *** as does a value of 0 for cap.<br><br>
384 ***
385 *** NOTE: The mtx_sparse_t is not a pointer to a struct, it IS the struct.
386 **/
387
388 typedef struct mtx_block_perm_structure *mtx_block_perm_t;
389 /**<
390 *** The mtx_block_perm_structure contains a copy of the information
391 *** needed to apply a previously derived and saved ordering to the
392 *** matrix in question. Precisely what this information is is nobody's
393 *** business. We reserve the right to change at a moment's notice what
394 *** the mtx_block_perm_structure actually is.<br><br>
395 ***
396 *** The utility of the mtx_block_perm_structure is as follows:
397 *** Say in a sequence of nonlinear solves, one wishes to solve
398 *** many matrices of the same nonzero pattern, but different
399 *** nonzero values. One wants to reorder a specified block
400 *** of the matrix one time and reuse that initial ordering at each
401 *** factorization. Now suppose one also wants to use a linear
402 *** solution method that messes with the ordering during factorization.
403 *** One can do the following:
404 *** Create the mtx in question, output assign and partition it.<br><br>
405 *** <pre>
406 *** Call bp=mtx_create_block_perm(mtx).
407 *** This returns a copy of the information we need to store.
408 *** Working with individual blocks is then as follows:
409 *** 1- Reorder the block in question in the mtx in any way you like.
410 *** 2- Call mtx_update_block_perm(mtx,block_number,bp)
411 *** This updates the saved information for rows and columns
412 *** *** within *** the block specified by copying from the mtx.
413 *** All sorts of inexpensive sanity checks will be performed to
414 *** ensure that this is a wise thing to do.
415 *** 3- Factor, or otherwise do anything you want to the permutation
416 *** of the original mtx, so long as you do not mess with more than
417 *** one block. Changing the order (size) of the mtx is not permitted.
418 *** 4- Call mtx_restore_block_perm(mtx,block_number,bp)
419 *** This takes information stored in bp and resets the permutation
420 *** of the mtx within block_number to match it.
421 *** All sorts of inexpensive sanity checks will be performed to
422 *** ensure that this is a wise thing to do.
423 *** 5- Repeat steps 3 and 4 ad nauseum, or loop back through 1 and 2
424 *** if moving from block to block.
425 *** 6- When done with the matrix, or at least the process in question,
426 *** don't forget to do an mtx_destroy_block_perm(bp).
427 *** </pre>
428 *** Other variations on this scheme are possible. See the headers for
429 *** details of the individual functions.
430 **/
431
432 #define mtx_ENTIRE_MATRIX NULL
433 /**<
434 *** Refers to the entire matrix in any function where a region is called for.
435 **/
436
437 #define mtx_ALL_ROWS NULL
438 #define mtx_ALL_COLS NULL
439 /**<
440 *** Refers to all rows/columns in any function where a row/column range
441 *** is called for.
442 **/
443
444 #define mtx_CREATE_SPARSE NULL
445 /**<
446 *** In calls that take an mtx_sparse_t as an argument AND return a
447 *** pointer to a sparse AND do not note otherwise in the header
448 *** mtx_CREATE_SPARSE can be passed the argument and the call will
449 *** allocate the sparse_t data structure it returns. In such cases,
450 *** it is the users responsibility to see that the sparse is
451 *** eventually destroyed.
452 **/
453
454 #define mtx_IGNORE_ZEROES (1)
455 #define mtx_SOFT_ZEROES (0)
456 /**<
457 *** Those functions which fill a sparse_t might optionally return
458 *** soft zeroes given mtx_SOFT_ZEROES.
459 **/
460
461 #define mtx_FIRST ((int32)(-1))
462 #define mtx_LAST ((int32)(-1))
463 #define mtx_NONE ((int32)(-1))
464 #define mtx_ALL_BLOCKS ((int32)(-1))
465
466 #include "mtx_basic.h"
467 /**
468 *** All operations involving accounting/memory
469 *** management and CHANGES to matrix data structure.
470 *** This includes data file i/o and floating point
471 *** exception recovery handling.
472 **/
473
474 #include "mtx_perms.h"
475 #include "mtx_reorder.h"
476 /***
477 *** Functions for permuting matrices, reordering, etc.
478 *** These probably ought to be combined. mtx_reorder is currently
479 *** thinly disguised linsol code that needs to be redone in the mtx idiom.
480 **/
481
482 #include "mtx_query.h"
483 /***
484 *** READ_ONLY functions that don't make more sense
485 *** in mtx2_perms or mtx2_basic. Most vector data
486 *** and basic math operations fall in this file.
487 **/
488
489 #include "mtx_linal.h"
490 /***
491 *** Functions which probably shouldn't exist at all,
492 *** but are simply too damned expensive to do without
493 *** knowing all the internals of the mtx. These tend
494 *** to be ultracritical path functions in linsol.
495 *** These functions generally have very, very
496 *** complex semantics and preconditions and assumptions.
497 *** You really DON'T want to add functions to this
498 *** file or use those that are in it.
499 **/
500
501 #endif /* ASC_MTX_H */

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