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

KLU Integration from scratch #3, new files

This commit is contained in:
Francesco Lannutti 2014-04-21 14:11:28 +02:00 committed by rlar
parent 4705569c2e
commit 17d9d61af6
6 changed files with 2057 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

151
src/include/ngspice/UFconfig.h Normal file
View File

@ -0,0 +1,151 @@
/* ========================================================================== */
/* === UFconfig.h =========================================================== */
/* ========================================================================== */
/* Configuration file for SuiteSparse: a Suite of Sparse matrix packages
* (AMD, COLAMD, CCOLAMD, CAMD, CHOLMOD, UMFPACK, CXSparse, and others).
*
* UFconfig.h provides the definition of the long integer. On most systems,
* a C program can be compiled in LP64 mode, in which long's and pointers are
* both 64-bits, and int's are 32-bits. Windows 64, however, uses the LLP64
* model, in which int's and long's are 32-bits, and long long's and pointers
* are 64-bits.
*
* SuiteSparse packages that include long integer versions are
* intended for the LP64 mode. However, as a workaround for Windows 64
* (and perhaps other systems), the long integer can be redefined.
*
* If _WIN64 is defined, then the __int64 type is used instead of long.
*
* The long integer can also be defined at compile time. For example, this
* could be added to UFconfig.mk:
*
* CFLAGS = -O -D'UF_long=long long' -D'UF_long_max=9223372036854775801' \
* -D'UF_long_idd="lld"'
*
* This file defines UF_long as either long (on all but _WIN64) or
* __int64 on Windows 64. The intent is that a UF_long is always a 64-bit
* integer in a 64-bit code. ptrdiff_t might be a better choice than long;
* it is always the same size as a pointer.
*
* This file also defines the SUITESPARSE_VERSION and related definitions.
*
* Copyright (c) 2007, University of Florida. No licensing restrictions
* apply to this file or to the UFconfig directory. Author: Timothy A. Davis.
*/
#ifndef _UFCONFIG_H
#define _UFCONFIG_H
#ifdef __cplusplus
extern "C" {
#endif
#include <limits.h>
#include <stdlib.h>
/* ========================================================================== */
/* === UF_long ============================================================== */
/* ========================================================================== */
#ifndef UF_long
#ifdef _WIN64
#define UF_long __int64
#define UF_long_max _I64_MAX
#define UF_long_idd "I64d"
#else
#define UF_long long
#define UF_long_max LONG_MAX
#define UF_long_idd "ld"
#endif
#define UF_long_id "%" UF_long_idd
#endif
/* ========================================================================== */
/* === UFconfig parameters and functions ==================================== */
/* ========================================================================== */
/* SuiteSparse-wide parameters will be placed in this struct. So far, they
are only used by RBio. */
typedef struct UFconfig_struct
{
void *(*malloc_memory) (size_t) ; /* pointer to malloc */
void *(*realloc_memory) (void *, size_t) ; /* pointer to realloc */
void (*free_memory) (void *) ; /* pointer to free */
void *(*calloc_memory) (size_t, size_t) ; /* pointer to calloc */
} UFconfig ;
void *UFmalloc /* pointer to allocated block of memory */
(
size_t nitems, /* number of items to malloc (>=1 is enforced) */
size_t size_of_item, /* sizeof each item */
int *ok, /* TRUE if successful, FALSE otherwise */
UFconfig *config /* SuiteSparse-wide configuration */
) ;
void *UFfree /* always returns NULL */
(
void *p, /* block to free */
UFconfig *config /* SuiteSparse-wide configuration */
) ;
/* ========================================================================== */
/* === SuiteSparse version ================================================== */
/* ========================================================================== */
/* SuiteSparse is not a package itself, but a collection of packages, some of
* which must be used together (UMFPACK requires AMD, CHOLMOD requires AMD,
* COLAMD, CAMD, and CCOLAMD, etc). A version number is provided here for the
* collection itself. The versions of packages within each version of
* SuiteSparse are meant to work together. Combining one packge from one
* version of SuiteSparse, with another package from another version of
* SuiteSparse, may or may not work.
*
* SuiteSparse Version 3.7.0 contains the following packages:
*
* UFconfig version 3.7.0 (version always the same as SuiteSparse)
* AMD version 2.2.3
* CAMD version 2.2.3
* CCOLAMD version 2.7.4
* COLAMD version 2.7.4
* BTF version 1.1.3
* CHOLMOD version 1.7.4
* CSparse3 version 3.0.2
* CSparse version 2.2.6
* CXSparse version 2.2.6
* KLU version 1.1.3
* LDL version 2.0.4
* RBio version 2.0.2
* SPQR version 1.2.3 (also called SuiteSparseQR)
* UFcollection version 1.6.0
* UMFPACK version 5.5.2
* SSMULT version 2.0.3
* spqr_rank version 1.0.0
* MATLAB_Tools various packages & M-files. No specific version number.
*
* Other package dependencies:
* BLAS required by CHOLMOD and UMFPACK
* LAPACK required by CHOLMOD
* METIS 4.0.1 required by CHOLMOD (optional) and KLU (optional)
*/
#define SUITESPARSE_DATE "Dec 15, 2011"
#define SUITESPARSE_VER_CODE(main,sub) ((main) * 1000 + (sub))
#define SUITESPARSE_MAIN_VERSION 3
#define SUITESPARSE_SUB_VERSION 7
#define SUITESPARSE_SUBSUB_VERSION 0
#define SUITESPARSE_VERSION \
SUITESPARSE_VER_CODE(SUITESPARSE_MAIN_VERSION,SUITESPARSE_SUB_VERSION)
#ifdef __cplusplus
}
#endif
#endif

412
src/include/ngspice/amd.h Normal file
View File

