diff options
Diffstat (limited to 'modules/sparse/includes')
| -rwxr-xr-x | modules/sparse/includes/dynlib_sparse.h | 28 | ||||
| -rwxr-xr-x | modules/sparse/includes/gw_sparse.h | 51 | ||||
| -rwxr-xr-x | modules/sparse/includes/sp.h | 52 | ||||
| -rwxr-xr-x | modules/sparse/includes/spConfig.h | 549 | ||||
| -rwxr-xr-x | modules/sparse/includes/spDefs.h | 828 | ||||
| -rwxr-xr-x | modules/sparse/includes/spmalloc.h | 33 | ||||
| -rwxr-xr-x | modules/sparse/includes/spmatrix.h | 325 |
7 files changed, 1866 insertions, 0 deletions
diff --git a/modules/sparse/includes/dynlib_sparse.h b/modules/sparse/includes/dynlib_sparse.h new file mode 100755 index 000000000..e35cc9b98 --- /dev/null +++ b/modules/sparse/includes/dynlib_sparse.h @@ -0,0 +1,28 @@ +/* +* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +* Copyright (C) DIGITEO - 2010 - Allan CORNET +* +* This file must be used under the terms of the CeCILL. +* This source file is licensed as described in the file COPYING, which +* you should have received as part of this distribution. The terms +* are also available at +* http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt +* +*/ + +/*--------------------------------------------------------------------------*/ +#ifndef __DYNLIB_SPARSE_H__ +#define __DYNLIB_SPARSE_H__ + +#ifdef _MSC_VER +#ifdef SPARSE_EXPORTS +#define SPARSE_IMPEXP __declspec(dllexport) +#else +#define SPARSE_IMPEXP __declspec(dllimport) +#endif +#else +#define SPARSE_IMPEXP +#endif + +#endif /* __DYNLIB_SPARSE_H__ */ +/*--------------------------------------------------------------------------*/ diff --git a/modules/sparse/includes/gw_sparse.h b/modules/sparse/includes/gw_sparse.h new file mode 100755 index 000000000..fc81c1ec2 --- /dev/null +++ b/modules/sparse/includes/gw_sparse.h @@ -0,0 +1,51 @@ +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2006 - INRIA - Allan CORNET + * Copyright (C) 2010 - DIGITEO - Allan CORNET + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + * + */ + +#ifndef __GW_SPARSE__ +#define __GW_SPARSE__ + +#include "machine.h" +#include "dynlib_sparse.h" +/*--------------------------------------------------------------------------*/ +SPARSE_IMPEXP int gw_sparse(void); +/*--------------------------------------------------------------------------*/ +int sci_sparsefunc (char *fname, unsigned long fname_len); +int sci_spget (char *fname, unsigned long fname_len); +int sci_full (char *fname, unsigned long fname_len); +int sci_lufact (char *fname, unsigned long fname_len); +int sci_lusolve (char *fname, unsigned long fname_len); +int sci_ludel (char *fname, unsigned long fname_len); +int sci_luget (char *fname, unsigned long fname_len); +int sci_spclean (char *fname, unsigned long fname_len); +int sci_nnz (char *fname, unsigned long fname_len); +int sci_spmax (char *fname, unsigned long fname_len); +int sci_spmin (char *fname, unsigned long fname_len); +int sci_spmatrix (char *fname, unsigned long fname_len); +int sci_spchol (char *fname, unsigned long fname_len); +int sci_fadj2sp (char *fname, unsigned long fname_len); +int sci_spcompa (char *fname, unsigned long fname_len); +int sci_ordmmd (char *fname, unsigned long fname_len); +int sci_blkfc1i (char *fname, unsigned long fname_len); +int sci_blkslvi (char *fname, unsigned long fname_len); +int sci_inpnvi (char *fname, unsigned long fname_len); +int sci_sfinit (char *fname, unsigned long fname_len); +int sci_symfcti (char *fname, unsigned long fname_len); +int sci_bfinit (char *fname, unsigned long fname_len); +int sci_msparse (char *fname, unsigned long fname_len); +int sci_mspget (char *fname, unsigned long fname_len); +int sci_mfull (char *fname, unsigned long fname_len); +extern int C2F(scita2lpd) (char *fname, unsigned long fname_len); +/*--------------------------------------------------------------------------*/ +#endif /* __GW_SPARSE__ */ +/*--------------------------------------------------------------------------*/ + diff --git a/modules/sparse/includes/sp.h b/modules/sparse/includes/sp.h new file mode 100755 index 000000000..b00d900dc --- /dev/null +++ b/modules/sparse/includes/sp.h @@ -0,0 +1,52 @@ +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) INRIA + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + * + */ + +#ifndef __SP_H__ +#define __SP_H__ + +#define SQR(x) ((x)*(x)) + +#define NB 32 /* block size for dgels */ +#define MINABSTOL 1e-8 +#define MAXITERS 100 +#define TOLC 1e-5 /* tolerance used for dual infeasibility */ +#define SIGTOL 1e-5 /* tolerance used for detecting zero steps +* dF or dZ */ +#define MINRCOND 1e-8 /* minimum rcond to declare F_i dependent */ + +/* BLAS 1 */ +double dnrm2_( ); +double ddot_( ); +void dcopy_( ); +void daxpy_( ); +void dscal_( ); + +/* BLAS 2 */ +void dgemv_( ); +void dspmv_( ); + +/* BLAS 3 */ +void dgemm_( ); + +/* LAPACK */ +void dgels_( ); +void dspgst_( ); +void dspev_( ); +void dspgv_( ); +void dtrcon_( ); + +int sp( /* int m, int L, double *F, int *blck_szs, double *c, + double *x, double *Z, double *ul, double nu, double abstol, + double reltol, double tv, int *iters, double *work, + int lwork, int *iwork, int *info */ ); + +#endif /* __SP_H__ */ diff --git a/modules/sparse/includes/spConfig.h b/modules/sparse/includes/spConfig.h new file mode 100755 index 000000000..d79b9e12d --- /dev/null +++ b/modules/sparse/includes/spConfig.h @@ -0,0 +1,549 @@ +/* + * CONFIGURATION MACRO DEFINITIONS for sparse matrix routines + * + * Author: Advising professor: + * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli + * U.C. Berkeley + * + * This file contains macros for the sparse matrix routines that are used + * to define the personality of the routines. The user is expected to + * modify this file to maximize the performance of the routines with + * his/her matrices. + * + * Macros are distinguished by using solely capital letters in their + * identifiers. This contrasts with C defined identifiers which are + * strictly lower case, and program variable and procedure names which use + * both upper and lower case. + */ + + +/* + * Revision and copyright information. + * + * Copyright (c) 1985,86,87,88 + * by Kenneth S. Kundert and the University of California. + * + * Permission to use, copy, modify, and distribute this software and + * its documentation for any purpose and without fee is hereby granted, + * provided that the copyright notices appear in all copies and + * supporting documentation and that the authors and the University of + * California are properly credited. The authors and the University of + * California make no representations as to the suitability of this + * software for any purpose. It is provided `as is', without express + * or implied warranty. + * + */ + + +#ifndef spCONFIG_DEFS +#define spCONFIG_DEFS + + + + +#ifdef spINSIDE_SPARSE +/* + * OPTIONS + * + * These are compiler options. Set each option to one to compile that + * section of the code. If a feature is not desired, set the macro + * to NO. Recommendations are given in brackets, [ignore them]. + * + * >>> Option descriptions: + * Arithmetic Precision + * The precision of the arithmetic used by Sparse can be set by + * changing changing the spREAL macro. This macro is +* contained in the file spMatrix.h. It is strongly suggested to + * used double precision with circuit simulators. Note that + * because C always performs arithmetic operations in double + * precision, the only benefit to using single precision is that + * less storage is required. There is often a noticeable speed + * penalty when using single precision. Sparse internally refers + * to a spREAL as a RealNumber. + * REAL + * This specifies that the routines are expected to handle real + * systems of equations. The routines can be compiled to handle + * both real and complex systems at the same time, but there is a + * slight speed and memory advantage if the routines are complied + * to handle only real systems of equations. + * spCOMPLEX + * This specifies that the routines will be complied to handle + * complex systems of equations. + * EXPANDABLE + * Setting this compiler flag true (1) makes the matrix + * expandable before it has been factored. If the matrix is + * expandable, then if an element is added that would be + * considered out of bounds in the current matrix, the size of + * the matrix is increased to hold that element. As a result, + * the size of the matrix need not be known before the matrix is + * built. The matrix can be allocated with size zero and + * expanded. + * TRANSLATE + * This option allows the set of external row and column numbers + * to be non-packed. In other words, the row and column numbers + * do not have to be contiguous. The priced paid for this + * flexibility is that when TRANSLATE is set true, the time + * required to initially build the matrix will be greater because + * the external row and column number must be translated into + * internal equivalents. This translation brings about other + * benefits though. First, the spGetElement() and + * spGetAdmittance() routines may be used after the matrix has + * been factored. Further, elements, and even rows and columns, + * may be added to the matrix, and row and columns may be deleted + * from the matrix, after it has been factored. Note that when + * the set of row and column number is not a packed set, neither + * are the RHS and Solution vectors. Thus the size of these + * vectors must be at least as large as the external size, which + * is the value of the largest given row or column numbers. + * INITIALIZE + * Causes the spInitialize(), spGetInitInfo(), and + * spInstallInitInfo() routines to be compiled. These routines + * allow the user to store and read one pointer in each nonzero + * element in the matrix. spInitialize() then calls a user + * specified function for each structural nonzero in the matrix, + * and includes this pointer as well as the external row and + * column numbers as arguments. This allows the user to write + * custom matrix initialization routines. + * DIAGONAL_PIVOTING + * Many matrices, and in particular node- and modified-node + * admittance matrices, tend to be nearly symmetric and nearly + * diagonally dominant. For these matrices, it is a good idea to + * select pivots from the diagonal. With this option enabled, + * this is exactly what happens, though if no satisfactory pivot + * can be found on the diagonal, an off-diagonal pivot will be + * used. If this option is disabled, Sparse does not + * preferentially search the diagonal. Because of this, Sparse + * has a wider variety of pivot candidates available, and so + * presumably fewer fill-ins will be created. However, the + * initial pivot selection process will take considerably longer. + * If working with node admittance matrices, or other matrices + * with a strong diagonal, it is probably best to use + * DIAGONAL_PIVOTING for two reasons. First, accuracy will be + * better because pivots will be chosen from the large diagonal + * elements, thus reducing the chance of growth. Second, a near + * optimal ordering will be chosen quickly. If the class of + * matrices you are working with does not have a strong diagonal, + * do not use DIAGONAL_PIVOTING, but consider using a larger + * threshold. When DIAGONAL_PIVOTING is turned off, the following + * options and constants are not used: MODIFIED_MARKOWITZ, + * MAX_MARKOWITZ_TIES, and TIES_MULTIPLIER. + * ARRAY_OFFSET + * This determines whether arrays start at an index of zero or one. + * This option is necessitated by the fact that standard C + * convention dictates that arrays begin with an index of zero but + * the standard mathematic convention states that arrays begin with + * an index of one. So if you prefer to start your arrays with + * zero, or your calling Sparse from FORTRAN, set ARRAY_OFFSET to + * NO or 0. Otherwise, set ARRAY_OFFSET to YES or 1. Note that if + * you use an offset of one, the arrays that you pass to Sparse + * must have an allocated length of one plus the size of the + * matrix. ARRAY_OFFSET must be either 0 or 1, no other offsets + * are valid. + * spSEPARATED_COMPLEX_VECTORS + * This specifies the format for complex vectors. If this is set + * false then a complex vector is made up of one double sized + * array of RealNumber's in which the real and imaginary numbers + * are placed in the alternately array in the array. In other + * words, the first entry would be Complex[1].Real, then comes + * Complex[1].Imag, then Complex[1].Real, etc. If + * spSEPARATED_COMPLEX_VECTORS is set true, then each complex + * vector is represented by two arrays of RealNumbers, one with + * the real terms, the other with the imaginary. [NO] + * MODIFIED_MARKOWITZ + * This specifies that the modified Markowitz method of pivot + * selection is to be used. The modified Markowitz method differs + * from standard Markowitz in two ways. First, under modified + * Markowitz, the search for a pivot can be terminated early if a + * adequate (in terms of sparsity) pivot candidate is found. + * Thus, when using modified Markowitz, the initial factorization + * can be faster, but at the expense of a suboptimal pivoting + * order that may slow subsequent factorizations. The second + * difference is in the way modified Markowitz breaks Markowitz + * ties. When two or more elements are pivot candidates and they + * all have the same Markowitz product, then the tie is broken by + * choosing the element that is best numerically. The numerically + * best element is the one with the largest ratio of its magnitude + * to the magnitude of the largest element in the same column, + * excluding itself. The modified Markowitz method results in + * marginally better accuracy. This option is most appropriate + * for use when working with very large matrices where the initial + * factor time represents an unacceptable burden. [NO] + * DELETE + * This specifies that the spDeleteRowAndCol() routine + * should be compiled. Note that for this routine to be + * compiled, both DELETE and TRANSLATE should be set true. + * STRIP + * This specifies that the spStripFills() routine should be compiled. + * MODIFIED_NODAL + * This specifies that the routine that preorders modified node + * admittance matrices should be compiled. This routine results + * in greater speed and accuracy if used with this type of + * matrix. + * QUAD_ELEMENT + * This specifies that the routines that allow four related + * elements to be entered into the matrix at once should be + * compiled. These elements are usually related to an + * admittance. The routines affected by QUAD_ELEMENT are the + * spGetAdmittance, spGetQuad and spGetOnes routines. + * TRANSPOSE + * This specifies that the routines that solve the matrix as if + * it was transposed should be compiled. These routines are + * useful when performing sensitivity analysis using the adjoint + * method. + * SCALING + * This specifies that the routine that performs scaling on the + * matrix should be complied. Scaling is not strongly + * supported. The routine to scale the matrix is provided, but + * no routines are provided to scale and descale the RHS and + * Solution vectors. It is suggested that if scaling is desired, + * it only be preformed when the pivot order is being chosen [in + * spOrderAndFactor()]. This is the only time scaling has + * an effect. The scaling may then either be removed from the + * solution by the user or the scaled factors may simply be + * thrown away. [NO] + * DOCUMENTATION + * This specifies that routines that are used to document the + * matrix, such as spPrint() and spFileMatrix(), should be + * compiled. + * DETERMINANT + * This specifies that the routine spDeterminant() should be complied. + * STABILITY + * This specifies that spLargestElement() and spRoundoff() should + * be compiled. These routines are used to check the stability (and + * hence the quality of the pivoting) of the factorization by + * computing a bound on the size of the element is the matrix E = + * A - LU. If this bound is very high after applying + * spOrderAndFactor(), then the pivot threshold should be raised. + * If the bound increases greatly after using spFactor(), then the + * matrix should probably be reordered. + * CONDITION + * This specifies that spCondition() and spNorm(), the code that + * computes a good estimate of the condition number of the matrix, + * should be compiled. + * PSEUDOCONDITION + * This specifies that spPseudoCondition(), the code that computes + * a crude and easily fooled indicator of ill-conditioning in the + * matrix, should be compiled. + * MULTIPLICATION + * This specifies that the routines to multiply the unfactored + * matrix by a vector should be compiled. + * FORTRAN + * This specifies that the FORTRAN interface routines should be + * compiled. When interfacing to FORTRAN programs, the ARRAY_OFFSET + * options should be set to NO. + * DEBUG + * This specifies that additional error checking will be compiled. + * The type of error checked are those that are common when the + * matrix routines are first integrated into a user's program. Once + * the routines have been integrated in and are running smoothly, this + * option should be turned off. + */ +/* Begin options. */ +#define REAL YES +#define EXPANDABLE YES +#define TRANSLATE YES +#define INITIALIZE YES +#define DIAGONAL_PIVOTING YES +#define ARRAY_OFFSET NO +#define MODIFIED_MARKOWITZ NO +#define SPARSEDELETE YES +#define STRIP YES +#define MODIFIED_NODAL YES +#define QUAD_ELEMENT YES +#define TRANSPOSE YES +#define SCALING YES +#define DOCUMENTATION YES +#define MULTIPLICATION YES +#define DETERMINANT YES +#define STABILITY YES +#define CONDITION YES +#define PSEUDOCONDITION YES +#define FORTRAN YES +#define DEBUG NO + +/* + * The following options affect Sparse exports and so are exported as a + * side effect. For this reason they use the `sp' prefix. The boolean + * constants YES an NO are not defined in spMatrix.h to avoid conflicts + * with user code, so use 0 for NO and 1 for YES. + */ +#endif /* spINSIDE_SPARSE */ +#define spCOMPLEX 1 +#define spSEPARATED_COMPLEX_VECTORS 0 +#ifdef spINSIDE_SPARSE + + + + + + + +/* + * MATRIX CONSTANTS + * + * These constants are used throughout the sparse matrix routines. They + * should be set to suit the type of matrix being solved. Recommendations + * are given in brackets. + * + * Some terminology should be defined. The Markowitz row count is the number + * of non-zero elements in a row excluding the one being considered as pivot. + * There is one Markowitz row count for every row. The Markowitz column + * is defined similarly for columns. The Markowitz product for an element + * is the product of its row and column counts. It is a measure of how much + * work would be required on the next step of the factorization if that + * element were chosen to be pivot. A small Markowitz product is desirable. + * + * >>> Constants descriptions: + * DEFAULT_THRESHOLD + * The relative threshold used if the user enters an invalid + * threshold. Also the threshold used by spFactor() when + * calling spOrderAndFactor(). The default threshold should + * not be less than or equal to zero nor larger than one. [0.001] + * DIAG_PIVOTING_AS_DEFAULT + * This indicates whether spOrderAndFactor() should use diagonal + * pivoting as default. This issue only arises when + * spOrderAndFactor() is called from spFactor(). + * SPACE_FOR_ELEMENTS + * This number multiplied by the size of the matrix equals the number + * of elements for which memory is initially allocated in + * spCreate(). [6] + * SPACE_FOR_FILL_INS + * This number multiplied by the size of the matrix equals the number + * of elements for which memory is initially allocated and specifically + * reserved for fill-ins in spCreate(). [4] + * ELEMENTS_PER_ALLOCATION + * The number of matrix elements requested from the malloc utility on + * each call to it. Setting this value greater than 1 reduces the + * amount of overhead spent in this system call. On a virtual memory + * machine, its good to allocate slightly less than a page worth of + * elements at a time (or some multiple thereof). + * [For the VAX, for real only use 41, otherwise use 31] + * MINIMUM_ALLOCATED_SIZE + * The minimum allocated size of a matrix. Note that this does not + * limit the minimum size of a matrix. This just prevents having to + * resize a matrix many times if the matrix is expandable, large and + * allocated with an estimated size of zero. This number should not + * be less than one. + * EXPANSION_FACTOR + * The amount the allocated size of the matrix is increased when it + * is expanded. + * MAX_MARKOWITZ_TIES + * This number is used for two slightly different things, both of which + * relate to the search for the best pivot. First, it is the maximum + * number of elements that are Markowitz tied that will be sifted + * through when trying to find the one that is numerically the best. + * Second, it creates an upper bound on how large a Markowitz product + * can be before it eliminates the possibility of early termination + * of the pivot search. In other words, if the product of the smallest + * Markowitz product yet found and TIES_MULTIPLIER is greater than + * MAX_MARKOWITZ_TIES, then no early termination takes place. + * Set MAX_MARKOWITZ_TIES to some small value if no early termination of + * the pivot search is desired. An array of RealNumbers is allocated + * of size MAX_MARKOWITZ_TIES so it must be positive and shouldn't + * be too large. Active when MODIFIED_MARKOWITZ is 1 (true). [100] + * TIES_MULTIPLIER + * Specifies the number of Markowitz ties that are allowed to occur + * before the search for the pivot is terminated early. Set to some + * large value if no early termination of the pivot search is desired. + * This number is multiplied times the Markowitz product to determine + * how many ties are required for early termination. This means that + * more elements will be searched before early termination if a large + * number of fill-ins could be created by accepting what is currently + * considered the best choice for the pivot. Active when + * MODIFIED_MARKOWITZ is 1 (true). Setting this number to zero + * effectively eliminates all pivoting, which should be avoided. + * This number must be positive. TIES_MULTIPLIER is also used when + * diagonal pivoting breaks down. [5] + * DEFAULT_PARTITION + * Which partition mode is used by spPartition() as default. + * Possibilities include + * spDIRECT_PARTITION -- each row used direct addressing, best for + * a few relatively dense matrices. + * spINDIRECT_PARTITION -- each row used indirect addressing, best + * for a few very sparse matrices. + * spAUTO_PARTITION -- direct or indirect addressing is chosen on + * a row-by-row basis, carries a large overhead, but speeds up + * both dense and sparse matrices, best if there is a large + * number of matrices that can use the same ordering. + */ + +/* Begin constants. */ +#define DEFAULT_THRESHOLD 1.0e-3 +#define DIAG_PIVOTING_AS_DEFAULT YES +#define SPACE_FOR_ELEMENTS 6 +#define SPACE_FOR_FILL_INS 4 +#define ELEMENTS_PER_ALLOCATION 31 +#define MINIMUM_ALLOCATED_SIZE 6 +#define EXPANSION_FACTOR 1.5 +#define MAX_MARKOWITZ_TIES 100 +#define TIES_MULTIPLIER 5 +#define DEFAULT_PARTITION spAUTO_PARTITION + + + + + + +/* + * PRINTER WIDTH + * + * This macro characterize the printer for the spPrint() routine. + * + * >>> Macros: + * PRINTER_WIDTH + * The number of characters per page width. Set to 80 for terminal, + * 132 for line printer. + */ + +/* Begin printer constants. */ +#define PRINTER_WIDTH 80 + + + + + + +/* + * MACHINE CONSTANTS + * + * These numbers must be updated when the program is ported to a new machine. + */ + +/* Begin machine constants. */ + +#ifdef notdef /* __STDC__ */ +/* + * This code is currently deleted because most ANSI standard C compilers + * do not provide the standard header files yet. + */ +# include <limits.h> +# include <float.h> +# define MACHINE_RESOLUTION DBL_EPSILON +# define LARGEST_REAL DBL_MAX +# define SMALLEST_REAL DBL_MIN +# define LARGEST_SHORT_INTEGER SHRT_MAX +# define LARGEST_LONG_INTEGER LONG_MAX +#else /* NOT defined(__STDC__) */ + +/* Apple MacOSX */ + +#ifdef __APPLE__ /* __STDC__ */ +# include <limits.h> +# include <float.h> +# define MACHINE_RESOLUTION DBL_EPSILON +# define LARGEST_REAL DBL_MAX +# define SMALLEST_REAL DBL_MIN +# define LARGEST_SHORT_INTEGER SHRT_MAX +# define LARGEST_LONG_INTEGER LONG_MAX +#endif /* NOT defined(__STDC__) */ + +/* VAX machine constants */ +#if (defined(vax) && !defined(netbsd)) +# define MACHINE_RESOLUTION 6.93889e-18 +# define LARGEST_REAL 1.70141e+38 +# define SMALLEST_REAL 2.938743e-39 +# define LARGEST_SHORT_INTEGER 32766 +# define LARGEST_LONG_INTEGER 2147483646 +#endif + +/* MIPS machine constants */ +#if (defined(mips) && !defined(netbsd)) +# define MACHINE_RESOLUTION 6.93889e-18 +# define LARGEST_REAL 1.70141e+38 +# define SMALLEST_REAL 2.938743e-39 +# define LARGEST_SHORT_INTEGER 32766 +# define LARGEST_LONG_INTEGER 2147483646 +#endif + + +/* hp9000 machine constants */ +#ifdef hpux +/* These values are correct for hp9000/300. Should be correct for others. */ +# define MACHINE_RESOLUTION 8.9e-15 +# define LARGEST_REAL 1.79769313486231e+308 +# define SMALLEST_REAL 2.22507385850721e-308 +# define LARGEST_SHORT_INTEGER 32766 +# define LARGEST_LONG_INTEGER 2147483646 +#endif + +/* IBM machine constants */ +#ifdef aix +/** The STDC option works on aix on gives the values that i've copied here */ +# define MACHINE_RESOLUTION 2.2204460492503131e-16 +# define LARGEST_REAL 1.7976931348623158e+308 +# define SMALLEST_REAL 2.2250738585072014e-308 +# define LARGEST_SHORT_INTEGER 32767 +# define LARGEST_LONG_INTEGER 2147483647 +#endif + +/* Sun machine constants */ +#if (defined(sun) && !defined(netbsd)) +/* These values are rumored to be the correct values. */ +# define MACHINE_RESOLUTION 8.9e-15 +# define LARGEST_REAL 1.79769313486231e+308 +# define SMALLEST_REAL 2.22507385850721e-308 +# define LARGEST_SHORT_INTEGER 32766 +# define LARGEST_LONG_INTEGER 2147483646 +#endif +/* DEC alpha machine constant*/ +#if (defined(__alpha) && !defined(netbsd)) +# include <limits.h> +# include <float.h> +# define MACHINE_RESOLUTION DBL_EPSILON +# define LARGEST_REAL DBL_MAX +# define SMALLEST_REAL DBL_MIN +# define LARGEST_SHORT_INTEGER SHRT_MAX +# define LARGEST_LONG_INTEGER LONG_MAX +#endif +#ifdef linux +# include <limits.h> +# include <float.h> +# define MACHINE_RESOLUTION DBL_EPSILON +# define LARGEST_REAL DBL_MAX +# define SMALLEST_REAL DBL_MIN +# define LARGEST_SHORT_INTEGER SHRT_MAX +# define LARGEST_LONG_INTEGER LONG_MAX +#endif +#if defined(netbsd) || defined(freebsd) +# include <limits.h> +# include <float.h> +# define MACHINE_RESOLUTION DBL_EPSILON +# define LARGEST_REAL DBL_MAX +# define SMALLEST_REAL DBL_MIN +# define LARGEST_SHORT_INTEGER SHRT_MAX +# define LARGEST_LONG_INTEGER LONG_MAX +#endif +#ifdef _MSC_VER +# include <limits.h> +# include <float.h> +# define MACHINE_RESOLUTION DBL_EPSILON +# define LARGEST_REAL DBL_MAX +# define SMALLEST_REAL DBL_MIN +/* XXXXX : a v'erifier */ +# define LARGEST_SHORT_INTEGER 32766 +# define LARGEST_LONG_INTEGER 2147483646 +#endif +#endif /* NOT defined(__STDC__) */ + +/* + * ANNOTATION + * + * This macro changes the amount of annotation produced by the matrix + * routines. The annotation is used as a debugging aid. Change the number + * associated with ANNOTATE to change the amount of annotation produced by + * the program. + */ + +/* Begin annotation definitions. */ +#define ANNOTATE NONE + +#define NONE 0 +#define ON_STRANGE_BEHAVIOR 1 +#define FULL 2 + +#endif /* spINSIDE_SPARSE */ + +#endif /* spCONFIG_DEFS */ + + + diff --git a/modules/sparse/includes/spDefs.h b/modules/sparse/includes/spDefs.h new file mode 100755 index 000000000..e2e57b6ea --- /dev/null +++ b/modules/sparse/includes/spDefs.h @@ -0,0 +1,828 @@ +#ifndef __SPDEFS_H__ +#define __SPDEFS_H__ +#include "core_math.h" +/* + * DATA STRUCTURE AND MACRO DEFINITIONS for Sparse. + * + * Author: Advising professor: + * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli + * UC Berkeley + * + * This file contains common type definitions and macros for the sparse + * matrix routines. These definitions are of no interest to the user. + */ + + +/* + * Revision and copyright information. + * + * Copyright (c) 1985,86,87,88 + * by Kenneth S. Kundert and the University of California. + * + * Permission to use, copy, modify, and distribute this software and + * its documentation for any purpose and without fee is hereby granted, + * provided that the copyright notices appear in all copies and + * supporting documentation and that the authors and the University of + * California are properly credited. The authors and the University of + * California make no representations as to the suitability of this + * software for any purpose. It is provided `as is', without express + * or implied warranty. + * + */ + + + + +/* + * IMPORTS + */ + +/* + * If running lint, change some of the compiler options to get a more + * complete inspection. + */ + +#ifdef lint +#undef REAL +#undef spCOMPLEX +#undef EXPANDABLE +#undef TRANSLATE +#undef INITIALIZE +#undef SPARSEDELETE +#undef STRIP +#undef MODIFIED_NODAL +#undef QUAD_ELEMENT +#undef TRANSPOSE +#undef SCALING +#undef DOCUMENTATION +#undef MULTIPLICATION +#undef DETERMINANT +#undef CONDITION +#undef PSEUDOCONDITION +#undef FORTRAN +#undef DEBUG + +#define REAL YES +#define spCOMPLEX YES +#define EXPANDABLE YES +#define TRANSLATE YES +#define INITIALIZE YES +#define SPARSEDELETE YES +#define STRIP YES +#define MODIFIED_NODAL YES +#define QUAD_ELEMENT YES +#define TRANSPOSE YES +#define SCALING YES +#define DOCUMENTATION YES +#define MULTIPLICATION YES +#define DETERMINANT YES +#define CONDITION YES +#define PSEUDOCONDITION YES +#define FORTRAN YES +#define DEBUG YES + +#define LINT YES +#else /* not lint */ +#define LINT NO +#endif /* not lint */ + + + +/* + * MACRO DEFINITIONS + * + * Macros are distinguished by using solely capital letters in their + * identifiers. This contrasts with C defined identifiers which are strictly + * lower case, and program variable and procedure names which use both upper + * and lower case. + */ + +/* Begin macros. */ + +/* Boolean data type */ + +#define SPBOOLEAN int +#define NO 0 +#define YES 1 +#define NOT ! +#define AND && +#define OR || + +/* NULL pointer */ +#ifndef NULL +#define NULL 0 +#endif + +#define SPARSE_ID 0x772773 /* Arbitrary (is Sparse on phone). */ +#define IS_SPARSE(matrix) ((matrix) != NULL && \ + (matrix)->ID == SPARSE_ID) +#define IS_VALID(matrix) ((matrix) != NULL && \ + (matrix)->ID == SPARSE_ID && \ + (matrix)->Error >= spOKAY && \ + (matrix)->Error < spFATAL) +#define IS_FACTORED(matrix) ((matrix)->Factored && !(matrix)->NeedsOrdering) + +/* Macro commands */ +/* Macro function that returns the square of a number. */ +#define SQR(a) ((a)*(a)) + +/* Macro procedure that swaps two entities. */ +#define SWAP(type, a, b) {type swapx; swapx = a; a = b; b = swapx;} + +/* Macro function that returns the approx absolute value of a complex number. */ +#if spCOMPLEX +#define ELEMENT_MAG(ptr) (Abs((ptr)->Real) + Abs((ptr)->Imag)) +#else +#define ELEMENT_MAG(ptr) ((ptr)->Real < 0.0 ? -(ptr)->Real : (ptr)->Real) +#endif + +/* Complex assignment statements. */ +#define CMPLX_ASSIGN(to,from) \ +{ (to).Real = (from).Real; \ + (to).Imag = (from).Imag; \ +} +#define CMPLX_CONJ_ASSIGN(to,from) \ +{ (to).Real = (from).Real; \ + (to).Imag = -(from).Imag; \ +} +#define CMPLX_NEGATE_ASSIGN(to,from) \ +{ (to).Real = -(from).Real; \ + (to).Imag = -(from).Imag; \ +} +#define CMPLX_CONJ_NEGATE_ASSIGN(to,from) \ +{ (to).Real = -(from).Real; \ + (to).Imag = (from).Imag; \ +} +#define CMPLX_CONJ(a) (a).Imag = -(a).Imag +#define CMPLX_NEGATE(a) \ +{ (a).Real = -(a).Real; \ + (a).Imag = -(a).Imag; \ +} + +/* Macro that returns the approx magnitude (L-1 norm) of a complex number. */ +#define CMPLX_1_NORM(a) (Abs((a).Real) + Abs((a).Imag)) + +/* Macro that returns the approx magnitude (L-infinity norm) of a complex. */ +#define CMPLX_INF_NORM(a) (Max (Abs((a).Real),Abs((a).Imag))) + +/* Macro function that returns the magnitude (L-2 norm) of a complex number. */ +#define CMPLX_2_NORM(a) (sqrt((a).Real*(a).Real + (a).Imag*(a).Imag)) + +/* Macro function that performs complex addition. */ +#define CMPLX_ADD(to,from_a,from_b) \ +{ (to).Real = (from_a).Real + (from_b).Real; \ + (to).Imag = (from_a).Imag + (from_b).Imag; \ +} + +/* Macro function that performs complex subtraction. */ +#define CMPLX_SUBT(to,from_a,from_b) \ +{ (to).Real = (from_a).Real - (from_b).Real; \ + (to).Imag = (from_a).Imag - (from_b).Imag; \ +} + +/* Macro function that is equivalent to += operator for complex numbers. */ +#define CMPLX_ADD_ASSIGN(to,from) \ +{ (to).Real += (from).Real; \ + (to).Imag += (from).Imag; \ +} + +/* Macro function that is equivalent to -= operator for complex numbers. */ +#define CMPLX_SUBT_ASSIGN(to,from) \ +{ (to).Real -= (from).Real; \ + (to).Imag -= (from).Imag; \ +} + +/* Macro function that multiplies a complex number by a scalar. */ +#define SCLR_MULT(to,sclr,cmplx) \ +{ (to).Real = (sclr) * (cmplx).Real; \ + (to).Imag = (sclr) * (cmplx).Imag; \ +} + +/* Macro function that multiply-assigns a complex number by a scalar. */ +#define SCLR_MULT_ASSIGN(to,sclr) \ +{ (to).Real *= (sclr); \ + (to).Imag *= (sclr); \ +} + +/* Macro function that multiplies two complex numbers. */ +#define CMPLX_MULT(to,from_a,from_b) \ +{ (to).Real = (from_a).Real * (from_b).Real - \ + (from_a).Imag * (from_b).Imag; \ + (to).Imag = (from_a).Real * (from_b).Imag + \ + (from_a).Imag * (from_b).Real; \ +} + +/* Macro function that implements to *= from for complex numbers. */ +#define CMPLX_MULT_ASSIGN(to,from) \ +{ RealNumber to_real_ = (to).Real; \ + (to).Real = to_real_ * (from).Real - \ + (to).Imag * (from).Imag; \ + (to).Imag = to_real_ * (from).Imag + \ + (to).Imag * (from).Real; \ +} + +/* Macro function that multiplies two complex numbers, the first of which is + * conjugated. */ +#define CMPLX_CONJ_MULT(to,from_a,from_b) \ +{ (to).Real = (from_a).Real * (from_b).Real + \ + (from_a).Imag * (from_b).Imag; \ + (to).Imag = (from_a).Real * (from_b).Imag - \ + (from_a).Imag * (from_b).Real; \ +} + +/* Macro function that multiplies two complex numbers and then adds them + * to another. to = add + mult_a * mult_b */ +#define CMPLX_MULT_ADD(to,mult_a,mult_b,add) \ +{ (to).Real = (mult_a).Real * (mult_b).Real - \ + (mult_a).Imag * (mult_b).Imag + (add).Real; \ + (to).Imag = (mult_a).Real * (mult_b).Imag + \ + (mult_a).Imag * (mult_b).Real + (add).Imag; \ +} + +/* Macro function that subtracts the product of two complex numbers from + * another. to = subt - mult_a * mult_b */ +#define CMPLX_MULT_SUBT(to,mult_a,mult_b,subt) \ +{ (to).Real = (subt).Real - (mult_a).Real * (mult_b).Real + \ + (mult_a).Imag * (mult_b).Imag; \ + (to).Imag = (subt).Imag - (mult_a).Real * (mult_b).Imag - \ + (mult_a).Imag * (mult_b).Real; \ +} + +/* Macro function that multiplies two complex numbers and then adds them + * to another. to = add + mult_a* * mult_b where mult_a* represents mult_a + * conjugate. */ +#define CMPLX_CONJ_MULT_ADD(to,mult_a,mult_b,add) \ +{ (to).Real = (mult_a).Real * (mult_b).Real + \ + (mult_a).Imag * (mult_b).Imag + (add).Real; \ + (to).Imag = (mult_a).Real * (mult_b).Imag - \ + (mult_a).Imag * (mult_b).Real + (add).Imag; \ +} + +/* Macro function that multiplies two complex numbers and then adds them + * to another. to += mult_a * mult_b */ +#define CMPLX_MULT_ADD_ASSIGN(to,from_a,from_b) \ +{ (to).Real += (from_a).Real * (from_b).Real - \ + (from_a).Imag * (from_b).Imag; \ + (to).Imag += (from_a).Real * (from_b).Imag + \ + (from_a).Imag * (from_b).Real; \ +} + +/* Macro function that multiplies two complex numbers and then subtracts them + * from another. */ +#define CMPLX_MULT_SUBT_ASSIGN(to,from_a,from_b) \ +{ (to).Real -= (from_a).Real * (from_b).Real - \ + (from_a).Imag * (from_b).Imag; \ + (to).Imag -= (from_a).Real * (from_b).Imag + \ + (from_a).Imag * (from_b).Real; \ +} + +/* Macro function that multiplies two complex numbers and then adds them + * to the destination. to += from_a* * from_b where from_a* represents from_a + * conjugate. */ +#define CMPLX_CONJ_MULT_ADD_ASSIGN(to,from_a,from_b) \ +{ (to).Real += (from_a).Real * (from_b).Real + \ + (from_a).Imag * (from_b).Imag; \ + (to).Imag += (from_a).Real * (from_b).Imag - \ + (from_a).Imag * (from_b).Real; \ +} + +/* Macro function that multiplies two complex numbers and then subtracts them + * from the destination. to -= from_a* * from_b where from_a* represents from_a + * conjugate. */ +#define CMPLX_CONJ_MULT_SUBT_ASSIGN(to,from_a,from_b) \ +{ (to).Real -= (from_a).Real * (from_b).Real + \ + (from_a).Imag * (from_b).Imag; \ + (to).Imag -= (from_a).Real * (from_b).Imag - \ + (from_a).Imag * (from_b).Real; \ +} + +/* + * Macro functions that provide complex division. + */ + +/* Complex division: to = num / den */ +#define CMPLX_DIV(to,num,den) \ +{ RealNumber r_, s_; \ + if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ + ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ + { r_ = (den).Imag / (den).Real; \ + s_ = (den).Real + r_*(den).Imag; \ + (to).Real = ((num).Real + r_*(num).Imag)/s_; \ + (to).Imag = ((num).Imag - r_*(num).Real)/s_; \ + } \ + else \ + { r_ = (den).Real / (den).Imag; \ + s_ = (den).Imag + r_*(den).Real; \ + (to).Real = (r_*(num).Real + (num).Imag)/s_; \ + (to).Imag = (r_*(num).Imag - (num).Real)/s_; \ + } \ +} + +/* Complex division and assignment: num /= den */ +#define CMPLX_DIV_ASSIGN(num,den) \ +{ RealNumber r_, s_, t_; \ + if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ + ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ + { r_ = (den).Imag / (den).Real; \ + s_ = (den).Real + r_*(den).Imag; \ + t_ = ((num).Real + r_*(num).Imag)/s_; \ + (num).Imag = ((num).Imag - r_*(num).Real)/s_; \ + (num).Real = t_; \ + } \ + else \ + { r_ = (den).Real / (den).Imag; \ + s_ = (den).Imag + r_*(den).Real; \ + t_ = (r_*(num).Real + (num).Imag)/s_; \ + (num).Imag = (r_*(num).Imag - (num).Real)/s_; \ + (num).Real = t_; \ + } \ +} + +/* Complex reciprocation: to = 1.0 / den */ +#define CMPLX_RECIPROCAL(to,den) \ +{ RealNumber r_; \ + if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ + ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ + { r_ = (den).Imag / (den).Real; \ + (to).Imag = -r_*((to).Real = 1.0/((den).Real + r_*(den).Imag)); \ + } \ + else \ + { r_ = (den).Real / (den).Imag; \ + (to).Real = -r_*((to).Imag = -1.0/((den).Imag + r_*(den).Real));\ + } \ +} + + + + + + +/* + * ASSERT and ABORT + * + * Macro used to assert that if the code is working correctly, then + * a condition must be true. If not, then execution is terminated + * and an error message is issued stating that there is an internal + * error and giving the file and line number. These assertions are + * not evaluated unless the DEBUG flag is true. + */ + +#if DEBUG +#define ASSERT(condition) if (NOT(condition)) ABORT() +#else +#define ASSERT(condition) +#endif + +#if DEBUG +#define ABORT() \ +{ (void)fflush(stdout); \ + (void)fprintf(stderr, _("sparse: panic in file `%s' at line %d.\n"), \ + __FILE__, __LINE__); \ + (void)fflush(stderr); \ + abort(); \ +} +#else +#define ABORT() +#endif + + + + + +/* + * IMAGINARY VECTORS + * + * The imaginary vectors iRHS and iSolution are only needed when the + * options spCOMPLEX and spSEPARATED_COMPLEX_VECTORS are set. The following + * macro makes it easy to include or exclude these vectors as needed. + */ + +#if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS +#define IMAG_VECTORS , iRHS, iSolution +#define IMAG_RHS , iRHS +#else +#define IMAG_VECTORS +#define IMAG_RHS +#endif + + +/* + * REAL NUMBER + */ + +/* Begin `RealNumber'. */ + +typedef spREAL RealNumber, *RealVector; + +/* + * COMPLEX NUMBER DATA STRUCTURE + * + * >>> Structure fields: + * Real (RealNumber) + * The real portion of the number. Real must be the first + * field in this structure. + * Imag (RealNumber) + * The imaginary portion of the number. This field must follow + * immediately after Real. + */ + +/* Begin `ComplexNumber'. */ + +typedef struct +{ + RealNumber Real; + RealNumber Imag; +} ComplexNumber, *ComplexVector; + + + + + + + + +/* + * MATRIX ELEMENT DATA STRUCTURE + * + * Every nonzero element in the matrix is stored in a dynamically allocated + * MatrixElement structure. These structures are linked together in an + * orthogonal linked list. Two different MatrixElement structures exist. + * One is used when only real matrices are expected, it is missing an entry + * for imaginary data. The other is used if complex matrices are expected. + * It contains an entry for imaginary data. + * + * >>> Structure fields: + * Real (RealNumber) + * The real portion of the value of the element. Real must be the first + * field in this structure. + * Imag (RealNumber) + * The imaginary portion of the value of the element. If the matrix + * routines are not compiled to handle complex matrices, then this + * field does not exist. If it exists, it must follow immediately after + * Real. + * Row (int) + * The row number of the element. + * Col (int) + * The column number of the element. + * NextInRow (struct MatrixElement *) + * NextInRow contains a pointer to the next element in the row to the + * right of this element. If this element is the last nonzero in the + * row then NextInRow contains NULL. + * NextInCol (struct MatrixElement *) + * NextInCol contains a pointer to the next element in the column below + * this element. If this element is the last nonzero in the column then + * NextInCol contains NULL. + * pInitInfo (char *) + * Pointer to user data used for initialization of the matrix element. + * Initialized to NULL. + * + * >>> Type definitions: + * ElementPtr + * A pointer to a MatrixElement. + * ArrayOfElementPtrs + * An array of ElementPtrs. Used for FirstInRow, FirstInCol and + * Diag pointer arrays. + */ + +/* Begin `MatrixElement'. */ + +struct MatrixElement +{ + RealNumber Real; +#if spCOMPLEX + RealNumber Imag; +#endif + int Row; + int Col; + struct MatrixElement *NextInRow; + struct MatrixElement *NextInCol; +#if INITIALIZE + char *pInitInfo; +#endif +}; + +typedef struct MatrixElement *ElementPtr; +typedef ElementPtr *ArrayOfElementPtrs; + + + + + + + + +/* + * SPALLOCATION DATA STRUCTURE + * + * The sparse matrix routines keep track of all memory that is allocated by + * the operating system so the memory can later be freed. This is done by + * saving the pointers to all the chunks of memory that are allocated to a + * particular matrix in an allocation list. That list is organized as a + * linked list so that it can grow without a priori bounds. + * + * >>> Structure fields: + * AllocatedPtr (char *) + * Pointer to chunk of memory that has been allocated for the matrix. + * NextRecord (struct AllocationRecord *) + * Pointer to the next allocation record. + */ + +/* Begin `AllocationRecord'. */ +struct AllocationRecord +{ + char *AllocatedPtr; + struct AllocationRecord *NextRecord; +}; + +typedef struct AllocationRecord *AllocationListPtr; + + + + + + + + + +/* + * FILL-IN LIST DATA STRUCTURE + * + * The sparse matrix routines keep track of all fill-ins separately from + * user specified elements so they may be removed by spStripFills(). Fill-ins + * are allocated in bunched in what is called a fill-in lists. The data + * structure defined below is used to organize these fill-in lists into a + * linked-list. + * + * >>> Structure fields: + * pFillinList (ElementPtr) + * Pointer to a fill-in list, or a bunch of fill-ins arranged contiguously + * in memory. + * NumberOfFillinsInList (int) + * Seems pretty self explanatory to me. + * Next (struct FillinListNodeStruct *) + * Pointer to the next fill-in list structures. + */ + +/* Begin `FillinListNodeStruct'. */ +struct FillinListNodeStruct +{ + ElementPtr pFillinList; + int NumberOfFillinsInList; + struct FillinListNodeStruct *Next; +}; + + + + + + + + + + +/* + * MATRIX FRAME DATA STRUCTURE + * + * This structure contains all the pointers that support the orthogonal + * linked list that contains the matrix elements. Also included in this + * structure are other numbers and pointers that are used globally by the + * sparse matrix routines and are associated with one particular matrix. + * + * >>> Type definitions: + * MatrixPtr + * A pointer to MatrixFrame. Essentially, a pointer to the matrix. + * + * >>> Structure fields: + * AbsThreshold (RealNumber) + * The absolute magnitude an element must have to be considered as a + * pivot candidate, except as a last resort. + * AllocatedExtSize (int) + * The allocated size of the arrays used to translate external row and + * column numbers to their internal values. + * AllocatedSize (int) + * The currently allocated size of the matrix; the size the matrix can + * grow to when EXPANDABLE is set true and AllocatedSize is the largest + * the matrix can get without requiring that the matrix frame be + * reallocated. + * Complex (SPBOOLEAN) + * The flag which indicates whether the matrix is complex (true) or + * real. + * CurrentSize (int) + * This number is used during the building of the matrix when the + * TRANSLATE option is set true. It indicates the number of internal + * rows and columns that have elements in them. + * Diag (ArrayOfElementPtrs) + * Array of pointers that points to the diagonal elements. + * DoCmplxDirect (SPBOOLEAN *) + * Array of flags, one for each column in matrix. If a flag is true + * then corresponding column in a complex matrix should be eliminated + * in spFactor() using direct addressing (rather than indirect + * addressing). + * DoRealDirect (SPBOOLEAN *) + * Array of flags, one for each column in matrix. If a flag is true + * then corresponding column in a real matrix should be eliminated + * in spFactor() using direct addressing (rather than indirect + * addressing). + * Elements (int) + * The number of original elements (total elements minus fill ins) + * present in matrix. + * Error (int) + * The error status of the sparse matrix package. + * ExtSize (int) + * The value of the largest external row or column number encountered. + * ExtToIntColMap (int []) + * An array that is used to convert external columns number to internal + * external column numbers. Present only if TRANSLATE option is set true. + * ExtToIntRowMap (int []) + * An array that is used to convert external row numbers to internal + * external row numbers. Present only if TRANSLATE option is set true. + * Factored (SPBOOLEAN) + * Indicates if matrix has been factored. This flag is set true in + * spFactor() and spOrderAndFactor() and set false in spCreate() + * and spClear(). + * Fillins (int) + * The number of fill-ins created during the factorization the matrix. + * FirstInCol (ArrayOfElementPtrs) + * Array of pointers that point to the first nonzero element of the + * column corresponding to the index. + * FirstInRow (ArrayOfElementPtrs) + * Array of pointers that point to the first nonzero element of the row + * corresponding to the index. + * ID (unsigned long int) + * A constant that provides the sparse data structure with a signature. + * When DEBUG is true, all externally available sparse routines check + * this signature to assure they are operating on a valid matrix. + * Intermediate (RealVector) + * Temporary storage used in the spSolve routines. Intermediate is an + * array used during forward and backward substitution. It is + * commonly called y when the forward and backward substitution process is + * denoted Ax = b => Ly = b and Ux = y. + * InternalVectorsAllocated (SPBOOLEAN) + * A flag that indicates whether the markowitz vectors and the + * Intermediate vector have been created. + * These vectors are created in CreateInternalVectors(). + * IntToExtColMap (int []) + * An array that is used to convert internal column numbers to external + * external column numbers. + * IntToExtRowMap (int []) + * An array that is used to convert internal row numbers to external + * external row numbers. + * MarkowitzCol (int []) + * An array that contains the count of the non-zero elements excluding + * the pivots for each column. Used to generate and update MarkowitzProd. + * MarkowitzProd (long []) + * The array of the products of the Markowitz row and column counts. The + * element with the smallest product is the best pivot to use to maintain + * sparsity. + * MarkowitzRow (int []) + * An array that contains the count of the non-zero elements excluding + * the pivots for each row. Used to generate and update MarkowitzProd. + * MaxRowCountInLowerTri (int) + * The maximum number of off-diagonal element in the rows of L, the + * lower triangular matrix. This quantity is used when computing an + * estimate of the roundoff error in the matrix. + * NeedsOrdering (SPBOOLEAN) + * This is a flag that signifies that the matrix needs to be ordered + * or reordered. NeedsOrdering is set true in spCreate() and + * spGetElement() or spGetAdmittance() if new elements are added to the + * matrix after it has been previously factored. It is set false in + * spOrderAndFactor(). + * NumberOfInterchangesIsOdd (SPBOOLEAN) + * Flag that indicates the sum of row and column interchange counts + * is an odd number. Used when determining the sign of the determinant. + * Partitioned (SPBOOLEAN) + * This flag indicates that the columns of the matrix have been + * partitioned into two groups. Those that will be addressed directly + * and those that will be addressed indirectly in spFactor(). + * PivotsOriginalCol (int) + * Column pivot was chosen from. + * PivotsOriginalRow (int) + * Row pivot was chosen from. + * PivotSelectionMethod (char) + * Character that indicates which pivot search method was successful. + * PreviousMatrixWasComplex (SPBOOLEAN) + * This flag in needed to determine how to clear the matrix. When + * dealing with real matrices, it is important that the imaginary terms + * in the matrix elements be zero. Thus, if the previous matrix was + * complex, then the current matrix will be cleared as if it were complex + * even if it is real. + * RelThreshold (RealNumber) + * The magnitude an element must have relative to others in its row + * to be considered as a pivot candidate, except as a last resort. + * Reordered (SPBOOLEAN) + * This flag signifies that the matrix has been reordered. It + * is cleared in spCreate(), set in spMNA_Preorder() and + * spOrderAndFactor() and is used in spPrint(). + * RowsLinked (SPBOOLEAN) + * A flag that indicates whether the row pointers exist. The AddByIndex + * routines do not generate the row pointers, which are needed by some + * of the other routines, such as spOrderAndFactor() and spScale(). + * The row pointers are generated in the function spcLinkRows(). + * SingularCol (int) + * Normally zero, but if matrix is found to be singular, SingularCol is + * assigned the external column number of pivot that was zero. + * SingularRow (int) + * Normally zero, but if matrix is found to be singular, SingularRow is + * assigned the external row number of pivot that was zero. + * Singletons (int) + * The number of singletons available for pivoting. Note that if row I + * and column I both contain singletons, only one of them is counted. + * Size (int) + * Number of rows and columns in the matrix. Does not change as matrix + * is factored. + * TrashCan (MatrixElement) + * This is a dummy MatrixElement that is used to by the user to stuff + * data related to the zero row or column. In other words, when the user + * adds an element in row zero or column zero, then the matrix returns + * a pointer to TrashCan. In this way the user can have a uniform way + * data into the matrix independent of whether a component is connected + * to ground. + * + * >>> The remaining fields are related to memory allocation. + * TopOfAllocationList (AllocationListPtr) + * Pointer which points to the top entry in a list. The list contains + * all the pointers to the segments of memory that have been allocated + * to this matrix. This is used when the memory is to be freed on + * deallocation of the matrix. + * RecordsRemaining (int) + * Number of slots left in the list of allocations. + * NextAvailElement (ElementPtr) + * Pointer to the next available element which has been allocated but as + * yet is unused. Matrix elements are allocated in groups of + * ELEMENTS_PER_ALLOCATION in order to speed element allocation and + * freeing. + * ElementsRemaining (int) + * Number of unused elements left in last block of elements allocated. + * NextAvailFillin (ElementPtr) + * Pointer to the next available fill-in which has been allocated but + * as yet is unused. Fill-ins are allocated in a group in order to keep + * them physically close in memory to the rest of the matrix. + * FillinsRemaining (int) + * Number of unused fill-ins left in the last block of fill-ins + * allocated. + * FirstFillinListNode (FillinListNodeStruct *) + * A pointer to the head of the linked-list that keeps track of the + * lists of fill-ins. + * LastFillinListNode (FillinListNodeStruct *) + * A pointer to the tail of the linked-list that keeps track of the + * lists of fill-ins. + */ + +/* Begin `MatrixFrame'. */ +struct MatrixFrame +{ + int NumRank ; /* the numerical Rank of the matrix */ + RealNumber AbsThreshold; + int AllocatedSize; + int AllocatedExtSize; + SPBOOLEAN Complex; + int CurrentSize; + ArrayOfElementPtrs Diag; + SPBOOLEAN *DoCmplxDirect; + SPBOOLEAN *DoRealDirect; + int Elements; + int Error; + int ExtSize; + int *ExtToIntColMap; + int *ExtToIntRowMap; + SPBOOLEAN Factored; + int Fillins; + ArrayOfElementPtrs FirstInCol; + ArrayOfElementPtrs FirstInRow; + unsigned long ID; + RealVector Intermediate; + SPBOOLEAN InternalVectorsAllocated; + int *IntToExtColMap; + int *IntToExtRowMap; + int *MarkowitzRow; + int *MarkowitzCol; + long *MarkowitzProd; + int MaxRowCountInLowerTri; + SPBOOLEAN NeedsOrdering; + SPBOOLEAN NumberOfInterchangesIsOdd; + SPBOOLEAN Partitioned; + int PivotsOriginalCol; + int PivotsOriginalRow; + char PivotSelectionMethod; + SPBOOLEAN PreviousMatrixWasComplex; + RealNumber RelThreshold; + SPBOOLEAN Reordered; + SPBOOLEAN RowsLinked; + int SingularCol; + int SingularRow; + int Singletons; + int Size; + struct MatrixElement TrashCan; + + AllocationListPtr TopOfAllocationList; + int RecordsRemaining; + ElementPtr NextAvailElement; + int ElementsRemaining; + ElementPtr NextAvailFillin; + int FillinsRemaining; + struct FillinListNodeStruct *FirstFillinListNode; + struct FillinListNodeStruct *LastFillinListNode; +}; +typedef struct MatrixFrame *MatrixPtr; + +#endif /* __SPDEFS_H__*/ diff --git a/modules/sparse/includes/spmalloc.h b/modules/sparse/includes/spmalloc.h new file mode 100755 index 000000000..1ac571ce0 --- /dev/null +++ b/modules/sparse/includes/spmalloc.h @@ -0,0 +1,33 @@ +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) INRIA + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + * + */ + +#ifndef __SPMALLOC_H__ +#define __SPMALLOC_H__ + +#include <stdlib.h> +#include "MALLOC.h" +#define SPMALLOC(x) MALLOC(((size_t) x)) +#define SPALLOC(type,number) ((type *)SPMALLOC((unsigned)(sizeof(type)*(number)))) +#define SPREALLOC(ptr,type,number) \ + ptr = (type *)REALLOC((char *)ptr,(unsigned)(sizeof(type)*(number))) + +#define SPFREE(ptr) { if ((ptr) != NULL) {FREE((void *)(ptr)); (ptr) = NULL;}} + + +/* Calloc that properly handles allocating a cleared vector. */ +#define SPCALLOC(ptr,type,number) \ +{ int i; ptr = SPALLOC(type, number); \ + if (ptr != (type *)NULL) \ + for(i=(number)-1;i>=0; i--) ptr[i] = (type) 0; \ +} + +#endif /*__SPMALLOC_H__*/ diff --git a/modules/sparse/includes/spmatrix.h b/modules/sparse/includes/spmatrix.h new file mode 100755 index 000000000..6b1e1fc35 --- /dev/null +++ b/modules/sparse/includes/spmatrix.h @@ -0,0 +1,325 @@ +#ifndef __SPMATRIX_H__ +#define __SPMATRIX_H__ + +/* + * EXPORTS for sparse matrix routines. + * + * Author: Advising professor: + * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli + * UC Berkeley + * + * This file contains definitions that are useful to the calling + * program. In particular, this file contains error keyword + * definitions, some macro functions that are used to quickly enter + * data into the matrix and the type definition of a data structure + * that acts as a template for entering admittances into the matrix. + * Also included is the type definitions for the various functions + * available to the user. + */ + + +/* + * Revision and copyright information. + * + * Copyright (c) 1985,86,87,88 + * by Kenneth S. Kundert and the University of California. + * + * Permission to use, copy, modify, and distribute this software and + * its documentation for any purpose and without fee is hereby granted, + * provided that the copyright notices appear in all copies and + * supporting documentation and that the authors and the University of + * California are properly credited. The authors and the University of + * California make no representations as to the suitability of this + * software for any purpose. It is provided `as is', without express + * or implied warranty. + * + */ + + + + +#ifndef spOKAY + +/* + * IMPORTS + * + * >>> Import descriptions: + * spConfig.h + * Macros that customize the sparse matrix routines. + */ + +#include "spConfig.h" + + + + + + +/* + * ERROR KEYWORDS + * + * The actual numbers used in the error codes are not sacred, they can be + * changed under the condition that the codes for the nonfatal errors are + * less than the code for spFATAL and similarly the codes for the fatal + * errors are greater than that for spFATAL. + * + * >>> Error descriptions: + * spOKAY + * No error has occurred. + * spSMALL_PIVOT + * When reordering the matrix, no element was found which satisfies the + * absolute threshold criteria. The largest element in the matrix was + * chosen as pivot. Non-fatal. + * spZERO_DIAG + * Fatal error. A zero was encountered on the diagonal the matrix. This + * does not necessarily imply that the matrix is singular. When this + * error occurs, the matrix should be reconstructed and factored using + * spOrderAndFactor(). + * spSINGULAR + * Fatal error. Matrix is singular, so no unique solution exists. + * spNO_MEMORY + * Fatal error. Indicates that not enough memory is available to handle + * the matrix. + * spPANIC + * Fatal error indicating that the routines are not prepared to + * handle the matrix that has been requested. This may occur when + * the matrix is specified to be real and the routines are not + * compiled for real matrices, or when the matrix is specified to + * be complex and the routines are not compiled to handle complex + * matrices. + * spFATAL + * Not an error flag, but rather the dividing line between fatal errors + * and warnings. + */ + +/* Begin error macros. */ +#define spOKAY 0 +#define spSMALL_PIVOT 1 +#define spZERO_DIAG 2 +#define spSINGULAR 3 +#define spNO_MEMORY 4 +#define spPANIC 5 + +#define spFATAL 2 + + + + + + + +/* + * KEYWORD DEFINITIONS + * + * Here we define what precision arithmetic Sparse will use. Double + * precision is suggested as being most appropriate for circuit + * simulation and for C. However, it is possible to change spREAL + * to a float for single precision arithmetic. Note that in C, single + * precision arithmetic is often slower than double precision. Sparse + * internally refers to spREALs as RealNumbers. + * + * Some C compilers, notably the old VMS compiler, do not handle the keyword + * "void" correctly. If this is true for your compiler, remove the + * comment delimiters from the redefinition of void to int below. + */ + +#define spREAL double +/* #define void int */ + + + + + +/* + * PARTITION TYPES + * + * When factoring a previously ordered matrix using spFactor(), Sparse + * operates on a row-at-a-time basis. For speed, on each step, the row + * being updated is copied into a full vector and the operations are + * performed on that vector. This can be done one of two ways, either + * using direct addressing or indirect addressing. Direct addressing + * is fastest when the matrix is relatively dense and indirect addressing + * is quite sparse. The user can select which partitioning mode is used. + * The following keywords are passed to spPartition() and indicate that + * Sparse should use only direct addressing, only indirect addressing, or + * that it should choose the best mode on a row-by-row basis. The time + * required to choose a partition is of the same order of the cost to factor + * the matrix. + * + * If you plan to factor a large number of matrices with the same structure, + * it is best to let Sparse choose the partition. Otherwise, you should + * choose the partition based on the predicted density of the matrix. + */ + +/* Begin partition keywords. */ + +#define spDEFAULT_PARTITION 0 +#define spDIRECT_PARTITION 1 +#define spINDIRECT_PARTITION 2 +#define spAUTO_PARTITION 3 + + + + + +/* + * MACRO FUNCTION DEFINITIONS + * + * >>> Macro descriptions: + * spADD_REAL_ELEMENT + * Macro function that adds data to a real element in the matrix by a + * pointer. + * spADD_IMAG_ELEMENT + * Macro function that adds data to a imaginary element in the matrix by + * a pointer. + * spADD_COMPLEX_ELEMENT + * Macro function that adds data to a complex element in the matrix by a + * pointer. + * spADD_REAL_QUAD + * Macro function that adds data to each of the four real matrix elements + * specified by the given template. + * spADD_IMAG_QUAD + * Macro function that adds data to each of the four imaginary matrix + * elements specified by the given template. + * spADD_COMPLEX_QUAD + * Macro function that adds data to each of the four complex matrix + * elements specified by the given template. + */ + +/* Begin Macros. */ +#define spADD_REAL_ELEMENT(element,real) *(element) += real + +#define spADD_IMAG_ELEMENT(element,imag) *(element+1) += imag + +#define spADD_COMPLEX_ELEMENT(element,real,imag) \ +{ *(element) += real; \ + *(element+1) += imag; \ +} + +#define spADD_REAL_QUAD(template,real) \ +{ *((template).Element1) += real; \ + *((template).Element2) += real; \ + *((template).Element3Negated) -= real; \ + *((template).Element4Negated) -= real; \ +} + +#define spADD_IMAG_QUAD(template,imag) \ +{ *((template).Element1+1) += imag; \ + *((template).Element2+1) += imag; \ + *((template).Element3Negated+1) -= imag; \ + *((template).Element4Negated+1) -= imag; \ +} + +#define spADD_COMPLEX_QUAD(template,real,imag) \ +{ *((template).Element1) += real; \ + *((template).Element2) += real; \ + *((template).Element3Negated) -= real; \ + *((template).Element4Negated) -= real; \ + *((template).Element1+1) += imag; \ + *((template).Element2+1) += imag; \ + *((template).Element3Negated+1) -= imag; \ + *((template).Element4Negated+1) -= imag; \ +} + + + + + + + +/* + * TYPE DEFINITION FOR COMPONENT TEMPLATE + * + * This data structure is used to hold pointers to four related elements in + * matrix. It is used in conjunction with the routines + * spGetAdmittance + * spGetQuad + * spGetOnes + * These routines stuff the structure which is later used by the spADD_QUAD + * macro functions above. It is also possible for the user to collect four + * pointers returned by spGetElement and stuff them into the template. + * The spADD_QUAD routines stuff data into the matrix in locations specified + * by Element1 and Element2 without changing the data. The data is negated + * before being placed in Element3 and Element4. + */ + +/* Begin `spTemplate'. */ +struct spTemplate +{ + spREAL *Element1 ; + spREAL *Element2 ; + spREAL *Element3Negated; + spREAL *Element4Negated; +}; + + + + + +/* + * FUNCTION TYPE DEFINITIONS + * + * The type of every user accessible function is declared here. + */ + +/* Begin function declarations. */ + +/* For compilers that understand function prototypes. */ + +extern void spClear( char* ); +extern spREAL spCondition( char*, spREAL, int* ); +extern char *spCreate( int, int, int* ); +extern void spDeleteRowAndCol( char*, int, int ); +extern void spDestroy( char* ); +extern int spElementCount( char* ); +extern int spError( char* ); +extern int spFactor( char* ); +extern int spFileMatrix( char*, char*, char*, int, int, int ); +extern int spFileStats( char*, char*, char* ); +extern int spFillinCount( char* ); +extern int spGetAdmittance( char*, int, int, struct spTemplate* ); +extern spREAL *spGetElement( char*, int, int ); +extern char *spGetInitInfo( spREAL* ); +extern int spGetOnes( char*, int, int, int, struct spTemplate* ); +extern int spGetQuad( char*, int, int, int, int, struct spTemplate* ); +extern int spGetSize( char*, int ); +extern int spInitialize( char*, int (*)() ); +extern void spInstallInitInfo( spREAL*, char* ); +extern spREAL spLargestElement( char* ); +extern void spMNA_Preorder( char* ); +extern spREAL spNorm( char* ); +extern int spOrderAndFactor( char*, spREAL[], spREAL, spREAL, int ); +extern void spPartition( char*, int ); +extern void spPrint( char*, int, int, int ); +extern spREAL spPseudoCondition( char* ); +extern spREAL spRoundoff( char*, spREAL ); +extern void spScale( char*, spREAL[], spREAL[] ); +extern void spSetComplex( char* ); +extern void spSetReal( char* ); +extern void spStripFills( char* ); +extern void spWhereSingular( char*, int*, int* ); + +/* Functions with argument lists that are dependent on options. */ + +#if spCOMPLEX +extern void spDeterminant ( char*, int*, spREAL*, spREAL* ); +#else /* NOT spCOMPLEX */ +extern void spDeterminant ( char*, int*, spREAL* ); +#endif /* NOT spCOMPLEX */ +#if spCOMPLEX && spSEPARATED_COMPLEX_VECTORS +extern int spFileVector( char*, char* , spREAL[], spREAL[]); +extern void spMultiply( char*, spREAL[], spREAL[], spREAL[], spREAL[] ); +extern void spMultTransposed(char*, spREAL[], spREAL[], spREAL[], spREAL[]); +extern void spSolve( char*, spREAL[], spREAL[], spREAL[], spREAL[] ); +extern void spSolveTransposed(char*, spREAL[], spREAL[], spREAL[], spREAL[]); +#else /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */ +extern int spFileVector( char*, char* , spREAL[] ); +extern void spMultiply( char*, spREAL[], spREAL[] ); +extern void spMultTransposed( char*, spREAL[], spREAL[] ); +extern void spSolve( char*, spREAL[], spREAL[] ); +extern void spSolveTransposed( char*, spREAL[], spREAL[] ); +#endif /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */ + +#endif /* spOKAY */ + +#endif /* __SPMATRIX_H__ */ |