From 89cf3cf5c6033d1dd84b985a5b54c5cb77b21959 Mon Sep 17 00:00:00 2001 From: dwarning <> Date: Mon, 2 Dec 2024 14:23:02 +0100 Subject: [PATCH] use ngspice specific utility functions --- src/include/ngspice/klu.h | 73 ++++++++++++++ src/maths/KLU/btf_maxtrans.c | 26 ++--- src/maths/KLU/btf_version.c | 2 +- src/maths/KLU/klu.c | 8 +- src/maths/KLU/klu_extract.c | 178 ----------------------------------- src/maths/KLU/klu_utils.c | 178 +++++++++++++++++++++++++++++++++++ 6 files changed, 271 insertions(+), 194 deletions(-) diff --git a/src/include/ngspice/klu.h b/src/include/ngspice/klu.h index 0bb93c176..a28adef2c 100644 --- a/src/include/ngspice/klu.h +++ b/src/include/ngspice/klu.h @@ -834,6 +834,7 @@ 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 *) ; +/* Francesco - Utilities */ int klu_print ( int *Ap, @@ -854,6 +855,78 @@ int klu_z_print int *IntToExtColMap ) ; +int klu_constant_multiply +( + int *Ap, + double *Ax, + int n, + klu_common *Common, + double constant +) ; + +int klu_z_constant_multiply +( + int *Ap, + double *Ax, + int n, + klu_common *Common, + double constant +) ; + +int klu_matrix_vector_multiply +( + int *Ap, /* CSR */ + int *Ai, /* CSR */ + double *Ax, /* CSR */ + double *RHS, + double *Solution, + int *IntToExtRowMap, + int *IntToExtColMap, + int n, + klu_common *Common +) ; + +int klu_z_matrix_vector_multiply +( + int *Ap, /* CSR */ + int *Ai, /* CSR */ + double *Ax, /* CSR */ + double *RHS, + double *Solution, + double *iRHS, + double *iSolution, + int *IntToExtRowMap, + int *IntToExtColMap, + int n, + klu_common *Common +) ; + +int klu_convert_matrix_in_CSR +( + int *Ap_CSC, /* CSC */ + int *Ai_CSC, /* CSC */ + double *Ax_CSC, /* CSC */ + int *Ap_CSR, /* CSR */ + int *Ai_CSR, /* CSR */ + double *Ax_CSR, /* CSR */ + int n, + int nz, + klu_common *Common +) ; + +int klu_z_convert_matrix_in_CSR +( + int *Ap_CSC, /* CSC */ + int *Ai_CSC, /* CSC */ + double *Ax_CSC, /* CSC */ + int *Ap_CSR, /* CSR */ + int *Ai_CSR, /* CSR */ + double *Ax_CSR, /* CSR */ + int n, + int nz, + klu_common *Common +) ; + typedef struct sBindElement { double *COO ; double *CSC ; diff --git a/src/maths/KLU/btf_maxtrans.c b/src/maths/KLU/btf_maxtrans.c index 87034c3d4..13ebcd4f1 100644 --- a/src/maths/KLU/btf_maxtrans.c +++ b/src/maths/KLU/btf_maxtrans.c @@ -1,6 +1,12 @@ -/* ========================================================================== */ -/* === BTF_MAXTRANS ========================================================= */ -/* ========================================================================== */ +//------------------------------------------------------------------------------ +// BTF/Source/btf_maxtrans: maximum transversal +//------------------------------------------------------------------------------ + +// BTF, Copyright (c) 2004-2022, University of Florida. All Rights Reserved. +// Author: Timothy A. Davis. +// SPDX-License-Identifier: LGPL-2.1+ + +//------------------------------------------------------------------------------ /* Finds a column permutation that maximizes the number of entries on the * diagonal of a sparse matrix. See btf.h for more information. @@ -40,10 +46,7 @@ * Thus, for general usage, cs_maxtrans is preferred. For square matrices that * are typically structurally non-singular, maxtrans is preferred. A partial * maxtrans can still be very useful when solving a sparse linear system. - * - * Copyright (c) 2004-2007. Tim Davis, University of Florida, - * with support from Sandia National Laboratories. All Rights Reserved. - */ + */ #include "ngspice/btf.h" #include "btf_internal.h" @@ -59,7 +62,7 @@ * * * column k is not matched to any row * * entries in the path are nonzero - * * the pairs (i1,j1), (i2,j2), (i3,j3) ..., (is,js) have been + * * the pairs (i1,j1), (i2,j2), (i3,j3) ..., (is,js) have been * previously matched to each other * * (i(s+1), js) is nonzero, and row i(s+1) is not matched to any column * @@ -130,7 +133,7 @@ * for (p = head ; ...) DO 90 K=1,JORD */ -static Int augment +static int augment ( Int k, /* which stage of the main loop we're in */ Int Ap [ ], /* column pointers, size n+1 */ @@ -314,7 +317,8 @@ Int BTF(maxtrans) /* returns # of columns in the matching */ ) { Int *Cheap, *Flag, *Istack, *Jstack, *Pstack ; - Int i, j, k, nmatch, work_limit_reached, result ; + Int i, j, k, nmatch, work_limit_reached ; + int result ; /* ---------------------------------------------------------------------- */ /* get workspace and initialize */ @@ -332,7 +336,7 @@ Int BTF(maxtrans) /* returns # of columns in the matching */ for (j = 0 ; j < ncol ; j++) { Cheap [j] = Ap [j] ; - Flag [j] = EMPTY ; + Flag [j] = EMPTY ; } /* all rows and columns are currently unmatched */ diff --git a/src/maths/KLU/btf_version.c b/src/maths/KLU/btf_version.c index e4ce6ebcc..151f24d79 100644 --- a/src/maths/KLU/btf_version.c +++ b/src/maths/KLU/btf_version.c @@ -8,7 +8,7 @@ //------------------------------------------------------------------------------ -#include "ngspice/btf.h" +#include "btf.h" void btf_version (int version [3]) { diff --git a/src/maths/KLU/klu.c b/src/maths/KLU/klu.c index 0251ac0aa..7791bf7f0 100644 --- a/src/maths/KLU/klu.c +++ b/src/maths/KLU/klu.c @@ -116,11 +116,11 @@ size_t KLU_kernel_factor /* 0 if failure, size of LU if OK */ { Lsize = -Lsize ; Lsize = MAX (Lsize, 1.0) ; - lsize = (int) Lsize * anz + n ; + lsize = Lsize * anz + n ; } else { - lsize = (int) Lsize ; + lsize = Lsize ; } usize = lsize ; @@ -130,8 +130,8 @@ size_t KLU_kernel_factor /* 0 if failure, size of LU if OK */ maxlnz = (((double) n) * ((double) n) + ((double) n)) / 2. ; maxlnz = MIN (maxlnz, ((double) Int_MAX)) ; - lsize = MIN ((int) maxlnz, lsize) ; - usize = MIN ((int) maxlnz, usize) ; + lsize = MIN (maxlnz, lsize) ; + usize = MIN (maxlnz, usize) ; PRINTF (("Welcome to klu: n %d anz %d k1 %d lsize %d usize %d maxlnz %g\n", n, anz, k1, lsize, usize, maxlnz)) ; diff --git a/src/maths/KLU/klu_extract.c b/src/maths/KLU/klu_extract.c index 57798312f..d393f18cd 100644 --- a/src/maths/KLU/klu_extract.c +++ b/src/maths/KLU/klu_extract.c @@ -294,181 +294,3 @@ int KLU_extract /* returns TRUE if successful, FALSE otherwise */ return (TRUE) ; } - -/* Francesco - Extract only Udiag */ -Int KLU_extract_Udiag /* returns TRUE if successful, FALSE otherwise */ -( - /* inputs: */ - KLU_numeric *Numeric, - KLU_symbolic *Symbolic, - - /* outputs, all of which must be allocated on input */ - - /* U */ - double *Ux, /* size nnz(U) */ -#ifdef COMPLEX - double *Uz, /* size nnz(U) for the complex case, ignored if real */ -#endif - - Int *P, - Int *Q, - double *Rs, - - KLU_common *Common -) -{ - Entry *Ukk ; - Int block, k1, k2, kk, i, n, nk, nblocks, nz ; - - if (Common == NULL) - { - return (FALSE) ; - } - - if (Common->status == KLU_EMPTY_MATRIX) - { - return (FALSE) ; - } - - if (Symbolic == NULL || Numeric == NULL) - { - Common->status = KLU_INVALID ; - return (FALSE) ; - } - - Common->status = KLU_OK ; - n = Symbolic->n ; - nblocks = Symbolic->nblocks ; - - /* ---------------------------------------------------------------------- */ - /* extract scale factors */ - /* ---------------------------------------------------------------------- */ - - if (Rs != NULL) - { - if (Numeric->Rs != NULL) - { - for (i = 0 ; i < n ; i++) - { - Rs [i] = Numeric->Rs [i] ; - } - } - else - { - /* no scaling */ - for (i = 0 ; i < n ; i++) - { - Rs [i] = 1 ; - } - } - } - - /* ---------------------------------------------------------------------- */ - /* extract final row permutation */ - /* ---------------------------------------------------------------------- */ - - if (P != NULL) - { - for (i = 0 ; i < n ; i++) - { - P [i] = Numeric->Pnum [i] ; - } - } - - /* ---------------------------------------------------------------------- */ - /* extract column permutation */ - /* ---------------------------------------------------------------------- */ - - if (Q != NULL) - { - for (i = 0 ; i < n ; i++) - { - Q [i] = Symbolic->Q [i] ; - } - } - - /* ---------------------------------------------------------------------- */ - /* extract each block of U */ - /* ---------------------------------------------------------------------- */ - - if (Ux != NULL -#ifdef COMPLEX - && Uz != NULL -#endif - ) - { - nz = 0 ; - for (block = 0 ; block < nblocks ; block++) - { - k1 = Symbolic->R [block] ; - k2 = Symbolic->R [block+1] ; - nk = k2 - k1 ; - Ukk = ((Entry *) Numeric->Udiag) + k1 ; - if (nk == 1) - { - /* singleton block */ - Ux [nz] = REAL (Ukk [0]) ; -#ifdef COMPLEX - Uz [nz] = IMAG (Ukk [0]) ; -#endif - nz++ ; - } - else - { - /* non-singleton block */ - for (kk = 0 ; kk < nk ; kk++) - { - /* add the diagonal entry */ - Ux [nz] = REAL (Ukk [kk]) ; -#ifdef COMPLEX - Uz [nz] = IMAG (Ukk [kk]) ; -#endif - nz++ ; - } - } - } - ASSERT (nz == Numeric->unz) ; - - } - - return (TRUE) ; -} - -Int KLU_print -( - Int *Ap, - Int *Ai, - double *Ax, - int n, - int *IntToExtRowMap, - int *IntToExtColMap -) -{ - Entry *Az ; - int i, j ; - - Az = (Entry *)Ax ; - for (i = 0 ; i < n ; i++) - { - for (j = Ap [i] ; j < Ap [i + 1] ; j++) - { - -#ifdef COMPLEX - if (IntToExtRowMap && IntToExtColMap) { - fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g j%-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j].Real, Az [j].Imag) ; - } else { - fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g j%-.9g\n", Ai [j] + 1, i + 1, Az [j].Real, Az [j].Imag) ; - } -#else - if (IntToExtRowMap && IntToExtColMap) { - fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j]) ; - } else { - fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g\n", Ai [j] + 1, i + 1, Az [j]) ; - } -#endif - - } - } - - return 0 ; -} diff --git a/src/maths/KLU/klu_utils.c b/src/maths/KLU/klu_utils.c index 2e1ad246f..becb7c446 100644 --- a/src/maths/KLU/klu_utils.c +++ b/src/maths/KLU/klu_utils.c @@ -171,3 +171,181 @@ KLU_convert_matrix_in_CSR /* return TRUE if successful, FALSE otherwise return (TRUE) ; } + +/* Francesco - Extract only Udiag */ +Int KLU_extract_Udiag /* returns TRUE if successful, FALSE otherwise */ +( + /* inputs: */ + KLU_numeric *Numeric, + KLU_symbolic *Symbolic, + + /* outputs, all of which must be allocated on input */ + + /* U */ + double *Ux, /* size nnz(U) */ +#ifdef COMPLEX + double *Uz, /* size nnz(U) for the complex case, ignored if real */ +#endif + + Int *P, + Int *Q, + double *Rs, + + KLU_common *Common +) +{ + Entry *Ukk ; + Int block, k1, k2, kk, i, n, nk, nblocks, nz ; + + if (Common == NULL) + { + return (FALSE) ; + } + + if (Common->status == KLU_EMPTY_MATRIX) + { + return (FALSE) ; + } + + if (Symbolic == NULL || Numeric == NULL) + { + Common->status = KLU_INVALID ; + return (FALSE) ; + } + + Common->status = KLU_OK ; + n = Symbolic->n ; + nblocks = Symbolic->nblocks ; + + /* ---------------------------------------------------------------------- */ + /* extract scale factors */ + /* ---------------------------------------------------------------------- */ + + if (Rs != NULL) + { + if (Numeric->Rs != NULL) + { + for (i = 0 ; i < n ; i++) + { + Rs [i] = Numeric->Rs [i] ; + } + } + else + { + /* no scaling */ + for (i = 0 ; i < n ; i++) + { + Rs [i] = 1 ; + } + } + } + + /* ---------------------------------------------------------------------- */ + /* extract final row permutation */ + /* ---------------------------------------------------------------------- */ + + if (P != NULL) + { + for (i = 0 ; i < n ; i++) + { + P [i] = Numeric->Pnum [i] ; + } + } + + /* ---------------------------------------------------------------------- */ + /* extract column permutation */ + /* ---------------------------------------------------------------------- */ + + if (Q != NULL) + { + for (i = 0 ; i < n ; i++) + { + Q [i] = Symbolic->Q [i] ; + } + } + + /* ---------------------------------------------------------------------- */ + /* extract each block of U */ + /* ---------------------------------------------------------------------- */ + + if (Ux != NULL +#ifdef COMPLEX + && Uz != NULL +#endif + ) + { + nz = 0 ; + for (block = 0 ; block < nblocks ; block++) + { + k1 = Symbolic->R [block] ; + k2 = Symbolic->R [block+1] ; + nk = k2 - k1 ; + Ukk = ((Entry *) Numeric->Udiag) + k1 ; + if (nk == 1) + { + /* singleton block */ + Ux [nz] = REAL (Ukk [0]) ; +#ifdef COMPLEX + Uz [nz] = IMAG (Ukk [0]) ; +#endif + nz++ ; + } + else + { + /* non-singleton block */ + for (kk = 0 ; kk < nk ; kk++) + { + /* add the diagonal entry */ + Ux [nz] = REAL (Ukk [kk]) ; +#ifdef COMPLEX + Uz [nz] = IMAG (Ukk [kk]) ; +#endif + nz++ ; + } + } + } + ASSERT (nz == Numeric->unz) ; + + } + + return (TRUE) ; +} + +Int KLU_print +( + Int *Ap, + Int *Ai, + double *Ax, + int n, + int *IntToExtRowMap, + int *IntToExtColMap +) +{ + Entry *Az ; + int i, j ; + + Az = (Entry *)Ax ; + for (i = 0 ; i < n ; i++) + { + for (j = Ap [i] ; j < Ap [i + 1] ; j++) + { + +#ifdef COMPLEX + if (IntToExtRowMap && IntToExtColMap) { + fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g j%-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j].Real, Az [j].Imag) ; + } else { + fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g j%-.9g\n", Ai [j] + 1, i + 1, Az [j].Real, Az [j].Imag) ; + } +#else + if (IntToExtRowMap && IntToExtColMap) { + fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g\n", IntToExtRowMap [Ai [j] + 1], IntToExtColMap [i + 1], Az [j]) ; + } else { + fprintf (stderr, "Row: %d\tCol: %d\tValue: %-.9g\n", Ai [j] + 1, i + 1, Az [j]) ; + } +#endif + + } + } + + return 0 ; +}