integrate OSDI with KLU

compile_linux_klu.sh

@@ -43,13 +43,13 @@ if test "$1" = "d"; then
    if [ $? -ne 0 ]; then  echo "cd debug failed"; exit 1 ; fi
   echo "configuring for 64 bit debug"
   echo
-  ../configure --with-x  --enable-cider --with-readline=yes --enable-openmp --enable-xspice --enable-klu CFLAGS="-g -m64 -O0 -Wall -Wno-unused-but-set-variable" LDFLAGS="-m64 -g"
+  ../configure --with-x  --enable-cider --with-readline=yes --enable-openmp --enable-xspice --enable-klu --enable-predictor --enable-osdi    CFLAGS="-g -m64 -O0 -Wall -Wno-unused-but-set-variable" LDFLAGS="-m64 -g"
 else
    cd release
    if [ $? -ne 0 ]; then  echo "cd release failed"; exit 1 ; fi
   echo "configuring for 64 bit release"
   echo 
-  ../configure --with-x  --enable-cider --with-readline=yes --enable-openmp --enable-xspice --enable-klu --disable-debug CFLAGS="-m64 -O2" LDFLAGS="-m64 -s"
+  ../configure --with-x  --enable-cider --with-readline=yes --enable-openmp --enable-xspice --enable-klu --enable-predictor --enable-osdi --disable-debug CFLAGS="-m64 -O2" LDFLAGS="-m64 -s"
 fi
 if [ $? -ne 0 ]; then  echo "../configure failed"; exit 1 ; fi
 

src/include/ngspice/osdiitf.h

@@ -11,6 +11,7 @@
 
 #pragma once
 
+#include "ngspice/config.h"
 #include "ngspice/devdefs.h"
 #include <stdint.h>
 
@@ -19,7 +20,13 @@ typedef struct OsdiRegistryEntry {
   uint32_t inst_offset;
   uint32_t dt;
   uint32_t temp;
+
   bool has_m;
+
+#ifdef KLU
+  uint32_t matrix_ptr_offset;
+#endif
+
 } OsdiRegistryEntry;
 
 typedef struct OsdiObjectFile {

src/osdi/osdidefs.h

@@ -73,9 +73,15 @@ typedef struct OsdiModelData {
 } OsdiModelData;
 
 extern size_t osdi_instance_data_off(const OsdiRegistryEntry *entry);
-extern void *osdi_instance_data(const OsdiRegistryEntry *entry, GENinstance *inst);
+extern void *osdi_instance_data(const OsdiRegistryEntry *entry,
-extern OsdiExtraInstData *osdi_extra_instance_data(const OsdiRegistryEntry *entry,
                                 GENinstance *inst);
+#ifdef KLU
+extern size_t osdi_instance_matrix_ptr_off(const OsdiRegistryEntry *entry);
+extern double **osdi_instance_matrix_ptr(const OsdiRegistryEntry *entry,
+                                         GENinstance *inst);
+#endif
+extern OsdiExtraInstData *
+osdi_extra_instance_data(const OsdiRegistryEntry *entry, GENinstance *inst);
 extern size_t osdi_model_data_off(void);
 extern void *osdi_model_data(GENmodel *model);
 extern void *osdi_model_data_from_inst(GENinstance *inst);

src/osdi/osdiext.h

@@ -26,7 +26,13 @@ extern int OSDIload(GENmodel *, CKTcircuit *);
 extern int OSDItemp(GENmodel *, CKTcircuit *);
 extern int OSDIacLoad(GENmodel *, CKTcircuit *);
 extern int OSDItrunc(GENmodel *, CKTcircuit *, double *);
-extern int OSDIpzLoad(GENmodel*, CKTcircuit*, SPcomplex*);
+extern int OSDIpzLoad(GENmodel *, CKTcircuit *, SPcomplex *);
+
+#ifdef KLU
+extern int OSDIbindCSC(GENmodel *inModel, CKTcircuit *ckt);
+extern int OSDIbindCSCComplexToReal(GENmodel *inModel, CKTcircuit *ckt);
+extern int OSDIbindCSCComplex(GENmodel *inModel, CKTcircuit *ckt);
+#endif
 
 /* extern int OSDIconvTest(GENmodel*,CKTcircuit*); */
 /* extern int OSDImDelete(GENmodel*); */

src/osdi/osdiinit.c

@@ -197,5 +197,11 @@ extern SPICEdev *osdi_create_spicedev(const OsdiRegistryEntry *entry) {
   OSDIinfo->DEVpzLoad = OSDIpzLoad;
   OSDIinfo->DEVtrunc = OSDItrunc;
 
+  #ifdef KLU
+  OSDIinfo->DEVbindCSC = OSDIbindCSC;
+  OSDIinfo->DEVbindCSCComplex = OSDIbindCSCComplex;
+  OSDIinfo->DEVbindCSCComplexToReal = OSDIbindCSCComplexToReal;
+  #endif
+
   return OSDIinfo;
 }

src/osdi/osdiregistry.c

@@ -9,11 +9,13 @@
  * Author: Pascal Kuthe <pascal.kuthe@semimod.de>
  */
 
+#include "ngspice/config.h"
 #include "ngspice/hash.h"
 #include "ngspice/memory.h"
 #include "ngspice/stringutil.h"
 #include "osdidefs.h"
 
+#include <stddef.h>
 #include <sys/stat.h>
 
 #include "osdi.h"
@@ -199,6 +201,13 @@ static char *resolve_input_path(const char *name) {
   return NULL;
 } /* end of function inp_pathresolve_at */
 
+static size_t pad_to_align(size_t alignment, size_t size) {
+  size_t padding = alignment - size % alignment;
+  if (padding == alignment) {
+    return 0;
+  }
+  return padding + size;
+}
 /**
  * Calculates the offset that the OSDI instance data has from the beginning of
  * the instance data allocated by ngspice. This offset is non trivial because
@@ -209,13 +218,21 @@ static char *resolve_input_path(const char *name) {
 static size_t calc_osdi_instance_data_off(const OsdiDescriptor *descr) {
   size_t res = sizeof(GENinstance) /* generic data */
                + descr->num_terminals * sizeof(int);
-  size_t padding = sizeof(max_align_t) - res % sizeof(max_align_t);
+#ifdef KLU
-  if (padding == sizeof(max_align_t)) {
+  res = pad_to_align(alignof(double *), res);
-    padding = 0;
+  res += ((size_t)descr->num_jacobian_entries) * 2 * sizeof(double *);
-  }
+#endif
-  return res + padding;
+  return pad_to_align(alignof(max_align_t), res);
 }
 
+#ifdef KLU
+static size_t calc_osdi_instance_matrix_off(const OsdiDescriptor *descr) {
+  size_t res = sizeof(GENinstance) /* generic data */
+               + descr->num_terminals * sizeof(int);
+  return pad_to_align(alignof(double *), res);
+}
+#endif
+
 #define INVALID_OBJECT                                                         \
   (OsdiObjectFile) { .num_entries = -1 }
 
@@ -386,6 +403,11 @@ extern OsdiObjectFile load_object_file(const char *input) {
         .dt = dt,
         .temp = temp,
         .has_m = has_m,
+
+#ifdef KLU
+        .matrix_ptr_offset = (uint32_t)calc_osdi_instance_matrix_off(descr),
+#endif
+
     };
   }
 
