/*
* 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 */