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

* src/include/spconfig.h: Removed spCOMPLEX, 

spSEPARATED_COMPLEX_VECTORS and spCOMPATIBILITY defines.  This
	made including this file from src/include/spmatrix.h unnecessary.
	Moved this file to src/maths/sparse/spconfig.h.

	* src/include/spmatrix.h: Removed include of
	src/include/spconfig.h.

	* src/maths/sparse/spalloc.c, src/maths/sparse/spbuild.c,
	src/maths/sparse/spcombin.c, src/maths/sparse/spdefs.h,
	src/maths/sparse/spfactor.c, src/maths/sparse/spoutput.c,
	src/maths/sparse/spsmp.c, src/maths/sparse/spsolve.c,
	src/maths/sparse/sputils.c: The other files affected by the
	removal of spCOMPLEX, spSEPARATED_COMPLEX_VECTORS and
	spCOMPATIBILITY defines.  Also: assertions are enabled by
	default.

	* src/include/smpdefs.h, src/maths/sparse/spsmp.c: SMPmatrix is
	now a typedef for void, instead of char.  Updated all function
	declarations to match this.  Also added function prototypes not
	previously mentioned in src/include/smpdefs.h.

	* src/include/complex.h: Updates of cast from char * to void *
This commit is contained in:
arno 2000-07-03 15:28:50 +00:00
parent 8aa1ceead5
commit 7ac9278389
15 changed files with 2199 additions and 2814 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/include/Makefile.am
View File

@ -47,7 +47,6 @@ noinst_HEADERS = \
sensgen.h \
sim.h \
smpdefs.h \
spconfig.h \
sperror.h \
spmatrix.h \
suffix.h \

11
src/include/complex.h
View File