@ -0,0 +1,412 @@
/* ========================================================================= */
/* === AMD: approximate minimum degree ordering =========================== */
/* ========================================================================= */
/* ------------------------------------------------------------------------- */
/* AMD Version 2.2, Copyright (c) 2007 by Timothy A. Davis, */
/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
/* email: davis at cise.ufl.edu CISE Department, Univ. of Florida. */
/* web: http://www.cise.ufl.edu/research/sparse/amd */
/* ------------------------------------------------------------------------- */
/* AMD finds a symmetric ordering P of a matrix A so that the Cholesky
* factorization of P*A*P' has fewer nonzeros and takes less work than the
* Cholesky factorization of A. If A is not symmetric, then it performs its
* ordering on the matrix A+A'. Two sets of user-callable routines are
* provided, one for int integers and the other for UF_long integers.
*
* The method is based on the approximate minimum degree algorithm, discussed
* in Amestoy, Davis, and Duff, "An approximate degree ordering algorithm",
* SIAM Journal of Matrix Analysis and Applications, vol. 17, no. 4, pp.
* 886-905, 1996. This package can perform both the AMD ordering (with
* aggressive absorption), and the AMDBAR ordering (without aggressive
* absorption) discussed in the above paper. This package differs from the
* Fortran codes discussed in the paper:
*
* (1) it can ignore "dense" rows and columns, leading to faster run times
* (2) it computes the ordering of A+A' if A is not symmetric
* (3) it is followed by a depth-first post-ordering of the assembly tree
* (or supernodal elimination tree)
*
* For historical reasons, the Fortran versions, amd.f and amdbar.f, have
* been left (nearly) unchanged. They compute the identical ordering as
* described in the above paper.
*/
#ifndef AMD_H
#define AMD_H
/* make it easy for C++ programs to include AMD */
#ifdef __cplusplus
extern "C" {
#endif
/* get the definition of size_t: */
#include <stddef.h>
/* define UF_long */
#include "UFconfig.h"
int amd_order /* returns AMD_OK, AMD_OK_BUT_JUMBLED,
* AMD_INVALID, or AMD_OUT_OF_MEMORY */
(
int n, /* A is n-by-n. n must be >= 0. */
const int Ap [ ], /* column pointers for A, of size n+1 */
const int Ai [ ], /* row indices of A, of size nz = Ap [n] */
int P [ ], /* output permutation, of size n */
double Control [ ], /* input Control settings, of size AMD_CONTROL */
double Info [ ] /* output Info statistics, of size AMD_INFO */
) ;
UF_long amd_l_order /* see above for description of arguments */
(
UF_long n,
const UF_long Ap [ ],
const UF_long Ai [ ],
UF_long P [ ],
double Control [ ],
double Info [ ]
) ;
/* Input arguments (not modified):
*
* n: the matrix A is n-by-n.
* Ap: an int/UF_long array of size n+1, containing column pointers of A.
* Ai: an int/UF_long array of size nz, containing the row indices of A,
* where nz = Ap [n].
* Control: a double array of size AMD_CONTROL, containing control
* parameters. Defaults are used if Control is NULL.
*
* Output arguments (not defined on input):
*
* P: an int/UF_long array of size n, containing the output permutation. If
* row i is the kth pivot row, then P [k] = i. In MATLAB notation,
* the reordered matrix is A (P,P).
* Info: a double array of size AMD_INFO, containing statistical
* information. Ignored if Info is NULL.
*
* On input, the matrix A is stored in column-oriented form. The row indices
* of nonzero entries in column j are stored in Ai [Ap [j] ... Ap [j+1]-1].
*
* If the row indices appear in ascending order in each column, and there
* are no duplicate entries, then amd_order is slightly more efficient in
* terms of time and memory usage. If this condition does not hold, a copy
* of the matrix is created (where these conditions do hold), and the copy is
* ordered. This feature is new to v2.0 (v1.2 and earlier required this
* condition to hold for the input matrix).
*
* Row indices must be in the range 0 to
* n-1. Ap [0] must be zero, and thus nz = Ap [n] is the number of nonzeros
* in A. The array Ap is of size n+1, and the array Ai is of size nz = Ap [n].
* The matrix does not need to be symmetric, and the diagonal does not need to
* be present (if diagonal entries are present, they are ignored except for
* the output statistic Info [AMD_NZDIAG]). The arrays Ai and Ap are not
* modified. This form of the Ap and Ai arrays to represent the nonzero
* pattern of the matrix A is the same as that used internally by MATLAB.
* If you wish to use a more flexible input structure, please see the
* umfpack_*_triplet_to_col routines in the UMFPACK package, at
* http://www.cise.ufl.edu/research/sparse/umfpack.
*
* Restrictions: n >= 0. Ap [0] = 0. Ap [j] <= Ap [j+1] for all j in the
* range 0 to n-1. nz = Ap [n] >= 0. Ai [0..nz-1] must be in the range 0
* to n-1. Finally, Ai, Ap, and P must not be NULL. If any of these
* restrictions are not met, AMD returns AMD_INVALID.
*
* AMD returns:
*
* AMD_OK if the matrix is valid and sufficient memory can be allocated to
* perform the ordering.
*
* AMD_OUT_OF_MEMORY if not enough memory can be allocated.
*
* AMD_INVALID if the input arguments n, Ap, Ai are invalid, or if P is
* NULL.
*
* AMD_OK_BUT_JUMBLED if the matrix had unsorted columns, and/or duplicate
* entries, but was otherwise valid.
*
* The AMD routine first forms the pattern of the matrix A+A', and then
* computes a fill-reducing ordering, P. If P [k] = i, then row/column i of
* the original is the kth pivotal row. In MATLAB notation, the permuted
* matrix is A (P,P), except that 0-based indexing is used instead of the
* 1-based indexing in MATLAB.
*
* The Control array is used to set various parameters for AMD. If a NULL
* pointer is passed, default values are used. The Control array is not
* modified.
*
* Control [AMD_DENSE]: controls the threshold for "dense" rows/columns.
* A dense row/column in A+A' can cause AMD to spend a lot of time in
* ordering the matrix. If Control [AMD_DENSE] >= 0, rows/columns
* with more than Control [AMD_DENSE] * sqrt (n) entries are ignored
* during the ordering, and placed last in the output order. The
* default value of Control [AMD_DENSE] is 10. If negative, no
* rows/columns are treated as "dense". Rows/columns with 16 or
* fewer off-diagonal entries are never considered "dense".
*
* Control [AMD_AGGRESSIVE]: controls whether or not to use aggressive
* absorption, in which a prior element is absorbed into the current
* element if is a subset of the current element, even if it is not
* adjacent to the current pivot element (refer to Amestoy, Davis,
* & Duff, 1996, for more details). The default value is nonzero,
* which means to perform aggressive absorption. This nearly always
* leads to a better ordering (because the approximate degrees are
* more accurate) and a lower execution time. There are cases where
* it can lead to a slightly worse ordering, however. To turn it off,
* set Control [AMD_AGGRESSIVE] to 0.
*
* Control [2..4] are not used in the current version, but may be used in
* future versions.
*
* The Info array provides statistics about the ordering on output. If it is
* not present, the statistics are not returned. This is not an error
* condition.
*
* Info [AMD_STATUS]: the return value of AMD, either AMD_OK,
* AMD_OK_BUT_JUMBLED, AMD_OUT_OF_MEMORY, or AMD_INVALID.
*
* Info [AMD_N]: n, the size of the input matrix
*
* Info [AMD_NZ]: the number of nonzeros in A, nz = Ap [n]
*
* Info [AMD_SYMMETRY]: the symmetry of the matrix A. It is the number
* of "matched" off-diagonal entries divided by the total number of
* off-diagonal entries. An entry A(i,j) is matched if A(j,i) is also
* an entry, for any pair (i,j) for which i != j. In MATLAB notation,
* S = spones (A) ;
* B = tril (S, -1) + triu (S, 1) ;
* symmetry = nnz (B & B') / nnz (B) ;
*
* Info [AMD_NZDIAG]: the number of entries on the diagonal of A.
*
* Info [AMD_NZ_A_PLUS_AT]: the number of nonzeros in A+A', excluding the
* diagonal. If A is perfectly symmetric (Info [AMD_SYMMETRY] = 1)
* with a fully nonzero diagonal, then Info [AMD_NZ_A_PLUS_AT] = nz-n
* (the smallest possible value). If A is perfectly unsymmetric
* (Info [AMD_SYMMETRY] = 0, for an upper triangular matrix, for
* example) with no diagonal, then Info [AMD_NZ_A_PLUS_AT] = 2*nz
* (the largest possible value).
*
* Info [AMD_NDENSE]: the number of "dense" rows/columns of A+A' that were
* removed from A prior to ordering. These are placed last in the
* output order P.
*
* Info [AMD_MEMORY]: the amount of memory used by AMD, in bytes. In the
* current version, this is 1.2 * Info [AMD_NZ_A_PLUS_AT] + 9*n
* times the size of an integer. This is at most 2.4nz + 9n. This
* excludes the size of the input arguments Ai, Ap, and P, which have
* a total size of nz + 2*n + 1 integers.
*
* Info [AMD_NCMPA]: the number of garbage collections performed.
*
* Info [AMD_LNZ]: the number of nonzeros in L (excluding the diagonal).
* This is a slight upper bound because mass elimination is combined
* with the approximate degree update. It is a rough upper bound if
* there are many "dense" rows/columns. The rest of the statistics,
* below, are also slight or rough upper bounds, for the same reasons.
* The post-ordering of the assembly tree might also not exactly
* correspond to a true elimination tree postordering.
*
* Info [AMD_NDIV]: the number of divide operations for a subsequent LDL'
* or LU factorization of the permuted matrix A (P,P).
*
* Info [AMD_NMULTSUBS_LDL]: the number of multiply-subtract pairs for a
* subsequent LDL' factorization of A (P,P).
*
* Info [AMD_NMULTSUBS_LU]: the number of multiply-subtract pairs for a
* subsequent LU factorization of A (P,P), assuming that no numerical
* pivoting is required.
*
* Info [AMD_DMAX]: the maximum number of nonzeros in any column of L,
* including the diagonal.
*
* Info [14..19] are not used in the current version, but may be used in
* future versions.
*/
/* ------------------------------------------------------------------------- */
/* direct interface to AMD */
/* ------------------------------------------------------------------------- */
/* amd_2 is the primary AMD ordering routine. It is not meant to be
* user-callable because of its restrictive inputs and because it destroys
* the user's input matrix. It does not check its inputs for errors, either.
* However, if you can work with these restrictions it can be faster than
* amd_order and use less memory (assuming that you can create your own copy
* of the matrix for AMD to destroy). Refer to AMD/Source/amd_2.c for a
* description of each parameter. */
void amd_2
(
int n,
int Pe [ ],
int Iw [ ],
int Len [ ],
int iwlen,
int pfree,
int Nv [ ],
int Next [ ],
int Last [ ],
int Head [ ],
int Elen [ ],
int Degree [ ],
int W [ ],
double Control [ ],
double Info [ ]
) ;
void amd_l2
(
UF_long n,
UF_long Pe [ ],
UF_long Iw [ ],
UF_long Len [ ],
UF_long iwlen,
UF_long pfree,
UF_long Nv [ ],
UF_long Next [ ],
UF_long Last [ ],
UF_long Head [ ],
UF_long Elen [ ],
UF_long Degree [ ],
UF_long W [ ],
double Control [ ],
double Info [ ]
) ;
/* ------------------------------------------------------------------------- */
/* amd_valid */
/* ------------------------------------------------------------------------- */
/* Returns AMD_OK or AMD_OK_BUT_JUMBLED if the matrix is valid as input to
* amd_order; the latter is returned if the matrix has unsorted and/or
* duplicate row indices in one or more columns. Returns AMD_INVALID if the
* matrix cannot be passed to amd_order. For amd_order, the matrix must also
* be square. The first two arguments are the number of rows and the number
* of columns of the matrix. For its use in AMD, these must both equal n.
*
* NOTE: this routine returned TRUE/FALSE in v1.2 and earlier.
*/
int amd_valid
(
int n_row, /* # of rows */
int n_col, /* # of columns */
const int Ap [ ], /* column pointers, of size n_col+1 */
const int Ai [ ] /* row indices, of size Ap [n_col] */
) ;
UF_long amd_l_valid
(
UF_long n_row,
UF_long n_col,
const UF_long Ap [ ],
const UF_long Ai [ ]
) ;
/* ------------------------------------------------------------------------- */
/* AMD memory manager and printf routines */
/* ------------------------------------------------------------------------- */
/* The user can redefine these to change the malloc, free, and printf routines
* that AMD uses. */
#ifndef EXTERN
#define EXTERN extern
#endif
EXTERN void *(*amd_malloc) (size_t) ; /* pointer to malloc */
EXTERN void (*amd_free) (void *) ; /* pointer to free */
EXTERN void *(*amd_realloc) (void *, size_t) ; /* pointer to realloc */
EXTERN void *(*amd_calloc) (size_t, size_t) ; /* pointer to calloc */
EXTERN int (*amd_printf) (const char *, ...) ; /* pointer to printf */
/* ------------------------------------------------------------------------- */
/* AMD Control and Info arrays */
/* ------------------------------------------------------------------------- */
/* amd_defaults: sets the default control settings */
void amd_defaults (double Control [ ]) ;
void amd_l_defaults (double Control [ ]) ;
/* amd_control: prints the control settings */
void amd_control (double Control [ ]) ;
void amd_l_control (double Control [ ]) ;
/* amd_info: prints the statistics */
void amd_info (double Info [ ]) ;
void amd_l_info (double Info [ ]) ;
#define AMD_CONTROL 5 /* size of Control array */
#define AMD_INFO 20 /* size of Info array */
/* contents of Control */
#define AMD_DENSE 0 /* "dense" if degree > Control [0] * sqrt (n) */
#define AMD_AGGRESSIVE 1 /* do aggressive absorption if Control [1] != 0 */
/* default Control settings */
#define AMD_DEFAULT_DENSE 10.0 /* default "dense" degree 10*sqrt(n) */
#define AMD_DEFAULT_AGGRESSIVE 1 /* do aggressive absorption by default */
/* contents of Info */
#define AMD_STATUS 0 /* return value of amd_order and amd_l_order */
#define AMD_N 1 /* A is n-by-n */
#define AMD_NZ 2 /* number of nonzeros in A */
#define AMD_SYMMETRY 3 /* symmetry of pattern (1 is sym., 0 is unsym.) */
#define AMD_NZDIAG 4 /* # of entries on diagonal */
#define AMD_NZ_A_PLUS_AT 5 /* nz in A+A' */
#define AMD_NDENSE 6 /* number of "dense" rows/columns in A */
#define AMD_MEMORY 7 /* amount of memory used by AMD */
#define AMD_NCMPA 8 /* number of garbage collections in AMD */
#define AMD_LNZ 9 /* approx. nz in L, excluding the diagonal */
#define AMD_NDIV 10 /* number of fl. point divides for LU and LDL' */
#define AMD_NMULTSUBS_LDL 11 /* number of fl. point (*,-) pairs for LDL' */
#define AMD_NMULTSUBS_LU 12 /* number of fl. point (*,-) pairs for LU */
#define AMD_DMAX 13 /* max nz. in any column of L, incl. diagonal */
/* ------------------------------------------------------------------------- */
/* return values of AMD */
/* ------------------------------------------------------------------------- */
#define AMD_OK 0 /* success */
#define AMD_OUT_OF_MEMORY -1 /* malloc failed, or problem too large */
#define AMD_INVALID -2 /* input arguments are not valid */
#define AMD_OK_BUT_JUMBLED 1 /* input matrix is OK for amd_order, but
* columns were not sorted, and/or duplicate entries were present. AMD had
* to do extra work before ordering the matrix. This is a warning, not an
* error. */
/* ========================================================================== */
/* === AMD version ========================================================== */
/* ========================================================================== */
/* AMD Version 1.2 and later include the following definitions.
* As an example, to test if the version you are using is 1.2 or later:
*
* #ifdef AMD_VERSION
* if (AMD_VERSION >= AMD_VERSION_CODE (1,2)) ...
* #endif
*
* This also works during compile-time:
*
* #if defined(AMD_VERSION) && (AMD_VERSION >= AMD_VERSION_CODE (1,2))
* printf ("This is version 1.2 or later\n") ;
* #else
* printf ("This is an early version\n") ;
* #endif
*
* Versions 1.1 and earlier of AMD do not include a #define'd version number.
*/
#define AMD_DATE "Dec 7, 2011"
#define AMD_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
#define AMD_MAIN_VERSION 2
#define AMD_SUB_VERSION 2
#define AMD_SUBSUB_VERSION 3
#define AMD_VERSION AMD_VERSION_CODE(AMD_MAIN_VERSION,AMD_SUB_VERSION)
#ifdef __cplusplus
}
#endif
#endif

