epilectrik
/
pyNgSpice
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

KLU Integration from scratch #1, new files

pre-master-46
Francesco Lannutti 12 years ago
committed by Holger Vogt
parent
60b4c3cda7
commit
cf5991ae81
4 changed files with 1010 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 20
      src/include/ngspice/klu-binding.h
  2. 97
      src/maths/KLU/Makefile.am
  3. 748
      src/maths/KLU/klusmp.c
  4. 145
      src/maths/sparse/spCSC.c

20
src/include/ngspice/klu-binding.h
View File

@ -0,0 +1,20 @@
#ifndef _KLU_BINDING_H
#define _KLU_BINDING_H
#define CREATE_KLU_BINDING_TABLE(ptr, binding, a, b) \
if ((here->a != 0) && (here->b != 0)) { \
i = here->ptr ; \
matched = (BindElement *) bsearch (&i, BindStruct, nz, sizeof(BindElement), BindCompare) ; \
here->binding = matched ; \
here->ptr = matched->CSC ; \
}
#define CONVERT_KLU_BINDING_TABLE_TO_COMPLEX(ptr, binding, a, b) \
if ((here->a != 0) && (here->b != 0)) \
here->ptr = here->binding->CSC_Complex ;
#define CONVERT_KLU_BINDING_TABLE_TO_REAL(ptr, binding, a, b) \
if ((here->a != 0) && (here->b != 0)) \
here->ptr = here->binding->CSC ;
#endif

97
src/maths/KLU/Makefile.am
View File

@ -0,0 +1,97 @@
## Process this file with automake to produce Makefile.in
noinst_LTLIBRARIES = libKLU_real.la libKLU_complex.la libKLU.la
libKLU_real_la_SOURCES = \
amd_1.c \
amd_2.c \
amd_aat.c \
amd_control.c \
amd_defaults.c \
amd_dump.c \
amd_global.c \
amd_info.c \
amd_order.c \
amd_postorder.c \
amd_post_tree.c \
amd_preprocess.c \
amd_valid.c \
btf_maxtrans.c \
btf_order.c \
btf_strongcomp.c \
colamd.c \
colamd_global.c \
klu.c \
klu_analyze.c \
klu_analyze_given.c \
klu_defaults.c \
klu_diagnostics.c \
klu_dump.c \
klu_extract.c \
klu_factor.c \
klu_free_numeric.c \
klu_free_symbolic.c \
klu_kernel.c \
klu_memory.c \
klu_refactor.c \
klu_scale.c \
klu_solve.c \
klu_sort.c \
klu_tsolve.c
libKLU_real_la_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/src/include
libKLU_complex_la_SOURCES = \
amd_1.c \
amd_2.c \
amd_aat.c \
amd_control.c \
amd_defaults.c \
amd_dump.c \
amd_global.c \
amd_info.c \
amd_order.c \
amd_postorder.c \
amd_post_tree.c \
amd_preprocess.c \
amd_valid.c \
btf_maxtrans.c \
btf_order.c \
btf_strongcomp.c \
colamd.c \
colamd_global.c \
klu.c \
klu_analyze.c \
klu_analyze_given.c \
klu_defaults.c \
klu_diagnostics.c \
klu_dump.c \
klu_extract.c \
klu_factor.c \
klu_free_numeric.c \
klu_free_symbolic.c \
klu_kernel.c \
klu_memory.c \
klu_refactor.c \
klu_scale.c \
klu_solve.c \
klu_sort.c \
klu_tsolve.c
libKLU_complex_la_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/src/include -DCOMPLEX
libKLU_la_SOURCES = \
klusmp.c
libKLU_la_LIBADD = \
libKLU_real.la \
libKLU_complex.la
libKLU_la_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/src/include
MAINTAINERCLEANFILES = Makefile.in

748
src/maths/KLU/klusmp.c
View File