@ -269,12 +269,9 @@ typedef struct
(A).imag = (B.imag) - (C.imag); \
}
/* Macro function that returns the approx absolute value of a complex number. */
#if spCOMPLEX
/* Macro function that returns the approx absolute value of a complex
number. */
#define ELEMENT_MAG(ptr) (ABS((ptr)->Real) + ABS((ptr)->Imag))
#else
#define ELEMENT_MAG(ptr) ((ptr)->Real < 0.0 ? -(ptr)->Real : (ptr)->Real)
#endif
/* Complex assignment statements. */
#define CMPLX_ASSIGN(to,from) \
@ -481,8 +478,8 @@ typedef struct
/* Complex reciprocation: to = 1.0 / den */
#define CMPLX_RECIPROCAL(to,den) \
{ RealNumber r_; \
if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \
((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \
if (((den).Real >= (den).Imag && (den).Real > -(den).Imag) || \
((den).Real < (den).Imag && (den).Real <= -(den).Imag)) \
{ r_ = (den).Imag / (den).Real; \
(to).Imag = -r_*((to).Real = 1.0/((den).Real + r_*(den).Imag)); \
} \

66
src/include/smpdefs.h
View File

@ -1,7 +1,7 @@
#ifndef SMP
#define SMP
typedef char SMPmatrix;
typedef void SMPmatrix;
typedef struct MatrixElement *SMPelement;
/**********
@ -12,52 +12,34 @@ Author: 1985 Thomas L. Quarles
#include "complex.h"
#include <stdio.h>
#ifdef __STDC__
int SMPaddElt( SMPmatrix *, int , int , double );
void SMPcClear( SMPmatrix *);
int SMPcLUfac( SMPmatrix *, double );
int SMPcProdDiag( SMPmatrix *, SPcomplex *, int *);
int SMPcReorder( SMPmatrix * , double , double , int *);
void SMPcSolve( SMPmatrix *, double [], double [], double [], double []);
void SMPclear( SMPmatrix *);
void SMPcolSwap( SMPmatrix * , int , int );
void SMPdestroy( SMPmatrix *);
int SMPfillup( SMPmatrix * );
SMPelement * SMPfindElt( SMPmatrix *, int , int , int );
void SMPgetError( SMPmatrix *, int *, int *);
int SMPluFac( SMPmatrix *, double , double );
double * SMPmakeElt( SMPmatrix * , int , int );
void SMPcClear( SMPmatrix *);
void SMPclear( SMPmatrix *);
int SMPcLUfac( SMPmatrix *, double );
int SMPluFac( SMPmatrix *, double , double );
int SMPcReorder( SMPmatrix * , double , double , int *);
int SMPreorder( SMPmatrix * , double , double , double );
void SMPcaSolve(SMPmatrix *Matrix, double RHS[], double iRHS[],
double Spare[], double iSpare[]);
void SMPcSolve( SMPmatrix *, double [], double [], double [], double []);
void SMPsolve( SMPmatrix *, double [], double []);
int SMPmatSize( SMPmatrix *);
int SMPnewMatrix( SMPmatrix ** );
int SMPnewNode( int , SMPmatrix *);
void SMPdestroy( SMPmatrix *);
int SMPpreOrder( SMPmatrix *);
void SMPprint( SMPmatrix * , FILE *);
int SMPreorder( SMPmatrix * , double , double , double );
void SMProwSwap( SMPmatrix * , int , int );
void SMPsolve( SMPmatrix *, double [], double []);
#else /* stdc */
int SMPaddElt();
void SMPcClear();
int SMPcLUfac();
int SMPcProdDiag();
int SMPcReorder();
void SMPcSolve();
void SMPclear();
void SMPcolSwap();
void SMPdestroy();
int SMPfillup();
SMPelement * SMPfindElt();
void SMPgetError();
int SMPluFac();
double * SMPmakeElt();
int SMPmatSize();
int SMPnewMatrix();
int SMPnewNode();
int SMPpreOrder();
void SMPprint();
int SMPreorder();
void SMProwSwap();
void SMPsolve();
#endif /* stdc */
void SMPgetError( SMPmatrix *, int *, int *);
int SMPcProdDiag( SMPmatrix *, SPcomplex *, int *);
int SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent);
SMPelement * SMPfindElt( SMPmatrix *, int , int , int );
int SMPcZeroCol(SMPmatrix *eMatrix, int Col);
int SMPcAddCol(SMPmatrix *eMatrix, int Accum_Col, int Addend_Col);
int SMPzeroRow(SMPmatrix *eMatrix, int Row);
#ifdef PARALLEL_ARCH
void SMPcombine(SMPmatrix *Matrix, double RHS[], double Spare[]);
void SMPcCombine(SMPmatrix *Matrix, double RHS[], double Spare[],
double iRHS[], double iSpare[]);
#endif
#endif /*SMP*/

172
src/include/spmatrix.h
View File

@ -42,17 +42,6 @@
#ifndef spOKAY
/*
* IMPORTS
*
* >>> Import descriptions:
* spConfig.h
* Macros that customize the sparse matrix routines.
*/
#include "spconfig.h"
@ -76,8 +65,7 @@
* 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(). In spCOMPATIBILITY mode, spZERO_DIAG is
* indistinguishable from spSINGULAR.
* spOrderAndFactor().
* spSINGULAR
* Fatal error. Matrix is singular, so no unique solution exists.
* spNO_MEMORY
@ -108,21 +96,6 @@
#define spFATAL E_BADMATRIX
#if spCOMPATIBILITY
#define NO_ERROR spOKAY
#define UNDER_FLOW spOKAY
#define OVER_FLOW spOKAY
#define ILL_CONDITIONED spSMALL_PIVOT
#define SINGULAR spSINGULAR
#define NO_MEMORY spNO_MEMORY
#define RANGE spPANIC
#define FATAL spFATAL
#undef spZERO_DIAG
#define spZERO_DIAG spSINGULAR
#endif /* spCOMPATIBILITY */
@ -145,10 +118,6 @@
#define spREAL double
/* #define void int */
#if spCOMPATIBILITY
#define SPARSE_REAL spREAL
#endif
/*
@ -242,15 +211,6 @@
*((template).Element4Negated+1) -= imag; \
}
#if spCOMPATIBILITY
#define ADD_REAL_ELEMENT_TO_MATRIX spADD_REAL_ELEMENT
#define ADD_IMAG_ELEMENT_TO_MATRIX spADD_IMAG_ELEMENT
#define ADD_COMPLEX_ELEMENT_TO_MATRIX spADD_COMPLEX_ELEMENT
#define ADD_REAL_QUAD_ELEMENT_TO_MATRIX spADD_REAL_ELEMENT
#define ADD_IMAG_QUAD_ELEMENT_TO_MATRIX spADD_IMAG_ELEMENT
#define ADD_COMPLEX_QUAD_ELEMENT_TO_MATRIX spADD_COMPLEX_ELEMENT
#endif
@ -272,10 +232,6 @@
* before being placed in Element3 and Element4.
*/
#if spCOMPATIBILITY
#define spTemplate TemplateStruct
#endif
/* Begin `spTemplate'. */
struct spTemplate
{ spREAL *Element1 ;
@ -296,97 +252,47 @@ struct spTemplate
/* Begin function declarations. */
extern void spClear( char* );
extern spREAL spCondition( char*, spREAL, int* );
extern char *spCreate( int, int, int* );
extern void spDeleteRowAndCol( char*, int, int );
extern void spDestroy( char* );
extern int spElementCount( char* );
extern int spError( char* );
extern int spFactor( char* );
extern int spFileMatrix( char*, char*, char*, int, int, int );
extern int spFileStats( char*, char*, char* );
extern int spFillinCount( char* );
extern int spGetAdmittance( char*, int, int, struct spTemplate* );
extern spREAL *spGetElement( char*, int, int );
extern void spClear( void * );
extern spREAL spCondition( void *, spREAL, int* );
extern void *spCreate( int, int, int* );
extern void spDeleteRowAndCol( void *, int, int );
extern void spDestroy(void *);
extern int spElementCount( void * );
extern int spError( void * );
extern int spFactor( void * );
extern int spFileMatrix( void *, char *, char *, int, int, int );
extern int spFileStats( void *, char *, char * );
extern int spFillinCount( void * );
extern int spGetAdmittance( void *, int, int, struct spTemplate* );
extern spREAL *spGetElement( void *, int, int );
extern char *spGetInitInfo( spREAL* );
extern int spGetOnes( char*, int, int, int, struct spTemplate* );
extern int spGetQuad( char*, int, int, int, int, struct spTemplate* );
extern int spGetSize( char*, int );
extern int spInitialize( char*, int (*)() );
extern void spInstallInitInfo( spREAL*, char* );
extern spREAL spLargestElement( char* );
extern void spMNA_Preorder( char* );
extern spREAL spNorm( char* );
extern int spOrderAndFactor( char*, spREAL*, spREAL, spREAL, int );
extern int spOriginalCount( char *);
extern void spPartition( char*, int );
extern void spPrint( char*, int, int, int );
extern spREAL spPseudoCondition( char* );
extern spREAL spRoundoff( char*, spREAL );
extern void spScale( char*, spREAL*, spREAL* );
extern void spSetComplex( char* );
extern void spSetReal( char* );
extern void spStripFills( char* );
extern void spWhereSingular( char*, int*, int* );
extern int spGetOnes( void *, int, int, int, struct spTemplate* );
extern int spGetQuad( void *, int, int, int, int, struct spTemplate* );
extern int spGetSize( void *, int );
extern int spInitialize( void *, int (*)() );
extern void spInstallInitInfo( spREAL*, void * );
extern spREAL spLargestElement( void * );
extern void spMNA_Preorder( void * );
extern spREAL spNorm( void * );
extern int spOrderAndFactor( void *, spREAL*, spREAL, spREAL, int );
extern int spOriginalCount( void *);
extern void spPartition( void *, int );
extern void spPrint( void *, int, int, int );
extern spREAL spPseudoCondition( void * );
extern spREAL spRoundoff( void *, spREAL );
extern void spScale( void *, spREAL*, spREAL* );
extern void spSetComplex( void * );
extern void spSetReal( void * );
extern void spStripFills( void * );
extern void spWhereSingular( void *, int*, int* );
/* Functions with argument lists that are dependent on options. */
#if spCOMPLEX
extern void spDeterminant ( char*, int*, spREAL*, spREAL* );
#else /* NOT spCOMPLEX */
extern void spDeterminant ( char*, int*, spREAL* );
#endif /* NOT spCOMPLEX */
#if spCOMPLEX && spSEPARATED_COMPLEX_VECTORS
extern int spFileVector( char*, char* , spREAL*, spREAL*);
extern void spMultiply( char*, spREAL*, spREAL*, spREAL*, spREAL* );
extern void spMultTransposed(char*,spREAL*,spREAL*,spREAL*,spREAL*);
extern void spSolve( char*, spREAL*, spREAL*, spREAL*, spREAL* );
extern void spSolveTransposed(char*,spREAL*,spREAL*,spREAL*,spREAL*);
#else /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */
extern int spFileVector( char*, char* , spREAL* );
extern void spMultiply( char*, spREAL*, spREAL* );
extern void spMultTransposed( char*, spREAL*, spREAL* );
extern void spSolve( char*, spREAL*, spREAL* );
extern void spSolveTransposed( char*, spREAL*, spREAL* );
#endif /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */
#if spCOMPATIBILITY
extern char *AllocateMatrix();
extern spREAL *AddElementToMatrix();
extern void AddRealElementToMatrix();
extern void AddImagElementToMatrix();
extern void AddComplexElementToMatrix();
extern void AddAdmittanceToMatrix();
extern void AddOnesToMatrix();
extern void AddQuadToMatrix();
extern void AddRealQuadElementToMatrix();
extern void AddImagQuadElementToMatrix();
extern void AddComplexQuadElementToMatrix();
extern void CleanMatrix();
extern void ClearMatrix();
extern int ClearMatrixError();
extern void DeallocateMatrix();
extern void DeleteRowAndColFromMatrix();
extern void Determinant();
extern int DecomposeMatrix();
extern int GetMatrixSize();
extern int MatrixElementCount();
extern int MatrixFillinCount();
extern void MatrixMultiply();
extern spREAL MatrixRoundoffError();
extern int MatrixError();
extern int OrderAndDecomposeMatrix();
extern void OutputMatrixToFile();
extern void OutputStatisticsToFile();
extern void OutputVectorToFile();
extern void PreorderForModifiedNodal();
extern void PrintMatrix();
extern void SetMatrixComplex();
extern void SetMatrixReal();
extern void SolveMatrix();
extern void SolveTransposedMatrix();
extern void ScaleMatrix();
#endif /* spCOMPATIBILITY */
extern void spDeterminant ( void *, int*, spREAL*, spREAL* );
extern int spFileVector( void *, char * , spREAL*, spREAL*);
extern void spMultiply( void *, spREAL*, spREAL*, spREAL*, spREAL* );
extern void spMultTransposed(void *,spREAL*,spREAL*,spREAL*,spREAL*);
extern void spSolve( void *, spREAL*, spREAL*, spREAL*, spREAL* );
extern void spSolveTransposed(void *,spREAL*,spREAL*,spREAL*,spREAL*);
#endif /* spOKAY */

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

@ -3,16 +3,17 @@
noinst_LIBRARIES = libsparse.a
libsparse_a_SOURCES = \
spalloc.c \
spbuild.c \
spcombin.c \
spdefs.h \
spextra.c \
spfactor.c \
spoutput.c \
spsmp.c \
spsolve.c \
sputils.c
spalloc.c \
spbuild.c \
spcombin.c \
spconfig.h \
spdefs.h \
spextra.c \
spfactor.c \
spoutput.c \
spsmp.c \
spsolve.c \
sputils.c

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

@ -58,6 +58,7 @@
* spDefs.h
* Matrix type and macro definitions for the sparse matrix routines.
*/
#include <assert.h>
#define spINSIDE_SPARSE
#include "spconfig.h"
@ -75,15 +76,9 @@
* Function declarations
*/
#ifdef __STDC__
static void InitializeElementBlocks( MatrixPtr, int, int );
static void RecordAllocation( MatrixPtr, char* );
static void RecordAllocation( MatrixPtr, void *);
static void AllocateBlockOfAllocationList( MatrixPtr );
#else /* __STDC__ */
static void InitializeElementBlocks();
static void RecordAllocation();
static void AllocateBlockOfAllocationList();
#endif /* __STDC__ */
@ -124,51 +119,42 @@ static void AllocateBlockOfAllocationList();
* Error is cleared in this routine.
*/
char *
spCreate( Size, Complex, pError )
int Size, *pError;
BOOLEAN Complex;
void *
spCreate(int Size, int Complex, int *pError)
{
register unsigned SizePlusOne;
register MatrixPtr Matrix;
register int I;
int AllocatedSize;
unsigned SizePlusOne;
MatrixPtr Matrix;
int I;
int AllocatedSize;
/* Begin `spCreate'. */
/* Clear error flag. */
/* Begin `spCreate'. */
/* Clear error flag. */
*pError = spOKAY;
/* Test for valid size. */
if ((Size < 0) OR (Size == 0 AND NOT EXPANDABLE))
{ *pError = spPANIC;
/* Test for valid size. */
if ((Size < 0) || (Size == 0 && !EXPANDABLE)) {
*pError = spPANIC;
return NULL;
}
/* Test for valid type. */
#if NOT spCOMPLEX
if (Complex)
{ *pError = spPANIC;
return NULL;
}
#endif
#if NOT REAL
if (NOT Complex)
{ *pError = spPANIC;
return NULL;
/* Test for valid type. */
#if !REAL
if (!Complex) {
*pError = spPANIC;
return NULL;
}
#endif
/* Create Matrix. */
/* Create Matrix. */
AllocatedSize = MAX( Size, MINIMUM_ALLOCATED_SIZE );
SizePlusOne = (unsigned)(AllocatedSize + 1);
if ((Matrix = ALLOC(struct MatrixFrame, 1)) == NULL)
{ *pError = spNO_MEMORY;
return NULL;
if ((Matrix = ALLOC(struct MatrixFrame, 1)) == NULL) {
*pError = spNO_MEMORY;
return NULL;
}
/* Initialize matrix */
/* Initialize matrix */
Matrix->ID = SPARSE_ID;
Matrix->Complex = Complex;
Matrix->PreviousMatrixWasComplex = Complex;
@ -208,14 +194,12 @@ int AllocatedSize;
Matrix->ElementsRemaining = 0;
Matrix->FillinsRemaining = 0;
RecordAllocation( Matrix, (char *)Matrix );
RecordAllocation( Matrix, (void *)Matrix );
if (Matrix->Error == spNO_MEMORY) goto MemoryError;
/* Take out the trash. */
/* 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;
@ -224,67 +208,68 @@ int AllocatedSize;
Matrix->TrashCan.pInitInfo = NULL;
#endif
/* Allocate space in memory for Diag pointer vector. */
/* 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. */
/* 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. */
/* 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. */
/* Allocate space in memory for IntToExtColMap vector. */
if (( Matrix->IntToExtColMap = ALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Allocate space in memory for IntToExtRowMap vector. */
/* Allocate space in memory for IntToExtRowMap vector. */
if (( Matrix->IntToExtRowMap = ALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Initialize MapIntToExt vectors. */
/* Initialize MapIntToExt vectors. */
for (I = 1; I <= AllocatedSize; I++)
{ Matrix->IntToExtRowMap[I] = I;
{
Matrix->IntToExtRowMap[I] = I;
Matrix->IntToExtColMap[I] = I;
}
#if TRANSLATE
/* Allocate space in memory for ExtToIntColMap vector. */
/* Allocate space in memory for ExtToIntColMap vector. */
if (( Matrix->ExtToIntColMap = ALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Allocate space in memory for ExtToIntRowMap vector. */
/* 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;
/* 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. */
/* Allocate space for fill-ins and initial set of elements. */
InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize,
SPACE_FOR_FILL_INS*AllocatedSize );
SPACE_FOR_FILL_INS*AllocatedSize );
if (Matrix->Error == spNO_MEMORY)
goto MemoryError;
return (char *)Matrix;
return (void *)Matrix;
MemoryError:
MemoryError:
/* Deallocate matrix and return no pointer to matrix if there is not enough
memory. */
/* Deallocate matrix and return no pointer to matrix if there is not enough
memory. */
*pError = spNO_MEMORY;
spDestroy( (char *)Matrix);
spDestroy( (void *)Matrix);
return NULL;
}
@ -320,48 +305,46 @@ MemoryError:
*/
ElementPtr
spcGetElement( Matrix )
MatrixPtr Matrix;
spcGetElement(MatrixPtr Matrix)
{
ElementPtr pElements;
ElementPtr pElements;
/* Begin `spcGetElement'. */
/* Begin `spcGetElement'. */
#if NOT COMBINE OR STRIP OR LINT
/* Allocate block of MatrixElements if necessary. */
if (Matrix->ElementsRemaining == 0)
{ pElements = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, (char *)pElements );
#if !COMBINE || STRIP || LINT
/* Allocate block of MatrixElements if necessary. */
if (Matrix->ElementsRemaining == 0) {
pElements = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, (void *)pElements );
if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailElement = pElements;
}
#endif
#if COMBINE OR STRIP OR LINT
#if COMBINE || STRIP || LINT
if (Matrix->ElementsRemaining == 0)
{ pListNode = Matrix->LastElementListNode;
{
pListNode = Matrix->LastElementListNode;
/* First see if there are any stripped elements left. */
if (pListNode->Next != NULL)
{ Matrix->LastElementListNode = pListNode = pListNode->Next;
/* 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. */
} else {
/* Allocate block of elements. */
pElements = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, (char *)pElements );
RecordAllocation( Matrix, (void *)pElements );
if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailElement = pElements;
/* Allocate an element list structure. */
/* Allocate an element list structure. */
pListNode->Next = ALLOC(struct ElementListNodeStruct,1);
RecordAllocation( Matrix, (char *)pListNode->Next );
if (Matrix->Error == spNO_MEMORY) return NULL;
RecordAllocation( Matrix, (void *)pListNode->Next );
if (Matrix->Error == spNO_MEMORY)
return NULL;
Matrix->LastElementListNode = pListNode = pListNode->Next;
pListNode->pElementList = pElements;
@ -371,7 +354,7 @@ ElementPtr pElements;
}
#endif
/* Update Element counter and return pointer to Element. */
/* Update Element counter and return pointer to Element. */
Matrix->ElementsRemaining--;
return Matrix->NextAvailElement++;
@ -387,11 +370,11 @@ ElementPtr pElements;
/*
* 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.
* 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)
@ -411,30 +394,26 @@ ElementPtr pElements;
* A pointer to the first element in the group of elements being allocated.
*
* >>> Possible errors:
* spNO_MEMORY
*/
* spNO_MEMORY */
static void
InitializeElementBlocks( Matrix, InitialNumberOfElements,
NumberOfFillinsExpected )
MatrixPtr Matrix;
int InitialNumberOfElements, NumberOfFillinsExpected;
InitializeElementBlocks(MatrixPtr Matrix, int InitialNumberOfElements,
int NumberOfFillinsExpected)
{
ElementPtr pElement;
ElementPtr pElement;
/* Begin `InitializeElementBlocks'. */
/* Begin `InitializeElementBlocks'. */
/* Allocate block of MatrixElements for elements. */
/* Allocate block of MatrixElements for elements. */
pElement = ALLOC(struct MatrixElement, InitialNumberOfElements);
RecordAllocation( Matrix, (char *)pElement );
RecordAllocation( Matrix, (void *)pElement );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->ElementsRemaining = InitialNumberOfElements;
Matrix->NextAvailElement = pElement;
/* Allocate an element list structure. */
/* Allocate an element list structure. */
Matrix->FirstElementListNode = ALLOC(struct ElementListNodeStruct,1);
RecordAllocation( Matrix, (char *)Matrix->FirstElementListNode );
RecordAllocation( Matrix, (void *)Matrix->FirstElementListNode );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->LastElementListNode = Matrix->FirstElementListNode;
@ -443,16 +422,16 @@ ElementPtr pElement;
InitialNumberOfElements;
Matrix->FirstElementListNode->Next = NULL;
/* Allocate block of MatrixElements for fill-ins. */
/* Allocate block of MatrixElements for fill-ins. */
pElement = ALLOC(struct MatrixElement, NumberOfFillinsExpected);
RecordAllocation( Matrix, (char *)pElement );
RecordAllocation( Matrix, (void *)pElement );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->FillinsRemaining = NumberOfFillinsExpected;
Matrix->NextAvailFillin = pElement;
/* Allocate a fill-in list structure. */
/* Allocate a fill-in list structure. */
Matrix->FirstFillinListNode = ALLOC(struct FillinListNodeStruct,1);
RecordAllocation( Matrix, (char *)Matrix->FirstFillinListNode );
RecordAllocation( Matrix, (void *)Matrix->FirstFillinListNode );
if (Matrix->Error == spNO_MEMORY) return;
Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
@ -475,10 +454,10 @@ ElementPtr pElement;
/*
* 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.
* 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.
@ -488,43 +467,38 @@ ElementPtr pElement;
* Pointer to matrix.
*
* >>> Possible errors:
* spNO_MEMORY
*/
* spNO_MEMORY */
ElementPtr
spcGetFillin( Matrix )
MatrixPtr Matrix;
spcGetFillin(MatrixPtr Matrix)
{
/* Begin `spcGetFillin'. */
/* Begin `spcGetFillin'. */
#if NOT STRIP OR LINT
#if !STRIP || LINT
if (Matrix->FillinsRemaining == 0)
return spcGetElement( Matrix );
#endif
#if STRIP OR LINT
#if STRIP || LINT
if (Matrix->FillinsRemaining == 0)
{ pListNode = Matrix->LastFillinListNode;
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;
/* 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. */
} else {
/* Allocate block of fill-ins. */
pFillins = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, (char *)pFillins );
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. */
/* Allocate a fill-in list structure. */
pListNode->Next = ALLOC(struct FillinListNodeStruct,1);
RecordAllocation( Matrix, (char *)pListNode->Next );
RecordAllocation( Matrix, (void *)pListNode->Next );
if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->LastFillinListNode = pListNode = pListNode->Next;
@ -535,7 +509,7 @@ MatrixPtr Matrix;
}
#endif
/* Update Fill-in counter and return pointer to Fill-in. */
/* Update Fill-in counter and return pointer to Fill-in. */
Matrix->FillinsRemaining--;
return Matrix->NextAvailFillin++;
}
@ -551,50 +525,43 @@ MatrixPtr Matrix;
/*
* RECORD A MEMORY ALLOCATION
*
* This routine is used to record all memory allocations so that the memory
* can be freed later.
* 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> (char *)
* The pointer returned by malloc or calloc. These pointers are saved in
* a list so that they can be easily freed.
* 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
*/
* spNO_MEMORY */
static void
RecordAllocation( Matrix, AllocatedPtr )
MatrixPtr Matrix;
char *AllocatedPtr;
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;
/* 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);
/* 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. */
/* Add Allocated pointer to Allocation List. */
(++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr;
Matrix->RecordsRemaining--;
return;
}
@ -615,42 +582,41 @@ char *AllocatedPtr;
*
* >>> 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.
* 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
*/
* spNO_MEMORY */
static void
AllocateBlockOfAllocationList( Matrix )
MatrixPtr Matrix;
AllocateBlockOfAllocationList(MatrixPtr Matrix)
{
register int I;
register AllocationListPtr ListPtr;
int I;
AllocationListPtr ListPtr;
/* Begin `AllocateBlockOfAllocationList'. */
/* Allocate block of records for allocation list. */
/* Begin `AllocateBlockOfAllocationList'. */
/* Allocate block of records for allocation list. */
ListPtr = ALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION+1));
if (ListPtr == NULL)
{ Matrix->Error = spNO_MEMORY;
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. */
/* 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--;
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 = (char *)ListPtr;
/* Record allocation of space for allocation list on allocation list. */
Matrix->TopOfAllocationList->AllocatedPtr = (void *)ListPtr;
Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;
return;
@ -669,7 +635,7 @@ register AllocationListPtr ListPtr;
* Deallocates pointers and elements of Matrix.
*
* >>> Arguments:
* Matrix <input> (char *)
* Matrix <input> (void *)
* Pointer to the matrix frame which is to be removed from memory.
*
* >>> Local variables:
@ -684,18 +650,16 @@ register AllocationListPtr ListPtr;
*/
void
spDestroy( eMatrix )
register char *eMatrix;
spDestroy(void *eMatrix)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register AllocationListPtr ListPtr, NextListPtr;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
AllocationListPtr ListPtr, NextListPtr;
/* Begin `spDestroy'. */
ASSERT( IS_SPARSE( Matrix ) );
/* Begin `spDestroy'. */
assert( IS_SPARSE( Matrix ) );
/* Deallocate the vectors that are located in the matrix frame. */
/* Deallocate the vectors that are located in the matrix frame. */
FREE( Matrix->IntToExtColMap );
FREE( Matrix->IntToExtRowMap );
FREE( Matrix->ExtToIntColMap );
@ -710,17 +674,14 @@ register AllocationListPtr ListPtr, NextListPtr;
FREE( Matrix->DoRealDirect );
FREE( Matrix->Intermediate );
/* Sequentially step through the list of allocated pointers freeing pointers
* along the way. */
/* Sequentially step through the list of allocated pointers
* freeing pointers along the way. */
ListPtr = Matrix->TopOfAllocationList;
while (ListPtr != NULL)
{ NextListPtr = ListPtr->NextRecord;
if ((char *) ListPtr == ListPtr->AllocatedPtr)
{
while (ListPtr != NULL) {
NextListPtr = ListPtr->NextRecord;
if ((void *) ListPtr == ListPtr->AllocatedPtr) {
FREE( ListPtr );
}
else
{
} else {
FREE( ListPtr->AllocatedPtr );
}
ListPtr = NextListPtr;
@ -737,29 +698,27 @@ register AllocationListPtr ListPtr, NextListPtr;
/*
* RETURN MATRIX ERROR STATUS
*
* This function is used to determine the error status of the given matrix.
* This function is used to determine the error status of the given
* matrix.
*
* >>> Returned:
* The error status of the given matrix.
*
* >>> Arguments:
* eMatrix <input> (char *)
* The matrix for which the error status is desired.
*/
* eMatrix <input> (void *)
* The matrix for which the error status is desired. */
int
spError( eMatrix )
char *eMatrix;
spError(void *eMatrix )
{
/* Begin `spError'. */
/* Begin `spError'. */
if (eMatrix != NULL)
{ ASSERT(((MatrixPtr)eMatrix)->ID == SPARSE_ID);
if (eMatrix != NULL) {
assert(((MatrixPtr)eMatrix)->ID == SPARSE_ID);
return ((MatrixPtr)eMatrix)->Error;
} else {
/* This error may actually be spPANIC, no way to tell. */
return spNO_MEMORY;
}
else return spNO_MEMORY; /* This error may actually be spPANIC,
* no way to tell. */
}
@ -777,7 +736,7 @@ char *eMatrix;
* detected as singular or where a zero was detected on the diagonal.
*
* >>> Arguments:
* eMatrix <input> (char *)
* eMatrix <input> (void *)
* The matrix for which the error status is desired.
* pRow <output> (int *)
* The row number.
@ -786,18 +745,16 @@ char *eMatrix;
*/
void
spWhereSingular( eMatrix, pRow, pCol )
char *eMatrix;
int *pRow, *pCol;
spWhereSingular(void *eMatrix, int *pRow, int *pCol)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
/* Begin `spWhereSingular'. */
ASSERT( IS_SPARSE( Matrix ) );
/* Begin `spWhereSingular'. */
assert( IS_SPARSE( Matrix ) );
if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG)
{ *pRow = Matrix->SingularRow;
if (Matrix->Error == spSINGULAR || Matrix->Error == spZERO_DIAG)
{
*pRow = Matrix->SingularRow;
*pCol = Matrix->SingularCol;
}
else *pRow = *pCol = 0;
@ -816,9 +773,9 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix;
* the matrix is returned.
*
* >>> Arguments:
* eMatrix <input> (char *)
* eMatrix <input> (void *)
* Pointer to matrix.
* External <input> (BOOLEAN)
* 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
@ -826,15 +783,12 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix;
*/
int
spGetSize( eMatrix, External )
char *eMatrix;
BOOLEAN External;
spGetSize(void *eMatrix, int External)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
/* Begin `spGetSize'. */
ASSERT( IS_SPARSE( Matrix ) );
/* Begin `spGetSize'. */
assert( IS_SPARSE( Matrix ) );
#if TRANSLATE
if (External)
@ -859,31 +813,27 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix;
* Forces matrix to be either real or complex.
*
* >>> Arguments:
* eMatrix <input> (char *)
* eMatrix <input> (void *)
* Pointer to matrix.
*/
void
spSetReal( eMatrix )
char *eMatrix;
spSetReal(void *eMatrix)
{
/* Begin `spSetReal'. */
/* Begin `spSetReal'. */
ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) AND REAL);
assert( IS_SPARSE( (MatrixPtr)eMatrix ));
((MatrixPtr)eMatrix)->Complex = NO;
return;
}
void
spSetComplex( eMatrix )
char *eMatrix;
spSetComplex(void *eMatrix)
{
/* Begin `spSetComplex'. */
/* Begin `spSetComplex'. */
ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) AND spCOMPLEX);
assert( IS_SPARSE( (MatrixPtr)eMatrix ));
((MatrixPtr)eMatrix)->Complex = YES;
return;
}
@ -904,40 +854,34 @@ char *eMatrix;
* of original elements can be returned.
*
* >>> Arguments:
* eMatrix <input> (char *)
* eMatrix <input> (void *)
* Pointer to matrix.
*/
int
spFillinCount( eMatrix )
char *eMatrix;
spFillinCount(void *eMatrix)
{
/* Begin `spFillinCount'. */
/* Begin `spFillinCount'. */
ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) );
assert( IS_SPARSE( (MatrixPtr)eMatrix ) );
return ((MatrixPtr)eMatrix)->Fillins;
}
int
spElementCount( eMatrix )
char *eMatrix;
spElementCount(void *eMatrix)
{
/* Begin `spElementCount'. */
/* Begin `spElementCount'. */
ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) );
assert( IS_SPARSE( (MatrixPtr)eMatrix ) );
return ((MatrixPtr)eMatrix)->Elements;
}
int
spOriginalCount( eMatrix )
char *eMatrix;
spOriginalCount(void *eMatrix)
{
/* Begin `spOriginalCount'. */
/* Begin `spOriginalCount'. */
ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) );
assert( IS_SPARSE( (MatrixPtr)eMatrix ) );
return ((MatrixPtr)eMatrix)->Originals;
}

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

@ -56,6 +56,7 @@
* spDefs.h
* Matrix type and macro definitions for the sparse matrix routines.
*/
#include <assert.h>
#define spINSIDE_SPARSE
#include "spconfig.h"
@ -69,17 +70,9 @@
/*
* Function declarations
*/
#ifdef __STDC__
static void Translate( MatrixPtr, int*, int* );
static void EnlargeMatrix( MatrixPtr, int );
static void ExpandTranslationArrays( MatrixPtr, int );
#else /* __STDC__ */
static void Translate();
static void EnlargeMatrix();
static void ExpandTranslationArrays();
#endif /* __STDC__ */
@ -100,45 +93,45 @@ static void ExpandTranslationArrays();
*/
void
spClear( eMatrix )
char *eMatrix;
spClear(void *eMatrix )
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register ElementPtr pElement;
register int I;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
ElementPtr pElement;
int I;
/* Begin `spClear'. */
ASSERT( IS_SPARSE( Matrix ) );
/* Begin `spClear'. */
assert( IS_SPARSE( Matrix ) );
/* Clear matrix. */
#if spCOMPLEX
if (Matrix->PreviousMatrixWasComplex OR Matrix->Complex)
{ for (I = Matrix->Size; I > 0; I--)
{ pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{ pElement->Real = 0.0;
pElement->Imag = 0.0;
pElement = pElement->NextInCol;
}
}
/* Clear matrix. */
if (Matrix->PreviousMatrixWasComplex || Matrix->Complex)
{
for (I = Matrix->Size; I > 0; I--)
{
pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{
pElement->Real = 0.0;
pElement->Imag = 0.0;
pElement = pElement->NextInCol;
}
}
}
else
#endif
{ for (I = Matrix->Size; I > 0; I--)
{ pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{ pElement->Real = 0.0;
pElement = pElement->NextInCol;
}
}
{
for (I = Matrix->Size; I > 0; I--)
{
pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{
pElement->Real = 0.0;
pElement = pElement->NextInCol;
}
}
}
/* Empty the trash. */
/* Empty the trash. */
Matrix->TrashCan.Real = 0.0;
#if spCOMPLEX
Matrix->TrashCan.Imag = 0.0;
#endif
Matrix->Error = spOKAY;
Matrix->Factored = NO;
@ -194,24 +187,21 @@ register int I;
*/
RealNumber *
spGetElement( eMatrix, Row, Col )
char *eMatrix;
int Row, Col;
spGetElement(void *eMatrix, int Row, int Col)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
RealNumber *pElement;
ElementPtr spcFindElementInCol();
void Translate();
MatrixPtr Matrix = (MatrixPtr)eMatrix;
RealNumber *pElement;
ElementPtr spcFindElementInCol();
void Translate();
/* Begin `spGetElement'. */
ASSERT( IS_SPARSE( Matrix ) AND Row >= 0 AND Col >= 0 );
/* Begin `spGetElement'. */
assert( IS_SPARSE( Matrix ) && Row >= 0 && Col >= 0 );
if ((Row == 0) OR (Col == 0))
return &Matrix->TrashCan.Real;
if ((Row == 0) || (Col == 0))
return &Matrix->TrashCan.Real;
#if NOT TRANSLATE
ASSERT(Matrix->NeedsOrdering);
#if !TRANSLATE
assert(Matrix->NeedsOrdering);
#endif
#if TRANSLATE
@ -219,40 +209,38 @@ void Translate();
if (Matrix->Error == spNO_MEMORY) return NULL;
#endif
#if NOT TRANSLATE
#if NOT EXPANDABLE
ASSERT(Row <= Matrix->Size AND Col <= Matrix->Size);
#if !TRANSLATE
#if !EXPANDABLE
assert(Row <= Matrix->Size && Col <= Matrix->Size);
#endif
#if EXPANDABLE
/* Re-size Matrix if necessary. */
if ((Row > Matrix->Size) OR (Col > Matrix->Size))
EnlargeMatrix( Matrix, MAX(Row, Col) );
/* Re-size Matrix if necessary. */
if ((Row > Matrix->Size) || (Col > Matrix->Size))
EnlargeMatrix( Matrix, MAX(Row, Col) );
if (Matrix->Error == spNO_MEMORY) return NULL;
#endif
#endif
/*
* The condition part of the following if statement tests to see if the
* element resides along the diagonal, if it does then it tests to see
* if the element has been created yet (Diag pointer not NULL). The
* pointer to the element is then assigned to Element after it is cast
* into a pointer to a RealNumber. This casting makes the pointer into
* a pointer to Real. This statement depends on the fact that Real
* is the first record in the MatrixElement structure.
*/
/* The condition part of the following if statement tests to see
* if the element resides along the diagonal, if it does then it
* tests to see if the element has been created yet (Diag pointer
* not NULL). The pointer to the element is then assigned to
* Element after it is cast into a pointer to a RealNumber. This
* casting makes the pointer into a pointer to Real. This
* statement depends on the fact that Real is the first record in
* the MatrixElement structure. */
if ((Row != Col) OR ((pElement = (RealNumber *)Matrix->Diag[Row]) == NULL))
if ((Row != Col) || ((pElement = (RealNumber *)Matrix->Diag[Row]) == NULL))
{
/*
* Element does not exist or does not reside along diagonal. Search
* column for element. As in the if statement above, the pointer to the
* element which is returned by spcFindElementInCol is cast into a
* pointer to Real, a RealNumber.
*/
pElement = (RealNumber*)spcFindElementInCol( Matrix,
&(Matrix->FirstInCol[Col]),
Row, Col, YES );
/* Element does not exist or does not reside along diagonal.
* Search column for element. As in the if statement above,
* the pointer to the element which is returned by
* spcFindElementInCol is cast into a pointer to Real, a
* RealNumber. */
pElement = (RealNumber*)spcFindElementInCol( Matrix,
&(Matrix->FirstInCol[Col]),
Row, Col, YES );
}
return pElement;
}
@ -290,7 +278,7 @@ void Translate();
* Row being searched for.
* Col (int)
* Column being searched.
* CreateIfMissing <input> (BOOLEAN)
* CreateIfMissing <input> (int)
* Indicates what to do if element is not found, create one or return a
* NULL pointer.
*
@ -300,40 +288,37 @@ void Translate();
*/
ElementPtr
spcFindElementInCol( Matrix, LastAddr, Row, Col, CreateIfMissing )
MatrixPtr Matrix;
register ElementPtr *LastAddr;
register int Row;
int Col;
BOOLEAN CreateIfMissing;
spcFindElementInCol(MatrixPtr Matrix, ElementPtr *LastAddr,
int Row, int Col, int CreateIfMissing)
{
register ElementPtr pElement;
ElementPtr spcCreateElement();
ElementPtr pElement;
ElementPtr spcCreateElement();
/* Begin `spcFindElementInCol'. */
/* Begin `spcFindElementInCol'. */
pElement = *LastAddr;
/* Search for element. */
/* Search for element. */
while (pElement != NULL)
{ if (pElement->Row < Row)
{
if (pElement->Row < Row)
{
/* Have not reached element yet. */
LastAddr = &(pElement->NextInCol);
pElement = pElement->NextInCol;
/* Have not reached element yet. */
LastAddr = &(pElement->NextInCol);
pElement = pElement->NextInCol;
}
else if (pElement->Row == Row)
else if (pElement->Row == Row)
{
/* Reached element. */
return pElement;
/* Reached element. */
return pElement;
}
else break; /* while loop */
else break; /* while loop */
}
/* Element does not exist and must be created. */
/* Element does not exist and must be created. */
if (CreateIfMissing)
return spcCreateElement( Matrix, Row, Col, LastAddr, NO );
else return NULL;
return spcCreateElement( Matrix, Row, Col, LastAddr, NO );
else
return NULL;
}
@ -378,41 +363,39 @@ ElementPtr spcCreateElement();
*/
static void
Translate( Matrix, Row, Col )
MatrixPtr Matrix;
int *Row, *Col;
Translate(MatrixPtr Matrix, int *Row, int *Col)
{
register int IntRow, IntCol, ExtRow, ExtCol;
int IntRow, IntCol, ExtRow, ExtCol;
/* Begin `Translate'. */
/* Begin `Translate'. */
ExtRow = *Row;
ExtCol = *Col;
/* Expand translation arrays if necessary. */
if ((ExtRow > Matrix->AllocatedExtSize) OR
/* Expand translation arrays if necessary. */
if ((ExtRow > Matrix->AllocatedExtSize) ||
(ExtCol > Matrix->AllocatedExtSize))
{
ExpandTranslationArrays( Matrix, MAX(ExtRow, ExtCol) );
if (Matrix->Error == spNO_MEMORY) return;
}
/* Set ExtSize if necessary. */
if ((ExtRow > Matrix->ExtSize) OR (ExtCol > Matrix->ExtSize))
/* Set ExtSize if necessary. */
if ((ExtRow > Matrix->ExtSize) || (ExtCol > Matrix->ExtSize))
Matrix->ExtSize = MAX(ExtRow, ExtCol);
/* Translate external row or node number to internal row or node number. */
/* Translate external row or node number to internal row or node number. */
if ((IntRow = Matrix->ExtToIntRowMap[ExtRow]) == -1)
{ Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize;
{
Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize;
Matrix->ExtToIntColMap[ExtRow] = Matrix->CurrentSize;
IntRow = Matrix->CurrentSize;
#if NOT EXPANDABLE
ASSERT(IntRow <= Matrix->Size);
#if !EXPANDABLE
assert(IntRow <= Matrix->Size);
#endif
#if EXPANDABLE
/* Re-size Matrix if necessary. */
/* Re-size Matrix if necessary. */
if (IntRow > Matrix->Size)
EnlargeMatrix( Matrix, IntRow );
if (Matrix->Error == spNO_MEMORY) return;
@ -422,18 +405,19 @@ register int IntRow, IntCol, ExtRow, ExtCol;
Matrix->IntToExtColMap[IntRow] = ExtRow;
}
/* Translate external column or node number to internal column or node number.*/
/* Translate external column or node number to internal column or node number.*/
if ((IntCol = Matrix->ExtToIntColMap[ExtCol]) == -1)
{ Matrix->ExtToIntRowMap[ExtCol] = ++Matrix->CurrentSize;
{
Matrix->ExtToIntRowMap[ExtCol] = ++Matrix->CurrentSize;
Matrix->ExtToIntColMap[ExtCol] = Matrix->CurrentSize;
IntCol = Matrix->CurrentSize;
#if NOT EXPANDABLE
ASSERT(IntCol <= Matrix->Size);
#if !EXPANDABLE
assert(IntCol <= Matrix->Size);
#endif
#if EXPANDABLE
/* Re-size Matrix if necessary. */
/* Re-size Matrix if necessary. */
if (IntCol > Matrix->Size)
EnlargeMatrix( Matrix, IntCol );
if (Matrix->Error == spNO_MEMORY) return;
@ -490,24 +474,19 @@ register int IntRow, IntCol, ExtRow, ExtCol;
*/
int
spGetAdmittance( Matrix, Node1, Node2, Template )
char *Matrix;
int Node1, Node2;
struct spTemplate *Template;
spGetAdmittance(void *Matrix, int Node1, int Node2,
struct spTemplate *Template)
{
/* Begin `spGetAdmittance'. */
/* Begin `spGetAdmittance'. */
Template->Element1 = spGetElement(Matrix, Node1, Node1 );
Template->Element2 = spGetElement(Matrix, Node2, Node2 );
Template->Element3Negated = spGetElement( Matrix, Node2, Node1 );
Template->Element4Negated = spGetElement( Matrix, Node1, Node2 );
if
( (Template->Element1 == NULL)
OR (Template->Element2 == NULL)
OR (Template->Element3Negated == NULL)
OR (Template->Element4Negated == NULL)
) return spNO_MEMORY;
if ((Template->Element1 == NULL) ||
(Template->Element2 == NULL) ||
(Template->Element3Negated == NULL) ||
(Template->Element4Negated == NULL))
return spNO_MEMORY;
if (Node1 == 0)
SWAP( RealNumber*, Template->Element1, Template->Element2 );
@ -578,23 +557,19 @@ struct spTemplate *Template;
*/
int
spGetQuad( Matrix, Row1, Row2, Col1, Col2, Template )
char *Matrix;
int Row1, Row2, Col1, Col2;
struct spTemplate *Template;
spGetQuad(void *Matrix, int Row1, int Row2, int Col1, int Col2,
struct spTemplate *Template)
{
/* Begin `spGetQuad'. */
/* Begin `spGetQuad'. */
Template->Element1 = spGetElement( Matrix, Row1, Col1);
Template->Element2 = spGetElement( Matrix, Row2, Col2 );
Template->Element3Negated = spGetElement( Matrix, Row2, Col1 );
Template->Element4Negated = spGetElement( Matrix, Row1, Col2 );
if
( (Template->Element1 == NULL)
OR (Template->Element2 == NULL)
OR (Template->Element3Negated == NULL)
OR (Template->Element4Negated == NULL)
) return spNO_MEMORY;
if ((Template->Element1 == NULL) ||
(Template->Element2 == NULL) ||
(Template->Element3Negated == NULL) ||
(Template->Element4Negated == NULL))
return spNO_MEMORY;
if (Template->Element1 == &((MatrixPtr)Matrix)->TrashCan.Real)
SWAP( RealNumber *, Template->Element1, Template->Element2 );
@ -653,23 +628,19 @@ struct spTemplate *Template;
*/
int
spGetOnes(Matrix, Pos, Neg, Eqn, Template)
char *Matrix;
int Pos, Neg, Eqn;
struct spTemplate *Template;
spGetOnes(void *Matrix, int Pos, int Neg, int Eqn,
struct spTemplate *Template)
{
/* Begin `spGetOnes'. */
/* Begin `spGetOnes'. */
Template->Element4Negated = spGetElement( Matrix, Neg, Eqn );
Template->Element3Negated = spGetElement( Matrix, Eqn, Neg );
Template->Element2 = spGetElement( Matrix, Pos, Eqn );
Template->Element1 = spGetElement( Matrix, Eqn, Pos );
if
( (Template->Element1 == NULL)
OR (Template->Element2 == NULL)
OR (Template->Element3Negated == NULL)
OR (Template->Element4Negated == NULL)
) return spNO_MEMORY;
if ((Template->Element1 == NULL) ||
(Template->Element2 == NULL) ||
(Template->Element3Negated == NULL) ||
(Template->Element4Negated == NULL))
return spNO_MEMORY;
spADD_REAL_QUAD( *Template, 1.0 );
return spOKAY;
@ -703,7 +674,7 @@ struct spTemplate *Template;
* This contains the address of the pointer to the element just above the
* one being created. It is used to speed the search and it is updated with
* address of the created element.
* Fillin <input> (BOOLEAN)
* Fillin <input> (int)
* Flag that indicates if created element is to be a fill-in.
*
* >>> Local variables:
@ -723,78 +694,73 @@ struct spTemplate *Template;
*/
ElementPtr
spcCreateElement( Matrix, Row, Col, LastAddr, Fillin )
MatrixPtr Matrix;
int Row;
register int Col;
register ElementPtr *LastAddr;
BOOLEAN Fillin;
spcCreateElement(MatrixPtr Matrix, int Row, int Col,
ElementPtr *LastAddr, int Fillin)
{
register ElementPtr pElement, pLastElement;
ElementPtr pCreatedElement, spcGetElement(), spcGetFillin();
ElementPtr pElement, pLastElement;
ElementPtr pCreatedElement, spcGetElement(), spcGetFillin();
/* Begin `spcCreateElement'. */
/* Begin `spcCreateElement'. */
if (Matrix->RowsLinked)
{
/* Row pointers cannot be ignored. */
/* Row pointers cannot be ignored. */
if (Fillin)
{ pElement = spcGetFillin( Matrix );
{
pElement = spcGetFillin( Matrix );
Matrix->Fillins++;
}
else
{ pElement = spcGetElement( Matrix );
{
pElement = spcGetElement( Matrix );
Matrix->Originals++;
Matrix->NeedsOrdering = YES;
}
if (pElement == NULL) return NULL;
/* If element is on diagonal, store pointer in Diag. */
/* If element is on diagonal, store pointer in Diag. */
if (Row == Col) Matrix->Diag[Row] = pElement;
/* Initialize Element. */
/* Initialize Element. */
pCreatedElement = pElement;
pElement->Row = Row;
pElement->Col = Col;
pElement->Real = 0.0;
#if spCOMPLEX
pElement->Imag = 0.0;
#endif
#if INITIALIZE
pElement->pInitInfo = NULL;
#endif
/* Splice element into column. */
/* Splice element into column. */
pElement->NextInCol = *LastAddr;
*LastAddr = pElement;
/* Search row for proper element position. */
/* Search row for proper element position. */
pElement = Matrix->FirstInRow[Row];
pLastElement = NULL;
while (pElement != NULL)
{
/* Search for element row position. */
/* Search for element row position. */
if (pElement->Col < Col)
{
/* Have not reached desired element. */
/* Have not reached desired element. */
pLastElement = pElement;
pElement = pElement->NextInRow;
}
else pElement = NULL;
}
/* Splice element into row. */
/* Splice element into row. */
pElement = pCreatedElement;
if (pLastElement == NULL)
{
/* Element is first in row. */
/* Element is first in row. */
pElement->NextInRow = Matrix->FirstInRow[Row];
Matrix->FirstInRow[Row] = pElement;
}
else
/* Element is not first in row. */
{
/* Element is not first in row. */
pElement->NextInRow = pLastElement->NextInRow;
pLastElement->NextInRow = pElement;
}
@ -802,34 +768,30 @@ ElementPtr pCreatedElement, spcGetElement(), spcGetFillin();
}
else
{
/*
* Matrix has not been factored yet. Thus get element rather than fill-in.
* Also, row pointers can be ignored.
*/
/* Matrix has not been factored yet. Thus get element rather
* than fill-in. Also, row pointers can be ignored. */
/* Allocate memory for Element. */
/* Allocate memory for Element. */
pElement = spcGetElement( Matrix );
Matrix->Originals++;
if (pElement == NULL) return NULL;
/* If element is on diagonal, store pointer in Diag. */
/* If element is on diagonal, store pointer in Diag. */
if (Row == Col) Matrix->Diag[Row] = pElement;
/* Initialize Element. */
/* Initialize Element. */
pCreatedElement = pElement;
pElement->Row = Row;
#if DEBUG
pElement->Col = Col;
#endif
pElement->Real = 0.0;
#if spCOMPLEX
pElement->Imag = 0.0;
#endif
#if INITIALIZE
pElement->pInitInfo = NULL;
#endif
/* Splice element into column. */
/* Splice element into column. */
pElement->NextInCol = *LastAddr;
*LastAddr = pElement;
}
@ -865,30 +827,29 @@ ElementPtr pCreatedElement, spcGetElement(), spcGetFillin();
* currently being operated upon.
* FirstInRowArray (ArrayOfElementPtrs)
* A pointer to the FirstInRow array. Same as Matrix->FirstInRow but
* resides in a register and requires less indirection so is faster to
* resides in a and requires less indirection so is faster to
* use.
* Col (int)
* Column currently being operated upon.
*/
void
spcLinkRows( Matrix )
MatrixPtr Matrix;
spcLinkRows(MatrixPtr Matrix)
{
register ElementPtr pElement, *FirstInRowEntry;
register ArrayOfElementPtrs FirstInRowArray;
register int Col;
ElementPtr pElement, *FirstInRowEntry;
ArrayOfElementPtrs FirstInRowArray;
int Col;
/* Begin `spcLinkRows'. */
/* Begin `spcLinkRows'. */
FirstInRowArray = Matrix->FirstInRow;
for (Col = Matrix->Size; Col >= 1; Col--)
{
/* Generate row links for the elements in the Col'th column. */
/* Generate row links for the elements in the Col'th column. */
pElement = Matrix->FirstInCol[Col];
while (pElement != NULL)
{ pElement->Col = Col;
{
pElement->Col = Col;
FirstInRowEntry = &FirstInRowArray[pElement->Row];
pElement->NextInRow = *FirstInRowEntry;
*FirstInRowEntry = pElement;
@ -923,48 +884,48 @@ register int Col;
*/
static void
EnlargeMatrix( Matrix, NewSize )
MatrixPtr Matrix;
register int NewSize;
EnlargeMatrix(MatrixPtr Matrix, int NewSize)
{
register int I, OldAllocatedSize = Matrix->AllocatedSize;
int I, OldAllocatedSize = Matrix->AllocatedSize;
/* Begin `EnlargeMatrix'. */
/* Begin `EnlargeMatrix'. */
Matrix->Size = NewSize;
if (NewSize <= OldAllocatedSize)
return;
/* Expand the matrix frame. */
/* Expand the matrix frame. */
NewSize = MAX( NewSize, EXPANSION_FACTOR * OldAllocatedSize );
Matrix->AllocatedSize = NewSize;
if (( REALLOC(Matrix->IntToExtColMap, int, NewSize+1)) == NULL)
{ Matrix->Error = spNO_MEMORY;
{
Matrix->Error = spNO_MEMORY;
return;
}
if (( REALLOC(Matrix->IntToExtRowMap, int, NewSize+1)) == NULL)
{ Matrix->Error = spNO_MEMORY;
{
Matrix->Error = spNO_MEMORY;
return;
}
if (( REALLOC(Matrix->Diag, ElementPtr, NewSize+1)) == NULL)
{ Matrix->Error = spNO_MEMORY;
{
Matrix->Error = spNO_MEMORY;
return;
}
if (( REALLOC(Matrix->FirstInCol, ElementPtr, NewSize+1)) == NULL)
{ Matrix->Error = spNO_MEMORY;
{
Matrix->Error = spNO_MEMORY;
return;
}
if (( REALLOC(Matrix->FirstInRow, ElementPtr, NewSize+1)) == NULL)
{ Matrix->Error = spNO_MEMORY;
{
Matrix->Error = spNO_MEMORY;
return;
}
/*
* Destroy the Markowitz and Intermediate vectors, they will be recreated
* in spOrderAndFactor().
*/
/* Destroy the Markowitz and Intermediate vectors, they will be
* recreated in spOrderAndFactor(). */
FREE( Matrix->MarkowitzRow );
FREE( Matrix->MarkowitzCol );
FREE( Matrix->MarkowitzProd );
@ -973,9 +934,10 @@ register int I, OldAllocatedSize = Matrix->AllocatedSize;
FREE( Matrix->Intermediate );
Matrix->InternalVectorsAllocated = NO;
/* Initialize the new portion of the vectors. */
/* Initialize the new portion of the vectors. */
for (I = OldAllocatedSize+1; I <= NewSize; I++)
{ Matrix->IntToExtColMap[I] = I;
{
Matrix->IntToExtColMap[I] = I;
Matrix->IntToExtRowMap[I] = I;
Matrix->Diag[I] = NULL;
Matrix->FirstInRow[I] = NULL;
@ -1012,35 +974,35 @@ register int I, OldAllocatedSize = Matrix->AllocatedSize;
*/
static void
ExpandTranslationArrays( Matrix, NewSize )
MatrixPtr Matrix;
register int NewSize;
ExpandTranslationArrays(MatrixPtr Matrix, int NewSize)
{
register int I, OldAllocatedSize = Matrix->AllocatedExtSize;
int I, OldAllocatedSize = Matrix->AllocatedExtSize;
/* Begin `ExpandTranslationArrays'. */
/* Begin `ExpandTranslationArrays'. */
Matrix->ExtSize = NewSize;
if (NewSize <= OldAllocatedSize)
return;
/* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */
/* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */
NewSize = MAX( NewSize, EXPANSION_FACTOR * OldAllocatedSize );
Matrix->AllocatedExtSize = NewSize;
if (( REALLOC(Matrix->ExtToIntRowMap, int, NewSize+1)) == NULL)
{ Matrix->Error = spNO_MEMORY;
{
Matrix->Error = spNO_MEMORY;
return;
}
if (( REALLOC(Matrix->ExtToIntColMap, int, NewSize+1)) == NULL)
{ Matrix->Error = spNO_MEMORY;
{
Matrix->Error = spNO_MEMORY;
return;
}
/* Initialize the new portion of the vectors. */
/* Initialize the new portion of the vectors. */
for (I = OldAllocatedSize+1; I <= NewSize; I++)
{ Matrix->ExtToIntRowMap[I] = -1;
{
Matrix->ExtToIntRowMap[I] = -1;
Matrix->ExtToIntColMap[I] = -1;
}
@ -1090,72 +1052,68 @@ register int I, OldAllocatedSize = Matrix->AllocatedExtSize;
*/
void
spInstallInitInfo( pElement, pInitInfo )
RealNumber *pElement;
char *pInitInfo;
spInstallInitInfo(RealNumber *pElement, char *pInitInfo)
{
/* Begin `spInstallInitInfo'. */
ASSERT(pElement != NULL);
/* Begin `spInstallInitInfo'. */
assert(pElement != NULL);
((ElementPtr)pElement)->pInitInfo = pInitInfo;
}
char *
spGetInitInfo( pElement )
RealNumber *pElement;
spGetInitInfo(RealNumber *pElement)
{
/* Begin `spGetInitInfo'. */
ASSERT(pElement != NULL);
/* Begin `spGetInitInfo'. */
assert(pElement != NULL);
return (char *)((ElementPtr)pElement)->pInitInfo;
}
int
spInitialize( eMatrix, pInit )
char *eMatrix;
int (*pInit)();
spInitialize(void *eMatrix, int (*pInit)())
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register ElementPtr pElement;
int J, Error, Col;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
ElementPtr pElement;
int J, Error, Col;
/* Begin `spInitialize'. */
ASSERT( IS_SPARSE( Matrix ) );
/* Begin `spInitialize'. */
assert( IS_SPARSE( Matrix ) );
#if spCOMPLEX
/* Clear imaginary part of matrix if matrix is real but was complex. */
if (Matrix->PreviousMatrixWasComplex AND NOT Matrix->Complex)
{ for (J = Matrix->Size; J > 0; J--)
{ pElement = Matrix->FirstInCol[J];
/* Clear imaginary part of matrix if matrix is real but was complex. */
if (Matrix->PreviousMatrixWasComplex && !Matrix->Complex)
{
for (J = Matrix->Size; J > 0; J--)
{
pElement = Matrix->FirstInCol[J];
while (pElement != NULL)
{ pElement->Imag = 0.0;
{
pElement->Imag = 0.0;
pElement = pElement->NextInCol;
}
}
}
#endif /* spCOMPLEX */
/* Initialize the matrix. */
/* Initialize the matrix. */
for (J = Matrix->Size; J > 0; J--)
{ pElement = Matrix->FirstInCol[J];
{
pElement = Matrix->FirstInCol[J];
Col = Matrix->IntToExtColMap[J];
while (pElement != NULL)
{ if (pElement->pInitInfo == NULL)
{ pElement->Real = 0.0;
# if spCOMPLEX
pElement->Imag = 0.0;
# endif
{
if (pElement->pInitInfo == NULL)
{
pElement->Real = 0.0;
pElement->Imag = 0.0;
}
else
{ Error = (*pInit)((RealNumber *)pElement, pElement->pInitInfo,
{
Error = (*pInit)((RealNumber *)pElement, pElement->pInitInfo,
Matrix->IntToExtRowMap[pElement->Row], Col);
if (Error)
{ Matrix->Error = spFATAL;
{
Matrix->Error = spFATAL;
return Error;
}
@ -1164,11 +1122,9 @@ int J, Error, Col;
}
}
/* Empty the trash. */
/* Empty the trash. */
Matrix->TrashCan.Real = 0.0;
#if spCOMPLEX
Matrix->TrashCan.Imag = 0.0;
#endif
Matrix->Error = spOKAY;
Matrix->Factored = NO;

44
src/maths/sparse/spcombin.c
View File

@ -58,20 +58,10 @@
#if COMBINE
#ifdef __STDC__
#if spSEPARATED_COMPLEX_VECTORS
static void CombineComplexMatrix( MatrixPtr,
RealVector, RealVector, RealVector, RealVector );
#else
static void CombineComplexMatrix( MatrixPtr, RealVector, RealVector );
#endif
static void ClearBuffer( MatrixPtr, int, int, ElementPtr );
static void ClearComplexBuffer( MatrixPtr, int, int, ElementPtr );
#else /* __STDC__ */
static void CombineComplexMatrix();
static void ClearBuffer();
static void ClearComplexBuffer();
#endif /* __STDC__ */
/*
* COMBINE MATRICES ON A MULTIPROCESSOR
@ -91,10 +81,10 @@ static void ClearComplexBuffer();
static double Buffer[SPBSIZE];
void
spCombine( eMatrix, RHS, Spare IMAG_VECTORS )
spCombine( eMatrix, RHS, Spare, iRHS, iSolution )
char *eMatrix;
RealVector RHS, Spare IMAG_VECTORS;
RealVector RHS, Spare, iRHS, iSolution;
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register ElementPtr pElement;
@ -105,18 +95,15 @@ struct ElementListNodeStruct *pListNode;
long type = MT_COMBINE, length = Matrix->Size + 1;
/* Begin `spCombine'. */
ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored );
if (NOT Matrix->InternalVectorsAllocated)
assert( IS_VALID(Matrix) && !Matrix->Factored );
if (!Matrix->InternalVectorsAllocated)
spcCreateInternalVectors( Matrix );
#if spCOMPLEX
if (Matrix->Complex) {
CombineComplexMatrix( Matrix, RHS, Spare IMAG_VECTORS );
CombineComplexMatrix( Matrix, RHS, Spare, iRHS, iSolution );
return;
}
#endif
#if REAL
Size = Matrix->Size;
/* Mark original non-zeroes. */
@ -167,15 +154,14 @@ long type = MT_COMBINE, length = Matrix->Size + 1;
DGOP_( &type, RHS, &length, "+" );
return;
#endif /* REAL */
}
#if spCOMPLEX
static void
CombineComplexMatrix( Matrix, RHS, Spare IMAG_VECTORS )
CombineComplexMatrix( Matrix, RHS, Spare, iRHS, iSolution )
MatrixPtr Matrix;
RealVector RHS, Spare IMAG_VECTORS;
RealVector RHS, Spare, iRHS, iSolution;
{
register ElementPtr pElement;
register int I, Size;
@ -185,7 +171,7 @@ struct ElementListNodeStruct *pListNode;
long type = MT_COMBINE, length = Matrix->Size + 1;
/* Begin `CombineComplexMatrix'. */
ASSERT(Matrix->Complex);
assert(Matrix->Complex);
Size = Matrix->Size;
/* Mark original non-zeroes. */
@ -235,19 +221,13 @@ long type = MT_COMBINE, length = Matrix->Size + 1;
}
/* Sum all RHS's together */
#if spSEPARATED_COMPLEX_VECTORS
DGOP_( &type, RHS, &length, "+" );
DGOP_( &type, iRHS, &length, "+" );
#else
length *= 2;
DGOP_( &type, RHS, &length, "+" );
#endif
return;
}
#endif /* spCOMPLEX */
#if REAL
static void
ClearBuffer( Matrix, NumElems, StartCol, StartElement )
@ -281,9 +261,8 @@ ElementPtr StartElement;
pElement = pElement->NextInCol;
}
}
#endif REAL
#if spCOMPLEX
static void
ClearComplexBuffer( Matrix, DataCount, StartCol, StartElement )
@ -318,5 +297,4 @@ ElementPtr StartElement;
pElement = pElement->NextInCol;
}
}
#endif /* spCOMPLEX */
#endif /* COMBINE */

54
src/include/spconfig.h → src/maths/sparse/spconfig.h
View File

@ -62,15 +62,6 @@
* less storage is required. There is often a noticeable speed
* penalty when using single precision. Sparse internally refers
* to a spREAL as a RealNumber.
* REAL
* This specifies that the routines are expected to handle real
* systems of equations. The routines can be compiled to handle
* both real and complex systems at the same time, but there is a
* slight speed and memory advantage if the routines are complied
* to handle only real systems of equations.
* spCOMPLEX
* This specifies that the routines will be complied to handle
* complex systems of equations.
* EXPANDABLE
* Setting this compiler flag true (1) makes the matrix
* expandable before it has been factored. If the matrix is
@ -129,28 +120,6 @@
* threshold. When DIAGONAL_PIVOTING is turned off, the following
* options and constants are not used: MODIFIED_MARKOWITZ,
* MAX_MARKOWITZ_TIES, and TIES_MULTIPLIER.
* ARRAY_OFFSET
* This determines whether arrays start at an index of zero or one.
* This option is necessitated by the fact that standard C
* convention dictates that arrays begin with an index of zero but
* the standard mathematic convention states that arrays begin with
* an index of one. So if you prefer to start your arrays with
* zero, or your calling Sparse from FORTRAN, set ARRAY_OFFSET to
* NO or 0. Otherwise, set ARRAY_OFFSET to YES or 1. Note that if
* you use an offset of one, the arrays that you pass to Sparse
* must have an allocated length of one plus the size of the
* matrix. ARRAY_OFFSET must be either 0 or 1, no other offsets
* are valid.
* spSEPARATED_COMPLEX_VECTORS
* This specifies the format for complex vectors. If this is set
* false then a complex vector is made up of one double sized
* array of RealNumber's in which the real and imaginary numbers
* are placed in the alternately array in the array. In other
* words, the first entry would be Complex[1].Real, then comes
* Complex[1].Imag, then Complex[1].Real, etc. If
* spSEPARATED_COMPLEX_VECTORS is set true, then each complex
* vector is represented by two arrays of RealNumbers, one with
* the real terms, the other with the imaginary. [NO]
* MODIFIED_MARKOWITZ
* This specifies that the modified Markowitz method of pivot
* selection is to be used. The modified Markowitz method differs
@ -229,33 +198,19 @@
* MULTIPLICATION
* This specifies that the routines to multiply the unfactored
* matrix by a vector should be compiled.
* FORTRAN
* This specifies that the FORTRAN interface routines should be
* compiled. When interfacing to FORTRAN programs, the ARRAY_OFFSET
* options should be set to NO.
* DEBUG
* This specifies that additional error checking will be compiled.
* The type of error checked are those that are common when the
* matrix routines are first integrated into a user's program. Once
* the routines have been integrated in and are running smoothly, this
* option should be turned off.
* spCOMPATIBILITY
* This specifies that Sparse1.3 should be configured to be upward
* compatible from Sparse1.2. This option is not suggested for use
* in new software. Sparse1.3, when configured to be compatible with
* Sparse1.2, is not completely compatible. In particular, if
* recompiling the calling program, it is necessary to change the
* of the Sparse include files. This option will go away in future
* versions of Sparse. [0]
*/
/* Begin options. */
#define REAL YES
#define EXPANDABLE YES
#define TRANSLATE YES
#define INITIALIZE NO
#define DIAGONAL_PIVOTING YES
#define ARRAY_OFFSET YES
#define MODIFIED_MARKOWITZ NO
#define DELETE NO
#define STRIP NO
@ -270,7 +225,6 @@
#define STABILITY NO
#define CONDITION NO
#define PSEUDOCONDITION NO
#define FORTRAN NO
#ifdef HAS_MINDATA
# define DEBUG NO
#else
@ -283,14 +237,6 @@
* constants YES an NO are not defined in spMatrix.h to avoid conflicts
* with user code, so use 0 for NO and 1 for YES.
*/
#endif /* spINSIDE_SPARSE */
#define spCOMPLEX 1
#define spSEPARATED_COMPLEX_VECTORS 1
#define spCOMPATIBILITY 0
#ifdef spINSIDE_SPARSE

264
src/maths/sparse/spdefs.h
View File

@ -1,3 +1,5 @@
#ifndef _SPDEFS_H
#define _SPDEFS_H
/*
* DATA STRUCTURE AND MACRO DEFINITIONS for Sparse.
*
@ -36,96 +38,29 @@
* IMPORTS
*/
#include <ngspice.h>
#include <stdio.h>
#include "ngspice.h"
#undef ABORT
#undef MALLOC
#undef FREE
#undef REALLOC
/*
* If running lint, change some of the compiler options to get a more
* complete inspection.
*/
#ifdef lint
#undef REAL
#undef spCOMPLEX
#undef EXPANDABLE
#undef TRANSLATE
#undef INITIALIZE
#undef 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
#undef spCOMPATIBILITY
#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 spCOMPATIBILITY YES
#define LINT YES
#else /* not lint */
#define LINT NO
#endif /* not lint */
/*
* MACRO DEFINITIONS
*
* Macros are distinguished by using solely capital letters in their
* identifiers. This contrasts with C defined identifiers which are strictly
* lower case, and program variable and procedure names which use both upper
* and lower case.
*/
* identifiers. This contrasts with C defined identifiers which are
* strictly lower case, and program variable and procedure names
* which use both upper and lower case. */
/* Begin macros. */
/* Boolean data type */
#define BOOLEAN int
#define NO 0
#define YES 1
#define NOT !
#define AND &&
#define OR ||
/* NULL pointer */
#ifndef NULL
#define NULL 0
#endif
#define SPARSE_ID 0x772773 /* Arbitrary (is Sparse on phone). */
#define IS_SPARSE(matrix) ((matrix) != NULL && \
@ -151,11 +86,7 @@
#define SWAP(type, a, b) {type swapx; swapx = a; a = b; b = swapx;}
/* Macro function that returns the approx absolute value of a complex number. */
#if spCOMPLEX
#define ELEMENT_MAG(ptr) (ABS((ptr)->Real) + ABS((ptr)->Imag))
#else
#define ELEMENT_MAG(ptr) ((ptr)->Real < 0.0 ? -(ptr)->Real : (ptr)->Real)
#endif
/* Complex assignment statements. */
#define CMPLX_ASSIGN(to,from) \
@ -317,114 +248,6 @@
(from_a).Imag * (from_b).Real; \
}
/*
* Macro functions that provide complex division.
*/
/* Complex division: to = num / den */
#define CMPLX_DIV(to,num,den) \
{ RealNumber r_, s_; \
if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \
((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \
{ r_ = (den).Imag / (den).Real; \
s_ = (den).Real + r_*(den).Imag; \
(to).Real = ((num).Real + r_*(num).Imag)/s_; \
(to).Imag = ((num).Imag - r_*(num).Real)/s_; \
} \
else \
{ r_ = (den).Real / (den).Imag; \
s_ = (den).Imag + r_*(den).Real; \
(to).Real = (r_*(num).Real + (num).Imag)/s_; \
(to).Imag = (r_*(num).Imag - (num).Real)/s_; \
} \
}
/* Complex division and assignment: num /= den */
#define CMPLX_DIV_ASSIGN(num,den) \
{ RealNumber r_, s_, t_; \
if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \
((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \
{ r_ = (den).Imag / (den).Real; \
s_ = (den).Real + r_*(den).Imag; \
t_ = ((num).Real + r_*(num).Imag)/s_; \
(num).Imag = ((num).Imag - r_*(num).Real)/s_; \
(num).Real = t_; \
} \
else \
{ r_ = (den).Real / (den).Imag; \
s_ = (den).Imag + r_*(den).Real; \
t_ = (r_*(num).Real + (num).Imag)/s_; \
(num).Imag = (r_*(num).Imag - (num).Real)/s_; \
(num).Real = t_; \
} \
}
/* Complex reciprocation: to = 1.0 / den */
#define CMPLX_RECIPROCAL(to,den) \
{ RealNumber r_; \
if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \
((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \
{ r_ = (den).Imag / (den).Real; \
(to).Imag = -r_*((to).Real = 1.0/((den).Real + r_*(den).Imag)); \
} \
else \
{ r_ = (den).Real / (den).Imag; \
(to).Real = -r_*((to).Imag = -1.0/((den).Imag + r_*(den).Real));\
} \
}
/*
* ASSERT and ABORT
*
* Macro used to assert that if the code is working correctly, then
* a condition must be true. If not, then execution is terminated
* and an error message is issued stating that there is an internal
* error and giving the file and line number. These assertions are
* not evaluated unless the DEBUG flag is true.
*/
#if DEBUG
#define ASSERT(condition) if (NOT(condition)) ABORT()
#else
#define ASSERT(condition)
#endif
#if DEBUG
#define ABORT() \
{ (void)fflush(stdout); \
(void)fprintf(stderr, "sparse: panic in file `%s' at line %d.\n", \
__FILE__, __LINE__); \
(void)fflush(stderr); \
abort(); \
}
#else
#define ABORT()
#endif
/*
* IMAGINARY VECTORS
*
* The imaginary vectors iRHS and iSolution are only needed when the
* options spCOMPLEX and spSEPARATED_COMPLEX_VECTORS are set. The following
* macro makes it easy to include or exclude these vectors as needed.
*/
#if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
#define IMAG_VECTORS , iRHS, iSolution
#define IMAG_RHS , iRHS
#else
#define IMAG_VECTORS
#define IMAG_RHS
#endif
#define ALLOC(type,number) ((type *)tmalloc((unsigned)(sizeof(type)*(number))))
#define REALLOC(ptr,type,number) \
@ -491,10 +314,9 @@
/* Begin `MatrixElement'. */
struct MatrixElement
{ RealNumber Real;
#if spCOMPLEX
{
RealNumber Real;
RealNumber Imag;
#endif
int Row;
int Col;
struct MatrixElement *NextInRow;
@ -532,7 +354,8 @@ typedef ElementPtr *ArrayOfElementPtrs;
/* Begin `AllocationRecord'. */
struct AllocationRecord
{ char *AllocatedPtr;
{
char *AllocatedPtr;
struct AllocationRecord *NextRecord;
};
@ -567,7 +390,8 @@ typedef struct AllocationRecord *AllocationListPtr;
/* Begin `FillinListNodeStruct'. */
struct FillinListNodeStruct
{ ElementPtr pFillinList;
{
ElementPtr pFillinList;
int NumberOfFillinsInList;
struct FillinListNodeStruct *Next;
};
@ -575,7 +399,8 @@ struct FillinListNodeStruct
/* Similar to above, but keeps track of the original Elements */
/* Begin `ElementListNodeStruct'. */
struct ElementListNodeStruct
{ ElementPtr pElementList;
{
ElementPtr pElementList;
int NumberOfElementsInList;
struct ElementListNodeStruct *Next;
};
@ -613,7 +438,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 (BOOLEAN)
* Complex (int)
* The flag which indicates whether the matrix is complex (true) or
* real.
* CurrentSize (int)
@ -622,12 +447,12 @@ struct ElementListNodeStruct
* rows and columns that have elements in them.
* Diag (ArrayOfElementPtrs)
* Array of pointers that points to the diagonal elements.
* DoCmplxDirect (BOOLEAN *)
* DoCmplxDirect (int *)
* 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 (BOOLEAN *)
* DoRealDirect (int *)
* 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
@ -644,7 +469,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 (BOOLEAN)
* Factored (int)
* Indicates if matrix has been factored. This flag is set true in
* spFactor() and spOrderAndFactor() and set false in spCreate()
* and spClear().
@ -665,7 +490,7 @@ 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 (BOOLEAN)
* InternalVectorsAllocated (int)
* A flag that indicates whether the Markowitz vectors and the
* Intermediate vector have been created.
* These vectors are created in spcCreateInternalVectors().
@ -689,19 +514,19 @@ 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 (BOOLEAN)
* NeedsOrdering (int)
* 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 (BOOLEAN)
* NumberOfInterchangesIsOdd (int)
* 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 (BOOLEAN)
* Partitioned (int)
* 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().
@ -711,7 +536,7 @@ struct ElementListNodeStruct
* Row pivot was chosen from.
* PivotSelectionMethod (char)
* Character that indicates which pivot search method was successful.
* PreviousMatrixWasComplex (BOOLEAN)
* PreviousMatrixWasComplex (int)
* 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
@ -720,11 +545,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 (BOOLEAN)
* Reordered (int)
* 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 (BOOLEAN)
* RowsLinked (int)
* 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().
@ -781,43 +606,44 @@ struct ElementListNodeStruct
/* Begin `MatrixFrame'. */
struct MatrixFrame
{ RealNumber AbsThreshold;
{
RealNumber AbsThreshold;
int AllocatedSize;
int AllocatedExtSize;
BOOLEAN Complex;
int Complex;
int CurrentSize;
ArrayOfElementPtrs Diag;
BOOLEAN *DoCmplxDirect;
BOOLEAN *DoRealDirect;
int *DoCmplxDirect;
int *DoRealDirect;
int Elements;
int Error;
int ExtSize;
int *ExtToIntColMap;
int *ExtToIntRowMap;
BOOLEAN Factored;
int Factored;
int Fillins;
ArrayOfElementPtrs FirstInCol;
ArrayOfElementPtrs FirstInRow;
unsigned long ID;
RealVector Intermediate;
BOOLEAN InternalVectorsAllocated;
int InternalVectorsAllocated;
int *IntToExtColMap;
int *IntToExtRowMap;
int *MarkowitzRow;
int *MarkowitzCol;
long *MarkowitzProd;
int MaxRowCountInLowerTri;
BOOLEAN NeedsOrdering;
BOOLEAN NumberOfInterchangesIsOdd;
int NeedsOrdering;
int NumberOfInterchangesIsOdd;
int Originals;
BOOLEAN Partitioned;
int Partitioned;
int PivotsOriginalCol;
int PivotsOriginalRow;
char PivotSelectionMethod;
BOOLEAN PreviousMatrixWasComplex;
int PreviousMatrixWasComplex;
RealNumber RelThreshold;
BOOLEAN Reordered;
BOOLEAN RowsLinked;
int Reordered;
int RowsLinked;
int SingularCol;
int SingularRow;
int Singletons;
@ -844,7 +670,6 @@ typedef struct MatrixFrame *MatrixPtr;
* Function declarations
*/
#ifdef __STDC__
extern ElementPtr spcGetElement( MatrixPtr );
extern ElementPtr spcGetFillin( MatrixPtr );
extern ElementPtr spcFindElementInCol( MatrixPtr, ElementPtr*, int, int, int );
@ -853,13 +678,6 @@ extern void spcCreateInternalVectors( MatrixPtr );
extern void spcLinkRows( MatrixPtr );
extern void spcColExchange( MatrixPtr, int, int );
extern void spcRowExchange( MatrixPtr, int, int );
#else /* __STDC__ */
extern ElementPtr spcGetElement();
extern ElementPtr spcGetFillin();
extern ElementPtr spcFindElementInCol();
extern ElementPtr spcCreateElement();
extern void spcCreateInternalVectors();
extern void spcLinkRows();
extern void spcColExchange();
extern void spcRowExchange();
#endif /* __STDC__ */
#endif

1227
src/maths/sparse/spfactor.c
View File

File diff suppressed because it is too large Load Diff

414
src/maths/sparse/spoutput.c
View File

@ -45,6 +45,7 @@
* spDefs.h
* Matrix type and macro definitions for the sparse matrix routines.
*/
#include <assert.h>
#define spINSIDE_SPARSE
#include "spconfig.h"
@ -130,30 +131,30 @@ int Printer_Width = PRINTER_WIDTH;
*/
void
spPrint( eMatrix, PrintReordered, Data, Header )
char *eMatrix;
int PrintReordered, Data, Header;
spPrint(void *eMatrix, int PrintReordered, int Data, int Header)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register int J = 0;
int I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns, ElementCount = 0;
double Magnitude, SmallestDiag, SmallestElement;
double LargestElement = 0.0, LargestDiag = 0.0;
ElementPtr pElement, *pImagElements;
int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
int J = 0;
int I, Row, Col, Size, Top;
int StartCol = 1, StopCol, Columns, ElementCount = 0;
double Magnitude, SmallestDiag, SmallestElement;
double LargestElement = 0.0, LargestDiag = 0.0;
ElementPtr pElement, *pImagElements;
int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
/* Begin `spPrint'. */
ASSERT( IS_SPARSE( Matrix ) );
/* Begin `spPrint'. */
assert( IS_SPARSE( Matrix ) );
Size = Matrix->Size;
CALLOC(pImagElements, ElementPtr, Printer_Width / 10 + 1);
if ( pImagElements == NULL)
{ Matrix->Error = spNO_MEMORY;
{
Matrix->Error = spNO_MEMORY;
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
@ -161,37 +162,43 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
#endif
CALLOC( PrintOrdToIntRowMap, int, Top + 1 );
if ( PrintOrdToIntRowMap == NULL)
{ Matrix->Error = spNO_MEMORY;
{
Matrix->Error = spNO_MEMORY;
FREE(pImagElements);
return;
}
CALLOC( PrintOrdToIntColMap, int, Top + 1 );
if (PrintOrdToIntColMap == NULL)
{ Matrix->Error = spNO_MEMORY;
{
Matrix->Error = spNO_MEMORY;
FREE(pImagElements);
FREE(PrintOrdToIntRowMap);
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 AND PrintReordered )
if ( Matrix->Reordered && PrintReordered )
printf("Matrix has been reordered.\n");
putchar('\n');
@ -204,29 +211,30 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
SmallestDiag = SmallestElement;
}
/* Determine how many columns to use. */
/* 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;
{
/* Calculatestrchr 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];
@ -235,64 +243,70 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
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 AND NOT Data)
{
if (PrintReordered && !Data)
printf("%4d", I);
else
printf("%4d", Matrix->IntToExtRowMap[ Row ]);
if (NOT Data) putchar(' ');
if (!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 AND pElement->Row != Row)
while(pElement != NULL && pElement->Row != Row)
pElement = pElement->NextInCol;
if (Data)
pImagElements[J - StartCol] = pElement;
if (pElement != NULL)
/* Case where element exists */
{ if (Data)
{
/* 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) AND (Magnitude != 0.0))
if ((Magnitude < SmallestElement) && (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('.');
@ -300,55 +314,61 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
}
putchar('\n');
#if spCOMPLEX
if (Matrix->Complex AND Data)
{ printf(" ");
if (Matrix->Complex && Data)
{
printf(" ");
for (J = StartCol; J <= StopCol; J++)
{ if (pImagElements[J - StartCol] != NULL)
{ printf(" %8.2gj",
{
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. */
/* Calculatestrchr 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. */
/* Calculate and print sparsity and number of fill-ins created. */
printf("\nDensity = %2.2f%%.\n", ((double)(ElementCount * 100)) /
((double)(Size * Size)));
((double)(Size * Size)));
printf("Number of originals = %1d.\n", Matrix->Originals);
if (NOT Matrix->NeedsOrdering)
if (!Matrix->NeedsOrdering)
printf("Number of fill-ins = %1d.\n", Matrix->Fillins);
}
putchar('\n');
@ -387,14 +407,14 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
* Name of file into which matrix is to be written.
* Label <input> (char *)
* String that is transferred to file and is used as a label.
* Reordered <input> (BOOLEAN)
* Reordered <input> (int)
* Specifies whether matrix should be output in reordered form,
* or in original order.
* Data <input> (BOOLEAN)
* Data <input> (int)
* Indicates that the element values should be output along with
* the indices for each element. This parameter must be TRUE if
* matrix is to be read by the sparse test program.
* Header <input> (BOOLEAN)
* Header <input> (int)
* Indicates that header is desired. This parameter must be TRUE if
* matrix is to be read by the sparse test program.
*
@ -412,113 +432,124 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
*/
int
spFileMatrix( eMatrix, File, Label, Reordered, Data, Header )
char *eMatrix, *Label, *File;
int Reordered, Data, Header;
spFileMatrix(void *eMatrix, char *File, char *Label, int Reordered,
int Data, int Header)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register int I, Size;
register ElementPtr pElement;
int Row, Col, Err;
FILE *pMatrixFile, *fopen();
MatrixPtr Matrix = (MatrixPtr)eMatrix;
int I, Size;
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 AND Data)
{ Err = fprintf
( pMatrixFile,
"Warning : The following matrix is factored in to LU form.\n"
);
if (Err < 0) return 0;
{
if (Matrix->Factored && 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;
}
/* Output matrix. */
if (NOT Data)
{ for (I = 1; I <= Size; I++)
{ pElement = Matrix->FirstInCol[I];
/* Output matrix. */
if (!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 spCOMPLEX
if (Data AND Matrix->Complex)
{ for (I = 1; I <= Size; I++)
{ pElement = Matrix->FirstInCol[I];
if (Data && 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, (double)pElement->Real, (double)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 REAL
if (Data AND NOT Matrix->Complex)
{ for (I = 1; I <= Size; I++)
{ pElement = Matrix->FirstInCol[I];
if (Data && !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, (double)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;
}
@ -547,11 +578,9 @@ FILE *pMatrixFile, *fopen();
* File <input> (char *)
* Name of file into which matrix is to be written.
* RHS <input> (RealNumber [])
* Right-hand side vector. This is only the real portion if
* spSEPARATED_COMPLEX_VECTORS is TRUE.
* Right-hand side vector, real portion
* iRHS <input> (RealNumber [])
* Right-hand side vector, imaginary portion. Not necessary if matrix
* is real or if spSEPARATED_COMPLEX_VECTORS is set FALSE.
* Right-hand side vector, imaginary portion.
*
* >>> Local variables:
* pMatrixFile (FILE *)
@ -559,86 +588,47 @@ FILE *pMatrixFile, *fopen();
* Size (int)
* The size of the matrix.
*
* >>> Obscure Macros
* IMAG_RHS
* Replaces itself with `, iRHS' if the options spCOMPLEX and
* spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
* without a trace.
*/
int
spFileVector( eMatrix, File, RHS IMAG_RHS )
char *eMatrix, *File;
RealVector RHS IMAG_RHS;
spFileVector(void *eMatrix, char *File, RealVector RHS, RealVector iRHS)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register int I, Size, Err;
FILE *pMatrixFile;
FILE *fopen();
MatrixPtr Matrix = (MatrixPtr)eMatrix;
int I, Size, Err;
FILE *pMatrixFile;
FILE *fopen();
/* Begin `spFileVector'. */
ASSERT( IS_SPARSE( Matrix ) AND RHS != NULL)
/* Begin `spFileVector'. */
assert( IS_SPARSE( Matrix ) && RHS != NULL);
/* Open File in append mode. */
if ((pMatrixFile = fopen(File,"a")) == NULL)
return 0;
/* Open File in append mode. */
pMatrixFile = fopen(File,"a");
if (pMatrixFile == NULL)
return 0;
/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spCOMPLEX
if (Matrix->Complex)
{
#if spSEPARATED_COMPLEX_VECTORS
ASSERT(iRHS != NULL)
--RHS;
--iRHS;
#else
RHS -= 2;
#endif
}
else
#endif /* spCOMPLEX */
--RHS;
#endif /* NOT ARRAY_OFFSET */
/* Output vector. */
/* Output vector. */
Size = Matrix->Size;
#if spCOMPLEX
if (Matrix->Complex)
{
#if spSEPARATED_COMPLEX_VECTORS
for (I = 1; I <= Size; I++)
{ Err = fprintf
( pMatrixFile, "%-.15g\t%-.15g\n",
(double)RHS[I], (double)iRHS[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)
{
for (I = 1; I <= Size; I++)
{
if (fprintf(pMatrixFile, "%-.15g\n", (double)RHS[I]) < 0)
return 0;
}
}
#endif /* REAL */
/* Close file. */
/* Close file. */
if (fclose(pMatrixFile) < 0) return 0;
return 1;
}
@ -688,27 +678,25 @@ FILE *fopen();
*/
int
spFileStats( eMatrix, File, Label )
char *eMatrix, *File, *Label;
spFileStats(void *eMatrix, char *File, char *Label)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register int Size, I;
register ElementPtr pElement;
int NumberOfElements;
RealNumber Data, LargestElement, SmallestElement;
FILE *pStatsFile, *fopen();
MatrixPtr Matrix = (MatrixPtr)eMatrix;
int Size, I;
ElementPtr pElement;
int NumberOfElements;
RealNumber Data, LargestElement, SmallestElement;
FILE *pStatsFile, *fopen();
/* 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 (NOT Matrix->Factored)
if (!Matrix->Factored)
fprintf(pStatsFile, "Matrix has not been factored.\n");
fprintf(pStatsFile, "||| Starting new matrix |||\n");
fprintf(pStatsFile, "%s\n", Label);
@ -718,19 +706,21 @@ FILE *pStatsFile, *fopen();
fprintf(pStatsFile, "Matrix is real.\n");
fprintf(pStatsFile," Size = %d\n",Size);
/* 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 AND Data != 0.0)
if (Data < SmallestElement && Data != 0.0)
SmallestElement = Data;
pElement = pElement->NextInCol;
}
@ -738,7 +728,7 @@ FILE *pStatsFile, *fopen();
SmallestElement = MIN( SmallestElement, LargestElement );
/* Output remaining statistics. */
/* Output remaining statistics. */
fprintf(pStatsFile, " Initial number of elements = %d\n",
NumberOfElements - Matrix->Fillins);
fprintf(pStatsFile,
@ -758,7 +748,7 @@ FILE *pStatsFile, *fopen();
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;
}

237
src/maths/sparse/spsmp.c
View File

@ -47,12 +47,10 @@
* 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
@ -66,11 +64,7 @@
* STABILITY NO
* CONDITION NO
* PSEUDOCONDITION NO
* FORTRAN NO
* DEBUG YES
* spCOMPLEX 1
* spSEPARATED_COMPLEX_VECTORS 1
* spCOMPATIBILITY 0
*
* spREAL double
*/
@ -100,63 +94,56 @@
* Spice3's matrix macro definitions.
*/
#include "ngspice.h"
#include <config.h>
#include <ngspice.h>
#include <assert.h>
#include <stdio.h>
#include "spmatrix.h"
#include "smpdefs.h"
#include <spmatrix.h>
#include <smpdefs.h>
#include "spdefs.h"
/* #define NO 0 */
/* #define YES 1 */
#ifdef __STDC__
static void LoadGmin( char * /*eMatrix*/, double /*Gmin*/ );
#else
static void LoadGmin( );
#endif
static void LoadGmin(SMPmatrix *eMatrix, double Gmin);
/*
* SMPaddElt()
*/
int
SMPaddElt( Matrix, Row, Col, Value )
SMPmatrix *Matrix;
int Row, Col;
double Value;
SMPaddElt(SMPmatrix *Matrix, int Row, int Col, double Value)
{
*spGetElement( (char *)Matrix, Row, Col ) = Value;
return spError( (char *)Matrix );
*spGetElement( (void *)Matrix, Row, Col ) = Value;
return spError( (void *)Matrix );
}
/*
* SMPmakeElt()
*/
double *
SMPmakeElt( Matrix, Row, Col )
SMPmatrix *Matrix;
int Row, Col;
SMPmakeElt(SMPmatrix *Matrix, int Row, int Col)
{
return spGetElement( (char *)Matrix, Row, Col );
return spGetElement( (void *)Matrix, Row, Col );
}
/*
* SMPcClear()
*/
void
SMPcClear( Matrix )
SMPmatrix *Matrix;
SMPcClear(SMPmatrix *Matrix)
{
spClear( (char *)Matrix );
spClear( (void *)Matrix );
}
/*
* SMPclear()
*/
void
SMPclear( Matrix )
SMPmatrix *Matrix;
SMPclear(SMPmatrix *Matrix)
{
spClear( (char *)Matrix );
spClear( (void *)Matrix );
}
/*
@ -164,12 +151,10 @@ SMPmatrix *Matrix;
*/
/*ARGSUSED*/
int
SMPcLUfac( Matrix, PivTol )
SMPmatrix *Matrix;
double PivTol;
SMPcLUfac(SMPmatrix *Matrix, double PivTol)
{
spSetComplex( (char *)Matrix );
return spFactor( (char *)Matrix );
spSetComplex( (void *)Matrix );
return spFactor( (void *)Matrix );
}
/*
@ -177,27 +162,23 @@ double PivTol;
*/
/*ARGSUSED*/
int
SMPluFac( Matrix, PivTol, Gmin )
SMPmatrix *Matrix;
double PivTol, Gmin;
SMPluFac(SMPmatrix *Matrix, double PivTol, double Gmin)
{
spSetReal( (char *)Matrix );
LoadGmin( (char *)Matrix, Gmin );
return spFactor( (char *)Matrix );
spSetReal( (void *)Matrix );
LoadGmin( (void *)Matrix, Gmin );
return spFactor( (void *)Matrix );
}
/*
* SMPcReorder()
*/
int
SMPcReorder( Matrix, PivTol, PivRel, NumSwaps )
SMPmatrix *Matrix;
double PivTol, PivRel;
int *NumSwaps;
SMPcReorder(SMPmatrix *Matrix, double PivTol, double PivRel,
int *NumSwaps)
{
*NumSwaps = 1;
spSetComplex( (char *)Matrix );
return spOrderAndFactor( (char *)Matrix, (spREAL*)NULL,
spSetComplex( (void *)Matrix );
return spOrderAndFactor( (void *)Matrix, (spREAL*)NULL,
(spREAL)PivRel, (spREAL)PivTol, YES );
}
@ -205,13 +186,11 @@ int *NumSwaps;
* SMPreorder()
*/
int
SMPreorder( Matrix, PivTol, PivRel, Gmin )
SMPmatrix *Matrix;
double PivTol, PivRel, Gmin;
SMPreorder(SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin)
{
spSetReal( (char *)Matrix );
LoadGmin( (char *)Matrix, Gmin );
return spOrderAndFactor( (char *)Matrix, (spREAL*)NULL,
spSetReal( (void *)Matrix );
LoadGmin( (void *)Matrix, Gmin );
return spOrderAndFactor( (void *)Matrix, (spREAL*)NULL,
(spREAL)PivRel, (spREAL)PivTol, YES );
}
@ -219,53 +198,47 @@ double PivTol, PivRel, Gmin;
* SMPcaSolve()
*/
void
SMPcaSolve( Matrix, RHS, iRHS, Spare, iSpare)
SMPmatrix *Matrix;
double RHS[], iRHS[], Spare[], iSpare[];
SMPcaSolve(SMPmatrix *Matrix, double RHS[], double iRHS[],
double Spare[], double iSpare[])
{
spSolveTransposed( (char *)Matrix, RHS, RHS, iRHS, iRHS );
spSolveTransposed( (void *)Matrix, RHS, RHS, iRHS, iRHS );
}
/*
* SMPcSolve()
*/
void
SMPcSolve( Matrix, RHS, iRHS, Spare, iSpare)
SMPmatrix *Matrix;
double RHS[], iRHS[], Spare[], iSpare[];
SMPcSolve(SMPmatrix *Matrix, double RHS[], double iRHS[],
double Spare[], double iSpare[])
{
spSolve( (char *)Matrix, RHS, RHS, iRHS, iRHS );
spSolve( (void *)Matrix, RHS, RHS, iRHS, iRHS );
}
/*
* SMPsolve()
*/
void
SMPsolve( Matrix, RHS, Spare )
SMPmatrix *Matrix;
double RHS[], Spare[];
SMPsolve(SMPmatrix *Matrix, double RHS[], double Spare[])
{
spSolve( (char *)Matrix, RHS, RHS, (spREAL*)NULL, (spREAL*)NULL );
spSolve( (void *)Matrix, RHS, RHS, (spREAL*)NULL, (spREAL*)NULL );
}
/*
* SMPmatSize()
*/
int
SMPmatSize( Matrix )
SMPmatrix *Matrix;
SMPmatSize(SMPmatrix *Matrix)
{
return spGetSize( (char *)Matrix, 1 );
return spGetSize( (void *)Matrix, 1 );
}
/*
* SMPnewMatrix()
*/
int
SMPnewMatrix( pMatrix )
SMPmatrix **pMatrix;
SMPnewMatrix(SMPmatrix **pMatrix)
{
int Error;
int Error;
*pMatrix = (SMPmatrix *)spCreate( 0, 1, &Error );
return Error;
}
@ -274,21 +247,19 @@ int Error;
* SMPdestroy()
*/
void
SMPdestroy( Matrix )
SMPmatrix *Matrix;
SMPdestroy(SMPmatrix *Matrix)
{
spDestroy( (char *)Matrix );
spDestroy( (void *)Matrix );
}
/*
* SMPpreOrder()
*/
int
SMPpreOrder( Matrix )
SMPmatrix *Matrix;
SMPpreOrder(SMPmatrix *Matrix)
{
spMNA_Preorder( (char *)Matrix );
return spError( (char *)Matrix );
spMNA_Preorder( (void *)Matrix );
return spError( (void *)Matrix );
}
/*
@ -296,22 +267,18 @@ SMPmatrix *Matrix;
*/
/*ARGSUSED*/
void
SMPprint( Matrix, File )
SMPmatrix *Matrix;
FILE *File;
SMPprint(SMPmatrix *Matrix, FILE *File)
{
spPrint( (char *)Matrix, 0, 1, 1 );
spPrint( (void *)Matrix, 0, 1, 1 );
}
/*
* SMPgetError()
*/
void
SMPgetError( Matrix, Col, Row)
SMPmatrix *Matrix;
int *Row, *Col;
SMPgetError(SMPmatrix *Matrix, int *Col, int *Row)
{
spWhereSingular( (char *)Matrix, Row, Col );
spWhereSingular( (void *)Matrix, Row, Col );
}
/*
@ -319,29 +286,23 @@ int *Row, *Col;
* note: obsolete for Spice3d2 and later
*/
int
SMPcProdDiag( Matrix, pMantissa, pExponent)
SMPmatrix *Matrix;
SPcomplex *pMantissa;
int *pExponent;
SMPcProdDiag(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
{
spDeterminant( (char *)Matrix, pExponent, &(pMantissa->real),
spDeterminant( (void *)Matrix, pExponent, &(pMantissa->real),
&(pMantissa->imag) );
return spError( (char *)Matrix );
return spError( (void *)Matrix );
}
/*
* SMPcDProd()
*/
int
SMPcDProd( Matrix, pMantissa, pExponent)
SMPmatrix *Matrix;
SPcomplex *pMantissa;
int *pExponent;
SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
{
double re, im, x, y, z;
int p;
spDeterminant( (char *)Matrix, &p, &re, &im);
spDeterminant( (void *)Matrix, &p, &re, &im);
#ifndef M_LN2
#define M_LN2 0.69314718055994530942
@ -411,7 +372,7 @@ int *pExponent;
printf("Determinant 10->2: (%20g,%20g)^%d\n", pMantissa->real,
pMantissa->imag, *pExponent);
#endif
return spError( (char *)Matrix );
return spError( (void *)Matrix );
}
@ -433,17 +394,15 @@ int *pExponent;
*/
static void
LoadGmin( eMatrix, Gmin )
char *eMatrix;
register double Gmin;
LoadGmin(SMPmatrix *eMatrix, double Gmin)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register int I;
register ArrayOfElementPtrs Diag;
register ElementPtr diag;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
int I;
ArrayOfElementPtrs Diag;
ElementPtr diag;
/* Begin `LoadGmin'. */
ASSERT( IS_SPARSE( Matrix ) );
/* Begin `LoadGmin'. */
assert( IS_SPARSE( Matrix ) );
if (Gmin != 0.0) {
Diag = Matrix->Diag;
@ -468,17 +427,13 @@ register ElementPtr diag;
*/
SMPelement *
SMPfindElt( eMatrix, Row, Col, CreateIfMissing )
SMPmatrix *eMatrix;
int Row, Col;
int CreateIfMissing;
SMPfindElt(SMPmatrix *eMatrix, int Row, int Col, int CreateIfMissing)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
ElementPtr Element;
MatrixPtr Matrix = (MatrixPtr)eMatrix;
ElementPtr Element;
/* Begin `SMPfindElt'. */
ASSERT( IS_SPARSE( Matrix ) );
/* Begin `SMPfindElt'. */
assert( IS_SPARSE( Matrix ) );
Row = Matrix->ExtToIntRowMap[Row];
Col = Matrix->ExtToIntColMap[Col];
Element = Matrix->FirstInCol[Col];
@ -492,10 +447,9 @@ ElementPtr Element;
* SMPcZeroCol()
*/
int
SMPcZeroCol( Matrix, Col )
MatrixPtr Matrix;
int Col;
SMPcZeroCol(SMPmatrix *eMatrix, int Col)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
ElementPtr Element;
Col = Matrix->ExtToIntColMap[Col];
@ -508,17 +462,16 @@ int Col;
Element->Imag = 0.0;
}
return spError( (char *)Matrix );
return spError( (void *)Matrix );
}
/*
* SMPcAddCol()
*/
int
SMPcAddCol( Matrix, Accum_Col, Addend_Col )
MatrixPtr Matrix;
int Accum_Col, Addend_Col;
SMPcAddCol(SMPmatrix *eMatrix, int Accum_Col, int Addend_Col)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
ElementPtr Accum, Addend, *Prev;
Accum_Col = Matrix->ExtToIntColMap[Accum_Col];
@ -541,17 +494,16 @@ int Accum_Col, Addend_Col;
Addend = Addend->NextInCol;
}
return spError( (char *)Matrix );
return spError( (void *)Matrix );
}
/*
* SMPzeroRow()
*/
int
SMPzeroRow( Matrix, Row )
MatrixPtr Matrix;
int Row;
SMPzeroRow(SMPmatrix *eMatrix, int Row)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
ElementPtr Element;
Row = Matrix->ExtToIntColMap[Row];
@ -559,8 +511,7 @@ int Row;
if (Matrix->RowsLinked == NO)
spcLinkRows(Matrix);
#if spCOMPLEX
if (Matrix->PreviousMatrixWasComplex OR Matrix->Complex) {
if (Matrix->PreviousMatrixWasComplex || Matrix->Complex) {
for (Element = Matrix->FirstInRow[Row];
Element != NULL;
Element = Element->NextInRow)
@ -568,9 +519,7 @@ int Row;
Element->Real = 0.0;
Element->Imag = 0.0;
}
} else
#endif
{
} else {
for (Element = Matrix->FirstInRow[Row];
Element != NULL;
Element = Element->NextInRow)
@ -579,7 +528,7 @@ int Row;
}
}
return spError( (char *)Matrix );
return spError( (void *)Matrix );
}
#ifdef PARALLEL_ARCH
@ -587,24 +536,20 @@ int Row;
* SMPcombine()
*/
void
SMPcombine( Matrix, RHS, Spare )
SMPmatrix *Matrix;
double RHS[], Spare[];
SMPcombine(SMPmatrix *Matrix, double RHS[], double Spare[])
{
spSetReal( (char *)Matrix );
spCombine( (char *)Matrix, RHS, Spare, (spREAL*)NULL, (spREAL*)NULL );
spSetReal( (void *)Matrix );
spCombine( (void *)Matrix, RHS, Spare, (spREAL*)NULL, (spREAL*)NULL );
}
/*
* SMPcCombine()
*/
void
SMPcCombine( Matrix, RHS, Spare, iRHS, iSpare )
SMPmatrix *Matrix;
double RHS[], Spare[];
double iRHS[], iSpare[];
SMPcCombine(SMPmatrix *Matrix, double RHS[], double Spare[],
double iRHS[], double iSpare[])
{
spSetComplex( (char *)Matrix );
spCombine( (char *)Matrix, RHS, Spare, iRHS, iSpare );
spSetComplex( (void *)Matrix );
spCombine( (void *)Matrix, RHS, Spare, iRHS, iSpare );
}
#endif /* PARALLEL_ARCH */

325
src/maths/sparse/spsolve.c
View File

@ -45,6 +45,7 @@
* spDefs.h
* Matrix type and macro definitions for the sparse matrix routines.
*/
#include <assert.h>
#define spINSIDE_SPARSE
#include "spconfig.h"
@ -58,20 +59,10 @@
* Function declarations
*/
#ifdef __STDC__
#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
#else /* __STDC__ */
static void SolveComplexMatrix();
static void SolveComplexTransposedMatrix();
#endif /* __STDC__ */
@ -102,13 +93,10 @@ static void SolveComplexTransposedMatrix();
* iRHS <input> (RealVector)
* iRHS is the imaginary portion of the input data array, the right
* hand side. This data is undisturbed and may be reused for other solves.
* This argument is only necessary if matrix is complex and if
* spSEPARATED_COMPLEX_VECTOR is set TRUE.
* iSolution <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. This argument is only necessary if matrix is complex
* and if spSEPARATED_COMPLEX_VECTOR is set TRUE.
* the same array.
*
* >>> Local variables:
* Intermediate (RealVector)
@ -131,89 +119,77 @@ static void SolveComplexTransposedMatrix();
* 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( eMatrix, RHS, Solution IMAG_VECTORS )
char *eMatrix;
RealVector RHS, Solution IMAG_VECTORS;
spSolve(void *eMatrix, RealVector RHS, RealVector Solution,
RealVector iRHS, RealVector iSolution)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register ElementPtr pElement;
register RealVector Intermediate;
register RealNumber Temp;
register int I, *pExtOrder, Size;
ElementPtr pPivot;
void SolveComplexMatrix();
MatrixPtr Matrix = (MatrixPtr)eMatrix;
ElementPtr pElement;
RealVector Intermediate;
RealNumber Temp;
int I, *pExtOrder, Size;
ElementPtr pPivot;
void SolveComplexMatrix();
/* Begin `spSolve'. */
ASSERT( IS_VALID(Matrix) AND IS_FACTORED(Matrix) );
/* Begin `spSolve'. */
assert( IS_VALID(Matrix) && IS_FACTORED(Matrix) );
#if spCOMPLEX
if (Matrix->Complex)
{ SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS );
{
SolveComplexMatrix( Matrix, RHS, Solution, iRHS, iSolution );
return;
}
#endif
#if REAL
Intermediate = Matrix->Intermediate;
Size = Matrix->Size;
/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
--RHS;
--Solution;
#endif
/* Initialize Intermediate vector. */
/* 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 */
}
@ -226,7 +202,6 @@ void SolveComplexMatrix();
#if spCOMPLEX
/*
* SOLVE COMPLEX MATRIX EQUATION
*
@ -280,69 +255,50 @@ void SolveComplexMatrix();
* Size of matrix. Made local to reduce indirection.
* Temp (ComplexNumber)
* 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.
*/
static void
SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS )
SolveComplexMatrix( Matrix, RHS, Solution , iRHS, iSolution )
MatrixPtr Matrix;
RealVector RHS, Solution IMAG_VECTORS;
RealVector RHS, Solution , iRHS, iSolution;
{
register ElementPtr pElement;
register ComplexVector Intermediate;
register int I, *pExtOrder, Size;
ElementPtr pPivot;
ComplexNumber Temp;
ElementPtr pElement;
ComplexVector Intermediate;
int I, *pExtOrder, Size;
ElementPtr pPivot;
ComplexNumber Temp;
/* Begin `SolveComplexMatrix'. */
/* Begin `SolveComplexMatrix'. */
Size = Matrix->Size;
Intermediate = (ComplexVector)Matrix->Intermediate;
/* 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. */
/* 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) OR (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) || (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;
@ -350,37 +306,32 @@ ComplexNumber Temp;
}
}
/* 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 */
@ -448,88 +399,77 @@ ComplexNumber Temp;
* 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
spSolveTransposed( eMatrix, RHS, Solution IMAG_VECTORS )
char *eMatrix;
RealVector RHS, Solution IMAG_VECTORS;
spSolveTransposed(void *eMatrix, RealVector RHS, RealVector Solution,
RealVector iRHS, RealVector iSolution)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register ElementPtr pElement;
register RealVector Intermediate;
register int I, *pExtOrder, Size;
ElementPtr pPivot;
RealNumber Temp;
void SolveComplexTransposedMatrix();
MatrixPtr Matrix = (MatrixPtr)eMatrix;
ElementPtr pElement;
RealVector Intermediate;
int I, *pExtOrder, Size;
ElementPtr pPivot;
RealNumber Temp;
void SolveComplexTransposedMatrix();
/* Begin `spSolveTransposed'. */
ASSERT( IS_VALID(Matrix) AND IS_FACTORED(Matrix) );
/* Begin `spSolveTransposed'. */
assert( IS_VALID(Matrix) && IS_FACTORED(Matrix) );
#if spCOMPLEX
if (Matrix->Complex)
{ SolveComplexTransposedMatrix( Matrix, RHS, Solution IMAG_VECTORS );
{
SolveComplexTransposedMatrix( Matrix, RHS, Solution , iRHS, iSolution );
return;
}
#endif
#if REAL
Size = Matrix->Size;
Intermediate = Matrix->Intermediate;
/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
--RHS;
--Solution;
#endif
/* Initialize Intermediate vector. */
/* 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 */
@ -542,7 +482,7 @@ void SolveComplexTransposedMatrix();
#if TRANSPOSE AND spCOMPLEX
#if TRANSPOSE
/*
* SOLVE COMPLEX TRANSPOSED MATRIX EQUATION
*
@ -560,23 +500,18 @@ void SolveComplexTransposedMatrix();
* RHS <input> (RealVector)
* 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 and
* spSEPARATED_COMPLEX_VECTORS is set TRUE.
* 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 and
* spSEPARATED_COMPLEX_VECTORS is set TRUE.
* 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
* hand side. This data is undisturbed and may be reused for other solves.
* If either spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set FALSE, 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 spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set
* FALSE, there is no need to supply this array.
* the same array.
*
* >>> Local variables:
* Intermediate (ComplexVector)
@ -599,65 +534,46 @@ void SolveComplexTransposedMatrix();
* Size of matrix. Made local to reduce indirection.
* Temp (ComplexNumber)
* 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.
*/
static void
SolveComplexTransposedMatrix(Matrix, RHS, Solution IMAG_VECTORS )
SolveComplexTransposedMatrix(Matrix, RHS, Solution , iRHS, iSolution )
MatrixPtr Matrix;
RealVector RHS, Solution IMAG_VECTORS;
RealVector RHS, Solution , iRHS, iSolution;
{
register ElementPtr pElement;
register ComplexVector Intermediate;
register int I, *pExtOrder, Size;
ElementPtr pPivot;
ComplexNumber Temp;
ElementPtr pElement;
ComplexVector Intermediate;
int I, *pExtOrder, Size;
ElementPtr pPivot;
ComplexNumber Temp;
/* Begin `SolveComplexTransposedMatrix'. */
/* Begin `SolveComplexTransposedMatrix'. */
Size = Matrix->Size;
Intermediate = (ComplexVector)Matrix->Intermediate;
/* 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. */
/* 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) OR (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) || (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;
@ -665,37 +581,34 @@ ComplexNumber Temp;
}
}
/* 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 AND spCOMPLEX */
#endif /* TRANSPOSE */

1229
src/maths/sparse/sputils.c
View File

File diff suppressed because it is too large Load Diff