263
src/include/ngspice/btf.h Normal file
View File

@ -0,0 +1,263 @@
/* ========================================================================== */
/* === BTF package ========================================================== */
/* ========================================================================== */
/* BTF_MAXTRANS: find a column permutation Q to give A*Q a zero-free diagonal
* BTF_STRONGCOMP: find a symmetric permutation P to put P*A*P' into block
* upper triangular form.
* BTF_ORDER: do both of the above (btf_maxtrans then btf_strongcomp).
*
* Copyright (c) 2004-2007. Tim Davis, University of Florida,
* with support from Sandia National Laboratories. All Rights Reserved.
*/
/* ========================================================================== */
/* === BTF_MAXTRANS ========================================================= */
/* ========================================================================== */
/* BTF_MAXTRANS: finds a permutation of the columns of a matrix so that it has a
* zero-free diagonal. The input is an m-by-n sparse matrix in compressed
* column form. The array Ap of size n+1 gives the starting and ending
* positions of the columns in the array Ai. Ap[0] must be zero. The array Ai
* contains the row indices of the nonzeros of the matrix A, and is of size
* Ap[n]. The row indices of column j are located in Ai[Ap[j] ... Ap[j+1]-1].
* Row indices must be in the range 0 to m-1. Duplicate entries may be present
* in any given column. The input matrix is not checked for validity (row
* indices out of the range 0 to m-1 will lead to an undeterminate result -
* possibly a core dump, for example). Row indices in any given column need
* not be in sorted order. However, if they are sorted and the matrix already
* has a zero-free diagonal, then the identity permutation is returned.
*
* The output of btf_maxtrans is an array Match of size n. If row i is matched
* with column j, then A(i,j) is nonzero, and then Match[i] = j. If the matrix
* is structurally nonsingular, all entries in the Match array are unique, and
* Match can be viewed as a column permutation if A is square. That is, column
* k of the original matrix becomes column Match[k] of the permuted matrix. In
* MATLAB, this can be expressed as (for non-structurally singular matrices):
*
* Match = maxtrans (A) ;
* B = A (:, Match) ;
*
* except of course here the A matrix and Match vector are all 0-based (rows
* and columns in the range 0 to n-1), not 1-based (rows/cols in range 1 to n).
* The MATLAB dmperm routine returns a row permutation. See the maxtrans
* mexFunction for more details.
*
* If row i is not matched to any column, then Match[i] is == -1. The
* btf_maxtrans routine returns the number of nonzeros on diagonal of the
* permuted matrix.
*
* In the MATLAB mexFunction interface to btf_maxtrans, 1 is added to the Match
* array to obtain a 1-based permutation. Thus, in MATLAB where A is m-by-n:
*
* q = maxtrans (A) ; % has entries in the range 0:n
* q % a column permutation (only if sprank(A)==n)
* B = A (:, q) ; % permuted matrix (only if sprank(A)==n)
* sum (q > 0) ; % same as "sprank (A)"
*
* This behaviour differs from p = dmperm (A) in MATLAB, which returns the
* matching as p(j)=i if row i and column j are matched, and p(j)=0 if column j
* is unmatched.
*
* p = dmperm (A) ; % has entries in the range 0:m
* p % a row permutation (only if sprank(A)==m)
* B = A (p, :) ; % permuted matrix (only if sprank(A)==m)
* sum (p > 0) ; % definition of sprank (A)
*
* This algorithm is based on the paper "On Algorithms for obtaining a maximum
* transversal" by Iain Duff, ACM Trans. Mathematical Software, vol 7, no. 1,
* pp. 315-330, and "Algorithm 575: Permutations for a zero-free diagonal",
* same issue, pp. 387-390. Algorithm 575 is MC21A in the Harwell Subroutine
* Library. This code is not merely a translation of the Fortran code into C.
* It is a completely new implementation of the basic underlying method (depth
* first search over a subgraph with nodes corresponding to columns matched so
* far, and cheap matching). This code was written with minimal observation of
* the MC21A/B code itself. See comments below for a comparison between the
* maxtrans and MC21A/B codes.
*
* This routine operates on a column-form matrix and produces a column
* permutation. MC21A uses a row-form matrix and produces a row permutation.
* The difference is merely one of convention in the comments and interpretation
* of the inputs and outputs. If you want a row permutation, simply pass a
* compressed-row sparse matrix to this routine and you will get a row
* permutation (just like MC21A). Similarly, you can pass a column-oriented
* matrix to MC21A and it will happily return a column permutation.
*/
#ifndef _BTF_H
#define _BTF_H
/* make it easy for C++ programs to include BTF */
#ifdef __cplusplus
extern "C" {
#endif
#include "UFconfig.h"
int btf_maxtrans /* returns # of columns matched */
(
/* --- input, not modified: --- */
int nrow, /* A is nrow-by-ncol in compressed column form */
int ncol,
int Ap [ ], /* size ncol+1 */
int Ai [ ], /* size nz = Ap [ncol] */
double maxwork, /* maximum amount of work to do is maxwork*nnz(A); no limit
* if <= 0 */
/* --- output, not defined on input --- */
double *work, /* work = -1 if maxwork > 0 and the total work performed
* reached the maximum of maxwork*nnz(A).
* Otherwise, work = the total work performed. */
int Match [ ], /* size nrow. Match [i] = j if column j matched to row i
* (see above for the singular-matrix case) */
/* --- workspace, not defined on input or output --- */
int Work [ ] /* size 5*ncol */
) ;
/* long integer version (all "int" parameters become "UF_long") */
UF_long btf_l_maxtrans (UF_long, UF_long, UF_long *, UF_long *, double,
double *, UF_long *, UF_long *) ;
/* ========================================================================== */
/* === BTF_STRONGCOMP ======================================================= */
/* ========================================================================== */
/* BTF_STRONGCOMP finds the strongly connected components of a graph, returning
* a symmetric permutation. The matrix A must be square, and is provided on
* input in compressed-column form (see BTF_MAXTRANS, above). The diagonal of
* the input matrix A (or A*Q if Q is provided on input) is ignored.
*
* If Q is not NULL on input, then the strongly connected components of A*Q are
* found. Q may be flagged on input, where Q[k] < 0 denotes a flagged column k.
* The permutation is j = BTF_UNFLIP (Q [k]). On output, Q is modified (the
* flags are preserved) so that P*A*Q is in block upper triangular form.
*
* If Q is NULL, then the permutation P is returned so that P*A*P' is in upper
* block triangular form.
*
* The vector R gives the block boundaries, where block b is in rows/columns
* R[b] to R[b+1]-1 of the permuted matrix, and where b ranges from 1 to the
* number of strongly connected components found.
*/
int btf_strongcomp /* return # of strongly connected components */
(
/* input, not modified: */
int n, /* A is n-by-n in compressed column form */
int Ap [ ], /* size n+1 */
int Ai [ ], /* size nz = Ap [n] */
/* optional input, modified (if present) on output: */
int Q [ ], /* size n, input column permutation */
/* output, not defined on input */
int P [ ], /* size n. P [k] = j if row and column j are kth row/col
* in permuted matrix. */
int R [ ], /* size n+1. block b is in rows/cols R[b] ... R[b+1]-1 */
/* workspace, not defined on input or output */
int Work [ ] /* size 4n */
) ;
UF_long btf_l_strongcomp (UF_long, UF_long *, UF_long *, UF_long *, UF_long *,
UF_long *, UF_long *) ;
/* ========================================================================== */
/* === BTF_ORDER ============================================================ */
/* ========================================================================== */
/* BTF_ORDER permutes a square matrix into upper block triangular form. It
* does this by first finding a maximum matching (or perhaps a limited matching
* if the work is limited), via the btf_maxtrans function. If a complete
* matching is not found, BTF_ORDER completes the permutation, but flags the
* columns of P*A*Q to denote which columns are not matched. If the matrix is
* structurally rank deficient, some of the entries on the diagonal of the
* permuted matrix will be zero. BTF_ORDER then calls btf_strongcomp to find
* the strongly-connected components.
*
* On output, P and Q are the row and column permutations, where i = P[k] if
* row i of A is the kth row of P*A*Q, and j = BTF_UNFLIP(Q[k]) if column j of
* A is the kth column of P*A*Q. If Q[k] < 0, then the (k,k)th entry in P*A*Q
* is structurally zero.
*
* The vector R gives the block boundaries, where block b is in rows/columns
* R[b] to R[b+1]-1 of the permuted matrix, and where b ranges from 1 to the
* number of strongly connected components found.
*/
int btf_order /* returns number of blocks found */
(
/* --- input, not modified: --- */
int n, /* A is n-by-n in compressed column form */
int Ap [ ], /* size n+1 */
int Ai [ ], /* size nz = Ap [n] */
double maxwork, /* do at most maxwork*nnz(A) work in the maximum
* transversal; no limit if <= 0 */
/* --- output, not defined on input --- */
double *work, /* return value from btf_maxtrans */
int P [ ], /* size n, row permutation */
int Q [ ], /* size n, column permutation */
int R [ ], /* size n+1. block b is in rows/cols R[b] ... R[b+1]-1 */
int *nmatch, /* # nonzeros on diagonal of P*A*Q */
/* --- workspace, not defined on input or output --- */
int Work [ ] /* size 5n */
) ;
UF_long btf_l_order (UF_long, UF_long *, UF_long *, double , double *,
UF_long *, UF_long *, UF_long *, UF_long *, UF_long *) ;
/* ========================================================================== */
/* === BTF marking of singular columns ====================================== */
/* ========================================================================== */
/* BTF_FLIP is a "negation about -1", and is used to mark an integer j
* that is normally non-negative. BTF_FLIP (-1) is -1. BTF_FLIP of
* a number > -1 is negative, and BTF_FLIP of a number < -1 is positive.
* BTF_FLIP (BTF_FLIP (j)) = j for all integers j. UNFLIP (j) acts
* like an "absolute value" operation, and is always >= -1. You can test
* whether or not an integer j is "flipped" with the BTF_ISFLIPPED (j)
* macro.
*/
#define BTF_FLIP(j) (-(j)-2)
#define BTF_ISFLIPPED(j) ((j) < -1)
#define BTF_UNFLIP(j) ((BTF_ISFLIPPED (j)) ? BTF_FLIP (j) : (j))
/* ========================================================================== */
/* === BTF version ========================================================== */
/* ========================================================================== */
/* All versions of BTF include these definitions.
* As an example, to test if the version you are using is 1.2 or later:
*
* if (BTF_VERSION >= BTF_VERSION_CODE (1,2)) ...
*
* This also works during compile-time:
*
* #if (BTF >= BTF_VERSION_CODE (1,2))
* printf ("This is version 1.2 or later\n") ;
* #else
* printf ("This is an early version\n") ;
* #endif
*/
#define BTF_DATE "Dec 7, 2011"
#define BTF_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
#define BTF_MAIN_VERSION 1
#define BTF_SUB_VERSION 1
#define BTF_SUBSUB_VERSION 3
#define BTF_VERSION BTF_VERSION_CODE(BTF_MAIN_VERSION,BTF_SUB_VERSION)
#ifdef __cplusplus
}
#endif
#endif

