diff --git a/src/include/ngspice/UFconfig.h b/src/include/ngspice/UFconfig.h new file mode 100644 index 000000000..4b36a9e31 --- /dev/null +++ b/src/include/ngspice/UFconfig.h @@ -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 +#include + +/* ========================================================================== */ +/* === 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 diff --git a/src/include/ngspice/amd.h b/src/include/ngspice/amd.h new file mode 100644 index 000000000..2fcaef117 --- /dev/null +++ b/src/include/ngspice/amd.h @@ -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 + +/* 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 diff --git a/src/include/ngspice/btf.h b/src/include/ngspice/btf.h new file mode 100644 index 000000000..32be79cf9 --- /dev/null +++ b/src/include/ngspice/btf.h @@ -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 diff --git a/src/include/ngspice/colamd.h b/src/include/ngspice/colamd.h new file mode 100644 index 000000000..91665cef7 --- /dev/null +++ b/src/include/ngspice/colamd.h @@ -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 + +/* ========================================================================== */ +/* === 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 */ diff --git a/src/include/ngspice/klu.h b/src/include/ngspice/klu.h new file mode 100644 index 000000000..1f05df595 --- /dev/null +++ b/src/include/ngspice/klu.h @@ -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 diff --git a/src/maths/sparse/spCSC.c b/src/maths/sparse/spCSC.c new file mode 100644 index 000000000..2d7cd704e --- /dev/null +++ b/src/maths/sparse/spCSC.c @@ -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) ; +}