generic/solver/mtx.h

/*
 *  mtx: Ascend Sparse Matrix Package
 *  by Benjamin Andrew Allan
 *  Derived from mtx by Karl Michael Westerberg
 *  Created: 5/3/90
 *  Version: $Revision: 1.11 $
 *  Version control file: $RCSfile: mtx.h,v $
 *  Date last modified: $Date: 1998/07/05 20:46:27 $
 *  Last modified by: $Author: ballan $
 *
 *  This file is part of the SLV solver.
 *
 *  Copyright (C) 1990 Karl Michael Westerberg
 *  Copyright (C) 1993 Joseph Zaher
 *  Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
 *  Copyright (C) 1995 Benjamin Andrew Allan, Kirk Andre' Abbott
 *
 *  The SLV solver is free software; you can redistribute
 *  it and/or modify it under the terms of the GNU General Public License as
 *  published by the Free Software Foundation; either version 2 of the
 *  License, or (at your option) any later version.
 *                                         
 *  The SLV solver is distributed in hope that it will be
 *  useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License along with
 *  the program; if not, write to the Free Software Foundation, Inc., 675
 *  Mass Ave, Cambridge, MA 02139 USA.  Check the file named COPYING.
 *  COPYING is found in ../compiler.
 */

/** @file
 *  mtx: Ascend Sparse Matrix Package.
 *  <pre>
 *  Contents:     Matrix modules module
 *
 *  Authors:      Karl Westerberg
 *                Joseph Zaher
 *                Ben Allan
 *
 *  Dates:    KW  06/90 - original version
 *            JZ  01/94 - added output assignment and partitioning
 *                        operations allowing them to be performed
 *                        individually.
 *            JZ  08/94 - reduced the internal practice of repeatedly
 *                        finding the existence of an element by re-
 *                        writing mtx_next_in_row and mtx_next_in_col
 *                        to use a static element pointer which can
 *                        advance and be rewound.  Dramatic cpu time
 *                        savings were obtained by having the values
 *                        of subsequent elements returned by these
 *                        functions to eliminate the need for calling
 *                        mtx_value on the obtained coordinate.  A new
 *                        function mtx_add_value was added to replace
 *                        mtx_set_value in applications where it is
 *                        known by the user that no nonzero currently
 *                        exists at the coordinate whose value is being
 *                        set.
 *            BA  08/94 - Protoized header file to reduce coding errors.
 *            BA  11/94 - Added mtx_write_region as it is done elsewhere
 *                        and seems a common lowlevel operation.
 *                        Made error output file specifiable in __debug.
 *                        Made functions explicitly void if needed.
 *                        For the following, see headers for explanation.
 *                        Added mtx_copy_complete and mtx_copy_wo_incidence.
 *                        Added mtx_drag from JZ.
 *                        Added mtx_write_region_plot from JZ for xgraphs.
 *                        Added mtx_free_reused_mem.
 *                        Added mtx_add_row_sparse and mtx_add_col_sparse.
 *                        Added mtx_nonzeros_in_region.
 *                        Added mtx_numbers_in_row.
 *                        Added mtx_numbers_in_col.
 *                        Added mtx_numbers_in_region.
 *                        Annotated next_in operator descriptions.
 *                        Renamed mtx_add_value to mtx_fill_value in
 *                        anticipation of a set of fill commands and to
 *                        keep the mtx_add namespace implying math only.
 *            BA  12/94 - Added a set of dense vector (double array)
 *                        handling routines:
 *                        mtx_org_row_vec \
 *                        mtx_org_col_vec | stuff an array which is cur or
 *                        mtx_cur_row_vec | org indexed from a row/col.
 *                        mtx_cur_col_vec /
 *                        mtx_set_org_row_vec \ UNIMPLEMENTED
 *                        mtx_set_org_col_vec | set a mtx row/col from
 *                        mtx_set_cur_row_vec | cur or org indexed array.
 *                        mtx_set_cur_col_vec /
 *                        mtx_fill_org_row_vec \
 *                        mtx_fill_org_col_vec | set empty mtx row/col from
 *                        mtx_fill_cur_row_vec | cur or org indexed array.
 *                        mtx_fill_cur_col_vec /
 *                        mtx_zr_org_vec_using_row \
 *                        mtx_zr_org_vec_using_col | zero org/cur vec entries
 *                        mtx_zr_cur_vec_using_row | matching row/col incid.
 *                        mtx_zr_cur_vec_using_col /
 *                        mtx_org_vec_add_row \
 *                        mtx_org_vec_add_col | vec += mtx row/col is a
 *                        mtx_cur_vec_add_row | cur/org indexed array.
 *                        mtx_cur_vec_add_col /
 *                        mtx_add_row_org_vec \ UNIMPLEMENTED
 *                        mtx_add_col_org_vec | mtx row/col += vec. vec a
 *                        mtx_add_row_cur_vec | cur/org indexed array.
 *                        mtx_add_col_cur_vec /
 *                        Added mtx_add_row/col_series for use in
 *                        common elimination schemes.
 *                        Added dot operators for dot products
 *                        mtx_row_dot_full_org_vec \
 *                        mtx_col_dot_full_org_vec | row/col dot vec indexed
 *                        mtx_row_dot_full_cur_vec | by org or cur row or col
 *                        mtx_col_dot_full_cur_vec /
 *
 *                        Declared a sparse vector type that can be mapped
 *                        directly to compressed row or column storage
 *                        called mtx_sparse_t.
 *
 *                        Added mtx_write_region_human.
 *                        Added mtx_read_region which is inverse of
 *                        mtx_write_region.
 *                        Revised mtx_row_max, mtx_col_max to return
 *                        signed value of element with max abs(value).
 *                        Consolidated mtx_copy* functions into one with
 *                        a set of macros for the various flavors.
 *
 *               6/95 ba  Added mtx_org_permute for 'natural' orderings.
 *                        Added mtx_reverse_diagonal which does that.
 *
 *               8/95 ba  Added slave matrices and attendant functions.
 *                        mtx_create_slave       adds a slave.
 *                        mtx_chattel_size       counts slave memory.
 *
 *                        Added map back functions with a drop tolerance.
 *                        mtx_dropfill_cur_row_vec \ cur indexed vector
 *                        mtx_dropfill_cur_col_vec / filled in mtx if big.
 *                        Added a coordinate list struct for passing lists
 *                        of elements around.
 *
 *              10/95 ba  Added sparse vector fetch routines and destroy.
 *                        mtx_org_row_sparse \
 *                        mtx_org_col_sparse | stuff a sparse which is cur or
 *                        mtx_cur_row_sparse | org indexed from a row/col.
 *                        mtx_cur_col_sparse /
 *                        mtx_destroy_sparse >  deallocate a sparse.
 *              11/95 ba  Added spec, but not code, for mtx_clear_rowlist
 *                        in anticipation of nonlinear solvers recomputing
 *                        the gradient of only the nonlinear rows.
 *                        Added also:
 *                        mtx_assemble, which does just that in much
 *                                      the finite element sense.
 *                        mtx_fill_org_value which takes org coord.
 *                        mtx_exception_recover  which resets internals
 *                          after floating point err has left things in
 *                          an illdefined state. Also in this process
 *                          unified several pieces of code that manage
 *                          internal buffers.
 *
 *              12/95 ba  mtx_steal_org_row_sparse \ clear incidence and
 *                        mtx_steal_org_col_sparse | stuff a sparse which is
 *                        mtx_steal_cur_row_sparse | cur or org indexed from
 *                        mtx_steal_cur_col_sparse / a row/col.
 *                        mtx_steal_org_row_vec \ clear incidence and
 *                        mtx_steal_org_col_vec | stuff an array which is
 *                        mtx_steal_cur_row_vec | cur or org indexed from
 *                        mtx_steal_cur_col_vec / a row/col.
 *                        mtx_get_pivot_col \  LU partial pivot selection
 *                        mtx_get_pivot_row /  with a sparsity tolerance.
 *               1/96 ba  mtx_fill_org_row_sparse \
 *                        mtx_fill_org_col_sparse | set empty row/col from
 *                        mtx_fill_cur_row_sparse | cur or org indexed sparse.
 *                        mtx_fill_cur_col_sparse /
 *                        mtx_write_region_human     \ Revised write human
 *                        mtx_write_region_human_rows| so we can have row
 *                        mtx_write_region_human_cols| or col oriented print.
 *                        mtx_write_region_human_f   / Supports old usage.
 *
 *               5/96 ba  Split mtx.[ch] into several files to make life
 *                        much easier for maintaining. All the mtx_*.c
 *                        files are intended to stand together as a
 *                        sparse matrix package. The split is for header
 *                        digestibility and for ease of maintenance.
 *                        The decomposition is as follows:
 *
 *                        mtx.h.
 *
 *                        This file, which includes the rest.
 *                        This describes general mtx concepts.
 *                        There is no mtx.c.
 *
 *                        mtx_basic.c mtx_basic.h Most memory management.
 *                        mtx_perms.c mtx_perms.h Permuation management.
 *                        mtx_query.c mtx_query.h Most queries, readonly math.
 *                        mtx_linal.c mtx_linal.h Unjustifiable pieces.
 *                        mtx_internal_use_only.[ch] Nobody's business.
 *
 *               4/97 ba  Added mtx_transpose, mtx_isa_transpose to perms.
 *                        Added mtx_duplicate_region to make a slave copy.
 *
 *               7/98 ba  Added mtx_write_region_matlab for harwell format.
 *
 *
 *  Description:  This module allows the user to create and manipulate
 *                matrices.  The following is list of what constitutes
 *                a "matrix":
 *
 *                   - A square nxn (where n is the order of the matrix)
 *                     rectangular region of numbers, indexed by
 *                     (row,col) pairs where 0<=row<n and 0<=col<n.
 *
 *                   - Row and column permutations which keep track of
 *                     where a given row/column "came from" originally.
 *
 *                The coefficient matrix (i.e. the nxn region of numbers)
 *                can be divided into two classes, non-zeros and zeros.
 *                Roughly speaking, a given element is a non-zero so long
 *                as its value has the POTENTIAL of not equaling zero.
 *                Thus, if the value of an element (r,c) is not equal to
 *                zero, then (r,c) must be classified as a non-zero:
 *                however the converse need not hold, unless either
 *                mtx_del_zr_in_row or mtx_del_zr_in_col has been
 *                recently called with r or c respectively.
 *
 *                The mtx_matrix_t includes a data structure for block
 *                information storage, a block being a matrix subregion.
 *                This feature supports partitioning solvers.
 *
 *                There are only 3 fundamental operations _on_ matrices
 *                  vector scaling  (mtx_mult_*)
 *                  vector addition (mtx_add_*)
 *                  permutation     (mtx_swap_*)
 *                Matrix elements are maintained as
 *                relatively unordered linked lists. Mtx_next_in_* is
 *                most generally useful operator for doing sparse math
 *                with this matrix package. There are several operations
 *                that take matrices/vectors and return vectors/scalars.
 *                Sparse matrix-matrix computations are best coded with
 *                vector primitives and knowledge of the specific matrix,
 *                so they are not provided here.
 *
 *                It's amazing for a package that only addresses the
 *                fundamentals, this still has over 100 calls and looks
 *                like an extremely well documented kitchen sink.
 *
 *                The user may grep on extern to get a semiformatted list
 *                of the operator names in this file. Please follow the
 *                formatting if you add functions to this module.
 *                i.e.  grep extern mtx2*h  cooks up a reduced header.
 **
 ** header conventions: (DON'T SKIP READING THIS SECTION!!)
 **   -$- in place of *** implies as described in compilation flags below.
 **   ! ! in place of *** implies a warning which should be read before
 **   ! ! using the function under discussion in certain ways.
 **
 ** compilation flags:
 **   MTX_DEBUG  - if TRUE, check_matrix will be invoked all over the place
 **                and -$- returns will be in effect.
 **                check_sparse is also used where appropriate if TRUE.
 **              - if FALSE, routines marked with -$- will assume that you
 **                are sending in a good mtx_matrix_t. Check accordingly.
 **                They will also assume any mtx_sparse_t * is good.
 **   MTX_DEBUG is wired in file mtx.c, but shouldn't be.
 **   MTX_DEBUG TRUE slows this code down incredibly and should not be
 **   used when compiling production code.
 **
 ** vocabulary for this header file.
 **   int32   a 32 bit signed integer (int usually, long on some archs.)
 **   real64  a 64 bit real (double on most systems.)
 **                 By using these typedefs, all the code remains intact
 **                 across different systems/compilers: only the typedefs
 **                 have to change if there is a difference in data size.
 **
 **   mtx, matrix:  a peculiar sparse matrix structure which should not
 **                 be manipulated except by means of functions in this
 **                 header. Attempts to go around by another path are
 **                 99.44% likely to fail, and 95% likely to fail
 **                 catastrophically.
 **
 **   cur_row,
 **   cur_col, or
 **   row,col:      an int32 denoting the ith(jth) row(column) as the
 **                 matrix is currently permuted (range 0..order-1).
 **   org_row, org_col:
 **                 an int32 denoting the ith(jth) row(col) as the matrix
 **                 is ordered in the unpermuted state (range 0..order-1).
 **
 **   vec:          an array of real64, generally of size
 **                 mtx->order or larger. Indexing scheme is indeterminate.
 **   org_vec:      a vec indexed by the org_row or org_col index.
 **   cur_vec:      a vec indexed by the (current) row or col index.
 **
 **   s,txxx:       source or target xxx, where xxx may be a  vec, row, col
 **                 org_row, or org_col.
 **
 **   order:        the working size of a matrix.
 **                 order N ==> cols/rows 0->N-1 exist.
 **   capacity:     the length, essentially, of the mtx edges.
 **                 This number is >= order, normally ==.
 **
 **   coord:        a pair of current indices (row,col).
 **   range:        a pair of current indices (low,high).
 **                 operations which take a range will do nothing if
 **                 low > high.
 **   region:       a pair of ranges (rowrange,colrange)
 **   Coords, ranges and regions are passed by pointer to mtx functions.
 **
 **   sparse:       a compressed vector, with indeterminate indexing scheme.
 **                 a sparse has a capacity (maximum amount of data)
 **                 a length (the current amount of data)
 **                 an array of real64 (the data) and
 **                 an array of int32 (where the data correspond to.)
 **                 The array of int32 may be cur or org and this
 **                 is determined by the context in which the sparse is used.
 **
 **                 Sparses are generally passed by pointer.
 *
 * Requires:      #include <stdio.h>
 *                #include <string.h>
 *                #include "utilities/ascConfig.h"
 *  </pre>
 */