255
src/include/ngspice/colamd.h Normal file
View File

@ -0,0 +1,255 @@
/* ========================================================================== */
/* === colamd/symamd prototypes and definitions ============================= */
/* ========================================================================== */
/* COLAMD / SYMAMD include file
You must include this file (colamd.h) in any routine that uses colamd,
symamd, or the related macros and definitions.
Authors:
The authors of the code itself are Stefan I. Larimore and Timothy A.
Davis (davis at cise.ufl.edu), University of Florida. The algorithm was
developed in collaboration with John Gilbert, Xerox PARC, and Esmond
Ng, Oak Ridge National Laboratory.
Acknowledgements:
This work was supported by the National Science Foundation, under
grants DMS-9504974 and DMS-9803599.
Notice:
Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
Permission is hereby granted to use, copy, modify, and/or distribute
this program, provided that the Copyright, this License, and the
Availability of the original version is retained on all copies and made
accessible to the end-user of any code or package that includes COLAMD
or any modified version of COLAMD.
Availability:
The colamd/symamd library is available at
http://www.cise.ufl.edu/research/sparse/colamd/
This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.h
file. It is required by the colamd.c, colamdmex.c, and symamdmex.c
files, and by any C code that calls the routines whose prototypes are
listed below, or that uses the colamd/symamd definitions listed below.
*/
#ifndef COLAMD_H
#define COLAMD_H
/* make it easy for C++ programs to include COLAMD */
#ifdef __cplusplus
extern "C" {
#endif
/* ========================================================================== */
/* === Include files ======================================================== */
/* ========================================================================== */
#include <stdlib.h>
/* ========================================================================== */
/* === COLAMD version ======================================================= */
/* ========================================================================== */
/* COLAMD Version 2.4 and later will include the following definitions.
* As an example, to test if the version you are using is 2.4 or later:
*
* #ifdef COLAMD_VERSION
* if (COLAMD_VERSION >= COLAMD_VERSION_CODE (2,4)) ...
* #endif
*
* This also works during compile-time:
*
* #if defined(COLAMD_VERSION) && (COLAMD_VERSION >= COLAMD_VERSION_CODE (2,4))
* printf ("This is version 2.4 or later\n") ;
* #else
* printf ("This is an early version\n") ;
* #endif
*
* Versions 2.3 and earlier of COLAMD do not include a #define'd version number.
*/
#define COLAMD_DATE "Dec 7, 2011"
#define COLAMD_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
#define COLAMD_MAIN_VERSION 2
#define COLAMD_SUB_VERSION 7
#define COLAMD_SUBSUB_VERSION 4
#define COLAMD_VERSION \
COLAMD_VERSION_CODE(COLAMD_MAIN_VERSION,COLAMD_SUB_VERSION)
/* ========================================================================== */
/* === Knob and statistics definitions ====================================== */
/* ========================================================================== */
/* size of the knobs [ ] array. Only knobs [0..1] are currently used. */
#define COLAMD_KNOBS 20
/* number of output statistics. Only stats [0..6] are currently used. */
#define COLAMD_STATS 20
/* knobs [0] and stats [0]: dense row knob and output statistic. */
#define COLAMD_DENSE_ROW 0
/* knobs [1] and stats [1]: dense column knob and output statistic. */
#define COLAMD_DENSE_COL 1
/* knobs [2]: aggressive absorption */
#define COLAMD_AGGRESSIVE 2
/* stats [2]: memory defragmentation count output statistic */
#define COLAMD_DEFRAG_COUNT 2
/* stats [3]: colamd status: zero OK, > 0 warning or notice, < 0 error */
#define COLAMD_STATUS 3
/* stats [4..6]: error info, or info on jumbled columns */
#define COLAMD_INFO1 4
#define COLAMD_INFO2 5
#define COLAMD_INFO3 6
/* error codes returned in stats [3]: */
#define COLAMD_OK (0)
#define COLAMD_OK_BUT_JUMBLED (1)
#define COLAMD_ERROR_A_not_present (-1)
#define COLAMD_ERROR_p_not_present (-2)
#define COLAMD_ERROR_nrow_negative (-3)
#define COLAMD_ERROR_ncol_negative (-4)
#define COLAMD_ERROR_nnz_negative (-5)
#define COLAMD_ERROR_p0_nonzero (-6)
#define COLAMD_ERROR_A_too_small (-7)
#define COLAMD_ERROR_col_length_negative (-8)
#define COLAMD_ERROR_row_index_out_of_bounds (-9)
#define COLAMD_ERROR_out_of_memory (-10)
#define COLAMD_ERROR_internal_error (-999)
/* ========================================================================== */
/* === Prototypes of user-callable routines ================================= */
/* ========================================================================== */
/* define UF_long */
#include "UFconfig.h"
size_t colamd_recommended /* returns recommended value of Alen, */
/* or 0 if input arguments are erroneous */
(
int nnz, /* nonzeros in A */
int n_row, /* number of rows in A */
int n_col /* number of columns in A */
) ;
size_t colamd_l_recommended /* returns recommended value of Alen, */
/* or 0 if input arguments are erroneous */
(
UF_long nnz, /* nonzeros in A */
UF_long n_row, /* number of rows in A */
UF_long n_col /* number of columns in A */
) ;
void colamd_set_defaults /* sets default parameters */
( /* knobs argument is modified on output */
double knobs [COLAMD_KNOBS] /* parameter settings for colamd */
) ;
void colamd_l_set_defaults /* sets default parameters */
( /* knobs argument is modified on output */
double knobs [COLAMD_KNOBS] /* parameter settings for colamd */
) ;
int colamd /* returns (1) if successful, (0) otherwise*/
( /* A and p arguments are modified on output */
int n_row, /* number of rows in A */
int n_col, /* number of columns in A */
int Alen, /* size of the array A */
int A [], /* row indices of A, of size Alen */
int p [], /* column pointers of A, of size n_col+1 */
double knobs [COLAMD_KNOBS],/* parameter settings for colamd */
int stats [COLAMD_STATS] /* colamd output statistics and error codes */
) ;
UF_long colamd_l /* returns (1) if successful, (0) otherwise*/
( /* A and p arguments are modified on output */
UF_long n_row, /* number of rows in A */
UF_long n_col, /* number of columns in A */
UF_long Alen, /* size of the array A */
UF_long A [], /* row indices of A, of size Alen */
UF_long p [], /* column pointers of A, of size n_col+1 */
double knobs [COLAMD_KNOBS],/* parameter settings for colamd */
UF_long stats [COLAMD_STATS]/* colamd output statistics and error codes */
) ;
int symamd /* return (1) if OK, (0) otherwise */
(
int n, /* number of rows and columns of A */
int A [], /* row indices of A */
int p [], /* column pointers of A */
int perm [], /* output permutation, size n_col+1 */
double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */
int stats [COLAMD_STATS], /* output statistics and error codes */
void * (*allocate) (size_t, size_t),
/* pointer to calloc (ANSI C) or */
/* mxCalloc (for MATLAB mexFunction) */
void (*release) (void *)
/* pointer to free (ANSI C) or */
/* mxFree (for MATLAB mexFunction) */
) ;
UF_long symamd_l /* return (1) if OK, (0) otherwise */
(
UF_long n, /* number of rows and columns of A */
UF_long A [], /* row indices of A */
UF_long p [], /* column pointers of A */
UF_long perm [], /* output permutation, size n_col+1 */
double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */
UF_long stats [COLAMD_STATS], /* output statistics and error codes */
void * (*allocate) (size_t, size_t),
/* pointer to calloc (ANSI C) or */
/* mxCalloc (for MATLAB mexFunction) */
void (*release) (void *)
/* pointer to free (ANSI C) or */
/* mxFree (for MATLAB mexFunction) */
) ;
void colamd_report
(
int stats [COLAMD_STATS]
) ;
void colamd_l_report
(
UF_long stats [COLAMD_STATS]
) ;
void symamd_report
(
int stats [COLAMD_STATS]
) ;
void symamd_l_report
(
UF_long stats [COLAMD_STATS]
) ;
#ifndef EXTERN
#define EXTERN extern
#endif
EXTERN int (*colamd_printf) (const char *, ...) ;
#ifdef __cplusplus
}
#endif
#endif /* COLAMD_H */