@@ -399,6 +421,15 @@ extern OsdiObjectFile load_object_file(const char *input) {
 inline size_t osdi_instance_data_off(const OsdiRegistryEntry *entry) {
   return entry->inst_offset;
 }
+#ifdef KLU
+inline size_t osdi_instance_matrix_ptr_off(const OsdiRegistryEntry *entry) {
+  return entry->matrix_ptr_offset;
+}
+inline double **osdi_instance_matrix_ptr(const OsdiRegistryEntry *entry,
+                                         GENinstance *inst) {
+  return (double **)(((char *)inst) + osdi_instance_matrix_ptr_off(entry));
+}
+#endif
 
 inline void *osdi_instance_data(const OsdiRegistryEntry *entry,
                                 GENinstance *inst) {

src/osdi/osdisetup.c

@@ -10,6 +10,7 @@
  */
 
 #include "ngspice/iferrmsg.h"
+#include "ngspice/klu.h"
 #include "ngspice/memory.h"
 #include "ngspice/ngspice.h"
 #include "ngspice/typedefs.h"
@@ -408,3 +409,123 @@ extern int OSDIunsetup(GENmodel *inModel, CKTcircuit *ckt) {
   }
   return (OK);
 }
+#ifdef KLU
+#include "ngspice/klu-binding.h"
+static int init_matrix_klu(SMPmatrix *matrix, const OsdiDescriptor *descr,
+                           void *inst, double **inst_matrix_ptrs) {
+  BindElement tmp;
+  BindElement *matched;
+  BindElement *bindings = matrix->SMPkluMatrix->KLUmatrixBindStructCOO;
+  size_t nz = (size_t)matrix->SMPkluMatrix->KLUmatrixLinkedListNZ;
+  uint32_t *node_mapping =
+      (uint32_t *)(((char *)inst) + descr->node_mapping_offset);
+  double **jacobian_ptr_resist =
+      (double **)(((char *)inst) + descr->jacobian_ptr_resist_offset);
+
+  for (uint32_t i = 0; i < descr->num_jacobian_entries; i++) {
+    uint32_t equation = descr->jacobian_entries[i].nodes.node_1;
+    uint32_t unkown = descr->jacobian_entries[i].nodes.node_2;
+    equation = node_mapping[equation];
+    unkown = node_mapping[unkown];
+    if (equation != 0 && unkown != 0) {
+      tmp.COO = jacobian_ptr_resist[i];
+      tmp.CSC = NULL;
+      tmp.CSC_Complex = NULL;
+      matched = (BindElement *)bsearch(&tmp, bindings, nz, sizeof(BindElement),
+                                       BindCompare);
+      if (matched == NULL) {
+        printf("Ptr %p not found in BindStruct Table\n",
+               jacobian_ptr_resist[i]);
+        return (E_PANIC);
+      }
+      uint32_t react_off = descr->jacobian_entries[i].react_ptr_off;
+      // complex number for ac analysis
+      if (react_off != UINT32_MAX) {
+        double **jacobian_ptr_react = (double **)(((char *)inst) + react_off);
+        *jacobian_ptr_react = matched->CSC_Complex + 1;
+      }
+      jacobian_ptr_resist[i] = matched->CSC;
+      inst_matrix_ptrs[2 * i] = matched->CSC;
+      inst_matrix_ptrs[2 * i + 1] = matched->CSC_Complex;
+    }
+  }
+  return (OK);
+}
+
+static int update_matrix_klu(const OsdiDescriptor *descr, void *inst,
+                             double **inst_matrix_ptrs, bool complex) {
+  uint32_t *node_mapping =
+      (uint32_t *)(((char *)inst) + descr->node_mapping_offset);
+  double **jacobian_ptr_resist =
+      (double **)(((char *)inst) + descr->jacobian_ptr_resist_offset);
+
+  for (uint32_t i = 0; i < descr->num_jacobian_entries; i++) {
+    uint32_t equation = descr->jacobian_entries[i].nodes.node_1;
+    uint32_t unkown = descr->jacobian_entries[i].nodes.node_2;
+    equation = node_mapping[equation];
+    unkown = node_mapping[unkown];
+    if (equation != 0 && unkown != 0) {
+      jacobian_ptr_resist[i] = inst_matrix_ptrs[2 * i + complex];
+    }
+  }
+  return (OK);
+}
+
+int OSDIbindCSC(GENmodel *inModel, CKTcircuit *ckt) {
+
+  OsdiRegistryEntry *entry = osdi_reg_entry_model(inModel);
+  const OsdiDescriptor *descr = entry->descriptor;
+  GENmodel *gen_model;
+  GENinstance *gen_inst;
+
+  /* setup a temporary buffer */
+  uint32_t *node_ids = TMALLOC(uint32_t, descr->num_nodes);
+
+  for (gen_model = inModel; gen_model; gen_model = gen_model->GENnextModel) {
+    void *model = osdi_model_data(gen_model);
+    for (gen_inst = gen_model->GENinstances; gen_inst;
+         gen_inst = gen_inst->GENnextInstance) {
+      void *inst = osdi_instance_data(entry, gen_inst);
+      double **matrix_ptrs = osdi_instance_matrix_ptr(entry, gen_inst);
+      int err = init_matrix_klu(ckt->CKTmatrix, descr, inst, matrix_ptrs);
+      if (err != (OK)) {
+        return err;
+      }
+    }
+  }
+
+  return (OK);
+}
+int OSDIupdateCSC(GENmodel *inModel, CKTcircuit *ckt, bool complex) {
+
+  OsdiRegistryEntry *entry = osdi_reg_entry_model(inModel);
+  const OsdiDescriptor *descr = entry->descriptor;
+  GENmodel *gen_model;
+  GENinstance *gen_inst;
+
+  /* setup a temporary buffer */
+  uint32_t *node_ids = TMALLOC(uint32_t, descr->num_nodes);
+
+  for (gen_model = inModel; gen_model; gen_model = gen_model->GENnextModel) {
+    void *model = osdi_model_data(gen_model);
+    for (gen_inst = gen_model->GENinstances; gen_inst;
+         gen_inst = gen_inst->GENnextInstance) {
+      void *inst = osdi_instance_data(entry, gen_inst);
+      double **matrix_ptrs = osdi_instance_matrix_ptr(entry, gen_inst);
+      int err = update_matrix_klu(descr, inst, matrix_ptrs, complex);
+      if (err != (OK)) {
+        return err;
+      }
+    }
+  }
+
+  return (OK);
+}
+int OSDIbindCSCComplexToReal(GENmodel *inModel, CKTcircuit *ckt) {
+  OSDIupdateCSC(inModel, ckt, false);
+}
+
+int OSDIbindCSCComplex(GENmodel *inModel, CKTcircuit *ckt) {
+  OSDIupdateCSC(inModel, ckt, true);
+}
+#endif