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

llvm format

This commit is contained in:
dwarning 2023-06-18 10:07:09 +02:00
parent 3638656ce6
commit 990d3a92d3
12 changed files with 3423 additions and 4508 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
src/include/ngspice/smpdefs.h
View File

@ -1,8 +1,8 @@
#ifndef ngspice_SMPDEFS_H
#define ngspice_SMPDEFS_H
typedef struct MatrixFrame SMPmatrix;
typedef struct MatrixElement SMPelement;
typedef struct MatrixFrame SMPmatrix;
typedef struct MatrixElement SMPelement;
/**********
Copyright 1990 Regents of the University of California. All rights reserved.
@ -10,32 +10,32 @@ Author: 1985 Thomas L. Quarles
Modified: 2000 AlansFixes
**********/
#include <stdio.h>
#include <math.h>
#include "ngspice/complex.h"
#include <math.h>
#include <stdio.h>
int SMPaddElt( SMPmatrix *, int , int , 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 );
void SMPdestroy( SMPmatrix *);
int SMPpreOrder( SMPmatrix *);
void SMPprint( SMPmatrix * , char *);
void SMPprintRHS( SMPmatrix * , char *, double*, double*);
void SMPgetError( SMPmatrix *, int *, int *);
int SMPcProdDiag( SMPmatrix *, SPcomplex *, int *);
int SMPaddElt(SMPmatrix *, int, int, 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);
void SMPdestroy(SMPmatrix *);
int SMPpreOrder(SMPmatrix *);
void SMPprint(SMPmatrix *, char *);
void SMPprintRHS(SMPmatrix *, char *, double *, double *);
void SMPgetError(SMPmatrix *, int *, int *);
int SMPcProdDiag(SMPmatrix *, SPcomplex *, int *);
int SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent);
SMPelement * SMPfindElt( SMPmatrix *, int , int , int );
SMPelement *SMPfindElt(SMPmatrix *, int, int, int);
int SMPcZeroCol(SMPmatrix *Matrix, int Col);
int SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col);
void SMPconstMult(SMPmatrix *, double);

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