831
src/include/ngspice/klu.h Normal file
View File

@ -0,0 +1,831 @@
/* ========================================================================== */
/* === klu include file ===================================================== */
/* ========================================================================== */
/* Include file for user programs that call klu_* routines */
#ifndef _KLU_H
#define _KLU_H
/* make it easy for C++ programs to include KLU */
#ifdef __cplusplus
extern "C" {
#endif
#include "amd.h"
#include "colamd.h"
#include "btf.h"
/* -------------------------------------------------------------------------- */
/* Symbolic object - contains the pre-ordering computed by klu_analyze */
/* -------------------------------------------------------------------------- */
typedef struct
{
/* A (P,Q) is in upper block triangular form. The kth block goes from
* row/col index R [k] to R [k+1]-1. The estimated number of nonzeros
* in the L factor of the kth block is Lnz [k].
*/
/* only computed if the AMD ordering is chosen: */
double symmetry ; /* symmetry of largest block */
double est_flops ; /* est. factorization flop count */
double lnz, unz ; /* estimated nz in L and U, including diagonals */
double *Lnz ; /* size n, but only Lnz [0..nblocks-1] is used */
/* computed for all orderings: */
int
n, /* input matrix A is n-by-n */
nz, /* # entries in input matrix */
*P, /* size n */
*Q, /* size n */
*R, /* size n+1, but only R [0..nblocks] is used */
nzoff, /* nz in off-diagonal blocks */
nblocks, /* number of blocks */
maxblock, /* size of largest block */
ordering, /* ordering used (AMD, COLAMD, or GIVEN) */
do_btf ; /* whether or not BTF preordering was requested */
/* only computed if BTF preordering requested */
int structural_rank ; /* 0 to n-1 if the matrix is structurally rank
* deficient. -1 if not computed. n if the matrix has
* full structural rank */
} klu_symbolic ;
typedef struct /* 64-bit version (otherwise same as above) */
{
double symmetry, est_flops, lnz, unz ;
double *Lnz ;
UF_long n, nz, *P, *Q, *R, nzoff, nblocks, maxblock, ordering, do_btf,
structural_rank ;
} klu_l_symbolic ;
/* -------------------------------------------------------------------------- */
/* Numeric object - contains the factors computed by klu_factor */
/* -------------------------------------------------------------------------- */
typedef struct
{
/* LU factors of each block, the pivot row permutation, and the
* entries in the off-diagonal blocks */
int n ; /* A is n-by-n */
int nblocks ; /* number of diagonal blocks */
int lnz ; /* actual nz in L, including diagonal */
int unz ; /* actual nz in U, including diagonal */
int max_lnz_block ; /* max actual nz in L in any one block, incl. diag */
int max_unz_block ; /* max actual nz in U in any one block, incl. diag */
int *Pnum ; /* size n. final pivot permutation */
int *Pinv ; /* size n. inverse of final pivot permutation */
/* LU factors of each block */
int *Lip ; /* size n. pointers into LUbx[block] for L */
int *Uip ; /* size n. pointers into LUbx[block] for U */
int *Llen ; /* size n. Llen [k] = # of entries in kth column of L */
int *Ulen ; /* size n. Ulen [k] = # of entries in kth column of U */
void **LUbx ; /* L and U indices and entries (excl. diagonal of U) */
size_t *LUsize ; /* size of each LUbx [block], in sizeof (Unit) */
void *Udiag ; /* diagonal of U */
/* scale factors; can be NULL if no scaling */
double *Rs ; /* size n. Rs [i] is scale factor for row i */
/* permanent workspace for factorization and solve */
size_t worksize ; /* size (in bytes) of Work */
void *Work ; /* workspace */
void *Xwork ; /* alias into Numeric->Work */
int *Iwork ; /* alias into Numeric->Work */
/* off-diagonal entries in a conventional compressed-column sparse matrix */
int *Offp ; /* size n+1, column pointers */
int *Offi ; /* size nzoff, row indices */
void *Offx ; /* size nzoff, numerical values */
int nzoff ;
} klu_numeric ;
typedef struct /* 64-bit version (otherwise same as above) */
{
UF_long n, nblocks, lnz, unz, max_lnz_block, max_unz_block, *Pnum, *Pinv,
*Lip, *Uip, *Llen, *Ulen ;
void **LUbx ;
size_t *LUsize ;
void *Udiag ;
double *Rs ;
size_t worksize ;
void *Work, *Xwork ;
UF_long *Iwork ;
UF_long *Offp, *Offi ;
void *Offx ;
UF_long nzoff ;
} klu_l_numeric ;
/* -------------------------------------------------------------------------- */
/* KLU control parameters and statistics */
/* -------------------------------------------------------------------------- */
/* Common->status values */
#define KLU_OK 0
#define KLU_SINGULAR (1) /* status > 0 is a warning, not an error */
#define KLU_OUT_OF_MEMORY (-2)
#define KLU_INVALID (-3)
#define KLU_TOO_LARGE (-4) /* integer overflow has occured */
typedef struct klu_common_struct
{
/* ---------------------------------------------------------------------- */
/* parameters */
/* ---------------------------------------------------------------------- */
double tol ; /* pivot tolerance for diagonal preference */
double memgrow ; /* realloc memory growth size for LU factors */
double initmem_amd ; /* init. memory size with AMD: c*nnz(L) + n */
double initmem ; /* init. memory size: c*nnz(A) + n */
double maxwork ; /* maxwork for BTF, <= 0 if no limit */
int btf ; /* use BTF pre-ordering, or not */
int ordering ; /* 0: AMD, 1: COLAMD, 2: user P and Q,
* 3: user function */
int scale ; /* row scaling: -1: none (and no error check),
* 0: none, 1: sum, 2: max */
/* memory management routines */
void *(*malloc_memory) (size_t) ; /* pointer to malloc */
void *(*realloc_memory) (void *, size_t) ; /* pointer to realloc */
void (*free_memory) (void *) ; /* pointer to free */
void *(*calloc_memory) (size_t, size_t) ; /* pointer to calloc */
/* pointer to user ordering function */
int (*user_order) (int, int *, int *, int *, struct klu_common_struct *) ;
/* pointer to user data, passed unchanged as the last parameter to the
* user ordering function (optional, the user function need not use this
* information). */
void *user_data ;
int halt_if_singular ; /* how to handle a singular matrix:
* FALSE: keep going. Return a Numeric object with a zero U(k,k). A
* divide-by-zero may occur when computing L(:,k). The Numeric object
* can be passed to klu_solve (a divide-by-zero will occur). It can
* also be safely passed to klu_refactor.
* TRUE: stop quickly. klu_factor will free the partially-constructed
* Numeric object. klu_refactor will not free it, but will leave the
* numerical values only partially defined. This is the default. */
/* ---------------------------------------------------------------------- */
/* statistics */
/* ---------------------------------------------------------------------- */
int status ; /* KLU_OK if OK, < 0 if error */
int nrealloc ; /* # of reallocations of L and U */
int structural_rank ; /* 0 to n-1 if the matrix is structurally rank
* deficient (as determined by maxtrans). -1 if not computed. n if the
* matrix has full structural rank. This is computed by klu_analyze
* if a BTF preordering is requested. */
int numerical_rank ; /* First k for which a zero U(k,k) was found,
* if the matrix was singular (in the range 0 to n-1). n if the matrix
* has full rank. This is not a true rank-estimation. It just reports
* where the first zero pivot was found. -1 if not computed.
* Computed by klu_factor and klu_refactor. */
int singular_col ; /* n if the matrix is not singular. If in the
* range 0 to n-1, this is the column index of the original matrix A that
* corresponds to the column of U that contains a zero diagonal entry.
* -1 if not computed. Computed by klu_factor and klu_refactor. */
int noffdiag ; /* # of off-diagonal pivots, -1 if not computed */
double flops ; /* actual factorization flop count, from klu_flops */
double rcond ; /* crude reciprocal condition est., from klu_rcond */
double condest ; /* accurate condition est., from klu_condest */
double rgrowth ; /* reciprocal pivot rgrowth, from klu_rgrowth */
double work ; /* actual work done in BTF, in klu_analyze */
size_t memusage ; /* current memory usage, in bytes */
size_t mempeak ; /* peak memory usage, in bytes */
} klu_common ;
typedef struct klu_l_common_struct /* 64-bit version (otherwise same as above)*/
{
double tol, memgrow, initmem_amd, initmem, maxwork ;
UF_long btf, ordering, scale ;
void *(*malloc_memory) (size_t) ;
void *(*realloc_memory) (void *, size_t) ;
void (*free_memory) (void *) ;
void *(*calloc_memory) (size_t, size_t) ;
UF_long (*user_order) (UF_long, UF_long *, UF_long *, UF_long *,
struct klu_l_common_struct *) ;
void *user_data ;
UF_long halt_if_singular ;
UF_long status, nrealloc, structural_rank, numerical_rank, singular_col,
noffdiag ;
double flops, rcond, condest, rgrowth, work ;
size_t memusage, mempeak ;
} klu_l_common ;
/* -------------------------------------------------------------------------- */
/* klu_defaults: sets default control parameters */
/* -------------------------------------------------------------------------- */
int klu_defaults
(
klu_common *Common
) ;
UF_long klu_l_defaults (klu_l_common *Common) ;
/* -------------------------------------------------------------------------- */
/* klu_analyze: orders and analyzes a matrix */
/* -------------------------------------------------------------------------- */
/* Order the matrix with BTF (or not), then order each block with AMD, COLAMD,
* a natural ordering, or with a user-provided ordering function */
klu_symbolic *klu_analyze
(
/* inputs, not modified */
int n, /* A is n-by-n */
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
klu_common *Common
) ;
klu_l_symbolic *klu_l_analyze (UF_long, UF_long *, UF_long *,
klu_l_common *Common) ;
/* -------------------------------------------------------------------------- */
/* klu_analyze_given: analyzes a matrix using given P and Q */
/* -------------------------------------------------------------------------- */
/* Order the matrix with BTF (or not), then use natural or given ordering
* P and Q on the blocks. P and Q are interpretted as identity
* if NULL. */
klu_symbolic *klu_analyze_given
(
/* inputs, not modified */
int n, /* A is n-by-n */
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
int P [ ], /* size n, user's row permutation (may be NULL) */
int Q [ ], /* size n, user's column permutation (may be NULL) */
klu_common *Common
) ;
klu_l_symbolic *klu_l_analyze_given (UF_long, UF_long *, UF_long *, UF_long *,
UF_long *, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_factor: factors a matrix using the klu_analyze results */
/* -------------------------------------------------------------------------- */
klu_numeric *klu_factor /* returns KLU_OK if OK, < 0 if error */
(
/* inputs, not modified */
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
double Ax [ ], /* size nz, numerical values */
klu_symbolic *Symbolic,
klu_common *Common
) ;
klu_numeric *klu_z_factor /* returns KLU_OK if OK, < 0 if error */
(
/* inputs, not modified */
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
double Ax [ ], /* size 2*nz, numerical values (real,imag pairs) */
klu_symbolic *Symbolic,
klu_common *Common
) ;
/* long / real version */
klu_l_numeric *klu_l_factor (UF_long *, UF_long *, double *, klu_l_symbolic *,
klu_l_common *) ;
/* long / complex version */
klu_l_numeric *klu_zl_factor (UF_long *, UF_long *, double *, klu_l_symbolic *,
klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_solve: solves Ax=b using the Symbolic and Numeric objects */
/* -------------------------------------------------------------------------- */
int klu_solve
(
/* inputs, not modified */
klu_symbolic *Symbolic,
klu_numeric *Numeric,
int ldim, /* leading dimension of B */
int nrhs, /* number of right-hand-sides */
/* right-hand-side on input, overwritten with solution to Ax=b on output */
double B [ ], /* size ldim*nrhs */
klu_common *Common
) ;
int klu_z_solve
(
/* inputs, not modified */
klu_symbolic *Symbolic,
klu_numeric *Numeric,
int ldim, /* leading dimension of B */
int nrhs, /* number of right-hand-sides */
/* right-hand-side on input, overwritten with solution to Ax=b on output */
double B [ ], /* size 2*ldim*nrhs */
klu_common *Common
) ;
UF_long klu_l_solve (klu_l_symbolic *, klu_l_numeric *, UF_long, UF_long,
double *, klu_l_common *) ;
UF_long klu_zl_solve (klu_l_symbolic *, klu_l_numeric *, UF_long, UF_long,
double *, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_tsolve: solves A'x=b using the Symbolic and Numeric objects */
/* -------------------------------------------------------------------------- */
int klu_tsolve
(
/* inputs, not modified */
klu_symbolic *Symbolic,
klu_numeric *Numeric,
int ldim, /* leading dimension of B */
int nrhs, /* number of right-hand-sides */
/* right-hand-side on input, overwritten with solution to Ax=b on output */
double B [ ], /* size ldim*nrhs */
klu_common *Common
) ;
int klu_z_tsolve
(
/* inputs, not modified */
klu_symbolic *Symbolic,
klu_numeric *Numeric,
int ldim, /* leading dimension of B */
int nrhs, /* number of right-hand-sides */
/* right-hand-side on input, overwritten with solution to Ax=b on output */
double B [ ], /* size 2*ldim*nrhs */
int conj_solve, /* TRUE: conjugate solve, FALSE: solve A.'x=b */
klu_common *Common
) ;
UF_long klu_l_tsolve (klu_l_symbolic *, klu_l_numeric *, UF_long, UF_long,
double *, klu_l_common *) ;
UF_long klu_zl_tsolve (klu_l_symbolic *, klu_l_numeric *, UF_long, UF_long,
double *, UF_long, klu_l_common * ) ;
/* -------------------------------------------------------------------------- */
/* klu_refactor: refactorizes matrix with same ordering as klu_factor */
/* -------------------------------------------------------------------------- */
int klu_refactor /* return TRUE if successful, FALSE otherwise */
(
/* inputs, not modified */
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
double Ax [ ], /* size nz, numerical values */
klu_symbolic *Symbolic,
/* input, and numerical values modified on output */
klu_numeric *Numeric,
klu_common *Common
) ;
int klu_z_refactor /* return TRUE if successful, FALSE otherwise */
(
/* inputs, not modified */
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
double Ax [ ], /* size 2*nz, numerical values */
klu_symbolic *Symbolic,
/* input, and numerical values modified on output */
klu_numeric *Numeric,
klu_common *Common
) ;
UF_long klu_l_refactor (UF_long *, UF_long *, double *, klu_l_symbolic *,
klu_l_numeric *, klu_l_common *) ;
UF_long klu_zl_refactor (UF_long *, UF_long *, double *, klu_l_symbolic *,
klu_l_numeric *, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_free_symbolic: destroys the Symbolic object */
/* -------------------------------------------------------------------------- */
int klu_free_symbolic
(
klu_symbolic **Symbolic,
klu_common *Common
) ;
UF_long klu_l_free_symbolic (klu_l_symbolic **, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_free_numeric: destroys the Numeric object */
/* -------------------------------------------------------------------------- */
/* Note that klu_free_numeric and klu_z_free_numeric are identical; each can
* free both kinds of Numeric objects (real and complex) */
int klu_free_numeric
(
klu_numeric **Numeric,
klu_common *Common
) ;
int klu_z_free_numeric
(
klu_numeric **Numeric,
klu_common *Common
) ;
UF_long klu_l_free_numeric (klu_l_numeric **, klu_l_common *) ;
UF_long klu_zl_free_numeric (klu_l_numeric **, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_sort: sorts the columns of the LU factorization */
/* -------------------------------------------------------------------------- */
/* this is not needed except for the MATLAB interface */
int klu_sort
(
/* inputs, not modified */
klu_symbolic *Symbolic,
/* input/output */
klu_numeric *Numeric,
klu_common *Common
) ;
int klu_z_sort
(
/* inputs, not modified */
klu_symbolic *Symbolic,
/* input/output */
klu_numeric *Numeric,
klu_common *Common
) ;
UF_long klu_l_sort (klu_l_symbolic *, klu_l_numeric *, klu_l_common *) ;
UF_long klu_zl_sort (klu_l_symbolic *, klu_l_numeric *, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_flops: determines # of flops performed in numeric factorzation */
/* -------------------------------------------------------------------------- */
int klu_flops
(
/* inputs, not modified */
klu_symbolic *Symbolic,
klu_numeric *Numeric,
/* input/output */
klu_common *Common
) ;
int klu_z_flops
(
/* inputs, not modified */
klu_symbolic *Symbolic,
klu_numeric *Numeric,
/* input/output */
klu_common *Common
) ;
UF_long klu_l_flops (klu_l_symbolic *, klu_l_numeric *, klu_l_common *) ;
UF_long klu_zl_flops (klu_l_symbolic *, klu_l_numeric *, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_rgrowth : compute the reciprocal pivot growth */
/* -------------------------------------------------------------------------- */
/* Pivot growth is computed after the input matrix is permuted, scaled, and
* off-diagonal entries pruned. This is because the LU factorization of each
* block takes as input the scaled diagonal blocks of the BTF form. The
* reciprocal pivot growth in column j of an LU factorization of a matrix C
* is the largest entry in C divided by the largest entry in U; then the overall
* reciprocal pivot growth is the smallest such value for all columns j. Note
* that the off-diagonal entries are not scaled, since they do not take part in
* the LU factorization of the diagonal blocks.
*
* In MATLAB notation:
*
* rgrowth = min (max (abs ((R \ A(p,q)) - F)) ./ max (abs (U))) */
int klu_rgrowth
(
int Ap [ ],
int Ai [ ],
double Ax [ ],
klu_symbolic *Symbolic,
klu_numeric *Numeric,
klu_common *Common /* Common->rgrowth = reciprocal pivot growth */
) ;
int klu_z_rgrowth
(
int Ap [ ],
int Ai [ ],
double Ax [ ],
klu_symbolic *Symbolic,
klu_numeric *Numeric,
klu_common *Common /* Common->rgrowth = reciprocal pivot growth */
) ;
UF_long klu_l_rgrowth (UF_long *, UF_long *, double *, klu_l_symbolic *,
klu_l_numeric *, klu_l_common *) ;
UF_long klu_zl_rgrowth (UF_long *, UF_long *, double *, klu_l_symbolic *,
klu_l_numeric *, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_condest */
/* -------------------------------------------------------------------------- */
/* Computes a reasonably accurate estimate of the 1-norm condition number, using
* Hager's method, as modified by Higham and Tisseur (same method as used in
* MATLAB's condest */
int klu_condest
(
int Ap [ ], /* size n+1, column pointers, not modified */
double Ax [ ], /* size nz = Ap[n], numerical values, not modified*/
klu_symbolic *Symbolic, /* symbolic analysis, not modified */
klu_numeric *Numeric, /* numeric factorization, not modified */
klu_common *Common /* result returned in Common->condest */
) ;
int klu_z_condest
(
int Ap [ ],
double Ax [ ], /* size 2*nz */
klu_symbolic *Symbolic,
klu_numeric *Numeric,
klu_common *Common /* result returned in Common->condest */
) ;
UF_long klu_l_condest (UF_long *, double *, klu_l_symbolic *, klu_l_numeric *,
klu_l_common *) ;
UF_long klu_zl_condest (UF_long *, double *, klu_l_symbolic *, klu_l_numeric *,
klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_rcond: compute min(abs(diag(U))) / max(abs(diag(U))) */
/* -------------------------------------------------------------------------- */
int klu_rcond
(
klu_symbolic *Symbolic, /* input, not modified */
klu_numeric *Numeric, /* input, not modified */
klu_common *Common /* result in Common->rcond */
) ;
int klu_z_rcond
(
klu_symbolic *Symbolic, /* input, not modified */
klu_numeric *Numeric, /* input, not modified */
klu_common *Common /* result in Common->rcond */
) ;
UF_long klu_l_rcond (klu_l_symbolic *, klu_l_numeric *, klu_l_common *) ;
UF_long klu_zl_rcond (klu_l_symbolic *, klu_l_numeric *, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_scale */
/* -------------------------------------------------------------------------- */
int klu_scale /* return TRUE if successful, FALSE otherwise */
(
/* inputs, not modified */
int scale, /* <0: none, no error check; 0: none, 1: sum, 2: max */
int n,
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
double Ax [ ],
/* outputs, not defined on input */
double Rs [ ],
/* workspace, not defined on input or output */
int W [ ], /* size n, can be NULL */
klu_common *Common
) ;
int klu_z_scale /* return TRUE if successful, FALSE otherwise */
(
/* inputs, not modified */
int scale, /* <0: none, no error check; 0: none, 1: sum, 2: max */
int n,
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
double Ax [ ],
/* outputs, not defined on input */
double Rs [ ],
/* workspace, not defined on input or output */
int W [ ], /* size n, can be NULL */
klu_common *Common
) ;
UF_long klu_l_scale (UF_long, UF_long, UF_long *, UF_long *, double *,
double *, UF_long *, klu_l_common *) ;
UF_long klu_zl_scale (UF_long, UF_long, UF_long *, UF_long *, double *,
double *, UF_long *, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* klu_extract */
/* -------------------------------------------------------------------------- */
int klu_extract /* returns TRUE if successful, FALSE otherwise */
(
/* inputs: */
klu_numeric *Numeric,
klu_symbolic *Symbolic,
/* outputs, either allocated on input, or ignored otherwise */
/* L */
int *Lp, /* size n+1 */
int *Li, /* size Numeric->lnz */
double *Lx, /* size Numeric->lnz */
/* U */
int *Up, /* size n+1 */
int *Ui, /* size Numeric->unz */
double *Ux, /* size Numeric->unz */
/* F */
int *Fp, /* size n+1 */
int *Fi, /* size Numeric->nzoff */
double *Fx, /* size Numeric->nzoff */
/* P, row permutation */
int *P, /* size n */
/* Q, column permutation */
int *Q, /* size n */
/* Rs, scale factors */
double *Rs, /* size n */
/* R, block boundaries */
int *R, /* size Symbolic->nblocks+1 (nblocks is at most n) */
klu_common *Common
) ;
int klu_z_extract /* returns TRUE if successful, FALSE otherwise */
(
/* inputs: */
klu_numeric *Numeric,
klu_symbolic *Symbolic,
/* outputs, all of which must be allocated on input */
/* L */
int *Lp, /* size n+1 */
int *Li, /* size nnz(L) */
double *Lx, /* size nnz(L) */
double *Lz, /* size nnz(L) for the complex case, ignored if real */
/* U */
int *Up, /* size n+1 */
int *Ui, /* size nnz(U) */
double *Ux, /* size nnz(U) */
double *Uz, /* size nnz(U) for the complex case, ignored if real */
/* F */
int *Fp, /* size n+1 */
int *Fi, /* size nnz(F) */
double *Fx, /* size nnz(F) */
double *Fz, /* size nnz(F) for the complex case, ignored if real */
/* P, row permutation */
int *P, /* size n */
/* Q, column permutation */
int *Q, /* size n */
/* Rs, scale factors */
double *Rs, /* size n */
/* R, block boundaries */
int *R, /* size Symbolic->nblocks+1 (nblocks is at most n) */
klu_common *Common
) ;
UF_long klu_l_extract (klu_l_numeric *, klu_l_symbolic *,
UF_long *, UF_long *, double *,
UF_long *, UF_long *, double *,
UF_long *, UF_long *, double *,
UF_long *, UF_long *, double *, UF_long *, klu_l_common *) ;
UF_long klu_zl_extract (klu_l_numeric *, klu_l_symbolic *,
UF_long *, UF_long *, double *, double *,
UF_long *, UF_long *, double *, double *,
UF_long *, UF_long *, double *, double *,
UF_long *, UF_long *, double *, UF_long *, klu_l_common *) ;
/* -------------------------------------------------------------------------- */
/* KLU memory management routines */
/* -------------------------------------------------------------------------- */
void *klu_malloc /* returns pointer to the newly malloc'd block */
(
/* ---- input ---- */
size_t n, /* number of items */
size_t size, /* size of each item */
/* --------------- */
klu_common *Common
) ;
void *klu_free /* always returns NULL */
(
/* ---- in/out --- */
void *p, /* block of memory to free */
size_t n, /* number of items */
size_t size, /* size of each item */
/* --------------- */
klu_common *Common
) ;
void *klu_realloc /* returns pointer to reallocated block */
(
/* ---- input ---- */
size_t nnew, /* requested # of items in reallocated block */
size_t nold, /* current size of block, in # of items */
size_t size, /* size of each item */
/* ---- in/out --- */
void *p, /* block of memory to realloc */
/* --------------- */
klu_common *Common
) ;
void *klu_l_malloc (size_t, size_t, klu_l_common *) ;
void *klu_l_free (void *, size_t, size_t, klu_l_common *) ;
void *klu_l_realloc (size_t, size_t, size_t, void *, klu_l_common *) ;
/* ========================================================================== */
/* === KLU version ========================================================== */
/* ========================================================================== */
/* All versions of KLU include these definitions.
* As an example, to test if the version you are using is 1.2 or later:
*
* if (KLU_VERSION >= KLU_VERSION_CODE (1,2)) ...
*
* This also works during compile-time:
*
* #if (KLU >= KLU_VERSION_CODE (1,2))
* printf ("This is version 1.2 or later\n") ;
* #else
* printf ("This is an early version\n") ;
* #endif
*/
#define KLU_DATE "Nov 30, 2009"
#define KLU_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
#define KLU_MAIN_VERSION 1
#define KLU_SUB_VERSION 1
#define KLU_SUBSUB_VERSION 1
#define KLU_VERSION KLU_VERSION_CODE(KLU_MAIN_VERSION,KLU_SUB_VERSION)
#ifdef __cplusplus
}
#endif
#endif

145
src/maths/sparse/spCSC.c Normal file
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) ;
}