@ -0,0 +1,748 @@
/*
* Spice3 COMPATIBILITY MODULE
*
* Author: Advising professor:
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
* UC Berkeley
*
* This module contains routines that make Sparse1.3 a direct
* replacement for the SMP sparse matrix package in Spice3c1 or Spice3d1.
* Sparse1.3 is in general a faster and more robust package than SMP.
* These advantages become significant on large circuits.
*
* >>> User accessible functions contained in this file:
* SMPaddElt
* SMPmakeElt
* SMPcClear
* SMPclear
* SMPcLUfac
* SMPluFac
* SMPcReorder
* SMPreorder
* SMPcaSolve
* SMPcSolve
* SMPsolve
* SMPmatSize
* SMPnewMatrix
* SMPdestroy
* SMPpreOrder
* SMPprint
* SMPgetError
* SMPcProdDiag
* LoadGmin
* SMPfindElt
* SMPcombine
* SMPcCombine
*/
/*
* To replace SMP with Sparse, rename the file spSpice3.h to
* spMatrix.h and place Sparse in a subdirectory of SPICE called
* `sparse'. Then on UNIX compile Sparse by executing `make spice'.
* If not on UNIX, after compiling Sparse and creating the sparse.a
* archive, compile this file (spSMP.c) and spSMP.o to the archive,
* then copy sparse.a into the SPICE main directory and rename it
* SMP.a. Finally link SPICE.
*
* To be compatible with SPICE, the following Sparse compiler options
* (in spConfig.h) should be set as shown below:
*
* EXPANDABLE YES
* TRANSLATE NO
* INITIALIZE NO or YES, YES for use with test prog.
* DIAGONAL_PIVOTING YES
* MODIFIED_MARKOWITZ NO
* DELETE NO
* STRIP NO
* MODIFIED_NODAL YES
* QUAD_ELEMENT NO
* TRANSPOSE YES
* SCALING NO
* DOCUMENTATION YES
* MULTIPLICATION NO
* DETERMINANT YES
* STABILITY NO
* CONDITION NO
* PSEUDOCONDITION NO
* DEBUG YES
*
* spREAL double
*/
/*
* Revision and copyright information.
*
* Copyright (c) 1985,86,87,88,89,90
* by Kenneth S. Kundert and the University of California.
*
* Permission to use, copy, modify, and distribute this software and its
* documentation for any purpose and without fee is hereby granted, provided
* that the above copyright notice appear in all copies and supporting
* documentation and that the authors and the University of California
* are properly credited. The authors and the University of California
* make no representations as to the suitability of this software for
* any purpose. It is provided `as is', without express or implied warranty.
*/
/*
* IMPORTS
*
* >>> Import descriptions:
* spMatrix.h
* Sparse macros and declarations.
* SMPdefs.h
* Spice3's matrix macro definitions.
*/
#include "ngspice/config.h"
#include <assert.h>
#include <stdio.h>
#include <math.h>
#include "ngspice/spmatrix.h"
#include "../sparse/spdefs.h"
#include "ngspice/smpdefs.h"
#if defined (_MSC_VER)
extern double scalbn(double, int);
#define logb _logb
extern double logb(double);
#endif
static void LoadGmin_CSC (double **diag, int n, double Gmin) ;
static void LoadGmin (SMPmatrix *eMatrix, double Gmin) ;
void
SMPmatrix_CSC (SMPmatrix *Matrix)
{
spMatrix_CSC (Matrix->SPmatrix, Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx,
Matrix->CKTkluAx_Complex, Matrix->CKTkluN, Matrix->CKTbindStruct, Matrix->CKTdiag_CSC) ;
return ;
}
void
SMPnnz (SMPmatrix *Matrix)
{
Matrix->CKTklunz = Matrix->SPmatrix->Elements ;
return ;
}
/*
* SMPaddElt()
*/
int
SMPaddElt (SMPmatrix *Matrix, int Row, int Col, double Value)
{
*spGetElement (Matrix->SPmatrix, Row, Col) = Value ;
return spError (Matrix->SPmatrix) ;
}
/*
* SMPmakeElt()
*/
double *
SMPmakeElt (SMPmatrix *Matrix, int Row, int Col)
{
return spGetElement (Matrix->SPmatrix, Row, Col) ;
}
/*
* SMPcClear()
*/
void
SMPcClear (SMPmatrix *Matrix)
{
int i ;
if (Matrix->CKTkluMODE)
{
spClear (Matrix->SPmatrix) ;
if (Matrix->CKTkluAx_Complex != NULL)
{
for (i = 0 ; i < 2 * Matrix->CKTklunz ; i++)
Matrix->CKTkluAx_Complex [i] = 0 ;
}
} else {
spClear (Matrix->SPmatrix) ;
}
}
/*
* SMPclear()
*/
void
SMPclear (SMPmatrix *Matrix)
{
int i ;
if (Matrix->CKTkluMODE)
{
spClear (Matrix->SPmatrix) ;
if (Matrix->CKTkluAx != NULL)
{
for (i = 0 ; i < Matrix->CKTklunz ; i++)
Matrix->CKTkluAx [i] = 0 ;
}
} else {
spClear (Matrix->SPmatrix) ;
}
}
#define NG_IGNORE(x) (void)x
/*
* SMPcLUfac()
*/
/*ARGSUSED*/
int
SMPcLUfac (SMPmatrix *Matrix, double PivTol)
{
int ret ;
NG_IGNORE (PivTol) ;
if (Matrix->CKTkluMODE)
{
spSetComplex (Matrix->SPmatrix) ;
ret = klu_z_refactor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx_Complex,
Matrix->CKTkluSymbolic, Matrix->CKTkluNumeric, Matrix->CKTkluCommon) ;
return (!ret) ;
} else {
spSetComplex (Matrix->SPmatrix) ;
return spFactor (Matrix->SPmatrix) ;
}
}
/*
* SMPluFac()
*/
/*ARGSUSED*/
int
SMPluFac (SMPmatrix *Matrix, double PivTol, double Gmin)
{
int ret ;
NG_IGNORE (PivTol) ;
if (Matrix->CKTkluMODE)
{
spSetReal (Matrix->SPmatrix) ;
LoadGmin_CSC (Matrix->CKTdiag_CSC, Matrix->CKTkluN, Gmin) ;
ret = klu_refactor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx,
Matrix->CKTkluSymbolic, Matrix->CKTkluNumeric, Matrix->CKTkluCommon) ;
return (!ret) ;
// if (ret == 1)
// return 0 ;
// else if (ret == 0)
// return (E_SINGULAR) ;
// else {
// fprintf (stderr, "KLU Error in re-factor!") ;
// return 1 ;
// }
} else {
spSetReal (Matrix->SPmatrix) ;
LoadGmin (Matrix, Gmin) ;
return spFactor (Matrix->SPmatrix) ;
}
}
/*
* SMPcReorder()
*/
int
SMPcReorder (SMPmatrix *Matrix, double PivTol, double PivRel, int *NumSwaps)
{
if (Matrix->CKTkluMODE)
{
*NumSwaps = 1 ;
spSetComplex (Matrix->SPmatrix) ;
if (Matrix->CKTkluNumeric != NULL)
{
klu_z_free_numeric (&(Matrix->CKTkluNumeric), Matrix->CKTkluCommon) ;
Matrix->CKTkluNumeric = klu_z_factor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx_Complex, Matrix->CKTkluSymbolic, Matrix->CKTkluCommon) ;
} else
Matrix->CKTkluNumeric = klu_z_factor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx_Complex, Matrix->CKTkluSymbolic, Matrix->CKTkluCommon) ;
if (Matrix->CKTkluNumeric == NULL)
return 1 ;
else
return 0 ;
} else {
*NumSwaps = 1 ;
spSetComplex (Matrix->SPmatrix) ;
return spOrderAndFactor (Matrix->SPmatrix, NULL, (spREAL)PivRel, (spREAL)PivTol, YES) ;
}
}
/*
* SMPreorder()
*/
int
SMPreorder (SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin)
{
if (Matrix->CKTkluMODE)
{
spSetReal (Matrix->SPmatrix) ;
LoadGmin_CSC (Matrix->CKTdiag_CSC, Matrix->CKTkluN, Gmin) ;
if (Matrix->CKTkluNumeric != NULL)
{
klu_free_numeric (&(Matrix->CKTkluNumeric), Matrix->CKTkluCommon) ;
Matrix->CKTkluNumeric = klu_factor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx, Matrix->CKTkluSymbolic, Matrix->CKTkluCommon) ;
} else
Matrix->CKTkluNumeric = klu_factor (Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluAx, Matrix->CKTkluSymbolic, Matrix->CKTkluCommon) ;
if (Matrix->CKTkluNumeric == NULL)
return 1 ;
else
return 0 ;
} else {
spSetReal (Matrix->SPmatrix) ;
LoadGmin (Matrix, Gmin) ;
return spOrderAndFactor (Matrix->SPmatrix, NULL, (spREAL)PivRel, (spREAL)PivTol, YES) ;
}
}
/*
* SMPcaSolve()
*/
void
SMPcaSolve (SMPmatrix *Matrix, double RHS[], double iRHS[], double Spare[], double iSpare[])
{
printf ("SMPcaSolve\n") ;
NG_IGNORE (iSpare) ;
NG_IGNORE (Spare) ;
spSolveTransposed (Matrix->SPmatrix, RHS, RHS, iRHS, iRHS) ;
}
/*
* SMPcSolve()
*/
void
SMPcSolve (SMPmatrix *Matrix, double RHS[], double iRHS[], double Spare[], double iSpare[])
{
int ret, i, *pExtOrder ;
NG_IGNORE (iSpare) ;
NG_IGNORE (Spare) ;
if (Matrix->CKTkluMODE)
{
pExtOrder = &Matrix->SPmatrix->IntToExtRowMap [Matrix->CKTkluN] ;
for (i = 2 * Matrix->CKTkluN - 1 ; i > 0 ; i -= 2)
{
Matrix->CKTkluIntermediate_Complex [i] = RHS [*(pExtOrder)] ;
Matrix->CKTkluIntermediate_Complex [i - 1] = iRHS [*(pExtOrder--)] ;
}
ret = klu_z_solve (Matrix->CKTkluSymbolic, Matrix->CKTkluNumeric, Matrix->CKTkluN, 1, Matrix->CKTkluIntermediate_Complex, Matrix->CKTkluCommon) ;
pExtOrder = &Matrix->SPmatrix->IntToExtColMap [Matrix->CKTkluN] ;
for (i = 2 * Matrix->CKTkluN - 1 ; i > 0 ; i -= 2)
{
RHS [*(pExtOrder)] = Matrix->CKTkluIntermediate_Complex [i] ;
iRHS [*(pExtOrder--)] = Matrix->CKTkluIntermediate_Complex [i - 1] ;
}
} else {
spSolve (Matrix->SPmatrix, RHS, RHS, iRHS, iRHS) ;
}
}
/*
* SMPsolve()
*/
void
SMPsolve (SMPmatrix *Matrix, double RHS[], double Spare[])
{
int ret, i, *pExtOrder ;
NG_IGNORE (Spare) ;
if (Matrix->CKTkluMODE) {
pExtOrder = &Matrix->SPmatrix->IntToExtRowMap [Matrix->CKTkluN] ;
for (i = Matrix->CKTkluN - 1 ; i >= 0 ; i--)
Matrix->CKTkluIntermediate [i] = RHS [*(pExtOrder--)] ;
ret = klu_solve (Matrix->CKTkluSymbolic, Matrix->CKTkluNumeric, Matrix->CKTkluN, 1, Matrix->CKTkluIntermediate, Matrix->CKTkluCommon) ;
pExtOrder = &Matrix->SPmatrix->IntToExtColMap [Matrix->CKTkluN] ;
for (i = Matrix->CKTkluN - 1 ; i >= 0 ; i--)
RHS [*(pExtOrder--)] = Matrix->CKTkluIntermediate [i] ;
} else {
spSolve (Matrix->SPmatrix, RHS, RHS, NULL, NULL) ;
}
}
/*
* SMPmatSize()
*/
int
SMPmatSize (SMPmatrix *Matrix)
{
return spGetSize (Matrix->SPmatrix, 1) ;
}
/*
* SMPnewMatrix()
*/
int
SMPnewMatrix (SMPmatrix *Matrix, int size)
{
int Error ;
Matrix->SPmatrix = spCreate (size, 1, &Error) ;
return Error ;
}
/*
* SMPdestroy()
*/
void
SMPdestroy (SMPmatrix *Matrix)
{
if (Matrix->CKTkluMODE)
{
spDestroy (Matrix->SPmatrix) ;
klu_free_numeric (&(Matrix->CKTkluNumeric), Matrix->CKTkluCommon) ;
klu_free_symbolic (&(Matrix->CKTkluSymbolic), Matrix->CKTkluCommon) ;
free (Matrix->CKTkluAp) ;
free (Matrix->CKTkluAi) ;
free (Matrix->CKTkluAx) ;
free (Matrix->CKTdiag_CSC) ;
free (Matrix->CKTkluIntermediate) ;
free (Matrix->CKTkluIntermediate_Complex) ;
} else {
spDestroy (Matrix->SPmatrix) ;
}
}
/*
* SMPpreOrder()
*/
int
SMPpreOrder (SMPmatrix *Matrix)
{
if (Matrix->CKTkluMODE)
{
Matrix->CKTkluSymbolic = klu_analyze (Matrix->CKTkluN, Matrix->CKTkluAp, Matrix->CKTkluAi, Matrix->CKTkluCommon) ;
return 0 ;
} else {
spMNA_Preorder (Matrix->SPmatrix) ;
return spError (Matrix->SPmatrix) ;
}
}
/*
* SMPprintRHS()
*/
void
SMPprintRHS (SMPmatrix *Matrix, char *Filename, RealVector RHS, RealVector iRHS)
{
if (!Matrix->CKTkluMODE)
spFileVector (Matrix->SPmatrix, Filename, RHS, iRHS) ;
}
/*
* SMPprint()
*/
void
SMPprint (SMPmatrix *Matrix, char *Filename)
{
if (!Matrix->CKTkluMODE)
{
if (Filename)
spFileMatrix (Matrix->SPmatrix, Filename, "Circuit Matrix", 0, 1, 1) ;
else
spPrint (Matrix->SPmatrix, 0, 1, 1) ;
}
}
/*
* SMPgetError()
*/
void
SMPgetError (SMPmatrix *Matrix, int *Col, int *Row)
{
if (Matrix->CKTkluMODE)
{
*Row = Matrix->SPmatrix->IntToExtRowMap [Matrix->CKTkluCommon->singular_col] ;
*Col = Matrix->SPmatrix->IntToExtColMap [Matrix->CKTkluCommon->singular_col] ;
} else {
spWhereSingular (Matrix->SPmatrix, Row, Col) ;
}
}
/*
* SMPcProdDiag()
* note: obsolete for Spice3d2 and later
*/
int
SMPcProdDiag (SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
{
spDeterminant (Matrix->SPmatrix, pExponent, &(pMantissa->real), &(pMantissa->imag)) ;
return spError (Matrix->SPmatrix) ;
}
/*
* SMPcDProd()
*/
int
SMPcDProd (SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
{
double re, im, x, y, z;
int p;
spDeterminant (Matrix->SPmatrix, &p, &re, &im) ;
#ifndef M_LN2
#define M_LN2 0.69314718055994530942
#endif
#ifndef M_LN10
#define M_LN10 2.30258509299404568402
#endif
#ifdef debug_print
printf ("Determinant 10: (%20g,%20g)^%d\n", re, im, p) ;
#endif
/* Convert base 10 numbers to base 2 numbers, for comparison */
y = p * M_LN10 / M_LN2;
x = (int) y;
y -= x;
/* ASSERT
* x = integral part of exponent, y = fraction part of exponent
*/
/* Fold in the fractional part */
#ifdef debug_print
printf (" ** base10 -> base2 int = %g, frac = %20g\n", x, y) ;
#endif
z = pow (2.0, y) ;
re *= z ;
im *= z ;
#ifdef debug_print
printf (" ** multiplier = %20g\n", z) ;
#endif
/* 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 ;
} else if (im != 0.0) {
z = logb (im) ;
y = 0 ;
} else {
/* 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 ;
*pExponent = (int)(x + y) ;
x = scalbn (re, (int) -y) ;
z = scalbn (im, (int) -y) ;
#ifdef debug_print
printf (" ** values are: re %g, im %g, y %g, re' %g, im' %g\n", re, im, y, x, z) ;
#endif
pMantissa->real = scalbn (re, (int) -y) ;
pMantissa->imag = scalbn (im, (int) -y) ;
#ifdef debug_print
printf ("Determinant 10->2: (%20g,%20g)^%d\n", pMantissa->real, pMantissa->imag, *pExponent) ;
#endif
return spError (Matrix->SPmatrix) ;
}
/*
* The following routines need internal knowledge of the Sparse data
* structures.
*/
/*
* LOAD GMIN
*
* This routine adds Gmin to each diagonal element. Because Gmin is
* added to the current diagonal, which may bear little relation to
* what the outside world thinks is a diagonal, and because the
* elements that are diagonals may change after calling spOrderAndFactor,
* use of this routine is not recommended. It is included here simply
* for compatibility with Spice3.
*/
static void
LoadGmin_CSC (double **diag, int n, double Gmin)
{
int i ;
if (Gmin != 0.0)
for (i = 0 ; i < n ; i++)
if (diag [i] != NULL)
*(diag [i]) += Gmin ;
}
static void
LoadGmin (SMPmatrix *eMatrix, double Gmin)
{
MatrixPtr Matrix = eMatrix->SPmatrix ;
int I ;
ArrayOfElementPtrs Diag ;
ElementPtr diag ;
/* Begin `LoadGmin'. */
assert (IS_SPARSE (Matrix)) ;
if (Gmin != 0.0) {
Diag = Matrix->Diag ;
for (I = Matrix->Size ; I > 0 ; I--)
{
if ((diag = Diag [I]) != NULL)
diag->Real += Gmin ;
}
}
return ;
}
/*
* FIND ELEMENT
*
* This routine finds an element in the matrix by row and column number.
* If the element exists, a pointer to it is returned. If not, then NULL
* is returned unless the CreateIfMissing flag is TRUE, in which case a
* pointer to the new element is returned.
*/
SMPelement *
SMPfindElt (SMPmatrix *eMatrix, int Row, int Col, int CreateIfMissing)
{
MatrixPtr Matrix = eMatrix->SPmatrix ;
ElementPtr Element ;
/* Begin `SMPfindElt'. */
assert (IS_SPARSE (Matrix)) ;
Row = Matrix->ExtToIntRowMap [Row] ;
Col = Matrix->ExtToIntColMap [Col] ;
Element = Matrix->FirstInCol [Col] ;
Element = spcFindElementInCol (Matrix, &Element, Row, Col, CreateIfMissing) ;
return (SMPelement *)Element ;
}
/* XXX The following should probably be implemented in spUtils */
/*
* SMPcZeroCol()
*/
int
SMPcZeroCol (SMPmatrix *eMatrix, int Col)
{
MatrixPtr Matrix = eMatrix->SPmatrix ;
ElementPtr Element ;
Col = Matrix->ExtToIntColMap [Col] ;
for (Element = Matrix->FirstInCol [Col] ; Element != NULL ; Element = Element->NextInCol)
{
Element->Real = 0.0 ;
Element->Imag = 0.0 ;
}
return spError (Matrix) ;
}
/*
* SMPcAddCol()
*/
int
SMPcAddCol (SMPmatrix *eMatrix, int Accum_Col, int Addend_Col)
{
MatrixPtr Matrix = eMatrix->SPmatrix ;
ElementPtr Accum, Addend, *Prev ;
Accum_Col = Matrix->ExtToIntColMap [Accum_Col] ;
Addend_Col = Matrix->ExtToIntColMap [Addend_Col] ;
Addend = Matrix->FirstInCol [Addend_Col] ;
Prev = &Matrix->FirstInCol [Accum_Col] ;
Accum = *Prev;
while (Addend != NULL)
{
while (Accum && Accum->Row < Addend->Row)
{
Prev = &Accum->NextInCol ;
Accum = *Prev ;
}
if (!Accum || Accum->Row > Addend->Row)
{
Accum = spcCreateElement (Matrix, Addend->Row, Accum_Col, Prev, 0) ;
}
Accum->Real += Addend->Real ;
Accum->Imag += Addend->Imag ;
Addend = Addend->NextInCol ;
}
return spError (Matrix) ;
}
/*
* SMPzeroRow()
*/
int
SMPzeroRow (SMPmatrix *eMatrix, int Row)
{
MatrixPtr Matrix = eMatrix->SPmatrix ;
ElementPtr Element ;
Row = Matrix->ExtToIntColMap [Row] ;
if (Matrix->RowsLinked == NO)
spcLinkRows (Matrix) ;
if (Matrix->PreviousMatrixWasComplex || Matrix->Complex)
{
for (Element = Matrix->FirstInRow[Row] ; Element != NULL; Element = Element->NextInRow)
{
Element->Real = 0.0 ;
Element->Imag = 0.0 ;
}
} else {
for (Element = Matrix->FirstInRow [Row] ; Element != NULL ; Element = Element->NextInRow)
{
Element->Real = 0.0 ;
}
}
return spError (Matrix) ;
}

145
src/maths/sparse/spCSC.c
View File

@ -0,0 +1,145 @@
/* Sparse Matrix to CSC Matrix Conversion Routines
* Including Dump Routines
*
* Author: Francesco Lannutti 2011-2012
*
* Instructions:
* spMatrix_CSC_dump and spRHS_CSC_dump are the dump routines;
* insert them in a point in your code after that the Sparse Matrix
* is filled in to dump the whole matrix in the CSC format.
* To solve correctly the resulting CSC linear system, it's crucial
* to perform another inversion of the Solution Vector following this code:
*
* pExtOrder = IntToExtColMap [n] ;
* for (i = n - 1 ; i >= 0 ; i--)
* RHS [*(pExtOrder--)] = Intermediate [i] ;
*/
/* Includes */
#include "ngspice/spmatrix.h"
#include "spdefs.h"
/* Body */
int
WriteCol_original (MatrixPtr Matrix, int Col, spREAL *CSC_Element, spREAL *CSC_Element_Complex, int *CSC_Row, BindElement *BindSparseCSC, spREAL **diag)
{
int i ;
ElementPtr current ;
i = 0 ;
current = Matrix->FirstInCol [Col] ;
while (current != NULL) {
BindSparseCSC [i].Sparse = (double *)current ;
BindSparseCSC [i].CSC = &(CSC_Element [i]) ;
BindSparseCSC [i].CSC_Complex = &(CSC_Element_Complex [2 * i]) ;
CSC_Row [i] = (current->Row) - 1 ;
if (CSC_Row [i] == Col - 1)
diag [0] = &(CSC_Element [i]) ;
i++ ;
current = current->NextInCol ;
}
return i ;
}
int
WriteCol_original_dump (MatrixPtr Matrix, int Col, spREAL *CSC_Element, int *CSC_Row)
{
int i ;
ElementPtr current ;
i = 0 ;
current = Matrix->FirstInCol [Col] ;
while (current != NULL) {
CSC_Element [i] = current->Real ;
CSC_Row [i] = (current->Row) - 1 ;
i++ ;
current = current->NextInCol ;
}
return i ;
}
void
spMatrix_CSC (MatrixPtr Matrix, int *Ap, int *Ai, double *Ax, double *Ax_Complex, int n, BindElement *BindSparseCSC, double **diag)
{
int offset, i ;
offset = 0 ;
Ap[0] = offset ;
for (i = 1 ; i <= n ; i++) {
offset += WriteCol_original (Matrix, i, (spREAL *)(Ax + offset), (spREAL *)(Ax_Complex + 2 * offset),
(int *)(Ai + offset), BindSparseCSC + offset, (spREAL **)(diag + (i - 1))) ;
Ap[i] = offset ;
}
}
void
spMatrix_CSC_dump (MatrixPtr Matrix, char *CSC_output)
{
FILE *output ;
int offset, i, j, *Ap, *Ai, n, nz ;
double *Ax ;
n = spGetSize (Matrix, 1) ;
nz = Matrix->Elements ;
Ap = (int *) SP_MALLOC (int, n + 1) ;
Ai = (int *) SP_MALLOC (int, nz) ;
Ax = (double *) SP_MALLOC (double, nz) ;
offset = 0 ;
Ap[0] = offset ;
for (i = 1 ; i <= n ; i++) {
offset += WriteCol_original_dump (Matrix, i, (spREAL *)(Ax + offset), (int *)(Ai + offset)) ;
Ap[i] = offset ;
}
output = fopen (CSC_output, "w") ;
fprintf (output, "%%%%MatrixMarket matrix coordinate real general\n") ;
fprintf (output, "%%-------------------------------------------------------------------------------\n") ;
fprintf (output, "%% Transient Matrix Dump\n%% Family: ISCAS Circuit\n") ;
fprintf (output, "%%-------------------------------------------------------------------------------\n") ;
fprintf (output, "%d %d %d\n", n, n, offset) ;
for (i = 0 ; i < n ; i++)
for (j = Ap [i] ; j < Ap [i + 1] ; j++)
fprintf (output, "%d %d %-.9g\n", Ai [j] + 1, i + 1, Ax [j]) ;
fclose (output) ;
output = fopen ("IntToExtColMap.txt", "w") ;
for (i = 1 ; i <= n ; i++)
fprintf (output, "%d\n", Matrix->IntToExtColMap [i]) ;
fclose (output) ;
SP_FREE (Ap) ;
SP_FREE (Ai) ;
SP_FREE (Ax) ;
}
void
spRHS_CSC_dump (RealNumber *RHS, char *CSC_output_b, MatrixPtr Matrix)
{
FILE *output ;
int i, n, *pExtOrder ;
double *Intermediate ;
n = spGetSize (Matrix, 1) ;
Intermediate = (double *) SP_MALLOC (double, n) ;
pExtOrder = &Matrix->IntToExtRowMap [n] ;
for (i = n - 1 ; i >= 0 ; i--)
Intermediate [i] = RHS [*(pExtOrder--)] ;
output = fopen (CSC_output_b, "w") ;
fprintf (output, "%%%%MatrixMarket matrix array real general\n") ;
fprintf (output, "%%-------------------------------------------------------------------------------\n") ;
fprintf (output, "%% Transient RHS Vector Dump\n%% Family: ISCAS Circuit\n") ;
fprintf (output, "%%-------------------------------------------------------------------------------\n") ;
fprintf (output, "%d %d\n", n, 1) ;
for (i = 1 ; i < n + 1 ; i++)
fprintf (output, "%-.9g\n", Intermediate [i]) ;
fclose (output) ;
SP_FREE (Intermediate) ;
}
Write Preview
Loading…
Cancel
Save