#ifndef __MTX_H_SEEN__
#define __MTX_H_SEEN__

/***********************************************************************\
  public mtx data structures
\***********************************************************************/
typedef struct mtx_header *mtx_matrix_t;
/**< Handle to the matrix */

/* Used to index rows and columns of matrices */

/**
 ***  Refers to (row,col) element of a matrix
 **/
typedef struct mtx_coord_struct {
   int32 row,col;
} mtx_coord_t;

/**
 *** provide a list type for giving functions a list of coordinates.
 **/
struct mtx_coord_list {
  mtx_coord_t coord;
  struct mtx_coord_list *next;
};

typedef struct mtx_range_struct {
  int32 low,high;
} mtx_range_t;
/**<
 ***  Range refered to is low..high, inclusive.  If low>high, range is empty
 **/

typedef struct mtx_region_struct {
  mtx_range_t row,col;
} mtx_region_t;
/**< 
 ***  Rectangular region of a matrix
 **/

typedef struct mtx_block_struct {
  int32 nblocks;
  mtx_region_t *block;      /**< array of regions denoting blocks */
} mtx_block_t;
/**< 
 *** Block structure of a matrix.
 **/

typedef struct mtx_sparse_vector_struct {
  real64 *data;
  int32 *idata;
  int32 cap;
  int32 len;
} mtx_sparse_t;
/**< 
 *** A structure for holding a compressed row or column vector.
 *** The semantics imposed for this structure are very limited
 *** so that one may be used (or reused) in many contexts.
 *** An mtx_sparse_t doesn't know where its data came from, or by
 *** what scheme the data is indexed.<br><br>
 ***
 *** The value of data[i] goes with mtx location somehow indexed idata[i].
 ***
 *** The values of idata[i] do not have to be in any particular
 *** order, and the vector knows nothing about whether the indices
 *** stored in idata are row indices, column indices, org_row indices
 *** or org_column indices. This information is determined by the context
 *** of the mtx_sparse_t.<br><br>
 ***
 *** cap is the size of data and idata.
 *** len is the length of data and idata (starting from location 0)
 *** which contains valid data.<br><br>
 ***
 *** A value of -1 for len indicates that pointers data and idata are invalid,
 *** as does a value of 0 for cap.<br><br>
 ***
 *** NOTE: The mtx_sparse_t is not a pointer to a struct, it IS the struct.
 **/

