luke
/
ngspice
mirror of https://git.code.sf.net/p/ngspice/ngspice
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

init sparse v1.4b

This commit is contained in:
dwarning 2023-06-05 14:07:56 +02:00
parent 10e86d72d3
commit 9f95c79537
15 changed files with 5971 additions and 5343 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

381
src/include/ngspice/spmatrix.h
View File

@ -1,9 +1,6 @@
/*
* EXPORTS for sparse matrix routines with SPICE3.
*
* Author: Advising professor:
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
* UC Berkeley
/* EXPORTS for sparse matrix routines. */
/*!
* \file
*
* This file contains definitions that are useful to the calling
* program. In particular, this file contains error keyword
@ -13,25 +10,21 @@
* Also included is the type definitions for the various functions
* available to the user.
*
* This file is a modified version of spMatrix.h that is used when
* interfacing to Spice3.
* Objects that begin with the \a spc prefix are considered private
* and should not be used.
*
* \author
* Kenneth S. Kundert <kundert@users.sourceforge.net>
*/
/*
* Revision and copyright information.
*
* Copyright (c) 1985,86,87,88,89,90
* by Kenneth S. Kundert and the University of California.
* Copyright (c) 1985-2003 by Kenneth S. Kundert
*
* 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.
* $Date: 2003/06/29 04:19:52 $
* $Revision: 1.2 $
*/
@ -39,6 +32,17 @@
#ifndef spOKAY
/*
* IMPORTS
*
* >>> Import descriptions:
* spConfig.h
* Macros that customize the sparse matrix routines.
*/
#include "../../maths/sparse/spConfig.h"
@ -50,47 +54,54 @@
* changed under the condition that the codes for the nonfatal errors are
* less than the code for spFATAL and similarly the codes for the fatal
* errors are greater than that for spFATAL.
*
* >>> Error descriptions:
* spOKAY
* No error has occurred.
* spSMALL_PIVOT
* When reordering the matrix, no element was found which satisfies the
* threshold criteria. The largest element in the matrix was chosen
* as pivot. Non-fatal.
* spZERO_DIAG
* Fatal error. A zero was encountered on the diagonal the matrix. This
* does not necessarily imply that the matrix is singular. When this
* error occurs, the matrix should be reconstructed and factored using
* spOrderAndFactor().
* spSINGULAR
* Fatal error. Matrix is singular, so no unique solution exists.
* spNO_MEMORY
* Fatal error. Indicates that not enough memory is available to handle
* the matrix.
* spPANIC
* Fatal error indicating that the routines are not prepared to
* handle the matrix that has been requested. This may occur when
* the matrix is specified to be real and the routines are not
* compiled for real matrices, or when the matrix is specified to
* be complex and the routines are not compiled to handle complex
* matrices.
* spFATAL
* Not an error flag, but rather the dividing line between fatal errors
* and warnings.
*/
#include "ngspice/sperror.h" /* Spice error definitions. */
/* Begin error macros. */
#define spOKAY OK
#define spSMALL_PIVOT OK
#define spZERO_DIAG E_SINGULAR
#define spSINGULAR E_SINGULAR
#define spNO_MEMORY E_NOMEM
#define spPANIC E_BADMATRIX
#define spFATAL E_BADMATRIX
#define spOKAY 0 /*!<
* Error code that indicates that no error has
* occurred.
*/
#define spSMALL_PIVOT 1 /*!<
* Non-fatal error code that indicates that, when
* reordering the matrix, no element was found that
* satisfies the absolute threshold criteria. The
* largest element in the matrix was chosen as pivot.
*/
#define spZERO_DIAG 2 /*!<
* Fatal error code that indicates that, a zero was
* encountered on the diagonal the matrix. This does
* not necessarily imply that the matrix is singular.
* When this error occurs, the matrix should be
* reconstructed and factored using
* spOrderAndFactor().
*/
#define spSINGULAR 3 /*!<
* Fatal error code that indicates that, matrix is
* singular, so no unique solution exists.
*/
#define spMANGLED 4 /*!<
* Fatal error code that indicates that, matrix has
* been mangled, results of requested operation are
* garbage.
*/
#define spNO_MEMORY 5 /*!<
* Fatal error code that indicates that not enough
* memory is available.
*/
#define spPANIC 6 /*!<
* Fatal error code that indicates that the routines
* are not prepared to handle the matrix that has
* been requested. This may occur when the matrix
* is specified to be real and the routines are not
* compiled for real matrices, or when the matrix is
* specified to be complex and the routines are not
* compiled to handle complex matrices.
*/
#define spFATAL 2 /*!<
* Error code that is not an error flag, but rather
* the dividing line between fatal errors and
* warnings.
*/
@ -99,21 +110,18 @@
/*
* KEYWORD DEFINITIONS
*
* Here we define what precision arithmetic Sparse will use. Double
* precision is suggested as being most appropriate for circuit
* simulation and for C. However, it is possible to change spREAL
* to a float for single precision arithmetic. Note that in C, single
* precision arithmetic is often slower than double precision. Sparse
* internally refers to spREALs as RealNumbers.
*
* Some C compilers, notably the old VMS compiler, do not handle the keyword
* "void" correctly. If this is true for your compiler, remove the
* comment delimiters from the redefinition of void to int below.
*/
#define spREAL double
/* #define void int */
#define spREAL double /*!<
* Defines the precision of the arithmetic used by
* \a Sparse will use. Double precision is suggested
* as being most appropriate for circuit simulation
* and for C. However, it is possible to change spREAL
* to a float for single precision arithmetic. Note
* that in C, single precision arithmetic is often
* slower than double precision. Sparse
* internally refers to spREALs as RealNumbers.
*/
@ -140,10 +148,32 @@
/* Begin partition keywords. */
#define spDEFAULT_PARTITION 0
#define spDIRECT_PARTITION 1
#define spINDIRECT_PARTITION 2
#define spAUTO_PARTITION 3
#define spDEFAULT_PARTITION 0 /*!<
* Partition code for spPartition().
* Indicates that the default partitioning
* mode should be used.
* \see spPartition()
*/
#define spDIRECT_PARTITION 1 /*!<
* Partition code for spPartition().
* Indicates that all rows should be placed
* in the direct addressing partition.
* \see spPartition()
*/
#define spINDIRECT_PARTITION 2 /*!<
* Partition code for spPartition().
* Indicates that all rows should be placed
* in the indirect addressing partition.
* \see spPartition()
*/
#define spAUTO_PARTITION 3 /*!<
* Partition code for spPartition().
* Indicates that \a Sparse should chose
* the best partition for each row based
* on some simple rules. This is generally
* preferred.
* \see spPartition()
*/
@ -151,38 +181,33 @@
/*
* MACRO FUNCTION DEFINITIONS
*
* >>> Macro descriptions:
* spADD_REAL_ELEMENT
* Macro function that adds data to a real element in the matrix by a
* pointer.
* spADD_IMAG_ELEMENT
* Macro function that adds data to a imaginary element in the matrix by
* a pointer.
* spADD_COMPLEX_ELEMENT
* Macro function that adds data to a complex element in the matrix by a
* pointer.
* spADD_REAL_QUAD
* Macro function that adds data to each of the four real matrix elements
* specified by the given template.
* spADD_IMAG_QUAD
* Macro function that adds data to each of the four imaginary matrix
* elements specified by the given template.
* spADD_COMPLEX_QUAD
* Macro function that adds data to each of the four complex matrix
* elements specified by the given template.
*/
/* Begin Macros. */
/*!
* Macro function that adds data to a real element in the matrix by a pointer.
*/
#define spADD_REAL_ELEMENT(element,real) *(element) += real
/*!
* Macro function that adds data to a imaginary element in the matrix by
* a pointer.
*/
#define spADD_IMAG_ELEMENT(element,imag) *(element+1) += imag
/*!
* Macro function that adds data to a complex element in the matrix by
* a pointer.
*/
#define spADD_COMPLEX_ELEMENT(element,real,imag) \
{ *(element) += real; \
*(element+1) += imag; \
}
/*!
* Macro function that adds data to each of the four real matrix elements
* specified by the given template.
*/
#define spADD_REAL_QUAD(template,real) \
{ *((template).Element1) += real; \
*((template).Element2) += real; \
@ -190,6 +215,10 @@
*((template).Element4Negated) -= real; \
}
/*!
* Macro function that adds data to each of the four imaginary matrix
* elements specified by the given template.
*/
#define spADD_IMAG_QUAD(template,imag) \
{ *((template).Element1+1) += imag; \
*((template).Element2+1) += imag; \
@ -197,6 +226,10 @@
*((template).Element4Negated+1) -= imag; \
}
/*!
* Macro function that adds data to each of the four complex matrix
* elements specified by the given template.
*/
#define spADD_COMPLEX_QUAD(template,real,imag) \
{ *((template).Element1) += real; \
*((template).Element2) += real; \
@ -212,33 +245,50 @@
/*
* TYPE DEFINITION FOR COMPONENT TEMPLATE
* TYPE DEFINITION FOR EXTERNAL MATRIX ELEMENT REFERENCES
*
* External type definitions for Sparse data objects.
*/
/*! Declares the type of the a pointer to a matrix. */
typedef spGenericPtr spMatrix;
/*! Declares the type of the a pointer to a matrix element. */
typedef spREAL spElement;
/*! Declares the type of the Sparse error codes. */
typedef int spError;
/* TYPE DEFINITION FOR COMPONENT TEMPLATE */
/*!
* This data structure is used to hold pointers to four related elements in
* matrix. It is used in conjunction with the routines
* spGetAdmittance
* spGetQuad
* spGetOnes
* These routines stuff the structure which is later used by the spADD_QUAD
* macro functions above. It is also possible for the user to collect four
* pointers returned by spGetElement and stuff them into the template.
* The spADD_QUAD routines stuff data into the matrix in locations specified
* by Element1 and Element2 without changing the data. The data is negated
* before being placed in Element3 and Element4.
* matrix. It is used in conjunction with the routines spGetAdmittance(),
* spGetQuad(), and spGetOnes(). These routines stuff the structure which
* is later used by the \a spADD_QUAD macro functions above. It is also
* possible for the user to collect four pointers returned by spGetElement()
* and stuff them into the template. The \a spADD_QUAD routines stuff data
* into the matrix in locations specified by \a Element1 and \a Element2
* without changing the data. The data is negated before being placed in
* \a Element3 and \a Element4.
*/
/* Begin `spTemplate'. */
struct spTemplate
{ spREAL *Element1 ;
spREAL *Element2 ;
spREAL *Element3Negated;
spREAL *Element4Negated;
{ spElement *Element1;
spElement *Element2;
spElement *Element3Negated;
spElement *Element4Negated;
};
typedef struct MatrixFrame *MatrixPtr;
/*
@ -249,49 +299,78 @@ typedef struct MatrixFrame *MatrixPtr;
/* Begin function declarations. */
extern void spClear( MatrixPtr );
extern spREAL spCondition( MatrixPtr, spREAL, int* );
extern MatrixPtr spCreate( int, int, int* );
extern void spDeleteRowAndCol( MatrixPtr, int, int );
extern void spDestroy( MatrixPtr);
extern int spElementCount( MatrixPtr );
extern int spError( MatrixPtr );
extern int spFactor( MatrixPtr );
extern int spFileMatrix( MatrixPtr, char *, char *, int, int, int );
extern int spFileStats( MatrixPtr, char *, char * );
extern int spFillinCount( MatrixPtr );
extern int spGetAdmittance( MatrixPtr, int, int, struct spTemplate* );
extern spREAL *spFindElement(MatrixPtr Matrix, int Row, int Col );
extern spREAL *spGetElement(MatrixPtr, int, int );
extern void *spGetInitInfo( spREAL* );
extern int spGetOnes( MatrixPtr, int, int, int, struct spTemplate* );
extern int spGetQuad( MatrixPtr, int, int, int, int, struct spTemplate* );
extern int spGetSize( MatrixPtr, int );
extern int spInitialize(MatrixPtr, int (*pInit)(spREAL*, void *InitInfo, int, int Col));
extern void spInstallInitInfo( spREAL*, void * );
extern spREAL spLargestElement( MatrixPtr );
extern void spMNA_Preorder( MatrixPtr );
extern spREAL spNorm( MatrixPtr );
extern int spOrderAndFactor(MatrixPtr, spREAL*, spREAL, spREAL, int );
extern int spOriginalCount( MatrixPtr);
extern void spPartition( MatrixPtr, int );
extern void spPrint(MatrixPtr, int, int, int );
extern spREAL spPseudoCondition( MatrixPtr );
extern spREAL spRoundoff( MatrixPtr, spREAL );
extern void spScale( MatrixPtr, spREAL*, spREAL* );
extern void spSetComplex( MatrixPtr );
extern void spSetReal( MatrixPtr );
extern void spStripFills( MatrixPtr );
extern void spWhereSingular(MatrixPtr, int*, int* );
extern void spConstMult(MatrixPtr, double);
spcEXTERN void spClear( spMatrix );
spcEXTERN spREAL spCondition( spMatrix, spREAL, int* );
spcEXTERN spMatrix spCreate( int, int, spError* );
spcEXTERN void spDeleteRowAndCol( spMatrix, int, int );
spcEXTERN void spDestroy( spMatrix );
spcEXTERN int spElementCount( spMatrix );
spcEXTERN int spOriginalCount( spMatrix );
spcEXTERN spError spErrorState( spMatrix );
#ifdef EOF
spcEXTERN void spErrorMessage( spMatrix, FILE*, char* );
#else
# define spErrorMessage(a,b,c) spcFUNC_NEEDS_FILE(_spErrorMessage,stdio)
#endif
spcEXTERN spError spFactor( spMatrix );
spcEXTERN int spFileMatrix( spMatrix, char*, char*, int, int, int );
spcEXTERN int spFileStats( spMatrix, char*, char* );
spcEXTERN int spFillinCount( spMatrix );
spcEXTERN spElement *spFindElement( spMatrix, int, int );
spcEXTERN spError spGetAdmittance( spMatrix, int, int,
struct spTemplate* );
spcEXTERN spElement *spGetElement( spMatrix, int, int );
spcEXTERN spGenericPtr spGetInitInfo( spElement* );
spcEXTERN spError spGetOnes( spMatrix, int, int, int,
struct spTemplate* );
spcEXTERN spError spGetQuad( spMatrix, int, int, int, int,
struct spTemplate* );
spcEXTERN int spGetSize( spMatrix, int );
spcEXTERN int spInitialize( spMatrix, int (*pInit)(spElement *, spGenericPtr, int, int) );
spcEXTERN void spInstallInitInfo( spElement*, spGenericPtr );
spcEXTERN spREAL spLargestElement( spMatrix );
spcEXTERN void spMNA_Preorder( spMatrix );
spcEXTERN spREAL spNorm( spMatrix );
spcEXTERN spError spOrderAndFactor( spMatrix, spREAL[], spREAL,
spREAL, int );
spcEXTERN void spPartition( spMatrix, int );
spcEXTERN void spPrint( spMatrix, int, int, int );
spcEXTERN spREAL spPseudoCondition( spMatrix );
spcEXTERN spREAL spRoundoff( spMatrix, spREAL );
spcEXTERN void spScale( spMatrix, spREAL[], spREAL[] );
spcEXTERN void spSetComplex( spMatrix );
spcEXTERN void spSetReal( spMatrix );
spcEXTERN void spStripFills( spMatrix );
spcEXTERN void spWhereSingular( spMatrix, int*, int* );
/* Functions with argument lists that are dependent on options. */
extern void spDeterminant ( MatrixPtr, int*, spREAL*, spREAL* );
extern int spFileVector( MatrixPtr, char * , spREAL*, spREAL*);
extern void spMultiply( MatrixPtr, spREAL*, spREAL*, spREAL*, spREAL* );
extern void spMultTransposed(MatrixPtr,spREAL*,spREAL*,spREAL*,spREAL*);
extern void spSolve( MatrixPtr, spREAL*, spREAL*, spREAL*, spREAL* );
extern void spSolveTransposed(MatrixPtr,spREAL*,spREAL*,spREAL*,spREAL*);
#if spCOMPLEX
spcEXTERN void spDeterminant( spMatrix, int*, spREAL*, spREAL* );
#else /* NOT spCOMPLEX */
spcEXTERN void spDeterminant( spMatrix, int*, spREAL* );
#endif /* NOT spCOMPLEX */
#if spCOMPLEX && spSEPARATED_COMPLEX_VECTORS
spcEXTERN int spFileVector( spMatrix, char* ,
spREAL[], spREAL[]);
spcEXTERN void spMultiply( spMatrix, spREAL[], spREAL[],
spREAL[], spREAL[] );
spcEXTERN void spMultTransposed( spMatrix, spREAL[], spREAL[],
spREAL[], spREAL[] );
spcEXTERN void spSolve( spMatrix, spREAL[], spREAL[], spREAL[],
spREAL[] );
spcEXTERN void spSolveTransposed( spMatrix, spREAL[], spREAL[],
spREAL[], spREAL[] );
#else /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */
spcEXTERN int spFileVector( spMatrix, char* , spREAL[] );
spcEXTERN void spMultiply( spMatrix, spREAL[], spREAL[] );
spcEXTERN void spMultTransposed( spMatrix,
spREAL[], spREAL[] );
spcEXTERN void spSolve( spMatrix, spREAL[], spREAL[] );
spcEXTERN void spSolveTransposed( spMatrix,
spREAL[], spREAL[] );
#endif /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */
#endif /* spOKAY */
//spcEXTERN void spConstMult(spMatrix, double);

18
src/maths/sparse/Makefile.am
View File

@ -3,16 +3,16 @@
noinst_LTLIBRARIES = libsparse.la
libsparse_la_SOURCES = \
spalloc.c \
spbuild.c \
spconfig.h \
spdefs.h \
spAllocate.c \
spBuild.c \
spConfig.h \
spDefs.h \
spextra.c \
spfactor.c \
spoutput.c \
spsmp.c \
spsolve.c \
sputils.c
spFactor.c \
spOutput.c \
spSMP.c \
spSolve.c \
spUtils.c

863
src/maths/sparse/spAllocate.c Normal file
View File

@ -0,0 +1,863 @@
/*
* MATRIX ALLOCATION MODULE
*
* Author: Advising professor:
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
* UC Berkeley
*/
/*!\file
* This file contains functions for allocating and freeing matrices, configuring them, and for
* accessing global information about the matrix (size, error status, etc.).
*
* Objects that begin with the \a spc prefix are considered private
* and should not be used.
*
* \author
* Kenneth S. Kundert <kundert@users.sourceforge.net>
*/
/* >>> User accessible functions contained in this file:
* spCreate
* spDestroy
* spErrorState
* spWhereSingular
* spGetSize
* spSetReal
* spSetComplex
* spFillinCount
* spElementCount
*
* >>> Other functions contained in this file:
* spcGetElement
* InitializeElementBlocks
* spcGetFillin
* RecordAllocation
* AllocateBlockOfAllocationList
* EnlargeMatrix
* ExpandTranslationArrays
*/
/*
* Revision and copyright information.
*
* Copyright (c) 1985-2003 by Kenneth S. Kundert
*/
/*
* 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 <stdio.h>
#include "spConfig.h"
#include "ngspice/spmatrix.h"
#include "spDefs.h"
/*
* Global strings
*/
char spcMatrixIsNotValid[] = "Matrix passed to Sparse is not valid";
char spcErrorsMustBeCleared[] = "Error not cleared";
char spcMatrixMustBeFactored[] = "Matrix must be factored";
char spcMatrixMustNotBeFactored[] = "Matrix must not be factored";
/*
* Function declarations
*/
//static spError ReserveElements( MatrixPtr, int );
static void InitializeElementBlocks( MatrixPtr, int, int );
static void RecordAllocation( MatrixPtr, void* );
static void AllocateBlockOfAllocationList( MatrixPtr );
/*!
* Allocates and initializes the data structures associated with a matrix.
*
* \return
* A pointer to the matrix is returned cast into \a spMatrix (typically a
* pointer to a void). This pointer is then passed and used by the other
* matrix routines to refer to a particular matrix. If an error occurs,
* the \a NULL pointer is returned.
*
* \param Size
* Size of matrix or estimate of size of matrix if matrix is \a EXPANDABLE.
* \param Complex
* Type of matrix. If \a 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 an \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.
* \param pError
* Returns error flag, needed because function \a spErrorState() will
* not work correctly if \a spCreate() returns \a NULL. Possible errors
* include \a spNO_MEMORY and \a spPANIC.
*/
/* >>> Local variables:
* AllocatedSize (int)
* The size of the matrix being allocated.
* Matrix (MatrixPtr)
* A pointer to the matrix frame being created.
*/
spMatrix
spCreate(
int Size,
int Complex,
int *pError
)
{
register unsigned SizePlusOne;
register MatrixPtr Matrix;
register int I;
int AllocatedSize;
/* Begin `spCreate'. */
/* Clear error flag. */
*pError = spOKAY;
/* Test for valid size. */
vASSERT( (Size >= 0) AND (Size != 0 OR EXPANDABLE), "Invalid size" );
/* Test for valid type. */
#if NOT spCOMPLEX
ASSERT( NOT Complex );
#endif
#if NOT REAL
ASSERT( Complex );
#endif
/* Create Matrix. */
AllocatedSize = MAX( Size, MINIMUM_ALLOCATED_SIZE );
SizePlusOne = (unsigned)(AllocatedSize + 1);
if ((Matrix = ALLOC(struct MatrixFrame, 1)) == NULL)
{ *pError = spNO_MEMORY;
return NULL;
}
/* Initialize matrix */
Matrix->ID = SPARSE_ID;
Matrix->Complex = Complex;
Matrix->PreviousMatrixWasComplex = Complex;
Matrix->Factored = NO;
Matrix->Elements = 0;
Matrix->Error = *pError;
Matrix->Fillins = 0;
Matrix->Reordered = NO;
Matrix->NeedsOrdering = YES;
Matrix->NumberOfInterchangesIsOdd = NO;
Matrix->Partitioned = NO;
Matrix->RowsLinked = NO;
Matrix->InternalVectorsAllocated = NO;
Matrix->SingularCol = 0;
Matrix->SingularRow = 0;
Matrix->Size = Size;
Matrix->AllocatedSize = AllocatedSize;
Matrix->ExtSize = Size;
Matrix->AllocatedExtSize = AllocatedSize;
Matrix->CurrentSize = 0;
Matrix->ExtToIntColMap = NULL;
Matrix->ExtToIntRowMap = NULL;
Matrix->IntToExtColMap = NULL;
Matrix->IntToExtRowMap = NULL;
Matrix->MarkowitzRow = NULL;
Matrix->MarkowitzCol = NULL;
Matrix->MarkowitzProd = NULL;
Matrix->DoCmplxDirect = NULL;
Matrix->DoRealDirect = NULL;
Matrix->Intermediate = NULL;
Matrix->RelThreshold = DEFAULT_THRESHOLD;
Matrix->AbsThreshold = 0.0;
Matrix->TopOfAllocationList = NULL;
Matrix->RecordsRemaining = 0;
Matrix->ElementsRemaining = 0;
Matrix->FillinsRemaining = 0;
RecordAllocation( Matrix, (void *)Matrix );
if (Matrix->Error == spNO_MEMORY) goto MemoryError;
/* Take out the trash. */
Matrix->TrashCan.Real = 0.0;
#if spCOMPLEX
Matrix->TrashCan.Imag = 0.0;
#endif
Matrix->TrashCan.Row = 0;
Matrix->TrashCan.Col = 0;
Matrix->TrashCan.NextInRow = NULL;
Matrix->TrashCan.NextInCol = NULL;
#if INITIALIZE
Matrix->TrashCan.pInitInfo = NULL;
#endif
/* Allocate space in memory for Diag pointer vector. */
CALLOC( Matrix->Diag, ElementPtr, SizePlusOne);
if (Matrix->Diag == NULL)
goto MemoryError;
/* Allocate space in memory for FirstInCol pointer vector. */
CALLOC( Matrix->FirstInCol, ElementPtr, SizePlusOne);
if (Matrix->FirstInCol == NULL)
goto MemoryError;
/* Allocate space in memory for FirstInRow pointer vector. */
CALLOC( Matrix->FirstInRow, ElementPtr, SizePlusOne);
if (Matrix->FirstInRow == NULL)
goto MemoryError;
/* Allocate space in memory for IntToExtColMap vector. */
if (( Matrix->IntToExtColMap = ALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Allocate space in memory for IntToExtRowMap vector. */
if (( Matrix->IntToExtRowMap = ALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Initialize MapIntToExt vectors. */
for (I = 1; I <= AllocatedSize; I++)
{ Matrix->IntToExtRowMap[I] = I;
Matrix->IntToExtColMap[I] = I;
}
#if TRANSLATE
/* Allocate space in memory for ExtToIntColMap vector. */
if (( Matrix->ExtToIntColMap = ALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Allocate space in memory for ExtToIntRowMap vector. */
if (( Matrix->ExtToIntRowMap = ALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Initialize MapExtToInt vectors. */
for (I = 1; I <= AllocatedSize; I++)
{ Matrix->ExtToIntColMap[I] = -1;
Matrix->ExtToIntRowMap[I] = -1;
}
Matrix->ExtToIntColMap[0] = 0;
Matrix->ExtToIntRowMap[0] = 0;
#endif
/* Allocate space for fill-ins and initial set of elements. */
InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize,
SPACE_FOR_FILL_INS*AllocatedSize );
if (Matrix->Error == spNO_MEMORY)
goto MemoryError;
return (char *)Matrix;
MemoryError:
/* Deallocate matrix and return no pointer to matrix if there is not enough
memory. */
*pError = spNO_MEMORY;
spDestroy( (char *)Matrix);
return NULL;
}
/*
* ELEMENT ALLOCATION
*
* This routine allocates space for matrix elements. It requests large blocks
* of storage from the system and doles out individual elements as required.
* This technique, as opposed to allocating elements individually, tends to
* speed the allocation process.
*
* >>> Returned:
* A pointer to an element.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to matrix.
*
* >>> Local variables:
* pElement (ElementPtr)
* A pointer to the first element in the group of elements being allocated.
*
* >>> Possible errors:
* spNO_MEMORY
*/
ElementPtr
spcGetElement( MatrixPtr Matrix )
{
ElementPtr pElement;
/* Begin `spcGetElement'. */
/* Allocate block of MatrixElements if necessary. */
if (Matrix->ElementsRemaining == 0)
{ pElement = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, (void *)pElement );
if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailElement = pElement;
}
/* Update Element counter and return pointer to Element. */
Matrix->ElementsRemaining--;
return Matrix->NextAvailElement++;
}
/*
* ELEMENT ALLOCATION INITIALIZATION
*
* This routine allocates space for matrix fill-ins and an initial set of
* elements. Besides being faster than allocating space for elements one
* at a time, it tends to keep the fill-ins physically close to the other
* matrix elements in the computer memory. This keeps virtual memory paging
* to a minimum.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to the matrix.
* InitialNumberOfElements <input> (int)
* This number is used as the size of the block of memory, in
* MatrixElements, reserved for elements. If more than this number of
* elements are generated, then more space is allocated later.
* NumberOfFillinsExpected <input> (int)
* This number is used as the size of the block of memory, in
* MatrixElements, reserved for fill-ins. If more than this number of
* fill-ins are generated, then more space is allocated, but they may
* not be physically close in computer's memory.
*
* >>> Local variables:
* pElement (ElementPtr)
* A pointer to the first element in the group of elements being allocated.
*
* >>> Possible errors:
* spNO_MEMORY
*/
static void
InitializeElementBlocks(
MatrixPtr Matrix,
int InitialNumberOfElements,
int NumberOfFillinsExpected
)
{
ElementPtr pElement;
/* Begin `InitializeElementBlocks'. */
/* Allocate block of MatrixElements for elements. */
pElement = ALLOC(struct MatrixElement, InitialNumberOfElements);
RecordAllocation( Matrix, (void *)pElement );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->ElementsRemaining = InitialNumberOfElements;
Matrix->NextAvailElement = pElement;
/* Allocate block of MatrixElements for fill-ins. */
pElement = ALLOC(struct MatrixElement, NumberOfFillinsExpected);
RecordAllocation( Matrix, (void *)pElement );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->FillinsRemaining = NumberOfFillinsExpected;
Matrix->NextAvailFillin = pElement;
/* Allocate a fill-in list structure. */
Matrix->FirstFillinListNode = ALLOC(struct FillinListNodeStruct,1);
RecordAllocation( Matrix, (void *)Matrix->FirstFillinListNode );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
Matrix->FirstFillinListNode->pFillinList = pElement;
Matrix->FirstFillinListNode->NumberOfFillinsInList =NumberOfFillinsExpected;
Matrix->FirstFillinListNode->Next = NULL;
return;
}
/*
* FILL-IN ALLOCATION
*
* This routine allocates space for matrix fill-ins. It requests large blocks
* of storage from the system and doles out individual elements as required.
* This technique, as opposed to allocating elements individually, tends to
* speed the allocation process.
*
* >>> Returned:
* A pointer to the fill-in.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to matrix.
*
* >>> Possible errors:
* spNO_MEMORY
*/
ElementPtr
spcGetFillin( MatrixPtr Matrix )
{
#if STRIP OR LINT
struct FillinListNodeStruct *pListNode;
ElementPtr pFillins;
#endif
/* Begin `spcGetFillin'. */
#if NOT STRIP OR LINT
if (Matrix->FillinsRemaining == 0)
return spcGetElement( Matrix );
#endif
#if STRIP OR LINT
if (Matrix->FillinsRemaining == 0)
{ pListNode = Matrix->LastFillinListNode;
/* First see if there are any stripped fill-ins left. */
if (pListNode->Next != NULL)
{ Matrix->LastFillinListNode = pListNode = pListNode->Next;
Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
Matrix->NextAvailFillin = pListNode->pFillinList;
}
else
{
/* Allocate block of fill-ins. */
pFillins = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, (void *)pFillins );
if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailFillin = pFillins;
/* Allocate a fill-in list structure. */
pListNode->Next = ALLOC(struct FillinListNodeStruct,1);
RecordAllocation( Matrix, (void *)pListNode->Next );
if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->LastFillinListNode = pListNode = pListNode->Next;
pListNode->pFillinList = pFillins;
pListNode->NumberOfFillinsInList = ELEMENTS_PER_ALLOCATION;
pListNode->Next = NULL;
}
}
#endif
/* Update Fill-in counter and return pointer to Fill-in. */
Matrix->FillinsRemaining--;
return Matrix->NextAvailFillin++;
}
/*
* RECORD A MEMORY ALLOCATION
*
* This routine is used to record all memory allocations so that the memory
* can be freed later.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to the matrix.
* AllocatedPtr <input> (void *)
* The pointer returned by malloc or calloc. These pointers are saved in
* a list so that they can be easily freed.
*
* >>> Possible errors:
* spNO_MEMORY
*/
static void
RecordAllocation(
MatrixPtr Matrix,
void *AllocatedPtr
)
{
/* Begin `RecordAllocation'. */
/*
* If Allocated pointer is NULL, assume that malloc returned a NULL pointer,
* which indicates a spNO_MEMORY error.
*/
if (AllocatedPtr == NULL)
{ Matrix->Error = spNO_MEMORY;
return;
}
/* Allocate block of MatrixElements if necessary. */
if (Matrix->RecordsRemaining == 0)
{ AllocateBlockOfAllocationList( Matrix );
if (Matrix->Error == spNO_MEMORY)
{ FREE(AllocatedPtr);
return;
}
}
/* Add Allocated pointer to Allocation List. */
(++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr;
Matrix->RecordsRemaining--;
return;
}
/*
* ADD A BLOCK OF SLOTS TO ALLOCATION LIST
*
* This routine increases the size of the allocation list.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to the matrix.
*
* >>> Local variables:
* ListPtr (AllocationListPtr)
* Pointer to the list that contains the pointers to segments of memory
* that were allocated by the operating system for the current matrix.
*
* >>> Possible errors:
* spNO_MEMORY
*/
static void
AllocateBlockOfAllocationList( MatrixPtr Matrix )
{
register int I;
register AllocationListPtr ListPtr;
/* Begin `AllocateBlockOfAllocationList'. */
/* Allocate block of records for allocation list. */
ListPtr = ALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION+1));
if (ListPtr == NULL)
{ Matrix->Error = spNO_MEMORY;
return;
}
/* String entries of allocation list into singly linked list. List is linked
such that any record points to the one before it. */
ListPtr->NextRecord = Matrix->TopOfAllocationList;
Matrix->TopOfAllocationList = ListPtr;
ListPtr += ELEMENTS_PER_ALLOCATION;
for (I = ELEMENTS_PER_ALLOCATION; I > 0; I--)
{ ListPtr->NextRecord = ListPtr - 1;
ListPtr--;
}
/* Record allocation of space for allocation list on allocation list. */
Matrix->TopOfAllocationList->AllocatedPtr = (void *)ListPtr;
Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;
return;
}
/*!
* Destroys a matrix and frees all memory associated with it.
*
* \param eMatrix
* Pointer to the matrix frame which is to be destroyed.
*/
/* >>> Local variables:
* ListPtr (AllocationListPtr)
* Pointer into the linked list of pointers to allocated data structures.
* Points to pointer to structure to be freed.
* NextListPtr (AllocationListPtr)
* Pointer into the linked list of pointers to allocated data structures.
* Points to the next pointer to structure to be freed. This is needed
* because the data structure to be freed could include the current node
* in the allocation list.
*/
void
spDestroy( spMatrix eMatrix )
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register AllocationListPtr ListPtr, NextListPtr;
/* Begin `spDestroy'. */
ASSERT_IS_SPARSE( Matrix );
/* Deallocate the vectors that are located in the matrix frame. */
FREE( Matrix->IntToExtColMap );
FREE( Matrix->IntToExtRowMap );
FREE( Matrix->ExtToIntColMap );
FREE( Matrix->ExtToIntRowMap );
FREE( Matrix->Diag );
FREE( Matrix->FirstInRow );
FREE( Matrix->FirstInCol );
FREE( Matrix->MarkowitzRow );
FREE( Matrix->MarkowitzCol );
FREE( Matrix->MarkowitzProd );
FREE( Matrix->DoCmplxDirect );
FREE( Matrix->DoRealDirect );
FREE( Matrix->Intermediate );
/* Sequentially step through the list of allocated pointers freeing pointers
* along the way. */
ListPtr = Matrix->TopOfAllocationList;
while (ListPtr != NULL)
{ NextListPtr = ListPtr->NextRecord;
free( ListPtr->AllocatedPtr );
ListPtr = NextListPtr;
}
return;
}
/*!
* This function returns the error status of the given matrix.
*
* \return
* The error status of the given matrix.
*
* \param eMatrix
* The pointer to the matrix for which the error status is desired.
*/
int
spErrorState( spMatrix eMatrix )
{
/* Begin `spErrorState'. */
if (eMatrix != NULL)
{ ASSERT_IS_SPARSE( (MatrixPtr)eMatrix );
return ((MatrixPtr)eMatrix)->Error;
}
else return spNO_MEMORY; /* This error may actually be spPANIC,
* no way to tell. */
}
/*!
* This function returns the row and column number where the matrix was
* detected as singular (if pivoting was allowed on the last factorization)
* or where a zero was detected on the diagonal (if pivoting was not
* allowed on the last factorization). Pivoting is performed only in
* spOrderAndFactor().
*
* \param eMatrix
* The matrix for which the error status is desired.
* \param pRow
* The row number.
* \param pCol
* The column number.
*/
void
spWhereSingular(
spMatrix eMatrix,
int *pRow,
int *pCol
)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
/* Begin `spWhereSingular'. */
ASSERT_IS_SPARSE( Matrix );
if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG)
{ *pRow = Matrix->SingularRow;
*pCol = Matrix->SingularCol;
}
else *pRow = *pCol = 0;
return;
}
/*!
* Returns the size of the matrix. Either the internal or external size of
* the matrix is returned.
*
* \param eMatrix
* Pointer to matrix.
* \param External
* If \a 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 \a TRANSLATE option is set true.
*/
int
spGetSize(
spMatrix eMatrix,
int External
)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
/* Begin `spGetSize'. */
ASSERT_IS_SPARSE( Matrix );
#if TRANSLATE
if (External)
return Matrix->ExtSize;
else
return Matrix->Size;
#else
return Matrix->Size;
#endif
}
/*!
* Forces matrix to be real.
*
* \param eMatrix
* Pointer to matrix.
*/
void
spSetReal( spMatrix eMatrix )
{
/* Begin `spSetReal'. */
ASSERT_IS_SPARSE( (MatrixPtr)eMatrix );
vASSERT( REAL, "Sparse not compiled to handle real matrices" );
((MatrixPtr)eMatrix)->Complex = NO;
return;
}
/*!
* Forces matrix to be complex.
*
* \param eMatrix
* Pointer to matrix.
*/
void
spSetComplex( spMatrix eMatrix )
{
/* Begin `spSetComplex'. */
ASSERT_IS_SPARSE( (MatrixPtr)eMatrix );
vASSERT( spCOMPLEX, "Sparse not compiled to handle complex matrices" );
((MatrixPtr)eMatrix)->Complex = YES;
return;
}
/*!
* This function returns the number of fill-ins that currently exists in a matrix.
*
* \param eMatrix
* Pointer to matrix.
*/
int
spFillinCount( spMatrix eMatrix )
{
/* Begin `spFillinCount'. */
ASSERT_IS_SPARSE( (MatrixPtr)eMatrix );
return ((MatrixPtr)eMatrix)->Fillins;
}
/*!
* This function returns the total number of elements (including fill-ins) that currently exists in a matrix.
*
* \param eMatrix
* Pointer to matrix.
*/
/* FIXME: Seems no different size entries available anymore */
int
spElementCount( spMatrix eMatrix )
{
/* Begin `spElementCount'. */
ASSERT_IS_SPARSE( (MatrixPtr)eMatrix );
return ((MatrixPtr)eMatrix)->Elements;
}
int
spOriginalCount( spMatrix eMatrix )
{
/* Begin `spOriginalCount'. */
ASSERT_IS_SPARSE( (MatrixPtr)eMatrix );
return ((MatrixPtr)eMatrix)->Elements;
}

1284
src/maths/sparse/spBuild.c Normal file
View File

File diff suppressed because it is too large Load Diff

572
src/maths/sparse/spConfig.h Normal file
View File

@ -0,0 +1,572 @@
/* CONFIGURATION MACRO DEFINITIONS for sparse matrix routines */
/*!
* \file
*
* This file contains macros for the sparse matrix routines that are used
* to define the personality of the routines. The user is expected to
* modify this file to maximize the performance of the routines with
* his/her matrices.
*
* Macros are distinguished by using solely capital letters in their
* identifiers. This contrasts with C defined identifiers which are
* strictly lower case, and program variable and procedure names which use
* both upper and lower case.
*
* Objects that begin with the \a spc prefix are considered private
* and should not be used.
*
* \author
* Kenneth S. Kundert <kundert@users.sourceforge.net>
*/
/*
* Revision and copyright information.
*
* Copyright (c) 1985-2003 by Kenneth S. Kundert
*
* $Date: 2003/06/30 19:41:29 $
* $Revision: 1.5 $
*/
#ifndef spCONFIG_DEFS
#define spCONFIG_DEFS
#ifdef spINSIDE_SPARSE
/*
* OPTIONS
*
* These are compiler options. Set each option to one to compile that
* section of the code. If a feature is not desired, set the macro
* to NO.
*/
/* Begin options. */
/* Arithmetic Precision
*
* The precision of the arithmetic used by Sparse can be set by
* changing changing the spREAL macro. This macro is
* contained in the file spMatrix.h. It is strongly suggested to
* used double precision with circuit simulators. Note that
* because C always performs arithmetic operations in double
* precision, the only benefit to using single precision is that
* less storage is required. There is often a noticeable speed
* penalty when using single precision. Sparse internally refers
* to a spREAL as a RealNumber.
*/
/*!
* This specifies that the routines are expected to handle real
* systems of equations. The routines can be compiled to handle
* both real and complex systems at the same time, but there is a
* slight speed and memory advantage if the routines are complied
* to handle only real systems of equations.
*/
#define REAL YES
/*!
* Setting this compiler flag true (1) makes the matrix
* expandable before it has been factored. If the matrix is
* expandable, then if an element is added that would be
* considered out of bounds in the current matrix, the size of
* the matrix is increased to hold that element. As a result,
* the size of the matrix need not be known before the matrix is
* built. The matrix can be allocated with size zero and expanded.
*/
#define EXPANDABLE YES
/*!
* This option allows the set of external row and column numbers
* to be non-packed. In other words, the row and column numbers
* do not have to be contiguous. The priced paid for this
* flexibility is that when \a TRANSLATE is set true, the time
* required to initially build the matrix will be greater because
* the external row and column number must be translated into
* internal equivalents. This translation brings about other
* benefits though. First, the spGetElement() and
* spGetAdmittance() routines may be used after the matrix has
* been factored. Further, elements, and even rows and columns,
* may be added to the matrix, and row and columns may be deleted
* from the matrix, after it has been factored. Note that when
* the set of row and column number is not a packed set, neither
* are the \a RHS and \a Solution vectors. Thus the size of these
* vectors must be at least as large as the external size, which
* is the value of the largest given row or column numbers.
*/
#define TRANSLATE YES
/*!
* Causes the spInitialize(), spGetInitInfo(), and
* spInstallInitInfo() routines to be compiled. These routines
* allow the user to store and read one pointer in each nonzero
* element in the matrix. spInitialize() then calls a user
* specified function for each structural nonzero in the matrix,
* and includes this pointer as well as the external row and
* column numbers as arguments. This allows the user to write
* custom matrix initialization routines.
*/
#define INITIALIZE NO
/*!
* Many matrices, and in particular node- and modified-node
* admittance matrices, tend to be nearly symmetric and nearly
* diagonally dominant. For these matrices, it is a good idea to
* select pivots from the diagonal. With this option enabled,
* this is exactly what happens, though if no satisfactory pivot
* can be found on the diagonal, an off-diagonal pivot will be
* used. If this option is disabled, Sparse does not
* preferentially search the diagonal. Because of this, Sparse
* has a wider variety of pivot candidates available, and so
* presumably fewer fill-ins will be created. However, the
* initial pivot selection process will take considerably longer.
* If working with node admittance matrices, or other matrices
* with a strong diagonal, it is probably best to use
* \a DIAGONAL_PIVOTING for two reasons. First, accuracy will be
* better because pivots will be chosen from the large diagonal
* elements, thus reducing the chance of growth. Second, a near
* optimal ordering will be chosen quickly. If the class of
* matrices you are working with does not have a strong diagonal,
* do not use \a DIAGONAL_PIVOTING, but consider using a larger
* threshold. When \a DIAGONAL_PIVOTING is turned off, the following
* options and constants are not used: \a MODIFIED_MARKOWITZ,
* \a MAX_MARKOWITZ_TIES, and \a TIES_MULTIPLIER.
*/
#define DIAGONAL_PIVOTING YES
/*!
* This determines whether arrays start at an index of zero or one.
* This option is necessitated by the fact that standard C
* convention dictates that arrays begin with an index of zero but
* the standard mathematic convention states that arrays begin with
* an index of one. So if you prefer to start your arrays with
* zero, or your calling Sparse from FORTRAN, set ARRAY_OFFSET to
* NO or 0. Otherwise, set ARRAY_OFFSET to YES or 1. Note that if
* you use an offset of one, the arrays that you pass to Sparse
* must have an allocated length of one plus the size of the
* matrix. ARRAY_OFFSET must be either 0 or 1, no other offsets
* are valid.
*/
#define ARRAY_OFFSET YES
/*!
* This specifies that the modified Markowitz method of pivot
* selection is to be used. The modified Markowitz method differs
* from standard Markowitz in two ways. First, under modified
* Markowitz, the search for a pivot can be terminated early if a
* adequate (in terms of sparsity) pivot candidate is found.
* Thus, when using modified Markowitz, the initial factorization
* can be faster, but at the expense of a suboptimal pivoting
* order that may slow subsequent factorizations. The second
* difference is in the way modified Markowitz breaks Markowitz
* ties. When two or more elements are pivot candidates and they
* all have the same Markowitz product, then the tie is broken by
* choosing the element that is best numerically. The numerically
* best element is the one with the largest ratio of its magnitude
* to the magnitude of the largest element in the same column,
* excluding itself. The modified Markowitz method results in
* marginally better accuracy. This option is most appropriate
* for use when working with very large matrices where the initial
* factor time represents an unacceptable burden. \a NO is recommended.
*/
#define MODIFIED_MARKOWITZ NO
/*!
* This specifies that the spDeleteRowAndCol() routine
* should be compiled. Note that for this routine to be
* compiled, both \a DELETE and \a TRANSLATE should be set true.
*/
#define DELETE NO
/*!
* This specifies that the spStripFills() routine should be compiled.
*/
#define STRIP NO
/*!
* This specifies that the routine that preorders modified node
* admittance matrices should be compiled. This routine results
* in greater speed and accuracy if used with this type of
* matrix.
*/
#define MODIFIED_NODAL YES
/*!
* This specifies that the routines that allow four related
* elements to be entered into the matrix at once should be
* compiled. These elements are usually related to an
* admittance. The routines affected by \a QUAD_ELEMENT are the
* spGetAdmittance(), spGetQuad() and spGetOnes() routines.
*/
#define QUAD_ELEMENT NO
/*!
* This specifies that the routines that solve the matrix as if
* it was transposed should be compiled. These routines are
* useful when performing sensitivity analysis using the adjoint
* method.
*/
#define TRANSPOSE YES
/*!
* This specifies that the routine that performs scaling on the
* matrix should be complied. Scaling is not strongly
* supported. The routine to scale the matrix is provided, but
* no routines are provided to scale and descale the RHS and
* Solution vectors. It is suggested that if scaling is desired,
* it only be preformed when the pivot order is being chosen [in
* spOrderAndFactor()]. This is the only time scaling has
* an effect. The scaling may then either be removed from the
* solution by the user or the scaled factors may simply be
* thrown away. \a NO is recommended.
*/
#define SCALING NO
/*!
* This specifies that routines that are used to document the
* matrix, such as spPrint() and spFileMatrix(), should be
* compiled.
*/
#define DOCUMENTATION YES
/*!
* This specifies that routines that are used to multily the
* matrix by a vector, such as spMultiply() and spMultTransposed(), should be
* compiled.
*/
#define MULTIPLICATION YES
/*!
* This specifies that the routine spDeterminant() should be complied.
*/
#define DETERMINANT YES
/*!
* This specifies that spLargestElement() and spRoundoff() should
* be compiled. These routines are used to check the stability (and
* hence the quality of the pivoting) of the factorization by
* computing a bound on the size of the element is the matrix
* \f$ E = A - LU \f$. If this bound is very high after applying
* spOrderAndFactor(), then the pivot threshold should be raised.
* If the bound increases greatly after using spFactor(), then the
* matrix should probably be reordered. Recomend \a NO.
*/
#define STABILITY NO
/*!
* This specifies that spCondition() and spNorm(), the code that
* computes a good estimate of the condition number of the matrix,
* should be compiled. Recomend \a NO.
*/
#define CONDITION NO
/*!
* This specifies that spPseudoCondition(), the code that computes
* a crude and easily fooled indicator of ill-conditioning in the
* matrix, should be compiled. Recomend \a NO.
*/
#define PSEUDOCONDITION NO
/*!
* This specifies that the \a FORTRAN interface routines should be
* compiled. When interfacing to \a FORTRAN programs, the \a ARRAY_OFFSET
* options should be set to NO.
*/
#define FORTRAN NO
/*!
* This specifies that additional error checking will be compiled.
* The type of error checked are those that are common when the
* matrix routines are first integrated into a user's program. Once
* the routines have been integrated in and are running smoothly, this
* option should be turned off. \a YES is recommended.
*/
#define DEBUG YES
#endif /* spINSIDE_SPARSE */
/*
* The following options affect Sparse exports and so are exported as a
* side effect. For this reason they use the `sp' prefix. The boolean
* constants YES an NO are not defined in spMatrix.h to avoid conflicts
* with user code, so use 0 for NO and 1 for YES.
*/
/*!
* This specifies that the routines will be complied to handle
* complex systems of equations.
*/
#define spCOMPLEX 1
/*!
* This specifies the format for complex vectors. If this is set
* false then a complex vector is made up of one double sized
* array of RealNumber's in which the real and imaginary numbers
* are placed alternately in the array. In other
* words, the first entry would be Complex[1].Real, then comes
* Complex[1].Imag, then Complex[2].Real, etc. If
* \a spSEPARATED_COMPLEX_VECTORS is set true, then each complex
* vector is represented by two arrays of \a spREALs, one with
* the real terms, the other with the imaginary. \a NO is recommended.
*/
#define spSEPARATED_COMPLEX_VECTORS 1
#ifdef spINSIDE_SPARSE
/*
* MATRIX CONSTANTS
*
* These constants are used throughout the sparse matrix routines. They
* should be set to suit the type of matrix being solved.
*/
/* Begin constants. */
/*!
* The relative threshold used if the user enters an invalid
* threshold. Also the threshold used by spFactor() when
* calling spOrderAndFactor(). The default threshold should
* not be less than or equal to zero nor larger than one.
* 0.001 is recommended.
*/
#define DEFAULT_THRESHOLD 1.0e-3
/*!
* This indicates whether spOrderAndFactor() should use diagonal
* pivoting as default. This issue only arises when
* spOrderAndFactor() is called from spFactor(). \a YES is recommended.
*/
#define DIAG_PIVOTING_AS_DEFAULT YES
/*!
* This number multiplied by the size of the matrix equals the number
* of elements for which memory is initially allocated in spCreate().
* 6 is recommended.
*/
#define SPACE_FOR_ELEMENTS 6
/*!
* This number multiplied by the size of the matrix equals the number
* of elements for which memory is initially allocated and specifically
* reserved for fill-ins in spCreate(). 4 is recommended.
*/
#define SPACE_FOR_FILL_INS 4
/*!
* The number of matrix elements requested from the malloc utility on
* each call to it. Setting this value greater than 1 reduces the
* amount of overhead spent in this system call. On a virtual memory
* machine, its good to allocate slightly less than a page worth of
* elements at a time (or some multiple thereof).
* 31 is recommended.
*/
#define ELEMENTS_PER_ALLOCATION 31
/*!
* The minimum allocated size of a matrix. Note that this does not
* limit the minimum size of a matrix. This just prevents having to
* resize a matrix many times if the matrix is expandable, large and
* allocated with an estimated size of zero. This number should not
* be less than one.
*/
#define MINIMUM_ALLOCATED_SIZE 6
/*!
* The amount the allocated size of the matrix is increased when it
* is expanded.
*/
#define EXPANSION_FACTOR 1.5
/*!
* Some terminology should be defined. The Markowitz row count is the number
* of non-zero elements in a row excluding the one being considered as pivot.
* There is one Markowitz row count for every row. The Markowitz column
* is defined similarly for columns. The Markowitz product for an element
* is the product of its row and column counts. It is a measure of how much
* work would be required on the next step of the factorization if that
* element were chosen to be pivot. A small Markowitz product is desirable.
*
* This number is used for two slightly different things, both of which
* relate to the search for the best pivot. First, it is the maximum
* number of elements that are Markowitz tied that will be sifted
* through when trying to find the one that is numerically the best.
* Second, it creates an upper bound on how large a Markowitz product
* can be before it eliminates the possibility of early termination
* of the pivot search. In other words, if the product of the smallest
* Markowitz product yet found and \a TIES_MULTIPLIER is greater than
* \a MAX_MARKOWITZ_TIES, then no early termination takes place.
* Set \a MAX_MARKOWITZ_TIES to some small value if no early termination of
* the pivot search is desired. An array of RealNumbers is allocated
* of size \a MAX_MARKOWITZ_TIES so it must be positive and shouldn't
* be too large. Active when MODIFIED_MARKOWITZ is 1 (YES).
* 100 is recommended.
* \see TIES_MULTIPLIER
*/
#define MAX_MARKOWITZ_TIES 100
/*!
* Specifies the number of Markowitz ties that are allowed to occur
* before the search for the pivot is terminated early. Set to some
* large value if no early termination of the pivot search is desired.
* This number is multiplied times the Markowitz product to determine
* how many ties are required for early termination. This means that
* more elements will be searched before early termination if a large
* number of fill-ins could be created by accepting what is currently
* considered the best choice for the pivot. Active when
* \a MODIFIED_MARKOWITZ is 1 (YES). Setting this number to zero
* effectively eliminates all pivoting, which should be avoided.
* This number must be positive. \a TIES_MULTIPLIER is also used when
* diagonal pivoting breaks down. 5 is recommended.
* \see MAX_MARKOWITZ_TIES
*/
#define TIES_MULTIPLIER 5
/*!
* Which partition mode is used by spPartition() as default.
* Possibilities include \a spDIRECT_PARTITION (each row used direct
* addressing, best for a few relatively dense matrices),
* \a spINDIRECT_PARTITION (each row used indirect addressing, best
* for a few very sparse matrices), and \a spAUTO_PARTITION (direct or
* indirect addressing is chosen on a row-by-row basis, carries a large
* overhead, but speeds up both dense and sparse matrices, best if there
* is a large number of matrices that can use the same ordering.
*/
#define DEFAULT_PARTITION spAUTO_PARTITION
/*!
* The number of characters per page width. Set to 80 for terminal,
* 132 for line printer. Controls how many columns printed by
* spPrint() per page width.
*/
#define PRINTER_WIDTH 80
#endif /* spINSIDE_SPARSE */
/*
* PORTABILITY MACROS
*/
#ifdef __STDC__
# define spcCONCAT(prefix,suffix) prefix ## suffix
# define spcQUOTE(x) # x
# define spcFUNC_NEEDS_FILE(func,file) \
func ## _requires_ ## file ## _to_be_included_
#else
# define spcCONCAT(prefix,suffix) prefix/**/suffix
# define spcQUOTE(x) "x"
# define spcFUNC_NEEDS_FILE(func,file) \
func/**/_requires_/**/file/**/_to_be_included_
#endif
#if defined(__cplusplus) || defined(c_plusplus)
/*
* Definitions for C++
*/
# define spcEXTERN extern "C"
# define spcNO_ARGS
# define spcCONST const
typedef void *spGenericPtr;
#else
#ifdef __STDC__
/*
* Definitions for ANSI C
*/
# define spcEXTERN extern
# define spcNO_ARGS void
# define spcCONST const
typedef void *spGenericPtr;
# else
/*
* Definitions for K&R C -- ignore function prototypes
*/
# define spcEXTERN extern
# define spcNO_ARGS
# define spcCONST
typedef char *spGenericPtr;
#endif
#endif
#ifdef spINSIDE_SPARSE
/*
* MACHINE CONSTANTS
*
* These numbers must be updated when the program is ported to a new machine.
*/
/* Begin machine constants. */
#include <limits.h>
#include <float.h>
/*! The resolution of spREAL. */
#define MACHINE_RESOLUTION DBL_EPSILON
/*! The largest possible value of spREAL. */
#define LARGEST_REAL DBL_MAX
/*! The smalles possible positive value of spREAL. */
#define SMALLEST_REAL DBL_MIN
/*! The largest possible value of shorts. */
#define LARGEST_SHORT_INTEGER SHRT_MAX
/*! The largest possible value of longs. */
#define LARGEST_LONG_INTEGER LONG_MAX
/* ANNOTATION */
/*!
* This macro changes the amount of annotation produced by the matrix
* routines. The annotation is used as a debugging aid. Change the number
* associated with \a ANNOTATE to change the amount of annotation produced by
* the program. Possible values include \a NONE, \a ON_STRANGE_BEHAVIOR, and
* \a FULL. \a NONE is recommended.
*/
#define ANNOTATE NONE
/*!
* A possible value for \a ANNOTATE. Disables all annotation.
*/
#define NONE 0
/*!
* A possible value for \a ANNOTATE. Causes annotation to be produce
* upon unusual occurances only.
*/
#define ON_STRANGE_BEHAVIOR 1
/*!
* A possible value for \a ANNOTATE. Enables full annotation.
*/
#define FULL 2
#endif /* spINSIDE_SPARSE */
#endif /* spCONFIG_DEFS */

519
src/maths/sparse/spdefs.h → src/maths/sparse/spDefs.h
View File

@ -1,5 +1,3 @@
#ifndef ngspice_SPDEFS_H
#define ngspice_SPDEFS_H
/*
* DATA STRUCTURE AND MACRO DEFINITIONS for Sparse.
*
@ -15,30 +13,67 @@
/*
* Revision and copyright information.
*
* Copyright (c) 1985,86,87,88,89,90
* by Kenneth S. Kundert and the University of California.
* Copyright (c) 1985-2003 by Kenneth S. Kundert
*
* 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.
* $Date: 2003/06/29 04:19:52 $
* $Revision: 1.2 $
*/
/*
* IMPORTS
* If running lint, change some of the compiler options to get a more
* complete inspection.
*/
#include <stdio.h>
#ifdef lint
#undef REAL
#undef spCOMPLEX
#undef EXPANDABLE
#undef TRANSLATE
#undef INITIALIZE
#undef DELETE
#undef STRIP
#undef MODIFIED_NODAL
#undef QUAD_ELEMENT
#undef TRANSPOSE
#undef SCALING
#undef DOCUMENTATION
#undef MULTIPLICATION
#undef DETERMINANT
#undef CONDITION
#undef PSEUDOCONDITION
#undef FORTRAN
#undef DEBUG
#define REAL YES
#define spCOMPLEX YES
#define EXPANDABLE YES
#define TRANSLATE YES
#define INITIALIZE YES
#define DELETE YES
#define STRIP YES
#define MODIFIED_NODAL YES
#define QUAD_ELEMENT YES
#define TRANSPOSE YES
#define SCALING YES
#define DOCUMENTATION YES
#define MULTIPLICATION YES
#define DETERMINANT YES
#define CONDITION YES
#define PSEUDOCONDITION YES
#define FORTRAN YES
#define DEBUG YES
#define LINT YES
#else /* not lint */
#define LINT NO
#endif /* not lint */
#undef ABORT
#undef FREE
@ -46,9 +81,10 @@
* MACRO DEFINITIONS
*
* Macros are distinguished by using solely capital letters in their
* identifiers. This contrasts with C defined identifiers which are
* strictly lower case, and program variable and procedure names
* which use both upper and lower case. */
* identifiers. This contrasts with C defined identifiers which are strictly
* lower case, and program variable and procedure names which use both upper
* and lower case.
*/
/* Begin macros. */
@ -56,15 +92,32 @@
#define BOOLEAN int
#define NO 0
#define YES 1
#define NOT !
#define AND &&
#define OR ||
#define SPARSE_ID 0x772773 /* Arbitrary (is Sparse on phone). */
#define IS_SPARSE(matrix) ((matrix) != NULL && \
(matrix)->ID == SPARSE_ID)
#define IS_VALID(matrix) ((matrix) != NULL && \
(matrix)->ID == SPARSE_ID && \
(matrix)->Error >= spOKAY && \
(matrix)->Error < spFATAL)
#define IS_FACTORED(matrix) ((matrix)->Factored && !(matrix)->NeedsOrdering)
/* NULL pointer */
#ifndef NULL
#define NULL 0
#endif
/* Define macros for validating matrix. */
#define SPARSE_ID 0xDeadBeef /* Arbitrary. */
#define IS_SPARSE(matrix) (((matrix) != NULL) AND \
((matrix)->ID == SPARSE_ID))
#define NO_ERRORS(matrix) (((matrix)->Error >= spOKAY) AND \
((matrix)->Error < spFATAL))
#define IS_FACTORED(matrix) ((matrix)->Factored AND \
NOT (matrix)->NeedsOrdering)
#define ASSERT_IS_SPARSE(matrix) vASSERT( IS_SPARSE(matrix), \
spcMatrixIsNotValid )
#define ASSERT_NO_ERRORS(matrix) vASSERT( NO_ERRORS(matrix), \
spcErrorsMustBeCleared )
#define ASSERT_IS_FACTORED(matrix) vASSERT( IS_FACTORED(matrix), \
spcMatrixMustBeFactored )
#define ASSERT_IS_NOT_FACTORED(matrix) vASSERT( NOT (matrix)->Factored, \
spcMatrixMustNotBeFactored )
/* Macro commands */
/* Macro functions that return the maximum or minimum independent of type. */
@ -78,36 +131,20 @@
#define SQR(a) ((a)*(a))
/* Macro procedure that swaps two entities. */
#define SWAP(type, a, b) \
do { \
type SWAP_macro_local = a; \
a = b; \
b = SWAP_macro_local; \
} while(0)
#define SWAP(type, a, b) {type swapx; swapx = a; a = b; b = swapx;}
/*
* COMPLEX OPERATION MACROS
*/
#ifndef ngspice_COMPLEX_H
/* Real and Complex numbers definition */
#define spREAL double
/* Begin `realNumber'. */
typedef spREAL RealNumber, *RealVector;
/* Begin `ComplexNumber'. */
typedef struct
{ RealNumber Real;
RealNumber Imag;
} ComplexNumber, *ComplexVector;
/* Macro function that returns the approx absolute value of a complex
number. */
/* Macro function that returns the approx absolute value of a complex number. */
#if spCOMPLEX
#define ELEMENT_MAG(ptr) (ABS((ptr)->Real) + ABS((ptr)->Imag))
#define CMPLX_ASSIGN_VALUE(cnum, vReal, vImag) \
{ (cnum).Real = vReal; \
(cnum).Imag = vImag; \
}
#else
#define ELEMENT_MAG(ptr) ((ptr)->Real < 0.0 ? -(ptr)->Real : (ptr)->Real)
#endif
/* Complex assignment statements. */
#define CMPLX_ASSIGN(to,from) \
@ -126,21 +163,12 @@ typedef struct
{ (to).Real = -(from).Real; \
(to).Imag = (from).Imag; \
}
#define CMPLX_CONJ(a) (a).Imag = -(a).Imag
#define CONJUGATE(a) (a).Imag = -(a).Imag
#define CMPLX_NEGATE(a) \
{ (a).Real = -(a).Real; \
(a).Imag = -(a).Imag; \
}
#define CMPLX_NEGATE_SELF(cnum) \
{ (cnum).Real = -(cnum).Real; \
(cnum).Imag = -(cnum).Imag; \
}
/* Macro that returns the approx magnitude (L-1 norm) of a complex number. */
#define CMPLX_1_NORM(a) (ABS((a).Real) + ABS((a).Imag))
@ -148,7 +176,7 @@ typedef struct
#define CMPLX_INF_NORM(a) (MAX (ABS((a).Real),ABS((a).Imag)))
/* Macro function that returns the magnitude (L-2 norm) of a complex number. */
#define CMPLX_2_NORM(a) (hypot((a).Real, (a).Imag))
#define CMPLX_2_NORM(a) (sqrt((a).Real*(a).Real + (a).Imag*(a).Imag))
/* Macro function that performs complex addition. */
#define CMPLX_ADD(to,from_a,from_b) \
@ -156,11 +184,6 @@ typedef struct
(to).Imag = (from_a).Imag + (from_b).Imag; \
}
/* Macro function that performs addition of a complex and a scalar. */
#define CMPLX_ADD_SELF_SCALAR(cnum, scalar) \
{ (cnum).Real += scalar; \
}
/* Macro function that performs complex subtraction. */
#define CMPLX_SUBT(to,from_a,from_b) \
{ (to).Real = (from_a).Real - (from_b).Real; \
@ -178,7 +201,7 @@ typedef struct
{ (to).Real -= (from).Real; \
(to).Imag -= (from).Imag; \
}
/* Macro function that multiplies a complex number by a scalar. */
#define SCLR_MULT(to,sclr,cmplx) \
{ (to).Real = (sclr) * (cmplx).Real; \
@ -198,32 +221,13 @@ typedef struct
(to).Imag = (from_a).Real * (from_b).Imag + \
(from_a).Imag * (from_b).Real; \
}
/* Macro function that multiplies a complex number and a scalar. */
#define CMPLX_MULT_SCALAR(to,from, scalar) \
{ (to).Real = (from).Real * scalar; \
(to).Imag = (from).Imag * scalar; \
}
/* Macro function that implements *= for a complex and a scalar number. */
#define CMPLX_MULT_SELF_SCALAR(cnum, scalar) \
{ (cnum).Real *= scalar; \
(cnum).Imag *= scalar; \
}
/* Macro function that multiply-assigns a complex number by a scalar. */
#define SCLR_MULT_ASSIGN(to,sclr) \
{ (to).Real *= (sclr); \
(to).Imag *= (sclr); \
}
/* Macro function that implements to *= from for complex numbers. */
#define CMPLX_MULT_ASSIGN(to,from) \
{ RealNumber to_Real_ = (to).Real; \
(to).Real = to_Real_ * (from).Real - \
{ RealNumber to_real_ = (to).Real; \
(to).Real = to_real_ * (from).Real - \
(to).Imag * (from).Imag; \
(to).Imag = to_Real_ * (from).Imag + \
(to).Imag = to_real_ * (from).Imag + \
(to).Imag * (from).Real; \
}
@ -309,8 +313,8 @@ typedef struct
/* Complex division: to = num / den */
#define CMPLX_DIV(to,num,den) \
{ RealNumber r_, s_; \
if (((den).Real >= (den).Imag && (den).Real > -(den).Imag) || \
((den).Real < (den).Imag && (den).Real <= -(den).Imag)) \
if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \
((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \
{ r_ = (den).Imag / (den).Real; \
s_ = (den).Real + r_*(den).Imag; \
(to).Real = ((num).Real + r_*(num).Imag)/s_; \
@ -327,8 +331,8 @@ typedef struct
/* Complex division and assignment: num /= den */
#define CMPLX_DIV_ASSIGN(num,den) \
{ RealNumber r_, s_, t_; \
if (((den).Real >= (den).Imag && (den).Real > -(den).Imag) || \
((den).Real < (den).Imag && (den).Real <= -(den).Imag)) \
if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \
((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \
{ r_ = (den).Imag / (den).Real; \
s_ = (den).Real + r_*(den).Imag; \
t_ = ((num).Real + r_*(num).Imag)/s_; \
@ -347,8 +351,8 @@ typedef struct
/* Complex reciprocation: to = 1.0 / den */
#define CMPLX_RECIPROCAL(to,den) \
{ RealNumber r_; \
if (((den).Real >= (den).Imag && (den).Real > -(den).Imag) || \
((den).Real < (den).Imag && (den).Real <= -(den).Imag)) \
if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \
((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \
{ r_ = (den).Imag / (den).Real; \
(to).Imag = -r_*((to).Real = 1.0/((den).Real + r_*(den).Imag)); \
} \
@ -360,40 +364,184 @@ typedef struct
/* Allocation */
extern void * tmalloc(size_t);
extern void txfree(const void *);
extern void * trealloc(const void *, size_t);
#define SP_MALLOC(type,number) (type *) tmalloc((size_t)(number) * sizeof(type))
#define SP_REALLOC(ptr,type,number) \
ptr = (type *) trealloc(ptr, (size_t)(number) * sizeof(type))
#define SP_FREE(ptr) { if ((ptr) != NULL) txfree(ptr); (ptr) = NULL; }
#include "ngspice/config.h"
/* A new calloc */
#ifndef HAVE_LIBGC
#define SP_CALLOC(ptr,type,number) \
{ ptr = (type *) calloc((size_t)(number), sizeof(type)); \
}
#else /* HAVE_LIBCG */
#define SP_CALLOC(ptr,type,number) \
{ ptr = (type *) tmalloc((size_t)(number) * sizeof(type)); \
}
#include <gc/gc.h>
#define tmalloc(m) GC_malloc(m)
#define trealloc(m, n) GC_realloc((m), (n))
#define tfree(m)
#define txfree(m)
#endif
#include "ngspice/defines.h"
/*
* ASSERT and ABORT
*
* Macro used to assert that if the code is working correctly, then
* a condition must be true. If not, then execution is terminated
* and an error message is issued stating that there is an internal
* error and giving the file and line number. These assertions are
* not evaluated unless the DEBUG flag is true.
*/
#if DEBUG
#define ASSERT(condition) \
{ if (NOT(condition)) \
{ (void)fflush(stdout); \
(void)fprintf(stderr, "sparse: internal error detected in file `%s' at line %d.\n assertion `%s' failed.\n",\
__FILE__, __LINE__, spcQUOTE(condition) ); \
(void)fflush(stderr); \
abort(); \
} \
}
#else
#define ASSERT(condition)
#endif
#if DEBUG
#define vASSERT(condition,message) \
{ if (NOT(condition)) \
vABORT(message); \
}
#else
#define vASSERT(condition,message)
#endif
#if DEBUG
#define vABORT(message) \
{ (void)fflush(stdout); \
(void)fprintf(stderr, "sparse: internal error detected in file `%s' at line %d.\n %s.\n", __FILE__, __LINE__, message );\
(void)fflush(stderr); \
abort(); \
}
#define ABORT() \
{ (void)fflush(stdout); \
(void)fprintf(stderr, "sparse: internal error detected in file `%s' at line %d.\n", __FILE__, __LINE__ ); \
(void)fflush(stderr); \
abort(); \
}
#else
#define vABORT(message) abort()
#define ABORT() abort()
#endif
/*
* IMAGINARY VECTORS
*
* The imaginary vectors iRHS and iSolution are only needed when the
* options spCOMPLEX and spSEPARATED_COMPLEX_VECTORS are set. The following
* macro makes it easy to include or exclude these vectors as needed.
*/
#if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
#define IMAG_VECTORS , iRHS, iSolution
#define IMAG_RHS , iRHS
#define IMAG_RHS_DECL , RealVector iRHS
#define IMAG_VECT_DECL , RealVector iRHS, RealVector iSolution
#else
#define IMAG_VECTORS
#define IMAG_RHS
#define IMAG_RHS_DECL
#define IMAG_VECT_DECL
#endif
/*
* MEMORY ALLOCATION
*/
#include <stddef.h>
spcEXTERN void *malloc(size_t size);
spcEXTERN void *calloc(size_t nmemb, size_t size);
spcEXTERN void *realloc(void *ptr, size_t size);
spcEXTERN void free(void *ptr);
spcEXTERN void abort(void);
#define ALLOC(type,number) ((type *)malloc((unsigned)(sizeof(type)*(number))))
#define REALLOC(ptr,type,number) \
ptr = (type *)realloc((char *)ptr,(unsigned)(sizeof(type)*(number)))
#define FREE(ptr) { if ((ptr) != NULL) free((char *)(ptr)); (ptr) = NULL; }
/* Calloc that properly handles allocating a cleared vector. */
#define CALLOC(ptr,type,number) \
{ int i; ptr = ALLOC(type, number); \
if (ptr != (type *)NULL) \
for(i=(number)-1;i>=0; i--) ptr[i] = (type) 0; \
}
/*
* Utility Functions
*/
/*
* Compute the product of two intergers while avoiding overflow.
* Used when computing Markowitz products.
*/
#define spcMarkoProd(product, op1, op2) \
if (( (op1) > LARGEST_SHORT_INTEGER AND (op2) != 0) OR \
( (op2) > LARGEST_SHORT_INTEGER AND (op1) != 0)) \
{ double fProduct = (double)(op1) * (double)(op2); \
if (fProduct >= LARGEST_LONG_INTEGER) \
(product) = LARGEST_LONG_INTEGER; \
else \
(product) = (long)fProduct; \
} \
else (product) = (op1)*(op2);
/*
* REAL NUMBER
*/
/* Begin `RealNumber'. */
typedef spREAL RealNumber, *RealVector;
/*
* COMPLEX NUMBER DATA STRUCTURE
*
* >>> Structure fields:
* Real (RealNumber)
* The real portion of the number. Real must be the first
* field in this structure.
* Imag (RealNumber)
* The imaginary portion of the number. This field must follow
* immediately after Real.
*/
/* Begin `ComplexNumber'. */
typedef struct
{ RealNumber Real;
RealNumber Imag;
} ComplexNumber, *ComplexVector;
/*
@ -427,7 +575,7 @@ extern void * trealloc(const void *, size_t);
* NextInCol contains a pointer to the next element in the column below
* this element. If this element is the last nonzero in the column then
* NextInCol contains NULL.
* pInitInfo (void *)
* pInitInfo (spGenericPtr)
* Pointer to user data used for initialization of the matrix element.
* Initialized to NULL.
*
@ -442,15 +590,16 @@ extern void * trealloc(const void *, size_t);
/* Begin `MatrixElement'. */
struct MatrixElement
{
RealNumber Real;
{ RealNumber Real;
#if spCOMPLEX
RealNumber Imag;
#endif
int Row;
int Col;
struct MatrixElement *NextInRow;
struct MatrixElement *NextInCol;
#if INITIALIZE
void *pInitInfo;
spGenericPtr pInitInfo;
#endif
};
@ -482,8 +631,7 @@ typedef ElementPtr *ArrayOfElementPtrs;
/* Begin `AllocationRecord'. */
struct AllocationRecord
{
void *AllocatedPtr;
{ void *AllocatedPtr;
struct AllocationRecord *NextRecord;
};
@ -518,21 +666,11 @@ typedef struct AllocationRecord *AllocationListPtr;
/* Begin `FillinListNodeStruct'. */
struct FillinListNodeStruct
{
ElementPtr pFillinList;
{ ElementPtr pFillinList;
int NumberOfFillinsInList;
struct FillinListNodeStruct *Next;
};
/* Similar to above, but keeps track of the original Elements */
/* Begin `ElementListNodeStruct'. */
struct ElementListNodeStruct
{
ElementPtr pElementList;
int NumberOfElementsInList;
struct ElementListNodeStruct *Next;
};
@ -566,7 +704,7 @@ struct ElementListNodeStruct
* grow to when EXPANDABLE is set true and AllocatedSize is the largest
* the matrix can get without requiring that the matrix frame be
* reallocated.
* Complex (int)
* Complex (BOOLEAN)
* The flag which indicates whether the matrix is complex (true) or
* real.
* CurrentSize (int)
@ -575,18 +713,19 @@ struct ElementListNodeStruct
* rows and columns that have elements in them.
* Diag (ArrayOfElementPtrs)
* Array of pointers that points to the diagonal elements.
* DoCmplxDirect (int *)
* DoCmplxDirect (BOOLEAN *)
* Array of flags, one for each column in matrix. If a flag is true
* then corresponding column in a complex matrix should be eliminated
* in spFactor() using direct addressing (rather than indirect
* addressing).
* DoRealDirect (int *)
* DoRealDirect (BOOLEAN *)
* Array of flags, one for each column in matrix. If a flag is true
* then corresponding column in a real matrix should be eliminated
* in spFactor() using direct addressing (rather than indirect
* addressing).
* Elements (int)
* The total number of elements present in matrix.
* The number of original elements (total elements minus fill ins)
* present in matrix.
* Error (int)
* The error status of the sparse matrix package.
* ExtSize (int)
@ -597,7 +736,7 @@ struct ElementListNodeStruct
* ExtToIntRowMap (int [])
* An array that is used to convert external row numbers to internal
* external row numbers. Present only if TRANSLATE option is set true.
* Factored (int)
* Factored (BOOLEAN)
* Indicates if matrix has been factored. This flag is set true in
* spFactor() and spOrderAndFactor() and set false in spCreate()
* and spClear().
@ -618,8 +757,8 @@ struct ElementListNodeStruct
* array used during forward and backward substitution. It is
* commonly called y when the forward and backward substitution process is
* denoted Ax = b => Ly = b and Ux = y.
* InternalVectorsAllocated (int)
* A flag that indicates whether the Markowitz vectors and the
* InternalVectorsAllocated (BOOLEAN)
* A flag that indicates whether theMmarkowitz vectors and the
* Intermediate vector have been created.
* These vectors are created in spcCreateInternalVectors().
* IntToExtColMap (int [])
@ -642,19 +781,16 @@ struct ElementListNodeStruct
* The maximum number of off-diagonal element in the rows of L, the
* lower triangular matrix. This quantity is used when computing an
* estimate of the roundoff error in the matrix.
* NeedsOrdering (int)
* NeedsOrdering (BOOLEAN)
* This is a flag that signifies that the matrix needs to be ordered
* or reordered. NeedsOrdering is set true in spCreate() and
* spGetElement() or spGetAdmittance() if new elements are added to the
* matrix after it has been previously factored. It is set false in
* spOrderAndFactor().
* NumberOfInterchangesIsOdd (int)
* NumberOfInterchangesIsOdd (BOOLEAN)
* Flag that indicates the sum of row and column interchange counts
* is an odd number. Used when determining the sign of the determinant.
* Originals (int)
* The number of original elements (total elements minus fill ins)
* present in matrix.
* Partitioned (int)
* Partitioned (BOOLEAN)
* This flag indicates that the columns of the matrix have been
* partitioned into two groups. Those that will be addressed directly
* and those that will be addressed indirectly in spFactor().
@ -664,7 +800,7 @@ struct ElementListNodeStruct
* Row pivot was chosen from.
* PivotSelectionMethod (char)
* Character that indicates which pivot search method was successful.
* PreviousMatrixWasComplex (int)
* PreviousMatrixWasComplex (BOOLEAN)
* This flag in needed to determine how to clear the matrix. When
* dealing with real matrices, it is important that the imaginary terms
* in the matrix elements be zero. Thus, if the previous matrix was
@ -673,11 +809,11 @@ struct ElementListNodeStruct
* RelThreshold (RealNumber)
* The magnitude an element must have relative to others in its row
* to be considered as a pivot candidate, except as a last resort.
* Reordered (int)
* Reordered (BOOLEAN)
* This flag signifies that the matrix has been reordered. It
* is cleared in spCreate(), set in spMNA_Preorder() and
* spOrderAndFactor() and is used in spPrint().
* RowsLinked (int)
* RowsLinked (BOOLEAN)
* A flag that indicates whether the row pointers exist. The AddByIndex
* routines do not generate the row pointers, which are needed by some
* of the other routines, such as spOrderAndFactor() and spScale().
@ -734,44 +870,42 @@ struct ElementListNodeStruct
/* Begin `MatrixFrame'. */
struct MatrixFrame
{
RealNumber AbsThreshold;
{ RealNumber AbsThreshold;
int AllocatedSize;
int AllocatedExtSize;
int Complex;
BOOLEAN Complex;
int CurrentSize;
ArrayOfElementPtrs Diag;
int *DoCmplxDirect;
int *DoRealDirect;
BOOLEAN *DoCmplxDirect;
BOOLEAN *DoRealDirect;
int Elements;
int Error;
int ExtSize;
int *ExtToIntColMap;
int *ExtToIntRowMap;
int Factored;
BOOLEAN Factored;
int Fillins;
ArrayOfElementPtrs FirstInCol;
ArrayOfElementPtrs FirstInRow;
unsigned long ID;
RealVector Intermediate;
int InternalVectorsAllocated;
BOOLEAN InternalVectorsAllocated;
int *IntToExtColMap;
int *IntToExtRowMap;
int *MarkowitzRow;
int *MarkowitzCol;
long *MarkowitzProd;
int MaxRowCountInLowerTri;
int NeedsOrdering;
int NumberOfInterchangesIsOdd;
int Originals;
int Partitioned;
BOOLEAN NeedsOrdering;
BOOLEAN NumberOfInterchangesIsOdd;
BOOLEAN Partitioned;
int PivotsOriginalCol;
int PivotsOriginalRow;
char PivotSelectionMethod;
int PreviousMatrixWasComplex;
BOOLEAN PreviousMatrixWasComplex;
RealNumber RelThreshold;
int Reordered;
int RowsLinked;
BOOLEAN Reordered;
BOOLEAN RowsLinked;
int SingularCol;
int SingularRow;
int Singletons;
@ -782,29 +916,32 @@ struct MatrixFrame
int RecordsRemaining;
ElementPtr NextAvailElement;
int ElementsRemaining;
struct ElementListNodeStruct *FirstElementListNode;
struct ElementListNodeStruct *LastElementListNode;
ElementPtr NextAvailFillin;
int FillinsRemaining;
struct FillinListNodeStruct *FirstFillinListNode;
struct FillinListNodeStruct *LastFillinListNode;
};
typedef struct MatrixFrame *MatrixPtr;
/*
* Function declarations
* Declarations
*/
extern ElementPtr spcGetElement( MatrixPtr );
extern ElementPtr spcGetFillin( MatrixPtr );
extern ElementPtr spcFindElementInCol( MatrixPtr, ElementPtr*, int, int, int );
extern ElementPtr spcCreateElement( MatrixPtr, int, int, ElementPtr*, int );
extern void spcCreateInternalVectors( MatrixPtr );
extern void spcLinkRows( MatrixPtr );
extern void spcColExchange( MatrixPtr, int, int );
extern void spcRowExchange( MatrixPtr, int, int );
void spErrorMessage(MatrixPtr, FILE *, char *);
#endif
spcEXTERN ElementPtr spcGetElement( MatrixPtr );
spcEXTERN ElementPtr spcGetFillin( MatrixPtr );
spcEXTERN ElementPtr spcFindElementInCol( MatrixPtr, ElementPtr*, int, int, int );
spcEXTERN ElementPtr spcFindDiag( MatrixPtr, int );
spcEXTERN ElementPtr spcCreateElement( MatrixPtr, int, int,
ElementPtr*, ElementPtr*, int );
spcEXTERN void spcCreateInternalVectors( MatrixPtr );
spcEXTERN void spcLinkRows( MatrixPtr );
spcEXTERN void spcColExchange( MatrixPtr, int, int );
spcEXTERN void spcRowExchange( MatrixPtr, int, int );
spcEXTERN char spcMatrixIsNotValid[];
spcEXTERN char spcErrorsMustBeCleared[];
spcEXTERN char spcMatrixMustBeFactored[];
spcEXTERN char spcMatrixMustNotBeFactored[];

2105
src/maths/sparse/spfactor.c → src/maths/sparse/spFactor.c
View File

File diff suppressed because it is too large Load Diff

593
src/maths/sparse/spoutput.c → src/maths/sparse/spOutput.c
View File

@ -4,11 +4,19 @@
* Author: Advisor:
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
* UC Berkeley
*
*/
/*! \file
*
* This file contains the output-to-file and output-to-screen routines for
* the matrix package.
*
* >>> User accessible functions contained in this file:
* Objects that begin with the \a spc prefix are considered private
* and should not be used.
*
* \author
* Kenneth S. Kundert <kundert@users.sourceforge.net>
*/
/* >>> User accessible functions contained in this file:
* spPrint
* spFileMatrix
* spFileVector
@ -21,19 +29,13 @@
/*
* Revision and copyright information.
*
* Copyright (c) 1985,86,87,88,89,90
* 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.
* Copyright (c) 1985-2003
* by Kenneth S. Kundert
*/
/*
* IMPORTS
*
@ -45,53 +47,44 @@
* spDefs.h
* Matrix type and macro definitions for the sparse matrix routines.
*/
#include <assert.h>
#include <stdlib.h>
#define spINSIDE_SPARSE
#include "spconfig.h"
#include <stdio.h>
#include "spConfig.h"
#include "ngspice/spmatrix.h"
#include "spdefs.h"
#include "spDefs.h"
int Printer_Width = PRINTER_WIDTH;
#include "ngspice/config.h"
#ifdef HAS_WINGUI
#include "ngspice/wstdio.h"
#endif
#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 <input> (char *)
* \param eMatrix
* Pointer to matrix.
* PrintReordered <input> (int)
* \param PrintReordered
* 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 <input> (int)
* Boolean flag that when FALSE indicates that output should be
* \param Data
* 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 <input> (int)
* \param Header
* Flag indicating that extra information should be given, such as row
* and column numbers.
*
* >>> Local variables:
*/
/* >>> Local variables:
* Col (int)
* Column being printed.
* ElementCount (int)
@ -135,75 +128,57 @@ int Printer_Width = PRINTER_WIDTH;
*/
void
spPrint(MatrixPtr Matrix, int PrintReordered, int Data, int Header)
spPrint(
spMatrix eMatrix,
int PrintReordered,
int Data,
int Header
)
{
int J = 0;
int I, Row, Col, Size, Top;
int StartCol = 1, StopCol, Columns, ElementCount = 0;
double Magnitude;
double SmallestDiag = 0;
double SmallestElement = 0;
double LargestElement = 0.0, LargestDiag = 0.0;
ElementPtr pElement, *pImagElements;
int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
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 ) );
/* Begin `spPrint'. */
ASSERT_IS_SPARSE( Matrix );
Size = Matrix->Size;
SP_CALLOC(pImagElements, ElementPtr, Printer_Width / 10 + 1);
if ( pImagElements == NULL)
{
Matrix->Error = spNO_MEMORY;
SP_FREE(pImagElements);
return;
}
/* Create a packed external to internal row and column translation
array. */
/* Create a packed external to internal row and column translation array. */
# if TRANSLATE
Top = Matrix->AllocatedExtSize;
#else
Top = Matrix->AllocatedSize;
#endif
SP_CALLOC( PrintOrdToIntRowMap, int, Top + 1 );
if ( PrintOrdToIntRowMap == NULL)
{
Matrix->Error = spNO_MEMORY;
SP_FREE(pImagElements);
return;
}
SP_CALLOC( PrintOrdToIntColMap, int, Top + 1 );
if (PrintOrdToIntColMap == NULL)
{
Matrix->Error = spNO_MEMORY;
SP_FREE(pImagElements);
SP_FREE(PrintOrdToIntRowMap);
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;
{ PrintOrdToIntRowMap[ Matrix->IntToExtRowMap[I] ] = I;
PrintOrdToIntColMap[ Matrix->IntToExtColMap[I] ] = I;
}
/* Pack the arrays. */
/* Pack the arrays. */
for (J = 1, I = 1; I <= Top; I++)
{
if (PrintOrdToIntRowMap[I] != 0)
{ if (PrintOrdToIntRowMap[I] != 0)
PrintOrdToIntRowMap[ J++ ] = PrintOrdToIntRowMap[ I ];
}
for (J = 1, I = 1; I <= Top; I++)
{
if (PrintOrdToIntColMap[I] != 0)
{ if (PrintOrdToIntColMap[I] != 0)
PrintOrdToIntColMap[ J++ ] = PrintOrdToIntColMap[ I ];
}
/* Print header. */
/* Print header. */
if (Header)
{
printf("MATRIX SUMMARY\n\n");
{ printf("MATRIX SUMMARY\n\n");
printf("Size of matrix = %1d x %1d.\n", Size, Size);
if ( Matrix->Reordered && PrintReordered )
if ( Matrix->Reordered AND PrintReordered )
printf("Matrix has been reordered.\n");
putchar('\n');
@ -215,31 +190,31 @@ spPrint(MatrixPtr Matrix, int PrintReordered, int Data, int Header)
SmallestElement = LARGEST_REAL;
SmallestDiag = SmallestElement;
}
if (Size == 0) return;
/* Determine how many columns to use. */
Columns = Printer_Width;
/* 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. */
/*
* Print matrix by printing groups of complete columns until all the columns
* are printed.
*/
J = 0;
while ( J <= Size )
{
/* Calculatestrchr of last column to printed in this group. */
StopCol = StartCol + Columns - 1;
/* Calculate index of last column to printed in this group. */
{ StopCol = StartCol + Columns - 1;
if (StopCol > Size)
StopCol = Size;
/* Label the columns. */
/* Label the columns. */
if (Header)
{
if (Data)
{
printf(" ");
{ if (Data)
{ printf(" ");
for (I = StartCol; I <= StopCol; I++)
{
if (PrintReordered)
{ if (PrintReordered)
Col = I;
else
Col = PrintOrdToIntColMap[I];
@ -248,70 +223,64 @@ spPrint(MatrixPtr Matrix, int PrintReordered, int Data, int Header)
printf("\n\n");
}
else
{
if (PrintReordered)
{ 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] ]);
{ printf("Columns %1d to %1d.\n",
Matrix->IntToExtColMap[ PrintOrdToIntColMap[StartCol] ],
Matrix->IntToExtColMap[ PrintOrdToIntColMap[StopCol] ]);
}
}
}
/* Print every row ... */
/* Print every row ... */
for (I = 1; I <= Size; I++)
{
if (PrintReordered)
{ if (PrintReordered)
Row = I;
else
Row = PrintOrdToIntRowMap[I];
if (Header)
{
if (PrintReordered && !Data)
{ if (PrintReordered AND NOT Data)
printf("%4d", I);
else
printf("%4d", Matrix->IntToExtRowMap[ Row ]);
if (!Data) putchar(' ');
if (NOT Data) putchar(' ');
}
/* ... in each column of the group. */
/* ... in each column of the group. */
for (J = StartCol; J <= StopCol; J++)
{
if (PrintReordered)
{ if (PrintReordered)
Col = J;
else
Col = PrintOrdToIntColMap[J];
pElement = Matrix->FirstInCol[Col];
while(pElement != NULL && pElement->Row != Row)
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.3g", pElement->Real);
/* Case where element exists */
{ if (Data)
printf(" %9.3g", (double)pElement->Real);
else
putchar('x');
/* Update status variables */
/* Update status variables */
if ( (Magnitude = ELEMENT_MAG(pElement)) > LargestElement )
LargestElement = Magnitude;
if ((Magnitude < SmallestElement) && (Magnitude != 0.0))
if ((Magnitude < SmallestElement) AND (Magnitude != 0.0))
SmallestElement = Magnitude;
ElementCount++;
}
/* Case where element is structurally zero */
/* Case where element is structurally zero */
else
{
if (Data)
{ if (Data)
printf(" ...");
else
putchar('.');
@ -319,68 +288,61 @@ spPrint(MatrixPtr Matrix, int PrintReordered, int Data, int Header)
}
putchar('\n');
if (Matrix->Complex && Data)
{
printf(" ");
#if spCOMPLEX
if (Matrix->Complex AND Data)
{ if (Header)
printf(" ");
for (J = StartCol; J <= StopCol; J++)
{
if (pImagElements[J - StartCol] != NULL)
{
printf(" %8.2gj",
pImagElements[J-StartCol]->Imag);
{ if (pImagElements[J - StartCol] != NULL)
{ printf(" %8.2gj",
(double)pImagElements[J-StartCol]->Imag);
}
else printf(" ");
}
putchar('\n');
}
#endif /* spCOMPLEX */
}
/* Calculatestrchr of first column in next group. */
/* Calculate index of first column in next group. */
StartCol = StopCol;
StartCol++;
putchar('\n');
}
if (Header)
{
printf("\nLargest element in matrix = %-1.4g.\n", LargestElement);
{ printf("\nLargest element in matrix = %-1.4g.\n", LargestElement);
printf("Smallest element in matrix = %-1.4g.\n", SmallestElement);
/* Search for largest and smallest diagonal values */
/* 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 (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 */
/* Print the largest and smallest diagonal values */
if ( Matrix->Factored )
{
printf("\nLargest diagonal element = %-1.4g.\n", LargestDiag);
{ printf("\nLargest diagonal element = %-1.4g.\n", LargestDiag);
printf("Smallest diagonal element = %-1.4g.\n", SmallestDiag);
}
else
{
printf("\nLargest pivot element = %-1.4g.\n", LargestDiag);
{ printf("\nLargest pivot element = %-1.4g.\n", LargestDiag);
printf("Smallest pivot element = %-1.4g.\n", SmallestDiag);
}
/* Calculate and print sparsity and number of fill-ins created. */
printf("\nDensity = %2.2f%%.\n", ((double)(ElementCount * 100)) /
((double)(Size * Size)));
printf("Number of originals = %1d.\n", Matrix->Originals);
if (!Matrix->NeedsOrdering)
/* Calculate and print sparsity and number of fill-ins created. */
printf("\nDensity = %2.2f%%.\n", ((double)ElementCount * 100.0)
/ (((double)Size * (double)Size)));
if (NOT Matrix->NeedsOrdering)
printf("Number of fill-ins = %1d.\n", Matrix->Fillins);
}
putchar('\n');
(void)fflush(stdout);
SP_FREE(PrintOrdToIntColMap);
SP_FREE(PrintOrdToIntRowMap);
FREE(PrintOrdToIntColMap);
FREE(PrintOrdToIntRowMap);
return;
}
@ -394,36 +356,33 @@ spPrint(MatrixPtr Matrix, int PrintReordered, int Data, int Header)
/*
* OUTPUT MATRIX TO FILE
*
/*!
* Writes matrix to file in format suitable to be read back in by the
* matrix test program.
*
* >>> Returns:
* \return
* One is returned if routine was successful, otherwise zero is returned.
* The calling function can query errno (the system global error variable)
* The calling function can query \a errno (the system global error variable)
* as to the reason why this routine failed.
*
* >>> Arguments:
* Matrix <input> (char *)
* \param eMatrix
* Pointer to matrix.
* File <input> (char *)
* \param File
* Name of file into which matrix is to be written.
* Label <input> (char *)
* \param Label
* String that is transferred to file and is used as a label.
* Reordered <input> (int)
* \param Reordered
* Specifies whether matrix should be output in reordered form,
* or in original order.
* Data <input> (int)
* \param Data
* Indicates that the element values should be output along with
* the indices for each element. This parameter must be TRUE if
* the indices for each element. This parameter must be true if
* matrix is to be read by the sparse test program.
* Header <input> (int)
* Indicates that header is desired. This parameter must be TRUE if
* \param Header
* Indicates that header is desired. This parameter must be true if
* matrix is to be read by the sparse test program.
*
* >>> Local variables:
*/
/* >>> Local variables:
* Col (int)
* The original column number of the element being output.
* pElement (ElementPtr)
@ -437,123 +396,118 @@ spPrint(MatrixPtr Matrix, int PrintReordered, int Data, int Header)
*/
int
spFileMatrix(MatrixPtr Matrix, char *File, char *Label, int Reordered,
int Data, int Header)
spFileMatrix(
spMatrix eMatrix,
char *File,
char *Label,
int Reordered,
int Data,
int Header
)
{
int I, Size;
ElementPtr pElement;
int Row, Col, Err;
FILE *pMatrixFile;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register int I, Size;
register ElementPtr pElement;
int Row, Col, Err;
FILE *pMatrixFile;
/* Begin `spFileMatrix'. */
assert( IS_SPARSE( Matrix ) );
/* Begin `spFileMatrix'. */
ASSERT_IS_SPARSE( Matrix );
/* Open file matrix file in write mode. */
/* Open file matrix file in write mode. */
if ((pMatrixFile = fopen(File, "w")) == NULL)
return 0;
/* Output header. */
/* Output header. */
Size = Matrix->Size;
if (Header)
{
if (Matrix->Factored && Data)
{
Err = fprintf(pMatrixFile,
"Warning : The following matrix is "
"factored in to LU form.\n");
if (Err < 0)
return 0;
{ 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;
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;
(Matrix->Complex ? "complex" : "real"));
if (Err < 0) return 0;
}
if (Size == 0) return 1;
/* Output matrix. */
if (!Data)
{
for (I = 1; I <= Size; I++)
{
pElement = Matrix->FirstInCol[I];
/* Output matrix. */
if (NOT Data)
{ for (I = 1; I <= Size; I++)
{ pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{
if (Reordered)
{
Row = pElement->Row;
{ if (Reordered)
{ Row = pElement->Row;
Col = I;
}
else
{
Row = Matrix->IntToExtRowMap[pElement->Row];
{ 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. */
/* Output terminator, a line of zeros. */
if (Header)
if (fprintf(pMatrixFile, "0\t0\n") < 0) return 0;
}
if (Data && Matrix->Complex)
{
for (I = 1; I <= Size; I++)
{
pElement = Matrix->FirstInCol[I];
#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;
{ if (Reordered)
{ Row = pElement->Row;
Col = I;
}
else
{
Row = Matrix->IntToExtRowMap[pElement->Row];
{ Row = Matrix->IntToExtRowMap[pElement->Row];
Col = Matrix->IntToExtColMap[I];
}
Err = fprintf
( pMatrixFile,"%d\t%d\t%-.15g\t%-.15g\n",
Row, Col, pElement->Real, pElement->Imag
);
( pMatrixFile,"%d\t%d\t%-.15g\t%-.15g\n",
Row, Col, (double)pElement->Real, (double)pElement->Imag
);
if (Err < 0) return 0;
pElement = pElement->NextInCol;
}
}
/* Output terminator, a line of zeros. */
/* Output terminator, a line of zeros. */
if (Header)
if (fprintf(pMatrixFile,"0\t0\t0.0\t0.0\n") < 0) return 0;
}
#endif /* spCOMPLEX */
if (Data && !Matrix->Complex)
{
for (I = 1; I <= Size; I++)
{
pElement = Matrix->FirstInCol[I];
#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];
{ Row = Matrix->IntToExtRowMap[pElement->Row];
Col = Matrix->IntToExtColMap[I];
Err = fprintf
( pMatrixFile,"%d\t%d\t%-.15g\n",
Row, Col, pElement->Real
);
( pMatrixFile,"%d\t%d\t%-.15g\n",
Row, Col, (double)pElement->Real
);
if (Err < 0) return 0;
pElement = pElement->NextInCol;
}
}
/* Output terminator, a line of zeros. */
/* Output terminator, a line of zeros. */
if (Header)
if (fprintf(pMatrixFile,"0\t0\t0.0\n") < 0) return 0;
}
#endif /* REAL */
/* Close file. */
/* Close file. */
if (fclose(pMatrixFile) < 0) return 0;
return 1;
}
@ -564,79 +518,117 @@ spFileMatrix(MatrixPtr Matrix, char *File, char *Label, int Reordered,
/*
* 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:
* \return
* One is returned if routine was successful, otherwise zero is returned.
* The calling function can query errno (the system global error variable)
* The calling function can query \a errno (the system global error variable)
* as to the reason why this routine failed.
*
* >>> Arguments:
* Matrix <input> (char *)
* \param eMatrix
* Pointer to matrix.
* File <input> (char *)
* \param File
* Name of file into which matrix is to be written.
* RHS <input> (RealNumber [])
* Right-hand side vector, real portion
* iRHS <input> (RealNumber [])
* Right-hand side vector, imaginary portion.
*
* >>> Local variables:
* \param RHS
* Right-hand side vector. This is only the real portion if
* \a spSEPARATED_COMPLEX_VECTORS is true.
* \param iRHS
* Right-hand side vector, imaginary portion. Not necessary if matrix
* is real or if \a spSEPARATED_COMPLEX_VECTORS is set false.
* \a iRHS is a macro that replaces itself with `, iRHS' if the options
* \a spCOMPLEX and \a spSEPARATED_COMPLEX_VECTORS are set, otherwise
* it disappears without a trace.
*/
/* >>> Local variables:
* pMatrixFile (FILE *)
* File pointer to the matrix file.
* Size (int)
* The size of the matrix.
*
*/
int
spFileVector(MatrixPtr Matrix, char *File, RealVector RHS, RealVector iRHS)
spFileVector(
spMatrix eMatrix,
char *File,
spREAL RHS[]
#if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
, spREAL iRHS[]
#endif
)
{
int I, Size, Err;
FILE *pMatrixFile;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register int I, Size, Err;
FILE *pMatrixFile;
/* Begin `spFileVector'. */
assert( IS_SPARSE( Matrix ) && RHS != NULL);
/* Begin `spFileVector'. */
ASSERT_IS_SPARSE( Matrix );
vASSERT( RHS != NULL, "Vector missing" );
if (File) {
/* Open File in write mode. */
pMatrixFile = fopen(File,"w");
if (pMatrixFile == NULL)
return 0;
}
else
pMatrixFile=stdout;
/* Open File in append mode. */
if ((pMatrixFile = fopen(File,"a")) == NULL)
return 0;
/* Output vector. */
Size = Matrix->Size;
/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spCOMPLEX
if (Matrix->Complex)
{
for (I = 1; I <= Size; I++)
{
Err = fprintf
( pMatrixFile, "%-.15g\t%-.15g\n",
RHS[I], iRHS[I]
);
if (Err < 0) return 0;
}
#if spSEPARATED_COMPLEX_VECTORS
vASSERT( iRHS != NULL, "Imaginary vector missing" );
--RHS;
--iRHS;
#else
RHS -= 2;
#endif
}
else
#endif /* spCOMPLEX */
--RHS;
#endif /* NOT ARRAY_OFFSET */
/* Output vector. */
Size = Matrix->Size;
if (Size == 0) return 1;
#if spCOMPLEX
if (Matrix->Complex)
{
for (I = 1; I <= Size; I++)
{
if (fprintf(pMatrixFile, "%-.15g\n", RHS[I]) < 0)
#if spSEPARATED_COMPLEX_VECTORS
for (I = 1; I <= Size; I++)
{ Err = fprintf
( pMatrixFile, "%-.15g\t%-.15g\n",
(double)RHS[I], (double)iRHS[I]
);
if (Err < 0) return 0;
}
#else
for (I = 1; I <= Size; I++)
{ Err = fprintf
( pMatrixFile, "%-.15g\t%-.15g\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, "%-.15g\n", (double)RHS[I]) < 0)
return 0;
}
}
#endif /* REAL */
/* Close file. */
if (File)
if (fclose(pMatrixFile) < 0) return 0;
/* Close file. */
if (fclose(pMatrixFile) < 0) return 0;
return 1;
}
@ -648,26 +640,23 @@ spFileVector(MatrixPtr Matrix, char *File, RealVector RHS, RealVector iRHS)
/*
* OUTPUT STATISTICS TO FILE
*
/*!
* Writes useful information concerning the matrix to a file. Should be
* executed after the matrix is factored.
*
* >>> Returns:
* \return
* One is returned if routine was successful, otherwise zero is returned.
* The calling function can query errno (the system global error variable)
* The calling function can query \a errno (the system global error variable)
* as to the reason why this routine failed.
*
* >>> Arguments:
* Matrix <input> (char *)
* \param eMatrix
* Pointer to matrix.
* File <input> (char *)
* \param File
* Name of file into which matrix is to be written.
* Label <input> (char *)
* \param Label
* String that is transferred to file and is used as a label.
*
* >>> Local variables:
*/
/* >>> Local variables:
* Data (RealNumber)
* The value of the matrix element being output.
* LargestElement (RealNumber)
@ -685,24 +674,29 @@ spFileVector(MatrixPtr Matrix, char *File, RealVector RHS, RealVector iRHS)
*/
int
spFileStats(MatrixPtr Matrix, char *File, char *Label)
spFileStats(
spMatrix eMatrix,
char *File,
char *Label
)
{
int Size, I;
ElementPtr pElement;
int NumberOfElements;
RealNumber Data, LargestElement, SmallestElement;
FILE *pStatsFile;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register int Size, I;
register ElementPtr pElement;
int NumberOfElements;
RealNumber Data, LargestElement, SmallestElement;
FILE *pStatsFile;
/* Begin `spFileStats'. */
assert( IS_SPARSE( Matrix ) );
/* Begin `spFileStats'. */
ASSERT_IS_SPARSE( Matrix );
/* Open File in append mode. */
/* Open File in append mode. */
if ((pStatsFile = fopen(File, "a")) == NULL)
return 0;
/* Output statistics. */
/* Output statistics. */
Size = Matrix->Size;
if (!Matrix->Factored)
if (NOT Matrix->Factored)
fprintf(pStatsFile, "Matrix has not been factored.\n");
fprintf(pStatsFile, "||| Starting new matrix |||\n");
fprintf(pStatsFile, "%s\n", Label);
@ -711,22 +705,21 @@ spFileStats(MatrixPtr Matrix, char *File, char *Label)
else
fprintf(pStatsFile, "Matrix is real.\n");
fprintf(pStatsFile," Size = %d\n",Size);
if (Size == 0) return 1;
/* Search matrix. */
/* Search matrix. */
NumberOfElements = 0;
LargestElement = 0.0;
SmallestElement = LARGEST_REAL;
for (I = 1; I <= Size; I++)
{
pElement = Matrix->FirstInCol[I];
{ pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{
NumberOfElements++;
{ NumberOfElements++;
Data = ELEMENT_MAG(pElement);
if (Data > LargestElement)
LargestElement = Data;
if (Data < SmallestElement && Data != 0.0)
if (Data < SmallestElement AND Data != 0.0)
SmallestElement = Data;
pElement = pElement->NextInCol;
}
@ -734,7 +727,7 @@ spFileStats(MatrixPtr Matrix, char *File, char *Label)
SmallestElement = MIN( SmallestElement, LargestElement );
/* Output remaining statistics. */
/* Output remaining statistics. */
fprintf(pStatsFile, " Initial number of elements = %d\n",
NumberOfElements - Matrix->Fillins);
fprintf(pStatsFile,
@ -748,13 +741,13 @@ spFileStats(MatrixPtr Matrix, char *File, char *Label)
fprintf(pStatsFile, " Average number of elements per row = %f\n",
(double)NumberOfElements / (double)Size);
fprintf(pStatsFile," Density = %f%%\n",
(double)(100.0*NumberOfElements)/(double)(Size*Size));
(100.0*(double)NumberOfElements)/((double)Size*(double)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. */
/* Close file. */
(void)fclose(pStatsFile);
return 1;
}

378
src/maths/sparse/spsmp.c → src/maths/sparse/spSMP.c
View File

@ -5,11 +5,14 @@
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
* UC Berkeley
*
* This module contains routines that make Sparse1.3 a direct
* replacement for the SMP sparse matrix package in Spice3c1 or Spice3d1.
* Sparse1.3 is in general a faster and more robust package than SMP.
* This module contains routines that make Sparse1.4 a direct
* replacement for the SMP sparse matrix package in Spice3c1 and Spice3d1.
* Sparse1.4 is in general a faster and more robust package than SMP.
* These advantages become significant on large circuits.
*
* This module is provided for convience only. It has not been tested
* with the recent version of Spice3 and is not supported.
*
* >>> User accessible functions contained in this file:
* SMPaddElt
* SMPmakeElt
@ -29,8 +32,6 @@
* SMPprint
* SMPgetError
* SMPcProdDiag
* LoadGmin
* SMPfindElt
*/
/*
@ -45,10 +46,12 @@
* To be compatible with SPICE, the following Sparse compiler options
* (in spConfig.h) should be set as shown below:
*
* REAL YES
* EXPANDABLE YES
* TRANSLATE NO
* INITIALIZE NO or YES, YES for use with test prog.
* DIAGONAL_PIVOTING YES
* ARRAY_OFFSET YES
* MODIFIED_MARKOWITZ NO
* DELETE NO
* STRIP NO
@ -62,7 +65,10 @@
* STABILITY NO
* CONDITION NO
* PSEUDOCONDITION NO
* FORTRAN NO
* DEBUG YES
* spCOMPLEX 1
* spSEPARATED_COMPLEX_VECTORS 1
*
* spREAL double
*/
@ -70,18 +76,10 @@
/*
* Revision and copyright information.
*
* Copyright (c) 1985,86,87,88,89,90
* 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.
* Copyright (c) 1985-2003 by Kenneth S. Kundert
*/
/*
* IMPORTS
*
@ -92,73 +90,71 @@
* Spice3's matrix macro definitions.
*/
#include "ngspice/config.h"
#include <assert.h>
#include <stdio.h>
#include <math.h>
#include "ngspice/spmatrix.h"
#include "spdefs.h"
#include "ngspice/smpdefs.h"
#if defined (_MSC_VER)
extern double scalbn(double, int);
#define logb _logb
extern double logb(double);
#endif
#define NO 0
#define YES 1
static void LoadGmin(SMPmatrix *Matrix, double Gmin);
typedef spREAL RealNumber, *RealVector;
static void LoadGmin(char *Matrix, double Gmin);
/*
* SMPaddElt()
*/
int
SMPaddElt(SMPmatrix *Matrix, int Row, int Col, double Value)
SMPaddElt( Matrix, Row, Col, Value )
SMPmatrix *Matrix;
int Row, Col;
double Value;
{
*spGetElement( Matrix, Row, Col ) = Value;
return spError( Matrix );
*spGetElement( (char *)Matrix, Row, Col ) = Value;
return spErrorState( (char *)Matrix );
}
/*
* SMPmakeElt()
*/
double *
SMPmakeElt(SMPmatrix *Matrix, int Row, int Col)
SMPmakeElt( Matrix, Row, Col )
SMPmatrix *Matrix;
int Row, Col;
{
return spGetElement( Matrix, Row, Col );
return spGetElement( (char *)Matrix, Row, Col );
}
/*
* SMPcClear()
*/
void
SMPcClear(SMPmatrix *Matrix)
SMPcClear( Matrix )
SMPmatrix *Matrix;
{
spClear( Matrix );
spClear( (char *)Matrix );
}
/*
* SMPclear()
*/
void
SMPclear(SMPmatrix *Matrix)
SMPclear( Matrix )
SMPmatrix *Matrix;
{
spClear( Matrix );
spClear( (char *)Matrix );
}
#define NG_IGNORE(x) (void)x
/*
* SMPcLUfac()
*/
/*ARGSUSED*/
int
SMPcLUfac(SMPmatrix *Matrix, double PivTol)
SMPcLUfac( Matrix, PivTol )
SMPmatrix *Matrix;
double PivTol;
{
NG_IGNORE(PivTol);
spSetComplex( Matrix );
return spFactor( Matrix );
spSetComplex( (char *)Matrix );
return spFactor( (char *)Matrix );
}
/*
@ -166,93 +162,97 @@ SMPcLUfac(SMPmatrix *Matrix, double PivTol)
*/
/*ARGSUSED*/
int
SMPluFac(SMPmatrix *Matrix, double PivTol, double Gmin)
SMPluFac( Matrix, PivTol, Gmin )
SMPmatrix *Matrix;
double PivTol, Gmin;
{
NG_IGNORE(PivTol);
spSetReal( Matrix );
LoadGmin( Matrix, Gmin );
return spFactor( Matrix );
spSetReal( (char *)Matrix );
LoadGmin( (char *)Matrix, Gmin );
return spFactor( (char *)Matrix );
}
/*
* SMPcReorder()
*/
int
SMPcReorder(SMPmatrix *Matrix, double PivTol, double PivRel,
int *NumSwaps)
SMPcReorder( Matrix, PivTol, PivRel, NumSwaps )
SMPmatrix *Matrix;
double PivTol, PivRel;
int *NumSwaps;
{
*NumSwaps = 1;
spSetComplex( Matrix );
return spOrderAndFactor( Matrix, NULL,
PivRel, PivTol, YES );
*NumSwaps = 0;
spSetComplex( (char *)Matrix );
return spOrderAndFactor( (char *)Matrix, (spREAL*)NULL,
(spREAL)PivRel, (spREAL)PivTol, YES );
}
/*
* SMPreorder()
*/
int
SMPreorder(SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin)
SMPreorder( Matrix, PivTol, PivRel, Gmin )
SMPmatrix *Matrix;
double PivTol, PivRel, Gmin;
{
spSetReal( Matrix );
LoadGmin( Matrix, Gmin );
return spOrderAndFactor( Matrix, NULL,
PivRel, PivTol, YES );
spSetComplex( (char *)Matrix );
LoadGmin( (char *)Matrix, Gmin );
return spOrderAndFactor( (char *)Matrix, (spREAL*)NULL,
(spREAL)PivRel, (spREAL)PivTol, YES );
}
/*
* SMPcaSolve()
*/
void
SMPcaSolve(SMPmatrix *Matrix, double RHS[], double iRHS[],
double Spare[], double iSpare[])
SMPcaSolve( Matrix, RHS, iRHS, Spare, iSpare)
SMPmatrix *Matrix;
double RHS[], iRHS[], Spare[], iSpare[];
{
NG_IGNORE(iSpare);
NG_IGNORE(Spare);
spSolveTransposed( Matrix, RHS, RHS, iRHS, iRHS );
spSolveTransposed( (char *)Matrix, RHS, RHS, iRHS, iRHS );
}
/*
* SMPcSolve()
*/
void
SMPcSolve(SMPmatrix *Matrix, double RHS[], double iRHS[],
double Spare[], double iSpare[])
SMPcSolve( Matrix, RHS, iRHS, Spare, iSpare)
SMPmatrix *Matrix;
double RHS[], iRHS[], Spare[], iSpare[];
{
NG_IGNORE(iSpare);
NG_IGNORE(Spare);
spSolve( Matrix, RHS, RHS, iRHS, iRHS );
spSolve( (char *)Matrix, RHS, RHS, iRHS, iRHS );
}
/*
* SMPsolve()
*/
void
SMPsolve(SMPmatrix *Matrix, double RHS[], double Spare[])
SMPsolve( Matrix, RHS, Spare )
SMPmatrix *Matrix;
double RHS[], Spare[];
{
NG_IGNORE(Spare);
spSolve( Matrix, RHS, RHS, NULL, NULL );
spSolve( (char *)Matrix, RHS, RHS, (spREAL*)NULL, (spREAL*)NULL );
}
/*
* SMPmatSize()
*/
int
SMPmatSize(SMPmatrix *Matrix)
SMPmatSize( Matrix )
SMPmatrix *Matrix;
{
return spGetSize( Matrix, 1 );
return spGetSize( (char *)Matrix, 1 );
}
/*
* SMPnewMatrix()
*/
int
SMPnewMatrix(SMPmatrix **pMatrix, int size)
SMPnewMatrix( pMatrix, dummy )
SMPmatrix **pMatrix;
int dummy;
{
int Error;
*pMatrix = spCreate( size, 1, &Error );
int Error;
*pMatrix = (SMPmatrix *)spCreate( 0, 1, &Error );
return Error;
}
@ -260,63 +260,67 @@ SMPnewMatrix(SMPmatrix **pMatrix, int size)
* SMPdestroy()
*/
void
SMPdestroy(SMPmatrix *Matrix)
SMPdestroy( Matrix )
SMPmatrix *Matrix;
{
spDestroy( Matrix );
spDestroy( (char *)Matrix );
}
/*
* SMPpreOrder()
*/
int
SMPpreOrder(SMPmatrix *Matrix)
SMPpreOrder( Matrix )
SMPmatrix *Matrix;
{
spMNA_Preorder( Matrix );
return spError( Matrix );
spMNA_Preorder( (char *)Matrix );
return spErrorState( (char *)Matrix );
}
/*
* SMPprint()
* SMPprintRHS()
*/
void
SMPprintRHS(SMPmatrix *Matrix, char *Filename, RealVector RHS, RealVector iRHS)
{
spFileVector( Matrix, Filename, RHS, iRHS );
spFileVector( Matrix, Filename, RHS, iRHS );
}
/*
* SMPprint()
*/
/*ARGSUSED*/
void
SMPprint(SMPmatrix *Matrix, char *Filename)
SMPprint( Matrix, File )
SMPmatrix *Matrix;
char *File;
{
if (Filename)
spFileMatrix(Matrix, Filename, "Circuit Matrix", 0, 1, 1 );
else
spPrint( Matrix, 0, 1, 1 );
spPrint( (char *)Matrix, 0, 1, 1 );
}
/*
* SMPgetError()
*/
void
SMPgetError(SMPmatrix *Matrix, int *Col, int *Row)
SMPgetError( Matrix, Col, Row)
SMPmatrix *Matrix;
int *Row, *Col;
{
spWhereSingular( Matrix, Row, Col );
spWhereSingular( (char *)Matrix, Row, Col );
}
/*
* SMPcProdDiag()
* note: obsolete for Spice3d2 and later
*/
int
SMPcProdDiag(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
SMPcProdDiag( Matrix, pMantissa, pExponent)
SMPmatrix *Matrix;
SPcomplex *pMantissa;
int *pExponent;
{
spDeterminant( Matrix, pExponent, &(pMantissa->real),
spDeterminant( (char *)Matrix, pExponent, &(pMantissa->real),
&(pMantissa->imag) );
return spError( Matrix );
return spErrorState( (char *)Matrix );
}
/*
@ -325,8 +329,8 @@ SMPcProdDiag(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
int
SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
{
double re, im, x, y, z;
int p;
double re, im, x, y, z;
int p;
spDeterminant( Matrix, &p, &re, &im);
@ -347,7 +351,7 @@ SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
y -= x;
/* ASSERT
* x = integral part of exponent, y = fraction part of exponent
* x = integral part of exponent, y = fraction part of exponent
*/
/* Fold in the fractional part */
@ -363,51 +367,45 @@ SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
/* Re-normalize (re or im may be > 2.0 or both < 1.0 */
if (re != 0.0) {
y = logb(re);
if (im != 0.0)
z = logb(im);
else
z = 0;
y = logb(re);
if (im != 0.0)
z = logb(im);
else
z = 0;
} else if (im != 0.0) {
z = logb(im);
y = 0;
z = logb(im);
y = 0;
} else {
/* Singular */
/*printf("10 -> singular\n");*/
y = 0;
z = 0;
/* Singular */
/*printf("10 -> singular\n");*/
y = 0;
z = 0;
}
#ifdef debug_print
printf(" ** renormalize changes = %g,%g\n", y, z);
#endif
if (y < z)
y = z;
y = z;
*pExponent = (int)(x + y);
x = scalbn(re, (int) -y);
z = scalbn(im, (int) -y);
#ifdef debug_print
printf(" ** values are: re %g, im %g, y %g, re' %g, im' %g\n",
re, im, y, x, z);
re, im, y, x, z);
#endif
pMantissa->real = scalbn(re, (int) -y);
pMantissa->imag = scalbn(im, (int) -y);
#ifdef debug_print
printf("Determinant 10->2: (%20g,%20g)^%d\n", pMantissa->real,
pMantissa->imag, *pExponent);
pMantissa->imag, *pExponent);
#endif
return spError( Matrix );
return spErrorState( Matrix );
}
/*
* The following routines need internal knowledge of the Sparse data
* structures.
*/
/*
* LOAD GMIN
*
@ -418,47 +416,68 @@ SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
* use of this routine is not recommended. It is included here simply
* for compatibility with Spice3.
*/
static void
LoadGmin(SMPmatrix *Matrix, double Gmin)
#include "spDefs.h"
void
LoadGmin( eMatrix, Gmin )
char *eMatrix;
register double Gmin;
{
int I;
ArrayOfElementPtrs Diag;
ElementPtr diag;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register int I;
register ArrayOfElementPtrs Diag;
/* Begin `LoadGmin'. */
assert( IS_SPARSE( Matrix ) );
/* Begin `spLoadGmin'. */
ASSERT_IS_SPARSE( Matrix );
if (Gmin != 0.0) {
Diag = Matrix->Diag;
for (I = Matrix->Size; I > 0; I--) {
if ((diag = Diag[I]) != NULL)
diag->Real += Gmin;
}
}
Diag = Matrix->Diag;
for (I = Matrix->Size; I > 0; I--)
Diag[I]->Real += Gmin;
return;
}
/*
* FIND ELEMENT
*
* This routine finds an element in the matrix by row and column number.
* If the element exists, a pointer to it is returned. If not, then NULL
* is returned unless the CreateIfMissing flag is TRUE, in which case a
* is returned unless the CreateIfMissing flag is true, in which case a
* pointer to the new element is returned.
*/
//SMPelement *
//SMPfindElt( Matrix, Row, Col, CreateIfMissing )
//
//SMPmatrix *Matrix;
//int Row, Col;
//int CreateIfMissing;
//{
////MatrixPtr Matrix = (MatrixPtr)eMatrix;
//spREAL *Element = (spREAL *)Matrix->FirstInCol[Col];
//
///* Begin `SMPfindElt'. */
// ASSERT_IS_SPARSE( Matrix );
// if (CreateIfMissing)
// { Element = spcCreateElement( Matrix, Row, Col,
// &Matrix->FirstInRow[Row],
// &Matrix->FirstInCol[Col], NO );
// }
// else Element = spcFindElement( Matrix, Row, Col );
// return (SMPelement *)Element;
//}
SMPelement *
SMPfindElt(SMPmatrix *Matrix, int Row, int Col, int CreateIfMissing)
{
ElementPtr Element;
/* Begin `SMPfindElt'. */
assert( IS_SPARSE( Matrix ) );
ASSERT_IS_SPARSE( Matrix );
Row = Matrix->ExtToIntRowMap[Row];
Col = Matrix->ExtToIntColMap[Col];
if (Col == -1)
@ -470,7 +489,6 @@ SMPfindElt(SMPmatrix *Matrix, int Row, int Col, int CreateIfMissing)
return Element;
}
/* XXX The following should probably be implemented in spUtils */
/*
* SMPcZeroCol()
@ -478,19 +496,19 @@ SMPfindElt(SMPmatrix *Matrix, int Row, int Col, int CreateIfMissing)
int
SMPcZeroCol(SMPmatrix *Matrix, int Col)
{
ElementPtr Element;
ElementPtr Element;
Col = Matrix->ExtToIntColMap[Col];
for (Element = Matrix->FirstInCol[Col];
Element != NULL;
Element = Element->NextInCol)
Element != NULL;
Element = Element->NextInCol)
{
Element->Real = 0.0;
Element->Imag = 0.0;
Element->Real = 0.0;
Element->Imag = 0.0;
}
return spError( Matrix );
return spErrorState( Matrix );
}
/*
@ -499,7 +517,7 @@ SMPcZeroCol(SMPmatrix *Matrix, int Col)
int
SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col)
{
ElementPtr Accum, Addend, *Prev;
ElementPtr Accum, Addend, *Prev;
Accum_Col = Matrix->ExtToIntColMap[Accum_Col];
Addend_Col = Matrix->ExtToIntColMap[Addend_Col];
@ -509,52 +527,29 @@ SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col)
Accum = *Prev;
while (Addend != NULL) {
while (Accum && Accum->Row < Addend->Row) {
Prev = &Accum->NextInCol;
Accum = *Prev;
}
if (!Accum || Accum->Row > Addend->Row) {
Accum = spcCreateElement(Matrix, Addend->Row, Accum_Col, Prev, 0);
}
Accum->Real += Addend->Real;
Accum->Imag += Addend->Imag;
Addend = Addend->NextInCol;
while (Accum && Accum->Row < Addend->Row) {
Prev = &Accum->NextInCol;
Accum = *Prev;
}
if (!Accum || Accum->Row > Addend->Row) {
Accum = spcCreateElement(Matrix, Addend->Row, Accum_Col, Prev, 0, 0);
}
Accum->Real += Addend->Real;
Accum->Imag += Addend->Imag;
Addend = Addend->NextInCol;
}
return spError( Matrix );
return spErrorState( Matrix );
}
/*
* SMPzeroRow()
* SMPmultiply()
*/
int
SMPzeroRow(SMPmatrix *Matrix, int Row)
void
SMPmultiply(SMPmatrix *Matrix, double *RHS, double *Solution, double *iRHS, double *iSolution)
{
ElementPtr Element;
Row = Matrix->ExtToIntColMap[Row];
if (Matrix->RowsLinked == NO)
spcLinkRows(Matrix);
if (Matrix->PreviousMatrixWasComplex || Matrix->Complex) {
for (Element = Matrix->FirstInRow[Row];
Element != NULL;
Element = Element->NextInRow)
{
Element->Real = 0.0;
Element->Imag = 0.0;
}
} else {
for (Element = Matrix->FirstInRow[Row];
Element != NULL;
Element = Element->NextInRow)
{
Element->Real = 0.0;
}
}
return spError( Matrix );
spMultiply(Matrix, RHS, Solution, iRHS, iSolution);
}
/*
@ -565,12 +560,3 @@ SMPconstMult(SMPmatrix *Matrix, double constant)
{
spConstMult(Matrix, constant);
}
/*
* SMPmultiply()
*/
void
SMPmultiply(SMPmatrix *Matrix, double *RHS, double *Solution, double *iRHS, double *iSolution)
{
spMultiply(Matrix, RHS, Solution, iRHS, iSolution);
}

502
src/maths/sparse/spsolve.c → src/maths/sparse/spSolve.c
View File

@ -4,11 +4,18 @@
* Author: Advising professor:
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
* UC Berkeley
*
*/
/*! \file
* This file contains the forward and backward substitution routines for
* the sparse matrix routines.
*
* >>> User accessible functions contained in this file:
* Objects that begin with the \a spc prefix are considered private
* and should not be used.
*
* \author
* Kenneth S. Kundert <kundert@users.sourceforge.net>
*/
/* >>> User accessible functions contained in this file:
* spSolve
* spSolveTransposed
*
@ -21,19 +28,12 @@
/*
* Revision and copyright information.
*
* Copyright (c) 1985,86,87,88,89,90
* 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.
* Copyright (c) 1985-2003
* by Kenneth S. Kundert
*/
/*
* IMPORTS
*
@ -45,12 +45,12 @@
* spDefs.h
* Matrix type and macro definitions for the sparse matrix routines.
*/
#include <assert.h>
#define spINSIDE_SPARSE
#include "spconfig.h"
#include <stdio.h>
#include "spConfig.h"
#include "ngspice/spmatrix.h"
#include "spdefs.h"
#include "spDefs.h"
@ -59,10 +59,16 @@
* Function declarations
*/
#if spSEPARATED_COMPLEX_VECTORS
static void SolveComplexMatrix( MatrixPtr,
RealVector, RealVector, RealVector, RealVector );
static void SolveComplexTransposedMatrix( MatrixPtr,
RealVector, RealVector, RealVector, RealVector );
#else
static void SolveComplexMatrix( MatrixPtr, RealVector, RealVector );
static void SolveComplexTransposedMatrix( MatrixPtr,
RealVector, RealVector );
#endif
@ -70,35 +76,35 @@ static void SolveComplexTransposedMatrix( MatrixPtr,
/*
* SOLVE MATRIX EQUATION
*
/*!
* Performs forward elimination and back substitution to find the
* unknown vector from the RHS vector and factored matrix. This
* unknown vector from the \a 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
* triangular matrix and that the diagonal of the upper triangular
* 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 <input> (char *)
* \param eMatrix
* Pointer to matrix.
* RHS <input> (RealVector)
* RHS is the input data array, the right hand side. This data is
* \param RHS
* \a RHS is the input data array, the right hand side. This data is
* undisturbed and may be reused for other solves.
* Solution <output> (RealVector)
* Solution is the output data array. This routine is constructed such that
* RHS and Solution can be the same array.
* iRHS <input> (RealVector)
* iRHS is the imaginary portion of the input data array, the right
* \param Solution
* \a Solution is the output data array. This routine is constructed
* such that \a RHS and \a Solution can be the same array.
* \param iRHS
* \a 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.
* iSolution <output> (RealVector)
* iSolution is the imaginary portion of the output data array. This
* routine is constructed such that iRHS and iSolution can be
* the same array.
*
* >>> Local variables:
* This argument is only necessary if matrix is complex and if
* \a spSEPARATED_COMPLEX_VECTOR is set true.
* \param iSolution
* \a iSolution is the imaginary portion of the output data array. This
* routine is constructed such that \a iRHS and \a iSolution can be
* the same array. This argument is only necessary if matrix is complex
* and if \a spSEPARATED_COMPLEX_VECTOR is set true.
*/
/* >>> Local variables:
* Intermediate (RealVector)
* Temporary storage for use in forward elimination and backward
* substitution. Commonly referred to as c, when the LU factorization
@ -119,75 +125,96 @@ static void SolveComplexTransposedMatrix( MatrixPtr,
* Size of matrix. Made local to reduce indirection.
* Temp (RealNumber)
* Temporary storage for entries in arrays.
*
* >>> 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
spSolve(MatrixPtr Matrix, RealVector RHS, RealVector Solution,
RealVector iRHS, RealVector iSolution)
spSolve(
spMatrix eMatrix,
spREAL RHS[],
spREAL Solution[]
# if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
, spREAL iRHS[]
, spREAL iSolution[]
# endif
)
{
ElementPtr pElement;
RealVector Intermediate;
RealNumber Temp;
int I, *pExtOrder, Size;
ElementPtr pPivot;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register ElementPtr pElement;
register RealVector Intermediate;
register RealNumber Temp;
register int I, *pExtOrder, Size;
ElementPtr pPivot;
void SolveComplexMatrix();
/* Begin `spSolve'. */
assert( IS_VALID(Matrix) && IS_FACTORED(Matrix) );
/* Begin `spSolve'. */
ASSERT_IS_SPARSE( Matrix );
ASSERT_NO_ERRORS( Matrix );
ASSERT_IS_FACTORED( Matrix );
#if spCOMPLEX
if (Matrix->Complex)
{
SolveComplexMatrix( Matrix, RHS, Solution, iRHS, iSolution );
{ SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS );
return;
}
#endif
#if REAL
Intermediate = Matrix->Intermediate;
Size = Matrix->Size;
/* Initialize Intermediate vector. */
/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
--RHS;
--Solution;
#endif
/* Initialize Intermediate vector. */
pExtOrder = &Matrix->IntToExtRowMap[Size];
for (I = Size; I > 0; I--)
Intermediate[I] = RHS[*(pExtOrder--)];
/* Forward elimination. Solves Lc = b.*/
/* Forward elimination. Solves Lc = b.*/
for (I = 1; I <= Size; I++)
{
/* This step of the elimination is skipped if Temp equals zero. */
{
/* This step of the elimination is skipped if Temp equals zero. */
if ((Temp = Intermediate[I]) != 0.0)
{
pPivot = Matrix->Diag[I];
{ pPivot = Matrix->Diag[I];
Intermediate[I] = (Temp *= pPivot->Real);
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
Intermediate[pElement->Row] -= Temp * pElement->Real;
{ Intermediate[pElement->Row] -= Temp * pElement->Real;
pElement = pElement->NextInCol;
}
}
}
/* Backward Substitution. Solves Ux = c.*/
/* Backward Substitution. Solves Ux = c.*/
for (I = Size; I > 0; I--)
{
Temp = Intermediate[I];
{ Temp = Intermediate[I];
pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{
Temp -= pElement->Real * Intermediate[pElement->Col];
{ Temp -= pElement->Real * Intermediate[pElement->Col];
pElement = pElement->NextInRow;
}
Intermediate[I] = Temp;
}
/* Unscramble Intermediate vector while placing data in to Solution vector. */
/* Unscramble Intermediate vector while placing data in to Solution vector. */
pExtOrder = &Matrix->IntToExtColMap[Size];
for (I = Size; I > 0; I--)
Solution[*(pExtOrder--)] = Intermediate[I];
return;
#endif /* REAL */
}
@ -200,39 +227,37 @@ spSolve(MatrixPtr Matrix, RealVector RHS, RealVector Solution,
/*
* SOLVE COMPLEX MATRIX EQUATION
*
#if spCOMPLEX
/*!
* 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
* triangular matrix and that the diagonal of the upper triangular
* 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 <input> (char *)
* \param Matrix
* Pointer to matrix.
* RHS <input> (RealVector)
* \param RHS
* RHS is the real portion of the input data array, the right hand
* side. This data is undisturbed and may be reused for other solves.
* Solution <output> (RealVector)
* \param Solution
* Solution is the real portion of the output data array. This routine
* is constructed such that RHS and Solution can be the same
* array.
* iRHS <input> (RealVector)
* \param iRHS
* 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, there is no need to
* If spSEPARATED_COMPLEX_VECTOR is set false, there is no need to
* supply this array.
* iSolution <output> (RealVector)
* \param iSolution
* 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, there is no
* the same array. If spSEPARATED_COMPLEX_VECTOR is set false, there is no
* need to supply this array.
*
* >>> Local variables:
*/
/* >>> Local variables:
* Intermediate (ComplexVector)
* Temporary storage for use in forward elimination and backward
* substitution. Commonly referred to as c, when the LU factorization
@ -256,44 +281,68 @@ spSolve(MatrixPtr Matrix, RealVector RHS, RealVector Solution,
*/
static void
SolveComplexMatrix( MatrixPtr Matrix, RealVector RHS, RealVector Solution , RealVector iRHS, RealVector iSolution )
SolveComplexMatrix(
MatrixPtr Matrix,
RealVector RHS,
RealVector Solution
# if spSEPARATED_COMPLEX_VECTORS
, RealVector iRHS
, RealVector iSolution
# endif
)
{
ElementPtr pElement;
ComplexVector Intermediate;
int I, *pExtOrder, Size;
ElementPtr pPivot;
ComplexNumber Temp;
register ElementPtr pElement;
register ComplexVector Intermediate;
register int I, *pExtOrder, Size;
ElementPtr pPivot;
#if NOT spSEPARATED_COMPLEX_VECTORS
register ComplexVector ExtVector;
#endif
ComplexNumber Temp;
/* Begin `SolveComplexMatrix'. */
/* Begin `SolveComplexMatrix'. */
Size = Matrix->Size;
Intermediate = (ComplexVector)Matrix->Intermediate;
/* Initialize Intermediate vector. */
/* 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. */
pExtOrder = &Matrix->IntToExtRowMap[Size];
#if spSEPARATED_COMPLEX_VECTORS
for (I = Size; I > 0; I--)
{
Intermediate[I].Real = RHS[*(pExtOrder)];
{ Intermediate[I].Real = RHS[*(pExtOrder)];
Intermediate[I].Imag = iRHS[*(pExtOrder--)];
}
#else
ExtVector = (ComplexVector)RHS;
for (I = Size; I > 0; I--)
Intermediate[I] = ExtVector[*(pExtOrder--)];
#endif
/* Forward substitution. Solves Lc = b.*/
/* Forward substitution. Solves Lc = b.*/
for (I = 1; I <= Size; I++)
{
Temp = Intermediate[I];
{ Temp = Intermediate[I];
/* This step of the substitution is skipped if Temp equals zero. */
if ((Temp.Real != 0.0) || (Temp.Imag != 0.0))
{
pPivot = Matrix->Diag[I];
/* Cmplx expr: Temp *= (1.0 / Pivot). */
/* This step of the substitution is skipped if Temp equals zero. */
if ((Temp.Real != 0.0) OR (Temp.Imag != 0.0))
{ pPivot = Matrix->Diag[I];
/* Cmplx expr: Temp *= (1.0 / Pivot). */
CMPLX_MULT_ASSIGN(Temp, *pPivot);
Intermediate[I] = Temp;
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
/* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */
/* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */
CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Row],
Temp, *pElement);
pElement = pElement->NextInCol;
@ -301,32 +350,37 @@ SolveComplexMatrix( MatrixPtr Matrix, RealVector RHS, RealVector Solution , Real
}
}
/* Backward Substitution. Solves Ux = c.*/
/* Backward Substitution. Solves Ux = c.*/
for (I = Size; I > 0; I--)
{
Temp = Intermediate[I];
{ Temp = Intermediate[I];
pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{
/* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */
/* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */
CMPLX_MULT_SUBT_ASSIGN(Temp, *pElement,Intermediate[pElement->Col]);
pElement = pElement->NextInRow;
}
Intermediate[I] = Temp;
}
/* Unscramble Intermediate vector while placing data in to Solution vector. */
/* Unscramble Intermediate vector while placing data in to Solution vector. */
pExtOrder = &Matrix->IntToExtColMap[Size];
#if spSEPARATED_COMPLEX_VECTORS
for (I = Size; I > 0; I--)
{
Solution[*(pExtOrder)] = Intermediate[I].Real;
{ Solution[*(pExtOrder)] = Intermediate[I].Real;
iSolution[*(pExtOrder--)] = Intermediate[I].Imag;
}
#else
ExtVector = (ComplexVector)Solution;
for (I = Size; I > 0; I--)
ExtVector[*(pExtOrder--)] = Intermediate[I];
#endif
return;
}
#endif /* spCOMPLEX */
@ -342,38 +396,35 @@ SolveComplexMatrix( MatrixPtr Matrix, RealVector RHS, RealVector Solution , Real
#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.
* matrix and that the diagonal of the untransposed upper
* triangular matrix consists of ones.
*
* >>> Arguments:
* Matrix <input> (char *)
* \param eMatrix
* Pointer to matrix.
* RHS <input> (RealVector)
* RHS is the input data array, the right hand side. This data is
* \param RHS
* \a RHS is the input data array, the right hand side. This data is
* undisturbed and may be reused for other solves.
* Solution <output> (RealVector)
* Solution is the output data array. This routine is constructed such that
* RHS and Solution can be the same array.
* iRHS <input> (RealVector)
* iRHS is the imaginary portion of the input data array, the right
* \param Solution
* \a Solution is the output data array. This routine is constructed
* such that \a RHS and \a Solution can be the same array.
* \param iRHS
* \a 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 <output> (RealVector)
* 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
* If \a spSEPARATED_COMPLEX_VECTOR is set false, or if matrix is real,
* there is no need to supply this array.
* \param iSolution
* \a iSolution is the imaginary portion of the output data array. This
* routine is constructed such that \a iRHS and \a iSolution can be
* the same array. If \a spSEPARATED_COMPLEX_VECTOR is set false, or if
* matrix is real, there is no need to supply this array.
*
* >>> Local variables:
*/
/* >>> Local variables:
* Intermediate (RealVector)
* Temporary storage for use in forward elimination and backward
* substitution. Commonly referred to as c, when the LU factorization
@ -399,70 +450,84 @@ SolveComplexMatrix( MatrixPtr Matrix, RealVector RHS, RealVector Solution , Real
/*VARARGS3*/
void
spSolveTransposed(MatrixPtr Matrix, RealVector RHS, RealVector Solution,
RealVector iRHS, RealVector iSolution)
spSolveTransposed(
spMatrix eMatrix,
spREAL RHS[],
spREAL Solution[]
# if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
, spREAL iRHS[]
, spREAL iSolution[]
# endif
)
{
ElementPtr pElement;
RealVector Intermediate;
int I, *pExtOrder, Size;
ElementPtr pPivot;
RealNumber Temp;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register ElementPtr pElement;
register RealVector Intermediate;
register int I, *pExtOrder, Size;
ElementPtr pPivot;
RealNumber Temp;
void SolveComplexTransposedMatrix();
/* Begin `spSolveTransposed'. */
assert( IS_VALID(Matrix) && IS_FACTORED(Matrix) );
/* Begin `spSolveTransposed'. */
ASSERT_IS_SPARSE( Matrix );
ASSERT_NO_ERRORS( Matrix );
ASSERT_IS_FACTORED( Matrix );
#if spCOMPLEX
if (Matrix->Complex)
{
SolveComplexTransposedMatrix( Matrix, RHS, Solution , iRHS, iSolution );
{ SolveComplexTransposedMatrix( Matrix, RHS, Solution IMAG_VECTORS );
return;
}
#endif
#if REAL
Size = Matrix->Size;
Intermediate = Matrix->Intermediate;
/* Initialize Intermediate vector. */
/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
--RHS;
--Solution;
#endif
/* Initialize Intermediate vector. */
pExtOrder = &Matrix->IntToExtColMap[Size];
for (I = Size; I > 0; I--)
Intermediate[I] = RHS[*(pExtOrder--)];
/* Forward elimination. */
/* Forward elimination. */
for (I = 1; I <= Size; I++)
{
/* This step of the elimination is skipped if Temp equals zero. */
{
/* This step of the elimination is skipped if Temp equals zero. */
if ((Temp = Intermediate[I]) != 0.0)
{
pElement = Matrix->Diag[I]->NextInRow;
{ pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{
Intermediate[pElement->Col] -= Temp * pElement->Real;
{ Intermediate[pElement->Col] -= Temp * pElement->Real;
pElement = pElement->NextInRow;
}
}
}
/* Backward Substitution. */
/* Backward Substitution. */
for (I = Size; I > 0; I--)
{
pPivot = Matrix->Diag[I];
{ pPivot = Matrix->Diag[I];
Temp = Intermediate[I];
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
Temp -= pElement->Real * Intermediate[pElement->Row];
{ Temp -= pElement->Real * Intermediate[pElement->Row];
pElement = pElement->NextInCol;
}
Intermediate[I] = Temp * pPivot->Real;
}
/* Unscramble Intermediate vector while placing data in to
Solution vector. */
/* Unscramble Intermediate vector while placing data in to Solution vector. */
pExtOrder = &Matrix->IntToExtRowMap[Size];
for (I = Size; I > 0; I--)
Solution[*(pExtOrder--)] = Intermediate[I];
return;
#endif /* REAL */
}
#endif /* TRANSPOSE */
@ -475,38 +540,40 @@ spSolveTransposed(MatrixPtr Matrix, RealVector RHS, RealVector Solution,
#if TRANSPOSE
/*
* SOLVE COMPLEX TRANSPOSED MATRIX EQUATION
*
#if TRANSPOSE AND spCOMPLEX
/*!
* 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.
* matrix and that the diagonal of the untransposed upper
* triangular matrix consists of ones.
*
* >>> Arguments:
* Matrix <input> (char *)
* \param Matrix
* Pointer to matrix.
* RHS <input> (RealVector)
* RHS is the input data array, the right hand
* \param RHS
* \a RHS is the input data array, the right hand
* side. This data is undisturbed and may be reused for other solves.
* This vector is only the real portion if the matrix is complex.
* Solution <output> (RealVector)
* Solution is the real portion of the output data array. This routine
* is constructed such that RHS and Solution can be the same array.
* This vector is only the real portion if the matrix is complex.
* iRHS <input> (RealVector)
* iRHS is the imaginary portion of the input data array, the right
* This vector is only the real portion if the matrix is complex and
* \a spSEPARATED_COMPLEX_VECTORS is set true.
* \param Solution
* \a Solution is the real portion of the output data array. This routine
* is constructed such that \a RHS and \a Solution can be the same array.
* This vector is only the real portion if the matrix is complex and
* \a spSEPARATED_COMPLEX_VECTORS is set true.
* \param iRHS
* \a 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.
* iSolution <output> (RealVector)
* iSolution is the imaginary portion of the output data array. This
* routine is constructed such that iRHS and iSolution can be
* the same array.
*
* >>> Local variables:
* If either \a spCOMPLEX or \a spSEPARATED_COMPLEX_VECTOR is set false,
* there is no need to supply this array.
* \param iSolution
* \a iSolution is the imaginary portion of the output data array. This
* routine is constructed such that \a iRHS and \a iSolution can be
* the same array. If \a spCOMPLEX or \a spSEPARATED_COMPLEX_VECTOR is set
* false, there is no need to supply this array.
*/
/* >>> Local variables:
* Intermediate (ComplexVector)
* Temporary storage for use in forward elimination and backward
* substitution. Commonly referred to as c, when the LU factorization
@ -530,40 +597,64 @@ spSolveTransposed(MatrixPtr Matrix, RealVector RHS, RealVector Solution,
*/
static void
SolveComplexTransposedMatrix(MatrixPtr Matrix, RealVector RHS, RealVector Solution , RealVector iRHS, RealVector iSolution )
SolveComplexTransposedMatrix(
MatrixPtr Matrix,
RealVector RHS,
RealVector Solution
# if spSEPARATED_COMPLEX_VECTORS
, RealVector iRHS
, RealVector iSolution
# endif
)
{
ElementPtr pElement;
ComplexVector Intermediate;
int I, *pExtOrder, Size;
ElementPtr pPivot;
ComplexNumber Temp;
register ElementPtr pElement;
register ComplexVector Intermediate;
register int I, *pExtOrder, Size;
#if NOT spSEPARATED_COMPLEX_VECTORS
register ComplexVector ExtVector;
#endif
ElementPtr pPivot;
ComplexNumber Temp;
/* Begin `SolveComplexTransposedMatrix'. */
/* Begin `SolveComplexTransposedMatrix'. */
Size = Matrix->Size;
Intermediate = (ComplexVector)Matrix->Intermediate;
/* Initialize Intermediate vector. */
/* 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. */
pExtOrder = &Matrix->IntToExtColMap[Size];
#if spSEPARATED_COMPLEX_VECTORS
for (I = Size; I > 0; I--)
{
Intermediate[I].Real = RHS[*(pExtOrder)];
{ Intermediate[I].Real = RHS[*(pExtOrder)];
Intermediate[I].Imag = iRHS[*(pExtOrder--)];
}
#else
ExtVector = (ComplexVector)RHS;
for (I = Size; I > 0; I--)
Intermediate[I] = ExtVector[*(pExtOrder--)];
#endif
/* Forward elimination. */
/* Forward elimination. */
for (I = 1; I <= Size; I++)
{
Temp = Intermediate[I];
{ Temp = Intermediate[I];
/* This step of the elimination is skipped if Temp equals zero. */
if ((Temp.Real != 0.0) || (Temp.Imag != 0.0))
{
pElement = Matrix->Diag[I]->NextInRow;
/* This step of the elimination is skipped if Temp equals zero. */
if ((Temp.Real != 0.0) OR (Temp.Imag != 0.0))
{ pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{
/* Cmplx expr: Intermediate[Element->Col] -= Temp * *Element. */
/* Cmplx expr: Intermediate[Element->Col] -= Temp * *Element. */
CMPLX_MULT_SUBT_ASSIGN( Intermediate[pElement->Col],
Temp, *pElement);
pElement = pElement->NextInRow;
@ -571,34 +662,37 @@ SolveComplexTransposedMatrix(MatrixPtr Matrix, RealVector RHS, RealVector Soluti
}
}
/* Backward Substitution. */
/* Backward Substitution. */
for (I = Size; I > 0; I--)
{
pPivot = Matrix->Diag[I];
{ pPivot = Matrix->Diag[I];
Temp = Intermediate[I];
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
/* Cmplx expr: Temp -= Intermediate[Element->Row] * *Element. */
/* Cmplx expr: Temp -= Intermediate[Element->Row] * *Element. */
CMPLX_MULT_SUBT_ASSIGN(Temp,Intermediate[pElement->Row],*pElement);
pElement = pElement->NextInCol;
}
/* Cmplx expr: Intermediate = Temp * (1.0 / *pPivot). */
/* Cmplx expr: Intermediate = Temp * (1.0 / *pPivot). */
CMPLX_MULT(Intermediate[I], Temp, *pPivot);
}
/* Unscramble Intermediate vector while placing data in to
Solution vector. */
/* Unscramble Intermediate vector while placing data in to Solution vector. */
pExtOrder = &Matrix->IntToExtRowMap[Size];
#if spSEPARATED_COMPLEX_VECTORS
for (I = Size; I > 0; I--)
{
Solution[*(pExtOrder)] = Intermediate[I].Real;
{ Solution[*(pExtOrder)] = Intermediate[I].Real;
iSolution[*(pExtOrder--)] = Intermediate[I].Imag;
}
#else
ExtVector = (ComplexVector)Solution;
for (I = Size; I > 0; I--)
ExtVector[*(pExtOrder--)] = Intermediate[I];
#endif
return;
}
#endif /* TRANSPOSE */
#endif /* TRANSPOSE AND spCOMPLEX */

1601
src/maths/sparse/sputils.c → src/maths/sparse/spUtils.c
View File

File diff suppressed because it is too large Load Diff

885
src/maths/sparse/spalloc.c
View File

@ -1,885 +0,0 @@
/*
* MATRIX ALLOCATION MODULE
*
* Author: Advising professor:
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
* UC Berkeley
*
* This file contains the allocation and deallocation routines for the
* sparse matrix routines.
*
* >>> User accessible functions contained in this file:
* spCreate
* spDestroy
* spError
* spWhereSingular
* spGetSize
* spSetReal
* spSetComplex
* spFillinCount
* spElementCount
* spOriginalCount
*
* >>> Other functions contained in this file:
* spcGetElement
* InitializeElementBlocks
* spcGetFillin
* RecordAllocation
* AllocateBlockOfAllocationList
* EnlargeMatrix
* ExpandTranslationArrays
*/
/*
* Revision and copyright information.
*
* Copyright (c) 1985,86,87,88,89,90
* by Kenneth S. Kundert and the University of California.
*
* Permission to use, copy, modify, and distribute this software and
* its documentation for any purpose and without fee is hereby granted,
* provided that the copyright notices appear in all copies and
* supporting documentation and that the authors and the University of
* California are properly credited. The authors and the University of
* California make no representations as to the suitability of this
* software for any purpose. It is provided `as is', without express
* or implied warranty.
*/
/*
* IMPORTS
*
* >>> 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.
*/
#include <assert.h>
#include <stdlib.h>
#define spINSIDE_SPARSE
#include "spconfig.h"
#include "ngspice/spmatrix.h"
#include "spdefs.h"
/*
* Function declarations
*/
static void InitializeElementBlocks( MatrixPtr, int, int );
static void RecordAllocation( MatrixPtr, void *);
static void AllocateBlockOfAllocationList( MatrixPtr );
/*
* MATRIX ALLOCATION
*
* Allocates and initializes the data structures associated with a matrix.
*
* >>> Returned:
* A pointer to the matrix is returned cast into the form of a pointer to
* a character. 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 <input> (int)
* Size of matrix or estimate of size of matrix if matrix is EXPANDABLE.
* Complex <input> (int)
* 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.
* pError <output> (int *)
* Returns error flag, needed because function spError() will not work
* correctly if spCreate() returns NULL.
*
* >>> Local variables:
* AllocatedSize (int)
* The size of the matrix being allocated.
* Matrix (MatrixPtr)
* A pointer to the matrix frame being created.
*
* >>> Possible errors:
* spNO_MEMORY
* spPANIC
* Error is cleared in this routine.
*/
MatrixPtr
spCreate(int Size, int Complex, int *pError)
{
unsigned SizePlusOne;
MatrixPtr Matrix;
int I;
int AllocatedSize;
/* Begin `spCreate'. */
/* Clear error flag. */
*pError = spOKAY;
/* Test for valid size. */
#if EXPANDABLE
if (Size < 0) {
*pError = spPANIC;
return NULL;
}
#else
if (Size <= 0) {
*pError = spPANIC;
return NULL;
}
#endif
#if 0 /* pn: skipped for cider */
/* Test for valid type. */
if (!Complex) {
*pError = spPANIC;
return NULL;
}
#endif
/* Create Matrix. */
AllocatedSize = MAX( Size, MINIMUM_ALLOCATED_SIZE );
SizePlusOne = (unsigned)(AllocatedSize + 1);
if ((Matrix = SP_MALLOC(struct MatrixFrame, 1)) == NULL) {
*pError = spNO_MEMORY;
return NULL;
}
/* Initialize matrix */
Matrix->ID = SPARSE_ID;
Matrix->Complex = Complex;
Matrix->PreviousMatrixWasComplex = Complex;
Matrix->Factored = NO;
Matrix->Elements = 0;
Matrix->Error = *pError;
Matrix->Originals = 0;
Matrix->Fillins = 0;
Matrix->Reordered = NO;
Matrix->NeedsOrdering = YES;
Matrix->NumberOfInterchangesIsOdd = NO;
Matrix->Partitioned = NO;
Matrix->RowsLinked = NO;
Matrix->InternalVectorsAllocated = NO;
Matrix->SingularCol = 0;
Matrix->SingularRow = 0;
Matrix->Size = Size;
Matrix->AllocatedSize = AllocatedSize;
Matrix->ExtSize = Size;
Matrix->AllocatedExtSize = AllocatedSize;
Matrix->CurrentSize = 0;
Matrix->ExtToIntColMap = NULL;
Matrix->ExtToIntRowMap = NULL;
Matrix->IntToExtColMap = NULL;
Matrix->IntToExtRowMap = NULL;
Matrix->MarkowitzRow = NULL;
Matrix->MarkowitzCol = NULL;
Matrix->MarkowitzProd = NULL;
Matrix->DoCmplxDirect = NULL;
Matrix->DoRealDirect = NULL;
Matrix->Intermediate = NULL;
Matrix->RelThreshold = DEFAULT_THRESHOLD;
Matrix->AbsThreshold = 0.0;
Matrix->TopOfAllocationList = NULL;
Matrix->RecordsRemaining = 0;
Matrix->ElementsRemaining = 0;
Matrix->FillinsRemaining = 0;
RecordAllocation( Matrix, Matrix );
if (Matrix->Error == spNO_MEMORY) goto MemoryError;
/* Take out the trash. */
Matrix->TrashCan.Real = 0.0;
Matrix->TrashCan.Imag = 0.0;
Matrix->TrashCan.Row = 0;
Matrix->TrashCan.Col = 0;
Matrix->TrashCan.NextInRow = NULL;
Matrix->TrashCan.NextInCol = NULL;
#if INITIALIZE
Matrix->TrashCan.pInitInfo = NULL;
#endif
/* Allocate space in memory for Diag pointer vector. */
SP_CALLOC( Matrix->Diag, ElementPtr, SizePlusOne);
if (Matrix->Diag == NULL)
goto MemoryError;
/* Allocate space in memory for FirstInCol pointer vector. */
SP_CALLOC( Matrix->FirstInCol, ElementPtr, SizePlusOne);
if (Matrix->FirstInCol == NULL)
goto MemoryError;
/* Allocate space in memory for FirstInRow pointer vector. */
SP_CALLOC( Matrix->FirstInRow, ElementPtr, SizePlusOne);
if (Matrix->FirstInRow == NULL)
goto MemoryError;
/* Allocate space in memory for IntToExtColMap vector. */
if (( Matrix->IntToExtColMap = SP_MALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Allocate space in memory for IntToExtRowMap vector. */
if (( Matrix->IntToExtRowMap = SP_MALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Initialize MapIntToExt vectors. */
for (I = 1; I <= AllocatedSize; I++)
{
Matrix->IntToExtRowMap[I] = I;
Matrix->IntToExtColMap[I] = I;
}
#if TRANSLATE
/* Allocate space in memory for ExtToIntColMap vector. */
if (( Matrix->ExtToIntColMap = SP_MALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Allocate space in memory for ExtToIntRowMap vector. */
if (( Matrix->ExtToIntRowMap = SP_MALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Initialize MapExtToInt vectors. */
for (I = 1; I <= AllocatedSize; I++) {
Matrix->ExtToIntColMap[I] = -1;
Matrix->ExtToIntRowMap[I] = -1;
}
Matrix->ExtToIntColMap[0] = 0;
Matrix->ExtToIntRowMap[0] = 0;
#endif
/* Allocate space for fill-ins and initial set of elements. */
InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize,
SPACE_FOR_FILL_INS*AllocatedSize );
if (Matrix->Error == spNO_MEMORY)
goto MemoryError;
return Matrix;
MemoryError:
/* Deallocate matrix and return no pointer to matrix if there is not enough
memory. */
*pError = spNO_MEMORY;
spDestroy(Matrix);
return NULL;
}
/*
* ELEMENT ALLOCATION
*
* This routine allocates space for matrix elements. It requests large blocks
* of storage from the system and doles out individual elements as required.
* This technique, as opposed to allocating elements individually, tends to
* speed the allocation process.
*
* >>> Returned:
* A pointer to an element.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to matrix.
*
* >>> Local variables:
* pElement (ElementPtr)
* A pointer to the first element in the group of elements being allocated.
*
* >>> Possible errors:
* spNO_MEMORY
*/
ElementPtr
spcGetElement(MatrixPtr Matrix)
{
ElementPtr pElements;
/* Begin `spcGetElement'. */
#if !COMBINE || STRIP || LINT
/* Allocate block of MatrixElements if necessary. */
if (Matrix->ElementsRemaining == 0) {
pElements = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, pElements );
if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailElement = pElements;
}
#endif
#if COMBINE || STRIP || LINT
if (Matrix->ElementsRemaining == 0)
{
pListNode = Matrix->LastElementListNode;
/* First see if there are any stripped elements left. */
if (pListNode->Next != NULL) {
Matrix->LastElementListNode = pListNode = pListNode->Next;
Matrix->ElementsRemaining = pListNode->NumberOfElementsInList;
Matrix->NextAvailElement = pListNode->pElementList;
} else {
/* Allocate block of elements. */
pElements = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, pElements );
if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailElement = pElements;
/* Allocate an element list structure. */
pListNode->Next = SP_MALLOC(struct ElementListNodeStruct,1);
RecordAllocation( Matrix, pListNode->Next );
if (Matrix->Error == spNO_MEMORY)
return NULL;
Matrix->LastElementListNode = pListNode = pListNode->Next;
pListNode->pElementList = pElements;
pListNode->NumberOfElementsInList = ELEMENTS_PER_ALLOCATION;
pListNode->Next = NULL;
}
}
#endif
/* Update Element counter and return pointer to Element. */
Matrix->ElementsRemaining--;
return Matrix->NextAvailElement++;
}
/*
* ELEMENT ALLOCATION INITIALIZATION
*
* This routine allocates space for matrix fill-ins and an initial
* set of elements. Besides being faster than allocating space for
* elements one at a time, it tends to keep the fill-ins physically
* close to the other matrix elements in the computer memory. This
* keeps virtual memory paging to a minimum.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to the matrix.
* InitialNumberOfElements <input> (int)
* This number is used as the size of the block of memory, in
* MatrixElements, reserved for elements. If more than this number of
* elements are generated, then more space is allocated later.
* NumberOfFillinsExpected <input> (int)
* This number is used as the size of the block of memory, in
* MatrixElements, reserved for fill-ins. If more than this number of
* fill-ins are generated, then more space is allocated, but they may
* not be physically close in computer's memory.
*
* >>> Local variables:
* pElement (ElementPtr)
* A pointer to the first element in the group of elements being allocated.
*
* >>> Possible errors:
* spNO_MEMORY */
static void
InitializeElementBlocks(MatrixPtr Matrix, int InitialNumberOfElements,
int NumberOfFillinsExpected)
{
ElementPtr pElement;
/* Begin `InitializeElementBlocks'. */
/* Allocate block of MatrixElements for elements. */
pElement = SP_MALLOC(struct MatrixElement, InitialNumberOfElements);
RecordAllocation( Matrix, pElement );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->ElementsRemaining = InitialNumberOfElements;
Matrix->NextAvailElement = pElement;
/* Allocate an element list structure. */
Matrix->FirstElementListNode = SP_MALLOC(struct ElementListNodeStruct,1);
RecordAllocation( Matrix, Matrix->FirstElementListNode );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->LastElementListNode = Matrix->FirstElementListNode;
Matrix->FirstElementListNode->pElementList = pElement;
Matrix->FirstElementListNode->NumberOfElementsInList =
InitialNumberOfElements;
Matrix->FirstElementListNode->Next = NULL;
/* Allocate block of MatrixElements for fill-ins. */
pElement = SP_MALLOC(struct MatrixElement, NumberOfFillinsExpected);
RecordAllocation( Matrix, pElement );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->FillinsRemaining = NumberOfFillinsExpected;
Matrix->NextAvailFillin = pElement;
/* Allocate a fill-in list structure. */
Matrix->FirstFillinListNode = SP_MALLOC(struct FillinListNodeStruct,1);
RecordAllocation( Matrix, Matrix->FirstFillinListNode );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
Matrix->FirstFillinListNode->pFillinList = pElement;
Matrix->FirstFillinListNode->NumberOfFillinsInList =NumberOfFillinsExpected;
Matrix->FirstFillinListNode->Next = NULL;
return;
}
/*
* FILL-IN ALLOCATION
*
* This routine allocates space for matrix fill-ins. It requests
* large blocks of storage from the system and doles out individual
* elements as required. This technique, as opposed to allocating
* elements individually, tends to speed the allocation process.
*
* >>> Returned:
* A pointer to the fill-in.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to matrix.
*
* >>> Possible errors:
* spNO_MEMORY */
ElementPtr
spcGetFillin(MatrixPtr Matrix)
{
/* Begin `spcGetFillin'. */
#if !STRIP || LINT
if (Matrix->FillinsRemaining == 0)
return spcGetElement( Matrix );
#endif
#if STRIP || LINT
if (Matrix->FillinsRemaining == 0) {
pListNode = Matrix->LastFillinListNode;
/* First see if there are any stripped fill-ins left. */
if (pListNode->Next != NULL) {
Matrix->LastFillinListNode = pListNode = pListNode->Next;
Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
Matrix->NextAvailFillin = pListNode->pFillinList;
} else {
/* Allocate block of fill-ins. */
pFillins = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, pFillins );
if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailFillin = pFillins;
/* Allocate a fill-in list structure. */
pListNode->Next = SP_MALLOC(struct FillinListNodeStruct,1);
RecordAllocation( Matrix, pListNode->Next );
if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->LastFillinListNode = pListNode = pListNode->Next;
pListNode->pFillinList = pFillins;
pListNode->NumberOfFillinsInList = ELEMENTS_PER_ALLOCATION;
pListNode->Next = NULL;
}
}
#endif
/* Update Fill-in counter and return pointer to Fill-in. */
Matrix->FillinsRemaining--;
return Matrix->NextAvailFillin++;
}
/*
* RECORD A MEMORY ALLOCATION
*
* This routine is used to record all memory allocations so that the
* memory can be freed later.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to the matrix.
* AllocatedPtr <input> (void *)
* The pointer returned by tmalloc or calloc. These pointers are
* saved in a list so that they can be easily freed.
*
* >>> Possible errors:
* spNO_MEMORY */
static void
RecordAllocation(MatrixPtr Matrix, void *AllocatedPtr )
{
/* Begin `RecordAllocation'. */
/* If Allocated pointer is NULL, assume that tmalloc returned a
* NULL pointer, which indicates a spNO_MEMORY error. */
if (AllocatedPtr == NULL) {
Matrix->Error = spNO_MEMORY;
return;
}
/* Allocate block of MatrixElements if necessary. */
if (Matrix->RecordsRemaining == 0) {
AllocateBlockOfAllocationList( Matrix );
if (Matrix->Error == spNO_MEMORY) {
SP_FREE(AllocatedPtr);
return;
}
}
/* Add Allocated pointer to Allocation List. */
(++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr;
Matrix->RecordsRemaining--;
return;
}
/*
* ADD A BLOCK OF SLOTS TO ALLOCATION LIST
*
* This routine increases the size of the allocation list.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to the matrix.
*
* >>> Local variables:
* ListPtr (AllocationListPtr)
* Pointer to the list that contains the pointers to segments of
* memory that were allocated by the operating system for the
* current matrix.
*
* >>> Possible errors:
* spNO_MEMORY */
static void
AllocateBlockOfAllocationList(MatrixPtr Matrix)
{
int I;
AllocationListPtr ListPtr;
/* Begin `AllocateBlockOfAllocationList'. */
/* Allocate block of records for allocation list. */
ListPtr = SP_MALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION+1));
if (ListPtr == NULL) {
Matrix->Error = spNO_MEMORY;
return;
}
/* String entries of allocation list into singly linked list.
List is linked such that any record points to the one before
it. */
ListPtr->NextRecord = Matrix->TopOfAllocationList;
Matrix->TopOfAllocationList = ListPtr;
ListPtr += ELEMENTS_PER_ALLOCATION;
for (I = ELEMENTS_PER_ALLOCATION; I > 0; I--) {
ListPtr->NextRecord = ListPtr - 1;
ListPtr--;
}
/* Record allocation of space for allocation list on allocation list. */
Matrix->TopOfAllocationList->AllocatedPtr = (void *)ListPtr;
Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;
return;
}
/*
* MATRIX DEALLOCATION
*
* Deallocates pointers and elements of Matrix.
*
* >>> Arguments:
* Matrix <input> (void *)
* Pointer to the matrix frame which is to be removed from memory.
*
* >>> Local variables:
* ListPtr (AllocationListPtr)
* Pointer into the linked list of pointers to allocated data structures.
* Points to pointer to structure to be freed.
* NextListPtr (AllocationListPtr)
* Pointer into the linked list of pointers to allocated data structures.
* Points to the next pointer to structure to be freed. This is needed
* because the data structure to be freed could include the current node
* in the allocation list.
*/
void
spDestroy(MatrixPtr Matrix)
{
AllocationListPtr ListPtr, NextListPtr;
/* Begin `spDestroy'. */
assert( IS_SPARSE( Matrix ) );
/* Deallocate the vectors that are located in the matrix frame. */
SP_FREE( Matrix->IntToExtColMap );
SP_FREE( Matrix->IntToExtRowMap );
SP_FREE( Matrix->ExtToIntColMap );
SP_FREE( Matrix->ExtToIntRowMap );
SP_FREE( Matrix->Diag );
SP_FREE( Matrix->FirstInRow );
SP_FREE( Matrix->FirstInCol );
SP_FREE( Matrix->MarkowitzRow );
SP_FREE( Matrix->MarkowitzCol );
SP_FREE( Matrix->MarkowitzProd );
SP_FREE( Matrix->DoCmplxDirect );
SP_FREE( Matrix->DoRealDirect );
SP_FREE( Matrix->Intermediate );
/* Sequentially step through the list of allocated pointers
* freeing pointers along the way. */
ListPtr = Matrix->TopOfAllocationList;
while (ListPtr != NULL) {
NextListPtr = ListPtr->NextRecord;
if ((void *) ListPtr == ListPtr->AllocatedPtr) {
SP_FREE( ListPtr );
} else {
SP_FREE( ListPtr->AllocatedPtr );
}
ListPtr = NextListPtr;
}
return;
}
/*
* RETURN MATRIX ERROR STATUS
*
* This function is used to determine the error status of the given
* matrix.
*
* >>> Returned:
* The error status of the given matrix.
*
* >>> Arguments:
* Matrix <input> (void *)
* The matrix for which the error status is desired. */
int
spError(MatrixPtr Matrix )
{
/* Begin `spError'. */
if (Matrix != NULL) {
assert(Matrix->ID == SPARSE_ID);
return Matrix->Error;
} else {
/* This error may actually be spPANIC, no way to tell. */
return spNO_MEMORY;
}
}
/*
* 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 <input> (void *)
* The matrix for which the error status is desired.
* pRow <output> (int *)
* The row number.
* pCol <output> (int *)
* The column number.
*/
void
spWhereSingular(MatrixPtr Matrix, int *pRow, int *pCol)
{
/* Begin `spWhereSingular'. */
assert( IS_SPARSE( Matrix ) );
if (Matrix->Error == spSINGULAR || Matrix->Error == spZERO_DIAG)
{
*pRow = Matrix->SingularRow;
*pCol = Matrix->SingularCol;
}
else *pRow = *pCol = 0;
return;
}
/*
* MATRIX SIZE
*
* Returns the size of the matrix. Either the internal or external size of
* the matrix is returned.
*
* >>> Arguments:
* Matrix <input> (void *)
* Pointer to matrix.
* External <input> (int)
* 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
spGetSize(MatrixPtr Matrix, int External)
{
/* Begin `spGetSize'. */
assert( IS_SPARSE( Matrix ) );
#if TRANSLATE
if (External)
return Matrix->ExtSize;
else
return Matrix->Size;
#else
return Matrix->Size;
#endif
}
/*
* SET MATRIX COMPLEX OR REAL
*
* Forces matrix to be either real or complex.
*
* >>> Arguments:
* Matrix <input> (void *)
* Pointer to matrix.
*/
void
spSetReal(MatrixPtr Matrix)
{
/* Begin `spSetReal'. */
assert( IS_SPARSE( Matrix ));
Matrix->Complex = NO;
return;
}
void
spSetComplex(MatrixPtr Matrix)
{
/* Begin `spSetComplex'. */
assert( IS_SPARSE( Matrix ));
Matrix->Complex = YES;
return;
}
/*
* ELEMENT, FILL-IN OR ORIGINAL COUNT
*
* Two functions used to return simple statistics. Either the number
* of total elements, or the number of fill-ins, or the number
* of original elements can be returned.
*
* >>> Arguments:
* Matrix <input> (void *)
* Pointer to matrix.
*/
int
spFillinCount(MatrixPtr Matrix)
{
/* Begin `spFillinCount'. */
assert( IS_SPARSE( Matrix ) );
return Matrix->Fillins;
}
int
spElementCount(MatrixPtr Matrix)
{
/* Begin `spElementCount'. */
assert( IS_SPARSE( Matrix ) );
return Matrix->Elements;
}
int
spOriginalCount(MatrixPtr Matrix)
{
/* Begin `spOriginalCount'. */
assert( IS_SPARSE( Matrix ) );
return Matrix->Originals;
}

1206
src/maths/sparse/spbuild.c
View File

File diff suppressed because it is too large Load Diff

402
src/maths/sparse/spconfig.h
View File

@ -1,402 +0,0 @@
/*
* CONFIGURATION MACRO DEFINITIONS for sparse matrix routines
*
* Author: Advising professor:
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
* U.C. Berkeley
*
* This file contains macros for the sparse matrix routines that are used
* to define the personality of the routines. The user is expected to
* modify this file to maximize the performance of the routines with
* his/her matrices.
*
* Macros are distinguished by using solely capital letters in their
* identifiers. This contrasts with C defined identifiers which are
* strictly lower case, and program variable and procedure names which use
* both upper and lower case.
*/
/*
* Revision and copyright information.
*
* Copyright (c) 1985,86,87,88,89,90
* 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 ngspice_SPCONFIG_H
#define ngspice_SPCONFIG_H
#ifdef spINSIDE_SPARSE
/*
* OPTIONS
*
* These are compiler options. Set each option to one to compile that
* section of the code. If a feature is not desired, set the macro
* to NO. Recommendations are given in brackets, [ignore them].
*
* >>> Option descriptions:
* Arithmetic Precision
* The precision of the arithmetic used by Sparse can be set by
* changing changing the spREAL macro. This macro is
* contained in the file spMatrix.h. It is strongly suggested to
* used double precision with circuit simulators. Note that
* because C always performs arithmetic operations in double
* precision, the only benefit to using single precision is that
* less storage is required. There is often a noticeable speed
* penalty when using single precision. Sparse internally refers
* to a spREAL as a RealNumber.
* EXPANDABLE
* Setting this compiler flag true (1) makes the matrix
* expandable before it has been factored. If the matrix is
* expandable, then if an element is added that would be
* considered out of bounds in the current matrix, the size of
* the matrix is increased to hold that element. As a result,
* the size of the matrix need not be known before the matrix is
* built. The matrix can be allocated with size zero and
* expanded.
* TRANSLATE
* This option allows the set of external row and column numbers
* to be non-packed. In other words, the row and column numbers
* do not have to be contiguous. The priced paid for this
* flexibility is that when TRANSLATE is set true, the time
* required to initially build the matrix will be greater because
* the external row and column number must be translated into
* internal equivalents. This translation brings about other
* benefits though. First, the spGetElement() and
* spGetAdmittance() routines may be used after the matrix has
* been factored. Further, elements, and even rows and columns,
* may be added to the matrix, and row and columns may be deleted
* from the matrix, after it has been factored. Note that when
* the set of row and column number is not a packed set, neither
* are the RHS and Solution vectors. Thus the size of these
* vectors must be at least as large as the external size, which
* is the value of the largest given row or column numbers.
* INITIALIZE
* Causes the spInitialize(), spGetInitInfo(), and
* spInstallInitInfo() routines to be compiled. These routines
* allow the user to store and read one pointer in each nonzero
* element in the matrix. spInitialize() then calls a user
* specified function for each structural nonzero in the matrix,
* and includes this pointer as well as the external row and
* column numbers as arguments. This allows the user to write
* custom matrix initialization routines.
* DIAGONAL_PIVOTING
* Many matrices, and in particular node- and modified-node
* admittance matrices, tend to be nearly symmetric and nearly
* diagonally dominant. For these matrices, it is a good idea to
* select pivots from the diagonal. With this option enabled,
* this is exactly what happens, though if no satisfactory pivot
* can be found on the diagonal, an off-diagonal pivot will be
* used. If this option is disabled, Sparse does not
* preferentially search the diagonal. Because of this, Sparse
* has a wider variety of pivot candidates available, and so
* presumably fewer fill-ins will be created. However, the
* initial pivot selection process will take considerably longer.
* If working with node admittance matrices, or other matrices
* with a strong diagonal, it is probably best to use
* DIAGONAL_PIVOTING for two reasons. First, accuracy will be
* better because pivots will be chosen from the large diagonal
* elements, thus reducing the chance of growth. Second, a near
* optimal ordering will be chosen quickly. If the class of
* matrices you are working with does not have a strong diagonal,
* do not use DIAGONAL_PIVOTING, but consider using a larger
* threshold. When DIAGONAL_PIVOTING is turned off, the following
* options and constants are not used: MODIFIED_MARKOWITZ,
* MAX_MARKOWITZ_TIES, and TIES_MULTIPLIER.
* MODIFIED_MARKOWITZ
* This specifies that the modified Markowitz method of pivot
* selection is to be used. The modified Markowitz method differs
* from standard Markowitz in two ways. First, under modified
* Markowitz, the search for a pivot can be terminated early if a
* adequate (in terms of sparsity) pivot candidate is found.
* Thus, when using modified Markowitz, the initial factorization
* can be faster, but at the expense of a suboptimal pivoting
* order that may slow subsequent factorizations. The second
* difference is in the way modified Markowitz breaks Markowitz
* ties. When two or more elements are pivot candidates and they
* all have the same Markowitz product, then the tie is broken by
* choosing the element that is best numerically. The numerically
* best element is the one with the largest ratio of its magnitude
* to the magnitude of the largest element in the same column,
* excluding itself. The modified Markowitz method results in
* marginally better accuracy. This option is most appropriate
* for use when working with very large matrices where the initial
* factor time represents an unacceptable burden. [NO]
* DELETE
* This specifies that the spDeleteRowAndCol() routine
* should be compiled. Note that for this routine to be
* compiled, both DELETE and TRANSLATE should be set true.
* STRIP
* This specifies that the spStripFills() routine should be compiled.
* MODIFIED_NODAL
* This specifies that the routine that preorders modified node
* admittance matrices should be compiled. This routine results
* in greater speed and accuracy if used with this type of
* matrix.
* QUAD_ELEMENT
* This specifies that the routines that allow four related
* elements to be entered into the matrix at once should be
* compiled. These elements are usually related to an
* admittance. The routines affected by QUAD_ELEMENT are the
* spGetAdmittance, spGetQuad and spGetOnes routines.
* TRANSPOSE
* This specifies that the routines that solve the matrix as if
* it was transposed should be compiled. These routines are
* useful when performing sensitivity analysis using the adjoint
* method.
* SCALING
* This specifies that the routine that performs scaling on the
* matrix should be complied. Scaling is not strongly
* supported. The routine to scale the matrix is provided, but
* no routines are provided to scale and descale the RHS and
* Solution vectors. It is suggested that if scaling is desired,
* it only be preformed when the pivot order is being chosen [in
* spOrderAndFactor()]. This is the only time scaling has
* an effect. The scaling may then either be removed from the
* solution by the user or the scaled factors may simply be
* thrown away. [NO]
* DOCUMENTATION
* This specifies that routines that are used to document the
* matrix, such as spPrint() and spFileMatrix(), should be
* compiled.
* DETERMINANT
* This specifies that the routine spDeterminant() should be complied.
* STABILITY
* This specifies that spLargestElement() and spRoundoff() should
* be compiled. These routines are used to check the stability (and
* hence the quality of the pivoting) of the factorization by
* computing a bound on the size of the element is the matrix E =
* A - LU. If this bound is very high after applying
* spOrderAndFactor(), then the pivot threshold should be raised.
* If the bound increases greatly after using spFactor(), then the
* matrix should probably be reordered.
* CONDITION
* This specifies that spCondition() and spNorm(), the code that
* computes a good estimate of the condition number of the matrix,
* should be compiled.
* PSEUDOCONDITION
* This specifies that spPseudoCondition(), the code that computes
* a crude and easily fooled indicator of ill-conditioning in the
* matrix, should be compiled.
* MULTIPLICATION
* This specifies that the routines to multiply the unfactored
* matrix by a vector should be compiled.
* DEBUG
* This specifies that additional error checking will be compiled.
* The type of error checked are those that are common when the
* matrix routines are first integrated into a user's program. Once
* the routines have been integrated in and are running smoothly, this
* option should be turned off.
*/
/* Begin options. */
#define EXPANDABLE YES
#define TRANSLATE YES
#define INITIALIZE NO
#define DIAGONAL_PIVOTING YES
#define MODIFIED_MARKOWITZ NO
#define DELETE NO
#define STRIP NO
#define MODIFIED_NODAL YES
#define QUAD_ELEMENT NO
#define TRANSPOSE YES
#define SCALING NO
#define DOCUMENTATION YES
#define MULTIPLICATION YES
#define DETERMINANT YES
#define DETERMINANT2 YES
#define STABILITY NO
#define CONDITION NO
#define PSEUDOCONDITION NO
#ifdef HAS_MINDATA
# define DEBUG NO
#else
# define DEBUG YES
#endif
/*
* The following options affect Sparse exports and so are exported as a
* side effect. For this reason they use the `sp' prefix. The boolean
* constants YES an NO are not defined in spMatrix.h to avoid conflicts
* with user code, so use 0 for NO and 1 for YES.
*/
/*
* MATRIX CONSTANTS
*
* These constants are used throughout the sparse matrix routines. They
* should be set to suit the type of matrix being solved. Recommendations
* are given in brackets.
*
* Some terminology should be defined. The Markowitz row count is the number
* of non-zero elements in a row excluding the one being considered as pivot.
* There is one Markowitz row count for every row. The Markowitz column
* is defined similarly for columns. The Markowitz product for an element
* is the product of its row and column counts. It is a measure of how much
* work would be required on the next step of the factorization if that
* element were chosen to be pivot. A small Markowitz product is desirable.
*
* >>> Constants descriptions:
* DEFAULT_THRESHOLD
* The relative threshold used if the user enters an invalid
* threshold. Also the threshold used by spFactor() when
* calling spOrderAndFactor(). The default threshold should
* not be less than or equal to zero nor larger than one. [0.001]
* DIAG_PIVOTING_AS_DEFAULT
* This indicates whether spOrderAndFactor() should use diagonal
* pivoting as default. This issue only arises when
* spOrderAndFactor() is called from spFactor().
* SPACE_FOR_ELEMENTS
* This number multiplied by the size of the matrix equals the number
* of elements for which memory is initially allocated in
* spCreate(). [6]
* SPACE_FOR_FILL_INS
* This number multiplied by the size of the matrix equals the number
* of elements for which memory is initially allocated and specifically
* reserved for fill-ins in spCreate(). [4]
* ELEMENTS_PER_ALLOCATION
* The number of matrix elements requested from the tmalloc utility on
* each call to it. Setting this value greater than 1 reduces the
* amount of overhead spent in this system call. On a virtual memory
* machine, its good to allocate slightly less than a page worth of
* elements at a time (or some multiple thereof).
* [For the VAX, for real only use 41, otherwise use 31]
* MINIMUM_ALLOCATED_SIZE
* The minimum allocated size of a matrix. Note that this does not
* limit the minimum size of a matrix. This just prevents having to
* resize a matrix many times if the matrix is expandable, large and
* allocated with an estimated size of zero. This number should not
* be less than one.
* EXPANSION_FACTOR
* The amount the allocated size of the matrix is increased when it
* is expanded.
* MAX_MARKOWITZ_TIES
* This number is used for two slightly different things, both of which
* relate to the search for the best pivot. First, it is the maximum
* number of elements that are Markowitz tied that will be sifted
* through when trying to find the one that is numerically the best.
* Second, it creates an upper bound on how large a Markowitz product
* can be before it eliminates the possibility of early termination
* of the pivot search. In other words, if the product of the smallest
* Markowitz product yet found and TIES_MULTIPLIER is greater than
* MAX_MARKOWITZ_TIES, then no early termination takes place.
* Set MAX_MARKOWITZ_TIES to some small value if no early termination of
* the pivot search is desired. An array of RealNumbers is allocated
* of size MAX_MARKOWITZ_TIES so it must be positive and shouldn't
* be too large. Active when MODIFIED_MARKOWITZ is 1 (true). [100]
* TIES_MULTIPLIER
* Specifies the number of Markowitz ties that are allowed to occur
* before the search for the pivot is terminated early. Set to some
* large value if no early termination of the pivot search is desired.
* This number is multiplied times the Markowitz product to determine
* how many ties are required for early termination. This means that
* more elements will be searched before early termination if a large
* number of fill-ins could be created by accepting what is currently
* considered the best choice for the pivot. Active when
* MODIFIED_MARKOWITZ is 1 (true). Setting this number to zero
* effectively eliminates all pivoting, which should be avoided.
* This number must be positive. TIES_MULTIPLIER is also used when
* diagonal pivoting breaks down. [5]
* DEFAULT_PARTITION
* Which partition mode is used by spPartition() as default.
* Possibilities include
* spDIRECT_PARTITION -- each row used direct addressing, best for
* a few relatively dense matrices.
* spINDIRECT_PARTITION -- each row used indirect addressing, best
* for a few very sparse matrices.
* spAUTO_PARTITION -- direct or indirect addressing is chosen on
* a row-by-row basis, carries a large overhead, but speeds up
* both dense and sparse matrices, best if there is a large
* number of matrices that can use the same ordering.
*/
/* Begin constants. */
#define DEFAULT_THRESHOLD 1.0e-3
#define DIAG_PIVOTING_AS_DEFAULT YES
#define SPACE_FOR_ELEMENTS 6
#define SPACE_FOR_FILL_INS 4
#define ELEMENTS_PER_ALLOCATION 31
#define MINIMUM_ALLOCATED_SIZE 6
#define EXPANSION_FACTOR 1.5
#define MAX_MARKOWITZ_TIES 100
#define TIES_MULTIPLIER 5
#define DEFAULT_PARTITION spAUTO_PARTITION
/*
* PRINTER WIDTH
*
* This macro characterize the printer for the spPrint() routine.
*
* >>> Macros:
* PRINTER_WIDTH
* The number of characters per page width. Set to 80 for terminal,
* 132 for line printer.
*/
/* Begin printer constants. */
#define PRINTER_WIDTH 80
/*
* MACHINE CONSTANTS
*
* These numbers must be updated when the program is ported to a new machine.
*/
/* Begin machine constants. */
/*
* Grab from Spice include files
*/
#define MACHINE_RESOLUTION DBL_EPSILON
#define LARGEST_REAL DBL_MAX
#define SMALLEST_REAL DBL_MIN
#define LARGEST_SHORT_INTEGER SHRT_MAX
#define LARGEST_LONG_INTEGER LONG_MAX
/*
* ANNOTATION
*
* This macro changes the amount of annotation produced by the matrix
* routines. The annotation is used as a debugging aid. Change the number
* associated with ANNOTATE to change the amount of annotation produced by
* the program.
*/
/* Begin annotation definitions. */
#define ANNOTATE NONE
#define NONE 0
#define ON_STRANGE_BEHAVIOR 1
#define FULL 2
#endif /* spINSIDE_SPARSE */
#endif

5
src/maths/sparse/spextra.c
View File

@ -23,10 +23,9 @@
*/
#define spINSIDE_SPARSE
#include "spconfig.h"
#include "spConfig.h"
#include "ngspice/spmatrix.h"
#include "spdefs.h"
#include "spDefs.h"
void
spConstMult(MatrixPtr matrix, double constant)