Squashed 'third_party/osqp/' content from commit 33454b3e23
Change-Id: I056df0582ca06664e86554c341a94c47ab932001
git-subtree-dir: third_party/osqp
git-subtree-split: 33454b3e236f1f44193bfbbb6b8c8e71f8f04e9a
Signed-off-by: Austin Schuh <austin.linux@gmail.com>
diff --git a/lin_sys/direct/CMakeLists.txt b/lin_sys/direct/CMakeLists.txt
new file mode 100644
index 0000000..c36c375
--- /dev/null
+++ b/lin_sys/direct/CMakeLists.txt
@@ -0,0 +1,44 @@
+# Add direct linear systems solvers
+
+# QDLDL (Default)
+# -------------------------
+add_subdirectory(qdldl)
+
+# Need to add qdldlobject only here because it cannot be
+# included in another object library such as linsys_qdldl
+set(direct_linsys_solvers $<TARGET_OBJECTS:linsys_qdldl> $<TARGET_OBJECTS:qdldlobject>)
+
+# NB. The second directory is added because we need to include the "qdldl_types.h" file in "qdldl_interface.h"
+set(direct_linsys_solvers_includes
+ "${CMAKE_CURRENT_SOURCE_DIR}/qdldl/"
+ "${CMAKE_CURRENT_SOURCE_DIR}/qdldl/qdldl_sources/include")
+
+
+# Add other solvers if embedded option is false
+if(NOT DEFINED EMBEDDED)
+
+# MKL Pardiso MKL
+# -----------
+# If MKL Pardiso is enabled, include pardiso directory
+if (ENABLE_MKL_PARDISO)
+ # Add Pardiso interface
+ add_subdirectory(pardiso)
+
+ add_library(linsys_pardiso OBJECT ${pardiso_sources})
+
+ # Add parent directory for library handler
+ target_include_directories(linsys_pardiso PRIVATE ${pardiso_includes} ${PROJECT_SOURCE_DIR}/include)
+
+ set(direct_linsys_solvers ${direct_linsys_solvers} $<TARGET_OBJECTS:linsys_pardiso>)
+
+ set(direct_linsys_solvers_includes "${direct_linsys_solvers_includes};${CMAKE_CURRENT_SOURCE_DIR}/pardiso/")
+endif()
+
+
+endif()
+
+# Make direct_linsys_solvers list visible from parent directory
+set(direct_linsys_solvers ${direct_linsys_solvers} PARENT_SCOPE)
+
+# Make external linsys solvers includes visible from parent directory
+set(direct_linsys_solvers_includes ${direct_linsys_solvers_includes} PARENT_SCOPE)
diff --git a/lin_sys/direct/pardiso/CMakeLists.txt b/lin_sys/direct/pardiso/CMakeLists.txt
new file mode 100644
index 0000000..83a3f86
--- /dev/null
+++ b/lin_sys/direct/pardiso/CMakeLists.txt
@@ -0,0 +1,14 @@
+set(pardiso_includes
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ PARENT_SCOPE
+)
+
+# Set source files
+set(
+ pardiso_sources
+ ${CMAKE_CURRENT_SOURCE_DIR}/pardiso_interface.h
+ ${CMAKE_CURRENT_SOURCE_DIR}/pardiso_interface.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/pardiso_loader.h
+ ${CMAKE_CURRENT_SOURCE_DIR}/pardiso_loader.c
+ PARENT_SCOPE
+)
diff --git a/lin_sys/direct/pardiso/pardiso_interface.c b/lin_sys/direct/pardiso/pardiso_interface.c
new file mode 100644
index 0000000..738334a
--- /dev/null
+++ b/lin_sys/direct/pardiso/pardiso_interface.c
@@ -0,0 +1,300 @@
+#include "pardiso_interface.h"
+
+#define MKL_INT c_int
+
+// Single Dynamic library interface
+#define MKL_INTERFACE_LP64 0x0
+#define MKL_INTERFACE_ILP64 0x1
+
+// Solver Phases
+#define PARDISO_SYMBOLIC (11)
+#define PARDISO_NUMERIC (22)
+#define PARDISO_SOLVE (33)
+#define PARDISO_CLEANUP (-1)
+
+
+// Prototypes for Pardiso functions
+void pardiso(void**, // pt
+ const c_int*, // maxfct
+ const c_int*, // mnum
+ const c_int*, // mtype
+ const c_int*, // phase
+ const c_int*, // n
+ const c_float*, // a
+ const c_int*, // ia
+ const c_int*, // ja
+ c_int*, // perm
+ const c_int*, //nrhs
+ c_int*, // iparam
+ const c_int*, //msglvl
+ c_float*, // b
+ c_float*, // x
+ c_int* // error
+ );
+c_int mkl_set_interface_layer(c_int);
+c_int mkl_get_max_threads();
+
+// Free LDL Factorization structure
+void free_linsys_solver_pardiso(pardiso_solver *s) {
+ if (s) {
+
+ // Free pardiso solver using internal function
+ s->phase = PARDISO_CLEANUP;
+ pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
+ &(s->nKKT), &(s->fdum), s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
+ s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
+
+ if ( s->error != 0 ){
+#ifdef PRINTING
+ c_eprint("Error during MKL Pardiso cleanup: %d", (int)s->error);
+#endif
+ }
+ // Check each attribute of the structure and free it if it exists
+ if (s->KKT) csc_spfree(s->KKT);
+ if (s->KKT_i) c_free(s->KKT_i);
+ if (s->KKT_p) c_free(s->KKT_p);
+ if (s->bp) c_free(s->bp);
+ if (s->sol) c_free(s->sol);
+ if (s->rho_inv_vec) c_free(s->rho_inv_vec);
+
+ // These are required for matrix updates
+ if (s->Pdiag_idx) c_free(s->Pdiag_idx);
+ if (s->PtoKKT) c_free(s->PtoKKT);
+ if (s->AtoKKT) c_free(s->AtoKKT);
+ if (s->rhotoKKT) c_free(s->rhotoKKT);
+
+ c_free(s);
+
+ }
+}
+
+
+// Initialize factorization structure
+c_int init_linsys_solver_pardiso(pardiso_solver ** sp, const csc * P, const csc * A, c_float sigma, const c_float * rho_vec, c_int polish){
+ c_int i; // loop counter
+ c_int nnzKKT; // Number of nonzeros in KKT
+ // Define Variables
+ c_int n_plus_m; // n_plus_m dimension
+
+
+ // Allocate private structure to store KKT factorization
+ pardiso_solver *s;
+ s = c_calloc(1, sizeof(pardiso_solver));
+ *sp = s;
+
+ // Size of KKT
+ s->n = P->n;
+ s->m = A->m;
+ n_plus_m = s->n + s->m;
+ s->nKKT = n_plus_m;
+
+ // Sigma parameter
+ s->sigma = sigma;
+
+ // Polishing flag
+ s->polish = polish;
+
+ // Link Functions
+ s->solve = &solve_linsys_pardiso;
+ s->free = &free_linsys_solver_pardiso;
+ s->update_matrices = &update_linsys_solver_matrices_pardiso;
+ s->update_rho_vec = &update_linsys_solver_rho_vec_pardiso;
+
+ // Assign type
+ s->type = MKL_PARDISO_SOLVER;
+
+ // Working vector
+ s->bp = (c_float *)c_malloc(sizeof(c_float) * n_plus_m);
+
+ // Solution vector
+ s->sol = (c_float *)c_malloc(sizeof(c_float) * n_plus_m);
+
+ // Parameter vector
+ s->rho_inv_vec = (c_float *)c_malloc(sizeof(c_float) * n_plus_m);
+
+ // Form KKT matrix
+ if (polish){ // Called from polish()
+ // Use s->rho_inv_vec for storing param2 = vec(delta)
+ for (i = 0; i < A->m; i++){
+ s->rho_inv_vec[i] = sigma;
+ }
+
+ s->KKT = form_KKT(P, A, 1, sigma, s->rho_inv_vec, OSQP_NULL, OSQP_NULL, OSQP_NULL, OSQP_NULL, OSQP_NULL);
+ }
+ else { // Called from ADMM algorithm
+
+ // Allocate vectors of indices
+ s->PtoKKT = c_malloc((P->p[P->n]) * sizeof(c_int));
+ s->AtoKKT = c_malloc((A->p[A->n]) * sizeof(c_int));
+ s->rhotoKKT = c_malloc((A->m) * sizeof(c_int));
+
+ // Use s->rho_inv_vec for storing param2 = rho_inv_vec
+ for (i = 0; i < A->m; i++){
+ s->rho_inv_vec[i] = 1. / rho_vec[i];
+ }
+
+ s->KKT = form_KKT(P, A, 1, sigma, s->rho_inv_vec,
+ s->PtoKKT, s->AtoKKT,
+ &(s->Pdiag_idx), &(s->Pdiag_n), s->rhotoKKT);
+ }
+
+ // Check if matrix has been created
+ if (!(s->KKT)) {
+#ifdef PRINTING
+ c_eprint("Error in forming KKT matrix");
+#endif
+ free_linsys_solver_pardiso(s);
+ return OSQP_LINSYS_SOLVER_INIT_ERROR;
+ } else {
+ // Adjust indexing for Pardiso
+ nnzKKT = s->KKT->p[s->KKT->m];
+ s->KKT_i = c_malloc((nnzKKT) * sizeof(c_int));
+ s->KKT_p = c_malloc((s->KKT->m + 1) * sizeof(c_int));
+
+ for(i = 0; i < nnzKKT; i++){
+ s->KKT_i[i] = s->KKT->i[i] + 1;
+ }
+ for(i = 0; i < n_plus_m+1; i++){
+ s->KKT_p[i] = s->KKT->p[i] + 1;
+ }
+
+ }
+
+ // Set MKL interface layer (Long integers if activated)
+#ifdef DLONG
+ mkl_set_interface_layer(MKL_INTERFACE_ILP64);
+#else
+ mkl_set_interface_layer(MKL_INTERFACE_LP64);
+#endif
+
+ // Set Pardiso variables
+ s->mtype = -2; // Real symmetric indefinite matrix
+ s->nrhs = 1; // Number of right hand sides
+ s->maxfct = 1; // Maximum number of numerical factorizations
+ s->mnum = 1; // Which factorization to use
+ s->msglvl = 0; // Do not print statistical information
+ s->error = 0; // Initialize error flag
+ for ( i = 0; i < 64; i++ ) {
+ s->iparm[i] = 0; // Setup Pardiso control parameters
+ s->pt[i] = 0; // Initialize the internal solver memory pointer
+ }
+ s->iparm[0] = 1; // No solver default
+ s->iparm[1] = 3; // Fill-in reordering from OpenMP
+ if (polish) {
+ s->iparm[5] = 1; // Write solution into b
+ } else {
+ s->iparm[5] = 0; // Do NOT write solution into b
+ }
+ /* s->iparm[7] = 2; // Max number of iterative refinement steps */
+ s->iparm[7] = 0; // Number of iterative refinement steps (auto, performs them only if perturbed pivots are obtained)
+ s->iparm[9] = 13; // Perturb the pivot elements with 1E-13
+ s->iparm[34] = 0; // Use Fortran-style indexing for indices
+ /* s->iparm[34] = 1; // Use C-style indexing for indices */
+
+ // Print number of threads
+ s->nthreads = mkl_get_max_threads();
+
+ // Reordering and symbolic factorization
+ s->phase = PARDISO_SYMBOLIC;
+ pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
+ &(s->nKKT), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
+ s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
+ if ( s->error != 0 ){
+#ifdef PRINTING
+ c_eprint("Error during symbolic factorization: %d", (int)s->error);
+#endif
+ free_linsys_solver_pardiso(s);
+ *sp = OSQP_NULL;
+ return OSQP_LINSYS_SOLVER_INIT_ERROR;
+ }
+
+ // Numerical factorization
+ s->phase = PARDISO_NUMERIC;
+ pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
+ &(s->nKKT), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
+ s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
+ if ( s->error ){
+#ifdef PRINTING
+ c_eprint("Error during numerical factorization: %d", (int)s->error);
+#endif
+ free_linsys_solver_pardiso(s);
+ *sp = OSQP_NULL;
+ return OSQP_LINSYS_SOLVER_INIT_ERROR;
+ }
+
+
+ // No error
+ return 0;
+}
+
+// Returns solution to linear system Ax = b with solution stored in b
+c_int solve_linsys_pardiso(pardiso_solver * s, c_float * b) {
+ c_int j;
+
+ // Back substitution and iterative refinement
+ s->phase = PARDISO_SOLVE;
+ pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
+ &(s->nKKT), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
+ s->iparm, &(s->msglvl), b, s->sol, &(s->error));
+ if ( s->error != 0 ){
+#ifdef PRINTING
+ c_eprint("Error during linear system solution: %d", (int)s->error);
+#endif
+ return 1;
+ }
+
+ if (!(s->polish)) {
+ /* copy x_tilde from s->sol */
+ for (j = 0 ; j < s->n ; j++) {
+ b[j] = s->sol[j];
+ }
+
+ /* compute z_tilde from b and s->sol */
+ for (j = 0 ; j < s->m ; j++) {
+ b[j + s->n] += s->rho_inv_vec[j] * s->sol[j + s->n];
+ }
+ }
+
+ return 0;
+}
+
+// Update solver structure with new P and A
+c_int update_linsys_solver_matrices_pardiso(pardiso_solver * s, const csc *P, const csc *A) {
+
+ // Update KKT matrix with new P
+ update_KKT_P(s->KKT, P, s->PtoKKT, s->sigma, s->Pdiag_idx, s->Pdiag_n);
+
+ // Update KKT matrix with new A
+ update_KKT_A(s->KKT, A, s->AtoKKT);
+
+ // Perform numerical factorization
+ s->phase = PARDISO_NUMERIC;
+ pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
+ &(s->nKKT), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
+ s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
+
+ // Return exit flag
+ return s->error;
+}
+
+
+c_int update_linsys_solver_rho_vec_pardiso(pardiso_solver * s, const c_float * rho_vec) {
+ c_int i;
+
+ // Update internal rho_inv_vec
+ for (i = 0; i < s->m; i++){
+ s->rho_inv_vec[i] = 1. / rho_vec[i];
+ }
+
+ // Update KKT matrix with new rho_vec
+ update_KKT_param2(s->KKT, s->rho_inv_vec, s->rhotoKKT, s->m);
+
+ // Perform numerical factorization
+ s->phase = PARDISO_NUMERIC;
+ pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
+ &(s->nKKT), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
+ s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
+
+ // Return exit flag
+ return s->error;
+}
diff --git a/lin_sys/direct/pardiso/pardiso_interface.h b/lin_sys/direct/pardiso/pardiso_interface.h
new file mode 100644
index 0000000..c0e334d
--- /dev/null
+++ b/lin_sys/direct/pardiso/pardiso_interface.h
@@ -0,0 +1,127 @@
+#ifndef PARDISO_H
+#define PARDISO_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include "lin_alg.h"
+#include "kkt.h"
+
+/**
+ * Pardiso solver structure
+ *
+ * NB: If we use Pardiso, we suppose that EMBEDDED is not enabled
+ */
+typedef struct pardiso pardiso_solver;
+
+struct pardiso {
+ enum linsys_solver_type type;
+
+ /**
+ * @name Functions
+ * @{
+ */
+ c_int (*solve)(struct pardiso * self, c_float * b);
+
+ void (*free)(struct pardiso * self); ///< Free workspace (only if desktop)
+
+ c_int (*update_matrices)(struct pardiso * self, const csc *P, const csc *A); ///< Update solver matrices
+ c_int (*update_rho_vec)(struct pardiso * self, const c_float * rho_vec); ///< Update rho_vec parameter
+
+ c_int nthreads;
+ /** @} */
+
+
+ /**
+ * @name Attributes
+ * @{
+ */
+ // Attributes
+ csc *KKT; ///< KKT matrix (in CSR format!)
+ c_int *KKT_i; ///< KKT column indices in 1-indexing for Pardiso
+ c_int *KKT_p; ///< KKT row pointers in 1-indexing for Pardiso
+ c_float *bp; ///< workspace memory for solves (rhs)
+ c_float *sol; ///< solution to the KKT system
+ c_float *rho_inv_vec; ///< parameter vector
+ c_float sigma; ///< scalar parameter
+ c_int polish; ///< polishing flag
+ c_int n; ///< number of QP variables
+ c_int m; ///< number of QP constraints
+
+ // Pardiso variables
+ void *pt[64]; ///< internal solver memory pointer pt
+ c_int iparm[64]; ///< Pardiso control parameters
+ c_int nKKT; ///< dimension of the linear system
+ c_int mtype; ///< matrix type (-2 for real and symmetric indefinite)
+ c_int nrhs; ///< number of right-hand sides (1 for our needs)
+ c_int maxfct; ///< maximum number of factors (1 for our needs)
+ c_int mnum; ///< indicates matrix for the solution phase (1 for our needs)
+ c_int phase; ///< control the execution phases of the solver
+ c_int error; ///< the error indicator (0 for no error)
+ c_int msglvl; ///< Message level information (0 for no output)
+ c_int idum; ///< dummy integer
+ c_float fdum; ///< dummy float
+
+ // These are required for matrix updates
+ c_int * Pdiag_idx, Pdiag_n; ///< index and number of diagonal elements in P
+ c_int * PtoKKT, * AtoKKT; ///< Index of elements from P and A to KKT matrix
+ c_int * rhotoKKT; ///< Index of rho places in KKT matrix
+
+ /** @} */
+};
+
+
+/**
+ * Initialize Pardiso Solver
+ *
+ * @param s Pointer to a private structure
+ * @param P Cost function matrix (upper triangular form)
+ * @param A Constraints matrix
+ * @param sigma Algorithm parameter. If polish, then sigma = delta.
+ * @param rho_vec Algorithm parameter. If polish, then rho_vec = OSQP_NULL.
+ * @param polish Flag whether we are initializing for polish or not
+ * @return Exitflag for error (0 if no errors)
+ */
+c_int init_linsys_solver_pardiso(pardiso_solver ** sp, const csc * P, const csc * A, c_float sigma, const c_float * rho_vec, c_int polish);
+
+
+/**
+ * Solve linear system and store result in b
+ * @param s Linear system solver structure
+ * @param b Right-hand side
+ * @return Exitflag
+ */
+c_int solve_linsys_pardiso(pardiso_solver * s, c_float * b);
+
+
+/**
+ * Update linear system solver matrices
+ * @param s Linear system solver structure
+ * @param P Matrix P
+ * @param A Matrix A
+ * @return Exitflag
+ */
+c_int update_linsys_solver_matrices_pardiso(pardiso_solver * s, const csc *P, const csc *A);
+
+
+/**
+ * Update rho parameter in linear system solver structure
+ * @param s Linear system solver structure
+ * @param rho_vec new rho_vec value
+ * @return exitflag
+ */
+c_int update_linsys_solver_rho_vec_pardiso(pardiso_solver * s, const c_float * rho_vec);
+
+
+/**
+ * Free linear system solver
+ * @param s linear system solver object
+ */
+void free_linsys_solver_pardiso(pardiso_solver * s);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/lin_sys/direct/pardiso/pardiso_loader.c b/lin_sys/direct/pardiso/pardiso_loader.c
new file mode 100644
index 0000000..f120e47
--- /dev/null
+++ b/lin_sys/direct/pardiso/pardiso_loader.c
@@ -0,0 +1,102 @@
+#include "lib_handler.h"
+#include "pardiso_loader.h"
+
+#include "glob_opts.h"
+#include "constants.h"
+
+#ifdef IS_WINDOWS
+#define PARDISOLIBNAME "mkl_rt." SHAREDLIBEXT
+#else
+#define PARDISOLIBNAME "libmkl_rt." SHAREDLIBEXT
+#endif
+
+typedef void (*voidfun)(void);
+
+voidfun lh_load_sym (soHandle_t h, const char *symName);
+
+
+// Interfaces for Pardiso functions
+typedef void (*pardiso_t)(void**, const c_int*, const c_int*, const c_int*,
+ const c_int*, const c_int*, const c_float*,
+ const c_int*, const c_int*, c_int*,
+ const c_int*, c_int*, const c_int*, c_float*,
+ c_float*, c_int*);
+typedef int (*mkl_set_ifl_t)(int);
+typedef int (*mkl_get_mt_t)();
+
+
+// Handlers are static variables
+static soHandle_t Pardiso_handle = OSQP_NULL;
+static pardiso_t func_pardiso = OSQP_NULL;
+static mkl_set_ifl_t func_mkl_set_interface_layer = OSQP_NULL;
+static mkl_get_mt_t func_mkl_get_max_threads = OSQP_NULL;
+
+// Wrappers for loaded Pardiso function handlers
+void pardiso(void** pt, const c_int* maxfct, const c_int* mnum,
+ const c_int* mtype, const c_int* phase, const c_int* n,
+ const c_float* a, const c_int* ia, const c_int* ja,
+ c_int* perm, const c_int* nrhs, c_int* iparm,
+ const c_int* msglvl, c_float* b, c_float* x,
+ c_int* error) {
+ if(func_pardiso){
+ // Call function pardiso only if it has been initialized
+ func_pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
+ perm, nrhs, iparm, msglvl, b, x, error);
+ }
+ else
+ {
+#ifdef PRINTING
+ c_eprint("Pardiso not loaded correctly");
+#endif
+ }
+}
+
+c_int mkl_set_interface_layer(c_int code) {
+ return (c_int)func_mkl_set_interface_layer((int)code);
+}
+
+c_int mkl_get_max_threads() {
+ return (c_int)func_mkl_get_max_threads();
+}
+
+
+c_int lh_load_pardiso(const char* libname) {
+ // DEBUG
+ // if (Pardiso_handle) return 0;
+
+ // Load Pardiso library
+ if (libname) {
+ Pardiso_handle = lh_load_lib(libname);
+ } else { /* try a default library name */
+ Pardiso_handle = lh_load_lib(PARDISOLIBNAME);
+ }
+ if (!Pardiso_handle) return 1;
+
+ // Load Pardiso functions
+ func_pardiso = (pardiso_t)lh_load_sym(Pardiso_handle, "pardiso");
+ if (!func_pardiso) return 1;
+
+ func_mkl_set_interface_layer = (mkl_set_ifl_t)lh_load_sym(Pardiso_handle,
+ "MKL_Set_Interface_Layer");
+ if (!func_mkl_set_interface_layer) return 1;
+
+ func_mkl_get_max_threads = (mkl_get_mt_t)lh_load_sym(Pardiso_handle,
+ "MKL_Get_Max_Threads");
+ if (!func_mkl_get_max_threads) return 1;
+
+ return 0;
+}
+
+c_int lh_unload_pardiso() {
+
+ if (Pardiso_handle == OSQP_NULL) return 0;
+
+ return lh_unload_lib(Pardiso_handle);
+
+ /* If multiple OSQP objects are laoded, the lines below cause a crash */
+ // Pardiso_handle = OSQP_NULL;
+ // func_pardiso = OSQP_NULL;
+ // func_mkl_set_interface_layer = OSQP_NULL;
+ // func_mkl_get_max_threads = OSQP_NULL;
+
+}
diff --git a/lin_sys/direct/pardiso/pardiso_loader.h b/lin_sys/direct/pardiso/pardiso_loader.h
new file mode 100644
index 0000000..aa96b24
--- /dev/null
+++ b/lin_sys/direct/pardiso/pardiso_loader.h
@@ -0,0 +1,29 @@
+#ifndef PARDISOLOADER_H
+#define PARDISOLOADER_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+
+/**
+ * Tries to load a shared library with Pardiso.
+ * Return a failure if the library cannot be loaded or not all Pardiso symbols are found.
+ * @param libname The name under which the Pardiso lib can be found, or OSQP_NULL to use a default name (mkl_rt.SHAREDLIBEXT).
+ * @return Zero on success, nonzero on failure.
+ */
+c_int lh_load_pardiso(const char* libname);
+
+/**
+ * Unloads the loaded Pardiso shared library.
+ * @return Zero on success, nonzero on failure.
+ */
+c_int lh_unload_pardiso();
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /*PARADISOLOADER_H*/
diff --git a/lin_sys/direct/qdldl/CMakeLists.txt b/lin_sys/direct/qdldl/CMakeLists.txt
new file mode 100644
index 0000000..d78c4a3
--- /dev/null
+++ b/lin_sys/direct/qdldl/CMakeLists.txt
@@ -0,0 +1,41 @@
+# Add qdldl
+add_subdirectory(qdldl_sources)
+
+
+if(NOT DEFINED EMBEDDED)
+set(
+ amd_sources
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/include/amd_internal.h
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/include/amd.h
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/include/SuiteSparse_config.h
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_1.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_2.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_aat.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_control.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_defaults.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_info.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_order.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_post_tree.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_postorder.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_preprocess.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/amd_valid.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/src/SuiteSparse_config.c
+)
+endif()
+
+
+set(qdldl_interface_includes
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ ${CMAKE_CURRENT_SOURCE_DIR}/amd/include
+ ${CMAKE_CURRENT_SOURCE_DIR}/qdldl_sources/include
+)
+
+set(qdldl_interface_src
+ ${amd_sources}
+ ${CMAKE_CURRENT_SOURCE_DIR}/qdldl_interface.h
+ ${CMAKE_CURRENT_SOURCE_DIR}/qdldl_interface.c
+)
+
+# Create object library for linear system solver interface
+add_library(linsys_qdldl OBJECT ${qdldl_interface_src})
+target_include_directories(linsys_qdldl PRIVATE ${qdldl_interface_includes} ${PROJECT_SOURCE_DIR}/include)
diff --git a/lin_sys/direct/qdldl/amd/LICENSE b/lin_sys/direct/qdldl/amd/LICENSE
new file mode 100644
index 0000000..36de6d8
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/LICENSE
@@ -0,0 +1,36 @@
+AMD, Copyright (c), 1996-2015, Timothy A. Davis,
+Patrick R. Amestoy, and Iain S. Duff. All Rights Reserved.
+
+Availability:
+
+ http://www.suitesparse.com
+
+-------------------------------------------------------------------------------
+AMD License: BSD 3-clause:
+-------------------------------------------------------------------------------
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions are met:
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+ * Neither the name of the organizations to which the authors are
+ affiliated, nor the names of its contributors may be used to endorse
+ or promote products derived from this software without specific prior
+ written permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
+ DAMAGE.
+
+
diff --git a/lin_sys/direct/qdldl/amd/include/SuiteSparse_config.h b/lin_sys/direct/qdldl/amd/include/SuiteSparse_config.h
new file mode 100644
index 0000000..abeb7b9
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/include/SuiteSparse_config.h
@@ -0,0 +1,246 @@
+/* ========================================================================== */
+/* === SuiteSparse_config =================================================== */
+/* ========================================================================== */
+
+/* Configuration file for SuiteSparse: a Suite of Sparse matrix packages
+ * (AMD, COLAMD, CCOLAMD, CAMD, CHOLMOD, UMFPACK, CXSparse, and others).
+ *
+ * SuiteSparse_config.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 SuiteSparse_config.mk:
+ *
+ * CFLAGS = -O -D'SuiteSparse_long=long long' \
+ * -D'SuiteSparse_long_max=9223372036854775801' -D'SuiteSparse_long_idd="lld"'
+ *
+ * This file defines SuiteSparse_long as either long (on all but _WIN64) or
+ * __int64 on Windows 64. The intent is that a SuiteSparse_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) 2012, Timothy A. Davis. No licensing restrictions apply
+ * to this file or to the SuiteSparse_config directory.
+ * Author: Timothy A. Davis.
+ */
+
+#ifndef SUITESPARSE_CONFIG_H
+#define SUITESPARSE_CONFIG_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include "glob_opts.h"
+#include <limits.h>
+#include <stdlib.h>
+
+/* ========================================================================== */
+/* === SuiteSparse_long ===================================================== */
+/* ========================================================================== */
+
+#ifndef SuiteSparse_long
+
+#define SuiteSparse_long long long
+#define SuiteSparse_long_max LONG_MAX
+#define SuiteSparse_long_idd "ld"
+
+#define SuiteSparse_long_id "%" SuiteSparse_long_idd
+#endif
+
+/* ========================================================================== */
+/* === SuiteSparse_config parameters and functions ========================== */
+/* ========================================================================== */
+
+/* SuiteSparse-wide parameters are placed in this struct. It is meant to be
+ an extern, globally-accessible struct. It is not meant to be updated
+ frequently by multiple threads. Rather, if an application needs to modify
+ SuiteSparse_config, it should do it once at the beginning of the application,
+ before multiple threads are launched.
+
+ The intent of these function pointers is that they not be used in your
+ application directly, except to assign them to the desired user-provided
+ functions. Rather, you should use the
+ */
+
+struct SuiteSparse_config_struct
+{
+ void *(*malloc_func) (size_t) ; /* pointer to malloc */
+ void *(*realloc_func) (void *, size_t) ; /* pointer to realloc */
+ void (*free_func) (void *) ; /* pointer to free */
+#if defined PYTHON
+ void (*printf_func) (const char *, ...) ; /* pointer to printf (in Python it returns void)*/
+#elif defined R_LANG
+ void (*printf_func) (const char *, ...) ; /* pointer to printf (in R it returns void)*/
+#else
+ int (*printf_func) (const char *, ...) ; /* pointer to printf */
+#endif
+ c_float (*hypot_func) (c_float, c_float) ; /* pointer to hypot */
+ int (*divcomplex_func) (c_float, c_float, c_float, c_float, c_float *, c_float *);
+} ;
+
+extern struct SuiteSparse_config_struct SuiteSparse_config ;
+
+void *SuiteSparse_malloc /* 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 */
+) ;
+
+void *SuiteSparse_calloc /* pointer to allocated block of memory */
+(
+ size_t nitems, /* number of items to calloc (>=1 is enforced) */
+ size_t size_of_item /* sizeof each item */
+) ;
+
+void *SuiteSparse_realloc /* pointer to reallocated block of memory, or
+ to original block if the realloc failed. */
+(
+ size_t nitems_new, /* new number of items in the object */
+ size_t nitems_old, /* old number of items in the object */
+ size_t size_of_item, /* sizeof each item */
+ void *p, /* old object to reallocate */
+ int *ok /* 1 if successful, 0 otherwise */
+) ;
+
+void *SuiteSparse_free /* always returns NULL */
+(
+ void *p /* block to free */
+) ;
+
+void SuiteSparse_tic /* start the timer */
+(
+ c_float tic [2] /* output, contents undefined on input */
+) ;
+
+c_float SuiteSparse_toc /* return time in seconds since last tic */
+(
+ c_float tic [2] /* input: from last call to SuiteSparse_tic */
+) ;
+
+c_float SuiteSparse_time /* returns current wall clock time in seconds */
+(
+ void
+) ;
+
+/* returns sqrt (x^2 + y^2), computed reliably */
+c_float SuiteSparse_hypot (c_float x, c_float y) ;
+
+/* complex division of c = a/b */
+int SuiteSparse_divcomplex
+(
+ c_float ar, c_float ai, /* real and imaginary parts of a */
+ c_float br, c_float bi, /* real and imaginary parts of b */
+ c_float *cr, c_float *ci /* real and imaginary parts of c */
+) ;
+
+/* OSQP: disabling this timing code */
+#define NTIMER
+
+/* OSQP: disabling Suitesparse printing */
+#define NPRINT
+
+/* determine which timer to use, if any */
+#ifndef NTIMER
+#ifdef _POSIX_C_SOURCE
+#if _POSIX_C_SOURCE >= 199309L
+#define SUITESPARSE_TIMER_ENABLED
+#endif
+#endif
+#endif
+
+/* SuiteSparse printf macro */
+#define SUITESPARSE_PRINTF(params) \
+{ \
+ if (SuiteSparse_config.printf_func != NULL) \
+ { \
+ (void) (SuiteSparse_config.printf_func) params ; \
+ } \
+}
+
+/* ========================================================================== */
+/* === 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 package from one
+ * version of SuiteSparse, with another package from another version of
+ * SuiteSparse, may or may not work.
+ *
+ * SuiteSparse contains the following packages:
+ *
+ * SuiteSparse_config version 4.5.3 (version always the same as SuiteSparse)
+ * AMD version 2.4.6
+ * BTF version 1.2.6
+ * CAMD version 2.4.6
+ * CCOLAMD version 2.9.6
+ * CHOLMOD version 3.0.11
+ * COLAMD version 2.9.6
+ * CSparse version 3.1.9
+ * CXSparse version 3.1.9
+ * GPUQREngine version 1.0.5
+ * KLU version 1.3.8
+ * LDL version 2.2.6
+ * RBio version 2.2.6
+ * SPQR version 2.0.7
+ * SuiteSparse_GPURuntime version 1.0.5
+ * UMFPACK version 5.7.6
+ * MATLAB_Tools various packages & M-files
+ * xerbla version 1.0.3
+ *
+ * Other package dependencies:
+ * BLAS required by CHOLMOD and UMFPACK
+ * LAPACK required by CHOLMOD
+ * METIS 5.1.0 required by CHOLMOD (optional) and KLU (optional)
+ * CUBLAS, CUDART NVIDIA libraries required by CHOLMOD and SPQR when
+ * they are compiled with GPU acceleration.
+ */
+
+int SuiteSparse_version /* returns SUITESPARSE_VERSION */
+(
+ /* output, not defined on input. Not used if NULL. Returns
+ the three version codes in version [0..2]:
+ version [0] is SUITESPARSE_MAIN_VERSION
+ version [1] is SUITESPARSE_SUB_VERSION
+ version [2] is SUITESPARSE_SUBSUB_VERSION
+ */
+ int version [3]
+) ;
+
+/* Versions prior to 4.2.0 do not have the above function. The following
+ code fragment will work with any version of SuiteSparse:
+
+ #ifdef SUITESPARSE_HAS_VERSION_FUNCTION
+ v = SuiteSparse_version (NULL) ;
+ #else
+ v = SUITESPARSE_VERSION ;
+ #endif
+*/
+#define SUITESPARSE_HAS_VERSION_FUNCTION
+
+#define SUITESPARSE_DATE "May 4, 2016"
+#define SUITESPARSE_VER_CODE(main,sub) ((main) * 1000 + (sub))
+#define SUITESPARSE_MAIN_VERSION 4
+#define SUITESPARSE_SUB_VERSION 5
+#define SUITESPARSE_SUBSUB_VERSION 3
+#define SUITESPARSE_VERSION \
+ SUITESPARSE_VER_CODE(SUITESPARSE_MAIN_VERSION,SUITESPARSE_SUB_VERSION)
+
+#ifdef __cplusplus
+}
+#endif
+#endif
diff --git a/lin_sys/direct/qdldl/amd/include/amd.h b/lin_sys/direct/qdldl/amd/include/amd.h
new file mode 100644
index 0000000..cd81e6e
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/include/amd.h
@@ -0,0 +1,401 @@
+/* ========================================================================= */
+/* === AMD: approximate minimum degree ordering =========================== */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD Version 2.4, Copyright (c) 1996-2013 by Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* 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 SuiteSparse_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
+
+#include "SuiteSparse_config.h"
+
+/* get the definition of size_t: */
+#include <stddef.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 */
+ c_float Control [ ], /* input Control settings, of size AMD_CONTROL */
+ c_float Info [ ] /* output Info statistics, of size AMD_INFO */
+) ;
+
+SuiteSparse_long amd_l_order /* see above for description of arguments */
+(
+ SuiteSparse_long n,
+ const SuiteSparse_long Ap [ ],
+ const SuiteSparse_long Ai [ ],
+ SuiteSparse_long P [ ],
+ c_float Control [ ],
+ c_float Info [ ]
+) ;
+
+/* Input arguments (not modified):
+ *
+ * n: the matrix A is n-by-n.
+ * Ap: an int/SuiteSparse_long array of size n+1, containing column
+ * pointers of A.
+ * Ai: an int/SuiteSparse_long array of size nz, containing the row
+ * indices of A, where nz = Ap [n].
+ * Control: a c_float array of size AMD_CONTROL, containing control
+ * parameters. Defaults are used if Control is NULL.
+ *
+ * Output arguments (not defined on input):
+ *
+ * P: an int/SuiteSparse_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 c_float 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.suitesparse.com.
+ *
+ * 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 [ ],
+ c_float Control [ ],
+ c_float Info [ ]
+) ;
+
+void amd_l2
+(
+ SuiteSparse_long n,
+ SuiteSparse_long Pe [ ],
+ SuiteSparse_long Iw [ ],
+ SuiteSparse_long Len [ ],
+ SuiteSparse_long iwlen,
+ SuiteSparse_long pfree,
+ SuiteSparse_long Nv [ ],
+ SuiteSparse_long Next [ ],
+ SuiteSparse_long Last [ ],
+ SuiteSparse_long Head [ ],
+ SuiteSparse_long Elen [ ],
+ SuiteSparse_long Degree [ ],
+ SuiteSparse_long W [ ],
+ c_float Control [ ],
+ c_float 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] */
+) ;
+
+SuiteSparse_long amd_l_valid
+(
+ SuiteSparse_long n_row,
+ SuiteSparse_long n_col,
+ const SuiteSparse_long Ap [ ],
+ const SuiteSparse_long Ai [ ]
+) ;
+
+/* ------------------------------------------------------------------------- */
+/* AMD memory manager and printf routines */
+/* ------------------------------------------------------------------------- */
+
+ /* moved to SuiteSparse_config.c */
+
+/* ------------------------------------------------------------------------- */
+/* AMD Control and Info arrays */
+/* ------------------------------------------------------------------------- */
+
+/* amd_defaults: sets the default control settings */
+void amd_defaults (c_float Control [ ]) ;
+void amd_l_defaults (c_float Control [ ]) ;
+
+/* amd_control: prints the control settings */
+void amd_control (c_float Control [ ]) ;
+void amd_l_control (c_float Control [ ]) ;
+
+/* amd_info: prints the statistics */
+void amd_info (c_float Info [ ]) ;
+void amd_l_info (c_float 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 "May 4, 2016"
+#define AMD_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
+#define AMD_MAIN_VERSION 2
+#define AMD_SUB_VERSION 4
+#define AMD_SUBSUB_VERSION 6
+#define AMD_VERSION AMD_VERSION_CODE(AMD_MAIN_VERSION,AMD_SUB_VERSION)
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/lin_sys/direct/qdldl/amd/include/amd_internal.h b/lin_sys/direct/qdldl/amd/include/amd_internal.h
new file mode 100644
index 0000000..7b7cea1
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/include/amd_internal.h
@@ -0,0 +1,328 @@
+/* ========================================================================= */
+/* === amd_internal.h ====================================================== */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* This file is for internal use in AMD itself, and does not normally need to
+ * be included in user code (it is included in UMFPACK, however). All others
+ * should use amd.h instead.
+ */
+
+/* ========================================================================= */
+/* === NDEBUG ============================================================== */
+/* ========================================================================= */
+
+/*
+ * Turning on debugging takes some work (see below). If you do not edit this
+ * file, then debugging is always turned off, regardless of whether or not
+ * -DNDEBUG is specified in your compiler options.
+ *
+ * If AMD is being compiled as a mexFunction, then MATLAB is defined,
+ * and mxAssert is used instead of assert. If debugging is not enabled, no
+ * MATLAB include files or functions are used. Thus, the AMD library libamd.a
+ * can be safely used in either a stand-alone C program or in another
+ * mexFunction, without any change.
+ */
+
+/*
+ AMD will be exceedingly slow when running in debug mode. The next three
+ lines ensure that debugging is turned off.
+*/
+#ifndef NDEBUG
+#define NDEBUG
+#endif
+
+#include "glob_opts.h"
+
+/*
+ To enable debugging, uncomment the following line:
+#undef NDEBUG
+*/
+
+/* ------------------------------------------------------------------------- */
+/* ANSI include files */
+/* ------------------------------------------------------------------------- */
+
+/* from stdlib.h: size_t, malloc, free, realloc, and calloc */
+#include <stdlib.h>
+
+#if !defined(NPRINT) || !defined(NDEBUG)
+/* from stdio.h: printf. Not included if NPRINT is defined at compile time.
+ * fopen and fscanf are used when debugging. */
+#include <stdio.h>
+#endif
+
+/* from limits.h: INT_MAX and LONG_MAX */
+#include <limits.h>
+
+/* from math.h: sqrt */
+#include <math.h>
+
+/* ------------------------------------------------------------------------- */
+/* MATLAB include files (only if being used in or via MATLAB) */
+/* ------------------------------------------------------------------------- */
+
+#ifdef MATLAB
+#include "matrix.h"
+#include "mex.h"
+#endif
+
+/* ------------------------------------------------------------------------- */
+/* basic definitions */
+/* ------------------------------------------------------------------------- */
+
+#ifdef FLIP
+#undef FLIP
+#endif
+
+#ifdef MAX
+#undef MAX
+#endif
+
+#ifdef MIN
+#undef MIN
+#endif
+
+#ifdef EMPTY
+#undef EMPTY
+#endif
+
+#ifdef GLOBAL
+#undef GLOBAL
+#endif
+
+#ifdef PRIVATE
+#undef PRIVATE
+#endif
+
+/* FLIP is a "negation about -1", and is used to mark an integer i that is
+ * normally non-negative. FLIP (EMPTY) is EMPTY. FLIP of a number > EMPTY
+ * is negative, and FLIP of a number < EMTPY is positive. FLIP (FLIP (i)) = i
+ * for all integers i. UNFLIP (i) is >= EMPTY. */
+#define EMPTY (-1)
+#define FLIP(i) (-(i)-2)
+#define UNFLIP(i) ((i < EMPTY) ? FLIP (i) : (i))
+
+/* for integer MAX/MIN, or for c_floats when we don't care how NaN's behave: */
+#define MAX(a,b) (((a) > (b)) ? (a) : (b))
+#define MIN(a,b) (((a) < (b)) ? (a) : (b))
+
+/* logical expression of p implies q: */
+#define IMPLIES(p,q) (!(p) || (q))
+
+/* Note that the IBM RS 6000 xlc predefines TRUE and FALSE in <types.h>. */
+/* The Compaq Alpha also predefines TRUE and FALSE. */
+#ifdef TRUE
+#undef TRUE
+#endif
+#ifdef FALSE
+#undef FALSE
+#endif
+
+#define TRUE (1)
+#define FALSE (0)
+#define PRIVATE static
+#define GLOBAL
+#define EMPTY (-1)
+
+/* Note that Linux's gcc 2.96 defines NULL as ((void *) 0), but other */
+/* compilers (even gcc 2.95.2 on Solaris) define NULL as 0 or (0). We */
+/* need to use the ANSI standard value of 0. */
+#ifdef NULL
+#undef NULL
+#endif
+
+#define NULL 0
+
+/* largest value of size_t */
+#ifndef SIZE_T_MAX
+#ifdef SIZE_MAX
+/* C99 only */
+#define SIZE_T_MAX SIZE_MAX
+#else
+#define SIZE_T_MAX ((size_t) (-1))
+#endif
+#endif
+
+/* ------------------------------------------------------------------------- */
+/* integer type for AMD: int or SuiteSparse_long */
+/* ------------------------------------------------------------------------- */
+
+#include "amd.h"
+
+#if defined (DLONG) || defined (ZLONG)
+
+#define Int SuiteSparse_long
+#define ID SuiteSparse_long_id
+#define Int_MAX SuiteSparse_long_max
+
+#define AMD_order amd_l_order
+#define AMD_defaults amd_l_defaults
+#define AMD_control amd_l_control
+#define AMD_info amd_l_info
+#define AMD_1 amd_l1
+#define AMD_2 amd_l2
+#define AMD_valid amd_l_valid
+#define AMD_aat amd_l_aat
+#define AMD_postorder amd_l_postorder
+#define AMD_post_tree amd_l_post_tree
+#define AMD_dump amd_l_dump
+#define AMD_debug amd_l_debug
+#define AMD_debug_init amd_l_debug_init
+#define AMD_preprocess amd_l_preprocess
+
+#else
+
+#define Int int
+#define ID "%d"
+#define Int_MAX INT_MAX
+
+#define AMD_order amd_order
+#define AMD_defaults amd_defaults
+#define AMD_control amd_control
+#define AMD_info amd_info
+#define AMD_1 amd_1
+#define AMD_2 amd_2
+#define AMD_valid amd_valid
+#define AMD_aat amd_aat
+#define AMD_postorder amd_postorder
+#define AMD_post_tree amd_post_tree
+#define AMD_dump amd_dump
+#define AMD_debug amd_debug
+#define AMD_debug_init amd_debug_init
+#define AMD_preprocess amd_preprocess
+
+#endif
+
+/* ------------------------------------------------------------------------- */
+/* AMD routine definitions (not user-callable) */
+/* ------------------------------------------------------------------------- */
+
+GLOBAL size_t AMD_aat
+(
+ Int n,
+ const Int Ap [ ],
+ const Int Ai [ ],
+ Int Len [ ],
+ Int Tp [ ],
+ c_float Info [ ]
+) ;
+
+GLOBAL void AMD_1
+(
+ Int n,
+ const Int Ap [ ],
+ const Int Ai [ ],
+ Int P [ ],
+ Int Pinv [ ],
+ Int Len [ ],
+ Int slen,
+ Int S [ ],
+ c_float Control [ ],
+ c_float Info [ ]
+) ;
+
+GLOBAL void AMD_postorder
+(
+ Int nn,
+ Int Parent [ ],
+ Int Npiv [ ],
+ Int Fsize [ ],
+ Int Order [ ],
+ Int Child [ ],
+ Int Sibling [ ],
+ Int Stack [ ]
+) ;
+
+GLOBAL Int AMD_post_tree
+(
+ Int root,
+ Int k,
+ Int Child [ ],
+ const Int Sibling [ ],
+ Int Order [ ],
+ Int Stack [ ]
+#ifndef NDEBUG
+ , Int nn
+#endif
+) ;
+
+GLOBAL void AMD_preprocess
+(
+ Int n,
+ const Int Ap [ ],
+ const Int Ai [ ],
+ Int Rp [ ],
+ Int Ri [ ],
+ Int W [ ],
+ Int Flag [ ]
+) ;
+
+/* ------------------------------------------------------------------------- */
+/* debugging definitions */
+/* ------------------------------------------------------------------------- */
+
+#ifndef NDEBUG
+
+/* from assert.h: assert macro */
+#include <assert.h>
+
+#ifndef EXTERN
+#define EXTERN extern
+#endif
+
+EXTERN Int AMD_debug ;
+
+GLOBAL void AMD_debug_init ( char *s ) ;
+
+GLOBAL void AMD_dump
+(
+ 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 [ ],
+ Int nel
+) ;
+
+#ifdef ASSERT
+#undef ASSERT
+#endif
+
+/* Use mxAssert if AMD is compiled into a mexFunction */
+#ifdef MATLAB
+#define ASSERT(expression) (mxAssert ((expression), ""))
+#else
+#define ASSERT(expression) (assert (expression))
+#endif
+
+#define AMD_DEBUG0(params) { SUITESPARSE_PRINTF (params) ; }
+#define AMD_DEBUG1(params) { if (AMD_debug >= 1) SUITESPARSE_PRINTF (params) ; }
+#define AMD_DEBUG2(params) { if (AMD_debug >= 2) SUITESPARSE_PRINTF (params) ; }
+#define AMD_DEBUG3(params) { if (AMD_debug >= 3) SUITESPARSE_PRINTF (params) ; }
+#define AMD_DEBUG4(params) { if (AMD_debug >= 4) SUITESPARSE_PRINTF (params) ; }
+
+#else
+
+/* no debugging */
+#define ASSERT(expression)
+#define AMD_DEBUG0(params)
+#define AMD_DEBUG1(params)
+#define AMD_DEBUG2(params)
+#define AMD_DEBUG3(params)
+#define AMD_DEBUG4(params)
+
+#endif
diff --git a/lin_sys/direct/qdldl/amd/src/SuiteSparse_config.c b/lin_sys/direct/qdldl/amd/src/SuiteSparse_config.c
new file mode 100644
index 0000000..cebdff0
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/SuiteSparse_config.c
@@ -0,0 +1,407 @@
+/* ========================================================================== */
+/* === SuiteSparse_config =================================================== */
+/* ========================================================================== */
+
+/* SuiteSparse configuration : memory manager and printf functions. */
+
+/* Copyright (c) 2013, Timothy A. Davis. No licensing restrictions
+ * apply to this file or to the SuiteSparse_config directory.
+ * Author: Timothy A. Davis.
+ */
+
+// Include OSQP Global options for memory management
+#include "glob_opts.h"
+
+#include <math.h>
+#include <stdlib.h>
+
+
+
+#ifndef NPRINT
+#include <stdio.h>
+#endif
+
+#ifdef MATLAB
+#include "mex.h"
+#include "matrix.h"
+#endif
+
+#ifndef NULL
+#define NULL ((void *) 0)
+#endif
+
+#include "SuiteSparse_config.h"
+
+/* -------------------------------------------------------------------------- */
+/* SuiteSparse_config : a global extern struct */
+/* -------------------------------------------------------------------------- */
+
+/* The SuiteSparse_config struct is available to all SuiteSparse functions and
+ to all applications that use those functions. It must be modified with
+ care, particularly in a multithreaded context. Normally, the application
+ will initialize this object once, via SuiteSparse_start, possibily followed
+ by application-specific modifications if the applications wants to use
+ alternative memory manager functions.
+
+ The user can redefine these global pointers at run-time to change the
+ memory manager and printf function used by SuiteSparse.
+
+ If -DNMALLOC is defined at compile-time, then no memory-manager is
+ specified. You must define them at run-time, after calling
+ SuiteSparse_start.
+
+ If -DPRINT is defined a compile time, then printf is disabled, and
+ SuiteSparse will not use printf.
+ */
+
+struct SuiteSparse_config_struct SuiteSparse_config =
+{
+ // Memory allocation from glob_opts.h in OSQP
+ c_malloc, c_realloc, c_free,
+ // Ignore printing
+ NULL,
+ SuiteSparse_hypot,
+ SuiteSparse_divcomplex
+
+} ;
+
+
+/* -------------------------------------------------------------------------- */
+/* SuiteSparse_malloc: malloc wrapper */
+/* -------------------------------------------------------------------------- */
+
+void *SuiteSparse_malloc /* pointer to allocated block of memory */
+(
+ size_t nitems, /* number of items to malloc */
+ size_t size_of_item /* sizeof each item */
+)
+{
+ void *p ;
+ size_t size ;
+ if (nitems < 1) nitems = 1 ;
+ if (size_of_item < 1) size_of_item = 1 ;
+ size = nitems * size_of_item ;
+
+ if (size != ((c_float) nitems) * size_of_item)
+ {
+ /* size_t overflow */
+ p = NULL ;
+ }
+ else
+ {
+ p = (void *) (SuiteSparse_config.malloc_func) (size) ;
+ // p = (void *) c_malloc(size) ;
+ }
+ return (p) ;
+}
+
+
+
+/* -------------------------------------------------------------------------- */
+/* SuiteSparse_realloc: realloc wrapper */
+/* -------------------------------------------------------------------------- */
+
+/* If p is non-NULL on input, it points to a previously allocated object of
+ size nitems_old * size_of_item. The object is reallocated to be of size
+ nitems_new * size_of_item. If p is NULL on input, then a new object of that
+ size is allocated. On success, a pointer to the new object is returned,
+ and ok is returned as 1. If the allocation fails, ok is set to 0 and a
+ pointer to the old (unmodified) object is returned.
+ */
+
+void *SuiteSparse_realloc /* pointer to reallocated block of memory, or
+ to original block if the realloc failed. */
+(
+ size_t nitems_new, /* new number of items in the object */
+ size_t nitems_old, /* old number of items in the object */
+ size_t size_of_item, /* sizeof each item */
+ void *p, /* old object to reallocate */
+ int *ok /* 1 if successful, 0 otherwise */
+)
+{
+ size_t size ;
+ if (nitems_old < 1) nitems_old = 1 ;
+ if (nitems_new < 1) nitems_new = 1 ;
+ if (size_of_item < 1) size_of_item = 1 ;
+ size = nitems_new * size_of_item ;
+
+ if (size != ((c_float) nitems_new) * size_of_item)
+ {
+ /* size_t overflow */
+ (*ok) = 0 ;
+ }
+ else if (p == NULL)
+ {
+ /* a fresh object is being allocated */
+ p = SuiteSparse_malloc (nitems_new, size_of_item) ;
+ (*ok) = (p != NULL) ;
+ }
+ else if (nitems_old == nitems_new)
+ {
+ /* the object does not change; do nothing */
+ (*ok) = 1 ;
+ }
+ else
+ {
+ /* change the size of the object from nitems_old to nitems_new */
+ void *pnew ;
+ // pnew = (void *) (SuiteSparse_config.realloc_func) (p, size) ;
+ pnew = (void *) c_realloc (p, size) ;
+ if (pnew == NULL)
+ {
+ if (nitems_new < nitems_old)
+ {
+ /* the attempt to reduce the size of the block failed, but
+ the old block is unchanged. So pretend to succeed. */
+ (*ok) = 1 ;
+ }
+ else
+ {
+ /* out of memory */
+ (*ok) = 0 ;
+ }
+ }
+ else
+ {
+ /* success */
+ p = pnew ;
+ (*ok) = 1 ;
+ }
+ }
+ return (p) ;
+}
+
+/* -------------------------------------------------------------------------- */
+/* SuiteSparse_free: free wrapper */
+/* -------------------------------------------------------------------------- */
+
+void *SuiteSparse_free /* always returns NULL */
+(
+ void *p /* block to free */
+)
+{
+ if (p)
+ {
+ (SuiteSparse_config.free_func) (p) ;
+ // c_free(p) ;
+ }
+ return (NULL) ;
+}
+
+
+/* -------------------------------------------------------------------------- */
+/* SuiteSparse_tic: return current wall clock time */
+/* -------------------------------------------------------------------------- */
+
+/* Returns the number of seconds (tic [0]) and nanoseconds (tic [1]) since some
+ * unspecified but fixed time in the past. If no timer is installed, zero is
+ * returned. A scalar c_float precision value for 'tic' could be used, but this
+ * might cause loss of precision because clock_getttime returns the time from
+ * some distant time in the past. Thus, an array of size 2 is used.
+ *
+ * The timer is enabled by default. To disable the timer, compile with
+ * -DNTIMER. If enabled on a POSIX C 1993 system, the timer requires linking
+ * with the -lrt library.
+ *
+ * example:
+ *
+ * c_float tic [2], r, s, t ;
+ * SuiteSparse_tic (tic) ; // start the timer
+ * // do some work A
+ * t = SuiteSparse_toc (tic) ; // t is time for work A, in seconds
+ * // do some work B
+ * s = SuiteSparse_toc (tic) ; // s is time for work A and B, in seconds
+ * SuiteSparse_tic (tic) ; // restart the timer
+ * // do some work C
+ * r = SuiteSparse_toc (tic) ; // s is time for work C, in seconds
+ *
+ * A c_float array of size 2 is used so that this routine can be more easily
+ * ported to non-POSIX systems. The caller does not rely on the POSIX
+ * <time.h> include file.
+ */
+
+#ifdef SUITESPARSE_TIMER_ENABLED
+
+#include <time.h>
+
+void SuiteSparse_tic
+(
+ c_float tic [2] /* output, contents undefined on input */
+)
+{
+ /* POSIX C 1993 timer, requires -librt */
+ struct timespec t ;
+ clock_gettime (CLOCK_MONOTONIC, &t) ;
+ tic [0] = (c_float) (t.tv_sec) ;
+ tic [1] = (c_float) (t.tv_nsec) ;
+}
+
+#else
+
+void SuiteSparse_tic
+(
+ c_float tic [2] /* output, contents undefined on input */
+)
+{
+ /* no timer installed */
+ tic [0] = 0 ;
+ tic [1] = 0 ;
+}
+
+#endif
+
+
+/* -------------------------------------------------------------------------- */
+/* SuiteSparse_toc: return time since last tic */
+/* -------------------------------------------------------------------------- */
+
+/* Assuming SuiteSparse_tic is accurate to the nanosecond, this function is
+ * accurate down to the nanosecond for 2^53 nanoseconds since the last call to
+ * SuiteSparse_tic, which is sufficient for SuiteSparse (about 104 days). If
+ * additional accuracy is required, the caller can use two calls to
+ * SuiteSparse_tic and do the calculations differently.
+ */
+
+c_float SuiteSparse_toc /* returns time in seconds since last tic */
+(
+ c_float tic [2] /* input, not modified from last call to SuiteSparse_tic */
+)
+{
+ c_float toc [2] ;
+ SuiteSparse_tic (toc) ;
+ return ((toc [0] - tic [0]) + 1e-9 * (toc [1] - tic [1])) ;
+}
+
+
+/* -------------------------------------------------------------------------- */
+/* SuiteSparse_time: return current wallclock time in seconds */
+/* -------------------------------------------------------------------------- */
+
+/* This function might not be accurate down to the nanosecond. */
+
+c_float SuiteSparse_time /* returns current wall clock time in seconds */
+(
+ void
+)
+{
+ c_float toc [2] ;
+ SuiteSparse_tic (toc) ;
+ return (toc [0] + 1e-9 * toc [1]) ;
+}
+
+
+/* -------------------------------------------------------------------------- */
+/* SuiteSparse_version: return the current version of SuiteSparse */
+/* -------------------------------------------------------------------------- */
+
+int SuiteSparse_version
+(
+ int version [3]
+)
+{
+ if (version != NULL)
+ {
+ version [0] = SUITESPARSE_MAIN_VERSION ;
+ version [1] = SUITESPARSE_SUB_VERSION ;
+ version [2] = SUITESPARSE_SUBSUB_VERSION ;
+ }
+ return (SUITESPARSE_VERSION) ;
+}
+
+/* -------------------------------------------------------------------------- */
+/* SuiteSparse_hypot */
+/* -------------------------------------------------------------------------- */
+
+/* There is an equivalent routine called hypot in <math.h>, which conforms
+ * to ANSI C99. However, SuiteSparse does not assume that ANSI C99 is
+ * available. You can use the ANSI C99 hypot routine with:
+ *
+ * #include <math.h>
+ *i SuiteSparse_config.hypot_func = hypot ;
+ *
+ * Default value of the SuiteSparse_config.hypot_func pointer is
+ * SuiteSparse_hypot, defined below.
+ *
+ * s = hypot (x,y) computes s = sqrt (x*x + y*y) but does so more accurately.
+ * The NaN cases for the c_float relops x >= y and x+y == x are safely ignored.
+ *
+ * Source: Algorithm 312, "Absolute value and square root of a complex number,"
+ * P. Friedland, Comm. ACM, vol 10, no 10, October 1967, page 665.
+ */
+
+c_float SuiteSparse_hypot (c_float x, c_float y)
+{
+ c_float s, r ;
+ x = fabs (x) ;
+ y = fabs (y) ;
+ if (x >= y)
+ {
+ if (x + y == x)
+ {
+ s = x ;
+ }
+ else
+ {
+ r = y / x ;
+ s = x * sqrt (1.0 + r*r) ;
+ }
+ }
+ else
+ {
+ if (y + x == y)
+ {
+ s = y ;
+ }
+ else
+ {
+ r = x / y ;
+ s = y * sqrt (1.0 + r*r) ;
+ }
+ }
+ return (s) ;
+}
+
+/* -------------------------------------------------------------------------- */
+/* SuiteSparse_divcomplex */
+/* -------------------------------------------------------------------------- */
+
+/* c = a/b where c, a, and b are complex. The real and imaginary parts are
+ * passed as separate arguments to this routine. The NaN case is ignored
+ * for the c_float relop br >= bi. Returns 1 if the denominator is zero,
+ * 0 otherwise.
+ *
+ * This uses ACM Algo 116, by R. L. Smith, 1962, which tries to avoid
+ * underflow and overflow.
+ *
+ * c can be the same variable as a or b.
+ *
+ * Default value of the SuiteSparse_config.divcomplex_func pointer is
+ * SuiteSparse_divcomplex.
+ */
+
+int SuiteSparse_divcomplex
+(
+ c_float ar, c_float ai, /* real and imaginary parts of a */
+ c_float br, c_float bi, /* real and imaginary parts of b */
+ c_float *cr, c_float *ci /* real and imaginary parts of c */
+)
+{
+ c_float tr, ti, r, den ;
+ if (fabs (br) >= fabs (bi))
+ {
+ r = bi / br ;
+ den = br + r * bi ;
+ tr = (ar + ai * r) / den ;
+ ti = (ai - ar * r) / den ;
+ }
+ else
+ {
+ r = br / bi ;
+ den = r * br + bi ;
+ tr = (ar * r + ai) / den ;
+ ti = (ai * r - ar) / den ;
+ }
+ *cr = tr ;
+ *ci = ti ;
+ return (den == 0.) ;
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_1.c b/lin_sys/direct/qdldl/amd/src/amd_1.c
new file mode 100644
index 0000000..06c8be6
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_1.c
@@ -0,0 +1,180 @@
+/* ========================================================================= */
+/* === AMD_1 =============================================================== */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* AMD_1: Construct A+A' for a sparse matrix A and perform the AMD ordering.
+ *
+ * The n-by-n sparse matrix A can be unsymmetric. It is stored in MATLAB-style
+ * compressed-column form, with sorted row indices in each column, and no
+ * duplicate entries. Diagonal entries may be present, but they are ignored.
+ * Row indices of column j of A are stored in Ai [Ap [j] ... Ap [j+1]-1].
+ * Ap [0] must be zero, and nz = Ap [n] is the number of entries in A. The
+ * size of the matrix, n, must be greater than or equal to zero.
+ *
+ * This routine must be preceded by a call to AMD_aat, which computes the
+ * number of entries in each row/column in A+A', excluding the diagonal.
+ * Len [j], on input, is the number of entries in row/column j of A+A'. This
+ * routine constructs the matrix A+A' and then calls AMD_2. No error checking
+ * is performed (this was done in AMD_valid).
+ */
+
+#include "amd_internal.h"
+
+GLOBAL void AMD_1
+(
+ Int n, /* n > 0 */
+ const Int Ap [ ], /* input of size n+1, not modified */
+ const Int Ai [ ], /* input of size nz = Ap [n], not modified */
+ Int P [ ], /* size n output permutation */
+ Int Pinv [ ], /* size n output inverse permutation */
+ Int Len [ ], /* size n input, undefined on output */
+ Int slen, /* slen >= sum (Len [0..n-1]) + 7n,
+ * ideally slen = 1.2 * sum (Len) + 8n */
+ Int S [ ], /* size slen workspace */
+ c_float Control [ ], /* input array of size AMD_CONTROL */
+ c_float Info [ ] /* output array of size AMD_INFO */
+)
+{
+ Int i, j, k, p, pfree, iwlen, pj, p1, p2, pj2, *Iw, *Pe, *Nv, *Head,
+ *Elen, *Degree, *s, *W, *Sp, *Tp ;
+
+ /* --------------------------------------------------------------------- */
+ /* construct the matrix for AMD_2 */
+ /* --------------------------------------------------------------------- */
+
+ ASSERT (n > 0) ;
+
+ iwlen = slen - 6*n ;
+ s = S ;
+ Pe = s ; s += n ;
+ Nv = s ; s += n ;
+ Head = s ; s += n ;
+ Elen = s ; s += n ;
+ Degree = s ; s += n ;
+ W = s ; s += n ;
+ Iw = s ; s += iwlen ;
+
+ ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
+
+ /* construct the pointers for A+A' */
+ Sp = Nv ; /* use Nv and W as workspace for Sp and Tp [ */
+ Tp = W ;
+ pfree = 0 ;
+ for (j = 0 ; j < n ; j++)
+ {
+ Pe [j] = pfree ;
+ Sp [j] = pfree ;
+ pfree += Len [j] ;
+ }
+
+ /* Note that this restriction on iwlen is slightly more restrictive than
+ * what is strictly required in AMD_2. AMD_2 can operate with no elbow
+ * room at all, but it will be very slow. For better performance, at
+ * least size-n elbow room is enforced. */
+ ASSERT (iwlen >= pfree + n) ;
+
+#ifndef NDEBUG
+ for (p = 0 ; p < iwlen ; p++) Iw [p] = EMPTY ;
+#endif
+
+ for (k = 0 ; k < n ; k++)
+ {
+ AMD_DEBUG1 (("Construct row/column k= "ID" of A+A'\n", k)) ;
+ p1 = Ap [k] ;
+ p2 = Ap [k+1] ;
+
+ /* construct A+A' */
+ for (p = p1 ; p < p2 ; )
+ {
+ /* scan the upper triangular part of A */
+ j = Ai [p] ;
+ ASSERT (j >= 0 && j < n) ;
+ if (j < k)
+ {
+ /* entry A (j,k) in the strictly upper triangular part */
+ ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
+ ASSERT (Sp [k] < (k == n-1 ? pfree : Pe [k+1])) ;
+ Iw [Sp [j]++] = k ;
+ Iw [Sp [k]++] = j ;
+ p++ ;
+ }
+ else if (j == k)
+ {
+ /* skip the diagonal */
+ p++ ;
+ break ;
+ }
+ else /* j > k */
+ {
+ /* first entry below the diagonal */
+ break ;
+ }
+ /* scan lower triangular part of A, in column j until reaching
+ * row k. Start where last scan left off. */
+ ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;
+ pj2 = Ap [j+1] ;
+ for (pj = Tp [j] ; pj < pj2 ; )
+ {
+ i = Ai [pj] ;
+ ASSERT (i >= 0 && i < n) ;
+ if (i < k)
+ {
+ /* A (i,j) is only in the lower part, not in upper */
+ ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ;
+ ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
+ Iw [Sp [i]++] = j ;
+ Iw [Sp [j]++] = i ;
+ pj++ ;
+ }
+ else if (i == k)
+ {
+ /* entry A (k,j) in lower part and A (j,k) in upper */
+ pj++ ;
+ break ;
+ }
+ else /* i > k */
+ {
+ /* consider this entry later, when k advances to i */
+ break ;
+ }
+ }
+ Tp [j] = pj ;
+ }
+ Tp [k] = p ;
+ }
+
+ /* clean up, for remaining mismatched entries */
+ for (j = 0 ; j < n ; j++)
+ {
+ for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)
+ {
+ i = Ai [pj] ;
+ ASSERT (i >= 0 && i < n) ;
+ /* A (i,j) is only in the lower part, not in upper */
+ ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ;
+ ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
+ Iw [Sp [i]++] = j ;
+ Iw [Sp [j]++] = i ;
+ }
+ }
+
+#ifndef NDEBUG
+ for (j = 0 ; j < n-1 ; j++) ASSERT (Sp [j] == Pe [j+1]) ;
+ ASSERT (Sp [n-1] == pfree) ;
+#endif
+
+ /* Tp and Sp no longer needed ] */
+
+ /* --------------------------------------------------------------------- */
+ /* order the matrix */
+ /* --------------------------------------------------------------------- */
+
+ AMD_2 (n, Pe, Iw, Len, iwlen, pfree,
+ Nv, Pinv, P, Head, Elen, Degree, W, Control, Info) ;
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_2.c b/lin_sys/direct/qdldl/amd/src/amd_2.c
new file mode 100644
index 0000000..fe06e75
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_2.c
@@ -0,0 +1,1842 @@
+/* ========================================================================= */
+/* === AMD_2 =============================================================== */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* AMD_2: performs the AMD ordering on a symmetric sparse matrix A, followed
+ * by a postordering (via depth-first search) of the assembly tree using the
+ * AMD_postorder routine.
+ */
+
+#include "amd_internal.h"
+
+/* ========================================================================= */
+/* === clear_flag ========================================================== */
+/* ========================================================================= */
+
+static Int clear_flag (Int wflg, Int wbig, Int W [ ], Int n)
+{
+ Int x ;
+ if (wflg < 2 || wflg >= wbig)
+ {
+ for (x = 0 ; x < n ; x++)
+ {
+ if (W [x] != 0) W [x] = 1 ;
+ }
+ wflg = 2 ;
+ }
+ /* at this point, W [0..n-1] < wflg holds */
+ return (wflg) ;
+}
+
+
+/* ========================================================================= */
+/* === AMD_2 =============================================================== */
+/* ========================================================================= */
+
+GLOBAL void AMD_2
+(
+ Int n, /* A is n-by-n, where n > 0 */
+ Int Pe [ ], /* Pe [0..n-1]: index in Iw of row i on input */
+ Int Iw [ ], /* workspace of size iwlen. Iw [0..pfree-1]
+ * holds the matrix on input */
+ Int Len [ ], /* Len [0..n-1]: length for row/column i on input */
+ Int iwlen, /* length of Iw. iwlen >= pfree + n */
+ Int pfree, /* Iw [pfree ... iwlen-1] is empty on input */
+
+ /* 7 size-n workspaces, not defined on input: */
+ Int Nv [ ], /* the size of each supernode on output */
+ Int Next [ ], /* the output inverse permutation */
+ Int Last [ ], /* the output permutation */
+ Int Head [ ],
+ Int Elen [ ], /* the size columns of L for each supernode */
+ Int Degree [ ],
+ Int W [ ],
+
+ /* control parameters and output statistics */
+ c_float Control [ ], /* array of size AMD_CONTROL */
+ c_float Info [ ] /* array of size AMD_INFO */
+)
+{
+
+/*
+ * Given a representation of the nonzero pattern of a symmetric matrix, A,
+ * (excluding the diagonal) perform an approximate minimum (UMFPACK/MA38-style)
+ * degree ordering to compute a pivot order such that the introduction of
+ * nonzeros (fill-in) in the Cholesky factors A = LL' is kept low. At each
+ * step, the pivot selected is the one with the minimum UMFAPACK/MA38-style
+ * upper-bound on the external degree. This routine can optionally perform
+ * aggresive absorption (as done by MC47B in the Harwell Subroutine
+ * Library).
+ *
+ * The approximate degree algorithm implemented here is the symmetric analog of
+ * the degree update algorithm in MA38 and UMFPACK (the Unsymmetric-pattern
+ * MultiFrontal PACKage, both by Davis and Duff). The routine is based on the
+ * MA27 minimum degree ordering algorithm by Iain Duff and John Reid.
+ *
+ * This routine is a translation of the original AMDBAR and MC47B routines,
+ * in Fortran, with the following modifications:
+ *
+ * (1) dense rows/columns are removed prior to ordering the matrix, and placed
+ * last in the output order. The presence of a dense row/column can
+ * increase the ordering time by up to O(n^2), unless they are removed
+ * prior to ordering.
+ *
+ * (2) the minimum degree ordering is followed by a postordering (depth-first
+ * search) of the assembly tree. Note that mass elimination (discussed
+ * below) combined with the approximate degree update can lead to the mass
+ * elimination of nodes with lower exact degree than the current pivot
+ * element. No additional fill-in is caused in the representation of the
+ * Schur complement. The mass-eliminated nodes merge with the current
+ * pivot element. They are ordered prior to the current pivot element.
+ * Because they can have lower exact degree than the current element, the
+ * merger of two or more of these nodes in the current pivot element can
+ * lead to a single element that is not a "fundamental supernode". The
+ * diagonal block can have zeros in it. Thus, the assembly tree used here
+ * is not guaranteed to be the precise supernodal elemination tree (with
+ * "funadmental" supernodes), and the postordering performed by this
+ * routine is not guaranteed to be a precise postordering of the
+ * elimination tree.
+ *
+ * (3) input parameters are added, to control aggressive absorption and the
+ * detection of "dense" rows/columns of A.
+ *
+ * (4) additional statistical information is returned, such as the number of
+ * nonzeros in L, and the flop counts for subsequent LDL' and LU
+ * factorizations. These are slight upper bounds, because of the mass
+ * elimination issue discussed above.
+ *
+ * (5) additional routines are added to interface this routine to MATLAB
+ * to provide a simple C-callable user-interface, to check inputs for
+ * errors, compute the symmetry of the pattern of A and the number of
+ * nonzeros in each row/column of A+A', to compute the pattern of A+A',
+ * to perform the assembly tree postordering, and to provide debugging
+ * ouput. Many of these functions are also provided by the Fortran
+ * Harwell Subroutine Library routine MC47A.
+ *
+ * (6) both int and SuiteSparse_long versions are provided. In the
+ * descriptions below and integer is and int or SuiteSparse_long depending
+ * on which version is being used.
+
+ **********************************************************************
+ ***** CAUTION: ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT. ******
+ **********************************************************************
+ ** If you want error checking, a more versatile input format, and a **
+ ** simpler user interface, use amd_order or amd_l_order instead. **
+ ** This routine is not meant to be user-callable. **
+ **********************************************************************
+
+ * ----------------------------------------------------------------------------
+ * References:
+ * ----------------------------------------------------------------------------
+ *
+ * [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern multifrontal
+ * method for sparse LU factorization", SIAM J. Matrix Analysis and
+ * Applications, vol. 18, no. 1, pp. 140-158. Discusses UMFPACK / MA38,
+ * which first introduced the approximate minimum degree used by this
+ * routine.
+ *
+ * [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, "An approximate
+ * minimum degree ordering algorithm," SIAM J. Matrix Analysis and
+ * Applications, vol. 17, no. 4, pp. 886-905, 1996. Discusses AMDBAR and
+ * MC47B, which are the Fortran versions of this routine.
+ *
+ * [3] Alan George and Joseph Liu, "The evolution of the minimum degree
+ * ordering algorithm," SIAM Review, vol. 31, no. 1, pp. 1-19, 1989.
+ * We list below the features mentioned in that paper that this code
+ * includes:
+ *
+ * mass elimination:
+ * Yes. MA27 relied on supervariable detection for mass elimination.
+ *
+ * indistinguishable nodes:
+ * Yes (we call these "supervariables"). This was also in the MA27
+ * code - although we modified the method of detecting them (the
+ * previous hash was the true degree, which we no longer keep track
+ * of). A supervariable is a set of rows with identical nonzero
+ * pattern. All variables in a supervariable are eliminated together.
+ * Each supervariable has as its numerical name that of one of its
+ * variables (its principal variable).
+ *
+ * quotient graph representation:
+ * Yes. We use the term "element" for the cliques formed during
+ * elimination. This was also in the MA27 code. The algorithm can
+ * operate in place, but it will work more efficiently if given some
+ * "elbow room."
+ *
+ * element absorption:
+ * Yes. This was also in the MA27 code.
+ *
+ * external degree:
+ * Yes. The MA27 code was based on the true degree.
+ *
+ * incomplete degree update and multiple elimination:
+ * No. This was not in MA27, either. Our method of degree update
+ * within MC47B is element-based, not variable-based. It is thus
+ * not well-suited for use with incomplete degree update or multiple
+ * elimination.
+ *
+ * Authors, and Copyright (C) 2004 by:
+ * Timothy A. Davis, Patrick Amestoy, Iain S. Duff, John K. Reid.
+ *
+ * Acknowledgements: This work (and the UMFPACK package) was supported by the
+ * National Science Foundation (ASC-9111263, DMS-9223088, and CCR-0203270).
+ * The UMFPACK/MA38 approximate degree update algorithm, the unsymmetric analog
+ * which forms the basis of AMD, was developed while Tim Davis was supported by
+ * CERFACS (Toulouse, France) in a post-doctoral position. This C version, and
+ * the etree postorder, were written while Tim Davis was on sabbatical at
+ * Stanford University and Lawrence Berkeley National Laboratory.
+
+ * ----------------------------------------------------------------------------
+ * INPUT ARGUMENTS (unaltered):
+ * ----------------------------------------------------------------------------
+
+ * n: The matrix order. Restriction: n >= 1.
+ *
+ * iwlen: The size of the Iw array. On input, the matrix is stored in
+ * Iw [0..pfree-1]. However, Iw [0..iwlen-1] should be slightly larger
+ * than what is required to hold the matrix, at least iwlen >= pfree + n.
+ * Otherwise, excessive compressions will take place. The recommended
+ * value of iwlen is 1.2 * pfree + n, which is the value used in the
+ * user-callable interface to this routine (amd_order.c). The algorithm
+ * will not run at all if iwlen < pfree. Restriction: iwlen >= pfree + n.
+ * Note that this is slightly more restrictive than the actual minimum
+ * (iwlen >= pfree), but AMD_2 will be very slow with no elbow room.
+ * Thus, this routine enforces a bare minimum elbow room of size n.
+ *
+ * pfree: On input the tail end of the array, Iw [pfree..iwlen-1], is empty,
+ * and the matrix is stored in Iw [0..pfree-1]. During execution,
+ * additional data is placed in Iw, and pfree is modified so that
+ * Iw [pfree..iwlen-1] is always the unused part of Iw.
+ *
+ * Control: A c_float array of size AMD_CONTROL containing input parameters
+ * that affect how the ordering is computed. If NULL, then default
+ * settings are used.
+ *
+ * Control [AMD_DENSE] is used to determine whether or not a given input
+ * row is "dense". A row is "dense" if the number of entries in the row
+ * exceeds Control [AMD_DENSE] times sqrt (n), except that rows with 16 or
+ * fewer entries are never considered "dense". To turn off the detection
+ * of dense rows, set Control [AMD_DENSE] to a negative number, or to a
+ * number larger than sqrt (n). The default value of Control [AMD_DENSE]
+ * is AMD_DEFAULT_DENSE, which is defined in amd.h as 10.
+ *
+ * Control [AMD_AGGRESSIVE] is used to determine whether or not aggressive
+ * absorption is to be performed. If nonzero, then aggressive absorption
+ * is performed (this is the default).
+
+ * ----------------------------------------------------------------------------
+ * INPUT/OUPUT ARGUMENTS:
+ * ----------------------------------------------------------------------------
+ *
+ * Pe: An integer array of size n. On input, Pe [i] is the index in Iw of
+ * the start of row i. Pe [i] is ignored if row i has no off-diagonal
+ * entries. Thus Pe [i] must be in the range 0 to pfree-1 for non-empty
+ * rows.
+ *
+ * During execution, it is used for both supervariables and elements:
+ *
+ * Principal supervariable i: index into Iw of the description of
+ * supervariable i. A supervariable represents one or more rows of
+ * the matrix with identical nonzero pattern. In this case,
+ * Pe [i] >= 0.
+ *
+ * Non-principal supervariable i: if i has been absorbed into another
+ * supervariable j, then Pe [i] = FLIP (j), where FLIP (j) is defined
+ * as (-(j)-2). Row j has the same pattern as row i. Note that j
+ * might later be absorbed into another supervariable j2, in which
+ * case Pe [i] is still FLIP (j), and Pe [j] = FLIP (j2) which is
+ * < EMPTY, where EMPTY is defined as (-1) in amd_internal.h.
+ *
+ * Unabsorbed element e: the index into Iw of the description of element
+ * e, if e has not yet been absorbed by a subsequent element. Element
+ * e is created when the supervariable of the same name is selected as
+ * the pivot. In this case, Pe [i] >= 0.
+ *
+ * Absorbed element e: if element e is absorbed into element e2, then
+ * Pe [e] = FLIP (e2). This occurs when the pattern of e (which we
+ * refer to as Le) is found to be a subset of the pattern of e2 (that
+ * is, Le2). In this case, Pe [i] < EMPTY. If element e is "null"
+ * (it has no nonzeros outside its pivot block), then Pe [e] = EMPTY,
+ * and e is the root of an assembly subtree (or the whole tree if
+ * there is just one such root).
+ *
+ * Dense variable i: if i is "dense", then Pe [i] = EMPTY.
+ *
+ * On output, Pe holds the assembly tree/forest, which implicitly
+ * represents a pivot order with identical fill-in as the actual order
+ * (via a depth-first search of the tree), as follows. If Nv [i] > 0,
+ * then i represents a node in the assembly tree, and the parent of i is
+ * Pe [i], or EMPTY if i is a root. If Nv [i] = 0, then (i, Pe [i])
+ * represents an edge in a subtree, the root of which is a node in the
+ * assembly tree. Note that i refers to a row/column in the original
+ * matrix, not the permuted matrix.
+ *
+ * Info: A c_float array of size AMD_INFO. If present, (that is, not NULL),
+ * then statistics about the ordering are returned in the Info array.
+ * See amd.h for a description.
+
+ * ----------------------------------------------------------------------------
+ * INPUT/MODIFIED (undefined on output):
+ * ----------------------------------------------------------------------------
+ *
+ * Len: An integer array of size n. On input, Len [i] holds the number of
+ * entries in row i of the matrix, excluding the diagonal. The contents
+ * of Len are undefined on output.
+ *
+ * Iw: An integer array of size iwlen. On input, Iw [0..pfree-1] holds the
+ * description of each row i in the matrix. The matrix must be symmetric,
+ * and both upper and lower triangular parts must be present. The
+ * diagonal must not be present. Row i is held as follows:
+ *
+ * Len [i]: the length of the row i data structure in the Iw array.
+ * Iw [Pe [i] ... Pe [i] + Len [i] - 1]:
+ * the list of column indices for nonzeros in row i (simple
+ * supervariables), excluding the diagonal. All supervariables
+ * start with one row/column each (supervariable i is just row i).
+ * If Len [i] is zero on input, then Pe [i] is ignored on input.
+ *
+ * Note that the rows need not be in any particular order, and there
+ * may be empty space between the rows.
+ *
+ * During execution, the supervariable i experiences fill-in. This is
+ * represented by placing in i a list of the elements that cause fill-in
+ * in supervariable i:
+ *
+ * Len [i]: the length of supervariable i in the Iw array.
+ * Iw [Pe [i] ... Pe [i] + Elen [i] - 1]:
+ * the list of elements that contain i. This list is kept short
+ * by removing absorbed elements.
+ * Iw [Pe [i] + Elen [i] ... Pe [i] + Len [i] - 1]:
+ * the list of supervariables in i. This list is kept short by
+ * removing nonprincipal variables, and any entry j that is also
+ * contained in at least one of the elements (j in Le) in the list
+ * for i (e in row i).
+ *
+ * When supervariable i is selected as pivot, we create an element e of
+ * the same name (e=i):
+ *
+ * Len [e]: the length of element e in the Iw array.
+ * Iw [Pe [e] ... Pe [e] + Len [e] - 1]:
+ * the list of supervariables in element e.
+ *
+ * An element represents the fill-in that occurs when supervariable i is
+ * selected as pivot (which represents the selection of row i and all
+ * non-principal variables whose principal variable is i). We use the
+ * term Le to denote the set of all supervariables in element e. Absorbed
+ * supervariables and elements are pruned from these lists when
+ * computationally convenient.
+ *
+ * CAUTION: THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION.
+ * The contents of Iw are undefined on output.
+
+ * ----------------------------------------------------------------------------
+ * OUTPUT (need not be set on input):
+ * ----------------------------------------------------------------------------
+ *
+ * Nv: An integer array of size n. During execution, ABS (Nv [i]) is equal to
+ * the number of rows that are represented by the principal supervariable
+ * i. If i is a nonprincipal or dense variable, then Nv [i] = 0.
+ * Initially, Nv [i] = 1 for all i. Nv [i] < 0 signifies that i is a
+ * principal variable in the pattern Lme of the current pivot element me.
+ * After element me is constructed, Nv [i] is set back to a positive
+ * value.
+ *
+ * On output, Nv [i] holds the number of pivots represented by super
+ * row/column i of the original matrix, or Nv [i] = 0 for non-principal
+ * rows/columns. Note that i refers to a row/column in the original
+ * matrix, not the permuted matrix.
+ *
+ * Elen: An integer array of size n. See the description of Iw above. At the
+ * start of execution, Elen [i] is set to zero for all rows i. During
+ * execution, Elen [i] is the number of elements in the list for
+ * supervariable i. When e becomes an element, Elen [e] = FLIP (esize) is
+ * set, where esize is the size of the element (the number of pivots, plus
+ * the number of nonpivotal entries). Thus Elen [e] < EMPTY.
+ * Elen (i) = EMPTY set when variable i becomes nonprincipal.
+ *
+ * For variables, Elen (i) >= EMPTY holds until just before the
+ * postordering and permutation vectors are computed. For elements,
+ * Elen [e] < EMPTY holds.
+ *
+ * On output, Elen [i] is the degree of the row/column in the Cholesky
+ * factorization of the permuted matrix, corresponding to the original row
+ * i, if i is a super row/column. It is equal to EMPTY if i is
+ * non-principal. Note that i refers to a row/column in the original
+ * matrix, not the permuted matrix.
+ *
+ * Note that the contents of Elen on output differ from the Fortran
+ * version (Elen holds the inverse permutation in the Fortran version,
+ * which is instead returned in the Next array in this C version,
+ * described below).
+ *
+ * Last: In a degree list, Last [i] is the supervariable preceding i, or EMPTY
+ * if i is the head of the list. In a hash bucket, Last [i] is the hash
+ * key for i.
+ *
+ * Last [Head [hash]] is also used as the head of a hash bucket if
+ * Head [hash] contains a degree list (see the description of Head,
+ * below).
+ *
+ * On output, Last [0..n-1] holds the permutation. That is, if
+ * i = Last [k], then row i is the kth pivot row (where k ranges from 0 to
+ * n-1). Row Last [k] of A is the kth row in the permuted matrix, PAP'.
+ *
+ * Next: Next [i] is the supervariable following i in a link list, or EMPTY if
+ * i is the last in the list. Used for two kinds of lists: degree lists
+ * and hash buckets (a supervariable can be in only one kind of list at a
+ * time).
+ *
+ * On output Next [0..n-1] holds the inverse permutation. That is, if
+ * k = Next [i], then row i is the kth pivot row. Row i of A appears as
+ * the (Next[i])-th row in the permuted matrix, PAP'.
+ *
+ * Note that the contents of Next on output differ from the Fortran
+ * version (Next is undefined on output in the Fortran version).
+
+ * ----------------------------------------------------------------------------
+ * LOCAL WORKSPACE (not input or output - used only during execution):
+ * ----------------------------------------------------------------------------
+ *
+ * Degree: An integer array of size n. If i is a supervariable, then
+ * Degree [i] holds the current approximation of the external degree of
+ * row i (an upper bound). The external degree is the number of nonzeros
+ * in row i, minus ABS (Nv [i]), the diagonal part. The bound is equal to
+ * the exact external degree if Elen [i] is less than or equal to two.
+ *
+ * We also use the term "external degree" for elements e to refer to
+ * |Le \ Lme|. If e is an element, then Degree [e] is |Le|, which is the
+ * degree of the off-diagonal part of the element e (not including the
+ * diagonal part).
+ *
+ * Head: An integer array of size n. Head is used for degree lists.
+ * Head [deg] is the first supervariable in a degree list. All
+ * supervariables i in a degree list Head [deg] have the same approximate
+ * degree, namely, deg = Degree [i]. If the list Head [deg] is empty then
+ * Head [deg] = EMPTY.
+ *
+ * During supervariable detection Head [hash] also serves as a pointer to
+ * a hash bucket. If Head [hash] >= 0, there is a degree list of degree
+ * hash. The hash bucket head pointer is Last [Head [hash]]. If
+ * Head [hash] = EMPTY, then the degree list and hash bucket are both
+ * empty. If Head [hash] < EMPTY, then the degree list is empty, and
+ * FLIP (Head [hash]) is the head of the hash bucket. After supervariable
+ * detection is complete, all hash buckets are empty, and the
+ * (Last [Head [hash]] = EMPTY) condition is restored for the non-empty
+ * degree lists.
+ *
+ * W: An integer array of size n. The flag array W determines the status of
+ * elements and variables, and the external degree of elements.
+ *
+ * for elements:
+ * if W [e] = 0, then the element e is absorbed.
+ * if W [e] >= wflg, then W [e] - wflg is the size of the set
+ * |Le \ Lme|, in terms of nonzeros (the sum of ABS (Nv [i]) for
+ * each principal variable i that is both in the pattern of
+ * element e and NOT in the pattern of the current pivot element,
+ * me).
+ * if wflg > W [e] > 0, then e is not absorbed and has not yet been
+ * seen in the scan of the element lists in the computation of
+ * |Le\Lme| in Scan 1 below.
+ *
+ * for variables:
+ * during supervariable detection, if W [j] != wflg then j is
+ * not in the pattern of variable i.
+ *
+ * The W array is initialized by setting W [i] = 1 for all i, and by
+ * setting wflg = 2. It is reinitialized if wflg becomes too large (to
+ * ensure that wflg+n does not cause integer overflow).
+
+ * ----------------------------------------------------------------------------
+ * LOCAL INTEGERS:
+ * ----------------------------------------------------------------------------
+ */
+
+ Int deg, degme, dext, lemax, e, elenme, eln, i, ilast, inext, j,
+ jlast, jnext, k, knt1, knt2, knt3, lenj, ln, me, mindeg, nel, nleft,
+ nvi, nvj, nvpiv, slenme, wbig, we, wflg, wnvi, ok, ndense, ncmpa,
+ dense, aggressive ;
+
+ unsigned Int hash ; /* unsigned, so that hash % n is well defined.*/
+
+/*
+ * deg: the degree of a variable or element
+ * degme: size, |Lme|, of the current element, me (= Degree [me])
+ * dext: external degree, |Le \ Lme|, of some element e
+ * lemax: largest |Le| seen so far (called dmax in Fortran version)
+ * e: an element
+ * elenme: the length, Elen [me], of element list of pivotal variable
+ * eln: the length, Elen [...], of an element list
+ * hash: the computed value of the hash function
+ * i: a supervariable
+ * ilast: the entry in a link list preceding i
+ * inext: the entry in a link list following i
+ * j: a supervariable
+ * jlast: the entry in a link list preceding j
+ * jnext: the entry in a link list, or path, following j
+ * k: the pivot order of an element or variable
+ * knt1: loop counter used during element construction
+ * knt2: loop counter used during element construction
+ * knt3: loop counter used during compression
+ * lenj: Len [j]
+ * ln: length of a supervariable list
+ * me: current supervariable being eliminated, and the current
+ * element created by eliminating that supervariable
+ * mindeg: current minimum degree
+ * nel: number of pivots selected so far
+ * nleft: n - nel, the number of nonpivotal rows/columns remaining
+ * nvi: the number of variables in a supervariable i (= Nv [i])
+ * nvj: the number of variables in a supervariable j (= Nv [j])
+ * nvpiv: number of pivots in current element
+ * slenme: number of variables in variable list of pivotal variable
+ * wbig: = (INT_MAX - n) for the int version, (SuiteSparse_long_max - n)
+ * for the SuiteSparse_long version. wflg is not allowed to
+ * be >= wbig.
+ * we: W [e]
+ * wflg: used for flagging the W array. See description of Iw.
+ * wnvi: wflg - Nv [i]
+ * x: either a supervariable or an element
+ *
+ * ok: true if supervariable j can be absorbed into i
+ * ndense: number of "dense" rows/columns
+ * dense: rows/columns with initial degree > dense are considered "dense"
+ * aggressive: true if aggressive absorption is being performed
+ * ncmpa: number of garbage collections
+
+ * ----------------------------------------------------------------------------
+ * LOCAL DOUBLES, used for statistical output only (except for alpha):
+ * ----------------------------------------------------------------------------
+ */
+
+ c_float f, r, ndiv, s, nms_lu, nms_ldl, dmax, alpha, lnz, lnzme ;
+
+/*
+ * f: nvpiv
+ * r: degme + nvpiv
+ * ndiv: number of divisions for LU or LDL' factorizations
+ * s: number of multiply-subtract pairs for LU factorization, for the
+ * current element me
+ * nms_lu number of multiply-subtract pairs for LU factorization
+ * nms_ldl number of multiply-subtract pairs for LDL' factorization
+ * dmax: the largest number of entries in any column of L, including the
+ * diagonal
+ * alpha: "dense" degree ratio
+ * lnz: the number of nonzeros in L (excluding the diagonal)
+ * lnzme: the number of nonzeros in L (excl. the diagonal) for the
+ * current element me
+
+ * ----------------------------------------------------------------------------
+ * LOCAL "POINTERS" (indices into the Iw array)
+ * ----------------------------------------------------------------------------
+*/
+
+ Int p, p1, p2, p3, p4, pdst, pend, pj, pme, pme1, pme2, pn, psrc ;
+
+/*
+ * Any parameter (Pe [...] or pfree) or local variable starting with "p" (for
+ * Pointer) is an index into Iw, and all indices into Iw use variables starting
+ * with "p." The only exception to this rule is the iwlen input argument.
+ *
+ * p: pointer into lots of things
+ * p1: Pe [i] for some variable i (start of element list)
+ * p2: Pe [i] + Elen [i] - 1 for some variable i
+ * p3: index of first supervariable in clean list
+ * p4:
+ * pdst: destination pointer, for compression
+ * pend: end of memory to compress
+ * pj: pointer into an element or variable
+ * pme: pointer into the current element (pme1...pme2)
+ * pme1: the current element, me, is stored in Iw [pme1...pme2]
+ * pme2: the end of the current element
+ * pn: pointer into a "clean" variable, also used to compress
+ * psrc: source pointer, for compression
+*/
+
+/* ========================================================================= */
+/* INITIALIZATIONS */
+/* ========================================================================= */
+
+ /* Note that this restriction on iwlen is slightly more restrictive than
+ * what is actually required in AMD_2. AMD_2 can operate with no elbow
+ * room at all, but it will be slow. For better performance, at least
+ * size-n elbow room is enforced. */
+ ASSERT (iwlen >= pfree + n) ;
+ ASSERT (n > 0) ;
+
+ /* initialize output statistics */
+ lnz = 0 ;
+ ndiv = 0 ;
+ nms_lu = 0 ;
+ nms_ldl = 0 ;
+ dmax = 1 ;
+ me = EMPTY ;
+
+ mindeg = 0 ;
+ ncmpa = 0 ;
+ nel = 0 ;
+ lemax = 0 ;
+
+ /* get control parameters */
+ if (Control != (c_float *) NULL)
+ {
+ alpha = Control [AMD_DENSE] ;
+ aggressive = (Control [AMD_AGGRESSIVE] != 0) ;
+ }
+ else
+ {
+ alpha = AMD_DEFAULT_DENSE ;
+ aggressive = AMD_DEFAULT_AGGRESSIVE ;
+ }
+ /* Note: if alpha is NaN, this is undefined: */
+ if (alpha < 0)
+ {
+ /* only remove completely dense rows/columns */
+ dense = n-2 ;
+ }
+ else
+ {
+ dense = (Int) (alpha * sqrt ((c_float) n)) ;
+ }
+ dense = MAX (16, dense) ;
+ dense = MIN (n, dense) ;
+ AMD_DEBUG1 (("\n\nAMD (debug), alpha %g, aggr. "ID"\n",
+ alpha, aggressive)) ;
+
+ for (i = 0 ; i < n ; i++)
+ {
+ Last [i] = EMPTY ;
+ Head [i] = EMPTY ;
+ Next [i] = EMPTY ;
+ /* if separate Hhead array is used for hash buckets: *
+ Hhead [i] = EMPTY ;
+ */
+ Nv [i] = 1 ;
+ W [i] = 1 ;
+ Elen [i] = 0 ;
+ Degree [i] = Len [i] ;
+ }
+
+#ifndef NDEBUG
+ AMD_DEBUG1 (("\n======Nel "ID" initial\n", nel)) ;
+ AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next, Last,
+ Head, Elen, Degree, W, -1) ;
+#endif
+
+ /* initialize wflg */
+ wbig = Int_MAX - n ;
+ wflg = clear_flag (0, wbig, W, n) ;
+
+ /* --------------------------------------------------------------------- */
+ /* initialize degree lists and eliminate dense and empty rows */
+ /* --------------------------------------------------------------------- */
+
+ ndense = 0 ;
+
+ for (i = 0 ; i < n ; i++)
+ {
+ deg = Degree [i] ;
+ ASSERT (deg >= 0 && deg < n) ;
+ if (deg == 0)
+ {
+
+ /* -------------------------------------------------------------
+ * we have a variable that can be eliminated at once because
+ * there is no off-diagonal non-zero in its row. Note that
+ * Nv [i] = 1 for an empty variable i. It is treated just
+ * the same as an eliminated element i.
+ * ------------------------------------------------------------- */
+
+ Elen [i] = FLIP (1) ;
+ nel++ ;
+ Pe [i] = EMPTY ;
+ W [i] = 0 ;
+
+ }
+ else if (deg > dense)
+ {
+
+ /* -------------------------------------------------------------
+ * Dense variables are not treated as elements, but as unordered,
+ * non-principal variables that have no parent. They do not take
+ * part in the postorder, since Nv [i] = 0. Note that the Fortran
+ * version does not have this option.
+ * ------------------------------------------------------------- */
+
+ AMD_DEBUG1 (("Dense node "ID" degree "ID"\n", i, deg)) ;
+ ndense++ ;
+ Nv [i] = 0 ; /* do not postorder this node */
+ Elen [i] = EMPTY ;
+ nel++ ;
+ Pe [i] = EMPTY ;
+
+ }
+ else
+ {
+
+ /* -------------------------------------------------------------
+ * place i in the degree list corresponding to its degree
+ * ------------------------------------------------------------- */
+
+ inext = Head [deg] ;
+ ASSERT (inext >= EMPTY && inext < n) ;
+ if (inext != EMPTY) Last [inext] = i ;
+ Next [i] = inext ;
+ Head [deg] = i ;
+
+ }
+ }
+
+/* ========================================================================= */
+/* WHILE (selecting pivots) DO */
+/* ========================================================================= */
+
+ while (nel < n)
+ {
+
+#ifndef NDEBUG
+ AMD_DEBUG1 (("\n======Nel "ID"\n", nel)) ;
+ if (AMD_debug >= 2)
+ {
+ AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next,
+ Last, Head, Elen, Degree, W, nel) ;
+ }
+#endif
+
+/* ========================================================================= */
+/* GET PIVOT OF MINIMUM DEGREE */
+/* ========================================================================= */
+
+ /* ----------------------------------------------------------------- */
+ /* find next supervariable for elimination */
+ /* ----------------------------------------------------------------- */
+
+ ASSERT (mindeg >= 0 && mindeg < n) ;
+ for (deg = mindeg ; deg < n ; deg++)
+ {
+ me = Head [deg] ;
+ if (me != EMPTY) break ;
+ }
+ mindeg = deg ;
+ ASSERT (me >= 0 && me < n) ;
+ AMD_DEBUG1 (("=================me: "ID"\n", me)) ;
+
+ /* ----------------------------------------------------------------- */
+ /* remove chosen variable from link list */
+ /* ----------------------------------------------------------------- */
+
+ inext = Next [me] ;
+ ASSERT (inext >= EMPTY && inext < n) ;
+ if (inext != EMPTY) Last [inext] = EMPTY ;
+ Head [deg] = inext ;
+
+ /* ----------------------------------------------------------------- */
+ /* me represents the elimination of pivots nel to nel+Nv[me]-1. */
+ /* place me itself as the first in this set. */
+ /* ----------------------------------------------------------------- */
+
+ elenme = Elen [me] ;
+ nvpiv = Nv [me] ;
+ ASSERT (nvpiv > 0) ;
+ nel += nvpiv ;
+
+/* ========================================================================= */
+/* CONSTRUCT NEW ELEMENT */
+/* ========================================================================= */
+
+ /* -----------------------------------------------------------------
+ * At this point, me is the pivotal supervariable. It will be
+ * converted into the current element. Scan list of the pivotal
+ * supervariable, me, setting tree pointers and constructing new list
+ * of supervariables for the new element, me. p is a pointer to the
+ * current position in the old list.
+ * ----------------------------------------------------------------- */
+
+ /* flag the variable "me" as being in Lme by negating Nv [me] */
+ Nv [me] = -nvpiv ;
+ degme = 0 ;
+ ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ;
+
+ if (elenme == 0)
+ {
+
+ /* ------------------------------------------------------------- */
+ /* construct the new element in place */
+ /* ------------------------------------------------------------- */
+
+ pme1 = Pe [me] ;
+ pme2 = pme1 - 1 ;
+
+ for (p = pme1 ; p <= pme1 + Len [me] - 1 ; p++)
+ {
+ i = Iw [p] ;
+ ASSERT (i >= 0 && i < n && Nv [i] >= 0) ;
+ nvi = Nv [i] ;
+ if (nvi > 0)
+ {
+
+ /* ----------------------------------------------------- */
+ /* i is a principal variable not yet placed in Lme. */
+ /* store i in new list */
+ /* ----------------------------------------------------- */
+
+ /* flag i as being in Lme by negating Nv [i] */
+ degme += nvi ;
+ Nv [i] = -nvi ;
+ Iw [++pme2] = i ;
+
+ /* ----------------------------------------------------- */
+ /* remove variable i from degree list. */
+ /* ----------------------------------------------------- */
+
+ ilast = Last [i] ;
+ inext = Next [i] ;
+ ASSERT (ilast >= EMPTY && ilast < n) ;
+ ASSERT (inext >= EMPTY && inext < n) ;
+ if (inext != EMPTY) Last [inext] = ilast ;
+ if (ilast != EMPTY)
+ {
+ Next [ilast] = inext ;
+ }
+ else
+ {
+ /* i is at the head of the degree list */
+ ASSERT (Degree [i] >= 0 && Degree [i] < n) ;
+ Head [Degree [i]] = inext ;
+ }
+ }
+ }
+ }
+ else
+ {
+
+ /* ------------------------------------------------------------- */
+ /* construct the new element in empty space, Iw [pfree ...] */
+ /* ------------------------------------------------------------- */
+
+ p = Pe [me] ;
+ pme1 = pfree ;
+ slenme = Len [me] - elenme ;
+
+ for (knt1 = 1 ; knt1 <= elenme + 1 ; knt1++)
+ {
+
+ if (knt1 > elenme)
+ {
+ /* search the supervariables in me. */
+ e = me ;
+ pj = p ;
+ ln = slenme ;
+ AMD_DEBUG2 (("Search sv: "ID" "ID" "ID"\n", me,pj,ln)) ;
+ }
+ else
+ {
+ /* search the elements in me. */
+ e = Iw [p++] ;
+ ASSERT (e >= 0 && e < n) ;
+ pj = Pe [e] ;
+ ln = Len [e] ;
+ AMD_DEBUG2 (("Search element e "ID" in me "ID"\n", e,me)) ;
+ ASSERT (Elen [e] < EMPTY && W [e] > 0 && pj >= 0) ;
+ }
+ ASSERT (ln >= 0 && (ln == 0 || (pj >= 0 && pj < iwlen))) ;
+
+ /* ---------------------------------------------------------
+ * search for different supervariables and add them to the
+ * new list, compressing when necessary. this loop is
+ * executed once for each element in the list and once for
+ * all the supervariables in the list.
+ * --------------------------------------------------------- */
+
+ for (knt2 = 1 ; knt2 <= ln ; knt2++)
+ {
+ i = Iw [pj++] ;
+ ASSERT (i >= 0 && i < n && (i == me || Elen [i] >= EMPTY));
+ nvi = Nv [i] ;
+ AMD_DEBUG2 ((": "ID" "ID" "ID" "ID"\n",
+ i, Elen [i], Nv [i], wflg)) ;
+
+ if (nvi > 0)
+ {
+
+ /* ------------------------------------------------- */
+ /* compress Iw, if necessary */
+ /* ------------------------------------------------- */
+
+ if (pfree >= iwlen)
+ {
+
+ AMD_DEBUG1 (("GARBAGE COLLECTION\n")) ;
+
+ /* prepare for compressing Iw by adjusting pointers
+ * and lengths so that the lists being searched in
+ * the inner and outer loops contain only the
+ * remaining entries. */
+
+ Pe [me] = p ;
+ Len [me] -= knt1 ;
+ /* check if nothing left of supervariable me */
+ if (Len [me] == 0) Pe [me] = EMPTY ;
+ Pe [e] = pj ;
+ Len [e] = ln - knt2 ;
+ /* nothing left of element e */
+ if (Len [e] == 0) Pe [e] = EMPTY ;
+
+ ncmpa++ ; /* one more garbage collection */
+
+ /* store first entry of each object in Pe */
+ /* FLIP the first entry in each object */
+ for (j = 0 ; j < n ; j++)
+ {
+ pn = Pe [j] ;
+ if (pn >= 0)
+ {
+ ASSERT (pn >= 0 && pn < iwlen) ;
+ Pe [j] = Iw [pn] ;
+ Iw [pn] = FLIP (j) ;
+ }
+ }
+
+ /* psrc/pdst point to source/destination */
+ psrc = 0 ;
+ pdst = 0 ;
+ pend = pme1 - 1 ;
+
+ while (psrc <= pend)
+ {
+ /* search for next FLIP'd entry */
+ j = FLIP (Iw [psrc++]) ;
+ if (j >= 0)
+ {
+ AMD_DEBUG2 (("Got object j: "ID"\n", j)) ;
+ Iw [pdst] = Pe [j] ;
+ Pe [j] = pdst++ ;
+ lenj = Len [j] ;
+ /* copy from source to destination */
+ for (knt3 = 0 ; knt3 <= lenj - 2 ; knt3++)
+ {
+ Iw [pdst++] = Iw [psrc++] ;
+ }
+ }
+ }
+
+ /* move the new partially-constructed element */
+ p1 = pdst ;
+ for (psrc = pme1 ; psrc <= pfree-1 ; psrc++)
+ {
+ Iw [pdst++] = Iw [psrc] ;
+ }
+ pme1 = p1 ;
+ pfree = pdst ;
+ pj = Pe [e] ;
+ p = Pe [me] ;
+
+ }
+
+ /* ------------------------------------------------- */
+ /* i is a principal variable not yet placed in Lme */
+ /* store i in new list */
+ /* ------------------------------------------------- */
+
+ /* flag i as being in Lme by negating Nv [i] */
+ degme += nvi ;
+ Nv [i] = -nvi ;
+ Iw [pfree++] = i ;
+ AMD_DEBUG2 ((" s: "ID" nv "ID"\n", i, Nv [i]));
+
+ /* ------------------------------------------------- */
+ /* remove variable i from degree link list */
+ /* ------------------------------------------------- */
+
+ ilast = Last [i] ;
+ inext = Next [i] ;
+ ASSERT (ilast >= EMPTY && ilast < n) ;
+ ASSERT (inext >= EMPTY && inext < n) ;
+ if (inext != EMPTY) Last [inext] = ilast ;
+ if (ilast != EMPTY)
+ {
+ Next [ilast] = inext ;
+ }
+ else
+ {
+ /* i is at the head of the degree list */
+ ASSERT (Degree [i] >= 0 && Degree [i] < n) ;
+ Head [Degree [i]] = inext ;
+ }
+ }
+ }
+
+ if (e != me)
+ {
+ /* set tree pointer and flag to indicate element e is
+ * absorbed into new element me (the parent of e is me) */
+ AMD_DEBUG1 ((" Element "ID" => "ID"\n", e, me)) ;
+ Pe [e] = FLIP (me) ;
+ W [e] = 0 ;
+ }
+ }
+
+ pme2 = pfree - 1 ;
+ }
+
+ /* ----------------------------------------------------------------- */
+ /* me has now been converted into an element in Iw [pme1..pme2] */
+ /* ----------------------------------------------------------------- */
+
+ /* degme holds the external degree of new element */
+ Degree [me] = degme ;
+ Pe [me] = pme1 ;
+ Len [me] = pme2 - pme1 + 1 ;
+ ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ;
+
+ Elen [me] = FLIP (nvpiv + degme) ;
+ /* FLIP (Elen (me)) is now the degree of pivot (including
+ * diagonal part). */
+
+#ifndef NDEBUG
+ AMD_DEBUG2 (("New element structure: length= "ID"\n", pme2-pme1+1)) ;
+ for (pme = pme1 ; pme <= pme2 ; pme++) AMD_DEBUG3 ((" "ID"", Iw[pme]));
+ AMD_DEBUG3 (("\n")) ;
+#endif
+
+ /* ----------------------------------------------------------------- */
+ /* make sure that wflg is not too large. */
+ /* ----------------------------------------------------------------- */
+
+ /* With the current value of wflg, wflg+n must not cause integer
+ * overflow */
+
+ wflg = clear_flag (wflg, wbig, W, n) ;
+
+/* ========================================================================= */
+/* COMPUTE (W [e] - wflg) = |Le\Lme| FOR ALL ELEMENTS */
+/* ========================================================================= */
+
+ /* -----------------------------------------------------------------
+ * Scan 1: compute the external degrees of previous elements with
+ * respect to the current element. That is:
+ * (W [e] - wflg) = |Le \ Lme|
+ * for each element e that appears in any supervariable in Lme. The
+ * notation Le refers to the pattern (list of supervariables) of a
+ * previous element e, where e is not yet absorbed, stored in
+ * Iw [Pe [e] + 1 ... Pe [e] + Len [e]]. The notation Lme
+ * refers to the pattern of the current element (stored in
+ * Iw [pme1..pme2]). If aggressive absorption is enabled, and
+ * (W [e] - wflg) becomes zero, then the element e will be absorbed
+ * in Scan 2.
+ * ----------------------------------------------------------------- */
+
+ AMD_DEBUG2 (("me: ")) ;
+ for (pme = pme1 ; pme <= pme2 ; pme++)
+ {
+ i = Iw [pme] ;
+ ASSERT (i >= 0 && i < n) ;
+ eln = Elen [i] ;
+ AMD_DEBUG3 ((""ID" Elen "ID": \n", i, eln)) ;
+ if (eln > 0)
+ {
+ /* note that Nv [i] has been negated to denote i in Lme: */
+ nvi = -Nv [i] ;
+ ASSERT (nvi > 0 && Pe [i] >= 0 && Pe [i] < iwlen) ;
+ wnvi = wflg - nvi ;
+ for (p = Pe [i] ; p <= Pe [i] + eln - 1 ; p++)
+ {
+ e = Iw [p] ;
+ ASSERT (e >= 0 && e < n) ;
+ we = W [e] ;
+ AMD_DEBUG4 ((" e "ID" we "ID" ", e, we)) ;
+ if (we >= wflg)
+ {
+ /* unabsorbed element e has been seen in this loop */
+ AMD_DEBUG4 ((" unabsorbed, first time seen")) ;
+ we -= nvi ;
+ }
+ else if (we != 0)
+ {
+ /* e is an unabsorbed element */
+ /* this is the first we have seen e in all of Scan 1 */
+ AMD_DEBUG4 ((" unabsorbed")) ;
+ we = Degree [e] + wnvi ;
+ }
+ AMD_DEBUG4 (("\n")) ;
+ W [e] = we ;
+ }
+ }
+ }
+ AMD_DEBUG2 (("\n")) ;
+
+/* ========================================================================= */
+/* DEGREE UPDATE AND ELEMENT ABSORPTION */
+/* ========================================================================= */
+
+ /* -----------------------------------------------------------------
+ * Scan 2: for each i in Lme, sum up the degree of Lme (which is
+ * degme), plus the sum of the external degrees of each Le for the
+ * elements e appearing within i, plus the supervariables in i.
+ * Place i in hash list.
+ * ----------------------------------------------------------------- */
+
+ for (pme = pme1 ; pme <= pme2 ; pme++)
+ {
+ i = Iw [pme] ;
+ ASSERT (i >= 0 && i < n && Nv [i] < 0 && Elen [i] >= 0) ;
+ AMD_DEBUG2 (("Updating: i "ID" "ID" "ID"\n", i, Elen[i], Len [i]));
+ p1 = Pe [i] ;
+ p2 = p1 + Elen [i] - 1 ;
+ pn = p1 ;
+ hash = 0 ;
+ deg = 0 ;
+ ASSERT (p1 >= 0 && p1 < iwlen && p2 >= -1 && p2 < iwlen) ;
+
+ /* ------------------------------------------------------------- */
+ /* scan the element list associated with supervariable i */
+ /* ------------------------------------------------------------- */
+
+ /* UMFPACK/MA38-style approximate degree: */
+ if (aggressive)
+ {
+ for (p = p1 ; p <= p2 ; p++)
+ {
+ e = Iw [p] ;
+ ASSERT (e >= 0 && e < n) ;
+ we = W [e] ;
+ if (we != 0)
+ {
+ /* e is an unabsorbed element */
+ /* dext = | Le \ Lme | */
+ dext = we - wflg ;
+ if (dext > 0)
+ {
+ deg += dext ;
+ Iw [pn++] = e ;
+ hash += e ;
+ AMD_DEBUG4 ((" e: "ID" hash = "ID"\n",e,hash)) ;
+ }
+ else
+ {
+ /* external degree of e is zero, absorb e into me*/
+ AMD_DEBUG1 ((" Element "ID" =>"ID" (aggressive)\n",
+ e, me)) ;
+ ASSERT (dext == 0) ;
+ Pe [e] = FLIP (me) ;
+ W [e] = 0 ;
+ }
+ }
+ }
+ }
+ else
+ {
+ for (p = p1 ; p <= p2 ; p++)
+ {
+ e = Iw [p] ;
+ ASSERT (e >= 0 && e < n) ;
+ we = W [e] ;
+ if (we != 0)
+ {
+ /* e is an unabsorbed element */
+ dext = we - wflg ;
+ ASSERT (dext >= 0) ;
+ deg += dext ;
+ Iw [pn++] = e ;
+ hash += e ;
+ AMD_DEBUG4 ((" e: "ID" hash = "ID"\n",e,hash)) ;
+ }
+ }
+ }
+
+ /* count the number of elements in i (including me): */
+ Elen [i] = pn - p1 + 1 ;
+
+ /* ------------------------------------------------------------- */
+ /* scan the supervariables in the list associated with i */
+ /* ------------------------------------------------------------- */
+
+ /* The bulk of the AMD run time is typically spent in this loop,
+ * particularly if the matrix has many dense rows that are not
+ * removed prior to ordering. */
+ p3 = pn ;
+ p4 = p1 + Len [i] ;
+ for (p = p2 + 1 ; p < p4 ; p++)
+ {
+ j = Iw [p] ;
+ ASSERT (j >= 0 && j < n) ;
+ nvj = Nv [j] ;
+ if (nvj > 0)
+ {
+ /* j is unabsorbed, and not in Lme. */
+ /* add to degree and add to new list */
+ deg += nvj ;
+ Iw [pn++] = j ;
+ hash += j ;
+ AMD_DEBUG4 ((" s: "ID" hash "ID" Nv[j]= "ID"\n",
+ j, hash, nvj)) ;
+ }
+ }
+
+ /* ------------------------------------------------------------- */
+ /* update the degree and check for mass elimination */
+ /* ------------------------------------------------------------- */
+
+ /* with aggressive absorption, deg==0 is identical to the
+ * Elen [i] == 1 && p3 == pn test, below. */
+ ASSERT (IMPLIES (aggressive, (deg==0) == (Elen[i]==1 && p3==pn))) ;
+
+ if (Elen [i] == 1 && p3 == pn)
+ {
+
+ /* --------------------------------------------------------- */
+ /* mass elimination */
+ /* --------------------------------------------------------- */
+
+ /* There is nothing left of this node except for an edge to
+ * the current pivot element. Elen [i] is 1, and there are
+ * no variables adjacent to node i. Absorb i into the
+ * current pivot element, me. Note that if there are two or
+ * more mass eliminations, fillin due to mass elimination is
+ * possible within the nvpiv-by-nvpiv pivot block. It is this
+ * step that causes AMD's analysis to be an upper bound.
+ *
+ * The reason is that the selected pivot has a lower
+ * approximate degree than the true degree of the two mass
+ * eliminated nodes. There is no edge between the two mass
+ * eliminated nodes. They are merged with the current pivot
+ * anyway.
+ *
+ * No fillin occurs in the Schur complement, in any case,
+ * and this effect does not decrease the quality of the
+ * ordering itself, just the quality of the nonzero and
+ * flop count analysis. It also means that the post-ordering
+ * is not an exact elimination tree post-ordering. */
+
+ AMD_DEBUG1 ((" MASS i "ID" => parent e "ID"\n", i, me)) ;
+ Pe [i] = FLIP (me) ;
+ nvi = -Nv [i] ;
+ degme -= nvi ;
+ nvpiv += nvi ;
+ nel += nvi ;
+ Nv [i] = 0 ;
+ Elen [i] = EMPTY ;
+
+ }
+ else
+ {
+
+ /* --------------------------------------------------------- */
+ /* update the upper-bound degree of i */
+ /* --------------------------------------------------------- */
+
+ /* the following degree does not yet include the size
+ * of the current element, which is added later: */
+
+ Degree [i] = MIN (Degree [i], deg) ;
+
+ /* --------------------------------------------------------- */
+ /* add me to the list for i */
+ /* --------------------------------------------------------- */
+
+ /* move first supervariable to end of list */
+ Iw [pn] = Iw [p3] ;
+ /* move first element to end of element part of list */
+ Iw [p3] = Iw [p1] ;
+ /* add new element, me, to front of list. */
+ Iw [p1] = me ;
+ /* store the new length of the list in Len [i] */
+ Len [i] = pn - p1 + 1 ;
+
+ /* --------------------------------------------------------- */
+ /* place in hash bucket. Save hash key of i in Last [i]. */
+ /* --------------------------------------------------------- */
+
+ /* NOTE: this can fail if hash is negative, because the ANSI C
+ * standard does not define a % b when a and/or b are negative.
+ * That's why hash is defined as an unsigned Int, to avoid this
+ * problem. */
+ hash = hash % n ;
+ ASSERT (((Int) hash) >= 0 && ((Int) hash) < n) ;
+
+ /* if the Hhead array is not used: */
+ j = Head [hash] ;
+ if (j <= EMPTY)
+ {
+ /* degree list is empty, hash head is FLIP (j) */
+ Next [i] = FLIP (j) ;
+ Head [hash] = FLIP (i) ;
+ }
+ else
+ {
+ /* degree list is not empty, use Last [Head [hash]] as
+ * hash head. */
+ Next [i] = Last [j] ;
+ Last [j] = i ;
+ }
+
+ /* if a separate Hhead array is used: *
+ Next [i] = Hhead [hash] ;
+ Hhead [hash] = i ;
+ */
+
+ Last [i] = hash ;
+ }
+ }
+
+ Degree [me] = degme ;
+
+ /* ----------------------------------------------------------------- */
+ /* Clear the counter array, W [...], by incrementing wflg. */
+ /* ----------------------------------------------------------------- */
+
+ /* make sure that wflg+n does not cause integer overflow */
+ lemax = MAX (lemax, degme) ;
+ wflg += lemax ;
+ wflg = clear_flag (wflg, wbig, W, n) ;
+ /* at this point, W [0..n-1] < wflg holds */
+
+/* ========================================================================= */
+/* SUPERVARIABLE DETECTION */
+/* ========================================================================= */
+
+ AMD_DEBUG1 (("Detecting supervariables:\n")) ;
+ for (pme = pme1 ; pme <= pme2 ; pme++)
+ {
+ i = Iw [pme] ;
+ ASSERT (i >= 0 && i < n) ;
+ AMD_DEBUG2 (("Consider i "ID" nv "ID"\n", i, Nv [i])) ;
+ if (Nv [i] < 0)
+ {
+ /* i is a principal variable in Lme */
+
+ /* ---------------------------------------------------------
+ * examine all hash buckets with 2 or more variables. We do
+ * this by examing all unique hash keys for supervariables in
+ * the pattern Lme of the current element, me
+ * --------------------------------------------------------- */
+
+ /* let i = head of hash bucket, and empty the hash bucket */
+ ASSERT (Last [i] >= 0 && Last [i] < n) ;
+ hash = Last [i] ;
+
+ /* if Hhead array is not used: */
+ j = Head [hash] ;
+ if (j == EMPTY)
+ {
+ /* hash bucket and degree list are both empty */
+ i = EMPTY ;
+ }
+ else if (j < EMPTY)
+ {
+ /* degree list is empty */
+ i = FLIP (j) ;
+ Head [hash] = EMPTY ;
+ }
+ else
+ {
+ /* degree list is not empty, restore Last [j] of head j */
+ i = Last [j] ;
+ Last [j] = EMPTY ;
+ }
+
+ /* if separate Hhead array is used: *
+ i = Hhead [hash] ;
+ Hhead [hash] = EMPTY ;
+ */
+
+ ASSERT (i >= EMPTY && i < n) ;
+ AMD_DEBUG2 (("----i "ID" hash "ID"\n", i, hash)) ;
+
+ while (i != EMPTY && Next [i] != EMPTY)
+ {
+
+ /* -----------------------------------------------------
+ * this bucket has one or more variables following i.
+ * scan all of them to see if i can absorb any entries
+ * that follow i in hash bucket. Scatter i into w.
+ * ----------------------------------------------------- */
+
+ ln = Len [i] ;
+ eln = Elen [i] ;
+ ASSERT (ln >= 0 && eln >= 0) ;
+ ASSERT (Pe [i] >= 0 && Pe [i] < iwlen) ;
+ /* do not flag the first element in the list (me) */
+ for (p = Pe [i] + 1 ; p <= Pe [i] + ln - 1 ; p++)
+ {
+ ASSERT (Iw [p] >= 0 && Iw [p] < n) ;
+ W [Iw [p]] = wflg ;
+ }
+
+ /* ----------------------------------------------------- */
+ /* scan every other entry j following i in bucket */
+ /* ----------------------------------------------------- */
+
+ jlast = i ;
+ j = Next [i] ;
+ ASSERT (j >= EMPTY && j < n) ;
+
+ while (j != EMPTY)
+ {
+ /* ------------------------------------------------- */
+ /* check if j and i have identical nonzero pattern */
+ /* ------------------------------------------------- */
+
+ AMD_DEBUG3 (("compare i "ID" and j "ID"\n", i,j)) ;
+
+ /* check if i and j have the same Len and Elen */
+ ASSERT (Len [j] >= 0 && Elen [j] >= 0) ;
+ ASSERT (Pe [j] >= 0 && Pe [j] < iwlen) ;
+ ok = (Len [j] == ln) && (Elen [j] == eln) ;
+ /* skip the first element in the list (me) */
+ for (p = Pe [j] + 1 ; ok && p <= Pe [j] + ln - 1 ; p++)
+ {
+ ASSERT (Iw [p] >= 0 && Iw [p] < n) ;
+ if (W [Iw [p]] != wflg) ok = 0 ;
+ }
+ if (ok)
+ {
+ /* --------------------------------------------- */
+ /* found it! j can be absorbed into i */
+ /* --------------------------------------------- */
+
+ AMD_DEBUG1 (("found it! j "ID" => i "ID"\n", j,i));
+ Pe [j] = FLIP (i) ;
+ /* both Nv [i] and Nv [j] are negated since they */
+ /* are in Lme, and the absolute values of each */
+ /* are the number of variables in i and j: */
+ Nv [i] += Nv [j] ;
+ Nv [j] = 0 ;
+ Elen [j] = EMPTY ;
+ /* delete j from hash bucket */
+ ASSERT (j != Next [j]) ;
+ j = Next [j] ;
+ Next [jlast] = j ;
+
+ }
+ else
+ {
+ /* j cannot be absorbed into i */
+ jlast = j ;
+ ASSERT (j != Next [j]) ;
+ j = Next [j] ;
+ }
+ ASSERT (j >= EMPTY && j < n) ;
+ }
+
+ /* -----------------------------------------------------
+ * no more variables can be absorbed into i
+ * go to next i in bucket and clear flag array
+ * ----------------------------------------------------- */
+
+ wflg++ ;
+ i = Next [i] ;
+ ASSERT (i >= EMPTY && i < n) ;
+
+ }
+ }
+ }
+ AMD_DEBUG2 (("detect done\n")) ;
+
+/* ========================================================================= */
+/* RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVARIABLES FROM ELEMENT */
+/* ========================================================================= */
+
+ p = pme1 ;
+ nleft = n - nel ;
+ for (pme = pme1 ; pme <= pme2 ; pme++)
+ {
+ i = Iw [pme] ;
+ ASSERT (i >= 0 && i < n) ;
+ nvi = -Nv [i] ;
+ AMD_DEBUG3 (("Restore i "ID" "ID"\n", i, nvi)) ;
+ if (nvi > 0)
+ {
+ /* i is a principal variable in Lme */
+ /* restore Nv [i] to signify that i is principal */
+ Nv [i] = nvi ;
+
+ /* --------------------------------------------------------- */
+ /* compute the external degree (add size of current element) */
+ /* --------------------------------------------------------- */
+
+ deg = Degree [i] + degme - nvi ;
+ deg = MIN (deg, nleft - nvi) ;
+ ASSERT (IMPLIES (aggressive, deg > 0) && deg >= 0 && deg < n) ;
+
+ /* --------------------------------------------------------- */
+ /* place the supervariable at the head of the degree list */
+ /* --------------------------------------------------------- */
+
+ inext = Head [deg] ;
+ ASSERT (inext >= EMPTY && inext < n) ;
+ if (inext != EMPTY) Last [inext] = i ;
+ Next [i] = inext ;
+ Last [i] = EMPTY ;
+ Head [deg] = i ;
+
+ /* --------------------------------------------------------- */
+ /* save the new degree, and find the minimum degree */
+ /* --------------------------------------------------------- */
+
+ mindeg = MIN (mindeg, deg) ;
+ Degree [i] = deg ;
+
+ /* --------------------------------------------------------- */
+ /* place the supervariable in the element pattern */
+ /* --------------------------------------------------------- */
+
+ Iw [p++] = i ;
+
+ }
+ }
+ AMD_DEBUG2 (("restore done\n")) ;
+
+/* ========================================================================= */
+/* FINALIZE THE NEW ELEMENT */
+/* ========================================================================= */
+
+ AMD_DEBUG2 (("ME = "ID" DONE\n", me)) ;
+ Nv [me] = nvpiv ;
+ /* save the length of the list for the new element me */
+ Len [me] = p - pme1 ;
+ if (Len [me] == 0)
+ {
+ /* there is nothing left of the current pivot element */
+ /* it is a root of the assembly tree */
+ Pe [me] = EMPTY ;
+ W [me] = 0 ;
+ }
+ if (elenme != 0)
+ {
+ /* element was not constructed in place: deallocate part of */
+ /* it since newly nonprincipal variables may have been removed */
+ pfree = p ;
+ }
+
+ /* The new element has nvpiv pivots and the size of the contribution
+ * block for a multifrontal method is degme-by-degme, not including
+ * the "dense" rows/columns. If the "dense" rows/columns are included,
+ * the frontal matrix is no larger than
+ * (degme+ndense)-by-(degme+ndense).
+ */
+
+ if (Info != (c_float *) NULL)
+ {
+ f = nvpiv ;
+ r = degme + ndense ;
+ dmax = MAX (dmax, f + r) ;
+
+ /* number of nonzeros in L (excluding the diagonal) */
+ lnzme = f*r + (f-1)*f/2 ;
+ lnz += lnzme ;
+
+ /* number of divide operations for LDL' and for LU */
+ ndiv += lnzme ;
+
+ /* number of multiply-subtract pairs for LU */
+ s = f*r*r + r*(f-1)*f + (f-1)*f*(2*f-1)/6 ;
+ nms_lu += s ;
+
+ /* number of multiply-subtract pairs for LDL' */
+ nms_ldl += (s + lnzme)/2 ;
+ }
+
+#ifndef NDEBUG
+ AMD_DEBUG2 (("finalize done nel "ID" n "ID"\n ::::\n", nel, n)) ;
+ for (pme = Pe [me] ; pme <= Pe [me] + Len [me] - 1 ; pme++)
+ {
+ AMD_DEBUG3 ((" "ID"", Iw [pme])) ;
+ }
+ AMD_DEBUG3 (("\n")) ;
+#endif
+
+ }
+
+/* ========================================================================= */
+/* DONE SELECTING PIVOTS */
+/* ========================================================================= */
+
+ if (Info != (c_float *) NULL)
+ {
+
+ /* count the work to factorize the ndense-by-ndense submatrix */
+ f = ndense ;
+ dmax = MAX (dmax, (c_float) ndense) ;
+
+ /* number of nonzeros in L (excluding the diagonal) */
+ lnzme = (f-1)*f/2 ;
+ lnz += lnzme ;
+
+ /* number of divide operations for LDL' and for LU */
+ ndiv += lnzme ;
+
+ /* number of multiply-subtract pairs for LU */
+ s = (f-1)*f*(2*f-1)/6 ;
+ nms_lu += s ;
+
+ /* number of multiply-subtract pairs for LDL' */
+ nms_ldl += (s + lnzme)/2 ;
+
+ /* number of nz's in L (excl. diagonal) */
+ Info [AMD_LNZ] = lnz ;
+
+ /* number of divide ops for LU and LDL' */
+ Info [AMD_NDIV] = ndiv ;
+
+ /* number of multiply-subtract pairs for LDL' */
+ Info [AMD_NMULTSUBS_LDL] = nms_ldl ;
+
+ /* number of multiply-subtract pairs for LU */
+ Info [AMD_NMULTSUBS_LU] = nms_lu ;
+
+ /* number of "dense" rows/columns */
+ Info [AMD_NDENSE] = ndense ;
+
+ /* largest front is dmax-by-dmax */
+ Info [AMD_DMAX] = dmax ;
+
+ /* number of garbage collections in AMD */
+ Info [AMD_NCMPA] = ncmpa ;
+
+ /* successful ordering */
+ Info [AMD_STATUS] = AMD_OK ;
+ }
+
+/* ========================================================================= */
+/* POST-ORDERING */
+/* ========================================================================= */
+
+/* -------------------------------------------------------------------------
+ * Variables at this point:
+ *
+ * Pe: holds the elimination tree. The parent of j is FLIP (Pe [j]),
+ * or EMPTY if j is a root. The tree holds both elements and
+ * non-principal (unordered) variables absorbed into them.
+ * Dense variables are non-principal and unordered.
+ *
+ * Elen: holds the size of each element, including the diagonal part.
+ * FLIP (Elen [e]) > 0 if e is an element. For unordered
+ * variables i, Elen [i] is EMPTY.
+ *
+ * Nv: Nv [e] > 0 is the number of pivots represented by the element e.
+ * For unordered variables i, Nv [i] is zero.
+ *
+ * Contents no longer needed:
+ * W, Iw, Len, Degree, Head, Next, Last.
+ *
+ * The matrix itself has been destroyed.
+ *
+ * n: the size of the matrix.
+ * No other scalars needed (pfree, iwlen, etc.)
+ * ------------------------------------------------------------------------- */
+
+ /* restore Pe */
+ for (i = 0 ; i < n ; i++)
+ {
+ Pe [i] = FLIP (Pe [i]) ;
+ }
+
+ /* restore Elen, for output information, and for postordering */
+ for (i = 0 ; i < n ; i++)
+ {
+ Elen [i] = FLIP (Elen [i]) ;
+ }
+
+/* Now the parent of j is Pe [j], or EMPTY if j is a root. Elen [e] > 0
+ * is the size of element e. Elen [i] is EMPTY for unordered variable i. */
+
+#ifndef NDEBUG
+ AMD_DEBUG2 (("\nTree:\n")) ;
+ for (i = 0 ; i < n ; i++)
+ {
+ AMD_DEBUG2 ((" "ID" parent: "ID" ", i, Pe [i])) ;
+ ASSERT (Pe [i] >= EMPTY && Pe [i] < n) ;
+ if (Nv [i] > 0)
+ {
+ /* this is an element */
+ e = i ;
+ AMD_DEBUG2 ((" element, size is "ID"\n", Elen [i])) ;
+ ASSERT (Elen [e] > 0) ;
+ }
+ AMD_DEBUG2 (("\n")) ;
+ }
+ AMD_DEBUG2 (("\nelements:\n")) ;
+ for (e = 0 ; e < n ; e++)
+ {
+ if (Nv [e] > 0)
+ {
+ AMD_DEBUG3 (("Element e= "ID" size "ID" nv "ID" \n", e,
+ Elen [e], Nv [e])) ;
+ }
+ }
+ AMD_DEBUG2 (("\nvariables:\n")) ;
+ for (i = 0 ; i < n ; i++)
+ {
+ Int cnt ;
+ if (Nv [i] == 0)
+ {
+ AMD_DEBUG3 (("i unordered: "ID"\n", i)) ;
+ j = Pe [i] ;
+ cnt = 0 ;
+ AMD_DEBUG3 ((" j: "ID"\n", j)) ;
+ if (j == EMPTY)
+ {
+ AMD_DEBUG3 ((" i is a dense variable\n")) ;
+ }
+ else
+ {
+ ASSERT (j >= 0 && j < n) ;
+ while (Nv [j] == 0)
+ {
+ AMD_DEBUG3 ((" j : "ID"\n", j)) ;
+ j = Pe [j] ;
+ AMD_DEBUG3 ((" j:: "ID"\n", j)) ;
+ cnt++ ;
+ if (cnt > n) break ;
+ }
+ e = j ;
+ AMD_DEBUG3 ((" got to e: "ID"\n", e)) ;
+ }
+ }
+ }
+#endif
+
+/* ========================================================================= */
+/* compress the paths of the variables */
+/* ========================================================================= */
+
+ for (i = 0 ; i < n ; i++)
+ {
+ if (Nv [i] == 0)
+ {
+
+ /* -------------------------------------------------------------
+ * i is an un-ordered row. Traverse the tree from i until
+ * reaching an element, e. The element, e, was the principal
+ * supervariable of i and all nodes in the path from i to when e
+ * was selected as pivot.
+ * ------------------------------------------------------------- */
+
+ AMD_DEBUG1 (("Path compression, i unordered: "ID"\n", i)) ;
+ j = Pe [i] ;
+ ASSERT (j >= EMPTY && j < n) ;
+ AMD_DEBUG3 ((" j: "ID"\n", j)) ;
+ if (j == EMPTY)
+ {
+ /* Skip a dense variable. It has no parent. */
+ AMD_DEBUG3 ((" i is a dense variable\n")) ;
+ continue ;
+ }
+
+ /* while (j is a variable) */
+ while (Nv [j] == 0)
+ {
+ AMD_DEBUG3 ((" j : "ID"\n", j)) ;
+ j = Pe [j] ;
+ AMD_DEBUG3 ((" j:: "ID"\n", j)) ;
+ ASSERT (j >= 0 && j < n) ;
+ }
+ /* got to an element e */
+ e = j ;
+ AMD_DEBUG3 (("got to e: "ID"\n", e)) ;
+
+ /* -------------------------------------------------------------
+ * traverse the path again from i to e, and compress the path
+ * (all nodes point to e). Path compression allows this code to
+ * compute in O(n) time.
+ * ------------------------------------------------------------- */
+
+ j = i ;
+ /* while (j is a variable) */
+ while (Nv [j] == 0)
+ {
+ jnext = Pe [j] ;
+ AMD_DEBUG3 (("j "ID" jnext "ID"\n", j, jnext)) ;
+ Pe [j] = e ;
+ j = jnext ;
+ ASSERT (j >= 0 && j < n) ;
+ }
+ }
+ }
+
+/* ========================================================================= */
+/* postorder the assembly tree */
+/* ========================================================================= */
+
+ AMD_postorder (n, Pe, Nv, Elen,
+ W, /* output order */
+ Head, Next, Last) ; /* workspace */
+
+/* ========================================================================= */
+/* compute output permutation and inverse permutation */
+/* ========================================================================= */
+
+ /* W [e] = k means that element e is the kth element in the new
+ * order. e is in the range 0 to n-1, and k is in the range 0 to
+ * the number of elements. Use Head for inverse order. */
+
+ for (k = 0 ; k < n ; k++)
+ {
+ Head [k] = EMPTY ;
+ Next [k] = EMPTY ;
+ }
+ for (e = 0 ; e < n ; e++)
+ {
+ k = W [e] ;
+ ASSERT ((k == EMPTY) == (Nv [e] == 0)) ;
+ if (k != EMPTY)
+ {
+ ASSERT (k >= 0 && k < n) ;
+ Head [k] = e ;
+ }
+ }
+
+ /* construct output inverse permutation in Next,
+ * and permutation in Last */
+ nel = 0 ;
+ for (k = 0 ; k < n ; k++)
+ {
+ e = Head [k] ;
+ if (e == EMPTY) break ;
+ ASSERT (e >= 0 && e < n && Nv [e] > 0) ;
+ Next [e] = nel ;
+ nel += Nv [e] ;
+ }
+ ASSERT (nel == n - ndense) ;
+
+ /* order non-principal variables (dense, & those merged into supervar's) */
+ for (i = 0 ; i < n ; i++)
+ {
+ if (Nv [i] == 0)
+ {
+ e = Pe [i] ;
+ ASSERT (e >= EMPTY && e < n) ;
+ if (e != EMPTY)
+ {
+ /* This is an unordered variable that was merged
+ * into element e via supernode detection or mass
+ * elimination of i when e became the pivot element.
+ * Place i in order just before e. */
+ ASSERT (Next [i] == EMPTY && Nv [e] > 0) ;
+ Next [i] = Next [e] ;
+ Next [e]++ ;
+ }
+ else
+ {
+ /* This is a dense unordered variable, with no parent.
+ * Place it last in the output order. */
+ Next [i] = nel++ ;
+ }
+ }
+ }
+ ASSERT (nel == n) ;
+
+ AMD_DEBUG2 (("\n\nPerm:\n")) ;
+ for (i = 0 ; i < n ; i++)
+ {
+ k = Next [i] ;
+ ASSERT (k >= 0 && k < n) ;
+ Last [k] = i ;
+ AMD_DEBUG2 ((" perm ["ID"] = "ID"\n", k, i)) ;
+ }
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_aat.c b/lin_sys/direct/qdldl/amd/src/amd_aat.c
new file mode 100644
index 0000000..cccb617
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_aat.c
@@ -0,0 +1,184 @@
+/* ========================================================================= */
+/* === AMD_aat ============================================================= */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* AMD_aat: compute the symmetry of the pattern of A, and count the number of
+ * nonzeros each column of A+A' (excluding the diagonal). Assumes the input
+ * matrix has no errors, with sorted columns and no duplicates
+ * (AMD_valid (n, n, Ap, Ai) must be AMD_OK, but this condition is not
+ * checked).
+ */
+
+#include "amd_internal.h"
+
+GLOBAL size_t AMD_aat /* returns nz in A+A' */
+(
+ Int n,
+ const Int Ap [ ],
+ const Int Ai [ ],
+ Int Len [ ], /* Len [j]: length of column j of A+A', excl diagonal*/
+ Int Tp [ ], /* workspace of size n */
+ c_float Info [ ]
+)
+{
+ Int p1, p2, p, i, j, pj, pj2, k, nzdiag, nzboth, nz ;
+ c_float sym ;
+ size_t nzaat ;
+
+#ifndef NDEBUG
+ AMD_debug_init ("AMD AAT") ;
+ for (k = 0 ; k < n ; k++) Tp [k] = EMPTY ;
+ ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
+#endif
+
+ if (Info != (c_float *) NULL)
+ {
+ /* clear the Info array, if it exists */
+ for (i = 0 ; i < AMD_INFO ; i++)
+ {
+ Info [i] = EMPTY ;
+ }
+ Info [AMD_STATUS] = AMD_OK ;
+ }
+
+ for (k = 0 ; k < n ; k++)
+ {
+ Len [k] = 0 ;
+ }
+
+ nzdiag = 0 ;
+ nzboth = 0 ;
+ nz = Ap [n] ;
+
+ for (k = 0 ; k < n ; k++)
+ {
+ p1 = Ap [k] ;
+ p2 = Ap [k+1] ;
+ AMD_DEBUG2 (("\nAAT Column: "ID" p1: "ID" p2: "ID"\n", k, p1, p2)) ;
+
+ /* construct A+A' */
+ for (p = p1 ; p < p2 ; )
+ {
+ /* scan the upper triangular part of A */
+ j = Ai [p] ;
+ if (j < k)
+ {
+ /* entry A (j,k) is in the strictly upper triangular part,
+ * add both A (j,k) and A (k,j) to the matrix A+A' */
+ Len [j]++ ;
+ Len [k]++ ;
+ AMD_DEBUG3 ((" upper ("ID","ID") ("ID","ID")\n", j,k, k,j));
+ p++ ;
+ }
+ else if (j == k)
+ {
+ /* skip the diagonal */
+ p++ ;
+ nzdiag++ ;
+ break ;
+ }
+ else /* j > k */
+ {
+ /* first entry below the diagonal */
+ break ;
+ }
+ /* scan lower triangular part of A, in column j until reaching
+ * row k. Start where last scan left off. */
+ ASSERT (Tp [j] != EMPTY) ;
+ ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;
+ pj2 = Ap [j+1] ;
+ for (pj = Tp [j] ; pj < pj2 ; )
+ {
+ i = Ai [pj] ;
+ if (i < k)
+ {
+ /* A (i,j) is only in the lower part, not in upper.
+ * add both A (i,j) and A (j,i) to the matrix A+A' */
+ Len [i]++ ;
+ Len [j]++ ;
+ AMD_DEBUG3 ((" lower ("ID","ID") ("ID","ID")\n",
+ i,j, j,i)) ;
+ pj++ ;
+ }
+ else if (i == k)
+ {
+ /* entry A (k,j) in lower part and A (j,k) in upper */
+ pj++ ;
+ nzboth++ ;
+ break ;
+ }
+ else /* i > k */
+ {
+ /* consider this entry later, when k advances to i */
+ break ;
+ }
+ }
+ Tp [j] = pj ;
+ }
+ /* Tp [k] points to the entry just below the diagonal in column k */
+ Tp [k] = p ;
+ }
+
+ /* clean up, for remaining mismatched entries */
+ for (j = 0 ; j < n ; j++)
+ {
+ for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)
+ {
+ i = Ai [pj] ;
+ /* A (i,j) is only in the lower part, not in upper.
+ * add both A (i,j) and A (j,i) to the matrix A+A' */
+ Len [i]++ ;
+ Len [j]++ ;
+ AMD_DEBUG3 ((" lower cleanup ("ID","ID") ("ID","ID")\n",
+ i,j, j,i)) ;
+ }
+ }
+
+ /* --------------------------------------------------------------------- */
+ /* compute the symmetry of the nonzero pattern of A */
+ /* --------------------------------------------------------------------- */
+
+ /* Given a matrix A, the symmetry of A is:
+ * B = tril (spones (A), -1) + triu (spones (A), 1) ;
+ * sym = nnz (B & B') / nnz (B) ;
+ * or 1 if nnz (B) is zero.
+ */
+
+ if (nz == nzdiag)
+ {
+ sym = 1 ;
+ }
+ else
+ {
+ sym = (2 * (c_float) nzboth) / ((c_float) (nz - nzdiag)) ;
+ }
+
+ nzaat = 0 ;
+ for (k = 0 ; k < n ; k++)
+ {
+ nzaat += Len [k] ;
+ }
+
+ AMD_DEBUG1 (("AMD nz in A+A', excluding diagonal (nzaat) = %g\n",
+ (c_float) nzaat)) ;
+ AMD_DEBUG1 ((" nzboth: "ID" nz: "ID" nzdiag: "ID" symmetry: %g\n",
+ nzboth, nz, nzdiag, sym)) ;
+
+ if (Info != (c_float *) NULL)
+ {
+ Info [AMD_STATUS] = AMD_OK ;
+ Info [AMD_N] = n ;
+ Info [AMD_NZ] = nz ;
+ Info [AMD_SYMMETRY] = sym ; /* symmetry of pattern of A */
+ Info [AMD_NZDIAG] = nzdiag ; /* nonzeros on diagonal of A */
+ Info [AMD_NZ_A_PLUS_AT] = nzaat ; /* nonzeros in A+A' */
+ }
+
+ return (nzaat) ;
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_control.c b/lin_sys/direct/qdldl/amd/src/amd_control.c
new file mode 100644
index 0000000..0cfd298
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_control.c
@@ -0,0 +1,64 @@
+/* ========================================================================= */
+/* === AMD_control ========================================================= */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* User-callable. Prints the control parameters for AMD. See amd.h
+ * for details. If the Control array is not present, the defaults are
+ * printed instead.
+ */
+
+#include "amd_internal.h"
+
+GLOBAL void AMD_control
+(
+ c_float Control [ ]
+)
+{
+ c_float alpha ;
+ Int aggressive ;
+
+ if (Control != (c_float *) NULL)
+ {
+ alpha = Control [AMD_DENSE] ;
+ aggressive = Control [AMD_AGGRESSIVE] != 0 ;
+ }
+ else
+ {
+ alpha = AMD_DEFAULT_DENSE ;
+ aggressive = AMD_DEFAULT_AGGRESSIVE ;
+ }
+
+ SUITESPARSE_PRINTF ((
+ "\nAMD version %d.%d.%d, %s: approximate minimum degree ordering\n"
+ " dense row parameter: %g\n", AMD_MAIN_VERSION, AMD_SUB_VERSION,
+ AMD_SUBSUB_VERSION, AMD_DATE, alpha)) ;
+
+ if (alpha < 0)
+ {
+ SUITESPARSE_PRINTF ((" no rows treated as dense\n")) ;
+ }
+ else
+ {
+ SUITESPARSE_PRINTF ((
+ " (rows with more than max (%g * sqrt (n), 16) entries are\n"
+ " considered \"dense\", and placed last in output permutation)\n",
+ alpha)) ;
+ }
+
+ if (aggressive)
+ {
+ SUITESPARSE_PRINTF ((" aggressive absorption: yes\n")) ;
+ }
+ else
+ {
+ SUITESPARSE_PRINTF ((" aggressive absorption: no\n")) ;
+ }
+
+ SUITESPARSE_PRINTF ((" size of AMD integer: %d\n\n", sizeof (Int))) ;
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_defaults.c b/lin_sys/direct/qdldl/amd/src/amd_defaults.c
new file mode 100644
index 0000000..50112f1
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_defaults.c
@@ -0,0 +1,37 @@
+/* ========================================================================= */
+/* === AMD_defaults ======================================================== */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* User-callable. Sets default control parameters for AMD. See amd.h
+ * for details.
+ */
+
+#include "amd_internal.h"
+
+/* ========================================================================= */
+/* === AMD defaults ======================================================== */
+/* ========================================================================= */
+
+GLOBAL void AMD_defaults
+(
+ c_float Control [ ]
+)
+{
+ Int i ;
+
+ if (Control != (c_float *) NULL)
+ {
+ for (i = 0 ; i < AMD_CONTROL ; i++)
+ {
+ Control [i] = 0 ;
+ }
+ Control [AMD_DENSE] = AMD_DEFAULT_DENSE ;
+ Control [AMD_AGGRESSIVE] = AMD_DEFAULT_AGGRESSIVE ;
+ }
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_info.c b/lin_sys/direct/qdldl/amd/src/amd_info.c
new file mode 100644
index 0000000..bf39378
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_info.c
@@ -0,0 +1,119 @@
+/* ========================================================================= */
+/* === AMD_info ============================================================ */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* User-callable. Prints the output statistics for AMD. See amd.h
+ * for details. If the Info array is not present, nothing is printed.
+ */
+
+#include "amd_internal.h"
+
+#define PRI(format,x) { if (x >= 0) { SUITESPARSE_PRINTF ((format, x)) ; }}
+
+GLOBAL void AMD_info
+(
+ c_float Info [ ]
+)
+{
+ c_float n, ndiv, nmultsubs_ldl, nmultsubs_lu, lnz, lnzd ;
+
+ SUITESPARSE_PRINTF (("\nAMD version %d.%d.%d, %s, results:\n",
+ AMD_MAIN_VERSION, AMD_SUB_VERSION, AMD_SUBSUB_VERSION, AMD_DATE)) ;
+
+ if (!Info)
+ {
+ return ;
+ }
+
+ n = Info [AMD_N] ;
+ ndiv = Info [AMD_NDIV] ;
+ nmultsubs_ldl = Info [AMD_NMULTSUBS_LDL] ;
+ nmultsubs_lu = Info [AMD_NMULTSUBS_LU] ;
+ lnz = Info [AMD_LNZ] ;
+ lnzd = (n >= 0 && lnz >= 0) ? (n + lnz) : (-1) ;
+
+ /* AMD return status */
+ SUITESPARSE_PRINTF ((" status: ")) ;
+ if (Info [AMD_STATUS] == AMD_OK)
+ {
+ SUITESPARSE_PRINTF (("OK\n")) ;
+ }
+ else if (Info [AMD_STATUS] == AMD_OUT_OF_MEMORY)
+ {
+ SUITESPARSE_PRINTF (("out of memory\n")) ;
+ }
+ else if (Info [AMD_STATUS] == AMD_INVALID)
+ {
+ SUITESPARSE_PRINTF (("invalid matrix\n")) ;
+ }
+ else if (Info [AMD_STATUS] == AMD_OK_BUT_JUMBLED)
+ {
+ SUITESPARSE_PRINTF (("OK, but jumbled\n")) ;
+ }
+ else
+ {
+ SUITESPARSE_PRINTF (("unknown\n")) ;
+ }
+
+ /* statistics about the input matrix */
+ PRI (" n, dimension of A: %.20g\n", n);
+ PRI (" nz, number of nonzeros in A: %.20g\n",
+ Info [AMD_NZ]) ;
+ PRI (" symmetry of A: %.4f\n",
+ Info [AMD_SYMMETRY]) ;
+ PRI (" number of nonzeros on diagonal: %.20g\n",
+ Info [AMD_NZDIAG]) ;
+ PRI (" nonzeros in pattern of A+A' (excl. diagonal): %.20g\n",
+ Info [AMD_NZ_A_PLUS_AT]) ;
+ PRI (" # dense rows/columns of A+A': %.20g\n",
+ Info [AMD_NDENSE]) ;
+
+ /* statistics about AMD's behavior */
+ PRI (" memory used, in bytes: %.20g\n",
+ Info [AMD_MEMORY]) ;
+ PRI (" # of memory compactions: %.20g\n",
+ Info [AMD_NCMPA]) ;
+
+ /* statistics about the ordering quality */
+ SUITESPARSE_PRINTF (("\n"
+ " The following approximate statistics are for a subsequent\n"
+ " factorization of A(P,P) + A(P,P)'. They are slight upper\n"
+ " bounds if there are no dense rows/columns in A+A', and become\n"
+ " looser if dense rows/columns exist.\n\n")) ;
+
+ PRI (" nonzeros in L (excluding diagonal): %.20g\n",
+ lnz) ;
+ PRI (" nonzeros in L (including diagonal): %.20g\n",
+ lnzd) ;
+ PRI (" # divide operations for LDL' or LU: %.20g\n",
+ ndiv) ;
+ PRI (" # multiply-subtract operations for LDL': %.20g\n",
+ nmultsubs_ldl) ;
+ PRI (" # multiply-subtract operations for LU: %.20g\n",
+ nmultsubs_lu) ;
+ PRI (" max nz. in any column of L (incl. diagonal): %.20g\n",
+ Info [AMD_DMAX]) ;
+
+ /* total flop counts for various factorizations */
+
+ if (n >= 0 && ndiv >= 0 && nmultsubs_ldl >= 0 && nmultsubs_lu >= 0)
+ {
+ SUITESPARSE_PRINTF (("\n"
+ " chol flop count for real A, sqrt counted as 1 flop: %.20g\n"
+ " LDL' flop count for real A: %.20g\n"
+ " LDL' flop count for complex A: %.20g\n"
+ " LU flop count for real A (with no pivoting): %.20g\n"
+ " LU flop count for complex A (with no pivoting): %.20g\n\n",
+ n + ndiv + 2*nmultsubs_ldl,
+ ndiv + 2*nmultsubs_ldl,
+ 9*ndiv + 8*nmultsubs_ldl,
+ ndiv + 2*nmultsubs_lu,
+ 9*ndiv + 8*nmultsubs_lu)) ;
+ }
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_order.c b/lin_sys/direct/qdldl/amd/src/amd_order.c
new file mode 100644
index 0000000..8b635a9
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_order.c
@@ -0,0 +1,199 @@
+/* ========================================================================= */
+/* === AMD_order =========================================================== */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* User-callable AMD minimum degree ordering routine. See amd.h for
+ * documentation.
+ */
+
+#include "amd_internal.h"
+
+/* ========================================================================= */
+/* === AMD_order =========================================================== */
+/* ========================================================================= */
+
+GLOBAL Int AMD_order
+(
+ Int n,
+ const Int Ap [ ],
+ const Int Ai [ ],
+ Int P [ ],
+ c_float Control [ ],
+ c_float Info [ ]
+)
+{
+ Int *Len, *S, nz, i, *Pinv, info, status, *Rp, *Ri, *Cp, *Ci, ok ;
+ size_t nzaat, slen ;
+ c_float mem = 0 ;
+
+#ifndef NDEBUG
+ AMD_debug_init ("amd") ;
+#endif
+
+ /* clear the Info array, if it exists */
+ info = Info != (c_float *) NULL ;
+ if (info)
+ {
+ for (i = 0 ; i < AMD_INFO ; i++)
+ {
+ Info [i] = EMPTY ;
+ }
+ Info [AMD_N] = n ;
+ Info [AMD_STATUS] = AMD_OK ;
+ }
+
+ /* make sure inputs exist and n is >= 0 */
+ if (Ai == (Int *) NULL || Ap == (Int *) NULL || P == (Int *) NULL || n < 0)
+ {
+ if (info) Info [AMD_STATUS] = AMD_INVALID ;
+ return (AMD_INVALID) ; /* arguments are invalid */
+ }
+
+ if (n == 0)
+ {
+ return (AMD_OK) ; /* n is 0 so there's nothing to do */
+ }
+
+ nz = Ap [n] ;
+ if (info)
+ {
+ Info [AMD_NZ] = nz ;
+ }
+ if (nz < 0)
+ {
+ if (info) Info [AMD_STATUS] = AMD_INVALID ;
+ return (AMD_INVALID) ;
+ }
+
+ /* check if n or nz will cause size_t overflow */
+ if (((size_t) n) >= SIZE_T_MAX / sizeof (Int)
+ || ((size_t) nz) >= SIZE_T_MAX / sizeof (Int))
+ {
+ if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
+ return (AMD_OUT_OF_MEMORY) ; /* problem too large */
+ }
+
+ /* check the input matrix: AMD_OK, AMD_INVALID, or AMD_OK_BUT_JUMBLED */
+ status = AMD_valid (n, n, Ap, Ai) ;
+
+ if (status == AMD_INVALID)
+ {
+ if (info) Info [AMD_STATUS] = AMD_INVALID ;
+ return (AMD_INVALID) ; /* matrix is invalid */
+ }
+
+ /* allocate two size-n integer workspaces */
+ Len = SuiteSparse_malloc (n, sizeof (Int)) ;
+ Pinv = SuiteSparse_malloc (n, sizeof (Int)) ;
+ mem += n ;
+ mem += n ;
+ if (!Len || !Pinv)
+ {
+ /* :: out of memory :: */
+ SuiteSparse_free (Len) ;
+ SuiteSparse_free (Pinv) ;
+ if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
+ return (AMD_OUT_OF_MEMORY) ;
+ }
+
+ if (status == AMD_OK_BUT_JUMBLED)
+ {
+ /* sort the input matrix and remove duplicate entries */
+ AMD_DEBUG1 (("Matrix is jumbled\n")) ;
+ Rp = SuiteSparse_malloc (n+1, sizeof (Int)) ;
+ Ri = SuiteSparse_malloc (nz, sizeof (Int)) ;
+ mem += (n+1) ;
+ mem += MAX (nz,1) ;
+ if (!Rp || !Ri)
+ {
+ /* :: out of memory :: */
+ SuiteSparse_free (Rp) ;
+ SuiteSparse_free (Ri) ;
+ SuiteSparse_free (Len) ;
+ SuiteSparse_free (Pinv) ;
+ if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
+ return (AMD_OUT_OF_MEMORY) ;
+ }
+ /* use Len and Pinv as workspace to create R = A' */
+ AMD_preprocess (n, Ap, Ai, Rp, Ri, Len, Pinv) ;
+ Cp = Rp ;
+ Ci = Ri ;
+ }
+ else
+ {
+ /* order the input matrix as-is. No need to compute R = A' first */
+ Rp = NULL ;
+ Ri = NULL ;
+ Cp = (Int *) Ap ;
+ Ci = (Int *) Ai ;
+ }
+
+ /* --------------------------------------------------------------------- */
+ /* determine the symmetry and count off-diagonal nonzeros in A+A' */
+ /* --------------------------------------------------------------------- */
+
+ nzaat = AMD_aat (n, Cp, Ci, Len, P, Info) ;
+ AMD_DEBUG1 (("nzaat: %g\n", (c_float) nzaat)) ;
+ ASSERT ((MAX (nz-n, 0) <= nzaat) && (nzaat <= 2 * (size_t) nz)) ;
+
+ /* --------------------------------------------------------------------- */
+ /* allocate workspace for matrix, elbow room, and 6 size-n vectors */
+ /* --------------------------------------------------------------------- */
+
+ S = NULL ;
+ slen = nzaat ; /* space for matrix */
+ ok = ((slen + nzaat/5) >= slen) ; /* check for size_t overflow */
+ slen += nzaat/5 ; /* add elbow room */
+ for (i = 0 ; ok && i < 7 ; i++)
+ {
+ ok = ((slen + n) > slen) ; /* check for size_t overflow */
+ slen += n ; /* size-n elbow room, 6 size-n work */
+ }
+ mem += slen ;
+ ok = ok && (slen < SIZE_T_MAX / sizeof (Int)) ; /* check for overflow */
+ ok = ok && (slen < Int_MAX) ; /* S[i] for Int i must be OK */
+ if (ok)
+ {
+ S = SuiteSparse_malloc (slen, sizeof (Int)) ;
+ }
+ AMD_DEBUG1 (("slen %g\n", (c_float) slen)) ;
+ if (!S)
+ {
+ /* :: out of memory :: (or problem too large) */
+ SuiteSparse_free (Rp) ;
+ SuiteSparse_free (Ri) ;
+ SuiteSparse_free (Len) ;
+ SuiteSparse_free (Pinv) ;
+ if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
+ return (AMD_OUT_OF_MEMORY) ;
+ }
+ if (info)
+ {
+ /* memory usage, in bytes. */
+ Info [AMD_MEMORY] = mem * sizeof (Int) ;
+ }
+
+ /* --------------------------------------------------------------------- */
+ /* order the matrix */
+ /* --------------------------------------------------------------------- */
+
+ AMD_1 (n, Cp, Ci, P, Pinv, Len, slen, S, Control, Info) ;
+
+ /* --------------------------------------------------------------------- */
+ /* free the workspace */
+ /* --------------------------------------------------------------------- */
+
+ SuiteSparse_free (Rp) ;
+ SuiteSparse_free (Ri) ;
+ SuiteSparse_free (Len) ;
+ SuiteSparse_free (Pinv) ;
+ SuiteSparse_free (S) ;
+ if (info) Info [AMD_STATUS] = status ;
+ return (status) ; /* successful ordering */
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_post_tree.c b/lin_sys/direct/qdldl/amd/src/amd_post_tree.c
new file mode 100644
index 0000000..516c95c
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_post_tree.c
@@ -0,0 +1,120 @@
+/* ========================================================================= */
+/* === AMD_post_tree ======================================================= */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* Post-ordering of a supernodal elimination tree. */
+
+#include "amd_internal.h"
+
+GLOBAL Int AMD_post_tree
+(
+ Int root, /* root of the tree */
+ Int k, /* start numbering at k */
+ Int Child [ ], /* input argument of size nn, undefined on
+ * output. Child [i] is the head of a link
+ * list of all nodes that are children of node
+ * i in the tree. */
+ const Int Sibling [ ], /* input argument of size nn, not modified.
+ * If f is a node in the link list of the
+ * children of node i, then Sibling [f] is the
+ * next child of node i.
+ */
+ Int Order [ ], /* output order, of size nn. Order [i] = k
+ * if node i is the kth node of the reordered
+ * tree. */
+ Int Stack [ ] /* workspace of size nn */
+#ifndef NDEBUG
+ , Int nn /* nodes are in the range 0..nn-1. */
+#endif
+)
+{
+ Int f, head, h, i ;
+
+#if 0
+ /* --------------------------------------------------------------------- */
+ /* recursive version (Stack [ ] is not used): */
+ /* --------------------------------------------------------------------- */
+
+ /* this is simple, but can caouse stack overflow if nn is large */
+ i = root ;
+ for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
+ {
+ k = AMD_post_tree (f, k, Child, Sibling, Order, Stack, nn) ;
+ }
+ Order [i] = k++ ;
+ return (k) ;
+#endif
+
+ /* --------------------------------------------------------------------- */
+ /* non-recursive version, using an explicit stack */
+ /* --------------------------------------------------------------------- */
+
+ /* push root on the stack */
+ head = 0 ;
+ Stack [0] = root ;
+
+ while (head >= 0)
+ {
+ /* get head of stack */
+ ASSERT (head < nn) ;
+ i = Stack [head] ;
+ AMD_DEBUG1 (("head of stack "ID" \n", i)) ;
+ ASSERT (i >= 0 && i < nn) ;
+
+ if (Child [i] != EMPTY)
+ {
+ /* the children of i are not yet ordered */
+ /* push each child onto the stack in reverse order */
+ /* so that small ones at the head of the list get popped first */
+ /* and the biggest one at the end of the list gets popped last */
+ for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
+ {
+ head++ ;
+ ASSERT (head < nn) ;
+ ASSERT (f >= 0 && f < nn) ;
+ }
+ h = head ;
+ ASSERT (head < nn) ;
+ for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
+ {
+ ASSERT (h > 0) ;
+ Stack [h--] = f ;
+ AMD_DEBUG1 (("push "ID" on stack\n", f)) ;
+ ASSERT (f >= 0 && f < nn) ;
+ }
+ ASSERT (Stack [h] == i) ;
+
+ /* delete child list so that i gets ordered next time we see it */
+ Child [i] = EMPTY ;
+ }
+ else
+ {
+ /* the children of i (if there were any) are already ordered */
+ /* remove i from the stack and order it. Front i is kth front */
+ head-- ;
+ AMD_DEBUG1 (("pop "ID" order "ID"\n", i, k)) ;
+ Order [i] = k++ ;
+ ASSERT (k <= nn) ;
+ }
+
+#ifndef NDEBUG
+ AMD_DEBUG1 (("\nStack:")) ;
+ for (h = head ; h >= 0 ; h--)
+ {
+ Int j = Stack [h] ;
+ AMD_DEBUG1 ((" "ID, j)) ;
+ ASSERT (j >= 0 && j < nn) ;
+ }
+ AMD_DEBUG1 (("\n\n")) ;
+ ASSERT (head < nn) ;
+#endif
+
+ }
+ return (k) ;
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_postorder.c b/lin_sys/direct/qdldl/amd/src/amd_postorder.c
new file mode 100644
index 0000000..e5aea7b
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_postorder.c
@@ -0,0 +1,206 @@
+/* ========================================================================= */
+/* === AMD_postorder ======================================================= */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* Perform a postordering (via depth-first search) of an assembly tree. */
+
+#include "amd_internal.h"
+
+GLOBAL void AMD_postorder
+(
+ /* inputs, not modified on output: */
+ Int nn, /* nodes are in the range 0..nn-1 */
+ Int Parent [ ], /* Parent [j] is the parent of j, or EMPTY if root */
+ Int Nv [ ], /* Nv [j] > 0 number of pivots represented by node j,
+ * or zero if j is not a node. */
+ Int Fsize [ ], /* Fsize [j]: size of node j */
+
+ /* output, not defined on input: */
+ Int Order [ ], /* output post-order */
+
+ /* workspaces of size nn: */
+ Int Child [ ],
+ Int Sibling [ ],
+ Int Stack [ ]
+)
+{
+ Int i, j, k, parent, frsize, f, fprev, maxfrsize, bigfprev, bigf, fnext ;
+
+ for (j = 0 ; j < nn ; j++)
+ {
+ Child [j] = EMPTY ;
+ Sibling [j] = EMPTY ;
+ }
+
+ /* --------------------------------------------------------------------- */
+ /* place the children in link lists - bigger elements tend to be last */
+ /* --------------------------------------------------------------------- */
+
+ for (j = nn-1 ; j >= 0 ; j--)
+ {
+ if (Nv [j] > 0)
+ {
+ /* this is an element */
+ parent = Parent [j] ;
+ if (parent != EMPTY)
+ {
+ /* place the element in link list of the children its parent */
+ /* bigger elements will tend to be at the end of the list */
+ Sibling [j] = Child [parent] ;
+ Child [parent] = j ;
+ }
+ }
+ }
+
+#ifndef NDEBUG
+ {
+ Int nels, ff, nchild ;
+ AMD_DEBUG1 (("\n\n================================ AMD_postorder:\n"));
+ nels = 0 ;
+ for (j = 0 ; j < nn ; j++)
+ {
+ if (Nv [j] > 0)
+ {
+ AMD_DEBUG1 (( ""ID" : nels "ID" npiv "ID" size "ID
+ " parent "ID" maxfr "ID"\n", j, nels,
+ Nv [j], Fsize [j], Parent [j], Fsize [j])) ;
+ /* this is an element */
+ /* dump the link list of children */
+ nchild = 0 ;
+ AMD_DEBUG1 ((" Children: ")) ;
+ for (ff = Child [j] ; ff != EMPTY ; ff = Sibling [ff])
+ {
+ AMD_DEBUG1 ((ID" ", ff)) ;
+ ASSERT (Parent [ff] == j) ;
+ nchild++ ;
+ ASSERT (nchild < nn) ;
+ }
+ AMD_DEBUG1 (("\n")) ;
+ parent = Parent [j] ;
+ if (parent != EMPTY)
+ {
+ ASSERT (Nv [parent] > 0) ;
+ }
+ nels++ ;
+ }
+ }
+ }
+ AMD_DEBUG1 (("\n\nGo through the children of each node, and put\n"
+ "the biggest child last in each list:\n")) ;
+#endif
+
+ /* --------------------------------------------------------------------- */
+ /* place the largest child last in the list of children for each node */
+ /* --------------------------------------------------------------------- */
+
+ for (i = 0 ; i < nn ; i++)
+ {
+ if (Nv [i] > 0 && Child [i] != EMPTY)
+ {
+
+#ifndef NDEBUG
+ Int nchild ;
+ AMD_DEBUG1 (("Before partial sort, element "ID"\n", i)) ;
+ nchild = 0 ;
+ for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
+ {
+ ASSERT (f >= 0 && f < nn) ;
+ AMD_DEBUG1 ((" f: "ID" size: "ID"\n", f, Fsize [f])) ;
+ nchild++ ;
+ ASSERT (nchild <= nn) ;
+ }
+#endif
+
+ /* find the biggest element in the child list */
+ fprev = EMPTY ;
+ maxfrsize = EMPTY ;
+ bigfprev = EMPTY ;
+ bigf = EMPTY ;
+ for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
+ {
+ ASSERT (f >= 0 && f < nn) ;
+ frsize = Fsize [f] ;
+ if (frsize >= maxfrsize)
+ {
+ /* this is the biggest seen so far */
+ maxfrsize = frsize ;
+ bigfprev = fprev ;
+ bigf = f ;
+ }
+ fprev = f ;
+ }
+ ASSERT (bigf != EMPTY) ;
+
+ fnext = Sibling [bigf] ;
+
+ AMD_DEBUG1 (("bigf "ID" maxfrsize "ID" bigfprev "ID" fnext "ID
+ " fprev " ID"\n", bigf, maxfrsize, bigfprev, fnext, fprev)) ;
+
+ if (fnext != EMPTY)
+ {
+ /* if fnext is EMPTY then bigf is already at the end of list */
+
+ if (bigfprev == EMPTY)
+ {
+ /* delete bigf from the element of the list */
+ Child [i] = fnext ;
+ }
+ else
+ {
+ /* delete bigf from the middle of the list */
+ Sibling [bigfprev] = fnext ;
+ }
+
+ /* put bigf at the end of the list */
+ Sibling [bigf] = EMPTY ;
+ ASSERT (Child [i] != EMPTY) ;
+ ASSERT (fprev != bigf) ;
+ ASSERT (fprev != EMPTY) ;
+ Sibling [fprev] = bigf ;
+ }
+
+#ifndef NDEBUG
+ AMD_DEBUG1 (("After partial sort, element "ID"\n", i)) ;
+ for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
+ {
+ ASSERT (f >= 0 && f < nn) ;
+ AMD_DEBUG1 ((" "ID" "ID"\n", f, Fsize [f])) ;
+ ASSERT (Nv [f] > 0) ;
+ nchild-- ;
+ }
+ ASSERT (nchild == 0) ;
+#endif
+
+ }
+ }
+
+ /* --------------------------------------------------------------------- */
+ /* postorder the assembly tree */
+ /* --------------------------------------------------------------------- */
+
+ for (i = 0 ; i < nn ; i++)
+ {
+ Order [i] = EMPTY ;
+ }
+
+ k = 0 ;
+
+ for (i = 0 ; i < nn ; i++)
+ {
+ if (Parent [i] == EMPTY && Nv [i] > 0)
+ {
+ AMD_DEBUG1 (("Root of assembly tree "ID"\n", i)) ;
+ k = AMD_post_tree (i, k, Child, Sibling, Order, Stack
+#ifndef NDEBUG
+ , nn
+#endif
+ ) ;
+ }
+ }
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_preprocess.c b/lin_sys/direct/qdldl/amd/src/amd_preprocess.c
new file mode 100644
index 0000000..a8139c3
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_preprocess.c
@@ -0,0 +1,118 @@
+/* ========================================================================= */
+/* === AMD_preprocess ====================================================== */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* Sorts, removes duplicate entries, and transposes from the nonzero pattern of
+ * a column-form matrix A, to obtain the matrix R. The input matrix can have
+ * duplicate entries and/or unsorted columns (AMD_valid (n,Ap,Ai) must not be
+ * AMD_INVALID).
+ *
+ * This input condition is NOT checked. This routine is not user-callable.
+ */
+
+#include "amd_internal.h"
+
+/* ========================================================================= */
+/* === AMD_preprocess ====================================================== */
+/* ========================================================================= */
+
+/* AMD_preprocess does not check its input for errors or allocate workspace.
+ * On input, the condition (AMD_valid (n,n,Ap,Ai) != AMD_INVALID) must hold.
+ */
+
+GLOBAL void AMD_preprocess
+(
+ Int n, /* input matrix: A is n-by-n */
+ const Int Ap [ ], /* size n+1 */
+ const Int Ai [ ], /* size nz = Ap [n] */
+
+ /* output matrix R: */
+ Int Rp [ ], /* size n+1 */
+ Int Ri [ ], /* size nz (or less, if duplicates present) */
+
+ Int W [ ], /* workspace of size n */
+ Int Flag [ ] /* workspace of size n */
+)
+{
+
+ /* --------------------------------------------------------------------- */
+ /* local variables */
+ /* --------------------------------------------------------------------- */
+
+ Int i, j, p, p2 ;
+
+ ASSERT (AMD_valid (n, n, Ap, Ai) != AMD_INVALID) ;
+
+ /* --------------------------------------------------------------------- */
+ /* count the entries in each row of A (excluding duplicates) */
+ /* --------------------------------------------------------------------- */
+
+ for (i = 0 ; i < n ; i++)
+ {
+ W [i] = 0 ; /* # of nonzeros in row i (excl duplicates) */
+ Flag [i] = EMPTY ; /* Flag [i] = j if i appears in column j */
+ }
+ for (j = 0 ; j < n ; j++)
+ {
+ p2 = Ap [j+1] ;
+ for (p = Ap [j] ; p < p2 ; p++)
+ {
+ i = Ai [p] ;
+ if (Flag [i] != j)
+ {
+ /* row index i has not yet appeared in column j */
+ W [i]++ ; /* one more entry in row i */
+ Flag [i] = j ; /* flag row index i as appearing in col j*/
+ }
+ }
+ }
+
+ /* --------------------------------------------------------------------- */
+ /* compute the row pointers for R */
+ /* --------------------------------------------------------------------- */
+
+ Rp [0] = 0 ;
+ for (i = 0 ; i < n ; i++)
+ {
+ Rp [i+1] = Rp [i] + W [i] ;
+ }
+ for (i = 0 ; i < n ; i++)
+ {
+ W [i] = Rp [i] ;
+ Flag [i] = EMPTY ;
+ }
+
+ /* --------------------------------------------------------------------- */
+ /* construct the row form matrix R */
+ /* --------------------------------------------------------------------- */
+
+ /* R = row form of pattern of A */
+ for (j = 0 ; j < n ; j++)
+ {
+ p2 = Ap [j+1] ;
+ for (p = Ap [j] ; p < p2 ; p++)
+ {
+ i = Ai [p] ;
+ if (Flag [i] != j)
+ {
+ /* row index i has not yet appeared in column j */
+ Ri [W [i]++] = j ; /* put col j in row i */
+ Flag [i] = j ; /* flag row index i as appearing in col j*/
+ }
+ }
+ }
+
+#ifndef NDEBUG
+ ASSERT (AMD_valid (n, n, Rp, Ri) == AMD_OK) ;
+ for (j = 0 ; j < n ; j++)
+ {
+ ASSERT (W [j] == Rp [j+1]) ;
+ }
+#endif
+}
diff --git a/lin_sys/direct/qdldl/amd/src/amd_valid.c b/lin_sys/direct/qdldl/amd/src/amd_valid.c
new file mode 100644
index 0000000..609abca
--- /dev/null
+++ b/lin_sys/direct/qdldl/amd/src/amd_valid.c
@@ -0,0 +1,92 @@
+/* ========================================================================= */
+/* === AMD_valid =========================================================== */
+/* ========================================================================= */
+
+/* ------------------------------------------------------------------------- */
+/* AMD, Copyright (c) Timothy A. Davis, */
+/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
+/* email: DrTimothyAldenDavis@gmail.com */
+/* ------------------------------------------------------------------------- */
+
+/* Check if a column-form matrix is valid or not. The matrix A is
+ * n_row-by-n_col. The row indices of entries in column j are in
+ * Ai [Ap [j] ... Ap [j+1]-1]. Required conditions are:
+ *
+ * n_row >= 0
+ * n_col >= 0
+ * nz = Ap [n_col] >= 0 number of entries in the matrix
+ * Ap [0] == 0
+ * Ap [j] <= Ap [j+1] for all j in the range 0 to n_col.
+ * Ai [0 ... nz-1] must be in the range 0 to n_row-1.
+ *
+ * If any of the above conditions hold, AMD_INVALID is returned. If the
+ * following condition holds, AMD_OK_BUT_JUMBLED is returned (a warning,
+ * not an error):
+ *
+ * row indices in Ai [Ap [j] ... Ap [j+1]-1] are not sorted in ascending
+ * order, and/or duplicate entries exist.
+ *
+ * Otherwise, AMD_OK is returned.
+ *
+ * In v1.2 and earlier, this function returned TRUE if the matrix was valid
+ * (now returns AMD_OK), or FALSE otherwise (now returns AMD_INVALID or
+ * AMD_OK_BUT_JUMBLED).
+ */
+
+#include "amd_internal.h"
+
+GLOBAL Int AMD_valid
+(
+ /* inputs, not modified on output: */
+ Int n_row, /* A is n_row-by-n_col */
+ Int n_col,
+ const Int Ap [ ], /* column pointers of A, of size n_col+1 */
+ const Int Ai [ ] /* row indices of A, of size nz = Ap [n_col] */
+)
+{
+ Int nz, j, p1, p2, ilast, i, p, result = AMD_OK ;
+
+ if (n_row < 0 || n_col < 0 || Ap == NULL || Ai == NULL)
+ {
+ return (AMD_INVALID) ;
+ }
+ nz = Ap [n_col] ;
+ if (Ap [0] != 0 || nz < 0)
+ {
+ /* column pointers must start at Ap [0] = 0, and Ap [n] must be >= 0 */
+ AMD_DEBUG0 (("column 0 pointer bad or nz < 0\n")) ;
+ return (AMD_INVALID) ;
+ }
+ for (j = 0 ; j < n_col ; j++)
+ {
+ p1 = Ap [j] ;
+ p2 = Ap [j+1] ;
+ AMD_DEBUG2 (("\nColumn: "ID" p1: "ID" p2: "ID"\n", j, p1, p2)) ;
+ if (p1 > p2)
+ {
+ /* column pointers must be ascending */
+ AMD_DEBUG0 (("column "ID" pointer bad\n", j)) ;
+ return (AMD_INVALID) ;
+ }
+ ilast = EMPTY ;
+ for (p = p1 ; p < p2 ; p++)
+ {
+ i = Ai [p] ;
+ AMD_DEBUG3 (("row: "ID"\n", i)) ;
+ if (i < 0 || i >= n_row)
+ {
+ /* row index out of range */
+ AMD_DEBUG0 (("index out of range, col "ID" row "ID"\n", j, i));
+ return (AMD_INVALID) ;
+ }
+ if (i <= ilast)
+ {
+ /* row index unsorted, or duplicate entry present */
+ AMD_DEBUG1 (("index unsorted/dupl col "ID" row "ID"\n", j, i));
+ result = AMD_OK_BUT_JUMBLED ;
+ }
+ ilast = i ;
+ }
+ }
+ return (result) ;
+}
diff --git a/lin_sys/direct/qdldl/qdldl_interface.c b/lin_sys/direct/qdldl/qdldl_interface.c
new file mode 100644
index 0000000..43f30ef
--- /dev/null
+++ b/lin_sys/direct/qdldl/qdldl_interface.c
@@ -0,0 +1,413 @@
+#include "glob_opts.h"
+
+#include "qdldl.h"
+#include "qdldl_interface.h"
+
+#ifndef EMBEDDED
+#include "amd.h"
+#endif
+
+#if EMBEDDED != 1
+#include "kkt.h"
+#endif
+
+#ifndef EMBEDDED
+
+// Free LDL Factorization structure
+void free_linsys_solver_qdldl(qdldl_solver *s) {
+ if (s) {
+ if (s->L) csc_spfree(s->L);
+ if (s->P) c_free(s->P);
+ if (s->Dinv) c_free(s->Dinv);
+ if (s->bp) c_free(s->bp);
+ if (s->sol) c_free(s->sol);
+ if (s->rho_inv_vec) c_free(s->rho_inv_vec);
+
+ // These are required for matrix updates
+ if (s->Pdiag_idx) c_free(s->Pdiag_idx);
+ if (s->KKT) csc_spfree(s->KKT);
+ if (s->PtoKKT) c_free(s->PtoKKT);
+ if (s->AtoKKT) c_free(s->AtoKKT);
+ if (s->rhotoKKT) c_free(s->rhotoKKT);
+
+ // QDLDL workspace
+ if (s->D) c_free(s->D);
+ if (s->etree) c_free(s->etree);
+ if (s->Lnz) c_free(s->Lnz);
+ if (s->iwork) c_free(s->iwork);
+ if (s->bwork) c_free(s->bwork);
+ if (s->fwork) c_free(s->fwork);
+ c_free(s);
+
+ }
+}
+
+
+/**
+ * Compute LDL factorization of matrix A
+ * @param A Matrix to be factorized
+ * @param p Private workspace
+ * @param nvar Number of QP variables
+ * @return exitstatus (0 is good)
+ */
+static c_int LDL_factor(csc *A, qdldl_solver * p, c_int nvar){
+
+ c_int sum_Lnz;
+ c_int factor_status;
+
+ // Compute elimination tree
+ sum_Lnz = QDLDL_etree(A->n, A->p, A->i, p->iwork, p->Lnz, p->etree);
+
+ if (sum_Lnz < 0){
+ // Error
+#ifdef PRINTING
+ c_eprint("Error in KKT matrix LDL factorization when computing the elimination tree.");
+ if(sum_Lnz == -1){
+ c_eprint("Matrix is not perfectly upper triangular.");
+ }
+ else if(sum_Lnz == -2){
+ c_eprint("Integer overflow in L nonzero count.");
+ }
+#endif
+ return sum_Lnz;
+ }
+
+ // Allocate memory for Li and Lx
+ p->L->i = (c_int *)c_malloc(sizeof(c_int)*sum_Lnz);
+ p->L->x = (c_float *)c_malloc(sizeof(c_float)*sum_Lnz);
+ p->L->nzmax = sum_Lnz;
+
+ // Factor matrix
+ factor_status = QDLDL_factor(A->n, A->p, A->i, A->x,
+ p->L->p, p->L->i, p->L->x,
+ p->D, p->Dinv, p->Lnz,
+ p->etree, p->bwork, p->iwork, p->fwork);
+
+
+ if (factor_status < 0){
+ // Error
+#ifdef PRINTING
+ c_eprint("Error in KKT matrix LDL factorization when computing the nonzero elements. There are zeros in the diagonal matrix");
+#endif
+ return factor_status;
+ } else if (factor_status < nvar) {
+ // Error: Number of positive elements of D should be equal to nvar
+#ifdef PRINTING
+ c_eprint("Error in KKT matrix LDL factorization when computing the nonzero elements. The problem seems to be non-convex");
+#endif
+ return -2;
+ }
+
+ return 0;
+
+}
+
+
+static c_int permute_KKT(csc ** KKT, qdldl_solver * p, c_int Pnz, c_int Anz, c_int m, c_int * PtoKKT, c_int * AtoKKT, c_int * rhotoKKT){
+ c_float *info;
+ c_int amd_status;
+ c_int * Pinv;
+ csc *KKT_temp;
+ c_int * KtoPKPt;
+ c_int i; // Indexing
+
+ info = (c_float *)c_malloc(AMD_INFO * sizeof(c_float));
+
+ // Compute permutation matrix P using AMD
+#ifdef DLONG
+ amd_status = amd_l_order((*KKT)->n, (*KKT)->p, (*KKT)->i, p->P, (c_float *)OSQP_NULL, info);
+#else
+ amd_status = amd_order((*KKT)->n, (*KKT)->p, (*KKT)->i, p->P, (c_float *)OSQP_NULL, info);
+#endif
+ if (amd_status < 0) {
+ // Free Amd info and return an error
+ c_free(info);
+ return amd_status;
+ }
+
+
+ // Inverse of the permutation vector
+ Pinv = csc_pinv(p->P, (*KKT)->n);
+
+ // Permute KKT matrix
+ if (!PtoKKT && !AtoKKT && !rhotoKKT){ // No vectors to be stored
+ // Assign values of mapping
+ KKT_temp = csc_symperm((*KKT), Pinv, OSQP_NULL, 1);
+ }
+ else {
+ // Allocate vector of mappings from unpermuted to permuted
+ KtoPKPt = c_malloc((*KKT)->p[(*KKT)->n] * sizeof(c_int));
+ KKT_temp = csc_symperm((*KKT), Pinv, KtoPKPt, 1);
+
+ // Update vectors PtoKKT, AtoKKT and rhotoKKT
+ if (PtoKKT){
+ for (i = 0; i < Pnz; i++){
+ PtoKKT[i] = KtoPKPt[PtoKKT[i]];
+ }
+ }
+ if (AtoKKT){
+ for (i = 0; i < Anz; i++){
+ AtoKKT[i] = KtoPKPt[AtoKKT[i]];
+ }
+ }
+ if (rhotoKKT){
+ for (i = 0; i < m; i++){
+ rhotoKKT[i] = KtoPKPt[rhotoKKT[i]];
+ }
+ }
+
+ // Cleanup vector of mapping
+ c_free(KtoPKPt);
+ }
+
+ // Cleanup
+ // Free previous KKT matrix and assign pointer to new one
+ csc_spfree((*KKT));
+ (*KKT) = KKT_temp;
+ // Free Pinv
+ c_free(Pinv);
+ // Free Amd info
+ c_free(info);
+
+ return 0;
+}
+
+
+// Initialize LDL Factorization structure
+c_int init_linsys_solver_qdldl(qdldl_solver ** sp, const csc * P, const csc * A, c_float sigma, const c_float * rho_vec, c_int polish){
+
+ // Define Variables
+ csc * KKT_temp; // Temporary KKT pointer
+ c_int i; // Loop counter
+ c_int n_plus_m; // Define n_plus_m dimension
+
+ // Allocate private structure to store KKT factorization
+ qdldl_solver *s;
+ s = c_calloc(1, sizeof(qdldl_solver));
+ *sp = s;
+
+ // Size of KKT
+ s->n = P->n;
+ s->m = A->m;
+ n_plus_m = s->n + s->m;
+
+ // Sigma parameter
+ s->sigma = sigma;
+
+ // Polishing flag
+ s->polish = polish;
+
+ // Link Functions
+ s->solve = &solve_linsys_qdldl;
+
+#ifndef EMBEDDED
+ s->free = &free_linsys_solver_qdldl;
+#endif
+
+#if EMBEDDED != 1
+ s->update_matrices = &update_linsys_solver_matrices_qdldl;
+ s->update_rho_vec = &update_linsys_solver_rho_vec_qdldl;
+#endif
+
+ // Assign type
+ s->type = QDLDL_SOLVER;
+
+ // Set number of threads to 1 (single threaded)
+ s->nthreads = 1;
+
+ // Sparse matrix L (lower triangular)
+ // NB: We don not allocate L completely (CSC elements)
+ // L will be allocated during the factorization depending on the
+ // resulting number of elements.
+ s->L = c_malloc(sizeof(csc));
+ s->L->m = n_plus_m;
+ s->L->n = n_plus_m;
+ s->L->nz = -1;
+
+ // Diagonal matrix stored as a vector D
+ s->Dinv = (QDLDL_float *)c_malloc(sizeof(QDLDL_float) * n_plus_m);
+ s->D = (QDLDL_float *)c_malloc(sizeof(QDLDL_float) * n_plus_m);
+
+ // Permutation vector P
+ s->P = (QDLDL_int *)c_malloc(sizeof(QDLDL_int) * n_plus_m);
+
+ // Working vector
+ s->bp = (QDLDL_float *)c_malloc(sizeof(QDLDL_float) * n_plus_m);
+
+ // Solution vector
+ s->sol = (QDLDL_float *)c_malloc(sizeof(QDLDL_float) * n_plus_m);
+
+ // Parameter vector
+ s->rho_inv_vec = (c_float *)c_malloc(sizeof(c_float) * s->m);
+
+ // Elimination tree workspace
+ s->etree = (QDLDL_int *)c_malloc(n_plus_m * sizeof(QDLDL_int));
+ s->Lnz = (QDLDL_int *)c_malloc(n_plus_m * sizeof(QDLDL_int));
+
+ // Preallocate L matrix (Lx and Li are sparsity dependent)
+ s->L->p = (c_int *)c_malloc((n_plus_m+1) * sizeof(QDLDL_int));
+
+ // Lx and Li are sparsity dependent, so set them to
+ // null initially so we don't try to free them prematurely
+ s->L->i = OSQP_NULL;
+ s->L->x = OSQP_NULL;
+
+ // Preallocate workspace
+ s->iwork = (QDLDL_int *)c_malloc(sizeof(QDLDL_int)*(3*n_plus_m));
+ s->bwork = (QDLDL_bool *)c_malloc(sizeof(QDLDL_bool)*n_plus_m);
+ s->fwork = (QDLDL_float *)c_malloc(sizeof(QDLDL_float)*n_plus_m);
+
+ // Form and permute KKT matrix
+ if (polish){ // Called from polish()
+ // Use s->rho_inv_vec for storing param2 = vec(delta)
+ for (i = 0; i < A->m; i++){
+ s->rho_inv_vec[i] = sigma;
+ }
+
+ KKT_temp = form_KKT(P, A, 0, sigma, s->rho_inv_vec, OSQP_NULL, OSQP_NULL, OSQP_NULL, OSQP_NULL, OSQP_NULL);
+
+ // Permute matrix
+ if (KKT_temp)
+ permute_KKT(&KKT_temp, s, OSQP_NULL, OSQP_NULL, OSQP_NULL, OSQP_NULL, OSQP_NULL, OSQP_NULL);
+ }
+ else { // Called from ADMM algorithm
+
+ // Allocate vectors of indices
+ s->PtoKKT = c_malloc((P->p[P->n]) * sizeof(c_int));
+ s->AtoKKT = c_malloc((A->p[A->n]) * sizeof(c_int));
+ s->rhotoKKT = c_malloc((A->m) * sizeof(c_int));
+
+ // Use p->rho_inv_vec for storing param2 = rho_inv_vec
+ for (i = 0; i < A->m; i++){
+ s->rho_inv_vec[i] = 1. / rho_vec[i];
+ }
+
+ KKT_temp = form_KKT(P, A, 0, sigma, s->rho_inv_vec,
+ s->PtoKKT, s->AtoKKT,
+ &(s->Pdiag_idx), &(s->Pdiag_n), s->rhotoKKT);
+
+ // Permute matrix
+ if (KKT_temp)
+ permute_KKT(&KKT_temp, s, P->p[P->n], A->p[A->n], A->m, s->PtoKKT, s->AtoKKT, s->rhotoKKT);
+ }
+
+ // Check if matrix has been created
+ if (!KKT_temp){
+#ifdef PRINTING
+ c_eprint("Error forming and permuting KKT matrix");
+#endif
+ free_linsys_solver_qdldl(s);
+ *sp = OSQP_NULL;
+ return OSQP_LINSYS_SOLVER_INIT_ERROR;
+ }
+
+ // Factorize the KKT matrix
+ if (LDL_factor(KKT_temp, s, P->n) < 0) {
+ csc_spfree(KKT_temp);
+ free_linsys_solver_qdldl(s);
+ *sp = OSQP_NULL;
+ return OSQP_NONCVX_ERROR;
+ }
+
+ if (polish){ // If KKT passed, assign it to KKT_temp
+ // Polish, no need for KKT_temp
+ csc_spfree(KKT_temp);
+ }
+ else { // If not embedded option 1 copy pointer to KKT_temp. Do not free it.
+ s->KKT = KKT_temp;
+ }
+
+
+ // No error
+ return 0;
+}
+
+#endif // EMBEDDED
+
+
+// Permute x = P*b using P
+void permute_x(c_int n, c_float * x, const c_float * b, const c_int * P) {
+ c_int j;
+ for (j = 0 ; j < n ; j++) x[j] = b[P[j]];
+}
+
+// Permute x = P'*b using P
+void permutet_x(c_int n, c_float * x, const c_float * b, const c_int * P) {
+ c_int j;
+ for (j = 0 ; j < n ; j++) x[P[j]] = b[j];
+}
+
+
+static void LDLSolve(c_float *x, c_float *b, const csc *L, const c_float *Dinv, const c_int *P, c_float *bp) {
+ /* solves P'LDL'P x = b for x */
+ permute_x(L->n, bp, b, P);
+ QDLDL_solve(L->n, L->p, L->i, L->x, Dinv, bp);
+ permutet_x(L->n, x, bp, P);
+
+}
+
+
+c_int solve_linsys_qdldl(qdldl_solver * s, c_float * b) {
+ c_int j;
+
+#ifndef EMBEDDED
+ if (s->polish) {
+ /* stores solution to the KKT system in b */
+ LDLSolve(b, b, s->L, s->Dinv, s->P, s->bp);
+ } else {
+#endif
+ /* stores solution to the KKT system in s->sol */
+ LDLSolve(s->sol, b, s->L, s->Dinv, s->P, s->bp);
+
+ /* copy x_tilde from s->sol */
+ for (j = 0 ; j < s->n ; j++) {
+ b[j] = s->sol[j];
+ }
+
+ /* compute z_tilde from b and s->sol */
+ for (j = 0 ; j < s->m ; j++) {
+ b[j + s->n] += s->rho_inv_vec[j] * s->sol[j + s->n];
+ }
+#ifndef EMBEDDED
+ }
+#endif
+
+ return 0;
+}
+
+
+#if EMBEDDED != 1
+// Update private structure with new P and A
+c_int update_linsys_solver_matrices_qdldl(qdldl_solver * s, const csc *P, const csc *A) {
+
+ // Update KKT matrix with new P
+ update_KKT_P(s->KKT, P, s->PtoKKT, s->sigma, s->Pdiag_idx, s->Pdiag_n);
+
+ // Update KKT matrix with new A
+ update_KKT_A(s->KKT, A, s->AtoKKT);
+
+ return (QDLDL_factor(s->KKT->n, s->KKT->p, s->KKT->i, s->KKT->x,
+ s->L->p, s->L->i, s->L->x, s->D, s->Dinv, s->Lnz,
+ s->etree, s->bwork, s->iwork, s->fwork) < 0);
+
+}
+
+
+c_int update_linsys_solver_rho_vec_qdldl(qdldl_solver * s, const c_float * rho_vec){
+ c_int i;
+
+ // Update internal rho_inv_vec
+ for (i = 0; i < s->m; i++){
+ s->rho_inv_vec[i] = 1. / rho_vec[i];
+ }
+
+ // Update KKT matrix with new rho_vec
+ update_KKT_param2(s->KKT, s->rho_inv_vec, s->rhotoKKT, s->m);
+
+ return (QDLDL_factor(s->KKT->n, s->KKT->p, s->KKT->i, s->KKT->x,
+ s->L->p, s->L->i, s->L->x, s->D, s->Dinv, s->Lnz,
+ s->etree, s->bwork, s->iwork, s->fwork) < 0);
+}
+
+
+#endif
diff --git a/lin_sys/direct/qdldl/qdldl_interface.h b/lin_sys/direct/qdldl/qdldl_interface.h
new file mode 100644
index 0000000..7995ef3
--- /dev/null
+++ b/lin_sys/direct/qdldl/qdldl_interface.h
@@ -0,0 +1,135 @@
+#ifndef QDLDL_INTERFACE_H
+#define QDLDL_INTERFACE_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include "types.h"
+#include "qdldl_types.h"
+
+/**
+ * QDLDL solver structure
+ */
+typedef struct qdldl qdldl_solver;
+
+struct qdldl {
+ enum linsys_solver_type type;
+
+ /**
+ * @name Functions
+ * @{
+ */
+ c_int (*solve)(struct qdldl * self, c_float * b);
+
+#ifndef EMBEDDED
+ void (*free)(struct qdldl * self); ///< Free workspace (only if desktop)
+#endif
+
+ // This used only in non embedded or embedded 2 version
+#if EMBEDDED != 1
+ c_int (*update_matrices)(struct qdldl * self, const csc *P, const csc *A); ///< Update solver matrices
+ c_int (*update_rho_vec)(struct qdldl * self, const c_float * rho_vec); ///< Update rho_vec parameter
+#endif
+
+#ifndef EMBEDDED
+ c_int nthreads;
+#endif
+ /** @} */
+
+ /**
+ * @name Attributes
+ * @{
+ */
+ csc *L; ///< lower triangular matrix in LDL factorization
+ c_float *Dinv; ///< inverse of diag matrix in LDL (as a vector)
+ c_int *P; ///< permutation of KKT matrix for factorization
+ c_float *bp; ///< workspace memory for solves
+ c_float *sol; ///< solution to the KKT system
+ c_float *rho_inv_vec; ///< parameter vector
+ c_float sigma; ///< scalar parameter
+#ifndef EMBEDDED
+ c_int polish; ///< polishing flag
+#endif
+ c_int n; ///< number of QP variables
+ c_int m; ///< number of QP constraints
+
+
+#if EMBEDDED != 1
+ // These are required for matrix updates
+ c_int * Pdiag_idx, Pdiag_n; ///< index and number of diagonal elements in P
+ csc * KKT; ///< Permuted KKT matrix in sparse form (used to update P and A matrices)
+ c_int * PtoKKT, * AtoKKT; ///< Index of elements from P and A to KKT matrix
+ c_int * rhotoKKT; ///< Index of rho places in KKT matrix
+ // QDLDL Numeric workspace
+ QDLDL_float *D;
+ QDLDL_int *etree;
+ QDLDL_int *Lnz;
+ QDLDL_int *iwork;
+ QDLDL_bool *bwork;
+ QDLDL_float *fwork;
+#endif
+
+ /** @} */
+};
+
+
+
+/**
+ * Initialize QDLDL Solver
+ *
+ * @param s Pointer to a private structure
+ * @param P Cost function matrix (upper triangular form)
+ * @param A Constraints matrix
+ * @param sigma Algorithm parameter. If polish, then sigma = delta.
+ * @param rho_vec Algorithm parameter. If polish, then rho_vec = OSQP_NULL.
+ * @param polish Flag whether we are initializing for polish or not
+ * @return Exitflag for error (0 if no errors)
+ */
+c_int init_linsys_solver_qdldl(qdldl_solver ** sp, const csc * P, const csc * A, c_float sigma, const c_float * rho_vec, c_int polish);
+
+/**
+ * Solve linear system and store result in b
+ * @param s Linear system solver structure
+ * @param b Right-hand side
+ * @return Exitflag
+ */
+c_int solve_linsys_qdldl(qdldl_solver * s, c_float * b);
+
+
+#if EMBEDDED != 1
+/**
+ * Update linear system solver matrices
+ * @param s Linear system solver structure
+ * @param P Matrix P
+ * @param A Matrix A
+ * @return Exitflag
+ */
+c_int update_linsys_solver_matrices_qdldl(qdldl_solver * s, const csc *P, const csc *A);
+
+
+
+
+/**
+ * Update rho_vec parameter in linear system solver structure
+ * @param s Linear system solver structure
+ * @param rho_vec new rho_vec value
+ * @return exitflag
+ */
+c_int update_linsys_solver_rho_vec_qdldl(qdldl_solver * s, const c_float * rho_vec);
+
+#endif
+
+#ifndef EMBEDDED
+/**
+ * Free linear system solver
+ * @param s linear system solver object
+ */
+void free_linsys_solver_qdldl(qdldl_solver * s);
+#endif
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/lin_sys/direct/qdldl/qdldl_sources b/lin_sys/direct/qdldl/qdldl_sources
new file mode 160000
index 0000000..7d16b70
--- /dev/null
+++ b/lin_sys/direct/qdldl/qdldl_sources
@@ -0,0 +1 @@
+Subproject commit 7d16b70a10a152682204d745d814b6eb63dc5cd2