# This is a shell archive. Remove anything before this line, # then unpack it by saving it in a file and typing "sh file". # # Wrapped by kundert at athos on Fri Jul 1 01:56:16 1988 # # This archive contains: # MAKE.COM spFortran.c spOutput.c spRevision # spTest.c spUtils.c mat0 mat1 # mat2 mat3 mat5 # LANG=""; export LANG echo x - MAKE.COM cat >MAKE.COM <<'@EOF' $ def lnk$library sys$library:vaxcrtl.olb $ cc spAllocate.c $ cc spBuild.c $ cc spFactor.c $ cc spFortran.c $ cc spOutput.c $ cc spSolve.c $ cc spTest.c $ cc spUtils.c $ library/create/object sparse.olb - spAllocate, - spBuild, - spFactor, - spFortran, - spOutput, - spSolve, - spUtils $ link/exe=sparse spTest.obj, sparse.olb/library $ sparse :== $$disk2:[kundert.sparse]sparse.exe @EOF chmod 644 MAKE.COM echo x - spFortran.c cat >spFortran.c <<'@EOF' /* * SPARSE FORTRAN MODULE * * Author: Advising professor: * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli * UC Berkeley * * This module contains routines that interface Sparse1.3 to a calling * program written in fortran. Almost every externally available Sparse1.3 * routine has a counterpart defined in this file, with the name the * same except the `sp' prefix is changed to `sf'. The spADD_ELEMENT * and spADD_QUAD macros are also replaced with the sfAdd1 and sfAdd4 * functions defined in this file. * * To ease porting this file to different operating systems, the names of * the functions can be easily redefined (search for `Routine Renaming'). * A simple example of a FORTRAN program that calls Sparse is included in * this file (search for Example). When interfacing to a FORTRAN program, * the ARRAY_OFFSET option should be set to NO (see spConfig.h). * * DISCLAIMER: * These interface routines were written by a C programmer who has little * experience with FORTRAN. The routines have had minimal testing. * Any interface between two languages is going to have portability * problems, this one is no exception. * * * >>> User accessible functions contained in this file: * sfCreate() * sfDestroy() * sfStripFills() * sfClear() * sfGetElement() * sfGetAdmittance() * sfGetQuad() * sfGetOnes() * sfAdd1Real() * sfAdd1Imag() * sfAdd1Complex() * sfAdd4Real() * sfAdd4Imag() * sfAdd4Complex() * sfOrderAndFactor() * sfFactor() * sfPartition() * sfSolve() * sfSolveTransposed() * sfPrint() * sfFileMatrix() * sfFileVector() * sfFileStats() * sfMNA_Preorder() * sfScale() * sfMultiply() * sfTransMultiply() * sfDeterminant() * sfError() * sfWhereSingular() * sfGetSize() * sfSetReal() * sfSetComplex() * sfFillinCount() * sfElementCount() * sfDeleteRowAndCol() * sfPseudoCondition() * sfCondition() * sfNorm() * sfLargestElement() * sfRoundoff() */ /* * FORTRAN -- C COMPATIBILITY * * Fortran and C data types correspond in the following way: * -- C -- -- FORTRAN -- * int INTEGER*4 or INTEGER*2 (machine dependent, usually int*4) * long INTEGER*4 * float REAL * double DOUBLE PRECISION (used by default in preference to float) * * The complex number format used by Sparse is compatible with that * used by FORTRAN. C pointers are passed to FORTRAN as longs, they should * not be used in any way in FORTRAN. */ /* * 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 above copyright notice 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 lint static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; static char RCSid[] = "@(#)$Header: spFortran.c,v 1.1 88/06/18 11:15:30 kundert Exp $"; #endif /* * IMPORTS * * >>> Import descriptions: * spConfig.h * Macros that customize the sparse matrix routines. * spMatrix.h * Macros and declarations to be imported by the user. * spDefs.h * Matrix type and macro definitions for the sparse matrix routines. */ #define spINSIDE_SPARSE #include "spConfig.h" #include "spMatrix.h" #include "spDefs.h" #if FORTRAN /* * Routine Renaming */ #define sfCreate sfcreate #define sfStripFills sfstripfills #define sfDestroy sfdestroy #define sfClear sfzero #define sfGetElement sfgetelement #define sfGetAdmittance sfgetadmittance #define sfGetQuad sfgetquad #define sfGetOnes sfgetones #define sfAdd1Real sfadd1real #define sfAdd1Imag sfadd1imag #define sfAdd1Complex sfadd1complex #define sfAdd4Real sfadd4real #define sfAdd4Imag sfadd4imag #define sfAdd4Complex sfadd4complex #define sfOrderAndFactor sforderandfactor #define sfFactor sffactor #define sfPartition sfpartition #define sfSolve sfsolve #define sfSolveTransposed sfsolvetransposed #define sfPrint sfprint #define sfFileMatrix sffilematrix #define sfFileVector sffilevector #define sfFileStats sffilestats #define sfMNA_Preorder sfmna_preorder #define sfScale sfscale #define sfMultiply sfmultiply #define sfDeterminant sfdeterminant #define sfError sferror #define sfWhereSingular sfwheresingular #define sfGetSize sfgetsize #define sfSetReal sfsetreal #define sfSetComplex sfsetcomplex #define sfFillinCount sffillincount #define sfElementCount sfelementcount #define sfDeleteRowAndCol sfdeleterowandcol #define sfPseudoCondition sfpseudocondition #define sfCondition sfcondition #define sfNorm sfnorm #define sfLargestElement sflargestelement #define sfRoundoff sfroundoff /* * Example of a FORTRAN Program Calling Sparse * integer matrix, error, sfCreate, sfGetElement, spFactor integer element(10) double precision rhs(4), solution(4) c matrix = sfCreate(4,0,error) element(1) = sfGetElement(matrix,1,1) element(2) = sfGetElement(matrix,1,2) element(3) = sfGetElement(matrix,2,1) element(4) = sfGetElement(matrix,2,2) element(5) = sfGetElement(matrix,2,3) element(6) = sfGetElement(matrix,3,2) element(7) = sfGetElement(matrix,3,3) element(8) = sfGetElement(matrix,3,4) element(9) = sfGetElement(matrix,4,3) element(10) = sfGetElement(matrix,4,4) call sfClear(matrix) call sfAdd1Real(element(1), 2d0) call sfAdd1Real(element(2), -1d0) call sfAdd1Real(element(3), -1d0) call sfAdd1Real(element(4), 3d0) call sfAdd1Real(element(5), -1d0) call sfAdd1Real(element(6), -1d0) call sfAdd1Real(element(7), 3d0) call sfAdd1Real(element(8), -1d0) call sfAdd1Real(element(9), -1d0) call sfAdd1Real(element(10), 3d0) call sfprint(matrix, .false., .false.) rhs(1) = 34d0 rhs(2) = 0d0 rhs(3) = 0d0 rhs(4) = 0d0 error = sfFactor(matrix) call sfSolve(matrix, rhs, solution) write (6, 10) rhs(1), rhs(2), rhs(3), rhs(4) 10 format (f 10.2) end * */ /* * MATRIX ALLOCATION * * Allocates and initializes the data structures associated with a matrix. * * >>> Returned: [INTEGER] * A pointer to the matrix is returned cast into an integer. This pointer * is then passed and used by the other matrix routines to refer to a * particular matrix. If an error occurs, the NULL pointer is returned. * * >>> Arguments: * Size (long *) [INTEGER] * Size of matrix or estimate of size of matrix if matrix is EXPANDABLE. * Complex (int *) [INTEGER or INTEGER*2] * Type of matrix. If Complex is 0 then the matrix is real, otherwise * the matrix will be complex. Note that if the routines are not set up * to handle the type of matrix requested, then a spPANIC error will occur. * Further note that if a matrix will be both real and complex, it must * be specified here as being complex. * Error (int *) [INTEGER or INTEGER*2] * Returns error flag, needed because function spError() will not work * correctly if spCreate() returns NULL. * * >>> Possible errors: * spNO_MEMORY * spPANIC * Error is cleared in this routine. */ long sfCreate( Size, Complex, Error ) int *Size, *Complex, *Error; { /* Begin `sfCreate'. */ return (long)spCreate(*Size, *Complex, Error ); } /* * MATRIX DEALLOCATION * * Deallocates pointers and elements of Matrix. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix frame which is to be removed from memory. */ void sfDestroy( Matrix ) long *Matrix; { /* Begin `sfDestroy'. */ spDestroy((char *)*Matrix); return; } #if STRIP /* * STRIP FILL-INS FROM MATRIX * * Strips the matrix of all fill-ins. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix to be stripped. */ void sfStripFills( Matrix ) long *Matrix; { /* Begin `sfStripFills'. */ spStripFills((char *)*Matrix); return; } #endif /* * CLEAR MATRIX * * Sets every element of the matrix to zero and clears the error flag. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix that is to be cleared. */ void sfClear( Matrix ) long *Matrix; { /* Begin `sfClear'. */ spClear((char *)*Matrix); return; } /* * SINGLE ELEMENT ADDITION TO MATRIX BY INDEX * * Finds element [Row,Col] and returns a pointer to it. If element is * not found then it is created and spliced into matrix. This routine * is only to be used after spCreate() and before spMNA_Preorder(), * spFactor() or spOrderAndFactor(). Returns a pointer to the * Real portion of a MatrixElement. This pointer is later used by * sfAddxxxxx() to directly access element. * * >>> Returns: [INTEGER] * Returns a pointer to the element. This pointer is then used to directly * access the element during successive builds. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix that the element is to be added to. * Row (int *) [INTEGER or INTEGER*2] * Row index for element. Must be in the range of [0..Size] unless * the options EXPANDABLE or TRANSLATE are used. Elements placed in * row zero are discarded. In no case may Row be less than zero. * Col (int *) [INTEGER or INTEGER*2] * Column index for element. Must be in the range of [0..Size] unless * the options EXPANDABLE or TRANSLATE are used. Elements placed in * column zero are discarded. In no case may Col be less than zero. * * >>> Possible errors: * spNO_MEMORY * Error is not cleared in this routine. */ long sfGetElement( Matrix, Row, Col ) long *Matrix; int *Row, *Col; { /* Begin `sfGetElement'. */ return (long)spGetElement((char *)*Matrix, *Row, *Col); } #if QUAD_ELEMENT /* * ADDITION OF ADMITTANCE TO MATRIX BY INDEX * * Performs same function as sfGetElement except rather than one * element, all four Matrix elements for a floating component are * added. This routine also works if component is grounded. Positive * elements are placed at [Node1,Node2] and [Node2,Node1]. This * routine is only to be used after sfCreate() and before * sfMNA_Preorder(), sfFactor() or sfOrderAndFactor(). * * >>> Returns: [INTEGER or INTEGER*2] * Error code. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix that component is to be entered in. * Node1 (int *) [INTEGER or INTEGER*2] * Row and column indices for elements. Must be in the range of [0..Size] * unless the options EXPANDABLE or TRANSLATE are used. Node zero is the * ground node. In no case may Node1 be less than zero. * Node2 (int *) [INTEGER or INTEGER*2] * Row and column indices for elements. Must be in the range of [0..Size] * unless the options EXPANDABLE or TRANSLATE are used. Node zero is the * ground node. In no case may Node2 be less than zero. * Template (long[4]) [INTEGER (4)] * Collection of pointers to four elements that are later used to directly * address elements. User must supply the template, this routine will * fill it. * * Possible errors: * spNO_MEMORY * Error is not cleared in this routine. */ int sfGetAdmittance( Matrix, Node1, Node2, Template ) long *Matrix, Template[4]; int *Node1, *Node2; { /* Begin `spGetAdmittance'. */ return ( spGetAdmittance((char *)*Matrix, *Node1, *Node2, (struct spTemplate *)Template ) ); } #endif /* QUAD_ELEMENT */ #if QUAD_ELEMENT /* * ADDITION OF FOUR ELEMENTS TO MATRIX BY INDEX * * Similar to sfGetAdmittance, except that sfGetAdmittance only * handles 2-terminal components, whereas sfGetQuad handles simple * 4-terminals as well. These 4-terminals are simply generalized * 2-terminals with the option of having the sense terminals different * from the source and sink terminals. sfGetQuad adds four * elements to the matrix. Positive elements occur at Row1,Col1 * Row2,Col2 while negative elements occur at Row1,Col2 and Row2,Col1. * The routine works fine if any of the rows and columns are zero. * This routine is only to be used after sfCreate() and before * sfMNA_Preorder(), sfFactor() or sfOrderAndFactor() * unless TRANSLATE is set true. * * >>> Returns: [INTEGER or INTEGER*2] * Error code. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix that component is to be entered in. * Row1 (int *) [INTEGER or INTEGER*2] * First row index for elements. Must be in the range of [0..Size] * unless the options EXPANDABLE or TRANSLATE are used. Zero is the * ground row. In no case may Row1 be less than zero. * Row2 (int *) [INTEGER or INTEGER*2] * Second row index for elements. Must be in the range of [0..Size] * unless the options EXPANDABLE or TRANSLATE are used. Zero is the * ground row. In no case may Row2 be less than zero. * Col1 (int *) [INTEGER or INTEGER*2] * First column index for elements. Must be in the range of [0..Size] * unless the options EXPANDABLE or TRANSLATE are used. Zero is the * ground column. In no case may Col1 be less than zero. * Col2 (int *) [INTEGER or INTEGER*2] * Second column index for elements. Must be in the range of [0..Size] * unless the options EXPANDABLE or TRANSLATE are used. Zero is the * ground column. In no case may Col2 be less than zero. * Template (long[4]) [INTEGER (4)] * Collection of pointers to four elements that are later used to directly * address elements. User must supply the template, this routine will * fill it. * * Possible errors: * spNO_MEMORY * Error is not cleared in this routine. */ int sfGetQuad( Matrix, Row1, Row2, Col1, Col2, Template ) long *Matrix, Template[4]; int *Row1, *Row2, *Col1, *Col2; { /* Begin `spGetQuad'. */ return ( spGetQuad( (char *)*Matrix, *Row1, *Row2, *Col1, *Col2, (struct spTemplate *)Template ) ); } #endif /* QUAD_ELEMENT */ #if QUAD_ELEMENT /* * ADDITION OF FOUR STRUCTURAL ONES TO MATRIX BY INDEX * * Performs similar function to sfGetQuad() except this routine is * meant for components that do not have an admittance representation. * * The following stamp is used: * Pos Neg Eqn * Pos [ . . 1 ] * Neg [ . . -1 ] * Eqn [ 1 -1 . ] * * >>> Returns: [INTEGER or INTEGER*2] * Error code. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix that component is to be entered in. * Pos (int *) [INTEGER or INTEGER*2] * See stamp above. Must be in the range of [0..Size] * unless the options EXPANDABLE or TRANSLATE are used. Zero is the * ground row. In no case may Pos be less than zero. * Neg (int *) [INTEGER or INTEGER*2] * See stamp above. Must be in the range of [0..Size] * unless the options EXPANDABLE or TRANSLATE are used. Zero is the * ground row. In no case may Neg be less than zero. * Eqn (int *) [INTEGER or INTEGER*2] * See stamp above. Must be in the range of [0..Size] * unless the options EXPANDABLE or TRANSLATE are used. Zero is the * ground row. In no case may Eqn be less than zero. * Template (long[4]) [INTEGER (4)] * Collection of pointers to four elements that are later used to directly * address elements. User must supply the template, this routine will * fill it. * * Possible errors: * spNO_MEMORY * Error is not cleared in this routine. */ int sfGetOnes(Matrix, Pos, Neg, Eqn, Template) long *Matrix, Template[4]; int *Pos, *Neg, *Eqn; { /* Begin `sfGetOnes'. */ return ( spGetOnes( (char *)*Matrix, *Pos, *Neg, *Eqn, (struct spTemplate *)Template ) ); } #endif /* QUAD_ELEMENT */ /* * ADD ELEMENT(S) DIRECTLY TO MATRIX * * Adds a value to an element or a set of four element in a matrix. * These elements are referenced by pointer, and so must already have * been created by spGetElement(), spGetAdmittance(), spGetQuad(), or * spGetOnes(). * * >>> Arguments: * Element (long *) [INTEGER] * Pointer to the element that is to be added to. * Template (long[4]) [INTEGER (4)] * Pointer to the element that is to be added to. * Real (spREAL *) [REAL or DOUBLE PRECISION] * Real portion of the number to be added to the element. * Imag (spREAL *) [REAL or DOUBLE PRECISION] * Imaginary portion of the number to be added to the element. */ void sfAdd1Real( Element, Real ) long *Element; RealNumber *Real; { /* Begin `sfAdd1Real'. */ *((RealNumber *)*Element) += *Real; } #if spCOMPLEX void sfAdd1Imag( Element, Imag ) long *Element; RealNumber *Imag; { /* Begin `sfAdd1Imag'. */ *(((RealNumber *)*Element)+1) += *Imag; } void sfAdd1Complex( Element, Real, Imag ) long *Element; RealNumber *Real, *Imag; { /* Begin `sfAdd1Complex'. */ *((RealNumber *)*Element) += *Real; *(((RealNumber *)*Element)+1) += *Imag; } #endif /* spCOMPLEX */ #if QUAD_ELEMENT void sfAdd4Real( Template, Real ) long Template[4]; RealNumber *Real; { /* Begin `sfAdd4Real'. */ *((RealNumber *)Template[0]) += *Real; *((RealNumber *)Template[1]) += *Real; *((RealNumber *)Template[2]) -= *Real; *((RealNumber *)Template[3]) -= *Real; } #if spCOMPLEX void sfAdd4Imag( Template, Imag ) long Template[4]; RealNumber *Imag; { /* Begin `sfAdd4Imag'. */ *(((RealNumber *)Template[0])+1) += *Imag; *(((RealNumber *)Template[1])+1) += *Imag; *(((RealNumber *)Template[2])+1) -= *Imag; *(((RealNumber *)Template[3])+1) -= *Imag; } void sfAdd4Complex( Template, Real, Imag ) long Template[4]; RealNumber *Real, *Imag; { /* Begin `sfAdd4Complex'. */ *((RealNumber *)Template[0]) += *Real; *((RealNumber *)Template[1]) += *Real; *((RealNumber *)Template[2]) -= *Real; *((RealNumber *)Template[3]) -= *Real; *(((RealNumber *)Template[0])+1) += *Imag; *(((RealNumber *)Template[1])+1) += *Imag; *(((RealNumber *)Template[2])+1) -= *Imag; *(((RealNumber *)Template[3])+1) -= *Imag; } #endif /* spCOMPLEX */ #endif /* QUAD_ELEMENT */ /* * ORDER AND FACTOR MATRIX * * This routine chooses a pivot order for the matrix and factors it * into LU form. It handles both the initial factorization and subsequent * factorizations when a reordering is desired. This is handled in a manner * that is transparent to the user. The routine uses a variation of * Gauss's method where the pivots are associated with L and the * diagonal terms of U are one. * * >>> Returned: [INTEGER of INTEGER*2] * The error code is returned. Possible errors are listed below. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix. * RHS (RealVector) [REAL (1) or DOUBLE PRECISION (1)] * Representative right-hand side vector that is used to determine * pivoting order when the right hand side vector is sparse. If * RHS is a NULL pointer then the RHS vector is assumed to * be full and it is not used when determining the pivoting * order. * RelThreshold (RealNumber *) [REAL or DOUBLE PRECISION] * This number determines what the pivot relative threshold will * be. It should be between zero and one. If it is one then the * pivoting method becomes complete pivoting, which is very slow * and tends to fill up the matrix. If it is set close to zero * the pivoting method becomes strict Markowitz with no * threshold. The pivot threshold is used to eliminate pivot * candidates that would cause excessive element growth if they * were used. Element growth is the cause of roundoff error. * Element growth occurs even in well-conditioned matrices. * Setting the RelThreshold large will reduce element growth and * roundoff error, but setting it too large will cause execution * time to be excessive and will result in a large number of * fill-ins. If this occurs, accuracy can actually be degraded * because of the large number of operations required on the * matrix due to the large number of fill-ins. A good value seems * to be 0.001. The default is chosen by giving a value larger * than one or less than or equal to zero. This value should be * increased and the matrix resolved if growth is found to be * excessive. Changing the pivot threshold does not improve * performance on matrices where growth is low, as is often the * case with ill-conditioned matrices. Once a valid threshold is * given, it becomes the new default. The default value of * RelThreshold was choosen for use with nearly diagonally * dominant matrices such as node- and modified-node admittance * matrices. For these matrices it is usually best to use * diagonal pivoting. For matrices without a strong diagonal, it * is usually best to use a larger threshold, such as 0.01 or * 0.1. * AbsThreshold (RealNumber *) [REAL or DOUBLE PRECISION] * The absolute magnitude an element must have to be considered * as a pivot candidate, except as a last resort. This number * should be set significantly smaller than the smallest diagonal * element that is is expected to be placed in the matrix. If * there is no reasonable prediction for the lower bound on these * elements, then AbsThreshold should be set to zero. * AbsThreshold is used to reduce the possibility of choosing as a * pivot an element that has suffered heavy cancellation and as a * result mainly consists of roundoff error. Once a valid * threshold is given, it becomes the new default. * DiagPivoting (long *) [LOGICAL] * A flag indicating that pivot selection should be confined to the * diagonal if possible. If DiagPivoting is nonzero and if * DIAGONAL_PIVOTING is enabled pivots will be chosen only from * the diagonal unless there are no diagonal elements that satisfy * the threshold criteria. Otherwise, the entire reduced * submatrix is searched when looking for a pivot. The diagonal * pivoting in Sparse is efficient and well refined, while the * off-diagonal pivoting is not. For symmetric and near symmetric * matrices, it is best to use diagonal pivoting because it * results in the best performance when reordering the matrix and * when factoring the matrix without ordering. If there is a * considerable amount of nonsymmetry in the matrix, then * off-diagonal pivoting may result in a better equation ordering * simply because there are more pivot candidates to choose from. * A better ordering results in faster subsequent factorizations. * However, the initial pivot selection process takes considerably * longer for off-diagonal pivoting. * * >>> Possible errors: * spNO_MEMORY * spSINGULAR * spSMALL_PIVOT * Error is cleared in this function. */ int sfOrderAndFactor( Matrix, RHS, RelThreshold, AbsThreshold, DiagPivoting ) long *Matrix, *DiagPivoting; RealNumber RHS[], *RelThreshold, *AbsThreshold; { /* Begin `sfOrderAndFactor'. */ return spOrderAndFactor( (char *)*Matrix, RHS, *RelThreshold, *AbsThreshold, (BOOLEAN)*DiagPivoting ); } /* * FACTOR MATRIX * * This routine is the companion routine to spOrderAndFactor(). * Unlike sfOrderAndFactor(), sfFactor() cannot change the ordering. * It is also faster than sfOrderAndFactor(). The standard way of * using these two routines is to first use sfOrderAndFactor() for the * initial factorization. For subsequent factorizations, sfFactor() * is used if there is some assurance that little growth will occur * (say for example, that the matrix is diagonally dominant). If * sfFactor() is called for the initial factorization of the matrix, * then sfOrderAndFactor() is automatically called with the default * threshold. This routine uses "row at a time" LU factorization. * Pivots are associated with the lower triangular matrix and the * diagonals of the upper triangular matrix are ones. * * >>> Returned: [INTEGER or INTEGER*2] * The error code is returned. Possible errors are listed below. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix. * * >>> Possible errors: * spNO_MEMORY * spSINGULAR * spZERO_DIAG * spSMALL_PIVOT * Error is cleared in this function. */ int sfFactor( Matrix ) long *Matrix; { /* Begin `sfFactor'. */ return spFactor((char *)*Matrix); } /* * PARTITION MATRIX * * This routine determines the cost to factor each row using both * direct and indirect addressing and decides, on a row-by-row basis, * which addressing mode is fastest. This information is used in * sfFactor() to speed the factorization. * * When factoring a previously ordered matrix using sfFactor(), 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 best when the matrix is quite sparse. The * user selects the type of partition used with Mode. If Mode is set * to spDIRECT_PARTITION, then the all rows are placed in the direct * addressing partition. Similarly, if Mode is set to * spINDIRECT_PARTITION, then the all rows are placed in the indirect * addressing partition. By setting Mode to spAUTO_PARTITION, the * user allows Sparse to select the partition for each row * individually. sfFactor() generally runs faster if Sparse is * allowed to choose its own partitioning, however choosing a * partition is expensive. 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. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix. * Mode (int *) [INTEGER or INTEGER*2] * Mode must be one of three special codes: spDIRECT_PARTITION, * spINDIRECT_PARTITION, or spAUTO_PARTITION. */ void sfPartition( Matrix, Mode ) long *Matrix; int *Mode; { /* Begin `sfPartition'. */ spPartition((char *)*Matrix, *Mode); } /* * SOLVE MATRIX EQUATION * * Performs forward elimination and back substitution to find the * unknown vector from the RHS vector and factored matrix. This * routine assumes that the pivots are associated with the lower * triangular (L) matrix and that the diagonal of the upper triangular * (U) matrix consists of ones. This routine arranges the computation * in different way than is traditionally used in order to exploit the * sparsity of the right-hand side. See the reference in spRevision. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix. * RHS (RealVector) [REAL (1) or DOUBLE PRECISION (1)] * RHS is the input data array, the right hand side. This data is * undisturbed and may be reused for other solves. * Solution (RealVector) [REAL (1) or DOUBLE PRECISION (1)] * Solution is the output data array. This routine is constructed such that * RHS and Solution can be the same array. * iRHS (RealVector) [REAL (1) or DOUBLE PRECISION (1)] * iRHS is the imaginary portion of the input data array, the right * hand side. This data is undisturbed and may be reused for other solves. * This argument is only necessary if matrix is complex and if * spSEPARATED_COMPLEX_VECTOR is set true. * iSolution (RealVector) [REAL (1) or DOUBLE PRECISION (1)] * iSolution is the imaginary portion of the output data array. This * routine is constructed such that iRHS and iSolution can be * the same array. This argument is only necessary if matrix is complex * and if spSEPARATED_COMPLEX_VECTOR is set true. * * >>> Obscure Macros * IMAG_VECTORS * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears * without a trace. */ /*VARARGS3*/ void sfSolve( Matrix, RHS, Solution IMAG_VECTORS ) long *Matrix; RealVector RHS, Solution IMAG_VECTORS; { /* Begin `sfSolve'. */ spSolve( (char *)*Matrix, RHS, Solution IMAG_VECTORS ); } #if TRANSPOSE /* * SOLVE TRANSPOSED MATRIX EQUATION * * Performs forward elimination and back substitution to find the * unknown vector from the RHS vector and transposed factored * matrix. This routine is useful when performing sensitivity analysis * on a circuit using the adjoint method. This routine assumes that * the pivots are associated with the untransposed lower triangular * (L) matrix and that the diagonal of the untransposed upper * triangular (U) matrix consists of ones. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix. * RHS (RealVector) [REAL (1) or DOUBLE PRECISION (1)] * RHS is the input data array, the right hand side. This data is * undisturbed and may be reused for other solves. * Solution (RealVector) [REAL (1) or DOUBLE PRECISION (1)] * Solution is the output data array. This routine is constructed such that * RHS and Solution can be the same array. * iRHS (RealVector) [REAL (1) or DOUBLE PRECISION (1)] * iRHS is the imaginary portion of the input data array, the right * hand side. This data is undisturbed and may be reused for other solves. * If spSEPARATED_COMPLEX_VECTOR is set false, or if matrix is real, there * is no need to supply this array. * iSolution (RealVector) [REAL (1) or DOUBLE PRECISION (1)] * iSolution is the imaginary portion of the output data array. This * routine is constructed such that iRHS and iSolution can be * the same array. If spSEPARATED_COMPLEX_VECTOR is set false, or if * matrix is real, there is no need to supply this array. * * >>> Obscure Macros * IMAG_VECTORS * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears * without a trace. */ /*VARARGS3*/ void sfSolveTransposed( Matrix, RHS, Solution IMAG_VECTORS ) long *Matrix; RealVector RHS, Solution IMAG_VECTORS; { /* Begin `sfSolveTransposed'. */ spSolveTransposed( (char *)*Matrix, RHS, Solution IMAG_VECTORS ); } #endif /* TRANSPOSE */ #if DOCUMENTATION /* * PRINT MATRIX * * Formats and send the matrix to standard output. Some elementary * statistics are also output. The matrix is output in a format that is * readable by people. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix. * PrintReordered (long *) [LOGICAL] * Indicates whether the matrix should be printed out in its original * form, as input by the user, or whether it should be printed in its * reordered form, as used by the matrix routines. A zero indicates that * the matrix should be printed as inputed, a one indicates that it * should be printed reordered. * Data (long *) [LOGICAL] * Boolean flag that when false indicates that output should be * compressed such that only the existence of an element should be * indicated rather than giving the actual value. Thus 10 times as many * can be printed on a row. A zero signifies that the matrix should * be printed compressed. A one indicates that the matrix should be * printed in all its glory. * Header (long *) [LOGICAL] * Flag indicating that extra information such as the row and column * numbers should be printed. */ void sfPrint( Matrix, Data, PrintReordered, Header ) long *Matrix, *PrintReordered, *Data, *Header; { /* Begin `sfPrint'. */ spPrint( (char *)*Matrix, (int)*PrintReordered, (int)*Data, (int)*Header ); } #endif /* DOCUMENTATION */ #if DOCUMENTATION /* * OUTPUT MATRIX TO FILE * * Writes matrix to file in format suitable to be read back in by the * matrix test program. Data is sent to a file with a fixed name because * it is impossible to pass strings from FORTRAN to C in a manner that is * portable. * * >>> Returns: * One is returned if routine was successful, otherwise zero is returned. * The calling function can query errno (the system global error variable) * as to the reason why this routine failed. * * >>> Arguments: [LOGICAL] * Matrix (long *) [INTEGER] * Pointer to matrix. * Reordered (long *) [LOGICAL] * Specifies whether matrix should be output in reordered form, * or in original order. * Data (long *) [LOGICAL] * Indicates that the element values should be output along with * the indices for each element. This parameter must be true if * matrix is to be read by the sparse test program. * Header (long *) [LOGICAL] * Indicates that header is desired. This parameter must be true if * matrix is to be read by the sparse test program. */ #define MATRIX_FILE_NAME "spMatrix" #define STATS_FILE_NAME "spStats" long sfFileMatrix( Matrix, Reordered, Data, Header ) long *Matrix, *Reordered, *Data, *Header; { /* Begin `sfFileMatrix'. */ return spFileMatrix( (char *)*Matrix, MATRIX_FILE_NAME, "", (int)*Reordered, (int)*Data, (int)*Header ); } #endif /* DOCUMENTATION */ #if DOCUMENTATION /* * OUTPUT SOURCE VECTOR TO FILE * * Writes vector to file in format suitable to be read back in by the * matrix test program. This routine should be executed after the function * sfFileMatrix. * * >>> Returns: * One is returned if routine was successful, otherwise zero is returned. * The calling function can query errno (the system global error variable) * as to the reason why this routine failed. * * >>> Arguments: * Matrix (long *) * Pointer to matrix. * RHS (RealNumber []) [REAL (1) or DOUBLE PRECISION (1)] * Right-hand side vector. This is only the real portion if * spSEPARATED_COMPLEX_VECTORS is true. * iRHS (RealNumber []) [REAL (1) or DOUBLE PRECISION (1)] * Right-hand side vector, imaginary portion. Not necessary if matrix * is real or if spSEPARATED_COMPLEX_VECTORS is set false. */ int sfFileVector( Matrix, RHS IMAG_RHS ) long *Matrix; RealVector RHS IMAG_RHS; { /* Begin `sfFileVector'. */ return spFileVector( (char *)*Matrix, MATRIX_FILE_NAME, RHS IMAG_RHS ); } #endif /* DOCUMENTATION */ #if DOCUMENTATION /* * OUTPUT STATISTICS TO FILE * * Writes useful information concerning the matrix to a file. Should be * executed after the matrix is factored. * * >>> Returns: [LOGICAL] * One is returned if routine was successful, otherwise zero is returned. * The calling function can query errno (the system global error variable) * as to the reason why this routine failed. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix. */ int sfFileStats( Matrix ) long *Matrix; { /* Begin `sfFileStats'. */ return spFileStats( (char *)*Matrix, STATS_FILE_NAME, "" ); } #endif /* DOCUMENTATION */ #if MODIFIED_NODAL /* * PREORDER MODIFIED NODE ADMITTANCE MATRIX TO REMOVE ZEROS FROM DIAGONAL * * This routine massages modified node admittance matrices to remove * zeros from the diagonal. It takes advantage of the fact that the * row and column associated with a zero diagonal usually have * structural ones placed symmetricly. This routine should be used * only on modified node admittance matrices and should be executed * after the matrix has been built but before the factorization * begins. It should be executed for the initial factorization only * and should be executed before the rows have been linked. Thus it * should be run before using spScale(), spMultiply(), * spDeleteRowAndCol(), or spNorm(). * * This routine exploits the fact that the structural one are placed * in the matrix in symmetric twins. For example, the stamps for * grounded and a floating voltage sources are * grounded: floating: * [ x x 1 ] [ x x 1 ] * [ x x ] [ x x -1 ] * [ 1 ] [ 1 -1 ] * Notice for the grounded source, there is one set of twins, and for * the grounded, there are two sets. We remove the zero from the diagonal * by swapping the rows associated with a set of twins. For example: * grounded: floating 1: floating 2: * [ 1 ] [ 1 -1 ] [ x x 1 ] * [ x x ] [ x x -1 ] [ 1 -1 ] * [ x x 1 ] [ x x 1 ] [ x x -1 ] * * It is important to deal with any zero diagonals that only have one * set of twins before dealing with those that have more than one because * swapping row destroys the symmetry of any twins in the rows being * swapped, which may limit future moves. Consider * [ x x 1 ] * [ x x -1 1 ] * [ 1 -1 ] * [ 1 ] * There is one set of twins for diagonal 4 and two for diagonal3. * Dealing with diagonal for first requires swapping rows 2 and 4. * [ x x 1 ] * [ 1 ] * [ 1 -1 ] * [ x x -1 1 ] * We can now deal with diagonal 3 by swapping rows 1 and 3. * [ 1 -1 ] * [ 1 ] * [ x x 1 ] * [ x x -1 1 ] * And we are done, there are no zeros left on the diagonal. However, if * we originally dealt with diagonal 3 first, we could swap rows 2 and 3 * [ x x 1 ] * [ 1 -1 ] * [ x x -1 1 ] * [ 1 ] * Diagonal 4 no longer has a symmetric twin and we cannot continue. * * So we always take care of lone twins first. When none remain, we * choose arbitrarily a set of twins for a diagonal with more than one set * and swap the rows corresponding to that twin. We then deal with any * lone twins that were created and repeat the procedure until no * zero diagonals with symmetric twins remain. * * In this particular implementation, columns are swapped rather than rows. * The algorithm used in this function was developed by Ken Kundert and * Tom Quarles. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix to be preordered. */ void sfMNA_Preorder( Matrix ) long *Matrix; { /* Begin `sfMNA_Preorder'. */ spMNA_Preorder( (char *)*Matrix ); } #endif /* MODIFIED_NODAL */ #if SCALING /* * SCALE MATRIX * * This function scales the matrix to enhance the possibility of * finding a good pivoting order. Note that scaling enhances accuracy * of the solution only if it affects the pivoting order, so it makes * no sense to scale the matrix before spFactor(). If scaling is * desired it should be done before spOrderAndFactor(). There * are several things to take into account when choosing the scale * factors. First, the scale factors are directly multiplied against * the elements in the matrix. To prevent roundoff, each scale factor * should be equal to an integer power of the number base of the * machine. Since most machines operate in base two, scale factors * should be a power of two. Second, the matrix should be scaled such * that the matrix of element uncertainties is equilibrated. Third, * this function multiplies the scale factors by the elements, so if * one row tends to have uncertainties 1000 times smaller than the * other rows, then its scale factor should be 1024, not 1/1024. * Fourth, to save time, this function does not scale rows or columns * if their scale factors are equal to one. Thus, the scale factors * should be normalized to the most common scale factor. Rows and * columns should be normalized separately. For example, if the size * of the matrix is 100 and 10 rows tend to have uncertainties near * 1e-6 and the remaining 90 have uncertainties near 1e-12, then the * scale factor for the 10 should be 1/1,048,576 and the scale factors * for the remaining 90 should be 1. Fifth, since this routine * directly operates on the matrix, it is necessary to apply the scale * factors to the RHS and Solution vectors. It may be easier to * simply use spOrderAndFactor() on a scaled matrix to choose the * pivoting order, and then throw away the matrix. Subsequent * factorizations, performed with spFactor(), will not need to have * the RHS and Solution vectors descaled. Lastly, this function * should not be executed before the function spMNA_Preorder. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix to be scaled. * SolutionScaleFactors (RealVector) [REAL(1) or DOUBLE PRECISION(1)] * The array of Solution scale factors. These factors scale the columns. * All scale factors are real valued. * RHS_ScaleFactors (RealVector) [REAL(1) or DOUBLE PRECISION(1)] * The array of RHS scale factors. These factors scale the rows. * All scale factors are real valued. */ void sfScale( Matrix, RHS_ScaleFactors, SolutionScaleFactors ) long *Matrix; RealVector RHS_ScaleFactors, SolutionScaleFactors; { /* Begin `sfScale'. */ spScale( (char *)*Matrix, RHS_ScaleFactors, SolutionScaleFactors ); } #endif /* SCALING */ #if MULTIPLICATION /* * MATRIX MULTIPLICATION * * Multiplies matrix by solution vector to find source vector. * Assumes matrix has not been factored. This routine can be used * as a test to see if solutions are correct. It should not be used * before PreorderFoModifiedNodal(). * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix. * RHS (RealVector) [REAL(1) or DOUBLE PRECISION(1)] * RHS is the right hand side. This is what is being solved for. * Solution (RealVector) [REAL(1) or DOUBLE PRECISION(1)] * Solution is the vector being multiplied by the matrix. * iRHS (RealVector) [REAL(1) or DOUBLE PRECISION(1)] * iRHS is the imaginary portion of the right hand side. This is * what is being solved for. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * iSolution (RealVector) [REAL(1) or DOUBLE PRECISION(1)] * iSolution is the imaginary portion of the vector being multiplied * by the matrix. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * * >>> Obscure Macros * IMAG_VECTORS * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears * without a trace. */ void sfMultiply( Matrix, RHS, Solution IMAG_VECTORS ) long *Matrix; RealVector Solution, RHS IMAG_VECTORS; { /* Begin `sfMultiply'. */ spMultiply( (char *)*Matrix, RHS, Solution IMAG_VECTORS ); } #endif /* MULTIPLICATION */ #if MULTIPLICATION AND TRANSPOSE /* * TRANSPOSED MATRIX MULTIPLICATION * * Multiplies transposed matrix by solution vector to find source vector. * Assumes matrix has not been factored. This routine can be used * as a test to see if solutions are correct. It should not be used * before PreorderFoModifiedNodal(). * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix. * RHS (RealVector) [REAL(1) or DOUBLE PRECISION(1)] * RHS is the right hand side. This is what is being solved for. * Solution (RealVector) [REAL(1) or DOUBLE PRECISION(1)] * Solution is the vector being multiplied by the matrix. * iRHS (RealVector) [REAL(1) or DOUBLE PRECISION(1)] * iRHS is the imaginary portion of the right hand side. This is * what is being solved for. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * iSolution (RealVector) [REAL(1) or DOUBLE PRECISION(1)] * iSolution is the imaginary portion of the vector being multiplied * by the matrix. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * * >>> Obscure Macros * IMAG_VECTORS * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears * without a trace. */ void sfMultTransposed( Matrix, RHS, Solution IMAG_VECTORS ) long *Matrix; RealVector Solution, RHS IMAG_VECTORS; { /* Begin `sfMultTransposed'. */ spMultTransposed( (char *)*Matrix, RHS, Solution IMAG_VECTORS ); } #endif /* MULTIPLICATION AND TRANSPOSE */ #if DETERMINANT /* * CALCULATE DETERMINANT * * This routine in capable of calculating the determinant of the * matrix once the LU factorization has been performed. Hence, only * use this routine after spFactor() and before spClear(). * The determinant equals the product of all the diagonal elements of * the lower triangular matrix L, except that this product may need * negating. Whether the product or the negative product equals the * determinant is determined by the number of row and column * interchanges performed. Note that the determinants of matrices can * be very large or very small. On large matrices, the determinant * can be far larger or smaller than can be represented by a floating * point number. For this reason the determinant is scaled to a * reasonable value and the logarithm of the scale factor is returned. * * >>> Arguments: * Matrix (long *) [INTEGER] * A pointer to the matrix for which the determinant is desired. * pExponent (int *) [INTEGER or INTEGER*2] * The logarithm base 10 of the scale factor for the determinant. To * find * the actual determinant, Exponent should be added to the exponent of * DeterminantReal. * pDeterminant (RealNumber *) [REAL or DOUBLE PRECISION] * The real portion of the determinant. This number is scaled to be * greater than or equal to 1.0 and less than 10.0. * piDeterminant (RealNumber *) [REAL or DOUBLE PRECISION] * The imaginary portion of the determinant. When the matrix is real * this pointer need not be supplied, nothing will be returned. This * number is scaled to be greater than or equal to 1.0 and less than 10.0. */ #if spCOMPLEX void sfDeterminant( Matrix, pExponent, pDeterminant, piDeterminant ) long *Matrix; RealNumber *pDeterminant, *piDeterminant; int *pExponent; { /* Begin `sfDeterminant'. */ spDeterminant( (char *)*Matrix, pExponent, pDeterminant, piDeterminant ); } #else /* spCOMPLEX */ void sfDeterminant( Matrix, pExponent, pDeterminant ) long *Matrix; RealNumber *pDeterminant; int *pExponent; { /* Begin `sfDeterminant'. */ spDeterminant( (char *)*Matrix, pExponent, pDeterminant ); } #endif /* spCOMPLEX */ #endif /* DETERMINANT */ /* * RETURN MATRIX ERROR STATUS * * This function is used to determine the error status of the given matrix. * * >>> Returned: [INTEGER or INTEGER*2] * The error status of the given matrix. * * >>> Arguments: * Matrix (long *) [INTEGER] * The matrix for which the error status is desired. */ int sfError( Matrix ) long *Matrix; { /* Begin `sfError'. */ return spError( (char *)*Matrix ); } /* * WHERE IS MATRIX SINGULAR * * This function returns the row and column number where the matrix was * detected as singular or where a zero was detected on the diagonal. * * >>> Arguments: * Matrix (long *) [INTEGER] * The matrix for which the error status is desired. * pRow (int *) [INTEGER or INTEGER*2] * The row number. * pCol (int *) [INTEGER or INTEGER*2] * The column number. */ void sfWhereSingular( Matrix, Row, Col ) long *Matrix; int *Row, *Col; { /* Begin `sfWhereSingular'. */ spWhereSingular( (char *)*Matrix, Row, Col ); } /* * MATRIX SIZE * * Returns the size of the matrix. Either the internal or external size of * the matrix is returned. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix. * External (BOOLEAN) [LOGICAL] * If External is set true, the external size , i.e., the value of the * largest external row or column number encountered is returned. * Otherwise the true size of the matrix is returned. These two sizes * may differ if the TRANSLATE option is set true. */ int sfGetSize( Matrix, External ) long *Matrix, *External; { /* Begin `sfGetSize'. */ return spGetSize( (char *)*Matrix, (BOOLEAN)*External ); } /* * SET MATRIX COMPLEX OR REAL * * Forces matrix to be either real or complex. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix. */ void sfSetReal( Matrix ) long *Matrix; { /* Begin `sfSetReal'. */ spSetReal( (char *)*Matrix ); } void sfSetComplex( Matrix ) long *Matrix; { /* Begin `sfSetComplex'. */ spSetComplex( (char *)*Matrix ); } /* * ELEMENT OR FILL-IN COUNT * * Two functions used to return simple statistics. Either the number * of total elements, or the number of fill-ins can be returned. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to matrix. */ int sfFillinCount( Matrix ) long *Matrix; { /* Begin `sfFillinCount'. */ return spFillinCount( (char *)*Matrix ); } int sfElementCount( Matrix ) long *Matrix; { /* Begin `sfElementCount'. */ return spElementCount( (char *)*Matrix ); } #if TRANSLATE AND DELETE /* * DELETE A ROW AND COLUMN FROM THE MATRIX * * Deletes a row and a column from a matrix. * * Sparse will abort if an attempt is made to delete a row or column that * doesn't exist. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix in which the row and column are to be deleted. * Row (int) [INTEGER or INTEGER*2] * Row to be deleted. * Col (int) [INTEGER or INTEGER*2] * Column to be deleted. */ void sfDeleteRowAndCol( Matrix, Row, Col ) long *Matrix; int *Row, *Col; { /* Begin `sfDeleteRowAndCol'. */ spDeleteRowAndCol( (char *)*Matrix, *Row, *Col ); } #endif #if PSEUDOCONDITION /* * CALCULATE PSEUDOCONDITION * * Computes the magnitude of the ratio of the largest to the smallest * pivots. This quantity is an indicator of ill-conditioning in the * matrix. If this ratio is large, and if the matrix is scaled such * that uncertainties in the RHS and the matrix entries are * equilibrated, then the matrix is ill-conditioned. However, a small * ratio does not necessarily imply that the matrix is * well-conditioned. This routine must only be used after a matrix * has been factored by sfOrderAndFactor() or sfFactor() and before it * is cleared by sfClear() or spInitialize(). The pseudocondition is faster * to compute than the condition number calculated by sfCondition(), but * is not as informative. * * >>> Returns: [REAL or DOUBLE PRECISION] * The magnitude of the ratio of the largest to smallest pivot used during * previous factorization. If the matrix was singular, zero is returned. * * >>> Arguments: * Matrix (long *) * Pointer to the matrix. */ RealNumber sfPseudoCondition( Matrix ) long *Matrix; { /* Begin `sfPseudoCondition'. */ return spPseudoCondition( (char *)Matrix ); } #endif #if CONDITION /* * ESTIMATE CONDITION NUMBER * * Computes an estimate of the condition number using a variation on * the LINPACK condition number estimation algorithm. This quantity is * an indicator of ill-conditioning in the matrix. To avoid problems * with overflow, the reciprocal of the condition number is returned. * If this number is small, and if the matrix is scaled such that * uncertainties in the RHS and the matrix entries are equilibrated, * then the matrix is ill-conditioned. If the this number is near * one, the matrix is well conditioned. This routine must only be * used after a matrix has been factored by sfOrderAndFactor() or * sfFactor() and before it is cleared by sfClear() or spInitialize(). * * Unlike the LINPACK condition number estimator, this routines * returns the L infinity condition number. This is an artifact of * Sparse placing ones on the diagonal of the upper triangular matrix * rather than the lower. This difference should be of no importance. * * References: * A.K. Cline, C.B. Moler, G.W. Stewart, J.H. Wilkinson. An estimate * for the condition number of a matrix. SIAM Journal on Numerical * Analysis. Vol. 16, No. 2, pages 368-375, April 1979. * * J.J. Dongarra, C.B. Moler, J.R. Bunch, G.W. Stewart. LINPACK * User's Guide. SIAM, 1979. * * Roger G. Grimes, John G. Lewis. Condition number estimation for * sparse matrices. SIAM Journal on Scientific and Statistical * Computing. Vol. 2, No. 4, pages 384-388, December 1981. * * Dianne Prost O'Leary. Estimating matrix condition numbers. SIAM * Journal on Scientific and Statistical Computing. Vol. 1, No. 2, * pages 205-209, June 1980. * * >>> Returns: [REAL or DOUBLE PRECISION] * The reciprocal of the condition number. If the matrix was singular, * zero is returned. * * >>> Arguments: * eMatrix (long *) * Pointer to the matrix. * NormOfMatrix (RealNumber *) [REAL or DOUBLE PRECISION] * The L-infinity norm of the unfactored matrix as computed by * spNorm(). * pError (int *) [INTEGER or INTEGER*2] * Used to return error code. * * >>> Possible errors: * spSINGULAR * spNO_MEMORY */ RealNumber sfCondition( Matrix, NormOfMatrix, pError ) long *Matrix; RealNumber *NormOfMatrix; int *pError; { /* Begin `sfCondition'. */ return spCondition( (char *)*Matrix, *NormOfMatrix, pError ); } /* * L-INFINITY MATRIX NORM * * Computes the L-infinity norm of an unfactored matrix. It is a fatal * error to pass this routine a factored matrix. * * One difficulty is that the rows may not be linked. * * >>> Returns: [REAL or DOUBLE PRECISION] * The largest absolute row sum of matrix. * * >>> Arguments: * Matrix (long *) * Pointer to the matrix. */ RealNumber sfNorm( Matrix ) long *Matrix; { /* Begin `sfNorm'. */ return spNorm( (char *)*Matrix ); } #endif /* CONDITION */ #if STABILITY /* * STABILITY OF FACTORIZATION * * The following routines are used to gauge the stability of a * factorization. If the factorization is determined to be too unstable, * then the matrix should be reordered. The routines compute quantities * that are needed in the computation of a bound on the error attributed * to any one element in the matrix during the factorization. In other * words, there is a matrix E = [e_ij] of error terms such that A+E = LU. * This routine finds a bound on |e_ij|. Erisman & Reid [1] showed that * |e_ij| < 3.01 u rho m_ij, where u is the machine rounding unit, * rho = max a_ij where the max is taken over every row i, column j, and * step k, and m_ij is the number of multiplications required in the * computation of l_ij if i > j or u_ij otherwise. Barlow [2] showed that * rho < max_i || l_i ||_p max_j || u_j ||_q where 1/p + 1/q = 1. * * The first routine finds the magnitude on the largest element in the * matrix. If the matrix has not yet been factored, the largest * element is found by direct search. If the matrix is factored, a * bound on the largest element in any of the reduced submatrices is * computed using Barlow with p = oo and q = 1. The ratio of these * two numbers is the growth, which can be used to determine if the * pivoting order is adequate. A large growth implies that * considerable error has been made in the factorization and that it * is probably a good idea to reorder the matrix. If a large growth * in encountered after using spFactor(), reconstruct the matrix and * refactor using spOrderAndFactor(). If a large growth is * encountered after using spOrderAndFactor(), refactor using * spOrderAndFactor() with the pivot threshold increased, say to 0.1. * * Using only the size of the matrix as an upper bound on m_ij and * Barlow's bound, the user can estimate the size of the matrix error * terms e_ij using the bound of Erisman and Reid. The second routine * computes a tighter bound (with more work) based on work by Gear * [3], |e_ij| < 1.01 u rho (t c^3 + (1 + t)c^2) where t is the * threshold and c is the maximum number of off-diagonal elements in * any row of L. The expensive part of computing this bound is * determining the maximum number of off-diagonals in L, which changes * only when the order of the matrix changes. This number is computed * and saved, and only recomputed if the matrix is reordered. * * [1] A. M. Erisman, J. K. Reid. Monitoring the stability of the * triangular factorization of a sparse matrix. Numerische * Mathematik. Vol. 22, No. 3, 1974, pp 183-186. * * [2] J. L. Barlow. A note on monitoring the stability of triangular * decomposition of sparse matrices. "SIAM Journal of Scientific * and Statistical Computing." Vol. 7, No. 1, January 1986, pp 166-168. * * [3] I. S. Duff, A. M. Erisman, J. K. Reid. "Direct Methods for Sparse * Matrices." Oxford 1986. pp 99. */ /* * LARGEST ELEMENT IN MATRIX * * >>> Returns: [REAL or DOUBLE PRECISION] * If matrix is not factored, returns the magnitude of the largest element in * the matrix. If the matrix is factored, a bound on the magnitude of the * largest element in any of the reduced submatrices is returned. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix. */ RealNumber sfLargestElement( Matrix ) long *Matrix; { /* Begin `sfLargestElement'. */ return spLargestElement( (char *)Matrix ); } /* * MATRIX ROUNDOFF ERROR * * >>> Returns: [REAL or DOUBLE PRECISION] * Returns a bound on the magnitude of the largest element in E = A - LU. * * >>> Arguments: * Matrix (long *) [INTEGER] * Pointer to the matrix. * Rho (RealNumber *) [REAL or DOUBLE PRECISION] * The bound on the magnitude of the largest element in any of the * reduced submatrices. This is the number computed by the function * spLargestElement() when given a factored matrix. If this number is * negative, the bound will be computed automatically. */ RealNumber sfRoundoff( Matrix, Rho ) long *Matrix; RealNumber *Rho; { /* Begin `sfRoundoff'. */ return spRoundoff( (char *)*Matrix, *Rho ); } #endif #endif /* FORTRAN */ @EOF chmod 644 spFortran.c echo x - spOutput.c cat >spOutput.c <<'@EOF' /* * MATRIX OUTPUT MODULE * * Author: Advisor: * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli * UC Berkeley * * This file contains the output-to-file and output-to-screen routines for * the matrix package. * * >>> User accessible functions contained in this file: * spPrint * spFileMatrix * spFileVector * spFileStats * * >>> Other functions contained in this file: */ /* * 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 lint static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; static char RCSid[] = "$Header: spOutput.c,v 1.2 88/06/18 11:15:36 kundert Exp $"; #endif /* * IMPORTS * * >>> Import descriptions: * spConfig.h * Macros that customize the sparse matrix routines. * spMatrix.h * Macros and declarations to be imported by the user. * spDefs.h * Matrix type and macro definitions for the sparse matrix routines. */ #define spINSIDE_SPARSE #include "spConfig.h" #include "spMatrix.h" #include "spDefs.h" #if DOCUMENTATION /* * PRINT MATRIX * * Formats and send the matrix to standard output. Some elementary * statistics are also output. The matrix is output in a format that is * readable by people. * * >>> Arguments: * Matrix (char *) * Pointer to matrix. * PrintReordered (int) * Indicates whether the matrix should be printed out in its original * form, as input by the user, or whether it should be printed in its * reordered form, as used by the matrix routines. A zero indicates that * the matrix should be printed as inputed, a one indicates that it * should be printed reordered. * Data (int) * Boolean flag that when false indicates that output should be * compressed such that only the existence of an element should be * indicated rather than giving the actual value. Thus 11 times as * many can be printed on a row. A zero signifies that the matrix * should be printed compressed. A one indicates that the matrix * should be printed in all its glory. * Header (int) * Flag indicating that extra information should be given, such as row * and column numbers. * * >>> Local variables: * Col (int) * Column being printed. * ElementCount (int) * Variable used to count the number of nonzero elements in the matrix. * LargestElement (RealNumber) * The magnitude of the largest element in the matrix. * LargestDiag (RealNumber) * The magnitude of the largest diagonal in the matrix. * Magnitude (RealNumber) * The absolute value of the matrix element being printed. * PrintOrdToIntColMap (int []) * A translation array that maps the order that columns will be * printed in (if not PrintReordered) to the internal column numbers. * PrintOrdToIntRowMap (int []) * A translation array that maps the order that rows will be * printed in (if not PrintReordered) to the internal row numbers. * pElement (ElementPtr) * Pointer to the element in the matrix that is to be printed. * pImagElements (ElementPtr [ ]) * Array of pointers to elements in the matrix. These pointers point * to the elements whose real values have just been printed. They are * used to quickly access those same elements so their imaginary values * can be printed. * Row (int) * Row being printed. * Size (int) * The size of the matrix. * SmallestDiag (RealNumber) * The magnitude of the smallest diagonal in the matrix. * SmallestElement (RealNumber) * The magnitude of the smallest element in the matrix excluding zero * elements. * StartCol (int) * The column number of the first column to be printed in the group of * columns currently being printed. * StopCol (int) * The column number of the last column to be printed in the group of * columns currently being printed. * Top (int) * The largest expected external row or column number. */ void spPrint( eMatrix, PrintReordered, Data, Header ) char *eMatrix; int PrintReordered, Data, Header; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register int J = 0; int I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns, ElementCount = 0; double Magnitude, SmallestDiag, SmallestElement; double LargestElement = 0.0, LargestDiag = 0.0; ElementPtr pElement, pImagElements[PRINTER_WIDTH/10+1]; int *PrintOrdToIntRowMap, *PrintOrdToIntColMap; /* Begin `spPrint'. */ ASSERT( IS_SPARSE( Matrix ) ); Size = Matrix->Size; /* Create a packed external to internal row and column translation array. */ # if TRANSLATE Top = Matrix->AllocatedExtSize; #else Top = Matrix->AllocatedSize; #endif CALLOC( PrintOrdToIntRowMap, int, Top + 1 ); CALLOC( PrintOrdToIntColMap, int, Top + 1 ); if ( PrintOrdToIntRowMap == NULL OR PrintOrdToIntColMap == NULL) { Matrix->Error = spNO_MEMORY; return; } for (I = 1; I <= Size; I++) { PrintOrdToIntRowMap[ Matrix->IntToExtRowMap[I] ] = I; PrintOrdToIntColMap[ Matrix->IntToExtColMap[I] ] = I; } /* Pack the arrays. */ for (J = 1, I = 1; I <= Top; I++) { if (PrintOrdToIntRowMap[I] != 0) PrintOrdToIntRowMap[ J++ ] = PrintOrdToIntRowMap[ I ]; } for (J = 1, I = 1; I <= Top; I++) { if (PrintOrdToIntColMap[I] != 0) PrintOrdToIntColMap[ J++ ] = PrintOrdToIntColMap[ I ]; } /* Print header. */ if (Header) { printf("MATRIX SUMMARY\n\n"); printf("Size of matrix = %1u x %1u.\n", Size, Size); if ( Matrix->Reordered AND PrintReordered ) printf("Matrix has been reordered.\n"); putchar('\n'); if ( Matrix->Factored ) printf("Matrix after factorization:\n"); else printf("Matrix before factorization:\n"); SmallestElement = LARGEST_REAL; SmallestDiag = SmallestElement; } /* Determine how many columns to use. */ Columns = PRINTER_WIDTH; if (Header) Columns -= 5; if (Data) Columns = (Columns+1) / 10; /* * Print matrix by printing groups of complete columns until all the columns * are printed. */ J = 0; while ( J <= Size ) /* Calculate index of last column to printed in this group. */ { StopCol = StartCol + Columns - 1; if (StopCol > Size) StopCol = Size; /* Label the columns. */ if (Header) { if (Data) { printf(" "); for (I = StartCol; I <= StopCol; I++) { if (PrintReordered) Col = I; else Col = PrintOrdToIntColMap[I]; printf(" %9d", Matrix->IntToExtColMap[ Col ]); } printf("\n\n"); } else { if (PrintReordered) printf("Columns %1d to %1d.\n",StartCol,StopCol); else { printf("Columns %1d to %1d.\n", Matrix->IntToExtColMap[ PrintOrdToIntColMap[StartCol] ], Matrix->IntToExtColMap[ PrintOrdToIntColMap[StopCol] ]); } } } /* Print every row ... */ for (I = 1; I <= Size; I++) { if (PrintReordered) Row = I; else Row = PrintOrdToIntRowMap[I]; if (Header) { if (PrintReordered AND NOT Data) printf("%4d", I); else printf("%4d", Matrix->IntToExtRowMap[ Row ]); if (NOT Data) putchar(' '); } /* ... in each column of the group. */ for (J = StartCol; J <= StopCol; J++) { if (PrintReordered) Col = J; else Col = PrintOrdToIntColMap[J]; pElement = Matrix->FirstInCol[Col]; while(pElement != NULL AND pElement->Row != Row) pElement = pElement->NextInCol; if (Data) pImagElements[J - StartCol] = pElement; if (pElement != NULL) /* Case where element exists */ { if (Data) printf(" %9.3lg", (double)pElement->Real); else putchar('x'); /* Update status variables */ if ( (Magnitude = ELEMENT_MAG(pElement)) > LargestElement ) LargestElement = Magnitude; if ((Magnitude < SmallestElement) AND (Magnitude != 0.0)) SmallestElement = Magnitude; ElementCount++; } /* Case where element is structurally zero */ else { if (Data) printf(" ..."); else putchar('.'); } } putchar('\n'); #if spCOMPLEX if (Matrix->Complex AND Data) { printf(" "); for (J = StartCol; J <= StopCol; J++) { if (pImagElements[J - StartCol] != NULL) { printf(" %8.2lgj", (double)pImagElements[J-StartCol]->Imag); } else printf(" "); } putchar('\n'); } #endif /* spCOMPLEX */ } /* Calculate index of first column in next group. */ StartCol = StopCol; StartCol++; putchar('\n'); } if (Header) { printf("\nLargest element in matrix = %-1.4lg.\n", LargestElement); printf("Smallest element in matrix = %-1.4lg.\n", SmallestElement); /* Search for largest and smallest diagonal values */ for (I = 1; I <= Size; I++) { if (Matrix->Diag[I] != NULL) { Magnitude = ELEMENT_MAG( Matrix->Diag[I] ); if ( Magnitude > LargestDiag ) LargestDiag = Magnitude; if ( Magnitude < SmallestDiag ) SmallestDiag = Magnitude; } } /* Print the largest and smallest diagonal values */ if ( Matrix->Factored ) { printf("\nLargest diagonal element = %-1.4lg.\n", LargestDiag); printf("Smallest diagonal element = %-1.4lg.\n", SmallestDiag); } else { printf("\nLargest pivot element = %-1.4lg.\n", LargestDiag); printf("Smallest pivot element = %-1.4lg.\n", SmallestDiag); } /* Calculate and print sparsity and number of fill-ins created. */ printf("\nDensity = %2.2lf%%.\n", ((double)(ElementCount * 100)) / ((double)(Size * Size))); if (NOT Matrix->NeedsOrdering) printf("Number of fill-ins = %1d.\n", Matrix->Fillins); } putchar('\n'); (void)fflush(stdout); FREE(PrintOrdToIntColMap); FREE(PrintOrdToIntRowMap); return; } /* * OUTPUT MATRIX TO FILE * * Writes matrix to file in format suitable to be read back in by the * matrix test program. * * >>> Returns: * One is returned if routine was successful, otherwise zero is returned. * The calling function can query errno (the system global error variable) * as to the reason why this routine failed. * * >>> Arguments: * Matrix (char *) * Pointer to matrix. * File (char *) * Name of file into which matrix is to be written. * Label (char *) * String that is transferred to file and is used as a label. * Reordered (BOOLEAN) * Specifies whether matrix should be output in reordered form, * or in original order. * Data (BOOLEAN) * Indicates that the element values should be output along with * the indices for each element. This parameter must be true if * matrix is to be read by the sparse test program. * Header (BOOLEAN) * Indicates that header is desired. This parameter must be true if * matrix is to be read by the sparse test program. * * >>> Local variables: * Col (int) * The original column number of the element being output. * pElement (ElementPtr) * Pointer to an element in the matrix. * pMatrixFile (FILE *) * File pointer to the matrix file. * Row (int) * The original row number of the element being output. * Size (int) * The size of the matrix. */ int spFileMatrix( eMatrix, File, Label, Reordered, Data, Header ) char *eMatrix, *Label, *File; int Reordered, Data, Header; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register int I, Size; register ElementPtr pElement; int Row, Col, Err; FILE *pMatrixFile, *fopen(); /* Begin `spFileMatrix'. */ ASSERT( IS_SPARSE( Matrix ) ); /* Open file matrix file in write mode. */ if ((pMatrixFile = fopen(File, "w")) == NULL) return 0; /* Output header. */ Size = Matrix->Size; if (Header) { if (Matrix->Factored AND Data) { Err = fprintf ( pMatrixFile, "Warning : The following matrix is factored in to LU form.\n" ); } if (Err < 0) return 0; if (fprintf(pMatrixFile, "%s\n", Label) < 0) return 0; Err = fprintf( pMatrixFile, "%d\t%s\n", Size, (Matrix->Complex ? "complex" : "real")); if (Err < 0) return 0; } /* Output matrix. */ if (NOT Data) { for (I = 1; I <= Size; I++) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { if (Reordered) { Row = pElement->Row; Col = I; } else { Row = Matrix->IntToExtRowMap[pElement->Row]; Col = Matrix->IntToExtColMap[I]; } pElement = pElement->NextInCol; if (fprintf(pMatrixFile, "%d\t%d\n", Row, Col) < 0) return 0; } } /* Output terminator, a line of zeros. */ if (Header) if (fprintf(pMatrixFile, "0\t0\n") < 0) return 0; } #if spCOMPLEX if (Data AND Matrix->Complex) { for (I = 1; I <= Size; I++) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { if (Reordered) { Row = pElement->Row; Col = I; } else { Row = Matrix->IntToExtRowMap[pElement->Row]; Col = Matrix->IntToExtColMap[I]; } Err = fprintf ( pMatrixFile,"%d\t%d\t%-.15lg\t%-.15lg\n", Row, Col, (double)pElement->Real, (double)pElement->Imag ); if (Err < 0) return 0; pElement = pElement->NextInCol; } } /* Output terminator, a line of zeros. */ if (Header) if (fprintf(pMatrixFile,"0\t0\t0.0\t0.0\n") < 0) return 0; } #endif /* spCOMPLEX */ #if REAL if (Data AND NOT Matrix->Complex) { for (I = 1; I <= Size; I++) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { Row = Matrix->IntToExtRowMap[pElement->Row]; Col = Matrix->IntToExtColMap[I]; Err = fprintf ( pMatrixFile,"%d\t%d\t%-.15lg\n", Row, Col, (double)pElement->Real ); if (Err < 0) return 0; pElement = pElement->NextInCol; } } /* Output terminator, a line of zeros. */ if (Header) if (fprintf(pMatrixFile,"0\t0\t0.0\n") < 0) return 0; } #endif /* REAL */ /* Close file. */ if (fclose(pMatrixFile) < 0) return 0; return 1; } /* * OUTPUT SOURCE VECTOR TO FILE * * Writes vector to file in format suitable to be read back in by the * matrix test program. This routine should be executed after the function * spFileMatrix. * * >>> Returns: * One is returned if routine was successful, otherwise zero is returned. * The calling function can query errno (the system global error variable) * as to the reason why this routine failed. * * >>> Arguments: * Matrix (char *) * Pointer to matrix. * File (char *) * Name of file into which matrix is to be written. * RHS (RealNumber []) * Right-hand side vector. This is only the real portion if * spSEPARATED_COMPLEX_VECTORS is true. * iRHS (RealNumber []) * Right-hand side vector, imaginary portion. Not necessary if matrix * is real or if spSEPARATED_COMPLEX_VECTORS is set false. * * >>> Local variables: * pMatrixFile (FILE *) * File pointer to the matrix file. * Size (int) * The size of the matrix. * * >>> Obscure Macros * IMAG_RHS * Replaces itself with `, iRHS' if the options spCOMPLEX and * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears * without a trace. */ int spFileVector( eMatrix, File, RHS IMAG_RHS ) char *eMatrix, *File; RealVector RHS IMAG_RHS; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register int I, Size, Err; FILE *pMatrixFile; FILE *fopen(); /* Begin `spFileVector'. */ ASSERT( IS_SPARSE( Matrix ) AND RHS != NULL) /* Open File in append mode. */ if ((pMatrixFile = fopen(File,"a")) == NULL) return 0; /* Correct array pointers for ARRAY_OFFSET. */ #if NOT ARRAY_OFFSET #if spCOMPLEX if (Matrix->Complex) { #if spSEPARATED_COMPLEX_VECTORS ASSERT(iRHS != NULL) --RHS; --iRHS; #else RHS -= 2; #endif } else #endif /* spCOMPLEX */ --RHS; #endif /* NOT ARRAY_OFFSET */ /* Output vector. */ Size = Matrix->Size; #if spCOMPLEX if (Matrix->Complex) { #if spSEPARATED_COMPLEX_VECTORS for (I = 1; I <= Size; I++) { Err = fprintf ( pMatrixFile, "%-.15lg\t%-.15lg\n", (double)RHS[I], (double)iRHS[I] ); if (Err < 0) return 0; } #else for (I = 1; I <= Size; I++) { Err = fprintf ( pMatrixFile, "%-.15lg\t%-.15lg\n", (double)RHS[2*I], (double)RHS[2*I+1] ); if (Err < 0) return 0; } #endif } #endif /* spCOMPLEX */ #if REAL AND spCOMPLEX else #endif #if REAL { for (I = 1; I <= Size; I++) { if (fprintf(pMatrixFile, "%-.15lg\n", (double)RHS[I]) < 0) return 0; } } #endif /* REAL */ /* Close file. */ if (fclose(pMatrixFile) < 0) return 0; return 1; } /* * OUTPUT STATISTICS TO FILE * * Writes useful information concerning the matrix to a file. Should be * executed after the matrix is factored. * * >>> Returns: * One is returned if routine was successful, otherwise zero is returned. * The calling function can query errno (the system global error variable) * as to the reason why this routine failed. * * >>> Arguments: * Matrix (char *) * Pointer to matrix. * File (char *) * Name of file into which matrix is to be written. * Label (char *) * String that is transferred to file and is used as a label. * * >>> Local variables: * Data (RealNumber) * The value of the matrix element being output. * LargestElement (RealNumber) * The largest element in the matrix. * NumberOfElements (int) * Number of nonzero elements in the matrix. * pElement (ElementPtr) * Pointer to an element in the matrix. * pStatsFile (FILE *) * File pointer to the statistics file. * Size (int) * The size of the matrix. * SmallestElement (RealNumber) * The smallest element in the matrix excluding zero elements. */ int spFileStats( eMatrix, File, Label ) char *eMatrix, *File, *Label; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register int Size, I; register ElementPtr pElement; int NumberOfElements; RealNumber Data, LargestElement, SmallestElement; FILE *pStatsFile, *fopen(); /* Begin `spFileStats'. */ ASSERT( IS_SPARSE( Matrix ) ); /* Open File in append mode. */ if ((pStatsFile = fopen(File, "a")) == NULL) return 0; /* Output statistics. */ Size = Matrix->Size; if (NOT Matrix->Factored) fprintf(pStatsFile, "Matrix has not been factored.\n"); fprintf(pStatsFile, "||| Starting new matrix |||\n"); fprintf(pStatsFile, "%s\n", Label); if (Matrix->Complex) fprintf(pStatsFile, "Matrix is complex.\n"); else fprintf(pStatsFile, "Matrix is real.\n"); fprintf(pStatsFile," Size = %d\n",Size); /* Search matrix. */ NumberOfElements = 0; LargestElement = 0.0; SmallestElement = LARGEST_REAL; for (I = 1; I <= Size; I++) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { NumberOfElements++; Data = ELEMENT_MAG(pElement); if (Data > LargestElement) LargestElement = Data; if (Data < SmallestElement AND Data != 0.0) SmallestElement = Data; pElement = pElement->NextInCol; } } SmallestElement = MIN( SmallestElement, LargestElement ); /* Output remaining statistics. */ fprintf(pStatsFile, " Initial number of elements = %d\n", NumberOfElements - Matrix->Fillins); fprintf(pStatsFile, " Initial average number of elements per row = %lf\n", (double)(NumberOfElements - Matrix->Fillins) / (double)Size); fprintf(pStatsFile, " Fill-ins = %d\n",Matrix->Fillins); fprintf(pStatsFile, " Average number of fill-ins per row = %lf%%\n", (double)Matrix->Fillins / (double)Size); fprintf(pStatsFile, " Total number of elements = %d\n", NumberOfElements); fprintf(pStatsFile, " Average number of elements per row = %lf\n", (double)NumberOfElements / (double)Size); fprintf(pStatsFile," Density = %lf%%\n", (double)(100.0*NumberOfElements)/(double)(Size*Size)); fprintf(pStatsFile," Relative Threshold = %e\n", Matrix->RelThreshold); fprintf(pStatsFile," Absolute Threshold = %e\n", Matrix->AbsThreshold); fprintf(pStatsFile," Largest Element = %e\n", LargestElement); fprintf(pStatsFile," Smallest Element = %e\n\n\n", SmallestElement); /* Close file. */ (void)fclose(pStatsFile); return 1; } #endif /* DOCUMENTATION */ @EOF chmod 644 spOutput.c echo x - spRevision cat >spRevision <<'@EOF' /* * REVISION HISTORY * * Author: * Kenneth S. Kundert * University of California at Berkeley * College of Engineering * Department of Electrical Engineering and Computer Science * * Advising professor: * Alberto Sangiovanni-Vincentelli * * * Sparse is available for a $150.00 charge. The package includes the * source code on tape, the user's guide and a paper that discusses techniques * for solving sparse systems of equations [1]. To obtain a copy of Sparse, * send a check or money order payable to the Regents of the University of * California to: * * EECS Industrial Liaison Program * Cindy Manly * University of California * Berkeley, CA 94720 * * Please allow four weeks for delivery. * * The University often does not have the resources to consult with * users on how to use or modify these programs. We would, however, like * to be notified of any problems or errors in the material provided and * appreciate copies on tape of any troublesome matrices. If the program * is enhanced or converted to run on other systems, we would like to receive * copies of the modified program so that it can be made available to the * public. * * References * [1] Kenneth S. Kundert. Sparse matrix techniques. In "Circuit Analysis, * Simulation and Design," vol. 3, pt. 1, Albert E. Ruehli (editor). * North-Holland, 1986. */ /* * 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 above copyright notice 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. */ /* * >>> Current revision information: * $Author: kundert $ * $Date: 85/12/10 21:03:02 $ * $Revision: 2.2 $ */ /* * >>> History: * Revision 1.1 January 1985 * Initial release. * * Revision 1.1a 20 March 1985 * Modified DecomposeMatrix() and OrderAndDecomposeMatrix() so that * the parameters Growth, PseudoCondition and LargestElement may be * given as NULL. * * Revision 1.1b 28 March 1985 * Corrected a bug that caused OrderAndDecomposeMatrix() to reorder * the matrix every time it was called. Also made many of the global * variables defined in MtrxDecom.c static. * * Revision 1.2 October 1985 * This new version of Sparse is meant to make it more compatible * with interactive circuit simulators. In it the TRANSLATE * option was added, along with the ability to access the matrix * with AddElementToMatrix() and AddAdmittanceToMatrix() after it * has been reordered. Also added were the DeleteRowAndColFromMatrix(), * CleanMatrix(), GetMatrixSize(), SetMatrixReal() and SetMatrixComplex() * routines. * * Revision 1.2a April 1986 * Fixed a bug that caused the matrix frame to get freed twice and one * in the test program that caused sparse to crash when a complex matrix * with no source vector was attempted. * * Revision 1.2b July 1986 * Modified the test routine so that it allocates vectors from the heap * rather than the stack. * * Revision 1.2c February 1987 * Fixed off-by-one error in PreorderForModifiedNodal(). * * Revision 1.2d * Modified the pivot selection algorithm so that singletons also meet * numerical threshold criterion. Deleted some global variables. Modified * test program to add command line options among other things. Made * thresholds sticky, so that once a valid threshold has been specified * for a particular matrix, it becomes the new default. * * Revision 1.3a July 1988 * Made numerous changes. Highlights are: * Routine names shortened and made unique to 8 characters. * Unordering factorization is faster. * Added self initialization feature. * Sparse now aborts on errors that should not occur. * Cleaned up test program. * Separated growth and pseudocondition calculations from factorization. * Added LINPACK condition number estimator. * Rewrote spMNA_Preorder, algorithm extended to fix inadequacies. * Eliminated all global variables. * Added DIAGONAL_PIVOTING option (only diagonal pivoting before). */ @EOF chmod 644 spRevision echo x - spTest.c cat >spTest.c <<'@EOF' /* * TEST MODULE for the sparse matrix routines * * Author: Advisor: * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli * UC Berkeley * * This file contains the test routine for the sparse matrix routines. * They are able to read matrices from files and solve them. * * >>> Functions contained in this file: * main * ReadMatrixFromFile * PrintMatrixErrorMessage * InterpretCommandLine * GetProgramName * Usage */ /* * 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 lint static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; static char RCSid[] = "@(#)$Header: spTest.c,v 1.2 88/06/18 11:16:24 kundert Exp $"; #endif /* * IMPORTS * * >>> Import descriptions: * stdio.h math.h ctype.h * Standard C libraries. * spConfig.h * Macros that customize the sparse matrix package. It is not normally * necessary, nor is normally particularly desirable to include this * file into the calling routines. Nor should spINSIDE_SPARSE be defined. * It is done in this test file so that the complex test routines may be * removed when they don't exist in Sparse. * spMatrix.h * Macros and declarations to be imported by the user. */ #define YES 1 #define NO 0 #include #include #include #define spINSIDE_SPARSE #include "spConfig.h" #undef spINSIDE_SPARSE #include "spMatrix.h" /* Declarations that should be in stdio.h. */ extern char *malloc(); #ifdef ultrix extern void free(); extern void exit(); extern void perror(); #else extern free(); extern exit(); extern perror(); #endif /* Declarations that should be in strings.h. */ extern int strcmp(), strncmp(), strlen(); #ifdef lint #undef MULTIPLICATION #undef DETERMINANT #define MULTIPLICATION YES #define DETERMINANT YES #define LINT YES #else /* not lint */ #define LINT NO #endif /* not lint */ /* * TIMING */ #ifndef vms # include #endif #ifdef notdef /* __STDC__ */ /* Deleted because some ANSI C compilers to not yet have complete h files. */ # include # define HZ CLK_TCK #endif #ifndef HZ # ifdef vax # ifdef unix # define HZ 60 # endif # ifdef vms # define HZ 100 # endif # endif # ifdef hpux # ifdef hp300 # define HZ 50 # endif # ifdef hp500 # define HZ 60 # endif # endif #endif /* Routine that queries the system to find the process time. */ double Time() { struct time {long user, system, childuser, childsystem;} time; (void)times(&time); return (double)time.user / (double)HZ; } /* * MACROS */ #define BOOLEAN int #define NOT ! #define AND && #define OR || #define ALLOC(type,number) ((type *)malloc((unsigned)(sizeof(type)*(number)))) #define ABS(a) ((a) < 0.0 ? -(a) : (a)) #ifdef MAX #undef MAX #endif #ifdef MIN #undef MIN #endif #define MAX(a,b) ((a) > (b) ? (a) : (b)) #define MIN(a,b) ((a) < (b) ? (a) : (b)) /* * 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 #else #define IMAG_VECTORS #endif /* * GLOBAL DECLARATIONS * * These variables, types, and macros are used throughout this file. * * >>> Macros * PRINT_LIMIT * The maximum number of terms to be printed form the solution vector. * May be overridden by the -n option. * * >>> Variable types: * ElementRecord * A structure for holding data for each matrix element. * * >>> Global variables: * ProgramName (char *) * The name of the test program with any path stripped off. * FileName (char *) * The name of the current file name. * Complex (BOOLEAN) * The flag that indicates whether the matrix is complex or real. * Element (ElementRecord[]) * Array used to hold information about all the elements. * Matrix (char *) * The pointer to the matrix. * Size (int) * The size of the matrix. * RHS (RealVector) * The right-hand side vector (b in Ax = b). * Solution (RealVector) * The solution vector (x in Ax = b). * RHS_Verif (RealVector) * The calculated RHS vector. */ /* Begin global declarations. */ #define PRINT_LIMIT 9 typedef spREAL RealNumber, *RealVector; typedef struct { RealNumber Real; RealNumber Imag; } ComplexNumber, *ComplexVector; static char *ProgramName; static char *FileName; static int Size, MaxSize = 0; static int PrintLimit = PRINT_LIMIT, Iterations = 1, ColumnAsRHS = 1; static BOOLEAN Complex, SolutionOnly = NO, RealAsComplex = NO, Transposed = NO; static BOOLEAN CreatePlotFiles = NO, UseColumnAsRHS = NO, PrintLimitSet = NO; static BOOLEAN ExpansionWarningGiven, DiagPivoting = DIAGONAL_PIVOTING; static char Message[BUFSIZ], *Matrix = NULL; static RealNumber AbsThreshold = (RealNumber)0, RelThreshold = (RealNumber)0; static RealVector RHS, Solution, RHS_Verif; static RealVector iRHS, iSolution, iRHS_Verif; static ComplexVector cRHS, cSolution, cRHS_Verif; int Init(); /* * MAIN TEST ROUTINE for sparse matrix routines * * This routine reads a matrix from a file and solves it several * times. The solution is printed along with some statistics. * The format expected for the input matrix is the same as what is output by * spFileMatrix and spFileVector. * * >>> Local variables: * Determinant (RealNumber) * Real portion of the determinant of the matrix. * iDeterminant (RealNumber) * Imaginary portion of the determinant of the matrix. * Error (int) * Holds error status. * Last (int) * The index of the last term to be printed of the solution. * Iterations (int) * The number of times that the matrix will be factored and solved. * MaxRHS (RealNumber) * The largest term in the given RHS vector. * Residual (RealNumber) * The sum of the magnitude of the differences in each corresponding * term of the given and calculated RHS vector. * Exponent (int) * Exponent for the determinant. * Growth (RealNumber) * The growth that has occurred in the matrix during the factorization. * * >>> Timing variables: * BuildTime (double) * The time required to build up the matrix including the time to clear * the matrix. * ConditionTime (double) * The time required to compute the condition number. * DeterminantTime (double) * The time required to compute the determinant. * FactorTime (double) * The time required to factor the matrix without ordering. * InitialFactorTime (double) * The time required to factor the matrix with ordering. * SolveTime (double) * The time required to perform the forward and backward elimination. * StartTime (double) * The time that a timing interval was started. * * >>> Global variables: * Complex * Matrix * Size * RHS * iRHS */ main(ArgCount, Args) int ArgCount; char *Args[]; { register long I; int Error, Last, Elements, Fillins; RealNumber MaxRHS, Residual; RealNumber Determinant, iDeterminant; int Exponent; RealNumber ConditionNumber, PseudoCondition; RealNumber LargestBefore, LargestAfter, Roundoff, InfNorm; BOOLEAN StandardInput; char j, PlotFile[BUFSIZ], ErrMsg[BUFSIZ]; double StartTime, BeginTime, BuildTime, FactorTime, SolveTime, PartitionTime; double InitialFactorTime, ConditionTime, DeterminantTime; extern double Time(); extern char *sbrk(); /* Begin `main'. */ BeginTime = Time(); ArgCount = InterpretCommandLine( ArgCount, Args ); /* Assure that the Sparse is compatible with this test program.*/ # if NOT EXPANDABLE OR NOT INITIALIZE OR NOT ARRAY_OFFSET fprintf(stderr, "%s: Sparse is configured inappropriately for test program.\n", ProgramName); fprintf(stderr, " Enable EXPANDABLE, INITIALIZE, and ARRAY_OFFSET in `spConfig.h'.\n"); exit(1); # endif do { /* Initialization. */ BuildTime = FactorTime = SolveTime = 0.0; ExpansionWarningGiven = NO; /* Create matrix. */ Matrix = spCreate( 0, spCOMPLEX, &Error ); if (Matrix == NULL) { fprintf(stderr, "%s: insufficient memory available.\n", ProgramName); exit(1); } if( Error >= spFATAL ) goto End; /* Read matrix. */ if (ArgCount == 0) FileName = NULL; else FileName = *(++Args); Error = ReadMatrixFromFile(); if( Error) goto End; StandardInput = (FileName == NULL); /* Clear solution vector if row and column numbers are not densely packed. */ if (spGetSize(Matrix, YES) != spGetSize(Matrix, NO)) { if (Complex OR RealAsComplex) { for (I = Size; I > 0; I--) cSolution[I].Real = cSolution[I].Imag = 0.0; } else { for (I = Size; I > 0; I--) Solution[I] = 0.0; } } /* Perform initial build, factor, and solve. */ (void)spInitialize(Matrix, Init); #if MODIFIED_NODAL spMNA_Preorder( Matrix ); #endif #if DOCUMENTATION if (CreatePlotFiles) { if (StandardInput) (void)sprintf(PlotFile, "bef"); else (void)sprintf(PlotFile, "%s.bef", FileName); if (NOT spFileMatrix( Matrix, PlotFile, FileName, NO, NO, NO )) { (void)sprintf(ErrMsg,"%s: plotfile `%s'",ProgramName,PlotFile); perror( ErrMsg ); } } #if NO spPrint( Matrix, NO /*reodered*/, NO /*data*/, NO /*header*/ ); #endif #endif /* DOCUMENTATION */ #if STABILITY if (NOT SolutionOnly) LargestBefore = spLargestElement(Matrix); #endif #if CONDITION if (NOT SolutionOnly) InfNorm = spNorm(Matrix); #endif StartTime = Time(); Error = spOrderAndFactor( Matrix, RHS, RelThreshold, AbsThreshold, DiagPivoting ); InitialFactorTime = Time() - StartTime; PrintMatrixErrorMessage( Error ); if( Error >= spFATAL ) goto End; #if DOCUMENTATION if (CreatePlotFiles) { if (StandardInput) (void)sprintf(PlotFile, "aft"); else (void)sprintf(PlotFile, "%s.aft", FileName); if (NOT spFileMatrix( Matrix, PlotFile, FileName, YES, NO, NO )) { (void)sprintf(ErrMsg,"%s: plotfile `%s'",ProgramName,PlotFile); perror( ErrMsg ); } } #if NO spFileStats( Matrix, FileName, "stats" ); #endif #endif /* DOCUMENTATION */ /* * IMAG_VECTORS is a macro that replaces itself with `, iRHS, iSolution' * if the options spCOMPLEX and spSEPARATED_COMPLEX_VECTORS are set, * otherwise it disappears without a trace. */ #if TRANSPOSE if (Transposed) spSolveTransposed( Matrix, RHS, Solution IMAG_VECTORS ); else #endif spSolve( Matrix, RHS, Solution IMAG_VECTORS ); if (SolutionOnly) Iterations = 0; else { #if STABILITY LargestAfter = spLargestElement(Matrix); Roundoff = spRoundoff(Matrix, LargestAfter); #endif #if CONDITION StartTime = Time(); ConditionNumber = spCondition(Matrix, InfNorm, &Error); ConditionTime = Time() - StartTime; PrintMatrixErrorMessage( Error ); #endif #if PSEUDOCONDITION PseudoCondition = spPseudoCondition(Matrix); #endif StartTime = Time(); spPartition( Matrix, spDEFAULT_PARTITION ); PartitionTime = Time() - StartTime; } /* Solve system of equations Iterations times. */ for(I = 1; I <= Iterations; I++) { StartTime = Time(); (void)spInitialize(Matrix, Init); BuildTime += Time() - StartTime; StartTime = Time(); Error = spFactor( Matrix ); FactorTime += Time() - StartTime; if( Error != spOKAY ) PrintMatrixErrorMessage( Error ); if( Error >= spFATAL ) goto End; StartTime = Time(); /* * IMAG_VECTORS is a macro that replaces itself with `, iRHS, iSolution' * if the options spCOMPLEX and spSEPARATED_COMPLEX_VECTORS are set, * otherwise it disappears without a trace. */ #if TRANSPOSE if (Transposed) spSolveTransposed(Matrix, RHS, Solution IMAG_VECTORS ); else #endif spSolve( Matrix, RHS, Solution IMAG_VECTORS ); SolveTime += Time() - StartTime; } /* Print Solution. */ if (SolutionOnly) { if (PrintLimitSet) Last = MIN( PrintLimit, Size ); else Last = Size; j = ' '; } else { Last = (PrintLimit > Size) ? Size : PrintLimit; if (Last > 0) printf("Solution:\n"); j = 'j'; } if (Complex OR RealAsComplex) { #if spSEPARATED_COMPLEX_VECTORS for (I = 1; I <= Last; I++) { printf("%-16.9lg %-.9lg%c\n", (double)Solution[I], (double)iSolution[I], j); } #else for (I = 1; I <= Last; I++) { printf("%-16.9lg %-.9lg%c\n", (double)cSolution[I].Real, (double)cSolution[I].Imag, j); } #endif } else { for (I = 1; I <= Last; I++) printf("%-.9lg\n", (double)Solution[I]); } if (Last < Size AND Last != 0) printf("Solution list truncated.\n"); printf("\n"); #if DETERMINANT /* Calculate determinant. */ if (NOT SolutionOnly) { StartTime = Time(); #if spCOMPLEX spDeterminant( Matrix, &Exponent, &Determinant, &iDeterminant ); #else spDeterminant( Matrix, &Exponent, &Determinant ); #endif DeterminantTime = Time() - StartTime; if (Complex OR RealAsComplex) { Determinant = hypot( Determinant, iDeterminant ); while (Determinant >= 10.0) { Determinant *= 0.1; Exponent++; } } } #else Determinant = 0.0; Exponent = 0; #endif #if MULTIPLICATION if (NOT SolutionOnly) { /* Calculate difference between actual RHS vector and RHS vector * calculated from solution. */ /* Find the largest element in the given RHS vector. */ MaxRHS = 0.0; if (Complex OR RealAsComplex) { #if spSEPARATED_COMPLEX_VECTORS for (I = 1; I <= Size; I++) { if (ABS(RHS[I]) > MaxRHS) MaxRHS = ABS(RHS[I]); if (ABS(iRHS[I]) > MaxRHS) MaxRHS = ABS(iRHS[I]); } #else for (I = 1; I <= Size; I++) { if (ABS(cRHS[I].Real) > MaxRHS) MaxRHS = ABS(cRHS[I].Real); if (ABS(cRHS[I].Imag) > MaxRHS) MaxRHS = ABS(cRHS[I].Imag); } #endif } else { for (I = 1; I <= Size; I++) { if (ABS(RHS[I]) > MaxRHS) MaxRHS = ABS(RHS[I]); } } /* Rebuild matrix. */ (void)spInitialize(Matrix, Init); #if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS if (Transposed) { spMultTransposed( Matrix, RHS_Verif, Solution, iRHS_Verif, iSolution); } else spMultiply(Matrix, RHS_Verif, Solution, iRHS_Verif, iSolution); #else if (Transposed) spMultTransposed( Matrix, RHS_Verif, Solution ); else spMultiply( Matrix, RHS_Verif, Solution ); #endif /* Calculate residual. */ Residual = 0.0; if (Complex OR RealAsComplex) { #if spSEPARATED_COMPLEX_VECTORS for (I = 1; I <= Size; I++) { Residual += ABS(RHS[I] - RHS_Verif[I]) + ABS(iRHS[I] - iRHS_Verif[I]); } #else for (I = 1; I <= Size; I++) { Residual += ABS(cRHS[I].Real - cRHS_Verif[I].Real) + ABS(cRHS[I].Imag - cRHS_Verif[I].Imag); } #endif } else { for (I = 1; I <= Size; I++) Residual += ABS(RHS[I] - RHS_Verif[I]); } } #endif /* Print statistics. */ if (NOT SolutionOnly) { Elements = spElementCount(Matrix); Fillins = spFillinCount(Matrix); printf("Initial factor time = %.2lf.\n", InitialFactorTime); printf("Partition time = %.2lf.\n", PartitionTime); if (Iterations > 0) { printf("Build time = %.3lf.\n", (BuildTime/Iterations)); printf("Factor time = %.3lf.\n",(FactorTime/Iterations)); printf("Solve time = %.3lf.\n", (SolveTime/Iterations)); } #if STABILITY printf("Condition time = %.2lf.\n", ConditionTime); #endif #if DETERMINANT printf("Determinant time = %.2lf.\n", DeterminantTime); #endif printf("\nTotal number of elements = %d.\n", Elements); printf("Average number of elements per row initially = %.2lf.\n", (double)(Elements - Fillins) / (double)spGetSize(Matrix, NO)); printf("Total number of fill-ins = %d.\n", Fillins); #if DETERMINANT OR MULTIPLICATION OR PSEUDOCONDITION OR CONDITION OR STABILITY putchar('\n'); #endif #if STABILITY if (LargestBefore != 0.0) printf("Growth = %.2lg.\n", LargestAfter / LargestBefore); printf("Max error in matrix = %.2lg.\n", Roundoff); #endif #if STABILITY if(ABS(ConditionNumber) > 10 * SMALLEST_REAL); printf("Condition number = %.2lg.\n", 1.0 / ConditionNumber); #endif #if CONDITION AND STABILITY printf("Estimated upper bound of error in solution = %.2lg.\n", Roundoff / ConditionNumber); #endif #if PSEUDOCONDITION printf("PseudoCondition = %.2lg.\n", PseudoCondition); #endif #if DETERMINANT printf("Determinant = %.3lg", (double)Determinant ); if (Determinant != 0.0 AND Exponent != 0) printf("e%d", Exponent); putchar('.'); putchar('\n'); #endif #if MULTIPLICATION if (MaxRHS != 0.0) printf("Normalized residual = %.2lg.\n", (Residual / MaxRHS)); #endif } End:; (void)fflush(stdout); CheckInAllComplexArrays(); spDestroy(Matrix); Matrix = NULL; } while( --ArgCount > 0 ); if (NOT SolutionOnly) { printf("\nAggregate resource usage:\n"); printf(" Time required = %.2lf seconds.\n", Time() - BeginTime); printf(" Virtual memory used = %d kBytes.\n\n", ((int)sbrk(0))/1000); } exit (0); } /* * READ MATRIX FROM FILE * * This function reads the input file for the matrix and the RHS vector. * If no RHS vector exists, one is created. If there is an error in the * file, the appropriate error messages are delivered to standard output. * * >>> Returned: * The error status is returned. If no error occurred, a zero is returned. * Otherwise, a one is returned. * * >>> Local variables: * pMatrixFile (FILE *) * The pointer to the file that holds the matrix. * InputString (char []) * String variable for holding input from the matrix file. * Message (char []) * String variable that contains a one line descriptor of the matrix. * Descriptor is taken from matrix file. * * >>> Global variables: * Complex * Size * Element * RHS * iRHS */ static int ReadMatrixFromFile() { long I, Reads; FILE *pMatrixFile; char sInput[BUFSIZ], sType[BUFSIZ], *p; int Error, Row, Col, Count = 0, LineNumber, RHS_Col, IntSize; double Real, Imag = 0.0; ComplexNumber *pValue, *pInitInfo, *CheckOutComplexArray(); RealNumber *pElement; static char *EndOfFile = "%s: unexpected end of file `%s' at line %d.\n"; static char *Syntax = "%s: syntax error in file `%s' at line %d.\n"; /* Begin `ReadMatrixFromFile'. */ /* Open matrix file in read mode. */ if (NOT SolutionOnly) putchar('\n'); if (FileName == NULL) { FileName = "standard input"; pMatrixFile = stdin; } else { pMatrixFile = fopen(FileName, "r"); if (pMatrixFile == NULL) { fprintf(stderr, "%s: file %s was not found.\n", ProgramName, FileName); return 1; } } Complex = NO; LineNumber = 1; /* Read and print label. */ if (NULL == fgets( Message, BUFSIZ, pMatrixFile )) { fprintf(stderr, EndOfFile, ProgramName, FileName, LineNumber); return 1; } /* For compatibility with the old file syntax. */ if (!strncmp( Message, "Starting", 8 )) { /* Test for complex matrix. */ if (strncmp( Message, "Starting complex", 15 ) == 0) Complex = YES; LineNumber++; if (NULL == fgets( Message, BUFSIZ, pMatrixFile )) { fprintf(stderr, EndOfFile, ProgramName, FileName, LineNumber); return 1; } } if (NOT SolutionOnly) printf("%-s\n", Message); /* Read size of matrix and determine type of matrix. */ LineNumber++; if (NULL == fgets( sInput, BUFSIZ, pMatrixFile )) { fprintf(stderr, EndOfFile, ProgramName, FileName, LineNumber); return 1; } if ((Reads = sscanf( sInput,"%d %s", &Size, sType )) < 1) { fprintf(stderr, Syntax, ProgramName, FileName, LineNumber); return 1; } if (Reads == 2) { for (p = sType; *p != '\0'; p++) if (isupper(*p)) *p += 'a'-'A'; if (strncmp( sType, "complex", 7 ) == 0) Complex = YES; else if (strncmp( sType, "real", 7 ) == 0) Complex = NO; else { fprintf(stderr, Syntax, ProgramName, FileName, LineNumber); return 1; } } EnlargeVectors( Size, YES ); RHS_Col = MIN( Size, ColumnAsRHS ); #if NOT spCOMPLEX if (Complex) { fprintf(stderr, "%s: Sparse is not configured to solve complex matrices.\n", ProgramName); fprintf(stderr," Enable spCOMPLEX in `spConfig.h'.\n"); return 1; } #endif #if NOT REAL if (NOT (Complex OR RealAsComplex)) { fprintf(stderr, "%s: Sparse is not configured to solve real matrices.\n", ProgramName); fprintf(stderr," Enable REAL in `spConfig.h'.\n"); return 1; } #endif /* Read matrix elements. */ do { if (Count == 0) pValue = CheckOutComplexArray( Count = 1000 ); LineNumber++; if (NULL == fgets( sInput, BUFSIZ, pMatrixFile )) { fprintf(stderr, "%s: unexpected end of file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } if (Complex) { Reads = sscanf( sInput,"%d%d%lf%lf", &Row, &Col, &Real, &Imag ); if (Reads != 4) { fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } } else { Reads = sscanf( sInput,"%d%d%lf", &Row, &Col, &Real ); if (Reads != 3) { fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } } if(Row < 0 OR Col < 0) { fprintf(stderr, "%s: index not positive in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } if(Row > Size OR Col > Size) { if (NOT ExpansionWarningGiven) { fprintf( stderr, "%s: computed and given matrix size differ in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); ExpansionWarningGiven = YES; } Size = MAX(Row, Col); EnlargeVectors( Size, NO ); } pElement = spGetElement( Matrix, Row, Col ); if (pElement == NULL) { fprintf(stderr, "%s: insufficient memory available.\n", ProgramName); exit(1); } pInitInfo = (ComplexNumber *)spGetInitInfo( pElement ); if (pInitInfo == NULL) { pValue[--Count].Real = Real; pValue[Count].Imag = Imag; spInstallInitInfo( pElement, (char *)(pValue + Count) ); } else { pInitInfo->Real += Real; pInitInfo->Imag += Imag; } /* Save into RHS vector if in desired column. */ if (Col == RHS_Col) { if (Complex OR RealAsComplex) { #if spSEPARATED_COMPLEX_VECTORS RHS[Row] = Real; iRHS[Row] = Imag; #else cRHS[Row].Real = Real; cRHS[Row].Imag = Imag; #endif } else RHS[Row] = Real; } } while (Row != 0 AND Col != 0); Size = spGetSize( Matrix, YES ); if (Error = spError( Matrix ) != spOKAY) { PrintMatrixErrorMessage( Error ); if( Error >= spFATAL ) return 1; } /* Read RHS vector. */ if (NOT UseColumnAsRHS AND (NULL != fgets( sInput, BUFSIZ, pMatrixFile ))) { /* RHS vector exists, read it. */ LineNumber++; for (I = 1; I <= Size; I++) { if (I != 1 OR (strncmp( sInput, "Beginning", 8 ) == 0)) { LineNumber++; if (NULL == fgets( sInput, BUFSIZ, pMatrixFile )) { fprintf(stderr, "%s: unexpected end of file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } } if (Complex) { #if spSEPARATED_COMPLEX_VECTORS Reads = sscanf( sInput,"%lf%lf", &RHS[I], &iRHS[I] ); #else Reads = sscanf( sInput, "%lf%lf", &cRHS[I].Real, &cRHS[I].Imag ); #endif if (Reads != 2) { fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } } else /* Not complex. */ { Reads = sscanf( sInput, "%lf", &RHS[I] ); if (Reads != 1) { fprintf(stderr, "%s: syntax error in file `%s' at line %d.\n", ProgramName, FileName, LineNumber); return 1; } } } if (RealAsComplex AND NOT Complex) { #if spSEPARATED_COMPLEX_VECTORS for (I = 1; I <= Size; I++) iRHS[I] = 0.0; #else for (I = Size; I > 0; I--) { cRHS[I].Real = RHS[I]; cRHS[I].Imag = 0.0; } #endif } } /* Print out the size and the type of the matrix. */ if (NOT SolutionOnly) { IntSize = spGetSize( Matrix, NO ); printf("Matrix is %d x %d ", IntSize, IntSize); if (IntSize != Size) printf("(external size is %d x %d) ", Size, Size); if (Complex OR RealAsComplex) printf("and complex.\n",Size,Size); else printf("and real.\n",Size,Size); } if (Complex OR RealAsComplex) spSetComplex( Matrix ); else spSetReal( Matrix ); (void)fclose(pMatrixFile); return 0; } /* * INITIALIZE MATRIX ELEMENT * * This function copys the InitInfo to the Real and Imag matrix element * values. * * >>> Returned: * A zero is returns, signifying that no error can be made. */ /*ARGSUSED*/ static int Init( pElement, pInitInfo, Row, Col) RealNumber *pElement; char *pInitInfo; int Row, Col; { /* Begin `Init'. */ *pElement = *(RealNumber *)pInitInfo; if (Complex OR RealAsComplex) *(pElement+1) = *((RealNumber *)pInitInfo+1); return 0; } /* * COMPLEX ARRAY ALLOCATION * * These routines are used to check out and in arrays of complex numbers. */ static struct Bin { ComplexVector Array; struct Bin *Next; } *FirstArray, *CurrentArray = NULL; static ComplexVector CheckOutComplexArray( Count ) int Count; { struct Bin Bin, *pBin; ComplexVector Temp; /* Begin `CheckOutComplexArray'. */ if (CurrentArray == NULL OR CurrentArray->Next == NULL) { pBin = ALLOC( struct Bin, 1); Temp = ALLOC( ComplexNumber, Count ); if (pBin == NULL OR Temp == NULL) { fprintf( stderr, "%s: insufficient memory available.\n", ProgramName); exit(1); } Bin.Array = Temp; Bin.Next = NULL; *pBin = Bin; if (CurrentArray == NULL) FirstArray = CurrentArray = pBin; else { CurrentArray->Next = pBin; CurrentArray = pBin; } } else { pBin = CurrentArray; CurrentArray = pBin->Next; } return pBin->Array; } static CheckInAllComplexArrays() { /* Begin `CheckInAllComplexArrays'. */ if (CurrentArray != NULL) CurrentArray = FirstArray; return; } /* * PRINT ERROR MESSAGE * * Prints error message for matrix routines if one is required. * * >>> Argument: * Error (int) * The error status of the matrix routines. */ static PrintMatrixErrorMessage( Error ) int Error; { int Row, Col; /* Begin `PrintMatrixErrorMessage'. */ if (Error == spOKAY) return; if (Error >= spFATAL) { fprintf(stderr, "%s: fatal error detected in file `%s'.\n", ProgramName, FileName ); } else fprintf(stderr, "%s: warning in `%s'.\n", ProgramName, FileName ); if (Error == spNO_MEMORY) fprintf(stderr, " Insufficient memory available.\n"); else if (Error == spSINGULAR) { spWhereSingular( Matrix, &Row, &Col ); printf(" Singular matrix (detected at row %d, column %d).\n", Row, Col ); } else if (Error == spZERO_DIAG) { spWhereSingular( Matrix, &Row, &Col ); printf(" Zero diagonal detected at row %d, column %d.\n", Row, Col ); } else if (Error == spPANIC) fprintf(stderr, " Matrix routines being used improperly.\n"); else if (Error == spSMALL_PIVOT) fprintf(stderr, " A pivot was chosen that is smaller than the absolute threshold.\n"); return; } /* * INTERPRET THE COMMAND LINE ARGUMENTS */ static int InterpretCommandLine( ArgCount, Args ) int ArgCount; char *Args[]; { int I, FileCount = 0; char *GetProgramName(); /* Begin `InterpretCommandLine'. */ /* Determine the name of the program. */ ProgramName = GetProgramName(Args[0]); /* Step through the argument list, interpreting and deleting the options. */ for (I = 1; I < ArgCount; I++) { if (!strcmp(Args[I], "-a")) { if (ArgCount == I OR (sscanf(Args[I+1],"%lf", &AbsThreshold) != 1)) { AbsThreshold = 0.0; Usage(Args[I]); } else I++; } else if (!strcmp(Args[I], "-r")) { if (ArgCount == I OR (sscanf(Args[I+1],"%lf", &RelThreshold) != 1)) { RelThreshold = 0.0; Usage(Args[I]); } else I++; } else if (!strcmp(Args[I], "-x")) { #if spCOMPLEX RealAsComplex = YES; #else fprintf(stderr, "%s: Sparse is not configured to solve complex matrices.\n", ProgramName); fprintf(stderr," Enable spCOMPLEX in `spConfig.h'.\n"); #endif } else if (!strcmp(Args[I], "-s")) SolutionOnly = YES; else if (!strcmp(Args[I], "-c")) DiagPivoting = NO; else if (!strcmp(Args[I], "-t")) { #if TRANSPOSE Transposed = YES; #else fprintf(stderr, "%s: Sparse is not configured to solve transposed system.\n", ProgramName); fprintf(stderr," Enable TRANSPOSE in `spConfig.h'.\n"); #endif } else if (!strcmp(Args[I], "-n")) { if (ArgCount == I OR (sscanf(Args[I+1],"%ld", &PrintLimit) != 1)) { PrintLimit = PRINT_LIMIT; Usage(Args[I]); } else { PrintLimitSet = YES; I++; } } else if (!strcmp(Args[I], "-i")) { if (ArgCount == I OR (sscanf(Args[I+1],"%ld", &Iterations) != 1)) { Iterations = 1; Usage(Args[I]); } else I++; } else if (!strcmp(Args[I], "-b")) { if (ArgCount == I OR (sscanf(Args[I+1],"%ld", &ColumnAsRHS) != 1)) { ColumnAsRHS = 1; UseColumnAsRHS = YES; } else { UseColumnAsRHS = YES; ColumnAsRHS = MAX( ColumnAsRHS, 1 ); I++; } } else if (!strcmp(Args[I], "-p")) { #if DOCUMENTATION CreatePlotFiles = YES; #else fprintf(stderr, "%s: Sparse is not configured to generate plot files.\n", ProgramName); fprintf(stderr," Enable DOCUMENTATION in `spConfig.h'.\n"); #endif } else if (!strcmp(Args[I], "-u")) Usage( (char *)NULL ); else if (Args[I][0] == '-') Usage(Args[I]); else Args[++FileCount] = Args[I]; } return FileCount; } /* * PROGRAM NAME * * Removes path from argv[0] and returns program name. * Assumes UNIX style path names. */ static char * GetProgramName( Arg0 ) char *Arg0; { char *pTail; /* Begin `GetProgramName'. */ pTail = Arg0 + strlen(Arg0)-1; while (pTail != Arg0 AND *(pTail-1) != '/') --pTail; return pTail; } /* * USAGE WARNING * * Sends a warning to the user when the command line syntax is wrong. */ static Usage( sBadOpt ) char *sBadOpt; { /* Begin `Usage'. */ if (sBadOpt != NULL) { fprintf(stderr, "%s: unknown or deformed option `%s'.\n", ProgramName, sBadOpt); } fprintf(stderr, "\nUsage: %s [options] [file names]\n\n", ProgramName); fprintf(stderr, "Options:\n"); fprintf(stderr, " -s Print solution rather than run statistics.\n"); fprintf( stderr, " -r x Use x as relative threshold.\n"); fprintf( stderr, " -a x Use x as absolute threshold.\n"); fprintf( stderr, " -n n Print first n terms of solution vector.\n"); fprintf( stderr, " -i n Repeat build/factor/solve n times.\n"); fprintf( stderr, " -b n Use n'th column of matrix as b in Ax=b.\n"); #if DIAGONAL_PIVOTING fprintf( stderr, " -c Use complete (as opposed to diagonal) pivoting.\n"); #endif #if DOCUMENTATION fprintf( stderr, " -p Create plot files `filename.bef' and `filename.aft'.\n"); #endif #if spCOMPLEX fprintf( stderr, " -x Treat real matrix as complex with imaginary part zero.\n"); #endif #if TRANSPOSE fprintf( stderr, " -t Solve transposed system.\n"); #endif fprintf( stderr, " -u Print usage message.\n"); exit(1); } /* * ENLARGE VECTORS * * Allocate or enlarge vectors. */ static EnlargeVectors( NewSize, Clear ) int NewSize, Clear; { int I, PrevSize = MaxSize; RealVector OldRHS = RHS, iOldRHS = iRHS; ComplexVector cOldRHS = cRHS; #if spCOMPLEX # define SCALE 2 #else # define SCALE 1 #endif /* Begin `EnlargeVectors'. */ if (NewSize > PrevSize) { if (MaxSize != 0) { free( (char *)Solution ); free( (char *)RHS_Verif ); } RHS = ALLOC( RealNumber, SCALE*(NewSize+1) ); Solution = ALLOC( RealNumber, SCALE*(NewSize+1) ); RHS_Verif = ALLOC( RealNumber, SCALE*(NewSize+1) ); if (NOT RHS OR NOT Solution OR NOT RHS_Verif) { fprintf(stderr, "%s: insufficient memory available.\n", ProgramName); exit(1); } cRHS = (ComplexVector)RHS; cSolution = (ComplexVector)Solution; cRHS_Verif = (ComplexVector)RHS_Verif; iRHS = RHS + NewSize + 1; iSolution = Solution + NewSize + 1; iRHS_Verif = RHS_Verif + NewSize + 1; /* Copy data from old RHS to new RHS. */ if (NOT Clear) { /* Copy old RHS vector to new. */ if (Complex OR RealAsComplex) { for (I = PrevSize; I > 0; I--) { #if spSEPARATED_COMPLEX_VECTORS OR LINT RHS[I] = OldRHS[I]; iRHS[I] = iOldRHS[I]; #endif #if spSEPARATED_COMPLEX_VECTORS OR LINT cRHS[I] = cOldRHS[I]; #endif } } else { for (I = PrevSize; I > 0; I--) RHS[I] = OldRHS[I]; } } if (MaxSize != 0) free( (char *)OldRHS ); MaxSize = NewSize; } /* Either completely clear or clear remaining portion of RHS vector. */ if ((NewSize > PrevSize) OR Clear) { if (Clear) { NewSize = MAX( NewSize, PrevSize ); PrevSize = 0; } if (Complex OR RealAsComplex) { for (I = NewSize; I > PrevSize; I--) { #if spSEPARATED_COMPLEX_VECTORS OR LINT RHS[I] = 0.0; iRHS[I] = 0.0; #endif #if NOT spSEPARATED_COMPLEX_VECTORS OR LINT cRHS[I].Real = 0.0; cRHS[I].Imag = 0.0; #endif } } else { for (I = NewSize; I > PrevSize; I--) RHS[I] = 0.0; } } return; } @EOF chmod 644 spTest.c echo x - spUtils.c cat >spUtils.c <<'@EOF' /* * MATRIX UTILITY MODULE * * Author: Advising professor: * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli * UC Berkeley * * This file contains various optional utility routines. * * >>> User accessible functions contained in this file: * spMNA_Preorder * spScale * spMultiply * spMultTransposed * spDeterminant * spStrip * spDeleteRowAndCol * spPseudoCondition * spCondition * spNorm * spLargestElement * spRoundoff * * >>> Other functions contained in this file: * CountTwins * SwapCols * ScaleComplexMatrix * ComplexMatrixMultiply * ComplexCondition */ /* * 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 lint static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; static char RCSid[] = "@(#)$Header: spUtils.c,v 1.2 88/06/18 11:16:48 kundert Exp $"; #endif /* * IMPORTS * * >>> Import descriptions: * spConfig.h * Macros that customize the sparse matrix routines. * spMatrix.h * Macros and declarations to be imported by the user. * spDefs.h * Matrix type and macro definitions for the sparse matrix routines. */ #define spINSIDE_SPARSE #include "spConfig.h" #include "spMatrix.h" #include "spDefs.h" #if MODIFIED_NODAL /* * PREORDER MODIFIED NODE ADMITTANCE MATRIX TO REMOVE ZEROS FROM DIAGONAL * * This routine massages modified node admittance matrices to remove * zeros from the diagonal. It takes advantage of the fact that the * row and column associated with a zero diagonal usually have * structural ones placed symmetricly. This routine should be used * only on modified node admittance matrices and should be executed * after the matrix has been built but before the factorization * begins. It should be executed for the initial factorization only * and should be executed before the rows have been linked. Thus it * should be run before using spScale(), spMultiply(), * spDeleteRowAndCol(), or spNorm(). * * This routine exploits the fact that the structural ones are placed * in the matrix in symmetric twins. For example, the stamps for * grounded and a floating voltage sources are * grounded: floating: * [ x x 1 ] [ x x 1 ] * [ x x ] [ x x -1 ] * [ 1 ] [ 1 -1 ] * Notice for the grounded source, there is one set of twins, and for * the floating, there are two sets. We remove the zero from the diagonal * by swapping the rows associated with a set of twins. For example: * grounded: floating 1: floating 2: * [ 1 ] [ 1 -1 ] [ x x 1 ] * [ x x ] [ x x -1 ] [ 1 -1 ] * [ x x 1 ] [ x x 1 ] [ x x -1 ] * * It is important to deal with any zero diagonals that only have one * set of twins before dealing with those that have more than one because * swapping row destroys the symmetry of any twins in the rows being * swapped, which may limit future moves. Consider * [ x x 1 ] * [ x x -1 1 ] * [ 1 -1 ] * [ 1 ] * There is one set of twins for diagonal 4 and two for diagonal 3. * Dealing with diagonal 4 first requires swapping rows 2 and 4. * [ x x 1 ] * [ 1 ] * [ 1 -1 ] * [ x x -1 1 ] * We can now deal with diagonal 3 by swapping rows 1 and 3. * [ 1 -1 ] * [ 1 ] * [ x x 1 ] * [ x x -1 1 ] * And we are done, there are no zeros left on the diagonal. However, if * we originally dealt with diagonal 3 first, we could swap rows 2 and 3 * [ x x 1 ] * [ 1 -1 ] * [ x x -1 1 ] * [ 1 ] * Diagonal 4 no longer has a symmetric twin and we cannot continue. * * So we always take care of lone twins first. When none remain, we * choose arbitrarily a set of twins for a diagonal with more than one set * and swap the rows corresponding to that twin. We then deal with any * lone twins that were created and repeat the procedure until no * zero diagonals with symmetric twins remain. * * In this particular implementation, columns are swapped rather than rows. * The algorithm used in this function was developed by Ken Kundert and * Tom Quarles. * * >>> Arguments: * eMatrix (char *) * Pointer to the matrix to be preordered. * * >>> Local variables; * J (int) * Column with zero diagonal being currently considered. * pTwin1 (ElementPtr) * Pointer to the twin found in the column belonging to the zero diagonal. * pTwin2 (ElementPtr) * Pointer to the twin found in the row belonging to the zero diagonal. * belonging to the zero diagonal. * AnotherPassNeeded (BOOLEAN) * Flag indicating that at least one zero diagonal with symmetric twins * remain. * StartAt (int) * Column number of first zero diagonal with symmetric twins. * Swapped (BOOLEAN) * Flag indicating that columns were swapped on this pass. * Twins (int) * Number of symmetric twins corresponding to current zero diagonal. */ void spMNA_Preorder( eMatrix ) char *eMatrix; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register int J, Size; ElementPtr pTwin1, pTwin2; int Twins, StartAt = 1; BOOLEAN Swapped, AnotherPassNeeded; /* Begin `spMNA_Preorder'. */ ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored ); if (Matrix->RowsLinked) return; Size = Matrix->Size; Matrix->Reordered = YES; do { AnotherPassNeeded = Swapped = NO; /* Search for zero diagonals with lone twins. */ for (J = StartAt; J <= Size; J++) { if (Matrix->Diag[J] == NULL) { Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 ); if (Twins == 1) { /* Lone twins found, swap rows. */ SwapCols( Matrix, pTwin1, pTwin2 ); Swapped = YES; } else if ((Twins > 1) AND NOT AnotherPassNeeded) { AnotherPassNeeded = YES; StartAt = J; } } } /* All lone twins are gone, look for zero diagonals with multiple twins. */ if (AnotherPassNeeded) { for (J = StartAt; NOT Swapped AND (J <= Size); J++) { if (Matrix->Diag[J] == NULL) { Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 ); SwapCols( Matrix, pTwin1, pTwin2 ); Swapped = YES; } } } } while (AnotherPassNeeded); return; } /* * COUNT TWINS * * This function counts the number of symmetric twins associated with * a zero diagonal and returns one set of twins if any exist. The * count is terminated early at two. */ static int CountTwins( Matrix, Col, ppTwin1, ppTwin2 ) MatrixPtr Matrix; int Col; ElementPtr *ppTwin1, *ppTwin2; { int Row, Twins = 0; ElementPtr pTwin1, pTwin2; /* Begin `CountTwins'. */ pTwin1 = Matrix->FirstInCol[Col]; while (pTwin1 != NULL) { if (ABS(pTwin1->Real) == 1.0) { Row = pTwin1->Row; pTwin2 = Matrix->FirstInCol[Row]; while ((pTwin2 != NULL) AND (pTwin2->Row != Col)) pTwin2 = pTwin2->NextInCol; if ((pTwin2 != NULL) AND (ABS(pTwin2->Real) == 1.0)) { /* Found symmetric twins. */ if (++Twins >= 2) return Twins; (*ppTwin1 = pTwin1)->Col = Col; (*ppTwin2 = pTwin2)->Col = Row; } } pTwin1 = pTwin1->NextInCol; } return Twins; } /* * SWAP COLUMNS * * This function swaps two columns and is applicable before the rows are * linked. */ static SwapCols( Matrix, pTwin1, pTwin2 ) MatrixPtr Matrix; ElementPtr pTwin1, pTwin2; { int Col1 = pTwin1->Col, Col2 = pTwin2->Col; /* Begin `SwapCols'. */ SWAP (ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]); SWAP (int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]); #if TRANSLATE Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col2]]=Col2; Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col1]]=Col1; #endif Matrix->Diag[Col1] = pTwin2; Matrix->Diag[Col2] = pTwin1; Matrix->NumberOfInterchangesIsOdd = NOT Matrix->NumberOfInterchangesIsOdd; return; } #endif /* MODIFIED_NODAL */ #if SCALING /* * SCALE MATRIX * * This function scales the matrix to enhance the possibility of * finding a good pivoting order. Note that scaling enhances accuracy * of the solution only if it affects the pivoting order, so it makes * no sense to scale the matrix before spFactor(). If scaling is * desired it should be done before spOrderAndFactor(). There * are several things to take into account when choosing the scale * factors. First, the scale factors are directly multiplied against * the elements in the matrix. To prevent roundoff, each scale factor * should be equal to an integer power of the number base of the * machine. Since most machines operate in base two, scale factors * should be a power of two. Second, the matrix should be scaled such * that the matrix of element uncertainties is equilibrated. Third, * this function multiplies the scale factors by the elements, so if * one row tends to have uncertainties 1000 times smaller than the * other rows, then its scale factor should be 1024, not 1/1024. * Fourth, to save time, this function does not scale rows or columns * if their scale factors are equal to one. Thus, the scale factors * should be normalized to the most common scale factor. Rows and * columns should be normalized separately. For example, if the size * of the matrix is 100 and 10 rows tend to have uncertainties near * 1e-6 and the remaining 90 have uncertainties near 1e-12, then the * scale factor for the 10 should be 1/1,048,576 and the scale factors * for the remaining 90 should be 1. Fifth, since this routine * directly operates on the matrix, it is necessary to apply the scale * factors to the RHS and Solution vectors. It may be easier to * simply use spOrderAndFactor() on a scaled matrix to choose the * pivoting order, and then throw away the matrix. Subsequent * factorizations, performed with spFactor(), will not need to have * the RHS and Solution vectors descaled. Lastly, this function * should not be executed before the function spMNA_Preorder. * * >>> Arguments: * eMatrix (char *) * Pointer to the matrix to be scaled. * SolutionScaleFactors (RealVector) * The array of Solution scale factors. These factors scale the columns. * All scale factors are real valued. * RHS_ScaleFactors (RealVector) * The array of RHS scale factors. These factors scale the rows. * All scale factors are real valued. * * >>> Local variables: * lSize (int) * Local version of the size of the matrix. * pElement (ElementPtr) * Pointer to an element in the matrix. * pExtOrder (int *) * Pointer into either IntToExtRowMap or IntToExtColMap vector. Used to * compensate for any row or column swaps that have been performed. * ScaleFactor (RealNumber) * The scale factor being used on the current row or column. */ void spScale( eMatrix, RHS_ScaleFactors, SolutionScaleFactors ) char *eMatrix; register RealVector RHS_ScaleFactors, SolutionScaleFactors; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register ElementPtr pElement; register int I, lSize, *pExtOrder; RealNumber ScaleFactor; void ScaleComplexMatrix(); /* Begin `spScale'. */ ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored ); if (NOT Matrix->RowsLinked) spcLinkRows( Matrix ); #if spCOMPLEX if (Matrix->Complex) { ScaleComplexMatrix( Matrix, RHS_ScaleFactors, SolutionScaleFactors ); return; } #endif #if REAL lSize = Matrix->Size; /* Correct pointers to arrays for ARRAY_OFFSET */ #if NOT ARRAY_OFFSET --RHS_ScaleFactors; --SolutionScaleFactors; #endif /* Scale Rows */ pExtOrder = &Matrix->IntToExtRowMap[1]; for (I = 1; I <= lSize; I++) { if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0) { pElement = Matrix->FirstInRow[I]; while (pElement != NULL) { pElement->Real *= ScaleFactor; pElement = pElement->NextInRow; } } } /* Scale Columns */ pExtOrder = &Matrix->IntToExtColMap[1]; for (I = 1; I <= lSize; I++) { if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { pElement->Real *= ScaleFactor; pElement = pElement->NextInCol; } } } return; #endif /* REAL */ } #endif /* SCALING */ #if spCOMPLEX AND SCALING /* * SCALE COMPLEX MATRIX * * This function scales the matrix to enhance the possibility of * finding a good pivoting order. Note that scaling enhances accuracy * of the solution only if it affects the pivoting order, so it makes * no sense to scale the matrix before spFactor(). If scaling is * desired it should be done before spOrderAndFactor(). There * are several things to take into account when choosing the scale * factors. First, the scale factors are directly multiplied against * the elements in the matrix. To prevent roundoff, each scale factor * should be equal to an integer power of the number base of the * machine. Since most machines operate in base two, scale factors * should be a power of two. Second, the matrix should be scaled such * that the matrix of element uncertainties is equilibrated. Third, * this function multiplies the scale factors by the elements, so if * one row tends to have uncertainties 1000 times smaller than the * other rows, then its scale factor should be 1024, not 1/1024. * Fourth, to save time, this function does not scale rows or columns * if their scale factors are equal to one. Thus, the scale factors * should be normalized to the most common scale factor. Rows and * columns should be normalized separately. For example, if the size * of the matrix is 100 and 10 rows tend to have uncertainties near * 1e-6 and the remaining 90 have uncertainties near 1e-12, then the * scale factor for the 10 should be 1/1,048,576 and the scale factors * for the remaining 90 should be 1. Fifth, since this routine * directly operates on the matrix, it is necessary to apply the scale * factors to the RHS and Solution vectors. It may be easier to * simply use spOrderAndFactor() on a scaled matrix to choose the * pivoting order, and then throw away the matrix. Subsequent * factorizations, performed with spFactor(), will not need to have * the RHS and Solution vectors descaled. Lastly, this function * should not be executed before the function spMNA_Preorder. * * >>> Arguments: * Matrix (char *) * Pointer to the matrix to be scaled. * SolutionScaleFactors (RealVector) * The array of Solution scale factors. These factors scale the columns. * All scale factors are real valued. * RHS_ScaleFactors (RealVector) * The array of RHS scale factors. These factors scale the rows. * All scale factors are real valued. * * >>> Local variables: * lSize (int) * Local version of the size of the matrix. * pElement (ElementPtr) * Pointer to an element in the matrix. * pExtOrder (int *) * Pointer into either IntToExtRowMap or IntToExtColMap vector. Used to * compensate for any row or column swaps that have been performed. * ScaleFactor (RealNumber) * The scale factor being used on the current row or column. */ static void ScaleComplexMatrix( Matrix, RHS_ScaleFactors, SolutionScaleFactors ) MatrixPtr Matrix; register RealVector RHS_ScaleFactors, SolutionScaleFactors; { register ElementPtr pElement; register int I, lSize, *pExtOrder; RealNumber ScaleFactor; /* Begin `ScaleComplexMatrix'. */ lSize = Matrix->Size; /* Correct pointers to arrays for ARRAY_OFFSET */ #if NOT ARRAY_OFFSET --RHS_ScaleFactors; --SolutionScaleFactors; #endif /* Scale Rows */ pExtOrder = &Matrix->IntToExtRowMap[1]; for (I = 1; I <= lSize; I++) { if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0) { pElement = Matrix->FirstInRow[I]; while (pElement != NULL) { pElement->Real *= ScaleFactor; pElement->Imag *= ScaleFactor; pElement = pElement->NextInRow; } } } /* Scale Columns */ pExtOrder = &Matrix->IntToExtColMap[1]; for (I = 1; I <= lSize; I++) { if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { pElement->Real *= ScaleFactor; pElement->Imag *= ScaleFactor; pElement = pElement->NextInCol; } } } return; } #endif /* SCALING AND spCOMPLEX */ #if MULTIPLICATION /* * MATRIX MULTIPLICATION * * Multiplies matrix by solution vector to find source vector. * Assumes matrix has not been factored. This routine can be used * as a test to see if solutions are correct. It should not be used * before spMNA_Preorder(). * * >>> Arguments: * eMatrix (char *) * Pointer to the matrix. * RHS (RealVector) * RHS is the right hand side. This is what is being solved for. * Solution (RealVector) * Solution is the vector being multiplied by the matrix. * iRHS (RealVector) * iRHS is the imaginary portion of the right hand side. This is * what is being solved for. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * iSolution (RealVector) * iSolution is the imaginary portion of the vector being multiplied * by the matrix. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * * >>> Obscure Macros * IMAG_VECTORS * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears * without a trace. */ void spMultiply( eMatrix, RHS, Solution IMAG_VECTORS ) char *eMatrix; RealVector RHS, Solution IMAG_VECTORS; { register ElementPtr pElement; register RealVector Vector; register RealNumber Sum; register int I, *pExtOrder; MatrixPtr Matrix = (MatrixPtr)eMatrix; extern void ComplexMatrixMultiply(); /* Begin `spMultiply'. */ ASSERT( IS_SPARSE( Matrix ) AND NOT Matrix->Factored ); if (NOT Matrix->RowsLinked) spcLinkRows(Matrix); #if spCOMPLEX if (Matrix->Complex) { ComplexMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS ); return; } #endif #if REAL #if NOT ARRAY_OFFSET /* Correct array pointers for ARRAY_OFFSET. */ --RHS; --Solution; #endif /* Initialize Intermediate vector with reordered Solution vector. */ Vector = Matrix->Intermediate; pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) Vector[I] = Solution[*(pExtOrder--)]; pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInRow[I]; Sum = 0.0; while (pElement != NULL) { Sum += pElement->Real * Vector[pElement->Col]; pElement = pElement->NextInRow; } RHS[*pExtOrder--] = Sum; } return; #endif /* REAL */ } #endif /* MULTIPLICATION */ #if spCOMPLEX AND MULTIPLICATION /* * COMPLEX MATRIX MULTIPLICATION * * Multiplies matrix by solution vector to find source vector. * Assumes matrix has not been factored. This routine can be used * as a test to see if solutions are correct. * * >>> Arguments: * Matrix (char *) * Pointer to the matrix. * RHS (RealVector) * RHS is the right hand side. This is what is being solved for. * This is only the real portion of the right-hand side if the matrix * is complex and spSEPARATED_COMPLEX_VECTORS is set true. * Solution (RealVector) * Solution is the vector being multiplied by the matrix. This is only * the real portion if the matrix is complex and * spSEPARATED_COMPLEX_VECTORS is set true. * iRHS (RealVector) * iRHS is the imaginary portion of the right hand side. This is * what is being solved for. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * iSolution (RealVector) * iSolution is the imaginary portion of the vector being multiplied * by the matrix. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * * >>> Obscure Macros * IMAG_VECTORS * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears * without a trace. */ static void ComplexMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS ) MatrixPtr Matrix; RealVector RHS, Solution IMAG_VECTORS; { register ElementPtr pElement; register ComplexVector Vector; register ComplexNumber Sum; register int I, *pExtOrder; /* Begin `ComplexMatrixMultiply'. */ /* Correct array pointers for ARRAY_OFFSET. */ #if NOT ARRAY_OFFSET #if spSEPARATED_COMPLEX_VECTORS --RHS; --iRHS; --Solution; --iSolution; #else RHS -= 2; Solution -= 2; #endif #endif /* Initialize Intermediate vector with reordered Solution vector. */ Vector = (ComplexVector)Matrix->Intermediate; pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; #if spSEPARATED_COMPLEX_VECTORS for (I = Matrix->Size; I > 0; I--) { Vector[I].Real = Solution[*pExtOrder]; Vector[I].Imag = iSolution[*(pExtOrder--)]; } #else for (I = Matrix->Size; I > 0; I--) Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)]; #endif pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInRow[I]; Sum.Real = Sum.Imag = 0.0; while (pElement != NULL) { /* Cmplx expression : Sum += Element * Vector[Col] */ CMPLX_MULT_ADD_ASSIGN( Sum, *pElement, Vector[pElement->Col] ); pElement = pElement->NextInRow; } #if spSEPARATED_COMPLEX_VECTORS RHS[*pExtOrder] = Sum.Real; iRHS[*pExtOrder--] = Sum.Imag; #else ((ComplexVector)RHS)[*pExtOrder--] = Sum; #endif } return; } #endif /* spCOMPLEX AND MULTIPLICATION */ #if MULTIPLICATION AND TRANSPOSE /* * TRANSPOSED MATRIX MULTIPLICATION * * Multiplies transposed matrix by solution vector to find source vector. * Assumes matrix has not been factored. This routine can be used * as a test to see if solutions are correct. It should not be used * before spMNA_Preorder(). * * >>> Arguments: * eMatrix (char *) * Pointer to the matrix. * RHS (RealVector) * RHS is the right hand side. This is what is being solved for. * Solution (RealVector) * Solution is the vector being multiplied by the matrix. * iRHS (RealVector) * iRHS is the imaginary portion of the right hand side. This is * what is being solved for. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * iSolution (RealVector) * iSolution is the imaginary portion of the vector being multiplied * by the matrix. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * * >>> Obscure Macros * IMAG_VECTORS * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears * without a trace. */ void spMultTransposed( eMatrix, RHS, Solution IMAG_VECTORS ) char *eMatrix; RealVector RHS, Solution IMAG_VECTORS; { register ElementPtr pElement; register RealVector Vector; register RealNumber Sum; register int I, *pExtOrder; MatrixPtr Matrix = (MatrixPtr)eMatrix; extern void ComplexTransposedMatrixMultiply(); /* Begin `spMultTransposed'. */ ASSERT( IS_SPARSE( Matrix ) AND NOT Matrix->Factored ); #if spCOMPLEX if (Matrix->Complex) { ComplexTransposedMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS ); return; } #endif #if REAL #if NOT ARRAY_OFFSET /* Correct array pointers for ARRAY_OFFSET. */ --RHS; --Solution; #endif /* Initialize Intermediate vector with reordered Solution vector. */ Vector = Matrix->Intermediate; pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) Vector[I] = Solution[*(pExtOrder--)]; pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInCol[I]; Sum = 0.0; while (pElement != NULL) { Sum += pElement->Real * Vector[pElement->Row]; pElement = pElement->NextInCol; } RHS[*pExtOrder--] = Sum; } return; #endif /* REAL */ } #endif /* MULTIPLICATION AND TRANSPOSE */ #if spCOMPLEX AND MULTIPLICATION AND TRANSPOSE /* * COMPLEX TRANSPOSED MATRIX MULTIPLICATION * * Multiplies transposed matrix by solution vector to find source vector. * Assumes matrix has not been factored. This routine can be used * as a test to see if solutions are correct. * * >>> Arguments: * Matrix (char *) * Pointer to the matrix. * RHS (RealVector) * RHS is the right hand side. This is what is being solved for. * This is only the real portion of the right-hand side if the matrix * is complex and spSEPARATED_COMPLEX_VECTORS is set true. * Solution (RealVector) * Solution is the vector being multiplied by the matrix. This is only * the real portion if the matrix is complex and * spSEPARATED_COMPLEX_VECTORS is set true. * iRHS (RealVector) * iRHS is the imaginary portion of the right hand side. This is * what is being solved for. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * iSolution (RealVector) * iSolution is the imaginary portion of the vector being multiplied * by the matrix. This is only necessary if the matrix is * complex and spSEPARATED_COMPLEX_VECTORS is true. * * >>> Obscure Macros * IMAG_VECTORS * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears * without a trace. */ static void ComplexTransposedMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS ) MatrixPtr Matrix; RealVector RHS, Solution IMAG_VECTORS; { register ElementPtr pElement; register ComplexVector Vector; register ComplexNumber Sum; register int I, *pExtOrder; /* Begin `ComplexMatrixMultiply'. */ /* Correct array pointers for ARRAY_OFFSET. */ #if NOT ARRAY_OFFSET #if spSEPARATED_COMPLEX_VECTORS --RHS; --iRHS; --Solution; --iSolution; #else RHS -= 2; Solution -= 2; #endif #endif /* Initialize Intermediate vector with reordered Solution vector. */ Vector = (ComplexVector)Matrix->Intermediate; pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size]; #if spSEPARATED_COMPLEX_VECTORS for (I = Matrix->Size; I > 0; I--) { Vector[I].Real = Solution[*pExtOrder]; Vector[I].Imag = iSolution[*(pExtOrder--)]; } #else for (I = Matrix->Size; I > 0; I--) Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)]; #endif pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInCol[I]; Sum.Real = Sum.Imag = 0.0; while (pElement != NULL) { /* Cmplx expression : Sum += Element * Vector[Row] */ CMPLX_MULT_ADD_ASSIGN( Sum, *pElement, Vector[pElement->Row] ); pElement = pElement->NextInCol; } #if spSEPARATED_COMPLEX_VECTORS RHS[*pExtOrder] = Sum.Real; iRHS[*pExtOrder--] = Sum.Imag; #else ((ComplexVector)RHS)[*pExtOrder--] = Sum; #endif } return; } #endif /* spCOMPLEX AND MULTIPLICATION AND TRANSPOSE */ #if DETERMINANT /* * CALCULATE DETERMINANT * * This routine in capable of calculating the determinant of the * matrix once the LU factorization has been performed. Hence, only * use this routine after spFactor() and before spClear(). * The determinant equals the product of all the diagonal elements of * the lower triangular matrix L, except that this product may need * negating. Whether the product or the negative product equals the * determinant is determined by the number of row and column * interchanges performed. Note that the determinants of matrices can * be very large or very small. On large matrices, the determinant * can be far larger or smaller than can be represented by a floating * point number. For this reason the determinant is scaled to a * reasonable value and the logarithm of the scale factor is returned. * * >>> Arguments: * eMatrix (char *) * A pointer to the matrix for which the determinant is desired. * pExponent (int *) * The logarithm base 10 of the scale factor for the determinant. To find * the actual determinant, Exponent should be added to the exponent of * Determinant. * pDeterminant (RealNumber *) * The real portion of the determinant. This number is scaled to be * greater than or equal to 1.0 and less than 10.0. * piDeterminant (RealNumber *) * The imaginary portion of the determinant. When the matrix is real * this pointer need not be supplied, nothing will be returned. This * number is scaled to be greater than or equal to 1.0 and less than 10.0. * * >>> Local variables: * Norm (RealNumber) * L-infinity norm of a complex number. * Size (int) * Local storage for Matrix->Size. Placed in a register for speed. * Temp (RealNumber) * Temporary storage for real portion of determinant. */ #if spCOMPLEX void spDeterminant( eMatrix, pExponent, pDeterminant, piDeterminant ) RealNumber *piDeterminant; #else void spDeterminant( eMatrix, pExponent, pDeterminant ) #endif char *eMatrix; register RealNumber *pDeterminant; int *pExponent; { register MatrixPtr Matrix = (MatrixPtr)eMatrix; register int I, Size; RealNumber Norm, nr, ni; ComplexNumber Pivot, cDeterminant; #define NORM(a) (nr = ABS((a).Real), ni = ABS((a).Imag), MAX (nr,ni)) /* Begin `spDeterminant'. */ ASSERT( IS_SPARSE( Matrix ) AND IS_FACTORED(Matrix) ); *pExponent = 0; if (Matrix->Error == spSINGULAR) { *pDeterminant = 0.0; #if spCOMPLEX if (Matrix->Complex) *piDeterminant = 0.0; #endif return; } Size = Matrix->Size; I = 0; #if spCOMPLEX if (Matrix->Complex) /* Complex Case. */ { cDeterminant.Real = 1.0; cDeterminant.Imag = 0.0; while (++I <= Size) { CMPLX_RECIPROCAL( Pivot, *Matrix->Diag[I] ); CMPLX_MULT_ASSIGN( cDeterminant, Pivot ); /* Scale Determinant. */ Norm = NORM( cDeterminant ); if (Norm != 0.0) { while (Norm >= 1.0e12) { cDeterminant.Real *= 1.0e-12; cDeterminant.Imag *= 1.0e-12; *pExponent += 12; Norm = NORM( cDeterminant ); } while (Norm < 1.0e-12) { cDeterminant.Real *= 1.0e12; cDeterminant.Imag *= 1.0e12; *pExponent -= 12; Norm = NORM( cDeterminant ); } } } /* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */ Norm = NORM( cDeterminant ); if (Norm != 0.0) { while (Norm >= 10.0) { cDeterminant.Real *= 0.1; cDeterminant.Imag *= 0.1; (*pExponent)++; Norm = NORM( cDeterminant ); } while (Norm < 1.0) { cDeterminant.Real *= 10.0; cDeterminant.Imag *= 10.0; (*pExponent)--; Norm = NORM( cDeterminant ); } } if (Matrix->NumberOfInterchangesIsOdd) CMPLX_NEGATE( cDeterminant ); *pDeterminant = cDeterminant.Real; *piDeterminant = cDeterminant.Imag; } #endif /* spCOMPLEX */ #if REAL AND spCOMPLEX else #endif #if REAL { /* Real Case. */ *pDeterminant = 1.0; while (++I <= Size) { *pDeterminant /= Matrix->Diag[I]->Real; /* Scale Determinant. */ if (*pDeterminant != 0.0) { while (ABS(*pDeterminant) >= 1.0e12) { *pDeterminant *= 1.0e-12; *pExponent += 12; } while (ABS(*pDeterminant) < 1.0e-12) { *pDeterminant *= 1.0e12; *pExponent -= 12; } } } /* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */ if (*pDeterminant != 0.0) { while (ABS(*pDeterminant) >= 10.0) { *pDeterminant *= 0.1; (*pExponent)++; } while (ABS(*pDeterminant) < 1.0) { *pDeterminant *= 10.0; (*pExponent)--; } } if (Matrix->NumberOfInterchangesIsOdd) *pDeterminant = -*pDeterminant; } #endif /* REAL */ } #endif /* DETERMINANT */ #if STRIP /* * STRIP FILL-INS FROM MATRIX * * Strips the matrix of all fill-ins. * * >>> Arguments: * Matrix (char *) * Pointer to the matrix to be stripped. * * >>> Local variables: * pElement (ElementPtr) * Pointer that is used to step through the matrix. * ppElement (ElementPtr *) * Pointer to the location of an ElementPtr. This location will be * updated if a fill-in is stripped from the matrix. * pFillin (ElementPtr) * Pointer used to step through the lists of fill-ins while marking them. * pLastFillin (ElementPtr) * A pointer to the last fill-in in the list. Used to terminate a loop. * pListNode (struct FillinListNodeStruct *) * A pointer to a node in the FillinList linked-list. */ void spStripFills( eMatrix ) char *eMatrix; { MatrixPtr Matrix = (MatrixPtr)eMatrix; struct FillinListNodeStruct *pListNode; /* Begin `spStripFills'. */ ASSERT( IS_SPARSE( Matrix ) ); if (Matrix->Fillins == 0) return; Matrix->NeedsOrdering = YES; Matrix->Elements -= Matrix->Fillins; Matrix->Fillins = 0; /* Mark the fill-ins. */ { register ElementPtr pFillin, pLastFillin; pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode; Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList; Matrix->NextAvailFillin = pListNode->pFillinList; while (pListNode != NULL) { pFillin = pListNode->pFillinList; pLastFillin = &(pFillin[ pListNode->NumberOfFillinsInList - 1 ]); while (pFillin <= pLastFillin) (pFillin++)->Row = 0; pListNode = pListNode->Next; } } /* Unlink fill-ins by searching for elements marked with Row = 0. */ { register ElementPtr pElement, *ppElement; register int I, Size = Matrix->Size; /* Unlink fill-ins in all columns. */ for (I = 1; I <= Size; I++) { ppElement = &(Matrix->FirstInCol[I]); while ((pElement = *ppElement) != NULL) { if (pElement->Row == 0) { *ppElement = pElement->NextInCol; /* Unlink fill-in. */ if (Matrix->Diag[pElement->Col] == pElement) Matrix->Diag[pElement->Col] = NULL; } else ppElement = &pElement->NextInCol; /* Skip element. */ } } /* Unlink fill-ins in all rows. */ for (I = 1; I <= Size; I++) { ppElement = &(Matrix->FirstInRow[I]); while ((pElement = *ppElement) != NULL) { if (pElement->Row == 0) *ppElement = pElement->NextInRow; /* Unlink fill-in. */ else ppElement = &pElement->NextInRow; /* Skip element. */ } } } return; } #endif #if TRANSLATE AND DELETE /* * DELETE A ROW AND COLUMN FROM THE MATRIX * * Deletes a row and a column from a matrix. * * Sparse will abort if an attempt is made to delete a row or column that * doesn't exist. * * >>> Arguments: * eMatrix (char *) * Pointer to the matrix in which the row and column are to be deleted. * Row (int) * Row to be deleted. * Col (int) * Column to be deleted. * * >>> Local variables: * ExtCol (int) * The external column that is being deleted. * ExtRow (int) * The external row that is being deleted. * pElement (ElementPtr) * Pointer to an element in the matrix. Used when scanning rows and * columns in order to eliminate elements from the last row or column. * ppElement (ElementPtr *) * Pointer to the location of an ElementPtr. This location will be * filled with a NULL pointer if it is the new last element in its row * or column. * pElement (ElementPtr) * Pointer to an element in the last row or column of the matrix. * Size (int) * The local version Matrix->Size, the size of the matrix. */ void spDeleteRowAndCol( eMatrix, Row, Col ) char *eMatrix; int Row, Col; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register ElementPtr pElement, *ppElement, pLastElement; int Size, ExtRow, ExtCol; ElementPtr spcFindElementInCol(); /* Begin `spDeleteRowAndCol'. */ ASSERT( IS_SPARSE(Matrix) AND Row > 0 AND Col > 0 ); ASSERT( Row <= Matrix->ExtSize AND Col <= Matrix->ExtSize ); Size = Matrix->Size; ExtRow = Row; ExtCol = Col; if (NOT Matrix->RowsLinked) spcLinkRows( Matrix ); Row = Matrix->ExtToIntRowMap[Row]; Col = Matrix->ExtToIntColMap[Col]; ASSERT( Row > 0 AND Col > 0 ); /* Move Row so that it is the last row in the matrix. */ if (Row != Size) spcRowExchange( Matrix, Row, Size ); /* Move Col so that it is the last column in the matrix. */ if (Col != Size) spcColExchange( Matrix, Col, Size ); /* Correct Diag pointers. */ if (Row == Col) SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Size] ) else { Matrix->Diag[Row] = spcFindElementInCol( Matrix, Matrix->FirstInCol+Row, Row, Row, NO ); Matrix->Diag[Col] = spcFindElementInCol( Matrix, Matrix->FirstInCol+Col, Col, Col, NO ); } /* * Delete last row and column of the matrix. */ /* Break the column links to every element in the last row. */ pLastElement = Matrix->FirstInRow[ Size ]; while (pLastElement != NULL) { ppElement = &(Matrix->FirstInCol[ pLastElement->Col ]); while ((pElement = *ppElement) != NULL) { if (pElement == pLastElement) *ppElement = NULL; /* Unlink last element in column. */ else ppElement = &pElement->NextInCol; /* Skip element. */ } pLastElement = pLastElement->NextInRow; } /* Break the row links to every element in the last column. */ pLastElement = Matrix->FirstInCol[ Size ]; while (pLastElement != NULL) { ppElement = &(Matrix->FirstInRow[ pLastElement->Row ]); while ((pElement = *ppElement) != NULL) { if (pElement == pLastElement) *ppElement = NULL; /* Unlink last element in row. */ else ppElement = &pElement->NextInRow; /* Skip element. */ } pLastElement = pLastElement->NextInCol; } /* Clean up some details. */ Matrix->Size = Size - 1; Matrix->Diag[Size] = NULL; Matrix->FirstInRow[Size] = NULL; Matrix->FirstInCol[Size] = NULL; Matrix->CurrentSize--; Matrix->ExtToIntRowMap[ExtRow] = -1; Matrix->ExtToIntColMap[ExtCol] = -1; Matrix->NeedsOrdering = YES; return; } #endif #if PSEUDOCONDITION /* * CALCULATE PSEUDOCONDITION * * Computes the magnitude of the ratio of the largest to the smallest * pivots. This quantity is an indicator of ill-conditioning in the * matrix. If this ratio is large, and if the matrix is scaled such * that uncertainties in the RHS and the matrix entries are * equilibrated, then the matrix is ill-conditioned. However, a small * ratio does not necessarily imply that the matrix is * well-conditioned. This routine must only be used after a matrix has * been factored by spOrderAndFactor() or spFactor() and before it is * cleared by spClear() or spInitialize(). The pseudocondition is * faster to compute than the condition number calculated by * spCondition(), but is not as informative. * * >>> Returns: * The magnitude of the ratio of the largest to smallest pivot used during * previous factorization. If the matrix was singular, zero is returned. * * >>> Arguments: * eMatrix (char *) * Pointer to the matrix. */ RealNumber spPseudoCondition( eMatrix ) char *eMatrix; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register int I; register ArrayOfElementPtrs Diag; RealNumber MaxPivot, MinPivot, Mag; /* Begin `spPseudoCondition'. */ ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) ); if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG) return 0.0; Diag = Matrix->Diag; MaxPivot = MinPivot = ELEMENT_MAG( Diag[1] ); for (I = 2; I <= Matrix->Size; I++) { Mag = ELEMENT_MAG( Diag[I] ); if (Mag > MaxPivot) MaxPivot = Mag; else if (Mag < MinPivot) MinPivot = Mag; } ASSERT( MaxPivot > 0.0); return MaxPivot / MinPivot; } #endif #if CONDITION /* * ESTIMATE CONDITION NUMBER * * Computes an estimate of the condition number using a variation on * the LINPACK condition number estimation algorithm. This quantity is * an indicator of ill-conditioning in the matrix. To avoid problems * with overflow, the reciprocal of the condition number is returned. * If this number is small, and if the matrix is scaled such that * uncertainties in the RHS and the matrix entries are equilibrated, * then the matrix is ill-conditioned. If the this number is near * one, the matrix is well conditioned. This routine must only be * used after a matrix has been factored by spOrderAndFactor() or * spFactor() and before it is cleared by spClear() or spInitialize(). * * Unlike the LINPACK condition number estimator, this routines * returns the L infinity condition number. This is an artifact of * Sparse placing ones on the diagonal of the upper triangular matrix * rather than the lower. This difference should be of no importance. * * References: * A.K. Cline, C.B. Moler, G.W. Stewart, J.H. Wilkinson. An estimate * for the condition number of a matrix. SIAM Journal on Numerical * Analysis. Vol. 16, No. 2, pages 368-375, April 1979. * * J.J. Dongarra, C.B. Moler, J.R. Bunch, G.W. Stewart. LINPACK * User's Guide. SIAM, 1979. * * Roger G. Grimes, John G. Lewis. Condition number estimation for * sparse matrices. SIAM Journal on Scientific and Statistical * Computing. Vol. 2, No. 4, pages 384-388, December 1981. * * Dianne Prost O'Leary. Estimating matrix condition numbers. SIAM * Journal on Scientific and Statistical Computing. Vol. 1, No. 2, * pages 205-209, June 1980. * * >>> Returns: * The reciprocal of the condition number. If the matrix was singular, * zero is returned. * * >>> Arguments: * eMatrix (char *) * Pointer to the matrix. * NormOfMatrix (RealNumber) * The L-infinity norm of the unfactored matrix as computed by * spNorm(). * pError (int *) * Used to return error code. * * >>> Possible errors: * spSINGULAR * spNO_MEMORY */ RealNumber spCondition( eMatrix, NormOfMatrix, pError ) char *eMatrix; RealNumber NormOfMatrix; int *pError; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register ElementPtr pElement; register RealVector T, Tm; register int I, K, Row; ElementPtr pPivot; int Size; RealNumber E, Em, Wp, Wm, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor; RealNumber Linpack, OLeary, InvNormOfInverse, ComplexCondition(); #define SLACK 1e4 /* Begin `spCondition'. */ ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) ); *pError = Matrix->Error; if (Matrix->Error >= spFATAL) return 0.0; if (NormOfMatrix == 0.0) { *pError = spSINGULAR; return 0.0; } #if spCOMPLEX if (Matrix->Complex) return ComplexCondition( Matrix, NormOfMatrix, pError ); #endif #if REAL Size = Matrix->Size; T = Matrix->Intermediate; #if spCOMPLEX Tm = Matrix->Intermediate + Size; #else Tm = ALLOC( RealNumber, Size+1 ); if (Tm == NULL) { *pError = spNO_MEMORY; return 0.0; } #endif for (I = Size; I > 0; I--) T[I] = 0.0; /* * Part 1. Ay = e. * Solve Ay = LUy = e where e consists of +1 and -1 terms with the sign * chosen to maximize the size of w in Lw = e. Since the terms in w can * get very large, scaling is used to avoid overflow. */ /* Forward elimination. Solves Lw = e while choosing e. */ E = 1.0; for (I = 1; I <= Size; I++) { pPivot = Matrix->Diag[I]; if (T[I] < 0.0) Em = -E; else Em = E; Wm = (Em + T[I]) * pPivot->Real; if (ABS(Wm) > SLACK) { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(Wm) ); for (K = Size; K > 0; K--) T[K] *= ScaleFactor; E *= ScaleFactor; Em *= ScaleFactor; Wm = (Em + T[I]) * pPivot->Real; } Wp = (T[I] - Em) * pPivot->Real; ASp = ABS(T[I] - Em); ASm = ABS(Em + T[I]); /* Update T for both values of W, minus value is placed in Tm. */ pElement = pPivot->NextInCol; while (pElement != NULL) { Row = pElement->Row; Tm[Row] = T[Row] - (Wm * pElement->Real); T[Row] -= (Wp * pElement->Real); ASp += ABS(T[Row]); ASm += ABS(Tm[Row]); pElement = pElement->NextInCol; } /* If minus value causes more growth, overwrite T with its values. */ if (ASm > ASp) { T[I] = Wm; pElement = pPivot->NextInCol; while (pElement != NULL) { T[pElement->Row] = Tm[pElement->Row]; pElement = pElement->NextInCol; } } else T[I] = Wp; } /* Compute 1-norm of T, which now contains w, and scale ||T|| to 1/SLACK. */ for (ASw = 0.0, I = Size; I > 0; I--) ASw += ABS(T[I]); ScaleFactor = 1.0 / (SLACK * ASw); if (ScaleFactor < 0.5) { for (I = Size; I > 0; I--) T[I] *= ScaleFactor; E *= ScaleFactor; } /* Backward Substitution. Solves Uy = w.*/ for (I = Size; I >= 1; I--) { pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) { T[I] -= pElement->Real * T[pElement->Col]; pElement = pElement->NextInRow; } if (ABS(T[I]) > SLACK) { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) ); for (K = Size; K > 0; K--) T[K] *= ScaleFactor; E *= ScaleFactor; } } /* Compute 1-norm of T, which now contains y, and scale ||T|| to 1/SLACK. */ for (ASy = 0.0, I = Size; I > 0; I--) ASy += ABS(T[I]); ScaleFactor = 1.0 / (SLACK * ASy); if (ScaleFactor < 0.5) { for (I = Size; I > 0; I--) T[I] *= ScaleFactor; ASy = 1.0 / SLACK; E *= ScaleFactor; } /* Compute infinity-norm of T for O'Leary's estimate. */ for (MaxY = 0.0, I = Size; I > 0; I--) if (MaxY < ABS(T[I])) MaxY = ABS(T[I]); /* * Part 2. A* z = y where the * represents the transpose. * Recall that A = LU implies A* = U* L*. */ /* Forward elimination, U* v = y. */ for (I = 1; I <= Size; I++) { pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) { T[pElement->Col] -= T[I] * pElement->Real; pElement = pElement->NextInRow; } if (ABS(T[I]) > SLACK) { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) ); for (K = Size; K > 0; K--) T[K] *= ScaleFactor; ASy *= ScaleFactor; } } /* Compute 1-norm of T, which now contains v, and scale ||T|| to 1/SLACK. */ for (ASv = 0.0, I = Size; I > 0; I--) ASv += ABS(T[I]); ScaleFactor = 1.0 / (SLACK * ASv); if (ScaleFactor < 0.5) { for (I = Size; I > 0; I--) T[I] *= ScaleFactor; ASy *= ScaleFactor; } /* Backward Substitution, L* z = v. */ for (I = Size; I >= 1; I--) { pPivot = Matrix->Diag[I]; pElement = pPivot->NextInCol; while (pElement != NULL) { T[I] -= pElement->Real * T[pElement->Row]; pElement = pElement->NextInCol; } T[I] *= pPivot->Real; if (ABS(T[I]) > SLACK) { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) ); for (K = Size; K > 0; K--) T[K] *= ScaleFactor; ASy *= ScaleFactor; } } /* Compute 1-norm of T, which now contains z. */ for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABS(T[I]); #if NOT spCOMPLEX FREE( Tm ); #endif Linpack = ASy / ASz; OLeary = E / MaxY; InvNormOfInverse = MIN( Linpack, OLeary ); return InvNormOfInverse / NormOfMatrix; #endif /* REAL */ } #if spCOMPLEX /* * ESTIMATE CONDITION NUMBER * * Complex version of spCondition(). * * >>> Returns: * The reciprocal of the condition number. * * >>> Arguments: * Matrix (MatrixPtr) * Pointer to the matrix. * NormOfMatrix (RealNumber) * The L-infinity norm of the unfactored matrix as computed by * spNorm(). * pError (int *) * Used to return error code. * * >>> Possible errors: * spNO_MEMORY */ static RealNumber ComplexCondition( Matrix, NormOfMatrix, pError ) MatrixPtr Matrix; RealNumber NormOfMatrix; int *pError; { register ElementPtr pElement; register ComplexVector T, Tm; register int I, K, Row; ElementPtr pPivot; int Size; RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor; RealNumber Linpack, OLeary, InvNormOfInverse; ComplexNumber Wp, Wm; /* Begin `ComplexCondition'. */ Size = Matrix->Size; T = (ComplexVector)Matrix->Intermediate; Tm = ALLOC( ComplexNumber, Size+1 ); if (Tm == NULL) { *pError = spNO_MEMORY; return 0.0; } for (I = Size; I > 0; I--) T[I].Real = T[I].Imag = 0.0; /* * Part 1. Ay = e. * Solve Ay = LUy = e where e consists of +1 and -1 terms with the sign * chosen to maximize the size of w in Lw = e. Since the terms in w can * get very large, scaling is used to avoid overflow. */ /* Forward elimination. Solves Lw = e while choosing e. */ E = 1.0; for (I = 1; I <= Size; I++) { pPivot = Matrix->Diag[I]; if (T[I].Real < 0.0) Em = -E; else Em = E; Wm = T[I]; Wm.Real += Em; ASm = CMPLX_1_NORM( Wm ); CMPLX_MULT_ASSIGN( Wm, *pPivot ); if (CMPLX_1_NORM(Wm) > SLACK) { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(Wm) ); for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor ); E *= ScaleFactor; Em *= ScaleFactor; ASm *= ScaleFactor; SCLR_MULT_ASSIGN( Wm, ScaleFactor ); } Wp = T[I]; Wp.Real -= Em; ASp = CMPLX_1_NORM( Wp ); CMPLX_MULT_ASSIGN( Wp, *pPivot ); /* Update T for both values of W, minus value is placed in Tm. */ pElement = pPivot->NextInCol; while (pElement != NULL) { Row = pElement->Row; /* Cmplx expr: Tm[Row] = T[Row] - (Wp * *pElement). */ CMPLX_MULT_SUBT( Tm[Row], Wm, *pElement, T[Row] ); /* Cmplx expr: T[Row] -= Wp * *pElement. */ CMPLX_MULT_SUBT_ASSIGN( T[Row], Wm, *pElement ); ASp += CMPLX_1_NORM(T[Row]); ASm += CMPLX_1_NORM(Tm[Row]); pElement = pElement->NextInCol; } /* If minus value causes more growth, overwrite T with its values. */ if (ASm > ASp) { T[I] = Wm; pElement = pPivot->NextInCol; while (pElement != NULL) { T[pElement->Row] = Tm[pElement->Row]; pElement = pElement->NextInCol; } } else T[I] = Wp; } /* Compute 1-norm of T, which now contains w, and scale ||T|| to 1/SLACK. */ for (ASw = 0.0, I = Size; I > 0; I--) ASw += CMPLX_1_NORM(T[I]); ScaleFactor = 1.0 / (SLACK * ASw); if (ScaleFactor < 0.5) { for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor ); E *= ScaleFactor; } /* Backward Substitution. Solves Uy = w.*/ for (I = Size; I >= 1; I--) { pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) { /* Cmplx expr: T[I] -= T[pElement->Col] * *pElement. */ CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Col], *pElement ); pElement = pElement->NextInRow; } if (CMPLX_1_NORM(T[I]) > SLACK) { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) ); for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor ); E *= ScaleFactor; } } /* Compute 1-norm of T, which now contains y, and scale ||T|| to 1/SLACK. */ for (ASy = 0.0, I = Size; I > 0; I--) ASy += CMPLX_1_NORM(T[I]); ScaleFactor = 1.0 / (SLACK * ASy); if (ScaleFactor < 0.5) { for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor ); ASy = 1.0 / SLACK; E *= ScaleFactor; } /* Compute infinity-norm of T for O'Leary's estimate. */ for (MaxY = 0.0, I = Size; I > 0; I--) if (MaxY < CMPLX_1_NORM(T[I])) MaxY = CMPLX_1_NORM(T[I]); /* * Part 2. A* z = y where the * represents the transpose. * Recall that A = LU implies A* = U* L*. */ /* Forward elimination, U* v = y. */ for (I = 1; I <= Size; I++) { pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) { /* Cmplx expr: T[pElement->Col] -= T[I] * *pElement. */ CMPLX_MULT_SUBT_ASSIGN( T[pElement->Col], T[I], *pElement ); pElement = pElement->NextInRow; } if (CMPLX_1_NORM(T[I]) > SLACK) { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) ); for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor ); ASy *= ScaleFactor; } } /* Compute 1-norm of T, which now contains v, and scale ||T|| to 1/SLACK. */ for (ASv = 0.0, I = Size; I > 0; I--) ASv += CMPLX_1_NORM(T[I]); ScaleFactor = 1.0 / (SLACK * ASv); if (ScaleFactor < 0.5) { for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor ); ASy *= ScaleFactor; } /* Backward Substitution, L* z = v. */ for (I = Size; I >= 1; I--) { pPivot = Matrix->Diag[I]; pElement = pPivot->NextInCol; while (pElement != NULL) { /* Cmplx expr: T[I] -= T[pElement->Row] * *pElement. */ CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Row], *pElement ); pElement = pElement->NextInCol; } CMPLX_MULT_ASSIGN( T[I], *pPivot ); if (CMPLX_1_NORM(T[I]) > SLACK) { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) ); for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor ); ASy *= ScaleFactor; } } /* Compute 1-norm of T, which now contains z. */ for (ASz = 0.0, I = Size; I > 0; I--) ASz += CMPLX_1_NORM(T[I]); FREE( Tm ); Linpack = ASy / ASz; OLeary = E / MaxY; InvNormOfInverse = MIN( Linpack, OLeary ); return InvNormOfInverse / NormOfMatrix; } #endif /* spCOMPLEX */ /* * L-INFINITY MATRIX NORM * * Computes the L-infinity norm of an unfactored matrix. It is a fatal * error to pass this routine a factored matrix. * * One difficulty is that the rows may not be linked. * * >>> Returns: * The largest absolute row sum of matrix. * * >>> Arguments: * eMatrix (char *) * Pointer to the matrix. */ RealNumber spNorm( eMatrix ) char *eMatrix; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register ElementPtr pElement; register int I; RealNumber Max = 0.0, AbsRowSum; /* Begin `spNorm'. */ ASSERT( IS_SPARSE(Matrix) AND NOT IS_FACTORED(Matrix) ); if (NOT Matrix->RowsLinked) spcLinkRows( Matrix ); /* Compute row sums. */ #if REAL if (NOT Matrix->Complex) { for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInRow[I]; AbsRowSum = 0.0; while (pElement != NULL) { AbsRowSum += ABS( pElement->Real ); pElement = pElement->NextInRow; } if (Max < AbsRowSum) Max = AbsRowSum; } } #endif #if spCOMPLEX if (Matrix->Complex) { for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInRow[I]; AbsRowSum = 0.0; while (pElement != NULL) { AbsRowSum += CMPLX_1_NORM( *pElement ); pElement = pElement->NextInRow; } if (Max < AbsRowSum) Max = AbsRowSum; } } #endif return Max; } #endif /* CONDITION */ #if STABILITY /* * STABILITY OF FACTORIZATION * * The following routines are used to gauge the stability of a * factorization. If the factorization is determined to be too unstable, * then the matrix should be reordered. The routines compute quantities * that are needed in the computation of a bound on the error attributed * to any one element in the matrix during the factorization. In other * words, there is a matrix E = [e_ij] of error terms such that A+E = LU. * This routine finds a bound on |e_ij|. Erisman & Reid [1] showed that * |e_ij| < 3.01 u rho m_ij, where u is the machine rounding unit, * rho = max a_ij where the max is taken over every row i, column j, and * step k, and m_ij is the number of multiplications required in the * computation of l_ij if i > j or u_ij otherwise. Barlow [2] showed that * rho < max_i || l_i ||_p max_j || u_j ||_q where 1/p + 1/q = 1. * * The first routine finds the magnitude on the largest element in the * matrix. If the matrix has not yet been factored, the largest * element is found by direct search. If the matrix is factored, a * bound on the largest element in any of the reduced submatrices is * computed using Barlow with p = oo and q = 1. The ratio of these * two numbers is the growth, which can be used to determine if the * pivoting order is adequate. A large growth implies that * considerable error has been made in the factorization and that it * is probably a good idea to reorder the matrix. If a large growth * in encountered after using spFactor(), reconstruct the matrix and * refactor using spOrderAndFactor(). If a large growth is * encountered after using spOrderAndFactor(), refactor using * spOrderAndFactor() with the pivot threshold increased, say to 0.1. * * Using only the size of the matrix as an upper bound on m_ij and * Barlow's bound, the user can estimate the size of the matrix error * terms e_ij using the bound of Erisman and Reid. The second routine * computes a tighter bound (with more work) based on work by Gear * [3], |e_ij| < 1.01 u rho (t c^3 + (1 + t)c^2) where t is the * threshold and c is the maximum number of off-diagonal elements in * any row of L. The expensive part of computing this bound is * determining the maximum number of off-diagonals in L, which changes * only when the order of the matrix changes. This number is computed * and saved, and only recomputed if the matrix is reordered. * * [1] A. M. Erisman, J. K. Reid. Monitoring the stability of the * triangular factorization of a sparse matrix. Numerische * Mathematik. Vol. 22, No. 3, 1974, pp 183-186. * * [2] J. L. Barlow. A note on monitoring the stability of triangular * decomposition of sparse matrices. "SIAM Journal of Scientific * and Statistical Computing." Vol. 7, No. 1, January 1986, pp 166-168. * * [3] I. S. Duff, A. M. Erisman, J. K. Reid. "Direct Methods for Sparse * Matrices." Oxford 1986. pp 99. */ /* * LARGEST ELEMENT IN MATRIX * * >>> Returns: * If matrix is not factored, returns the magnitude of the largest element in * the matrix. If the matrix is factored, a bound on the magnitude of the * largest element in any of the reduced submatrices is returned. * * >>> Arguments: * eMatrix (char *) * Pointer to the matrix. */ RealNumber spLargestElement( eMatrix ) char *eMatrix; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register int I; RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0; RealNumber Pivot; ComplexNumber cPivot; register ElementPtr pElement, pDiag; /* Begin `spLargestElement'. */ ASSERT( IS_SPARSE(Matrix) ); #if REAL if (Matrix->Factored AND NOT Matrix->Complex) { if (Matrix->Error == spSINGULAR) return 0.0; /* Find the bound on the size of the largest element over all factorization. */ for (I = 1; I <= Matrix->Size; I++) { pDiag = Matrix->Diag[I]; /* Lower triangular matrix. */ Pivot = 1.0 / pDiag->Real; Mag = ABS( Pivot ); if (Mag > MaxRow) MaxRow = Mag; pElement = Matrix->FirstInRow[I]; while (pElement != pDiag) { Mag = ABS( pElement->Real ); if (Mag > MaxRow) MaxRow = Mag; pElement = pElement->NextInRow; } /* Upper triangular matrix. */ pElement = Matrix->FirstInCol[I]; AbsColSum = 1.0; /* Diagonal of U is unity. */ while (pElement != pDiag) { AbsColSum += ABS( pElement->Real ); pElement = pElement->NextInCol; } if (AbsColSum > MaxCol) MaxCol = AbsColSum; } } else if (NOT Matrix->Complex) { for (I = 1; I <= Matrix->Size; I++) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { Mag = ABS( pElement->Real ); if (Mag > Max) Max = Mag; pElement = pElement->NextInCol; } } return Max; } #endif #if spCOMPLEX if (Matrix->Factored AND Matrix->Complex) { if (Matrix->Error == spSINGULAR) return 0.0; /* Find the bound on the size of the largest element over all factorization. */ for (I = 1; I <= Matrix->Size; I++) { pDiag = Matrix->Diag[I]; /* Lower triangular matrix. */ CMPLX_RECIPROCAL( cPivot, *pDiag ); Mag = CMPLX_1_NORM( cPivot ); if (Mag > MaxRow) MaxRow = Mag; pElement = Matrix->FirstInRow[I]; while (pElement != pDiag) { Mag = CMPLX_1_NORM( *pElement ); if (Mag > MaxRow) MaxRow = Mag; pElement = pElement->NextInRow; } /* Upper triangular matrix. */ pElement = Matrix->FirstInCol[I]; AbsColSum = 1.0; /* Diagonal of U is unity. */ while (pElement != pDiag) { AbsColSum += CMPLX_1_NORM( *pElement ); pElement = pElement->NextInCol; } if (AbsColSum > MaxCol) MaxCol = AbsColSum; } } else if (Matrix->Complex) { for (I = 1; I <= Matrix->Size; I++) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { Mag = CMPLX_1_NORM( *pElement ); if (Mag > Max) Max = Mag; pElement = pElement->NextInCol; } } return Max; } #endif return MaxRow * MaxCol; } /* * MATRIX ROUNDOFF ERROR * * >>> Returns: * Returns a bound on the magnitude of the largest element in E = A - LU. * * >>> Arguments: * eMatrix (char *) * Pointer to the matrix. * Rho (RealNumber) * The bound on the magnitude of the largest element in any of the * reduced submatrices. This is the number computed by the function * spLargestElement() when given a factored matrix. If this number is * negative, the bound will be computed automatically. */ RealNumber spRoundoff( eMatrix, Rho ) char *eMatrix; RealNumber Rho; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register ElementPtr pElement; register int Count, I, MaxCount = 0; RealNumber Reid, Gear; /* Begin `spRoundoff'. */ ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) ); /* Compute Barlow's bound if it is not given. */ if (Rho < 0.0) Rho = spLargestElement( eMatrix ); /* Find the maximum number of off-diagonals in L if not previously computed. */ if (Matrix->MaxRowCountInLowerTri < 0) { for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInRow[I]; Count = 0; while (pElement->Col < I) { Count++; pElement = pElement->NextInRow; } if (Count > MaxCount) MaxCount = Count; } Matrix->MaxRowCountInLowerTri = MaxCount; } else MaxCount = Matrix->MaxRowCountInLowerTri; /* Compute error bound. */ Gear = 1.01*((MaxCount + 1) * Matrix->RelThreshold + 1.0) * SQR(MaxCount); Reid = 3.01 * Matrix->Size; if (Gear < Reid) return (MACHINE_RESOLUTION * Rho * Gear); else return (MACHINE_RESOLUTION * Rho * Reid); } #endif @EOF chmod 644 spUtils.c echo x - mat0 cat >mat0 <<'@EOF' mat0 -- Simple matrix. 4 real 1 1 2.0 1 2 -1.0 2 1 -1.0 2 2 3.0 2 3 -1.0 3 2 -1.0 3 3 3.0 3 4 -1.0 4 3 -1.0 4 4 3.0 0 0 0.0 34.0 0.0 0.0 0.0 @EOF chmod 644 mat0 echo x - mat1 cat >mat1 <<'@EOF' mat1 -- Typical circuit matrix, somewhat ill-conditioned. 23 real 1 1 0.0201303 1 6 -0.0108325 1 19 -0.00992221 2 2 0.00196826 2 6 -7.00819e-05 2 16 -0.000149258 3 3 0.00233598 3 6 -6.89538e-05 3 13 -3.81701e-05 3 14 -0.000161645 4 4 0.00181074 4 6 -7.11072e-05 5 5 0.00193524 5 6 -8.35603e-05 5 8 -0.000160605 6 1 -1.72443e-05 6 2 -9.41688e-05 6 3 -4.87603e-05 6 4 -9.88527e-05 6 5 -1.01501e-05 6 6 0.0237867 7 7 0.0020023 7 20 -0.002 8 5 -6.08416e-05 8 8 0.000306666 8 11 3.52493e-05 9 9 0.00012008 10 10 0.00200237 10 22 -0.002 11 11 9.74511e-05 12 12 0.000182758 13 3 -0.000458012 13 13 0.00205898 13 21 -0.002 14 3 -4.65368e-05 14 14 0.000305252 14 17 3.7642e-05 15 15 0.00012008 16 2 -0.000156934 16 16 0.00217081 16 23 -0.002 17 17 9.5915e-05 18 18 0.000182758 19 1 -0.00413682 19 19 0.0189592 20 7 -0.002 20 20 10.0037 20 21 -10 21 13 -0.002 21 20 -10 21 21 10.002 22 10 -0.002 22 22 10.0037 22 23 -10 23 16 -0.002 23 22 -10 23 23 10.002 0 0 0.0 -2.57565e-05 -4.48769e-07 -4.13483e-06 -1.01714e-06 -1.37342e-05 2.90276e-07 1.51576e-09 -1.56322e-07 0 5.2198e-08 0 0 -4.63536e-07 1.99061e-06 0 1.95645e-06 0 0 3.74634e-05 6.29115e-06 -4.69183e-06 -3.57397e-06 1.88248e-06 @EOF chmod 644 mat1 echo x - mat2 cat >mat2 <<'@EOF' mat2 -- Typical circuit matrix. 8 real 1 1 0.0010686 1 2 -1.66e-05 1 5 -0.001 1 6 -4.2e-05 2 1 5.9554e-05 2 2 0.000131298 2 6 -9.16434e-05 2 7 0.000148333 3 3 0.00213033 3 5 -0.002 3 7 0.000190213 3 8 -4.2e-05 4 4 0.000275862 4 6 -4.6936e-05 5 1 -0.001 5 3 -0.002 5 5 0.00308858 5 8 0.000666678 6 1 -0.000118154 6 2 -6.68894e-08 6 4 -4.71518e-05 6 6 0.000366644 7 2 -0.000200677 7 3 -1.8e-06 7 7 0.000333767 8 3 -0.000198585 8 5 -7.2e-06 8 8 0.000431315 0 0 0.0 4.03934e-07 -3.044e-07 -1.53475e-05 5.28676e-08 1.5958e-05 -9.9608e-08 -2.13155e-07 4.15313e-09 @EOF chmod 644 mat2 echo x - mat3 cat >mat3 <<'@EOF' mat3 -- Ill-conditioned matrix. 22 real 1 1 0.0201303 1 6 -0.0108325 1 19 -0.00992221 2 2 0.00196826 2 6 -7.00819e-05 2 16 -0.000149258 3 3 0.00233598 3 6 -6.89538e-05 3 13 -3.81701e-05 3 14 -0.000161645 4 4 0.00181074 4 6 -7.11072e-05 5 5 0.00193524 5 6 -8.35603e-05 5 8 -0.000160605 6 1 -1.72443e-05 6 2 -9.41688e-05 6 3 -4.87603e-05 6 4 -9.88527e-05 6 5 -1.01501e-05 6 6 0.0237867 7 5 -6.08416e-05 7 8 0.000306666 7 11 3.52493e-05 8 9 0.00012008 9 10 0.00200237 9 22 -0.002 10 11 9.74511e-05 11 12 0.000182758 12 3 -0.000458012 12 13 0.00205898 12 21 -0.002 13 3 -4.65368e-05 13 14 0.000305252 13 17 3.7642e-05 14 15 0.00012008 15 2 -0.000156934 15 16 0.00217081 16 17 9.5915e-05 17 18 0.000182758 18 1 -0.00413682 18 19 0.0189592 19 7 -0.002 19 20 10.0037 19 21 -10 20 13 -0.002 20 20 -10 20 21 10.002 21 10 -0.002 21 22 10.0037 22 16 -0.002 22 22 -10 0 0 0.0 -2.57565e-05 -4.48769e-07 -4.13483e-06 -1.01714e-06 -1.37342e-05 2.90276e-07 1.51576e-09 -1.56322e-07 0 5.2198e-08 0 0 -4.63536e-07 1.99061e-06 0 1.95645e-06 0 0 3.74634e-05 6.29115e-06 -4.69183e-06 -3.57397e-06 @EOF chmod 644 mat3 echo x - mat5 cat >mat5 <<'@EOF' mat5 -- Test of complete pivoting capability. 3 real 1 2 1.0 2 3 1.0 3 1 1.0 0 0 0.0 2.0 3.0 1.0 @EOF chmod 644 mat5 exit 0 .