typedef struct mtx_block_perm_structure *mtx_block_perm_t;
/**< 
 ***  The mtx_block_perm_structure contains a copy of the information
 ***  needed to apply a previously derived and saved ordering to the
 ***  matrix in question. Precisely what this information is is nobody's
 ***  business. We reserve the right to change at a moment's notice what
 ***  the mtx_block_perm_structure actually is.<br><br>
 ***
 ***  The utility of the mtx_block_perm_structure is as follows:
 ***  Say in a sequence of nonlinear solves, one wishes to solve
 ***  many matrices of the same nonzero pattern, but different
 ***  nonzero values. One wants to reorder a specified block
 ***  of the matrix one time and reuse that initial ordering at each
 ***  factorization. Now suppose one also wants to use a linear
 ***  solution method that messes with the ordering during factorization.
 ***  One can do the following:
 ***  Create the mtx in question, output assign and partition it.<br><br>
 ***  <pre>
 ***  Call bp=mtx_create_block_perm(mtx).
 ***   This returns a copy of the information we need to store.
 ***  Working with individual blocks is then as follows:
 *** 1- Reorder the block in question in the mtx in any way you like.
 *** 2- Call mtx_update_block_perm(mtx,block_number,bp)
 ***      This updates the saved information for rows and columns
 ***      *** within *** the block specified by copying from the mtx.
 ***      All sorts of inexpensive sanity checks will be performed to
 ***      ensure that this is a wise thing to do.
 *** 3- Factor, or otherwise do anything you want to the permutation
 ***    of the original mtx, so long as you do not mess with more than
 ***    one block. Changing the order (size) of the mtx is not permitted.
 *** 4- Call mtx_restore_block_perm(mtx,block_number,bp)
 ***     This takes information stored in bp and resets the permutation
 ***     of the mtx within block_number to match it.
 ***     All sorts of inexpensive sanity checks will be performed to
 ***     ensure that this is a wise thing to do.
 *** 5- Repeat steps 3 and 4 ad nauseum, or loop back through 1 and 2
 ***    if moving from block to block.
 *** 6- When done with the matrix, or at least the process in question,
 ***    don't forget to do an mtx_destroy_block_perm(bp).
 ***  </pre>
 ***  Other variations on this scheme are possible. See the headers for
 ***  details of the individual functions.
 **/