@ -17,7 +17,6 @@
* Kenneth S. Kundert <kundert@users.sourceforge.net>
*/
/*
* Revision and copyright information.
*
@ -27,10 +26,7 @@
* $Revision: 1.2 $
*/
#ifndef spOKAY
#ifndef spOKAY
/*
* IMPORTS
@ -40,13 +36,9 @@
* Macros that customize the sparse matrix routines.
*/
#include <stdio.h>
#include "../../maths/sparse/spConfig.h"
/*
* ERROR KEYWORDS
*
@ -103,38 +95,34 @@
// * warnings.
// */
#include "ngspice/sperror.h" /* Spice error definitions. */
#include "ngspice/sperror.h" /* Spice error definitions. */
/* Begin error macros. */
#define spOKAY OK
#define spSMALL_PIVOT OK
#define spZERO_DIAG E_SINGULAR
#define spSINGULAR E_SINGULAR
#define spNO_MEMORY E_NOMEM
#define spPANIC E_BADMATRIX
#define spMANGLED E_BADMATRIX
#define spFATAL E_BADMATRIX
#define spOKAY OK
#define spSMALL_PIVOT OK
#define spZERO_DIAG E_SINGULAR
#define spSINGULAR E_SINGULAR
#define spNO_MEMORY E_NOMEM
#define spPANIC E_BADMATRIX
#define spMANGLED E_BADMATRIX
#define spFATAL E_BADMATRIX
/*
* KEYWORD DEFINITIONS
*/
#define spREAL double /*!<
* Defines the precision of the arithmetic used by
* \a Sparse will use. Double precision is suggested
* as being most appropriate for circuit simulation
* and for C. However, it is possible to change spREAL
* to a float for single precision arithmetic. Note
* that in C, single precision arithmetic is often
* slower than double precision. Sparse
* internally refers to spREALs as RealNumbers.
*/
#define spREAL \
double /*!< \
* Defines the precision of the arithmetic used by \
* \a Sparse will use. Double precision is suggested \
* as being most appropriate for circuit simulation \
* and for C. However, it is possible to change spREAL \
* to a float for single precision arithmetic. Note \
* that in C, single precision arithmetic is often \
* slower than double precision. Sparse \
* internally refers to spREALs as RealNumbers. \
*/
/*
* PARTITION TYPES
@ -159,37 +147,37 @@
/* Begin partition keywords. */
#define spDEFAULT_PARTITION 0 /*!<
* Partition code for spPartition().
* Indicates that the default partitioning
* mode should be used.
* \see spPartition()
*/
#define spDIRECT_PARTITION 1 /*!<
* Partition code for spPartition().
* Indicates that all rows should be placed
* in the direct addressing partition.
* \see spPartition()
*/
#define spINDIRECT_PARTITION 2 /*!<
* Partition code for spPartition().
* Indicates that all rows should be placed
* in the indirect addressing partition.
* \see spPartition()
*/
#define spAUTO_PARTITION 3 /*!<
* Partition code for spPartition().
* Indicates that \a Sparse should chose
* the best partition for each row based
* on some simple rules. This is generally
* preferred.
* \see spPartition()
*/
#define spDEFAULT_PARTITION \
0 /*!< \
* Partition code for spPartition(). \
* Indicates that the default partitioning \
* mode should be used. \
* \see spPartition() \
*/
#define spDIRECT_PARTITION \
1 /*!< \
* Partition code for spPartition(). \
* Indicates that all rows should be placed \
* in the direct addressing partition. \
* \see spPartition() \
*/
#define spINDIRECT_PARTITION \
2 /*!< \
* Partition code for spPartition(). \
* Indicates that all rows should be placed \
* in the indirect addressing partition. \
* \see spPartition() \
*/
#define spAUTO_PARTITION \
3 /*!< \
* Partition code for spPartition(). \
* Indicates that \a Sparse should chose \
* the best partition for each row based \
* on some simple rules. This is generally \
* preferred. \
* \see spPartition() \
*/
/*
* MACRO FUNCTION DEFINITIONS
*/
@ -198,66 +186,64 @@
/*!
* Macro function that adds data to a real element in the matrix by a pointer.
*/
#define spADD_REAL_ELEMENT(element,real) *(element) += real
#define spADD_REAL_ELEMENT(element, real) *(element) += real
/*!
* Macro function that adds data to a imaginary element in the matrix by
* a pointer.
*/
#define spADD_IMAG_ELEMENT(element,imag) *(element+1) += imag
#define spADD_IMAG_ELEMENT(element, imag) *(element + 1) += imag
/*!
* Macro function that adds data to a complex element in the matrix by
* a pointer.
*/
#define spADD_COMPLEX_ELEMENT(element,real,imag) \
{ *(element) += real; \
*(element+1) += imag; \
}
#define spADD_COMPLEX_ELEMENT(element, real, imag) \
{ \
*(element) += real; \
*(element + 1) += imag; \
}
/*!
* Macro function that adds data to each of the four real matrix elements
* specified by the given template.
*/
#define spADD_REAL_QUAD(template,real) \
{ *((template).Element1) += real; \
*((template).Element2) += real; \
*((template).Element3Negated) -= real; \
*((template).Element4Negated) -= real; \
}
#define spADD_REAL_QUAD(template, real) \
{ \
*((template).Element1) += real; \
*((template).Element2) += real; \
*((template).Element3Negated) -= real; \
*((template).Element4Negated) -= real; \
}
/*!
* Macro function that adds data to each of the four imaginary matrix
* elements specified by the given template.
*/
#define spADD_IMAG_QUAD(template,imag) \
{ *((template).Element1+1) += imag; \
*((template).Element2+1) += imag; \
*((template).Element3Negated+1) -= imag; \
*((template).Element4Negated+1) -= imag; \
}
#define spADD_IMAG_QUAD(template, imag) \
{ \
*((template).Element1 + 1) += imag; \
*((template).Element2 + 1) += imag; \
*((template).Element3Negated + 1) -= imag; \
*((template).Element4Negated + 1) -= imag; \
}
/*!
* Macro function that adds data to each of the four complex matrix
* elements specified by the given template.
*/
#define spADD_COMPLEX_QUAD(template,real,imag) \
{ *((template).Element1) += real; \
*((template).Element2) += real; \
*((template).Element3Negated) -= real; \
*((template).Element4Negated) -= real; \
*((template).Element1+1) += imag; \
*((template).Element2+1) += imag; \
*((template).Element3Negated+1) -= imag; \
*((template).Element4Negated+1) -= imag; \
}
#define spADD_COMPLEX_QUAD(template, real, imag) \
{ \
*((template).Element1) += real; \
*((template).Element2) += real; \
*((template).Element3Negated) -= real; \
*((template).Element4Negated) -= real; \
*((template).Element1 + 1) += imag; \
*((template).Element2 + 1) += imag; \
*((template).Element3Negated + 1) -= imag; \
*((template).Element4Negated + 1) -= imag; \
}
/*
* TYPE DEFINITION FOR EXTERNAL MATRIX ELEMENT REFERENCES
*
@ -271,12 +257,8 @@ typedef spGenericPtr spMatrix;
typedef spREAL spElement;
/*! Declares the type of the Sparse error codes. */
//typedef int spError;
// typedef int spError;
/* TYPE DEFINITION FOR COMPONENT TEMPLATE */
/*!
* This data structure is used to hold pointers to four related elements in
@ -291,18 +273,15 @@ typedef spREAL spElement;
*/
/* Begin `spTemplate'. */
struct spTemplate
{ spElement *Element1;
spElement *Element2;
spElement *Element3Negated;
spElement *Element4Negated;
struct spTemplate {
spElement *Element1;
spElement *Element2;
spElement *Element3Negated;
spElement *Element4Negated;
};
typedef struct MatrixFrame *MatrixPtr;
/*
* FUNCTION TYPE DEFINITIONS
*
@ -311,78 +290,69 @@ typedef struct MatrixFrame *MatrixPtr;
/* Begin function declarations. */
spcEXTERN void spClear( MatrixPtr );
spcEXTERN spREAL spCondition( MatrixPtr, spREAL, int* );
spcEXTERN MatrixPtr spCreate( int, int, int* );
spcEXTERN void spDeleteRowAndCol( MatrixPtr, int, int );
spcEXTERN void spDestroy( MatrixPtr );
spcEXTERN int spElementCount( MatrixPtr );
spcEXTERN int spOriginalCount( MatrixPtr );
spcEXTERN int spError( MatrixPtr );
spcEXTERN void spClear(MatrixPtr);
spcEXTERN spREAL spCondition(MatrixPtr, spREAL, int *);
spcEXTERN MatrixPtr spCreate(int, int, int *);
spcEXTERN void spDeleteRowAndCol(MatrixPtr, int, int);
spcEXTERN void spDestroy(MatrixPtr);
spcEXTERN int spElementCount(MatrixPtr);
spcEXTERN int spOriginalCount(MatrixPtr);
spcEXTERN int spError(MatrixPtr);
#ifdef EOF
spcEXTERN void spErrorMessage( MatrixPtr, FILE*, char* );
spcEXTERN void spErrorMessage(MatrixPtr, FILE *, char *);
#else
# define spErrorMessage(a,b,c) spcFUNC_NEEDS_FILE(_spErrorMessage,stdio)
#define spErrorMessage(a, b, c) spcFUNC_NEEDS_FILE(_spErrorMessage, stdio)
#endif
spcEXTERN int spFactor( MatrixPtr );
spcEXTERN int spFileMatrix( MatrixPtr, char*, char*, int, int, int );
spcEXTERN int spFileStats( MatrixPtr, char*, char* );
spcEXTERN int spFillinCount( MatrixPtr );
spcEXTERN spElement *spFindElement( MatrixPtr, int, int );
spcEXTERN int spGetAdmittance( MatrixPtr, int, int,
struct spTemplate* );
spcEXTERN spElement *spGetElement( MatrixPtr, int, int );
spcEXTERN spGenericPtr spGetInitInfo( spElement* );
spcEXTERN int spGetOnes( MatrixPtr, int, int, int,
struct spTemplate* );
spcEXTERN int spGetQuad( MatrixPtr, int, int, int, int,
struct spTemplate* );
spcEXTERN int spGetSize( MatrixPtr, int );
spcEXTERN int spInitialize( MatrixPtr, int (*pInit)(spElement *, spGenericPtr, int, int) );
spcEXTERN void spInstallInitInfo( spElement*, spGenericPtr );
spcEXTERN spREAL spLargestElement( MatrixPtr );
spcEXTERN void spMNA_Preorder( MatrixPtr );
spcEXTERN spREAL spNorm( MatrixPtr );
spcEXTERN int spOrderAndFactor( MatrixPtr, spREAL[], spREAL,
spREAL, int );
spcEXTERN void spPartition( MatrixPtr, int );
spcEXTERN void spPrint( MatrixPtr, int, int, int );
spcEXTERN spREAL spPseudoCondition( MatrixPtr );
spcEXTERN spREAL spRoundoff( MatrixPtr, spREAL );
spcEXTERN void spScale( MatrixPtr, spREAL[], spREAL[] );
spcEXTERN void spSetComplex( MatrixPtr );
spcEXTERN void spSetReal( MatrixPtr );
spcEXTERN void spStripFills( MatrixPtr );
spcEXTERN void spWhereSingular( MatrixPtr, int*, int* );
spcEXTERN void spConstMult( MatrixPtr, double );
spcEXTERN int spFactor(MatrixPtr);
spcEXTERN int spFileMatrix(MatrixPtr, char *, char *, int, int, int);
spcEXTERN int spFileStats(MatrixPtr, char *, char *);
spcEXTERN int spFillinCount(MatrixPtr);
spcEXTERN spElement *spFindElement(MatrixPtr, int, int);
spcEXTERN int spGetAdmittance(MatrixPtr, int, int, struct spTemplate *);
spcEXTERN spElement *spGetElement(MatrixPtr, int, int);
spcEXTERN spGenericPtr spGetInitInfo(spElement *);
spcEXTERN int spGetOnes(MatrixPtr, int, int, int, struct spTemplate *);
spcEXTERN int spGetQuad(MatrixPtr, int, int, int, int, struct spTemplate *);
spcEXTERN int spGetSize(MatrixPtr, int);
spcEXTERN int spInitialize(MatrixPtr,
int (*pInit)(spElement *, spGenericPtr, int, int));
spcEXTERN void spInstallInitInfo(spElement *, spGenericPtr);
spcEXTERN spREAL spLargestElement(MatrixPtr);
spcEXTERN void spMNA_Preorder(MatrixPtr);
spcEXTERN spREAL spNorm(MatrixPtr);
spcEXTERN int spOrderAndFactor(MatrixPtr, spREAL[], spREAL, spREAL, int);
spcEXTERN void spPartition(MatrixPtr, int);
spcEXTERN void spPrint(MatrixPtr, int, int, int);
spcEXTERN spREAL spPseudoCondition(MatrixPtr);
spcEXTERN spREAL spRoundoff(MatrixPtr, spREAL);
spcEXTERN void spScale(MatrixPtr, spREAL[], spREAL[]);
spcEXTERN void spSetComplex(MatrixPtr);
spcEXTERN void spSetReal(MatrixPtr);
spcEXTERN void spStripFills(MatrixPtr);
spcEXTERN void spWhereSingular(MatrixPtr, int *, int *);
spcEXTERN void spConstMult(MatrixPtr, double);
/* Functions with argument lists that are dependent on options. */
#if spCOMPLEX
spcEXTERN void spDeterminant( MatrixPtr, int*, spREAL*, spREAL* );
#else /* NOT spCOMPLEX */
spcEXTERN void spDeterminant( MatrixPtr, int*, spREAL* );
spcEXTERN void spDeterminant(MatrixPtr, int *, spREAL *, spREAL *);
#else /* NOT spCOMPLEX */
spcEXTERN void spDeterminant(MatrixPtr, int *, spREAL *);
#endif /* NOT spCOMPLEX */
#if spCOMPLEX && spSEPARATED_COMPLEX_VECTORS
spcEXTERN int spFileVector( MatrixPtr, char* ,
spREAL[], spREAL[]);
spcEXTERN void spMultiply( MatrixPtr, spREAL[], spREAL[],
spREAL[], spREAL[] );
spcEXTERN void spMultTransposed( MatrixPtr, spREAL[], spREAL[],
spREAL[], spREAL[] );
spcEXTERN void spSolve( MatrixPtr, spREAL[], spREAL[], spREAL[],
spREAL[] );
spcEXTERN void spSolveTransposed( MatrixPtr, spREAL[], spREAL[],
spREAL[], spREAL[] );
#else /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */
spcEXTERN int spFileVector( MatrixPtr, char* , spREAL[] );
spcEXTERN void spMultiply( MatrixPtr, spREAL[], spREAL[] );
spcEXTERN void spMultTransposed( MatrixPtr,
spREAL[], spREAL[] );
spcEXTERN void spSolve( MatrixPtr, spREAL[], spREAL[] );
spcEXTERN void spSolveTransposed( MatrixPtr,
spREAL[], spREAL[] );
spcEXTERN int spFileVector(MatrixPtr, char *, spREAL[], spREAL[]);
spcEXTERN void spMultiply(MatrixPtr, spREAL[], spREAL[], spREAL[], spREAL[]);
spcEXTERN void spMultTransposed(MatrixPtr, spREAL[], spREAL[], spREAL[],
spREAL[]);
spcEXTERN void spSolve(MatrixPtr, spREAL[], spREAL[], spREAL[], spREAL[]);
spcEXTERN void spSolveTransposed(MatrixPtr, spREAL[], spREAL[], spREAL[],
spREAL[]);
#else /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */
spcEXTERN int spFileVector(MatrixPtr, char *, spREAL[]);
spcEXTERN void spMultiply(MatrixPtr, spREAL[], spREAL[]);
spcEXTERN void spMultTransposed(MatrixPtr, spREAL[], spREAL[]);
spcEXTERN void spSolve(MatrixPtr, spREAL[], spREAL[]);
spcEXTERN void spSolveTransposed(MatrixPtr, spREAL[], spREAL[]);
#endif /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */
#endif /* spOKAY */
#endif /* spOKAY */

535
src/maths/sparse/spAllocate.c
View File

@ -6,8 +6,9 @@
* UC Berkeley
*/
/*!\file
* This file contains functions for allocating and freeing matrices, configuring them, and for
* accessing global information about the matrix (size, error status, etc.).
* This file contains functions for allocating and freeing matrices,
* configuring them, and for accessing global information about the matrix
* (size, error status, etc.).
*
* Objects that begin with the \a spc prefix are considered private
* and should not be used.
@ -36,7 +37,6 @@
* ExpandTranslationArrays
*/
/*
* Revision and copyright information.
*
@ -46,12 +46,10 @@
#ifdef notdef
static char copyright[] =
"Sparse1.4: Copyright (c) 1985-2003 by Kenneth S. Kundert";
static char RCSid[] =
"@(#)$Header: /cvsroot/sparse/src/spAllocate.c,v 1.3 2003/06/29 04:19:52 kundert Exp $";
static char RCSid[] = "@(#)$Header: /cvsroot/sparse/src/spAllocate.c,v 1.3 "
"2003/06/29 04:19:52 kundert Exp $";
#endif
/*
* IMPORTS
*
@ -65,14 +63,10 @@ static char RCSid[] =
*/
#define spINSIDE_SPARSE
#include <stdio.h>
#include "spConfig.h"
#include "ngspice/spmatrix.h"
#include "spConfig.h"
#include "spDefs.h"
#include <stdio.h>
/*
* Global strings
@ -83,20 +77,15 @@ char spcErrorsMustBeCleared[] = "Error not cleared";
char spcMatrixMustBeFactored[] = "Matrix must be factored";
char spcMatrixMustNotBeFactored[] = "Matrix must not be factored";
/*
* Function declarations
*/
//static spError ReserveElements( MatrixPtr, int );
static void InitializeElementBlocks( MatrixPtr, int, int );
static void RecordAllocation( MatrixPtr, void * );
static void AllocateBlockOfAllocationList( MatrixPtr );
// static spError ReserveElements( MatrixPtr, int );
static void InitializeElementBlocks(MatrixPtr, int, int);
static void RecordAllocation(MatrixPtr, void *);
static void AllocateBlockOfAllocationList(MatrixPtr);
/*!
* Allocates and initializes the data structures associated with a matrix.
*
@ -111,8 +100,8 @@ static void AllocateBlockOfAllocationList( MatrixPtr );
* \param Complex
* Type of matrix. If \a Complex is 0 then the matrix is real, otherwise
* the matrix will be complex. Note that if the routines are not set up
* to handle the type of matrix requested, then an \a spPANIC error will occur.
* Further note that if a matrix will be both real and complex, it must
* to handle the type of matrix requested, then an \a spPANIC error will
* occur. Further note that if a matrix will be both real and complex, it must
* be specified here as being complex.
* \param pError
* Returns error flag, needed because function \a spError() will
@ -126,43 +115,37 @@ static void AllocateBlockOfAllocationList( MatrixPtr );
* A pointer to the matrix frame being created.
*/
MatrixPtr
spCreate(
int Size,
int Complex,
int *pError
)
{
unsigned SizePlusOne;
MatrixPtr Matrix;
int I;
int AllocatedSize;
MatrixPtr spCreate(int Size, int Complex, int *pError) {
unsigned SizePlusOne;
MatrixPtr Matrix;
int I;
int AllocatedSize;
/* Begin `spCreate'. */
/* Clear error flag. */
/* Begin `spCreate'. */
/* Clear error flag. */
*pError = spOKAY;
/* Test for valid size. */
vASSERT( (Size >= 0) AND (Size != 0 OR EXPANDABLE), "Invalid size" );
/* Test for valid size. */
vASSERT((Size >= 0) AND(Size != 0 OR EXPANDABLE), "Invalid size");
/* Test for valid type. */
#if NOT spCOMPLEX
ASSERT( NOT Complex );
ASSERT(NOT Complex);
#endif
#if NOT REAL
ASSERT( Complex );
ASSERT(Complex);
#endif
/* Create Matrix. */
AllocatedSize = MAX( Size, MINIMUM_ALLOCATED_SIZE );
/* Create Matrix. */
AllocatedSize = MAX(Size, MINIMUM_ALLOCATED_SIZE);
SizePlusOne = (unsigned)(AllocatedSize + 1);
if ((Matrix = SP_MALLOC(struct MatrixFrame, 1)) == NULL)
{ *pError = spNO_MEMORY;
if ((Matrix = SP_MALLOC(struct MatrixFrame, 1)) == NULL) {
*pError = spNO_MEMORY;
return NULL;
}
/* Initialize matrix */
/* Initialize matrix */
Matrix->ID = SPARSE_ID;
Matrix->Complex = Complex;
Matrix->PreviousMatrixWasComplex = Complex;
@ -201,10 +184,11 @@ int AllocatedSize;
Matrix->ElementsRemaining = 0;
Matrix->FillinsRemaining = 0;
RecordAllocation( Matrix, Matrix );
if (Matrix->Error == spNO_MEMORY) goto MemoryError; /* FIXME: Use of memory after free */
RecordAllocation(Matrix, Matrix);
if (Matrix->Error == spNO_MEMORY)
goto MemoryError; /* FIXME: Use of memory after free */
/* Take out the trash. */
/* Take out the trash. */
Matrix->TrashCan.Real = 0.0;
#if spCOMPLEX
Matrix->TrashCan.Imag = 0.0;
@ -217,56 +201,56 @@ int AllocatedSize;
Matrix->TrashCan.pInitInfo = NULL;
#endif
/* Allocate space in memory for Diag pointer vector. */
SP_CALLOC( Matrix->Diag, ElementPtr, SizePlusOne);
/* Allocate space in memory for Diag pointer vector. */
SP_CALLOC(Matrix->Diag, ElementPtr, SizePlusOne);
if (Matrix->Diag == NULL)
goto MemoryError;
/* Allocate space in memory for FirstInCol pointer vector. */
SP_CALLOC( Matrix->FirstInCol, ElementPtr, SizePlusOne);
/* Allocate space in memory for FirstInCol pointer vector. */
SP_CALLOC(Matrix->FirstInCol, ElementPtr, SizePlusOne);
if (Matrix->FirstInCol == NULL)
goto MemoryError;
/* Allocate space in memory for FirstInRow pointer vector. */
SP_CALLOC( Matrix->FirstInRow, ElementPtr, SizePlusOne);
/* Allocate space in memory for FirstInRow pointer vector. */
SP_CALLOC(Matrix->FirstInRow, ElementPtr, SizePlusOne);
if (Matrix->FirstInRow == NULL)
goto MemoryError;
/* Allocate space in memory for IntToExtColMap vector. */
if (( Matrix->IntToExtColMap = SP_MALLOC(int, SizePlusOne)) == NULL)
/* Allocate space in memory for IntToExtColMap vector. */
if ((Matrix->IntToExtColMap = SP_MALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Allocate space in memory for IntToExtRowMap vector. */
if (( Matrix->IntToExtRowMap = SP_MALLOC(int, SizePlusOne)) == NULL)
/* Allocate space in memory for IntToExtRowMap vector. */
if ((Matrix->IntToExtRowMap = SP_MALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Initialize MapIntToExt vectors. */
for (I = 1; I <= AllocatedSize; I++)
{ Matrix->IntToExtRowMap[I] = I;
/* Initialize MapIntToExt vectors. */
for (I = 1; I <= AllocatedSize; I++) {
Matrix->IntToExtRowMap[I] = I;
Matrix->IntToExtColMap[I] = I;
}
#if TRANSLATE
/* Allocate space in memory for ExtToIntColMap vector. */
if (( Matrix->ExtToIntColMap = SP_MALLOC(int, SizePlusOne)) == NULL)
/* Allocate space in memory for ExtToIntColMap vector. */
if ((Matrix->ExtToIntColMap = SP_MALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Allocate space in memory for ExtToIntRowMap vector. */
if (( Matrix->ExtToIntRowMap = SP_MALLOC(int, SizePlusOne)) == NULL)
/* Allocate space in memory for ExtToIntRowMap vector. */
if ((Matrix->ExtToIntRowMap = SP_MALLOC(int, SizePlusOne)) == NULL)
goto MemoryError;
/* Initialize MapExtToInt vectors. */
for (I = 1; I <= AllocatedSize; I++)
{ Matrix->ExtToIntColMap[I] = -1;
/* Initialize MapExtToInt vectors. */
for (I = 1; I <= AllocatedSize; I++) {
Matrix->ExtToIntColMap[I] = -1;
Matrix->ExtToIntRowMap[I] = -1;
}
Matrix->ExtToIntColMap[0] = 0;
Matrix->ExtToIntRowMap[0] = 0;
#endif
/* Allocate space for fill-ins and initial set of elements. */
InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize,
SPACE_FOR_FILL_INS*AllocatedSize );
/* Allocate space for fill-ins and initial set of elements. */
InitializeElementBlocks(Matrix, SPACE_FOR_ELEMENTS * AllocatedSize,
SPACE_FOR_FILL_INS * AllocatedSize);
if (Matrix->Error == spNO_MEMORY)
goto MemoryError;
@ -274,21 +258,13 @@ int AllocatedSize;
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( Matrix );
spDestroy(Matrix);
return NULL;
}
/*
* ELEMENT ALLOCATION
*
@ -312,34 +288,26 @@ MemoryError:
* spNO_MEMORY
*/
ElementPtr
spcGetElement( MatrixPtr Matrix )
{
ElementPtr pElement;
ElementPtr spcGetElement(MatrixPtr Matrix) {
ElementPtr pElement;
/* Begin `spcGetElement'. */
/* Begin `spcGetElement'. */
/* Allocate block of MatrixElements if necessary. */
if (Matrix->ElementsRemaining == 0)
{ pElement = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, pElement );
if (Matrix->Error == spNO_MEMORY) return NULL;
/* Allocate block of MatrixElements if necessary. */
if (Matrix->ElementsRemaining == 0) {
pElement = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation(Matrix, pElement);
if (Matrix->Error == spNO_MEMORY)
return NULL;
Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailElement = pElement;
}
/* Update Element counter and return pointer to Element. */
/* Update Element counter and return pointer to Element. */
Matrix->ElementsRemaining--;
return Matrix->NextAvailElement++;
}
/*
* ELEMENT ALLOCATION INITIALIZATION
*
@ -370,53 +338,44 @@ ElementPtr pElement;
* spNO_MEMORY
*/
static void
InitializeElementBlocks(
MatrixPtr Matrix,
int InitialNumberOfElements,
int NumberOfFillinsExpected
)
{
ElementPtr pElement;
static void InitializeElementBlocks(MatrixPtr Matrix,
int InitialNumberOfElements,
int NumberOfFillinsExpected) {
ElementPtr pElement;
/* Begin `InitializeElementBlocks'. */
/* Begin `InitializeElementBlocks'. */
/* Allocate block of MatrixElements for elements. */
/* Allocate block of MatrixElements for elements. */
pElement = SP_MALLOC(struct MatrixElement, InitialNumberOfElements);
RecordAllocation( Matrix, pElement );
if (Matrix->Error == spNO_MEMORY) return;
RecordAllocation(Matrix, pElement);
if (Matrix->Error == spNO_MEMORY)
return;
Matrix->ElementsRemaining = InitialNumberOfElements;
Matrix->NextAvailElement = pElement;
/* Allocate block of MatrixElements for fill-ins. */
/* Allocate block of MatrixElements for fill-ins. */
pElement = SP_MALLOC(struct MatrixElement, NumberOfFillinsExpected);
RecordAllocation( Matrix, pElement );
if (Matrix->Error == spNO_MEMORY) return;
RecordAllocation(Matrix, pElement);
if (Matrix->Error == spNO_MEMORY)
return;
Matrix->FillinsRemaining = NumberOfFillinsExpected;
Matrix->NextAvailFillin = pElement;
/* Allocate a fill-in list structure. */
Matrix->FirstFillinListNode = SP_MALLOC(struct FillinListNodeStruct,1);
RecordAllocation( Matrix, Matrix->FirstFillinListNode );
if (Matrix->Error == spNO_MEMORY) return;
/* Allocate a fill-in list structure. */
Matrix->FirstFillinListNode = SP_MALLOC(struct FillinListNodeStruct, 1);
RecordAllocation(Matrix, Matrix->FirstFillinListNode);
if (Matrix->Error == spNO_MEMORY)
return;
Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
Matrix->FirstFillinListNode->pFillinList = pElement;
Matrix->FirstFillinListNode->NumberOfFillinsInList =NumberOfFillinsExpected;
Matrix->FirstFillinListNode->NumberOfFillinsInList =
NumberOfFillinsExpected;
Matrix->FirstFillinListNode->Next = NULL;
return;
}
/*
* FILL-IN ALLOCATION
*
@ -436,44 +395,42 @@ ElementPtr pElement;
* spNO_MEMORY
*/
ElementPtr
spcGetFillin( MatrixPtr Matrix )
{
ElementPtr spcGetFillin(MatrixPtr Matrix) {
#if STRIP OR LINT
struct FillinListNodeStruct *pListNode;
ElementPtr pFillins;
struct FillinListNodeStruct *pListNode;
ElementPtr pFillins;
#endif
/* Begin `spcGetFillin'. */
/* Begin `spcGetFillin'. */
#if NOT STRIP OR LINT
if (Matrix->FillinsRemaining == 0)
return spcGetElement( Matrix );
return spcGetElement(Matrix);
#endif
#if STRIP OR 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 = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, pFillins );
if (Matrix->Error == spNO_MEMORY) return NULL;
RecordAllocation(Matrix, pFillins);
if (Matrix->Error == spNO_MEMORY)
return NULL;
Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailFillin = pFillins;
/* Allocate a fill-in list structure. */
pListNode->Next = SP_MALLOC(struct FillinListNodeStruct,1);
RecordAllocation( Matrix, pListNode->Next );
if (Matrix->Error == spNO_MEMORY) return NULL;
/* Allocate a fill-in list structure. */
pListNode->Next = SP_MALLOC(struct FillinListNodeStruct, 1);
RecordAllocation(Matrix, pListNode->Next);
if (Matrix->Error == spNO_MEMORY)
return NULL;
Matrix->LastFillinListNode = pListNode = pListNode->Next;
pListNode->pFillinList = pFillins;
@ -483,19 +440,11 @@ ElementPtr pFillins;
}
#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++;
}
/*
* RECORD A MEMORY ALLOCATION
*
@ -505,7 +454,7 @@ ElementPtr pFillins;
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to the matrix.
* AllocatedPtr <input>
* AllocatedPtr <input>
* The pointer returned by malloc or calloc. These pointers are saved in
* a list so that they can be easily freed.
*
@ -513,45 +462,32 @@ ElementPtr pFillins;
* spNO_MEMORY
*/
static void
RecordAllocation(
MatrixPtr Matrix,
void *AllocatedPtr
)
{
/* Begin `RecordAllocation'. */
/*
* If Allocated pointer is NULL, assume that malloc returned a NULL pointer,
* which indicates a spNO_MEMORY error.
*/
if (AllocatedPtr == NULL)
{ Matrix->Error = spNO_MEMORY;
static void RecordAllocation(MatrixPtr Matrix, void *AllocatedPtr) {
/* Begin `RecordAllocation'. */
/*
* If Allocated pointer is NULL, assume that malloc returned a NULL pointer,
* which indicates a spNO_MEMORY error.
*/
if (AllocatedPtr == NULL) {
Matrix->Error = spNO_MEMORY;
return;
}
/* Allocate block of MatrixElements if necessary. */
if (Matrix->RecordsRemaining == 0)
{ AllocateBlockOfAllocationList( Matrix );
if (Matrix->Error == spNO_MEMORY)
{ SP_FREE(AllocatedPtr);
/* Allocate block of MatrixElements if necessary. */
if (Matrix->RecordsRemaining == 0) {
AllocateBlockOfAllocationList(Matrix);
if (Matrix->Error == spNO_MEMORY) {
SP_FREE(AllocatedPtr);
return;
}
}
/* Add Allocated pointer to Allocation List. */
/* Add Allocated pointer to Allocation List. */
(++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr;
Matrix->RecordsRemaining--;
return;
}
/*
* ADD A BLOCK OF SLOTS TO ALLOCATION LIST
*
@ -570,45 +506,36 @@ RecordAllocation(
* spNO_MEMORY
*/
static void
AllocateBlockOfAllocationList( MatrixPtr Matrix )
{
int I;
AllocationListPtr ListPtr;
static void AllocateBlockOfAllocationList(MatrixPtr Matrix) {
int I;
AllocationListPtr ListPtr;
/* Begin `AllocateBlockOfAllocationList'. */
/* Allocate block of records for allocation list. */
ListPtr = SP_MALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION+1));
if (ListPtr == NULL)
{ Matrix->Error = spNO_MEMORY;
/* Begin `AllocateBlockOfAllocationList'. */
/* Allocate block of records for allocation list. */
ListPtr = SP_MALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION + 1));
if (ListPtr == NULL) {
Matrix->Error = spNO_MEMORY;
return;
}
/* String entries of allocation list into singly linked list. List is linked
such that any record points to the one before it. */
/* 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. */
/* Record allocation of space for allocation list on allocation list. */
Matrix->TopOfAllocationList->AllocatedPtr = (void *)ListPtr;
Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;
return;
}
/*!
* Destroys a matrix and frees all memory associated with it.
*
@ -626,50 +553,42 @@ AllocationListPtr ListPtr;
* in the allocation list.
*/
void
spDestroy( MatrixPtr Matrix )
{
AllocationListPtr ListPtr, NextListPtr;
void spDestroy(MatrixPtr Matrix) {
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. */
SP_FREE( Matrix->IntToExtColMap );
SP_FREE( Matrix->IntToExtRowMap );
SP_FREE( Matrix->ExtToIntColMap );
SP_FREE( Matrix->ExtToIntRowMap );
SP_FREE( Matrix->Diag );
SP_FREE( Matrix->FirstInRow );
SP_FREE( Matrix->FirstInCol );
SP_FREE( Matrix->MarkowitzRow );
SP_FREE( Matrix->MarkowitzCol );
SP_FREE( Matrix->MarkowitzProd );
SP_FREE( Matrix->DoCmplxDirect );
SP_FREE( Matrix->DoRealDirect );
SP_FREE( Matrix->Intermediate );
/* Deallocate the vectors that are located in the matrix frame. */
SP_FREE(Matrix->IntToExtColMap);
SP_FREE(Matrix->IntToExtRowMap);
SP_FREE(Matrix->ExtToIntColMap);
SP_FREE(Matrix->ExtToIntRowMap);
SP_FREE(Matrix->Diag);
SP_FREE(Matrix->FirstInRow);
SP_FREE(Matrix->FirstInCol);
SP_FREE(Matrix->MarkowitzRow);
SP_FREE(Matrix->MarkowitzCol);
SP_FREE(Matrix->MarkowitzProd);
SP_FREE(Matrix->DoCmplxDirect);
SP_FREE(Matrix->DoRealDirect);
SP_FREE(Matrix->Intermediate);
/* Sequentially step through the list of allocated pointers freeing pointers
* along the way. */
/* Sequentially step through the list of allocated pointers freeing pointers
* along the way. */
ListPtr = Matrix->TopOfAllocationList;
while (ListPtr != NULL)
{ NextListPtr = ListPtr->NextRecord;
if ((void *) ListPtr == ListPtr->AllocatedPtr) {
SP_FREE( ListPtr );
while (ListPtr != NULL) {
NextListPtr = ListPtr->NextRecord;
if ((void *)ListPtr == ListPtr->AllocatedPtr) {
SP_FREE(ListPtr);
} else {
SP_FREE( ListPtr->AllocatedPtr );
SP_FREE(ListPtr->AllocatedPtr);
}
ListPtr = NextListPtr;
}
return;
}
/*!
* This function returns the error status of the given matrix.
*
@ -680,27 +599,17 @@ AllocationListPtr ListPtr, NextListPtr;
* The pointer to the matrix for which the error status is desired.
*/
int
spError( MatrixPtr Matrix )
{
/* Begin `spError'. */
int spError(MatrixPtr Matrix) {
/* Begin `spError'. */
if (Matrix != NULL)
{ ASSERT_IS_SPARSE( Matrix );
if (Matrix != NULL) {
ASSERT_IS_SPARSE(Matrix);
return Matrix->Error;
}
else return spNO_MEMORY; /* This error may actually be spPANIC,
* no way to tell. */
} else
return spNO_MEMORY; /* This error may actually be spPANIC,
* no way to tell. */
}
/*!
* This function returns the row and column number where the matrix was
* detected as singular (if pivoting was allowed on the last factorization)
@ -716,29 +625,18 @@ spError( MatrixPtr Matrix )
* The column number.
*/
void
spWhereSingular(
MatrixPtr Matrix,
int *pRow,
int *pCol
)
{
/* Begin `spWhereSingular'. */
ASSERT_IS_SPARSE( Matrix );
void spWhereSingular(MatrixPtr Matrix, int *pRow, int *pCol) {
/* Begin `spWhereSingular'. */
ASSERT_IS_SPARSE(Matrix);
if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG)
{ *pRow = Matrix->SingularRow;
if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG) {
*pRow = Matrix->SingularRow;
*pCol = Matrix->SingularCol;
}
else *pRow = *pCol = 0;
} else
*pRow = *pCol = 0;
return;
}
/*!
* Returns the size of the matrix. Either the internal or external size of
* the matrix is returned.
@ -752,14 +650,9 @@ spWhereSingular(
* may differ if the \a TRANSLATE option is set true.
*/
int
spGetSize(
MatrixPtr Matrix,
int External
)
{
/* Begin `spGetSize'. */
ASSERT_IS_SPARSE( Matrix );
int spGetSize(MatrixPtr Matrix, int External) {
/* Begin `spGetSize'. */
ASSERT_IS_SPARSE(Matrix);
#if TRANSLATE
if (External)
@ -771,13 +664,6 @@ spGetSize(
#endif
}
/*!
* Forces matrix to be real.
*
@ -785,18 +671,15 @@ spGetSize(
* Pointer to matrix.
*/
void
spSetReal( MatrixPtr Matrix )
{
/* Begin `spSetReal'. */
void spSetReal(MatrixPtr Matrix) {
/* Begin `spSetReal'. */
ASSERT_IS_SPARSE( Matrix );
vASSERT( REAL, "Sparse not compiled to handle real matrices" );
ASSERT_IS_SPARSE(Matrix);
vASSERT(REAL, "Sparse not compiled to handle real matrices");
Matrix->Complex = NO;
return;
}
/*!
* Forces matrix to be complex.
*
@ -804,44 +687,33 @@ spSetReal( MatrixPtr Matrix )
* Pointer to matrix.
*/
void
spSetComplex( MatrixPtr Matrix )
{
/* Begin `spSetComplex'. */
void spSetComplex(MatrixPtr Matrix) {
/* Begin `spSetComplex'. */
ASSERT_IS_SPARSE( Matrix );
vASSERT( spCOMPLEX, "Sparse not compiled to handle complex matrices" );
ASSERT_IS_SPARSE(Matrix);
vASSERT(spCOMPLEX, "Sparse not compiled to handle complex matrices");
Matrix->Complex = YES;
return;
}
/*!
* This function returns the number of fill-ins that currently exists in a matrix.
* This function returns the number of fill-ins that currently exists in a
* matrix.
*
* \param Matrix
* Pointer to matrix.
*/
int
spFillinCount( MatrixPtr Matrix )
{
/* Begin `spFillinCount'. */
int spFillinCount(MatrixPtr Matrix) {
/* Begin `spFillinCount'. */
ASSERT_IS_SPARSE( Matrix );
ASSERT_IS_SPARSE(Matrix);
return Matrix->Fillins;
}
/*!
* This function returns the total number of elements (including fill-ins) that currently exists in a matrix.
* This function returns the total number of elements (including fill-ins) that
* currently exists in a matrix.
*
* \param Matrix
* Pointer to matrix.
@ -849,21 +721,16 @@ spFillinCount( MatrixPtr Matrix )
/* FIXME: Seems no different size entries available anymore */
int
spElementCount( MatrixPtr Matrix )
{
/* Begin `spElementCount'. */
int spElementCount(MatrixPtr Matrix) {
/* Begin `spElementCount'. */
ASSERT_IS_SPARSE( Matrix );
ASSERT_IS_SPARSE(Matrix);
return Matrix->Elements;
}
int
spOriginalCount( MatrixPtr Matrix )
{
int spOriginalCount(MatrixPtr Matrix) {
/* Begin `spOriginalCount'. */
ASSERT_IS_SPARSE( Matrix );
ASSERT_IS_SPARSE(Matrix);
return Matrix->Elements;
}

957
src/maths/sparse/spBuild.c
View File

File diff suppressed because it is too large Load Diff

178
src/maths/sparse/spConfig.h
View File

@ -19,7 +19,6 @@
* Kenneth S. Kundert <kundert@users.sourceforge.net>
*/
/*
* Revision and copyright information.
*
@ -29,13 +28,9 @@
* $Revision: 1.5 $
*/
#ifndef spCONFIG_DEFS
#define spCONFIG_DEFS
#ifdef spINSIDE_SPARSE
/*
* OPTIONS
@ -67,7 +62,7 @@
* slight speed and memory advantage if the routines are complied
* to handle only real systems of equations.
*/
#define REAL YES
#define REAL YES
/*!
* Setting this compiler flag true (1) makes the matrix
@ -78,7 +73,7 @@
* the size of the matrix need not be known before the matrix is
* built. The matrix can be allocated with size zero and expanded.
*/
#define EXPANDABLE YES
#define EXPANDABLE YES
/*!
* This option allows the set of external row and column numbers
@ -98,7 +93,7 @@
* vectors must be at least as large as the external size, which
* is the value of the largest given row or column numbers.
*/
#define TRANSLATE YES
#define TRANSLATE YES
/*!
* Causes the spInitialize(), spGetInitInfo(), and
@ -110,7 +105,7 @@
* column numbers as arguments. This allows the user to write
* custom matrix initialization routines.
*/
#define INITIALIZE NO
#define INITIALIZE NO
/*!
* Many matrices, and in particular node- and modified-node
@ -136,7 +131,7 @@
* options and constants are not used: \a MODIFIED_MARKOWITZ,
* \a MAX_MARKOWITZ_TIES, and \a TIES_MULTIPLIER.
*/
#define DIAGONAL_PIVOTING YES
#define DIAGONAL_PIVOTING YES
/*!
* This determines whether arrays start at an index of zero or one.
@ -151,7 +146,7 @@
* matrix. ARRAY_OFFSET must be either 0 or 1, no other offsets
* are valid.
*/
#define ARRAY_OFFSET YES
#define ARRAY_OFFSET YES
/*!
* This specifies that the modified Markowitz method of pivot
@ -173,19 +168,19 @@
* for use when working with very large matrices where the initial
* factor time represents an unacceptable burden. \a NO is recommended.
*/
#define MODIFIED_MARKOWITZ NO
#define MODIFIED_MARKOWITZ NO
/*!
* This specifies that the spDeleteRowAndCol() routine
* should be compiled. Note that for this routine to be
* compiled, both \a DELETE and \a TRANSLATE should be set true.
*/
#define DELETE NO
#define DELETE NO
/*!
* This specifies that the spStripFills() routine should be compiled.
*/
#define STRIP NO
#define STRIP NO
/*!
* This specifies that the routine that preorders modified node
@ -193,7 +188,7 @@
* in greater speed and accuracy if used with this type of
* matrix.
*/
#define MODIFIED_NODAL YES
#define MODIFIED_NODAL YES
/*!
* This specifies that the routines that allow four related
@ -202,7 +197,7 @@
* admittance. The routines affected by \a QUAD_ELEMENT are the
* spGetAdmittance(), spGetQuad() and spGetOnes() routines.
*/
#define QUAD_ELEMENT NO
#define QUAD_ELEMENT NO
/*!
* This specifies that the routines that solve the matrix as if
@ -210,7 +205,7 @@
* useful when performing sensitivity analysis using the adjoint
* method.
*/
#define TRANSPOSE YES
#define TRANSPOSE YES
/*!
* This specifies that the routine that performs scaling on the
@ -224,26 +219,26 @@
* solution by the user or the scaled factors may simply be
* thrown away. \a NO is recommended.
*/
#define SCALING NO
#define SCALING NO
/*!
* This specifies that routines that are used to document the
* matrix, such as spPrint() and spFileMatrix(), should be
* compiled.
*/
#define DOCUMENTATION YES
#define DOCUMENTATION YES
/*!
* This specifies that routines that are used to multiply the
* matrix by a vector, such as spMultiply() and spMultTransposed(), should be
* compiled.
*/
#define MULTIPLICATION YES
#define MULTIPLICATION YES
/*!
* This specifies that the routine spDeterminant() should be complied.
*/
#define DETERMINANT YES
#define DETERMINANT YES
/*!
* This specifies that spLargestElement() and spRoundoff() should
@ -255,28 +250,28 @@
* If the bound increases greatly after using spFactor(), then the
* matrix should probably be reordered. Recomend \a NO.
*/
#define STABILITY NO
#define STABILITY NO
/*!
* This specifies that spCondition() and spNorm(), the code that
* computes a good estimate of the condition number of the matrix,
* should be compiled. Recomend \a NO.
*/
#define CONDITION NO
#define CONDITION NO
/*!
* This specifies that spPseudoCondition(), the code that computes
* a crude and easily fooled indicator of ill-conditioning in the
* matrix, should be compiled. Recomend \a NO.
*/
#define PSEUDOCONDITION NO
#define PSEUDOCONDITION NO
/*!
* This specifies that the \a FORTRAN interface routines should be
* compiled. When interfacing to \a FORTRAN programs, the \a ARRAY_OFFSET
* options should be set to NO.
*/
#define FORTRAN NO
#define FORTRAN NO
/*!
* This specifies that additional error checking will be compiled.
@ -285,7 +280,7 @@
* the routines have been integrated in and are running smoothly, this
* option should be turned off. \a YES is recommended.
*/
#define DEBUG YES
#define DEBUG YES
#endif /* spINSIDE_SPARSE */
@ -300,7 +295,7 @@
* This specifies that the routines will be complied to handle
* complex systems of equations.
*/
#define spCOMPLEX 1
#define spCOMPLEX 1
/*!
* This specifies the format for complex vectors. If this is set
@ -313,16 +308,10 @@
* vector is represented by two arrays of \a spREALs, one with
* the real terms, the other with the imaginary. \a NO is recommended.
*/
#define spSEPARATED_COMPLEX_VECTORS 1
#define spSEPARATED_COMPLEX_VECTORS 1
#ifdef spINSIDE_SPARSE
/*
* MATRIX CONSTANTS
*
@ -339,28 +328,28 @@
* not be less than or equal to zero nor larger than one.
* 0.001 is recommended.
*/
#define DEFAULT_THRESHOLD 1.0e-3
#define DEFAULT_THRESHOLD 1.0e-3
/*!
* This indicates whether spOrderAndFactor() should use diagonal
* pivoting as default. This issue only arises when
* spOrderAndFactor() is called from spFactor(). \a YES is recommended.
*/
#define DIAG_PIVOTING_AS_DEFAULT YES
#define DIAG_PIVOTING_AS_DEFAULT YES
/*!
* This number multiplied by the size of the matrix equals the number
* of elements for which memory is initially allocated in spCreate().
* 6 is recommended.
*/
#define SPACE_FOR_ELEMENTS 6
#define SPACE_FOR_ELEMENTS 6
/*!
* This number multiplied by the size of the matrix equals the number
* of elements for which memory is initially allocated and specifically
* reserved for fill-ins in spCreate(). 4 is recommended.
*/
#define SPACE_FOR_FILL_INS 4
#define SPACE_FOR_FILL_INS 4
/*!
* The number of matrix elements requested from the malloc utility on
@ -370,7 +359,7 @@
* elements at a time (or some multiple thereof).
* 31 is recommended.
*/
#define ELEMENTS_PER_ALLOCATION 31
#define ELEMENTS_PER_ALLOCATION 31
/*!
* The minimum allocated size of a matrix. Note that this does not
@ -379,13 +368,13 @@
* allocated with an estimated size of zero. This number should not
* be less than one.
*/
#define MINIMUM_ALLOCATED_SIZE 6
#define MINIMUM_ALLOCATED_SIZE 6
/*!
* The amount the allocated size of the matrix is increased when it
* is expanded.
*/
#define EXPANSION_FACTOR 1.5
#define EXPANSION_FACTOR 1.5
/*!
* Some terminology should be defined. The Markowitz row count is the number
@ -412,7 +401,7 @@
* 100 is recommended.
* \see TIES_MULTIPLIER
*/
#define MAX_MARKOWITZ_TIES 100
#define MAX_MARKOWITZ_TIES 100
/*!
* Specifies the number of Markowitz ties that are allowed to occur
@ -429,7 +418,7 @@
* diagonal pivoting breaks down. 5 is recommended.
* \see MAX_MARKOWITZ_TIES
*/
#define TIES_MULTIPLIER 5
#define TIES_MULTIPLIER 5
/*!
* Which partition mode is used by spPartition() as default.
@ -441,77 +430,61 @@
* overhead, but speeds up both dense and sparse matrices, best if there
* is a large number of matrices that can use the same ordering.
*/
#define DEFAULT_PARTITION spAUTO_PARTITION
#define DEFAULT_PARTITION spAUTO_PARTITION
/*!
* The number of characters per page width. Set to 80 for terminal,
* 132 for line printer. Controls how many columns printed by
* spPrint() per page width.
*/
#define PRINTER_WIDTH 80
#define PRINTER_WIDTH 80
#endif /* spINSIDE_SPARSE */
/*
* PORTABILITY MACROS
*/
#ifdef __STDC__
# define spcCONCAT(prefix,suffix) prefix ## suffix
# define spcQUOTE(x) # x
# define spcFUNC_NEEDS_FILE(func,file) \
func ## _requires_ ## file ## _to_be_included_
#define spcCONCAT(prefix, suffix) prefix##suffix
#define spcQUOTE(x) #x
#define spcFUNC_NEEDS_FILE(func, file) func##_requires_##file##_to_be_included_
#else
# define spcCONCAT(prefix,suffix) prefix/**/suffix
# define spcQUOTE(x) "x"
# define spcFUNC_NEEDS_FILE(func,file) \
func/**/_requires_/**/file/**/_to_be_included_
#define spcCONCAT(prefix, suffix) prefix /**/ suffix
#define spcQUOTE(x) "x"
#define spcFUNC_NEEDS_FILE(func, file) \
func /**/ _requires_ /**/ file /**/ _to_be_included_
#endif
#if defined(__cplusplus) || defined(c_plusplus)
/*
* Definitions for C++
*/
# define spcEXTERN extern "C"
# define spcNO_ARGS
# define spcCONST const
typedef void *spGenericPtr;
/*
* Definitions for C++
*/
#define spcEXTERN extern "C"
#define spcNO_ARGS
#define spcCONST const
typedef void *spGenericPtr;
#else
#ifdef __STDC__
/*
* Definitions for ANSI C
*/
# define spcEXTERN extern
# define spcNO_ARGS void
# define spcCONST const
typedef void *spGenericPtr;
# else
/*
* Definitions for K&R C -- ignore function prototypes
*/
# define spcEXTERN extern
# define spcNO_ARGS
# define spcCONST
typedef char *spGenericPtr;
/*
* Definitions for ANSI C
*/
#define spcEXTERN extern
#define spcNO_ARGS void
#define spcCONST const
typedef void *spGenericPtr;
#else
/*
* Definitions for K&R C -- ignore function prototypes
*/
#define spcEXTERN extern
#define spcNO_ARGS
#define spcCONST
typedef char *spGenericPtr;
#endif
#endif
#ifdef spINSIDE_SPARSE
/*
* MACHINE CONSTANTS
*
@ -519,29 +492,24 @@
*/
/* Begin machine constants. */
#include <limits.h>
#include <float.h>
#include <limits.h>
/*! The resolution of spREAL. */
#define MACHINE_RESOLUTION DBL_EPSILON
#define MACHINE_RESOLUTION DBL_EPSILON
/*! The largest possible value of spREAL. */
#define LARGEST_REAL DBL_MAX
#define LARGEST_REAL DBL_MAX
/*! The smalles possible positive value of spREAL. */
#define SMALLEST_REAL DBL_MIN
#define SMALLEST_REAL DBL_MIN
/*! The largest possible value of shorts. */
#define LARGEST_SHORT_INTEGER SHRT_MAX
#define LARGEST_SHORT_INTEGER SHRT_MAX
/*! The largest possible value of longs. */
#define LARGEST_LONG_INTEGER LONG_MAX
#define LARGEST_LONG_INTEGER LONG_MAX
/* ANNOTATION */
/*!
* This macro changes the amount of annotation produced by the matrix
@ -550,23 +518,23 @@
* the program. Possible values include \a NONE, \a ON_STRANGE_BEHAVIOR, and
* \a FULL. \a NONE is recommended.
*/
#define ANNOTATE NONE
#define ANNOTATE NONE
/*!
* A possible value for \a ANNOTATE. Disables all annotation.
*/
#define NONE 0
#define NONE 0
/*!
* A possible value for \a ANNOTATE. Causes annotation to be produce
* upon unusual occurances only.
*/
#define ON_STRANGE_BEHAVIOR 1
#define ON_STRANGE_BEHAVIOR 1
/*!
* A possible value for \a ANNOTATE. Enables full annotation.
*/
#define FULL 2
#define FULL 2
#endif /* spINSIDE_SPARSE */
#endif /* spCONFIG_DEFS */

834
src/maths/sparse/spDefs.h
View File

File diff suppressed because it is too large Load Diff

23
src/maths/sparse/spExtra.c
View File

@ -9,7 +9,6 @@
* >>> Other functions contained in this file:
*/
/*
* IMPORTS
*
@ -23,25 +22,21 @@
*/
#define spINSIDE_SPARSE
#include <stdio.h>
#include "spConfig.h"
#include "ngspice/spmatrix.h"
#include "spConfig.h"
#include "spDefs.h"
#include <stdio.h>
void
spConstMult(
MatrixPtr Matrix,
double constant
)
{
ElementPtr pElement;
int I;
int size = Matrix->Size;
void spConstMult(MatrixPtr Matrix, double constant) {
ElementPtr pElement;
int I;
int size = Matrix->Size;
ASSERT_IS_SPARSE( Matrix );
ASSERT_IS_SPARSE(Matrix);
for (I = 1; I <= size; I++) {
for (pElement = Matrix->FirstInCol[I]; pElement; pElement = pElement->NextInCol) {
for (pElement = Matrix->FirstInCol[I]; pElement;
pElement = pElement->NextInCol) {
pElement->Real *= constant;
#if spCOMPLEX
pElement->Imag *= constant;

2228
src/maths/sparse/spFactor.c
View File

File diff suppressed because it is too large Load Diff

541
src/maths/sparse/spOutput.c
View File

@ -25,7 +25,6 @@
* >>> Other functions contained in this file:
*/
/*
* Revision and copyright information.
*
@ -36,13 +35,10 @@
#ifdef notdef
static char copyright[] =
"Sparse1.4: Copyright (c) 1985-2003 by Kenneth S. Kundert";
static char RCSid[] =
"$Header: /cvsroot/sparse/src/spOutput.c,v 1.3 2003/06/29 04:19:52 kundert Exp $";
static char RCSid[] = "$Header: /cvsroot/sparse/src/spOutput.c,v 1.3 "
"2003/06/29 04:19:52 kundert Exp $";
#endif
/*
* IMPORTS
*
@ -56,17 +52,13 @@ static char RCSid[] =
*/
#define spINSIDE_SPARSE
#include <stdio.h>
#include "spConfig.h"
#include "ngspice/spmatrix.h"
#include "spConfig.h"
#include "spDefs.h"
#include <stdio.h>
#if DOCUMENTATION
/*!
* Formats and send the matrix to standard output. Some elementary
* statistics are also output. The matrix is output in a format that is
@ -134,66 +126,60 @@ static char RCSid[] =
* The largest expected external row or column number.
*/
void
spPrint(
MatrixPtr Matrix,
int PrintReordered,
int Data,
int Header
)
{
int J = 0;
int I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns, ElementCount = 0;
double Magnitude, SmallestDiag = 0.0, SmallestElement = 0.0;
double LargestElement = 0.0, LargestDiag = 0.0;
ElementPtr pElement;
void spPrint(MatrixPtr Matrix, int PrintReordered, int Data, int Header) {
int J = 0;
int I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns,
ElementCount = 0;
double Magnitude, SmallestDiag = 0.0, SmallestElement = 0.0;
double LargestElement = 0.0, LargestDiag = 0.0;
ElementPtr pElement;
#if spCOMPLEX
ElementPtr pImagElements[PRINTER_WIDTH/10+1];
ElementPtr pImagElements[PRINTER_WIDTH / 10 + 1];
#endif
int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
/* Begin `spPrint'. */
ASSERT_IS_SPARSE( Matrix );
/* Begin `spPrint'. */
ASSERT_IS_SPARSE(Matrix);
Size = Matrix->Size;
/* Create a packed external to internal row and column translation array. */
# if TRANSLATE
#if TRANSLATE
Top = Matrix->AllocatedExtSize;
#else
Top = Matrix->AllocatedSize;
#endif
SP_CALLOC( PrintOrdToIntRowMap, int, Top + 1 );
SP_CALLOC( PrintOrdToIntColMap, int, Top + 1 );
if ( PrintOrdToIntRowMap == NULL OR PrintOrdToIntColMap == NULL)
{ Matrix->Error = spNO_MEMORY;
SP_CALLOC(PrintOrdToIntRowMap, int, Top + 1);
SP_CALLOC(PrintOrdToIntColMap, int, Top + 1);
if (PrintOrdToIntRowMap == NULL OR PrintOrdToIntColMap == NULL) {
Matrix->Error = spNO_MEMORY;
SP_FREE(PrintOrdToIntColMap);
SP_FREE(PrintOrdToIntRowMap);
return;
}
for (I = 1; I <= Size; I++)
{ PrintOrdToIntRowMap[ Matrix->IntToExtRowMap[I] ] = I;
PrintOrdToIntColMap[ Matrix->IntToExtColMap[I] ] = I;
for (I = 1; I <= Size; I++) {
PrintOrdToIntRowMap[Matrix->IntToExtRowMap[I]] = I;
PrintOrdToIntColMap[Matrix->IntToExtColMap[I]] = I;
}
/* Pack the arrays. */
for (J = 1, I = 1; I <= Top; I++)
{ if (PrintOrdToIntRowMap[I] != 0)
PrintOrdToIntRowMap[ J++ ] = PrintOrdToIntRowMap[ I ];
/* Pack the arrays. */
for (J = 1, I = 1; I <= Top; I++) {
if (PrintOrdToIntRowMap[I] != 0)
PrintOrdToIntRowMap[J++] = PrintOrdToIntRowMap[I];
}
for (J = 1, I = 1; I <= Top; I++)
{ if (PrintOrdToIntColMap[I] != 0)
PrintOrdToIntColMap[ J++ ] = PrintOrdToIntColMap[ I ];
for (J = 1, I = 1; I <= Top; I++) {
if (PrintOrdToIntColMap[I] != 0)
PrintOrdToIntColMap[J++] = PrintOrdToIntColMap[I];
}
/* Print header. */
if (Header)
{ printf("MATRIX SUMMARY\n\n");
/* Print header. */
if (Header) {
printf("MATRIX SUMMARY\n\n");
printf("Size of matrix = %1d x %1d.\n", Size, Size);
if ( Matrix->Reordered AND PrintReordered )
if (Matrix->Reordered AND PrintReordered)
printf("Matrix has been reordered.\n");
putchar('\n');
if ( Matrix->Factored )
if (Matrix->Factored)
printf("Matrix after factorization:\n");
else
printf("Matrix before factorization:\n");
@ -201,73 +187,78 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
SmallestElement = LARGEST_REAL;
SmallestDiag = SmallestElement;
}
if (Size == 0) return;
if (Size == 0)
return;
/* 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;
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 )
while (J <= Size)
/* Calculate index of last column to printed in this group. */
{ StopCol = StartCol + Columns - 1;
/* Calculate index of last column to printed in this group. */
{
StopCol = StartCol + Columns - 1;
if (StopCol > Size)
StopCol = Size;
/* Label the columns. */
if (Header)
{ if (Data)
{ printf(" ");
for (I = StartCol; I <= StopCol; I++)
{ if (PrintReordered)
/* Label the columns. */
if (Header) {
if (Data) {
printf(" ");
for (I = StartCol; I <= StopCol; I++) {
if (PrintReordered)
Col = I;
else
Col = PrintOrdToIntColMap[I];
printf(" %9d", Matrix->IntToExtColMap[ Col ]);
printf(" %9d", Matrix->IntToExtColMap[Col]);
}
printf("\n\n");
}
else
{ if (PrintReordered)
printf("Columns %1d to %1d.\n",StartCol,StopCol);
else
{ printf("Columns %1d to %1d.\n",
Matrix->IntToExtColMap[ PrintOrdToIntColMap[StartCol] ],
Matrix->IntToExtColMap[ PrintOrdToIntColMap[StopCol] ]);
} else {
if (PrintReordered)
printf("Columns %1d to %1d.\n", StartCol, StopCol);
else {
printf(
"Columns %1d to %1d.\n",
Matrix->IntToExtColMap[PrintOrdToIntColMap[StartCol]],
Matrix->IntToExtColMap[PrintOrdToIntColMap[StopCol]]);
}
}
}
/* Print every row ... */
for (I = 1; I <= Size; I++)
{ if (PrintReordered)
/* Print every row ... */
for (I = 1; I <= Size; I++) {
if (PrintReordered)
Row = I;
else
Row = PrintOrdToIntRowMap[I];
if (Header)
{ if (PrintReordered AND NOT Data)
if (Header) {
if (PrintReordered AND NOT Data)
printf("%4d", I);
else
printf("%4d", Matrix->IntToExtRowMap[ Row ]);
if (NOT Data) putchar(' ');
printf("%4d", Matrix->IntToExtRowMap[Row]);
if (NOT Data)
putchar(' ');
}
/* ... in each column of the group. */
for (J = StartCol; J <= StopCol; J++)
{ if (PrintReordered)
/* ... in each column of the group. */
for (J = StartCol; J <= StopCol; J++) {
if (PrintReordered)
Col = J;
else
Col = PrintOrdToIntColMap[J];
pElement = Matrix->FirstInCol[Col];
while(pElement != NULL AND pElement->Row != Row)
while (pElement != NULL AND pElement->Row != Row)
pElement = pElement->NextInCol;
#if spCOMPLEX
@ -277,23 +268,24 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
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 */
if ( (Magnitude = ELEMENT_MAG(pElement)) > LargestElement )
/* Update status variables */
if ((Magnitude = ELEMENT_MAG(pElement)) > LargestElement)
LargestElement = Magnitude;
if ((Magnitude < SmallestElement) AND (Magnitude != 0.0))
if ((Magnitude < SmallestElement) AND(Magnitude != 0.0))
SmallestElement = Magnitude;
ElementCount++;
}
/* Case where element is structurally zero */
else
{ if (Data)
/* Case where element is structurally zero */
else {
if (Data)
printf(" ...");
else
putchar('.');
@ -302,52 +294,53 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
putchar('\n');
#if spCOMPLEX
if (Matrix->Complex AND Data)
{ if (Header)
printf(" ");
for (J = StartCol; J <= StopCol; J++)
{ if (pImagElements[J - StartCol] != NULL)
{ printf(" %8.2gj",
(double)pImagElements[J-StartCol]->Imag);
}
else printf(" ");
if (Matrix->Complex AND Data) {
if (Header)
printf(" ");
for (J = StartCol; J <= StopCol; J++) {
if (pImagElements[J - StartCol] != NULL) {
printf(" %8.2gj",
(double)pImagElements[J - StartCol]->Imag);
} else
printf(" ");
}
putchar('\n');
}
#endif /* spCOMPLEX */
}
/* Calculate index of first column in next group. */
/* Calculate index of first column in next group. */
StartCol = StopCol;
StartCol++;
putchar('\n');
}
if (Header)
{ printf("\nLargest element in matrix = %-1.4g.\n", LargestElement);
if (Header) {
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 */
for (I = 1; I <= Size; I++)
{ if (Matrix->Diag[I] != NULL)
{ Magnitude = ELEMENT_MAG( Matrix->Diag[I] );
if ( Magnitude > LargestDiag ) LargestDiag = Magnitude;
if ( Magnitude < SmallestDiag ) SmallestDiag = Magnitude;
/* Search for largest and smallest diagonal values */
for (I = 1; I <= Size; I++) {
if (Matrix->Diag[I] != NULL) {
Magnitude = ELEMENT_MAG(Matrix->Diag[I]);
if (Magnitude > LargestDiag)
LargestDiag = Magnitude;
if (Magnitude < SmallestDiag)
SmallestDiag = Magnitude;
}
}
/* Print the largest and smallest diagonal values */
if ( Matrix->Factored )
{ printf("\nLargest diagonal element = %-1.4g.\n", LargestDiag);
/* Print the largest and smallest diagonal values */
if (Matrix->Factored) {
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);
} else {
printf("\nLargest pivot element = %-1.4g.\n", LargestDiag);
printf("Smallest pivot element = %-1.4g.\n", SmallestDiag);
}
/* Calculate and print sparsity and number of fill-ins created. */
printf("\nDensity = %2.2f%%.\n", ((double)ElementCount * 100.0)
/ (((double)Size * (double)Size)));
/* Calculate and print sparsity and number of fill-ins created. */
printf("\nDensity = %2.2f%%.\n", ((double)ElementCount * 100.0) /
(((double)Size * (double)Size)));
if (NOT Matrix->NeedsOrdering)
printf("Number of fill-ins = %1d.\n", Matrix->Fillins);
}
@ -359,16 +352,6 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
return;
}
/*!
* Writes matrix to file in format suitable to be read back in by the
* matrix test program.
@ -408,128 +391,116 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
* The size of the matrix.
*/
int
spFileMatrix(
MatrixPtr Matrix,
char *File,
char *Label,
int Reordered,
int Data,
int Header
)
{
int I, Size;
ElementPtr pElement;
int Row, Col, Err;
FILE *pMatrixFile;
int spFileMatrix(MatrixPtr Matrix, char *File, char *Label, int Reordered,
int Data, int Header) {
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 (Header) {
if (Matrix->Factored AND Data) {
Err = fprintf(
pMatrixFile,
"Warning : The following matrix is factored in to LU form.\n");
if (Err < 0)
return 0;
}
if (fprintf(pMatrixFile, "%s\n", Label) < 0) return 0;
Err = fprintf( pMatrixFile, "%d\t%s\n", Size,
(Matrix->Complex ? "complex" : "real"));
if (Err < 0) return 0;
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;
}
if (Size == 0) return 1;
if (Size == 0)
return 1;
/* Output matrix. */
if (NOT Data)
{ for (I = 1; I <= Size; I++)
{ pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{ if (Reordered)
{ Row = pElement->Row;
/* Output matrix. */
if (NOT Data) {
for (I = 1; I <= Size; I++) {
pElement = Matrix->FirstInCol[I];
while (pElement != NULL) {
if (Reordered) {
Row = pElement->Row;
Col = I;
}
else
{ Row = Matrix->IntToExtRowMap[pElement->Row];
} else {
Row = Matrix->IntToExtRowMap[pElement->Row];
Col = Matrix->IntToExtColMap[I];
}
pElement = pElement->NextInCol;
if (fprintf(pMatrixFile, "%d\t%d\n", Row, Col) < 0) return 0;
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 (fprintf(pMatrixFile, "0\t0\n") < 0)
return 0;
}
#if spCOMPLEX
if (Data AND Matrix->Complex)
{ for (I = 1; I <= Size; I++)
{ pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{ if (Reordered)
{ Row = pElement->Row;
if (Data AND Matrix->Complex) {
for (I = 1; I <= Size; I++) {
pElement = Matrix->FirstInCol[I];
while (pElement != NULL) {
if (Reordered) {
Row = pElement->Row;
Col = I;
}
else
{ Row = Matrix->IntToExtRowMap[pElement->Row];
} else {
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
);
if (Err < 0) return 0;
Err = fprintf(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;
if (fprintf(pMatrixFile, "0\t0\t0.0\t0.0\n") < 0)
return 0;
}
#endif /* spCOMPLEX */
#if REAL
if (Data AND NOT Matrix->Complex)
{ for (I = 1; I <= Size; I++)
{ pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{ Row = Matrix->IntToExtRowMap[pElement->Row];
if (Data AND NOT Matrix->Complex) {
for (I = 1; I <= Size; I++) {
pElement = Matrix->FirstInCol[I];
while (pElement != NULL) {
Row = Matrix->IntToExtRowMap[pElement->Row];
Col = Matrix->IntToExtColMap[I];
Err = fprintf
( pMatrixFile,"%d\t%d\t%-.15g\n",
Row, Col, (double)pElement->Real
);
if (Err < 0) return 0;
Err = fprintf(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;
if (fprintf(pMatrixFile, "0\t0\t0.0\n") < 0)
return 0;
}
#endif /* REAL */
/* Close file. */
if (fclose(pMatrixFile) < 0) return 0;
/* Close file. */
if (fclose(pMatrixFile) < 0)
return 0;
return 1;
}
/*!
* Writes vector to file in format suitable to be read back in by the
* matrix test program. This routine should be executed after the function
@ -561,29 +532,25 @@ FILE *pMatrixFile;
* The size of the matrix.
*/
int
spFileVector(
MatrixPtr Matrix,
char *File,
spREAL RHS[]
int spFileVector(MatrixPtr Matrix, char *File, spREAL RHS[]
#if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
, spREAL iRHS[]
,
spREAL iRHS[]
#endif
)
{
int I, Size;
) {
int I, Size;
#if spCOMPLEX
int Err;
int Err;
#endif
FILE *pMatrixFile;
FILE *pMatrixFile;
/* Begin `spFileVector'. */
ASSERT_IS_SPARSE( Matrix );
vASSERT( RHS != NULL, "Vector missing" );
/* Begin `spFileVector'. */
ASSERT_IS_SPARSE(Matrix);
vASSERT(RHS != NULL, "Vector missing");
if (File) {
/* Open File in write mode. */
pMatrixFile = fopen(File,"w");
pMatrixFile = fopen(File, "w");
if (pMatrixFile == NULL)
return 0;
}
@ -591,44 +558,39 @@ FILE *pMatrixFile;
/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spCOMPLEX
if (Matrix->Complex)
{
if (Matrix->Complex) {
#if spSEPARATED_COMPLEX_VECTORS
vASSERT( iRHS != NULL, "Imaginary vector missing" );
vASSERT(iRHS != NULL, "Imaginary vector missing");
--RHS;
--iRHS;
#else
RHS -= 2;
#endif
}
else
} else
#endif /* spCOMPLEX */
--RHS;
#endif /* NOT ARRAY_OFFSET */
/* Output vector. */
/* Output vector. */
Size = Matrix->Size;
if (Size == 0) return 1;
if (Size == 0)
return 1;
#if spCOMPLEX
if (Matrix->Complex)
{
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]
);
if (Err < 0) return 0;
for (I = 1; I <= Size; I++) {
Err = fprintf(pMatrixFile, "%-.15g\t%-.15g\n", (double)RHS[I],
(double)iRHS[I]);
if (Err < 0)
return 0;
}
#else
for (I = 1; I <= Size; I++)
{ Err = fprintf
( pMatrixFile, "%-.15g\t%-.15g\n",
(double)RHS[2*I], (double)RHS[2*I+1]
);
if (Err < 0) return 0;
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
}
@ -637,26 +599,20 @@ FILE *pMatrixFile;
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. */
if (fclose(pMatrixFile) < 0) return 0;
/* Close file. */
if (fclose(pMatrixFile) < 0)
return 0;
return 1;
}
/*!
* Writes useful information concerning the matrix to a file. Should be
* executed after the matrix is factored.
@ -690,27 +646,21 @@ FILE *pMatrixFile;
* The smallest element in the matrix excluding zero elements.
*/
int
spFileStats(
MatrixPtr Matrix,
char *File,
char *Label
)
{
int Size, I;
ElementPtr pElement;
int NumberOfElements;
RealNumber Data, LargestElement, SmallestElement;
FILE *pStatsFile;
int spFileStats(MatrixPtr Matrix, char *File, char *Label) {
int Size, I;
ElementPtr pElement;
int NumberOfElements;
RealNumber Data, LargestElement, SmallestElement;
FILE *pStatsFile;
/* Begin `spFileStats'. */
ASSERT_IS_SPARSE( Matrix );
/* Begin `spFileStats'. */
ASSERT_IS_SPARSE(Matrix);
/* Open File in append mode. */
/* Open File in append mode. */
if ((pStatsFile = fopen(File, "a")) == NULL)
return 0;
/* Output statistics. */
/* Output statistics. */
Size = Matrix->Size;
if (NOT Matrix->Factored)
fprintf(pStatsFile, "Matrix has not been factored.\n");
@ -720,18 +670,19 @@ FILE *pStatsFile;
fprintf(pStatsFile, "Matrix is complex.\n");
else
fprintf(pStatsFile, "Matrix is real.\n");
fprintf(pStatsFile," Size = %d\n",Size);
if (Size == 0) return 1;
fprintf(pStatsFile, " Size = %d\n", Size);
if (Size == 0)
return 1;
/* Search matrix. */
/* Search matrix. */
NumberOfElements = 0;
LargestElement = 0.0;
SmallestElement = LARGEST_REAL;
for (I = 1; I <= Size; I++)
{ pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{ NumberOfElements++;
for (I = 1; I <= Size; I++) {
pElement = Matrix->FirstInCol[I];
while (pElement != NULL) {
NumberOfElements++;
Data = ELEMENT_MAG(pElement);
if (Data > LargestElement)
LargestElement = Data;
@ -741,29 +692,29 @@ FILE *pStatsFile;
}
}
SmallestElement = MIN( SmallestElement, LargestElement );
SmallestElement = MIN(SmallestElement, LargestElement);
/* Output remaining statistics. */
/* Output remaining statistics. */
fprintf(pStatsFile, " Initial number of elements = %d\n",
NumberOfElements - Matrix->Fillins);
fprintf(pStatsFile,
" Initial average number of elements per row = %f\n",
(double)(NumberOfElements - Matrix->Fillins) / (double)Size);
fprintf(pStatsFile, " Fill-ins = %d\n",Matrix->Fillins);
fprintf(pStatsFile, " Fill-ins = %d\n", Matrix->Fillins);
fprintf(pStatsFile, " Average number of fill-ins per row = %f%%\n",
(double)Matrix->Fillins / (double)Size);
fprintf(pStatsFile, " Total number of elements = %d\n",
NumberOfElements);
fprintf(pStatsFile, " Average number of elements per row = %f\n",
(double)NumberOfElements / (double)Size);
fprintf(pStatsFile," Density = %f%%\n",
(100.0*(double)NumberOfElements)/((double)Size*(double)Size));
fprintf(pStatsFile," Relative Threshold = %e\n", Matrix->RelThreshold);
fprintf(pStatsFile," Absolute Threshold = %e\n", Matrix->AbsThreshold);
fprintf(pStatsFile," Largest Element = %e\n", LargestElement);
fprintf(pStatsFile," Smallest Element = %e\n\n\n", SmallestElement);
fprintf(pStatsFile, " Density = %f%%\n",
(100.0 * (double)NumberOfElements) / ((double)Size * (double)Size));
fprintf(pStatsFile, " Relative Threshold = %e\n", Matrix->RelThreshold);
fprintf(pStatsFile, " Absolute Threshold = %e\n", Matrix->AbsThreshold);
fprintf(pStatsFile, " Largest Element = %e\n", LargestElement);
fprintf(pStatsFile, " Smallest Element = %e\n\n\n", SmallestElement);
/* Close file. */
/* Close file. */
(void)fclose(pStatsFile);
return 1;
}

331
src/maths/sparse/spSMP.c
View File

@ -34,7 +34,7 @@
* SMPcProdDiag
* LoadGmin
* SMPfindElt
*/
*/
/*
* To replace SMP with Sparse, rename the file spSpice3.h to
@ -80,17 +80,14 @@
*
* Copyright (c) 1985-2003 by Kenneth S. Kundert
*/
#ifdef notdef
static char copyright[] =
"Sparse1.4: Copyright (c) 1985-2003 by Kenneth S. Kundert";
static char RCSid[] =
"@(#)$Header: /cvsroot/sparse/src/spSMP.c,v 1.3 2003/06/30 19:40:51 kundert Exp $";
static char RCSid[] = "@(#)$Header: /cvsroot/sparse/src/spSMP.c,v 1.3 "
"2003/06/30 19:40:51 kundert Exp $";
#endif
/*
* IMPORTS
*
@ -101,229 +98,164 @@ static char RCSid[] =
* Spice3's matrix macro definitions.
*/
#include "ngspice/spmatrix.h"
#include "ngspice/smpdefs.h"
#include "ngspice/spmatrix.h"
#include "spDefs.h"
#define NO 0
#define YES 1
#define NO 0
#define YES 1
#define NG_IGNORE(x) (void)x
#define NG_IGNORE(x) (void)x
static void LoadGmin(MatrixPtr Matrix, double Gmin);
/*
* SMPaddElt()
*/
int
SMPaddElt(
SMPmatrix *Matrix,
int Row, int Col,
double Value)
{
*spGetElement( Matrix, Row, Col ) = Value;
return spError( Matrix );
int SMPaddElt(SMPmatrix *Matrix, int Row, int Col, double Value) {
*spGetElement(Matrix, Row, Col) = Value;
return spError(Matrix);
}
/*
* SMPmakeElt()
*/
double *
SMPmakeElt(
SMPmatrix *Matrix,
int Row, int Col)
{
return spGetElement( Matrix, Row, Col );
double *SMPmakeElt(SMPmatrix *Matrix, int Row, int Col) {
return spGetElement(Matrix, Row, Col);
}
/*
* SMPcClear()
*/
void
SMPcClear(
SMPmatrix *Matrix)
{
spClear( Matrix );
}
void SMPcClear(SMPmatrix *Matrix) { spClear(Matrix); }
/*
* SMPclear()
*/
void
SMPclear(
SMPmatrix *Matrix)
{
spClear( Matrix );
}
void SMPclear(SMPmatrix *Matrix) { spClear(Matrix); }
/*
* SMPcLUfac()
*/
/*ARGSUSED*/
int
SMPcLUfac(
SMPmatrix *Matrix,
double PivTol)
{
int SMPcLUfac(SMPmatrix *Matrix, double PivTol) {
NG_IGNORE(PivTol);
spSetComplex( Matrix );
return spFactor( Matrix );
spSetComplex(Matrix);
return spFactor(Matrix);
}
/*
* SMPluFac()
*/
/*ARGSUSED*/
int
SMPluFac(
SMPmatrix *Matrix,
double PivTol, double Gmin)
{
int SMPluFac(SMPmatrix *Matrix, double PivTol, double Gmin) {
NG_IGNORE(PivTol);
spSetReal( Matrix );
LoadGmin( Matrix, Gmin );
return spFactor( Matrix );
spSetReal(Matrix);
LoadGmin(Matrix, Gmin);
return spFactor(Matrix);
}
/*
* SMPcReorder()
*/
int
SMPcReorder(
SMPmatrix *Matrix,
double PivTol, double PivRel,
int *NumSwaps)
{
int SMPcReorder(SMPmatrix *Matrix, double PivTol, double PivRel,
int *NumSwaps) {
*NumSwaps = 1;
spSetComplex( Matrix );
return spOrderAndFactor( Matrix, NULL,
PivRel, PivTol, YES );
spSetComplex(Matrix);
return spOrderAndFactor(Matrix, NULL, PivRel, PivTol, YES);
}
/*
* SMPreorder()
*/
int
SMPreorder(
SMPmatrix *Matrix,
double PivTol, double PivRel, double Gmin)
{
spSetReal( Matrix );
LoadGmin( Matrix, Gmin );
return spOrderAndFactor( Matrix, NULL,
PivRel, PivTol, YES );
int SMPreorder(SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin) {
spSetReal(Matrix);
LoadGmin(Matrix, Gmin);
return spOrderAndFactor(Matrix, NULL, PivRel, PivTol, YES);
}
/*
* SMPcaSolve()
*/
void
SMPcaSolve(
SMPmatrix *Matrix,
double RHS[], double iRHS[], double Spare[], double iSpare[])
{
void SMPcaSolve(SMPmatrix *Matrix, double RHS[], double iRHS[], double Spare[],
double iSpare[]) {
NG_IGNORE(Spare);
NG_IGNORE(iSpare);
#if spCOMPLEX
spSolveTransposed( Matrix, RHS, RHS, iRHS, iRHS );
spSolveTransposed(Matrix, RHS, RHS, iRHS, iRHS);
#else
spSolveTransposed( Matrix, RHS, RHS );
spSolveTransposed(Matrix, RHS, RHS);
#endif
}
/*
* SMPcSolve()
*/
void
SMPcSolve(
SMPmatrix *Matrix,
double RHS[], double iRHS[], double Spare[], double iSpare[])
{
void SMPcSolve(SMPmatrix *Matrix, double RHS[], double iRHS[], double Spare[],
double iSpare[]) {
NG_IGNORE(Spare);
NG_IGNORE(iSpare);
#if spCOMPLEX
spSolve( Matrix, RHS, RHS, iRHS, iRHS );
spSolve(Matrix, RHS, RHS, iRHS, iRHS);
#else
spSolve( Matrix, RHS, RHS );
spSolve(Matrix, RHS, RHS);
#endif
}
/*
* SMPsolve()
*/
void
SMPsolve(
SMPmatrix *Matrix,
double RHS[], double Spare[])
{
void SMPsolve(SMPmatrix *Matrix, double RHS[], double Spare[]) {
NG_IGNORE(Spare);
#if spCOMPLEX
spSolve( Matrix, RHS, RHS, NULL, NULL );
spSolve(Matrix, RHS, RHS, NULL, NULL);
#else
spSolve( Matrix, RHS, RHS );
spSolve(Matrix, RHS, RHS);
#endif
}
/*
* SMPmatSize()
*/
int
SMPmatSize(
SMPmatrix *Matrix)
{
return spGetSize( Matrix, 1 );
}
int SMPmatSize(SMPmatrix *Matrix) { return spGetSize(Matrix, 1); }
/*
* SMPnewMatrix()
*/
int
SMPnewMatrix(
SMPmatrix **pMatrix,
int size)
{
int SMPnewMatrix(SMPmatrix **pMatrix, int size) {
int Error;
*pMatrix = (SMPmatrix *)spCreate( size, 1, &Error );
*pMatrix = (SMPmatrix *)spCreate(size, 1, &Error);
return Error;
}
/*
* SMPdestroy()
*/
void
SMPdestroy(
SMPmatrix *Matrix)
{
spDestroy( Matrix );
}
void SMPdestroy(SMPmatrix *Matrix) { spDestroy(Matrix); }
/*
* SMPpreOrder()
*/
int
SMPpreOrder(
SMPmatrix *Matrix)
{
spMNA_Preorder( Matrix );
return spError( Matrix );
int SMPpreOrder(SMPmatrix *Matrix) {
spMNA_Preorder(Matrix);
return spError(Matrix);
}
/*
* SMPprintRHS()
*/
void
SMPprintRHS(SMPmatrix *Matrix, char *Filename, RealVector RHS, RealVector iRHS)
{
void SMPprintRHS(SMPmatrix *Matrix, char *Filename, RealVector RHS,
RealVector iRHS) {
#if spCOMPLEX
spFileVector( Matrix, Filename, RHS, iRHS );
spFileVector(Matrix, Filename, RHS, iRHS);
#else
spFileVector( Matrix, Filename, RHS );
spFileVector(Matrix, Filename, RHS);
#endif
}
@ -331,66 +263,50 @@ SMPprintRHS(SMPmatrix *Matrix, char *Filename, RealVector RHS, RealVector iRHS)
* SMPprint()
*/
/*ARGSUSED*/
void
SMPprint(
SMPmatrix *Matrix,
char *Filename)
{
void SMPprint(SMPmatrix *Matrix, char *Filename) {
if (Filename)
spFileMatrix(Matrix, Filename, "Circuit Matrix", 0, 1, 1 );
spFileMatrix(Matrix, Filename, "Circuit Matrix", 0, 1, 1);
else
spPrint( Matrix, 0, 1, 1 );
spPrint(Matrix, 0, 1, 1);
}
/*
* SMPgetError()
*/
void
SMPgetError(
SMPmatrix *Matrix,
int *Row, int *Col)
{
spWhereSingular( Matrix, Row, Col );
void SMPgetError(SMPmatrix *Matrix, int *Row, int *Col) {
spWhereSingular(Matrix, Row, Col);
}
/*
* SMPcProdDiag()
*/
int
SMPcProdDiag(
SMPmatrix *Matrix,
SPcomplex *pMantissa,
int *pExponent)
{
int SMPcProdDiag(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent) {
#if spCOMPLEX
spDeterminant( Matrix, pExponent, &(pMantissa->real),
&(pMantissa->imag) );
spDeterminant(Matrix, pExponent, &(pMantissa->real), &(pMantissa->imag));
#else
spDeterminant( Matrix, pExponent, &(pMantissa->real) );
spDeterminant(Matrix, pExponent, &(pMantissa->real));
#endif
return spError( Matrix );
return spError(Matrix);
}
/*
* SMPcDProd()
*/
int
SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
{
double re, im, x, y, z;
int p;
int SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent) {
double re, im, x, y, z;
int p;
#if spCOMPLEX
spDeterminant( Matrix, &p, &re, &im);
spDeterminant(Matrix, &p, &re, &im);
#else
spDeterminant( Matrix, &p, &re );
spDeterminant(Matrix, &p, &re);
im = 0.0;
#endif
#ifndef M_LN2
#define M_LN2 0.69314718055994530942
#define M_LN2 0.69314718055994530942
#endif
#ifndef M_LN10
#define M_LN10 2.30258509299404568402
#define M_LN10 2.30258509299404568402
#endif
#ifdef debug_print
@ -399,7 +315,7 @@ SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
/* Convert base 10 numbers to base 2 numbers, for comparison */
y = p * M_LN10 / M_LN2;
x = (int) y;
x = (int)y;
y -= x;
/* ASSERT
@ -419,45 +335,44 @@ SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
/* Re-normalize (re or im may be > 2.0 or both < 1.0 */
if (re != 0.0) {
y = logb(re);
if (im != 0.0)
z = logb(im);
else
z = 0;
y = logb(re);
if (im != 0.0)
z = logb(im);
else
z = 0;
} else if (im != 0.0) {
z = logb(im);
y = 0;
z = logb(im);
y = 0;
} else {
/* Singular */
/*printf("10 -> singular\n");*/
y = 0;
z = 0;
/* Singular */
/*printf("10 -> singular\n");*/
y = 0;
z = 0;
}
#ifdef debug_print
printf(" ** renormalize changes = %g,%g\n", y, z);
#endif
if (y < z)
y = z;
y = z;
#ifdef debug_print
*pExponent = (int)(x + y);
x = scalbn(re, (int) -y);
z = scalbn(im, (int) -y);
printf(" ** values are: re %g, im %g, y %g, re' %g, im' %g\n",
re, im, y, x, z);
x = scalbn(re, (int)-y);
z = scalbn(im, (int)-y);
printf(" ** values are: re %g, im %g, y %g, re' %g, im' %g\n", re, im, y, x,
z);
#endif
pMantissa->real = scalbn(re, (int) -y);
pMantissa->imag = scalbn(im, (int) -y);
pMantissa->real = scalbn(re, (int)-y);
pMantissa->imag = scalbn(im, (int)-y);
#ifdef debug_print
printf("Determinant 10->2: (%20g,%20g)^%d\n", pMantissa->real,
pMantissa->imag, *pExponent);
pMantissa->imag, *pExponent);
#endif
return spError( Matrix );
return spError(Matrix);
}
/*
* LOAD GMIN
*
@ -468,17 +383,13 @@ SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
* use of this routine is not recommended. It is included here simply
* for compatibility with Spice3.
*/
static void
LoadGmin(
MatrixPtr Matrix,
double Gmin)
{
int I;
ArrayOfElementPtrs Diag;
ElementPtr diag;
static void LoadGmin(MatrixPtr Matrix, double Gmin) {
int I;
ArrayOfElementPtrs Diag;
ElementPtr diag;
/* Begin `spLoadGmin'. */
ASSERT_IS_SPARSE( Matrix );
/* Begin `spLoadGmin'. */
ASSERT_IS_SPARSE(Matrix);
if (Gmin != 0.0) {
Diag = Matrix->Diag;
@ -490,9 +401,6 @@ ElementPtr diag;
return;
}
/*
* FIND ELEMENT
*
@ -502,18 +410,17 @@ ElementPtr diag;
* pointer to the new element is returned.
*/
SMPelement *
SMPfindElt(SMPmatrix *Matrix, int Row, int Col, int CreateIfMissing)
{
SMPelement *SMPfindElt(SMPmatrix *Matrix, int Row, int Col,
int CreateIfMissing) {
ElementPtr Element;
/* Begin `SMPfindElt'. */
ASSERT_IS_SPARSE( Matrix );
ASSERT_IS_SPARSE(Matrix);
Row = Matrix->ExtToIntRowMap[Row];
Col = Matrix->ExtToIntColMap[Col];
if (Col == -1)
/* No element available */
/* No element available */
return NULL;
Element = Matrix->FirstInCol[Col];
@ -521,37 +428,30 @@ SMPfindElt(SMPmatrix *Matrix, int Row, int Col, int CreateIfMissing)
return Element;
}
/*
* SMPcZeroCol()
*/
int
SMPcZeroCol(SMPmatrix *Matrix, int Col)
{
ElementPtr Element;
int SMPcZeroCol(SMPmatrix *Matrix, int Col) {
ElementPtr Element;
Col = Matrix->ExtToIntColMap[Col];
for (Element = Matrix->FirstInCol[Col];
Element != NULL;
Element = Element->NextInCol)
{
Element->Real = 0.0;
for (Element = Matrix->FirstInCol[Col]; Element != NULL;
Element = Element->NextInCol) {
Element->Real = 0.0;
#if spCOMPLEX
Element->Imag = 0.0;
Element->Imag = 0.0;
#endif
}
return spError( Matrix );
return spError(Matrix);
}
/*
* SMPcAddCol()
*/
int
SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col)
{
ElementPtr Accum, Addend, *Prev;
int SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col) {
ElementPtr Accum, Addend, *Prev;
Accum_Col = Matrix->ExtToIntColMap[Accum_Col];
Addend_Col = Matrix->ExtToIntColMap[Addend_Col];
@ -566,7 +466,8 @@ SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col)
Accum = *Prev;
}
if (!Accum || Accum->Row > Addend->Row) {
Accum = spcCreateElement(Matrix, Addend->Row, Accum_Col, Prev, 0, NO);
Accum =
spcCreateElement(Matrix, Addend->Row, Accum_Col, Prev, 0, NO);
}
Accum->Real += Addend->Real;
#if spCOMPLEX
@ -575,28 +476,24 @@ SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col)
Addend = Addend->NextInCol;
}
return spError( Matrix );
return spError(Matrix);
}
/*
* SMPconstMult()
*/
void
SMPconstMult(SMPmatrix *Matrix, double constant)
{
void SMPconstMult(SMPmatrix *Matrix, double constant) {
spConstMult(Matrix, constant);
}
/*
* SMPmultiply()
*/
void
SMPmultiply(SMPmatrix *Matrix, double *RHS, double *Solution, double *iRHS, double *iSolution)
{
void SMPmultiply(SMPmatrix *Matrix, double *RHS, double *Solution, double *iRHS,
double *iSolution) {
#if spCOMPLEX
spMultiply(Matrix, RHS, Solution, iRHS, iSolution);
#else
spMultiply(Matrix, RHS, Solution);
#endif
}

390
src/maths/sparse/spSolve.c
View File

@ -24,7 +24,6 @@
* SolveComplexTransposedMatrix
*/
/*
* Revision and copyright information.
*
@ -35,12 +34,10 @@
#ifdef notdef
static char copyright[] =
"Sparse1.4: Copyright (c) 1985-2003 by Kenneth S. Kundert";
static char RCSid[] =
"@(#)$Header: /cvsroot/sparse/src/spSolve.c,v 1.3 2003/06/29 04:19:52 kundert Exp $";
static char RCSid[] = "@(#)$Header: /cvsroot/sparse/src/spSolve.c,v 1.3 "
"2003/06/29 04:19:52 kundert Exp $";
#endif
/*
* IMPORTS
*
@ -54,13 +51,10 @@ static char RCSid[] =
*/
#define spINSIDE_SPARSE
#include <stdio.h>
#include "spConfig.h"
#include "ngspice/spmatrix.h"
#include "spConfig.h"
#include "spDefs.h"
#include <stdio.h>
/*
* Function declarations
@ -68,23 +62,16 @@ static char RCSid[] =
#if spSEPARATED_COMPLEX_VECTORS
#if spCOMPLEX
static void SolveComplexMatrix( MatrixPtr,
RealVector, RealVector, RealVector, RealVector );
static void SolveComplexTransposedMatrix( MatrixPtr,
RealVector, RealVector, RealVector, RealVector );
static void SolveComplexMatrix(MatrixPtr, RealVector, RealVector, RealVector,
RealVector);
static void SolveComplexTransposedMatrix(MatrixPtr, RealVector, RealVector,
RealVector, RealVector);
#endif
#else
static void SolveComplexMatrix( MatrixPtr, RealVector, RealVector );
static void SolveComplexTransposedMatrix( MatrixPtr,
RealVector, RealVector );
static void SolveComplexMatrix(MatrixPtr, RealVector, RealVector);
static void SolveComplexTransposedMatrix(MatrixPtr, RealVector, RealVector);
#endif
/*!
* Performs forward elimination and back substitution to find the
* unknown vector from the \a RHS vector and factored matrix. This
@ -144,33 +131,28 @@ static void SolveComplexTransposedMatrix( MatrixPtr,
/*VARARGS3*/
void
spSolve(
MatrixPtr Matrix,
spREAL RHS[],
spREAL Solution[]
# if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
, spREAL iRHS[]
, spREAL iSolution[]
# endif
)
{
void spSolve(MatrixPtr Matrix, spREAL RHS[], spREAL Solution[]
#if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
,
spREAL iRHS[], spREAL iSolution[]
#endif
) {
#if REAL
ElementPtr pElement;
RealVector Intermediate;
RealNumber Temp;
int I, *pExtOrder, Size;
ElementPtr pPivot;
ElementPtr pElement;
RealVector Intermediate;
RealNumber Temp;
int I, *pExtOrder, Size;
ElementPtr pPivot;
#endif
/* Begin `spSolve'. */
ASSERT_IS_SPARSE( Matrix );
ASSERT_NO_ERRORS( Matrix );
ASSERT_IS_FACTORED( Matrix );
/* Begin `spSolve'. */
ASSERT_IS_SPARSE(Matrix);
ASSERT_NO_ERRORS(Matrix);
ASSERT_IS_FACTORED(Matrix);
#if spCOMPLEX
if (Matrix->Complex)
{ SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS );
if (Matrix->Complex) {
SolveComplexMatrix(Matrix, RHS, Solution IMAG_VECTORS);
return;
}
#endif
@ -185,39 +167,39 @@ ElementPtr pPivot;
--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.*/
for (I = 1; I <= Size; I++)
{
/* This step of the elimination is skipped if Temp equals zero. */
if ((Temp = Intermediate[I]) != 0.0)
{ pPivot = Matrix->Diag[I];
/* Forward elimination. Solves Lc = b.*/
for (I = 1; I <= Size; I++) {
/* This step of the elimination is skipped if Temp equals zero. */
if ((Temp = Intermediate[I]) != 0.0) {
pPivot = Matrix->Diag[I];
Intermediate[I] = (Temp *= pPivot->Real);
pElement = pPivot->NextInCol;
while (pElement != NULL)
{ Intermediate[pElement->Row] -= Temp * pElement->Real;
while (pElement != NULL) {
Intermediate[pElement->Row] -= Temp * pElement->Real;
pElement = pElement->NextInCol;
}
}
}
/* Backward Substitution. Solves Ux = c.*/
for (I = Size; I > 0; I--)
{ Temp = Intermediate[I];
/* Backward Substitution. Solves Ux = c.*/
for (I = Size; I > 0; I--) {
Temp = Intermediate[I];
pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{ Temp -= pElement->Real * Intermediate[pElement->Col];
while (pElement != NULL) {
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];
@ -226,16 +208,6 @@ ElementPtr pPivot;
#endif /* REAL */
}
#if spCOMPLEX
/*!
* Performs forward elimination and back substitution to find the
@ -289,27 +261,23 @@ ElementPtr pPivot;
* Temporary storage for entries in arrays.
*/
static void
SolveComplexMatrix(
MatrixPtr Matrix,
RealVector RHS,
RealVector Solution
# if spSEPARATED_COMPLEX_VECTORS
, RealVector iRHS
, RealVector iSolution
# endif
)
{
ElementPtr pElement;
ComplexVector Intermediate;
int I, *pExtOrder, Size;
ElementPtr pPivot;
#if NOT spSEPARATED_COMPLEX_VECTORS
ComplexVector ExtVector;
static void SolveComplexMatrix(MatrixPtr Matrix, RealVector RHS,
RealVector Solution
#if spSEPARATED_COMPLEX_VECTORS
,
RealVector iRHS, RealVector iSolution
#endif
ComplexNumber Temp;
) {
ElementPtr pElement;
ComplexVector Intermediate;
int I, *pExtOrder, Size;
ElementPtr pPivot;
#if NOT spSEPARATED_COMPLEX_VECTORS
ComplexVector ExtVector;
#endif
ComplexNumber Temp;
/* Begin `SolveComplexMatrix'. */
/* Begin `SolveComplexMatrix'. */
Size = Matrix->Size;
Intermediate = (ComplexVector)Matrix->Intermediate;
@ -317,19 +285,22 @@ ComplexNumber Temp;
/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spSEPARATED_COMPLEX_VECTORS
--RHS; --iRHS;
--Solution; --iSolution;
--RHS;
--iRHS;
--Solution;
--iSolution;
#else
RHS -= 2; Solution -= 2;
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)];
for (I = Size; I > 0; I--) {
Intermediate[I].Real = RHS[*(pExtOrder)];
Intermediate[I].Imag = iRHS[*(pExtOrder--)];
}
#else
@ -338,47 +309,47 @@ ComplexNumber Temp;
Intermediate[I] = ExtVector[*(pExtOrder--)];
#endif
/* Forward substitution. Solves Lc = b.*/
for (I = 1; I <= Size; I++)
{ Temp = Intermediate[I];
/* Forward substitution. Solves Lc = b.*/
for (I = 1; I <= Size; 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) OR(Temp.Imag != 0.0)) {
pPivot = Matrix->Diag[I];
/* Cmplx expr: Temp *= (1.0 / Pivot). */
CMPLX_MULT_ASSIGN(Temp, *pPivot);
Intermediate[I] = Temp;
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
/* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */
CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Row],
Temp, *pElement);
while (pElement != NULL) {
/* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */
CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Row], Temp,
*pElement);
pElement = pElement->NextInCol;
}
}
}
/* Backward Substitution. Solves Ux = c.*/
for (I = Size; I > 0; I--)
{ Temp = Intermediate[I];
/* Backward Substitution. Solves Ux = c.*/
for (I = Size; I > 0; I--) {
Temp = Intermediate[I];
pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{
/* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */
CMPLX_MULT_SUBT_ASSIGN(Temp, *pElement,Intermediate[pElement->Col]);
while (pElement != NULL) {
/* 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;
for (I = Size; I > 0; I--) {
Solution[*(pExtOrder)] = Intermediate[I].Real;
iSolution[*(pExtOrder--)] = Intermediate[I].Imag;
}
#else
@ -391,19 +362,6 @@ ComplexNumber Temp;
}
#endif /* spCOMPLEX */
#if TRANSPOSE
/*!
* Performs forward elimination and back substitution to find the
@ -458,33 +416,28 @@ ComplexNumber Temp;
/*VARARGS3*/
void
spSolveTransposed(
MatrixPtr Matrix,
spREAL RHS[],
spREAL Solution[]
# if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
, spREAL iRHS[]
, spREAL iSolution[]
# endif
)
{
void spSolveTransposed(MatrixPtr Matrix, spREAL RHS[], spREAL Solution[]
#if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
,
spREAL iRHS[], spREAL iSolution[]
#endif
) {
#if REAL
ElementPtr pElement;
RealVector Intermediate;
int I, *pExtOrder, Size;
ElementPtr pPivot;
RealNumber Temp;
ElementPtr pElement;
RealVector Intermediate;
int I, *pExtOrder, Size;
ElementPtr pPivot;
RealNumber Temp;
#endif
/* Begin `spSolveTransposed'. */
ASSERT_IS_SPARSE( Matrix );
ASSERT_NO_ERRORS( Matrix );
ASSERT_IS_FACTORED( Matrix );
/* Begin `spSolveTransposed'. */
ASSERT_IS_SPARSE(Matrix);
ASSERT_NO_ERRORS(Matrix);
ASSERT_IS_FACTORED(Matrix);
#if spCOMPLEX
if (Matrix->Complex)
{ SolveComplexTransposedMatrix( Matrix, RHS, Solution IMAG_VECTORS );
if (Matrix->Complex) {
SolveComplexTransposedMatrix(Matrix, RHS, Solution IMAG_VECTORS);
return;
}
#endif
@ -499,38 +452,37 @@ RealNumber Temp;
--Solution;
#endif
/* Initialize Intermediate vector. */
/* Initialize Intermediate vector. */
pExtOrder = &Matrix->IntToExtColMap[Size];
for (I = Size; I > 0; I--)
Intermediate[I] = RHS[*(pExtOrder--)];
/* Forward elimination. */
for (I = 1; I <= Size; I++)
{
/* This step of the elimination is skipped if Temp equals zero. */
if ((Temp = Intermediate[I]) != 0.0)
{ pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{ Intermediate[pElement->Col] -= Temp * pElement->Real;
/* Forward elimination. */
for (I = 1; I <= Size; I++) {
/* This step of the elimination is skipped if Temp equals zero. */
if ((Temp = Intermediate[I]) != 0.0) {
pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL) {
Intermediate[pElement->Col] -= Temp * pElement->Real;
pElement = pElement->NextInRow;
}
}
}
/* Backward Substitution. */
for (I = Size; I > 0; I--)
{ pPivot = Matrix->Diag[I];
/* Backward Substitution. */
for (I = Size; I > 0; I--) {
pPivot = Matrix->Diag[I];
Temp = Intermediate[I];
pElement = pPivot->NextInCol;
while (pElement != NULL)
{ Temp -= pElement->Real * Intermediate[pElement->Row];
while (pElement != NULL) {
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];
@ -540,15 +492,6 @@ RealNumber Temp;
}
#endif /* TRANSPOSE */
#if TRANSPOSE AND spCOMPLEX
/*!
* Performs forward elimination and back substitution to find the
@ -605,27 +548,23 @@ RealNumber Temp;
* Temporary storage for entries in arrays.
*/
static void
SolveComplexTransposedMatrix(
MatrixPtr Matrix,
RealVector RHS,
RealVector Solution
# if spSEPARATED_COMPLEX_VECTORS
, RealVector iRHS
, RealVector iSolution
# endif
)
{
ElementPtr pElement;
ComplexVector Intermediate;
int I, *pExtOrder, Size;
#if NOT spSEPARATED_COMPLEX_VECTORS
ComplexVector ExtVector;
static void SolveComplexTransposedMatrix(MatrixPtr Matrix, RealVector RHS,
RealVector Solution
#if spSEPARATED_COMPLEX_VECTORS
,
RealVector iRHS, RealVector iSolution
#endif
ElementPtr pPivot;
ComplexNumber Temp;
) {
ElementPtr pElement;
ComplexVector Intermediate;
int I, *pExtOrder, Size;
#if NOT spSEPARATED_COMPLEX_VECTORS
ComplexVector ExtVector;
#endif
ElementPtr pPivot;
ComplexNumber Temp;
/* Begin `SolveComplexTransposedMatrix'. */
/* Begin `SolveComplexTransposedMatrix'. */
Size = Matrix->Size;
Intermediate = (ComplexVector)Matrix->Intermediate;
@ -633,19 +572,22 @@ ComplexNumber Temp;
/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spSEPARATED_COMPLEX_VECTORS
--RHS; --iRHS;
--Solution; --iSolution;
--RHS;
--iRHS;
--Solution;
--iSolution;
#else
RHS -= 2; Solution -= 2;
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)];
for (I = Size; I > 0; I--) {
Intermediate[I].Real = RHS[*(pExtOrder)];
Intermediate[I].Imag = iRHS[*(pExtOrder--)];
}
#else
@ -654,46 +596,46 @@ ComplexNumber Temp;
Intermediate[I] = ExtVector[*(pExtOrder--)];
#endif
/* Forward elimination. */
for (I = 1; I <= Size; I++)
{ Temp = Intermediate[I];
/* Forward elimination. */
for (I = 1; I <= Size; 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;
while (pElement != NULL)
{
/* Cmplx expr: Intermediate[Element->Col] -= Temp * *Element. */
CMPLX_MULT_SUBT_ASSIGN( Intermediate[pElement->Col],
Temp, *pElement);
/* This step of the elimination is skipped if Temp equals zero. */
if ((Temp.Real != 0.0) OR(Temp.Imag != 0.0)) {
pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL) {
/* Cmplx expr: Intermediate[Element->Col] -= Temp * *Element. */
CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Col], Temp,
*pElement);
pElement = pElement->NextInRow;
}
}
}
/* Backward Substitution. */
for (I = Size; I > 0; I--)
{ pPivot = Matrix->Diag[I];
/* Backward Substitution. */
for (I = Size; I > 0; I--) {
pPivot = Matrix->Diag[I];
Temp = Intermediate[I];
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
/* Cmplx expr: Temp -= Intermediate[Element->Row] * *Element. */
CMPLX_MULT_SUBT_ASSIGN(Temp,Intermediate[pElement->Row],*pElement);
while (pElement != NULL) {
/* 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;
for (I = Size; I > 0; I--) {
Solution[*(pExtOrder)] = Intermediate[I].Real;
iSolution[*(pExtOrder--)] = Intermediate[I].Imag;
}
#else

1544
src/maths/sparse/spUtils.c
View File

File diff suppressed because it is too large Load Diff