#define mtx_ENTIRE_MATRIX NULL
/**< 
 ***  Refers to the entire matrix in any function where a region is called for.
 **/

#define mtx_ALL_ROWS NULL
#define mtx_ALL_COLS NULL
/**< 
 ***  Refers to all rows/columns in any function where a row/column range
 ***  is called for.
 **/

#define mtx_CREATE_SPARSE NULL
/**< 
 ***  In calls that take an mtx_sparse_t as an argument AND return a
 ***  pointer to a sparse AND do not note otherwise in the header
 ***  mtx_CREATE_SPARSE can be passed the argument and the call will
 ***  allocate the sparse_t data structure it returns. In such cases,
 ***  it is the users responsibility to see that the sparse is
 ***  eventually destroyed.
 **/

#define mtx_IGNORE_ZEROES (1)
#define mtx_SOFT_ZEROES (0)
/**< 
 *** Those functions which fill a sparse_t might optionally return
 *** soft zeroes given mtx_SOFT_ZEROES.
 **/

#define mtx_FIRST      ((int32)(-1))
#define mtx_LAST       ((int32)(-1))
#define mtx_NONE       ((int32)(-1))
#define mtx_ALL_BLOCKS ((int32)(-1))

#include "mtx_basic.h"
 /**
 ***  All operations involving accounting/memory
 ***  management and CHANGES to matrix data structure.
 ***  This includes data file i/o and floating point
 ***  exception recovery handling.
 **/

#include "mtx_perms.h"
#include "mtx_reorder.h"
/***                        
 ***  Functions for permuting matrices, reordering, etc.
 ***  These probably ought to be combined. mtx_reorder is currently
 ***  thinly disguised linsol code that needs to be redone in the mtx idiom.
 **/

#include "mtx_query.h"
/***
 ***  READ_ONLY functions that don't make more sense
 ***  in mtx2_perms or mtx2_basic. Most vector data
 ***  and basic math operations fall in this file.
 **/

#include "mtx_linal.h"
/***                        
 ***  Functions which probably shouldn't exist at all,
 ***  but are simply too damned expensive to do without
 ***  knowing all the internals of the mtx. These tend
 ***  to be ultracritical path functions in linsol.
 ***  These functions generally have very, very
 ***  complex semantics and preconditions and assumptions.
 ***  You really DON'T want to add functions to this
 ***  file or use those that are in it.
 **/

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