[med-svn] [Git][med-team/sofa-framework][master] 4 commits: Do not exclude csparse which is not in Debian; remove SuiteSparse from...
Andreas Tille
gitlab at salsa.debian.org
Thu Sep 19 08:15:22 BST 2019
Andreas Tille pushed to branch master at Debian Med / sofa-framework
Commits:
1256191b by Andreas Tille at 2019-09-19T06:25:39Z
Do not exclude csparse which is not in Debian; remove SuiteSparse from Files-Excludes since its not in the upstream source any more
- - - - -
f24dbaa8 by Andreas Tille at 2019-09-19T06:27:05Z
New upstream version 19.06.01+dfsg2
- - - - -
34a0346a by Andreas Tille at 2019-09-19T06:30:22Z
Updated version 19.06.01+dfsg2 from 'upstream/19.06.01+dfsg2'
with Debian dir 241d2a26c4aa95a7d7aaa2a38e2dfe5ad4731ca0
- - - - -
bbf9b000 by Andreas Tille at 2019-09-19T06:35:52Z
Freshly packaged upstream tarball, use dversionmangle=auto in watch file
- - - - -
11 changed files:
- debian/changelog
- debian/copyright
- debian/watch
- + extlibs/csparse/CMakeLists.txt
- + extlibs/csparse/COPYING.txt
- + extlibs/csparse/CSparseConfig.cmake.in
- + extlibs/csparse/UFconfig.h
- + extlibs/csparse/csparse.c
- + extlibs/csparse/csparse.h
- + extlibs/csparse/ldl.c
- + extlibs/csparse/ldl.h
Changes:
=====================================
debian/changelog
=====================================
@@ -1,4 +1,4 @@
-sofa-framework (19.06.01+dfsg-1) UNRELEASED; urgency=medium
+sofa-framework (19.06.01+dfsg2-1) UNRELEASED; urgency=medium
* New upstream version
* Point watch file to Github
=====================================
debian/copyright
=====================================
@@ -4,7 +4,6 @@ Debianized-By: Dominique Belhachemi <domibel at debian.org>
Debianized-Date: Tue, 24 Jun 2008 20:51:56 +0200
Source: https://github.com/sofa-framework/sofa/releases
Files-Excluded: */libQGLViewer*
- */csparse
*/gtest
*/qwt*
*/eigen-3*
@@ -13,7 +12,6 @@ Files-Excluded: */libQGLViewer*
*/tinyxml
*/*.dll
*/*.lib
- */SuiteSparse
Files: *
Copyright: 2006-2008 MGH, INRIA, USTL, UJF, CNRS
=====================================
debian/watch
=====================================
@@ -1,4 +1,4 @@
version=4
-opts="repacksuffix=+dfsg,dversionmangle=s/\+dfsg//g,repack,compression=xz" \
+opts="repacksuffix=+dfsg,dversionmangle=auto,repack,compression=xz" \
https://github.com/sofa-framework/sofa/releases .*/archive/v at ANY_VERSION@@ARCHIVE_EXT@
=====================================
extlibs/csparse/CMakeLists.txt
=====================================
@@ -0,0 +1,24 @@
+cmake_minimum_required(VERSION 3.1)
+project(csparse VERSION 1.2.0)
+
+set(HEADER_FILES
+ csparse.h
+ ldl.h
+ UFconfig.h)
+
+set(SOURCE_FILES
+ csparse.c
+ ldl.c)
+
+# The code must be relocatable if we want to link a shared library against it.
+if("x${CMAKE_CXX_COMPILER_ID}" STREQUAL "xGNU" OR "x${CMAKE_CXX_COMPILER_ID}" STREQUAL "xClang")
+ add_compile_options(-fPIC)
+endif()
+
+add_library(${PROJECT_NAME} STATIC ${HEADER_FILES} ${SOURCE_FILES})
+target_include_directories(${PROJECT_NAME} PUBLIC "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>")
+target_include_directories(${PROJECT_NAME} PUBLIC "$<INSTALL_INTERFACE:include/extlibs/CSparse>")
+set_target_properties(${PROJECT_NAME} PROPERTIES PUBLIC_HEADER "${HEADER_FILES}")
+
+include(${SOFA_KERNEL_SOURCE_DIR}/SofaFramework/SofaMacros.cmake)
+sofa_generate_package(NAME CSparse VERSION ${PROJECT_VERSION} TARGETS ${PROJECT_NAME} INCLUDE_INSTALL_DIR "extlibs/CSparse")
=====================================
extlibs/csparse/COPYING.txt
=====================================
@@ -0,0 +1,3 @@
+CSPARSE: a Concise Sparse matrix package.
+Copyright (c) 2006, Timothy A. Davis.
+http://www.cise.ufl.edu/research/sparse/CSparse
=====================================
extlibs/csparse/CSparseConfig.cmake.in
=====================================
@@ -0,0 +1,12 @@
+# CMake package configuration file for the csparse library
+
+ at PACKAGE_INIT@
+
+if(NOT TARGET csparse)
+ include("${CMAKE_CURRENT_LIST_DIR}/CSparseTargets.cmake")
+endif()
+
+check_required_components(csparse)
+
+set(CSparse_LIBRARIES csparse)
+set(CSparse_INCLUDE_DIRS @PACKAGE_CSPARSE_INCLUDE_DIR@)
=====================================
extlibs/csparse/UFconfig.h
=====================================
@@ -0,0 +1,118 @@
+/* ========================================================================== */
+/* === UFconfig.h =========================================================== */
+/* ========================================================================== */
+
+/* Configuration file for SuiteSparse: a Suite of Sparse matrix packages
+ * (AMD, COLAMD, CCOLAMD, CAMD, CHOLMOD, UMFPACK, CXSparse, and others).
+ *
+ * UFconfig.h provides the definition of the long integer. On most systems,
+ * a C program can be compiled in LP64 mode, in which long's and pointers are
+ * both 64-bits, and int's are 32-bits. Windows 64, however, uses the LLP64
+ * model, in which int's and long's are 32-bits, and long long's and pointers
+ * are 64-bits.
+ *
+ * SuiteSparse packages that include long integer versions are
+ * intended for the LP64 mode. However, as a workaround for Windows 64
+ * (and perhaps other systems), the long integer can be redefined.
+ *
+ * If _WIN64 is defined, then the __int64 type is used instead of long.
+ *
+ * The long integer can also be defined at compile time. For example, this
+ * could be added to UFconfig.mk:
+ *
+ * CFLAGS = -O -D'UF_long=long long' -D'UF_long_max=9223372036854775801' \
+ * -D'UF_long_id="%lld"'
+ *
+ * This file defines UF_long as either long (on all but _WIN64) or
+ * __int64 on Windows 64. The intent is that a UF_long is always a 64-bit
+ * integer in a 64-bit code. ptrdiff_t might be a better choice than long;
+ * it is always the same size as a pointer.
+ *
+ * This file also defines the SUITESPARSE_VERSION and related definitions.
+ *
+ * Copyright (c) 2007, University of Florida. No licensing restrictions
+ * apply to this file or to the UFconfig directory. Author: Timothy A. Davis.
+ */
+
+#ifndef _UFCONFIG_H
+#define _UFCONFIG_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include <limits.h>
+
+/* ========================================================================== */
+/* === UF_long ============================================================== */
+/* ========================================================================== */
+
+#ifndef UF_long
+
+#ifdef _WIN64
+
+#define UF_long __int64
+#define UF_long_max _I64_MAX
+#define UF_long_id "%I64d"
+
+#else
+
+#define UF_long long
+#define UF_long_max LONG_MAX
+#define UF_long_id "%ld"
+
+#endif
+#endif
+
+/* ========================================================================== */
+/* === SuiteSparse version ================================================== */
+/* ========================================================================== */
+
+/* SuiteSparse is not a package itself, but a collection of packages, some of
+ * which must be used together (UMFPACK requires AMD, CHOLMOD requires AMD,
+ * COLAMD, CAMD, and CCOLAMD, etc). A version number is provided here for the
+ * collection itself. The versions of packages within each version of
+ * SuiteSparse are meant to work together. Combining one packge from one
+ * version of SuiteSparse, with another package from another version of
+ * SuiteSparse, may or may not work.
+ *
+ * SuiteSparse Version 3.2.0 contains the following packages:
+ *
+ * AMD version 2.2.0
+ * CAMD version 2.2.0
+ * COLAMD version 2.7.1
+ * CCOLAMD version 2.7.1
+ * CHOLMOD version 1.7.0
+ * CSparse version 2.2.1
+ * CXSparse version 2.2.1
+ * KLU version 1.0.1
+ * BTF version 1.0.1
+ * LDL version 2.0.1
+ * UFconfig version number is the same as SuiteSparse
+ * UMFPACK version 5.2.0
+ * RBio version 1.1.1
+ * UFcollection version 1.1.1
+ * LINFACTOR version 1.1.0
+ * MESHND version 1.1.0
+ * SSMULT version 1.1.0
+ * MATLAB_Tools no specific version number
+ * SuiteSparseQR version 1.1.0
+ *
+ * Other package dependencies:
+ * BLAS required by CHOLMOD and UMFPACK
+ * LAPACK required by CHOLMOD
+ * METIS 4.0.1 required by CHOLMOD (optional) and KLU (optional)
+ */
+
+#define SUITESPARSE_DATE "Sept 20, 2008"
+#define SUITESPARSE_VER_CODE(main,sub) ((main) * 1000 + (sub))
+#define SUITESPARSE_MAIN_VERSION 3
+#define SUITESPARSE_SUB_VERSION 2
+#define SUITESPARSE_SUBSUB_VERSION 0
+#define SUITESPARSE_VERSION \
+ SUITESPARSE_VER_CODE(SUITESPARSE_MAIN_VERSION,SUITESPARSE_SUB_VERSION)
+
+#ifdef __cplusplus
+}
+#endif
+#endif
=====================================
extlibs/csparse/csparse.c
=====================================
@@ -0,0 +1,2122 @@
+# include <stdlib.h>
+# include <limits.h>
+# include <math.h>
+# include <stdio.h>
+
+# include "csparse.h"
+
+cs *cs_add ( const cs *A, const cs *B, double alpha, double beta )
+/*
+ Purpose:
+
+ CS_ADD computes C = alpha*A + beta*B for sparse A and B.
+
+ Reference:
+
+ Timothy Davis,
+ Direct Methods for Sparse Linear Systems,
+ SIAM, Philadelphia, 2006.
+*/
+{
+ int p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values ;
+ double *x, *Bx, *Cx ;
+ cs *C ;
+ if (!A || !B) return (NULL) ; /* check inputs */
+ m = A->m ; anz = A->p [A->n] ;
+ n = B->n ; Bp = B->p ; Bx = B->x ; bnz = Bp [n] ;
+ w = cs_calloc (m, sizeof (int)) ;
+ values = (A->x != NULL) && (Bx != NULL) ;
+ x = values ? cs_malloc (m, sizeof (double)) : NULL ;
+ C = cs_spalloc (m, n, anz + bnz, values, 0) ;
+ if (!C || !w || (values && !x)) return (cs_done (C, w, x, 0)) ;
+ Cp = C->p ; Ci = C->i ; Cx = C->x ;
+ for (j = 0 ; j < n ; j++)
+ {
+ Cp [j] = nz ; /* column j of C starts here */
+ nz = cs_scatter (A, j, alpha, w, x, j+1, C, nz) ; /* alpha*A(:,j)*/
+ nz = cs_scatter (B, j, beta, w, x, j+1, C, nz) ; /* beta*B(:,j) */
+ if (values) for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Ci [p]] ;
+ }
+ Cp [n] = nz ; /* finalize the last column of C */
+ cs_sprealloc (C, 0) ; /* remove extra space from C */
+ return (cs_done (C, w, x, 1)) ; /* success; free workspace, return C */
+}
+static int cs_wclear (int mark, int lemax, int *w, int n)
+/*
+ Purpose:
+
+ CS_WCLEAR clears W.
+
+ Reference:
+
+ Timothy Davis,
+ Direct Methods for Sparse Linear Systems,
+ SIAM, Philadelphia, 2006.
+*/
+{
+ int k ;
+ if (mark < 2 || (mark + lemax < 0))
+ {
+ for (k = 0 ; k < n ; k++) if (w [k] != 0) w [k] = 1 ;
+ mark = 2 ;
+ }
+ return (mark) ; /* at this point, w [0..n-1] < mark holds */
+}
+
+/* keep off-diagonal entries; drop diagonal entries */
+static int cs_diag (int i, int j, double aij, void * other)
+{
+ return (i != j);
+ (void)aij; (void)other; // unused parameters
+}
+
+/* p = amd(A+A') if symmetric is true, or amd(A'A) otherwise */
+int *cs_amd ( const cs *A, int order )
+/*
+ Purpose:
+
+ CS_AMD carries out the approximate minimum degree algorithm.
+
+ Reference:
+
+ Timothy Davis,
+ Direct Methods for Sparse Linear Systems,
+ SIAM, Philadelphia, 2006.
+
+ Parameters:
+
+ Input, int ORDER:
+ -1:natural,
+ 0:Cholesky,
+ 1:LU,
+ 2:QR
+*/
+{
+ cs *C, *A2, *AT ;
+ int *Cp, *Ci, *last, *ww, *len, *nv, *next, *P, *head, *elen, *degree, *w,
+ *hhead, *ATp, *ATi, d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
+ k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
+ ok, cnz, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, n, m ;
+ unsigned int h ;
+ /* --- Construct matrix C ----------------------------------------------- */
+ if (!A || order < 0) return (NULL) ; /* check inputs; quick return */
+ AT = cs_transpose (A, 0) ; /* compute A' */
+ if (!AT) return (NULL) ;
+ m = A->m ; n = A->n ;
+ dense = (int)CS_MAX (16, 10 * sqrt ((double) n)) ; /* find dense threshold */
+ dense = CS_MIN (n-2, dense) ;
+ if (order == 0 && n == m)
+ {
+ C = cs_add (A, AT, 0, 0) ; /* C = A+A' */
+ }
+ else if (order == 1)
+ {
+ ATp = AT->p ; /* drop dense columns from AT */
+ ATi = AT->i ;
+ for (p2 = 0, j = 0 ; j < m ; j++)
+ {
+ p = ATp [j] ; /* column j of AT starts here */
+ ATp [j] = p2 ; /* new column j starts here */
+ if (ATp [j+1] - p > dense) continue ; /* skip dense col j */
+ for ( ; p < ATp [j+1] ; p++) ATi [p2++] = ATi [p] ;
+ }
+ ATp [m] = p2 ; /* finalize AT */
+ A2 = cs_transpose (AT, 0) ; /* A2 = AT' */
+ C = A2 ? cs_multiply (AT, A2) : NULL ; /* C=A'*A with no dense rows */
+ cs_spfree (A2) ;
+ }
+ else
+ {
+ C = cs_multiply (AT, A) ; /* C=A'*A */
+ }
+ cs_spfree (AT) ;
+ if (!C) return (NULL) ;
+ P = cs_malloc (n+1, sizeof (int)) ; /* allocate result */
+ ww = cs_malloc (8*(n+1), sizeof (int)) ;/* get workspace */
+ len = ww ; nv = ww + (n+1) ; next = ww + 2*(n+1) ;
+ head = ww + 3*(n+1) ; elen = ww + 4*(n+1) ; degree = ww + 5*(n+1) ;
+ w = ww + 6*(n+1) ; hhead = ww + 7*(n+1) ;
+ last = P ; /* use P as workspace for last */
+ cs_fkeep (C, &cs_diag, NULL) ; /* drop diagonal entries */
+ Cp = C->p ;
+ cnz = Cp [n] ;
+ if (!cs_sprealloc (C, cnz+cnz/5+2*n)) return (cs_idone (P, C, ww, 0)) ;
+ /* --- Initialize quotient graph ---------------------------------------- */
+ for (k = 0 ; k < n ; k++) len [k] = Cp [k+1] - Cp [k] ;
+ len [n] = 0 ;
+ nzmax = C->nzmax ;
+ Ci = C->i ;
+ for (i = 0 ; i <= n ; i++)
+ {
+ head [i] = -1 ; /* degree list i is empty */
+ last [i] = -1 ;
+ next [i] = -1 ;
+ hhead [i] = -1 ; /* hash list i is empty */
+ nv [i] = 1 ; /* node i is just one node */
+ w [i] = 1 ; /* node i is alive */
+ elen [i] = 0 ; /* Ek of node i is empty */
+ degree [i] = len [i] ; /* degree of node i */
+ }
+ mark = cs_wclear (0, 0, w, n) ; /* clear w */
+ elen [n] = -2 ; /* n is a dead element */
+ Cp [n] = -1 ; /* n is a root of assembly tree */
+ w [n] = 0 ; /* n is a dead element */
+ /* --- Initialize degree lists ------------------------------------------ */
+ for (i = 0 ; i < n ; i++)
+ {
+ d = degree [i] ;
+ if (d == 0) /* node i is empty */
+ {
+ elen [i] = -2 ; /* element i is dead */
+ nel++ ;
+ Cp [i] = -1 ; /* i is a root of assemby tree */
+ w [i] = 0 ;
+ }
+ else if (d > dense) /* node i is dense */
+ {
+ nv [i] = 0 ; /* absorb i into element n */
+ elen [i] = -1 ; /* node i is dead */
+ nel++ ;
+ Cp [i] = CS_FLIP (n) ;
+ nv [n]++ ;
+ }
+ else
+ {
+ if (head [d] != -1) last [head [d]] = i ;
+ next [i] = head [d] ; /* put node i in degree list d */
+ head [d] = i ;
+ }
+ }
+ while (nel < n) /* while (selecting pivots) do */
+ {
+ /* --- Select node of minimum approximate degree -------------------- */
+ for (k = -1 ; mindeg < n && (k = head [mindeg]) == -1 ; mindeg++) ;
+ if (next [k] != -1) last [next [k]] = -1 ;
+ head [mindeg] = next [k] ; /* remove k from degree list */
+ elenk = elen [k] ; /* elenk = |Ek| */
+ nvk = nv [k] ; /* # of nodes k represents */
+ nel += nvk ; /* nv[k] nodes of A eliminated */
+ /* --- Garbage collection ------------------------------------------- */
+ if (elenk > 0 && cnz + mindeg >= nzmax)
+ {
+ for (j = 0 ; j < n ; j++)
+ {
+ if ((p = Cp [j]) >= 0) /* j is a live node or element */
+ {
+ Cp [j] = Ci [p] ; /* save first entry of object */
+ Ci [p] = CS_FLIP (j) ; /* first entry is now CS_FLIP(j) */
+ }
+ }
+ for (q = 0, p = 0 ; p < cnz ; ) /* scan all of memory */
+ {
+ if ((j = CS_FLIP (Ci [p++])) >= 0) /* found object j */
+ {
+ Ci [q] = Cp [j] ; /* restore first entry of object */
+ Cp [j] = q++ ; /* new pointer to object j */
+ for (k3 = 0 ; k3 < len [j]-1 ; k3++) Ci [q++] = Ci [p++] ;
+ }
+ }
+ cnz = q ; /* Ci [cnz...nzmax-1] now free */
+ }
+ /* --- Construct new element ---------------------------------------- */
+ dk = 0 ;
+ nv [k] = -nvk ; /* flag k as in Lk */
+ p = Cp [k] ;
+ pk1 = (elenk == 0) ? p : cnz ; /* do in place if elen[k] == 0 */
+ pk2 = pk1 ;
+ for (k1 = 1 ; k1 <= elenk + 1 ; k1++)
+ {
+ if (k1 > elenk)
+ {
+ e = k ; /* search the nodes in k */
+ pj = p ; /* list of nodes starts at Ci[pj]*/
+ ln = len [k] - elenk ; /* length of list of nodes in k */
+ }
+ else
+ {
+ e = Ci [p++] ; /* search the nodes in e */
+ pj = Cp [e] ;
+ ln = len [e] ; /* length of list of nodes in e */
+ }
+ for (k2 = 1 ; k2 <= ln ; k2++)
+ {
+ i = Ci [pj++] ;
+ if ((nvi = nv [i]) <= 0) continue ; /* node i dead, or seen */
+ dk += nvi ; /* degree[Lk] += size of node i */
+ nv [i] = -nvi ; /* negate nv[i] to denote i in Lk*/
+ Ci [pk2++] = i ; /* place i in Lk */
+ if (next [i] != -1) last [next [i]] = last [i] ;
+ if (last [i] != -1) /* remove i from degree list */
+ {
+ next [last [i]] = next [i] ;
+ }
+ else
+ {
+ head [degree [i]] = next [i] ;
+ }
+ }
+ if (e != k)
+ {
+ Cp [e] = CS_FLIP (k) ; /* absorb e into k */
+ w [e] = 0 ; /* e is now a dead element */
+ }
+ }
+ if (elenk != 0) cnz = pk2 ; /* Ci [cnz...nzmax] is free */
+ degree [k] = dk ; /* external degree of k - |Lk\i| */
+ Cp [k] = pk1 ; /* element k is in Ci[pk1..pk2-1] */
+ len [k] = pk2 - pk1 ;
+ elen [k] = -2 ; /* k is now an element */
+ /* --- Find set differences ----------------------------------------- */
+ mark = cs_wclear (mark, lemax, w, n) ; /* clear w if necessary */
+ for (pk = pk1 ; pk < pk2 ; pk++) /* scan 1: find |Le\Lk| */
+ {
+ i = Ci [pk] ;
+ if ((eln = elen [i]) <= 0) continue ;/* skip if elen[i] empty */
+ nvi = -nv [i] ; /* nv [i] was negated */
+ wnvi = mark - nvi ;
+ for (p = Cp [i] ; p <= Cp [i] + eln - 1 ; p++) /* scan Ei */
+ {
+ e = Ci [p] ;
+ if (w [e] >= mark)
+ {
+ w [e] -= nvi ; /* decrement |Le\Lk| */
+ }
+ else if (w [e] != 0) /* ensure e is a live element */
+ {
+ w [e] = degree [e] + wnvi ; /* 1st time e seen in scan 1 */
+ }
+ }
+ }
+ /* --- Degree update ------------------------------------------------ */
+ for (pk = pk1 ; pk < pk2 ; pk++) /* scan2: degree update */
+ {
+ i = Ci [pk] ; /* consider node i in Lk */
+ p1 = Cp [i] ;
+ p2 = p1 + elen [i] - 1 ;
+ pn = p1 ;
+ for (h = 0, d = 0, p = p1 ; p <= p2 ; p++) /* scan Ei */
+ {
+ e = Ci [p] ;
+ if (w [e] != 0) /* e is an unabsorbed element */
+ {
+ dext = w [e] - mark ; /* dext = |Le\Lk| */
+ if (dext > 0)
+ {
+ d += dext ; /* sum up the set differences */
+ Ci [pn++] = e ; /* keep e in Ei */
+ h += e ; /* compute the hash of node i */
+ }
+ else
+ {
+ Cp [e] = CS_FLIP (k) ; /* aggressive absorb. e->k */
+ w [e] = 0 ; /* e is a dead element */
+ }
+ }
+ }
+ elen [i] = pn - p1 + 1 ; /* elen[i] = |Ei| */
+ p3 = pn ;
+ p4 = p1 + len [i] ;
+ for (p = p2 + 1 ; p < p4 ; p++) /* prune edges in Ai */
+ {
+ j = Ci [p] ;
+ if ((nvj = nv [j]) <= 0) continue ; /* node j dead or in Lk */
+ d += nvj ; /* degree(i) += |j| */
+ Ci [pn++] = j ; /* place j in node list of i */
+ h += j ; /* compute hash for node i */
+ }
+ if (d == 0) /* check for mass elimination */
+ {
+ Cp [i] = CS_FLIP (k) ; /* absorb i into k */
+ nvi = -nv [i] ;
+ dk -= nvi ; /* |Lk| -= |i| */
+ nvk += nvi ; /* |k| += nv[i] */
+ nel += nvi ;
+ nv [i] = 0 ;
+ elen [i] = -1 ; /* node i is dead */
+ }
+ else
+ {
+ degree [i] = CS_MIN (degree [i], d) ; /* update degree(i) */
+ Ci [pn] = Ci [p3] ; /* move first node to end */
+ Ci [p3] = Ci [p1] ; /* move 1st el. to end of Ei */
+ Ci [p1] = k ; /* add k as 1st element in of Ei */
+ len [i] = pn - p1 + 1 ; /* new len of adj. list of node i */
+ h %= n ; /* finalize hash of i */
+ next [i] = hhead [h] ; /* place i in hash bucket */
+ hhead [h] = i ;
+ last [i] = h ; /* save hash of i in last[i] */
+ }
+ } /* scan2 is done */
+ degree [k] = dk ; /* finalize |Lk| */
+ lemax = CS_MAX (lemax, dk) ;
+ mark = cs_wclear (mark+lemax, lemax, w, n) ; /* clear w */
+ /* --- Supernode detection ------------------------------------------ */
+ for (pk = pk1 ; pk < pk2 ; pk++)
+ {
+ i = Ci [pk] ;
+ if (nv [i] >= 0) continue ; /* skip if i is dead */
+ h = last [i] ; /* scan hash bucket of node i */
+ i = hhead [h] ;
+ hhead [h] = -1 ; /* hash bucket will be empty */
+ for ( ; i != -1 && next [i] != -1 ; i = next [i], mark++)
+ {
+ ln = len [i] ;
+ eln = elen [i] ;
+ for (p = Cp[i]+1 ; p <= Cp[i]+ln-1 ; p++) w [Ci [p]] = mark ;
+ jlast = i ;
+ for (j = next [i] ; j != -1 ; ) /* compare i with all j */
+ {
+ ok = (len [j] == ln) && (elen [j] == eln) ;
+ for (p = Cp [j] + 1 ; ok && p <= Cp [j] + ln - 1 ; p++)
+ {
+ if (w [Ci [p]] != mark) ok = 0 ; /* compare i and j*/
+ }
+ if (ok) /* i and j are identical */
+ {
+ Cp [j] = CS_FLIP (i) ; /* absorb j into i */
+ nv [i] += nv [j] ;
+ nv [j] = 0 ;
+ elen [j] = -1 ; /* node j is dead */
+ j = next [j] ; /* delete j from hash bucket */
+ next [jlast] = j ;
+ }
+ else
+ {
+ jlast = j ; /* j and i are different */
+ j = next [j] ;
+ }
+ }
+ }
+ }
+ /* --- Finalize new element------------------------------------------ */
+ for (p = pk1, pk = pk1 ; pk < pk2 ; pk++) /* finalize Lk */
+ {
+ i = Ci [pk] ;
+ if ((nvi = -nv [i]) <= 0) continue ;/* skip if i is dead */
+ nv [i] = nvi ; /* restore nv[i] */
+ d = degree [i] + dk - nvi ; /* compute external degree(i) */
+ d = CS_MIN (d, n - nel - nvi) ;
+ if (head [d] != -1) last [head [d]] = i ;
+ next [i] = head [d] ; /* put i back in degree list */
+ last [i] = -1 ;
+ head [d] = i ;
+ mindeg = CS_MIN (mindeg, d) ; /* find new minimum degree */
+ degree [i] = d ;
+ Ci [p++] = i ; /* place i in Lk */
+ }
+ nv [k] = nvk ; /* # nodes absorbed into k */
+ if ((len [k] = p-pk1) == 0) /* length of adj list of element k*/
+ {
+ Cp [k] = -1 ; /* k is a root of the tree */
+ w [k] = 0 ; /* k is now a dead element */
+ }
+ if (elenk != 0) cnz = p ; /* free unused space in Lk */
+ }
+ /* --- Postordering ----------------------------------------------------- */
+ for (i = 0 ; i < n ; i++) Cp [i] = CS_FLIP (Cp [i]) ;/* fix assembly tree */
+ for (j = 0 ; j <= n ; j++) head [j] = -1 ;
+ for (j = n ; j >= 0 ; j--) /* place unordered nodes in lists */
+ {
+ if (nv [j] > 0) continue ; /* skip if j is an element */
+ next [j] = head [Cp [j]] ; /* place j in list of its parent */
+ head [Cp [j]] = j ;
+ }
+ for (e = n ; e >= 0 ; e--) /* place elements in lists */
+ {
+ if (nv [e] <= 0) continue ; /* skip unless e is an element */
+ if (Cp [e] != -1)
+ {
+ next [e] = head [Cp [e]] ; /* place e in list of its parent */
+ head [Cp [e]] = e ;
+ }
+ }
+ for (k = 0, i = 0 ; i <= n ; i++) /* postorder the assembly tree */
+ {
+ if (Cp [i] == -1) k = cs_tdfs (i, k, head, next, P, w) ;
+ }
+ return (cs_idone (P, C, ww, 1)) ;
+}
+
+/* compute nonzero pattern of L(k,:) */
+static
+int cs_ereach (const cs *A, int k, const int *parent, int *s, int *w,
+ double *x, int top)
+{
+ int i, p, len, *Ap = A->p, *Ai = A->i ;
+ double *Ax = A->x ;
+ for (p = Ap [k] ; p < Ap [k+1] ; p++) /* get pattern of L(k,:) */
+ {
+ i = Ai [p] ; /* A(i,k) is nonzero */
+ if (i > k) continue ; /* only use upper triangular part of A */
+ x [i] = Ax [p] ; /* x(i) = A(i,k) */
+ for (len = 0 ; w [i] != k ; i = parent [i]) /* traverse up etree */
+ {
+ s [len++] = i ; /* L(k,i) is nonzero */
+ w [i] = k ; /* mark i as visited */
+ }
+ while (len > 0) s [--top] = s [--len] ; /* push path onto stack */
+ }
+ return (top) ; /* s [top..n-1] contains pattern of L(k,:)*/
+}
+
+/* L = chol (A, [Pinv parent cp]), Pinv is optional */
+csn *cs_chol (const cs *A, const css *S)
+{
+ double d, lki, *Lx, *x ;
+ int top, i, p, k, n, *Li, *Lp, *cp, *Pinv, *w, *s, *c, *parent ;
+ cs *L, *C, *E ;
+ csn *N ;
+ if (!A || !S || !S->cp || !S->parent) return (NULL) ; /* check inputs */
+ n = A->n ;
+ N = cs_calloc (1, sizeof (csn)) ;
+ w = cs_malloc (3*n, sizeof (int)) ; s = w + n, c = w + 2*n ;
+ x = cs_malloc (n, sizeof (double)) ;
+ cp = S->cp ; Pinv = S->Pinv ; parent = S->parent ;
+ C = Pinv ? cs_symperm (A, Pinv, 1) : ((cs *) A) ;
+ E = Pinv ? C : NULL ;
+ if (!N || !w || !x || !C) return (cs_ndone (N, E, w, x, 0)) ;
+ N->L = L = cs_spalloc (n, n, cp [n], 1, 0) ;
+ if (!L) return (cs_ndone (N, E, w, x, 0)) ;
+ Lp = L->p ; Li = L->i ; Lx = L->x ;
+ for (k = 0 ; k < n ; k++)
+ {
+ /* --- Nonzero pattern of L(k,:) ------------------------------------ */
+ Lp [k] = c [k] = cp [k] ; /* column k of L starts here */
+ x [k] = 0 ; /* x (0:k) is now zero */
+ w [k] = k ; /* mark node k as visited */
+ top = cs_ereach (C, k, parent, s, w, x, n) ; /* find row k of L*/
+ d = x [k] ; /* d = C(k,k) */
+ x [k] = 0 ; /* clear workspace for k+1st iteration */
+ /* --- Triangular solve --------------------------------------------- */
+ for ( ; top < n ; top++) /* solve L(0:k-1,0:k-1) * x = C(:,k) */
+ {
+ i = s [top] ; /* s [top..n-1] is pattern of L(k,:) */
+ lki = x [i] / Lx [Lp [i]] ; /* L(k,i) = x (i) / L(i,i) */
+ x [i] = 0 ; /* clear workspace for k+1st iteration */
+ for (p = Lp [i] + 1 ; p < c [i] ; p++)
+ {
+ x [Li [p]] -= Lx [p] * lki ;
+ }
+ d -= lki * lki ; /* d = d - L(k,i)*L(k,i) */
+ p = c [i]++ ;
+ Li [p] = k ; /* store L(k,i) in column i */
+ Lx [p] = lki ;
+ }
+ /* --- Compute L(k,k) ----------------------------------------------- */
+ if (d <= 0) return (cs_ndone (N, E, w, x, 0)) ; /* not pos def */
+ p = c [k]++ ;
+ Li [p] = k ; /* store L(k,k) = sqrt (d) in column k */
+ Lx [p] = sqrt (d) ;
+ }
+ Lp [n] = cp [n] ; /* finalize L */
+ return (cs_ndone (N, E, w, x, 1)) ; /* success: free E,w,x; return N */
+}
+
+
+/* x=A\b where A is symmetric positive definite; b overwritten with solution */
+int cs_cholsol (const cs *A, double *b, int order)
+{
+ double *x ;
+ css *S ;
+ csn *N ;
+ int n, ok ;
+ if (!A || !b) return (0) ; /* check inputs */
+ n = A->n ;
+ S = cs_schol (A, order) ; /* ordering and symbolic analysis */
+ N = cs_chol (A, S) ; /* numeric Cholesky factorization */
+ x = cs_malloc (n, sizeof (double)) ;
+ ok = (S && N && x) ;
+ if (ok)
+ {
+ cs_ipvec (n, S->Pinv, b, x) ; /* x = P*b */
+ cs_lsolve (N->L, x) ; /* x = L\x */
+ cs_ltsolve (N->L, x) ; /* x = L'\x */
+ cs_pvec (n, S->Pinv, x, b) ; /* b = P'*x */
+ }
+ cs_free (x) ;
+ cs_sfree (S) ;
+ cs_nfree (N) ;
+ return (ok) ;
+}
+
+/* process edge (j,i) of the matrix */
+static void cs_cedge (int j, int i, const int *first, int *maxfirst, int *delta,
+ int *prevleaf, int *ancestor)
+{
+ int q, s, sparent, jprev ;
+ if (i <= j || first [j] <= maxfirst [i]) return ;
+ maxfirst [i] = first [j] ; /* update max first[j] seen so far */
+ jprev = prevleaf [i] ; /* j is a leaf of the ith subtree */
+ delta [j]++ ; /* A(i,j) is in the skeleton matrix */
+ if (jprev != -1)
+ {
+ /* q = least common ancestor of jprev and j */
+ for (q = jprev ; q != ancestor [q] ; q = ancestor [q]) ;
+ for (s = jprev ; s != q ; s = sparent)
+ {
+ sparent = ancestor [s] ; /* path compression */
+ ancestor [s] = q ;
+ }
+ delta [q]-- ; /* decrement to account for overlap in q */
+ }
+ prevleaf [i] = j ; /* j is now previous leaf of ith subtree */
+}
+
+/* colcount = column counts of LL'=A or LL'=A'A, given parent & post ordering */
+int *cs_counts (const cs *A, const int *parent, const int *post, int ata)
+{
+ int i, j, k, p, n, m, ii, s, *ATp, *ATi, *maxfirst, *prevleaf, *ancestor,
+ *head = NULL, *next = NULL, *colcount, *w, *first, *delta ;
+ cs *AT ;
+ if (!A || !parent || !post) return (NULL) ; /* check inputs */
+ m = A->m ; n = A->n ;
+ s = 4*n + (ata ? (n+m+1) : 0) ;
+ w = cs_malloc (s, sizeof (int)) ; first = w+3*n ; /* get workspace */
+ ancestor = w ; maxfirst = w+n ; prevleaf = w+2*n ;
+ delta = colcount = cs_malloc (n, sizeof (int)) ; /* allocate result */
+ AT = cs_transpose (A, 0) ;
+ if (!AT || !colcount || !w) return (cs_idone (colcount, AT, w, 1)) ;
+ for (k = 0 ; k < s ; k++) w [k] = -1 ; /* clear workspace w [0..s-1] */
+ for (k = 0 ; k < n ; k++) /* find first [j] */
+ {
+ j = post [k] ;
+ delta [j] = (first [j] == -1) ? 1 : 0 ; /* delta[j]=1 if j is a leaf */
+ for ( ; j != -1 && first [j] == -1 ; j = parent [j]) first [j] = k ;
+ }
+ ATp = AT->p ; ATi = AT->i ;
+ if (ata)
+ {
+ head = w+4*n ; next = w+5*n+1 ;
+ for (k = 0 ; k < n ; k++) w [post [k]] = k ; /* invert post */
+ for (i = 0 ; i < m ; i++)
+ {
+ k = n ; /* k = least postordered column in row i */
+ for (p = ATp [i] ; p < ATp [i+1] ; p++) k = CS_MIN (k, w [ATi [p]]);
+ next [i] = head [k] ; /* place row i in link list k */
+ head [k] = i ;
+ }
+ }
+ for (i = 0 ; i < n ; i++) ancestor [i] = i ; /* each node in its own set */
+ for (k = 0 ; k < n ; k++)
+ {
+ j = post [k] ; /* j is the kth node in postordered etree */
+ if (parent [j] != -1) delta [parent [j]]-- ; /* j is not a root */
+ if (ata)
+ {
+ for (ii = head [k] ; ii != -1 ; ii = next [ii])
+ {
+ for (p = ATp [ii] ; p < ATp [ii+1] ; p++)
+ cs_cedge (j, ATi [p], first, maxfirst, delta, prevleaf,
+ ancestor) ;
+ }
+ }
+ else
+ {
+ for (p = ATp [j] ; p < ATp [j+1] ; p++)
+ cs_cedge (j, ATi [p], first, maxfirst, delta, prevleaf,
+ ancestor) ;
+ }
+ if (parent [j] != -1) ancestor [j] = parent [j] ;
+ }
+ for (j = 0 ; j < n ; j++) /* sum up delta's of each child */
+ {
+ if (parent [j] != -1) colcount [parent [j]] += colcount [j] ;
+ }
+ return (cs_idone (colcount, AT, w, 1)) ; /* success: free workspace */
+}
+
+/* p [0..n] = cumulative sum of c [0..n-1], and then copy p [0..n-1] into c */
+int cs_cumsum (int *p, int *c, int n)
+{
+ int i, nz = 0 ;
+ if (!p || !c) return (-1) ; /* check inputs */
+ for (i = 0 ; i < n ; i++)
+ {
+ p [i] = nz ;
+ nz += c [i] ;
+ c [i] = p [i] ;
+ }
+ p [n] = nz ;
+ return (nz) ; /* return sum (c [0..n-1]) */
+}
+
+/* depth-first-search of the graph of a matrix, starting at node j */
+int cs_dfs (int j, cs *L, int top, int *xi, int *pstack, const int *Pinv)
+{
+ int i, p, p2, done, jnew, head = 0, *Lp, *Li ;
+ if (!L || !xi || !pstack) return (-1) ;
+ Lp = L->p ; Li = L->i ;
+ xi [0] = j ; /* initialize the recursion stack */
+ while (head >= 0)
+ {
+ j = xi [head] ; /* get j from the top of the recursion stack */
+ jnew = Pinv ? (Pinv [j]) : j ;
+ if (!CS_MARKED(Lp,j))
+ {
+ CS_MARK (Lp,j) ; /* mark node j as visited */
+ pstack [head] = (jnew < 0) ? 0 : CS_UNFLIP (Lp [jnew]) ;
+ }
+ done = 1 ; /* node j done if no unvisited neighbors */
+ p2 = (jnew < 0) ? 0 : CS_UNFLIP (Lp [jnew+1]) ;
+ for (p = pstack [head] ; p < p2 ; p++) /* examine all neighbors of j */
+ {
+ i = Li [p] ; /* consider neighbor node i */
+ if (CS_MARKED (Lp,i)) continue ; /* skip visited node i */
+ pstack [head] = p ; /* pause depth-first search of node j */
+ xi [++head] = i ; /* start dfs at node i */
+ done = 0 ; /* node j is not done */
+ break ; /* break, to start dfs (i) */
+ }
+ if (done) /* depth-first search at node j is done */
+ {
+ head-- ; /* remove j from the recursion stack */
+ xi [--top] = j ; /* and place in the output stack */
+ }
+ }
+ return (top) ;
+}
+
+/* breadth-first search for coarse decomposition (C0,C1,R1 or R0,R3,C3) */
+static int cs_bfs (const cs *A, int n, int *wi, int *wj, int *queue,
+ const int *imatch, const int *jmatch, int mark)
+{
+ int *Ap, *Ai, head = 0, tail = 0, j, i, p, j2 ;
+ cs *C ;
+ for (j = 0 ; j < n ; j++) /* place all unmatched nodes in queue */
+ {
+ if (imatch [j] >= 0) continue ; /* skip j if matched */
+ wj [j] = 0 ; /* j in set C0 (R0 if transpose) */
+ queue [tail++] = j ; /* place unmatched col j in queue */
+ }
+ if (tail == 0) return (1) ; /* quick return if no unmatched nodes */
+ C = (mark == 1) ? ((cs *) A) : cs_transpose (A, 0) ;
+ if (!C) return (0) ; /* bfs of C=A' to find R0,R3,C3 */
+ Ap = C->p ; Ai = C->i ;
+ while (head < tail) /* while queue is not empty */
+ {
+ j = queue [head++] ; /* get the head of the queue */
+ for (p = Ap [j] ; p < Ap [j+1] ; p++)
+ {
+ i = Ai [p] ;
+ if (wi [i] >= 0) continue ; /* skip if i is marked */
+ wi [i] = mark ; /* i in set R1 (C3 if transpose) */
+ j2 = jmatch [i] ; /* traverse alternating path to j2 */
+ if (wj [j2] >= 0) continue ;/* skip j2 if it is marked */
+ wj [j2] = mark ; /* j2 in set C1 (R3 if transpose) */
+ queue [tail++] = j2 ; /* add j2 to queue */
+ }
+ }
+ if (mark != 1) cs_spfree (C) ; /* free A' if it was created */
+ return (1) ;
+}
+
+/* collect matched rows and columns into P and Q */
+static void cs_matched (int m, const int *wi, const int *jmatch, int *P, int *Q,
+ int *cc, int *rr, int set, int mark)
+{
+ int kc = cc [set], i ;
+ int kr = rr [set-1] ;
+ for (i = 0 ; i < m ; i++)
+ {
+ if (wi [i] != mark) continue ; /* skip if i is not in R set */
+ P [kr++] = i ;
+ Q [kc++] = jmatch [i] ;
+ }
+ cc [set+1] = kc ;
+ rr [set] = kr ;
+}
+
+
+static void cs_unmatched (int m, const int *wi, int *P, int *rr, int set)
+/*
+ Purpose:
+
+ CS_UNMATCHED collects unmatched rows into the permutation vector P.
+*/
+{
+ int i, kr = rr [set] ;
+ for (i = 0 ; i < m ; i++) if (wi [i] == 0) P [kr++] = i ;
+ rr [set+1] = kr ;
+}
+
+/* return 1 if row i is in R2 */
+static int cs_rprune (int i, int j, double aij, void *other)
+{
+ int *rr = (int *) other ;
+ return (i >= rr [1] && i < rr [2]) ;
+ (void)j; (void)aij; // unused parameters
+}
+
+/* Given A, find coarse dmperm */
+csd *cs_dmperm (const cs *A)
+{
+ int m, n, i, j, k, p, cnz, nc, *jmatch, *imatch, *wi, *wj, *Pinv, *Cp, *Ci,
+ *Ps, *Rs, nb1, nb2, *P, *Q, *cc, *rr, *R, *S, ok ;
+ cs *C ;
+ csd *D, *scc ;
+ /* --- Maximum matching ------------------------------------------------- */
+ if (!A) return (NULL) ; /* check inputs */
+ m = A->m ; n = A->n ;
+ D = cs_dalloc (m, n) ; /* allocate result */
+ if (!D) return (NULL) ;
+ P = D->P ; Q = D->Q ; R = D->R ; S = D->S ; cc = D->cc ; rr = D->rr ;
+ jmatch = cs_maxtrans (A) ; /* max transversal */
+ imatch = jmatch + m ; /* imatch = inverse of jmatch */
+ if (!jmatch) return (cs_ddone (D, NULL, jmatch, 0)) ;
+ /* --- Coarse decomposition --------------------------------------------- */
+ wi = R ; wj = S ; /* use R and S as workspace */
+ for (j = 0 ; j < n ; j++) wj [j] = -1 ; /* unmark all cols for bfs */
+ for (i = 0 ; i < m ; i++) wi [i] = -1 ; /* unmark all rows for bfs */
+ cs_bfs (A, n, wi, wj, Q, imatch, jmatch, 1) ; /* find C0, C1, R1 */
+ ok = cs_bfs (A, m, wj, wi, P, jmatch, imatch, 3) ; /* find R0, R3, C3 */
+ if (!ok) return (cs_ddone (D, NULL, jmatch, 0)) ;
+ cs_unmatched (n, wj, Q, cc, 0) ; /* unmatched set C0 */
+ cs_matched (m, wi, jmatch, P, Q, cc, rr, 1, 1) ; /* set R1 and C1 */
+ cs_matched (m, wi, jmatch, P, Q, cc, rr, 2, -1) ; /* set R2 and C2 */
+ cs_matched (m, wi, jmatch, P, Q, cc, rr, 3, 3) ; /* set R3 and C3 */
+ cs_unmatched (m, wi, P, rr, 3) ; /* unmatched set R0 */
+ cs_free (jmatch) ;
+ /* --- Fine decomposition ----------------------------------------------- */
+ Pinv = cs_pinv (P, m) ; /* Pinv=P' */
+ if (!Pinv) return (cs_ddone (D, NULL, NULL, 0)) ;
+ C = cs_permute (A, Pinv, Q, 0) ;/* C=A(P,Q) (it will hold A(R2,C2)) */
+ cs_free (Pinv) ;
+ if (!C) return (cs_ddone (D, NULL, NULL, 0)) ;
+ Cp = C->p ; Ci = C->i ;
+ nc = cc [3] - cc [2] ; /* delete cols C0, C1, and C3 from C */
+ if (cc [2] > 0) for (j = cc [2] ; j <= cc [3] ; j++) Cp [j-cc[2]] = Cp [j] ;
+ C->n = nc ;
+ if (rr [2] - rr [1] < m) /* delete rows R0, R1, and R3 from C */
+ {
+ cs_fkeep (C, cs_rprune, rr) ;
+ cnz = Cp [nc] ;
+ if (rr [1] > 0) for (p = 0 ; p < cnz ; p++) Ci [p] -= rr [1] ;
+ }
+ C->m = nc ;
+ scc = cs_scc (C) ; /* find strongly-connected components of C*/
+ if (!scc) return (cs_ddone (D, C, NULL, 0)) ;
+ /* --- Combine coarse and fine decompositions --------------------------- */
+ Ps = scc->P ; /* C(Ps,Ps) is the permuted matrix */
+ Rs = scc->R ; /* kth block is Rs[k]..Rs[k+1]-1 */
+ nb1 = scc->nb ; /* # of blocks of A(*/
+ for (k = 0 ; k < nc ; k++) wj [k] = Q [Ps [k] + cc [2]] ; /* combine */
+ for (k = 0 ; k < nc ; k++) Q [k + cc [2]] = wj [k] ;
+ for (k = 0 ; k < nc ; k++) wi [k] = P [Ps [k] + rr [1]] ;
+ for (k = 0 ; k < nc ; k++) P [k + rr [1]] = wi [k] ;
+ nb2 = 0 ; /* create the fine block partitions */
+ R [0] = 0 ;
+ S [0] = 0 ;
+ if (cc [2] > 0) nb2++ ; /* leading coarse block A (R1, [C0 C1]) */
+ for (k = 0 ; k < nb1 ; k++) /* coarse block A (R2,C2) */
+ {
+ R [nb2] = Rs [k] + rr [1] ; /* A (R2,C2) splits into nb1 fine blocks */
+ S [nb2] = Rs [k] + cc [2] ;
+ nb2++ ;
+ }
+ if (rr [2] < m)
+ {
+ R [nb2] = rr [2] ; /* trailing coarse block A ([R3 R0], C3) */
+ S [nb2] = cc [3] ;
+ nb2++ ;
+ }
+ R [nb2] = m ;
+ S [nb2] = n ;
+ D->nb = nb2 ;
+ cs_dfree (scc) ;
+ return (cs_ddone (D, C, NULL, 1)) ;
+}
+
+static int cs_tol (int i, int j, double aij, void *tol)
+{
+ return (fabs (aij) > *((double *) tol)) ;
+ (void)i; (void)j; // unused parameters
+}
+int cs_droptol (cs *A, double tol)
+{
+ return (cs_fkeep (A, &cs_tol, &tol)) ; /* keep all large entries */
+}
+
+static int cs_nonzero (int i, int j, double aij, void *other)
+{
+ return (aij != 0) ;
+ (void)i; (void)j; (void)other; // unused parameters
+}
+int cs_dropzeros (cs *A)
+{
+ return (cs_fkeep (A, &cs_nonzero, NULL)) ; /* keep all nonzero entries */
+}
+int cs_dupl (cs *A)
+/*
+ Purpose:
+
+ CS_DUPL removes duplicate entries from A.
+
+ Reference:
+
+ Timothy Davis,
+ Direct Methods for Sparse Linear Systems,
+ SIAM, Philadelphia, 2006.
+*/
+{
+ int i, j, p, q, nz = 0, n, m, *Ap, *Ai, *w ;
+ double *Ax ;
+ if (!A) return (0) ; /* check inputs */
+ m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
+ w = cs_malloc (m, sizeof (int)) ; /* get workspace */
+ if (!w) return (0) ; /* out of memory */
+ for (i = 0 ; i < m ; i++) w [i] = -1 ; /* row i not yet seen */
+ for (j = 0 ; j < n ; j++)
+ {
+ q = nz ; /* column j will start at q */
+ for (p = Ap [j] ; p < Ap [j+1] ; p++)
+ {
+ i = Ai [p] ; /* A(i,j) is nonzero */
+ if (w [i] >= q)
+ {
+ Ax [w [i]] += Ax [p] ; /* A(i,j) is a duplicate */
+ }
+ else
+ {
+ w [i] = nz ; /* record where row i occurs */
+ Ai [nz] = i ; /* keep A(i,j) */
+ Ax [nz++] = Ax [p] ;
+ }
+ }
+ Ap [j] = q ; /* record start of column j */
+ }
+ Ap [n] = nz ; /* finalize A */
+ cs_free (w) ; /* free workspace */
+ return (cs_sprealloc (A, 0)) ; /* remove extra space from A */
+}
+
+/* add an entry to a triplet matrix; return 1 if ok, 0 otherwise */
+int cs_entry (cs *T, int i, int j, double x)
+{
+ if (!T || (T->nz >= T->nzmax && !cs_sprealloc (T, 2*(T->nzmax)))) return(0);
+ if (T->x) T->x [T->nz] = x ;
+ T->i [T->nz] = i ;
+ T->p [T->nz++] = j ;
+ T->m = CS_MAX (T->m, i+1) ;
+ T->n = CS_MAX (T->n, j+1) ;
+ return (1) ;
+}
+
+/* compute the etree of A (using triu(A), or A'A without forming A'A */
+int *cs_etree (const cs *A, int ata)
+{
+ int i, k, p, m, n, inext, *Ap, *Ai, *w, *parent, *ancestor, *prev ;
+ if (!A) return (NULL) ; /* check inputs */
+ m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ;
+ parent = cs_malloc (n, sizeof (int)) ;
+ w = cs_malloc (n + (ata ? m : 0), sizeof (int)) ;
+ ancestor = w ; prev = w + n ;
+ if (!w || !parent) return (cs_idone (parent, NULL, w, 0)) ;
+ if (ata) for (i = 0 ; i < m ; i++) prev [i] = -1 ;
+ for (k = 0 ; k < n ; k++)
+ {
+ parent [k] = -1 ; /* node k has no parent yet */
+ ancestor [k] = -1 ; /* nor does k have an ancestor */
+ for (p = Ap [k] ; p < Ap [k+1] ; p++)
+ {
+ i = ata ? (prev [Ai [p]]) : (Ai [p]) ;
+ for ( ; i != -1 && i < k ; i = inext) /* traverse from i to k */
+ {
+ inext = ancestor [i] ; /* inext = ancestor of i */
+ ancestor [i] = k ; /* path compression */
+ if (inext == -1) parent [i] = k ; /* no anc., parent is k */
+ }
+ if (ata) prev [Ai [p]] = k ;
+ }
+ }
+ return (cs_idone (parent, NULL, w, 1)) ;
+}
+
+/* drop entries for which fkeep(A(i,j)) is false; return nz if OK, else -1 */
+int cs_fkeep (cs *A, int (*fkeep) (int, int, double, void *), void *other)
+{
+ int j, p, nz = 0, n, *Ap, *Ai ;
+ double *Ax ;
+ if (!A || !fkeep) return (-1) ; /* check inputs */
+ n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
+ for (j = 0 ; j < n ; j++)
+ {
+ p = Ap [j] ; /* get current location of col j */
+ Ap [j] = nz ; /* record new location of col j */
+ for ( ; p < Ap [j+1] ; p++)
+ {
+ if (fkeep (Ai [p], j, Ax ? Ax [p] : 1, other))
+ {
+ if (Ax) Ax [nz] = Ax [p] ; /* keep A(i,j) */
+ Ai [nz++] = Ai [p] ;
+ }
+ }
+ }
+ return (Ap [n] = nz) ; /* finalize A and return nnz(A) */
+}
+
+/* y = A*x+y */
+int cs_gaxpy (const cs *A, const double *x, double *y)
+{
+ int p, j, n, *Ap, *Ai ;
+ double *Ax ;
+ if (!A || !x || !y) return (0) ; /* check inputs */
+ n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
+ for (j = 0 ; j < n ; j++)
+ {
+ for (p = Ap [j] ; p < Ap [j+1] ; p++)
+ {
+ y [Ai [p]] += Ax [p] * x [j] ;
+ }
+ }
+ return (1) ;
+}
+
+/* apply the ith Householder vector to x */
+int cs_happly (const cs *V, int i, double beta, double *x)
+{
+ int p, *Vp, *Vi ;
+ double *Vx, tau = 0 ;
+ if (!V || !x) return (0) ; /* check inputs */
+ Vp = V->p ; Vi = V->i ; Vx = V->x ;
+ for (p = Vp [i] ; p < Vp [i+1] ; p++) /* tau = v'*x */
+ {
+ tau += Vx [p] * x [Vi [p]] ;
+ }
+ tau *= beta ; /* tau = beta*(v'*x) */
+ for (p = Vp [i] ; p < Vp [i+1] ; p++) /* x = x - v*tau */
+ {
+ x [Vi [p]] -= Vx [p] * tau ;
+ }
+ return (1) ;
+}
+
+/* create a Householder reflection [v,beta,s]=house(x), overwrite x with v,
+ * where (I-beta*v*v')*x = s*x. See Algo 5.1.1, Golub & Van Loan, 3rd ed. */
+double cs_house (double *x, double *beta, int n)
+{
+ double s, sigma = 0 ;
+ int i ;
+ if (!x || !beta) return (-1) ; /* check inputs */
+ for (i = 1 ; i < n ; i++) sigma += x [i] * x [i] ;
+ if (sigma == 0)
+ {
+ s = fabs (x [0]) ; /* s = |x(0)| */
+ (*beta) = (x [0] <= 0) ? 2 : 0 ;
+ x [0] = 1 ;
+ }
+ else
+ {
+ s = sqrt (x [0] * x [0] + sigma) ; /* s = norm (x) */
+ x [0] = (x [0] <= 0) ? (x [0] - s) : (-sigma / (x [0] + s)) ;
+ (*beta) = -1. / (s * x [0]) ;
+ }
+ return (s) ;
+}
+
+/* x(P) = b, for dense vectors x and b; P=NULL denotes identity */
+int cs_ipvec (int n, const int *P, const double *b, double *x)
+{
+ int k ;
+ if (!x || !b) return (0) ; /* check inputs */
+ for (k = 0 ; k < n ; k++) x [P ? P [k] : k] = b [k] ;
+ return (1) ;
+}
+cs *cs_load ( FILE *f )
+/*
+ Purpose:
+
+ CS_LOAD loads a triplet matrix from a file.
+
+ Reference:
+
+ Timothy Davis,
+ Direct Methods for Sparse Linear Systems,
+ SIAM, Philadelphia, 2006.
+*/
+{
+ int i, j ;
+ double x ;
+ cs *T ;
+ if (!f) return (NULL) ;
+ T = cs_spalloc (0, 0, 1, 1, 1) ;
+ while (fscanf (f, "%d %d %lg\n", &i, &j, &x) == 3)
+ {
+ if (!cs_entry (T, i, j, x)) return (cs_spfree (T)) ;
+ }
+ return (T) ;
+}
+int cs_lsolve ( const cs *L, double *x )
+/*
+ Purpose:
+
+ CS_LSOLVE solves L*x=b.
+
+ Discussion:
+
+ On input, X contains the right hand side, and on output, the solution.
+
+ Reference:
+
+ Timothy Davis,
+ Direct Methods for Sparse Linear Systems,
+ SIAM, Philadelphia, 2006.
+*/
+{
+ int p, j, n, *Lp, *Li ;
+ double *Lx ;
+ if (!L || !x) return (0) ; /* check inputs */
+ n = L->n ; Lp = L->p ; Li = L->i ; Lx = L->x ;
+ for (j = 0 ; j < n ; j++)
+ {
+ x [j] /= Lx [Lp [j]] ;
+ for (p = Lp [j]+1 ; p < Lp [j+1] ; p++)
+ {
+ x [Li [p]] -= Lx [p] * x [j] ;
+ }
+ }
+ return (1) ;
+}
+int cs_ltsolve ( const cs *L, double *x )
+/*
+ Purpose:
+
+ CS_LTSOLVE solves L'*x=b.
+
+ Discussion:
+
+ On input, X contains the right hand side, and on output, the solution.
+
+ Reference:
+
+ Timothy Davis,
+ Direct Methods for Sparse Linear Systems,
+ SIAM, Philadelphia, 2006.
+*/
+{
+ int p, j, n, *Lp, *Li ;
+ double *Lx ;
+ if (!L || !x) return (0) ; /* check inputs */
+ n = L->n ; Lp = L->p ; Li = L->i ; Lx = L->x ;
+ for (j = n-1 ; j >= 0 ; j--)
+ {
+ for (p = Lp [j]+1 ; p < Lp [j+1] ; p++)
+ {
+ x [j] -= Lx [p] * x [Li [p]] ;
+ }
+ x [j] /= Lx [Lp [j]] ;
+ }
+ return (1) ;
+}
+
+/* [L,U,Pinv]=lu(A, [Q lnz unz]). lnz and unz can be guess */
+csn *cs_lu (const cs *A, const css *S, double tol)
+{
+ cs *L, *U ;
+ csn *N ;
+ double pivot, *Lx, *Ux, *x, a, t ;
+ int *Lp, *Li, *Up, *Ui, *Pinv, *xi, *Q, n, ipiv, k, top, p, i, col, lnz,unz;
+ if (!A || !S) return (NULL) ; /* check inputs */
+ n = A->n ;
+ Q = S->Q ; lnz = S->lnz ; unz = S->unz ;
+ x = cs_malloc (n, sizeof (double)) ;
+ xi = cs_malloc (2*n, sizeof (int)) ;
+ N = cs_calloc (1, sizeof (csn)) ;
+ if (!x || !xi || !N) return (cs_ndone (N, NULL, xi, x, 0)) ;
+ N->L = L = cs_spalloc (n, n, lnz, 1, 0) ; /* initial L and U */
+ N->U = U = cs_spalloc (n, n, unz, 1, 0) ;
+ N->Pinv = Pinv = cs_malloc (n, sizeof (int)) ;
+ if (!L || !U || !Pinv) return (cs_ndone (N, NULL, xi, x, 0)) ;
+ Lp = L->p ; Up = U->p ;
+ for (i = 0 ; i < n ; i++) x [i] = 0 ; /* clear workspace */
+ for (i = 0 ; i < n ; i++) Pinv [i] = -1 ; /* no rows pivotal yet */
+ for (k = 0 ; k <= n ; k++) Lp [k] = 0 ; /* no cols of L yet */
+ lnz = unz = 0 ;
+ for (k = 0 ; k < n ; k++) /* compute L(:,k) and U(:,k) */
+ {
+ /* --- Triangular solve --------------------------------------------- */
+ Lp [k] = lnz ; /* L(:,k) starts here */
+ Up [k] = unz ; /* U(:,k) starts here */
+ if ((lnz + n > L->nzmax && !cs_sprealloc (L, 2*L->nzmax + n)) ||
+ (unz + n > U->nzmax && !cs_sprealloc (U, 2*U->nzmax + n)))
+ {
+ return (cs_ndone (N, NULL, xi, x, 0)) ;
+ }
+ Li = L->i ; Lx = L->x ; Ui = U->i ; Ux = U->x ;
+ col = Q ? (Q [k]) : k ;
+ top = cs_splsolve (L, A, col, xi, x, Pinv) ; /* x = L\A(:,col) */
+ /* --- Find pivot --------------------------------------------------- */
+ ipiv = -1 ;
+ a = -1 ;
+ for (p = top ; p < n ; p++)
+ {
+ i = xi [p] ; /* x(i) is nonzero */
+ if (Pinv [i] < 0) /* row i is not pivotal */
+ {
+ if ((t = fabs (x [i])) > a)
+ {
+ a = t ; /* largest pivot candidate so far */
+ ipiv = i ;
+ }
+ }
+ else /* x(i) is the entry U(Pinv[i],k) */
+ {
+ Ui [unz] = Pinv [i] ;
+ Ux [unz++] = x [i] ;
+ }
+ }
+ if (ipiv == -1 || a <= 0) return (cs_ndone (N, NULL, xi, x, 0)) ;
+ if (Pinv [col] < 0 && fabs (x [col]) >= a*tol) ipiv = col ;
+ /* --- Divide by pivot ---------------------------------------------- */
+ pivot = x [ipiv] ; /* the chosen pivot */
+ Ui [unz] = k ; /* last entry in U(:,k) is U(k,k) */
+ Ux [unz++] = pivot ;
+ Pinv [ipiv] = k ; /* ipiv is the kth pivot row */
+ Li [lnz] = ipiv ; /* first entry in L(:,k) is L(k,k) = 1 */
+ Lx [lnz++] = 1 ;
+ for (p = top ; p < n ; p++) /* L(k+1:n,k) = x / pivot */
+ {
+ i = xi [p] ;
+ if (Pinv [i] < 0) /* x(i) is an entry in L(:,k) */
+ {
+ Li [lnz] = i ; /* save unpermuted row in L */
+ Lx [lnz++] = x [i] / pivot ; /* scale pivot column */
+ }
+ x [i] = 0 ; /* x [0..n-1] = 0 for next k */
+ }
+ }
+ /* --- Finalize L and U ------------------------------------------------- */
+ Lp [n] = lnz ;
+ Up [n] = unz ;
+ Li = L->i ; /* fix row indices of L for final Pinv */
+ for (p = 0 ; p < lnz ; p++) Li [p] = Pinv [Li [p]] ;
+ cs_sprealloc (L, 0) ; /* remove extra space from L and U */
+ cs_sprealloc (U, 0) ;
+ return (cs_ndone (N, NULL, xi, x, 1)) ; /* success */
+}
+
+/* x=A\b where A is unsymmetric; b overwritten with solution */
+int cs_lusol (const cs *A, double *b, int order, double tol)
+{
+ double *x ;
+ css *S ;
+ csn *N ;
+ int n, ok ;
+ if (!A || !b) return (0) ; /* check inputs */
+ n = A->n ;
+ S = cs_sqr (A, order, 0) ; /* ordering and symbolic analysis */
+ N = cs_lu (A, S, tol) ; /* numeric LU factorization */
+ x = cs_malloc (n, sizeof (double)) ;
+ ok = (S && N && x) ;
+ if (ok)
+ {
+ cs_ipvec (n, N->Pinv, b, x) ; /* x = P*b */
+ cs_lsolve (N->L, x) ; /* x = L\x */
+ cs_usolve (N->U, x) ; /* x = U\x */
+ cs_ipvec (n, S->Q, x, b) ; /* b = Q*x */
+ }
+ cs_free (x) ;
+ cs_sfree (S) ;
+ cs_nfree (N) ;
+ return (ok) ;
+}
+
+#ifdef MATLAB_MEX_FILE
+#define malloc mxMalloc
+#define free mxFree
+#define realloc mxRealloc
+#define calloc mxCalloc
+#endif
+
+/* wrapper for malloc */
+void *cs_malloc (int n, size_t size)
+{
+ return (CS_OVERFLOW (n,size) ? NULL : malloc (CS_MAX (n,1) * size)) ;
+}
+
+/* wrapper for calloc */
+void *cs_calloc (int n, size_t size)
+{
+ return (CS_OVERFLOW (n,size) ? NULL : calloc (CS_MAX (n,1), size)) ;
+}
+
+/* wrapper for free */
+void *cs_free (void *p)
+{
+ if (p) free (p) ; /* free p if it is not already NULL */
+ return (NULL) ; /* return NULL to simplify the use of cs_free */
+}
+
+/* wrapper for realloc */
+void *cs_realloc (void *p, int n, size_t size, int *ok)
+{
+ void *p2 ;
+ *ok = !CS_OVERFLOW (n,size) ; /* guard against int overflow */
+ if (!(*ok)) return (p) ; /* p unchanged if n too large */
+ p2 = realloc (p, CS_MAX (n,1) * size) ; /* realloc the block */
+ *ok = (p2 != NULL) ;
+ return ((*ok) ? p2 : p) ; /* return original p if failure */
+}
+
+/* find an augmenting path starting at column k and extend the match if found */
+static void cs_augment (int k, const cs *A, int *jmatch, int *cheap, int *w,
+ int *js, int *is, int *ps)
+{
+ int found = 0, p, i = -1, *Ap = A->p, *Ai = A->i, head = 0, j ;
+ js [0] = k ; /* start with just node k in jstack */
+ while (head >= 0)
+ {
+ /* --- Start (or continue) depth-first-search at node j ------------- */
+ j = js [head] ; /* get j from top of jstack */
+ if (w [j] != k) /* 1st time j visited for kth path */
+ {
+ w [j] = k ; /* mark j as visited for kth path */
+ for (p = cheap [j] ; p < Ap [j+1] && !found ; p++)
+ {
+ i = Ai [p] ; /* try a cheap assignment (i,j) */
+ found = (jmatch [i] == -1) ;
+ }
+ cheap [j] = p ; /* start here next time j is traversed*/
+ if (found)
+ {
+ is [head] = i ; /* column j matched with row i */
+ break ; /* end of augmenting path */
+ }
+ ps [head] = Ap [j] ; /* no cheap match: start dfs for j */
+ }
+ /* --- Depth-first-search of neighbors of j ------------------------- */
+ for (p = ps [head] ; p < Ap [j+1] ; p++)
+ {
+ i = Ai [p] ; /* consider row i */
+ if (w [jmatch [i]] == k) continue ; /* skip jmatch [i] if marked */
+ ps [head] = p + 1 ; /* pause dfs of node j */
+ is [head] = i ; /* i will be matched with j if found */
+ js [++head] = jmatch [i] ; /* start dfs at column jmatch [i] */
+ break ;
+ }
+ if (p == Ap [j+1]) head-- ; /* node j is done; pop from stack */
+ } /* augment the match if path found: */
+ if (found) for (p = head ; p >= 0 ; p--) jmatch [is [p]] = js [p] ;
+}
+
+/* find a maximum transveral */
+int *cs_maxtrans (const cs *A) /* returns jmatch [0..m-1]; imatch [0..n-1] */
+{
+ int i, j, k, n, m, p, n2 = 0, m2 = 0, *Ap, *jimatch, *w, *cheap, *js, *is,
+ *ps, *Ai, *Cp, *jmatch, *imatch ;
+ cs *C ;
+ if (!A) return (NULL) ; /* check inputs */
+ n = A->n ; m = A->m ; Ap = A->p ; Ai = A->i ;
+ w = jimatch = cs_calloc (m+n, sizeof (int)) ; /* allocate result */
+ if (!jimatch) return (NULL) ;
+ for (j = 0 ; j < n ; j++) /* count non-empty rows and columns */
+ {
+ n2 += (Ap [j] < Ap [j+1]) ;
+ for (p = Ap [j] ; p < Ap [j+1] ; p++) w [Ai [p]] = 1 ;
+ }
+ for (i = 0 ; i < m ; i++) m2 += w [i] ;
+ C = (m2 < n2) ? cs_transpose (A,0) : ((cs *) A) ; /* transpose if needed */
+ if (!C) return (cs_idone (jimatch, (m2 < n2) ? C : NULL, NULL, 0)) ;
+ n = C->n ; m = C->m ; Cp = C->p ;
+ jmatch = (m2 < n2) ? jimatch + n : jimatch ;
+ imatch = (m2 < n2) ? jimatch : jimatch + m ;
+ w = cs_malloc (5*n, sizeof (int)) ; /* allocate workspace */
+ if (!w) return (cs_idone (jimatch, (m2 < n2) ? C : NULL, w, 0)) ;
+ cheap = w + n ; js = w + 2*n ; is = w + 3*n ; ps = w + 4*n ;
+ for (j = 0 ; j < n ; j++) cheap [j] = Cp [j] ; /* for cheap assignment */
+ for (j = 0 ; j < n ; j++) w [j] = -1 ; /* all columns unflagged */
+ for (i = 0 ; i < m ; i++) jmatch [i] = -1 ; /* nothing matched yet */
+ for (k = 0 ; k < n ; k++) cs_augment (k, C, jmatch, cheap, w, js, is, ps) ;
+ for (j = 0 ; j < n ; j++) imatch [j] = -1 ; /* find row match */
+ for (i = 0 ; i < m ; i++) if (jmatch [i] >= 0) imatch [jmatch [i]] = i ;
+ return (cs_idone (jimatch, (m2 < n2) ? C : NULL, w, 1)) ;
+}
+
+/* C = A*B */
+cs *cs_multiply (const cs *A, const cs *B)
+{
+ int p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values, *Bi ;
+ double *x, *Bx, *Cx ;
+ cs *C ;
+ if (!A || !B) return (NULL) ; /* check inputs */
+ m = A->m ; anz = A->p [A->n] ;
+ n = B->n ; Bp = B->p ; Bi = B->i ; Bx = B->x ; bnz = Bp [n] ;
+ w = cs_calloc (m, sizeof (int)) ;
+ values = (A->x != NULL) && (Bx != NULL) ;
+ x = values ? cs_malloc (m, sizeof (double)) : NULL ;
+ C = cs_spalloc (m, n, anz + bnz, values, 0) ;
+ if (!C || !w || (values && !x)) return (cs_done (C, w, x, 0)) ;
+ Cp = C->p ;
+ for (j = 0 ; j < n ; j++)
+ {
+ if (nz + m > C->nzmax && !cs_sprealloc (C, 2*(C->nzmax)+m))
+ {
+ return (cs_done (C, w, x, 0)) ; /* out of memory */
+ }
+ Ci = C->i ; Cx = C->x ; /* C may have been reallocated */
+ Cp [j] = nz ; /* column j of C starts here */
+ for (p = Bp [j] ; p < Bp [j+1] ; p++)
+ {
+ nz = cs_scatter (A, Bi [p], Bx ? Bx [p] : 1, w, x, j+1, C, nz) ;
+ }
+ if (values) for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Ci [p]] ;
+ }
+ Cp [n] = nz ; /* finalize the last column of C */
+ cs_sprealloc (C, 0) ; /* remove extra space from C */
+ return (cs_done (C, w, x, 1)) ; /* success; free workspace, return C */
+}
+
+/* 1-norm of a sparse matrix = max (sum (abs (A))), largest column sum */
+double cs_norm (const cs *A)
+{
+ int p, j, n, *Ap ;
+ double *Ax, norm = 0, s ;
+ if (!A || !A->x) return (-1) ; /* check inputs */
+ n = A->n ; Ap = A->p ; Ax = A->x ;
+ for (j = 0 ; j < n ; j++)
+ {
+ for (s = 0, p = Ap [j] ; p < Ap [j+1] ; p++) s += fabs (Ax [p]) ;
+ norm = CS_MAX (norm, s) ;
+ }
+ return (norm) ;
+}
+
+/* C = A(P,Q) where P and Q are permutations of 0..m-1 and 0..n-1. */
+cs *cs_permute (const cs *A, const int *Pinv, const int *Q, int values)
+{
+ int p, j, k, nz = 0, m, n, *Ap, *Ai, *Cp, *Ci ;
+ double *Cx, *Ax ;
+ cs *C ;
+ if (!A) return (NULL) ; /* check inputs */
+ m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
+ C = cs_spalloc (m, n, Ap [n], values && Ax != NULL, 0) ;
+ if (!C) return (cs_done (C, NULL, NULL, 0)) ; /* out of memory */
+ Cp = C->p ; Ci = C->i ; Cx = C->x ;
+ for (k = 0 ; k < n ; k++)
+ {
+ Cp [k] = nz ; /* column k of C is column Q[k] of A */
+ j = Q ? (Q [k]) : k ;
+ for (p = Ap [j] ; p < Ap [j+1] ; p++)
+ {
+ if (Cx) Cx [nz] = Ax [p] ; /* row i of A is row Pinv[i] of C */
+ Ci [nz++] = Pinv ? (Pinv [Ai [p]]) : Ai [p] ;
+ }
+ }
+ Cp [n] = nz ; /* finalize the last column of C */
+ return (cs_done (C, NULL, NULL, 1)) ;
+}
+
+/* Pinv = P', or P = Pinv' */
+int *cs_pinv (int const *P, int n)
+{
+ int k, *Pinv ;
+ if (!P) return (NULL) ; /* P = NULL denotes identity */
+ Pinv = cs_malloc (n, sizeof (int)) ; /* allocate resuult */
+ if (!Pinv) return (NULL) ; /* out of memory */
+ for (k = 0 ; k < n ; k++) Pinv [P [k]] = k ;/* invert the permutation */
+ return (Pinv) ; /* return result */
+}
+
+/* post order a forest */
+int *cs_post (int n, const int *parent)
+{
+ int j, k = 0, *post, *w, *head, *next, *stack ;
+ if (!parent) return (NULL) ; /* check inputs */
+ post = cs_malloc (n, sizeof (int)) ; /* allocate result */
+ w = cs_malloc (3*n, sizeof (int)) ; /* 3*n workspace */
+ head = w ; next = w + n ; stack = w + 2*n ;
+ if (!w || !post) return (cs_idone (post, NULL, w, 0)) ;
+ for (j = 0 ; j < n ; j++) head [j] = -1 ; /* empty link lists */
+ for (j = n-1 ; j >= 0 ; j--) /* traverse nodes in reverse order*/
+ {
+ if (parent [j] == -1) continue ; /* j is a root */
+ next [j] = head [parent [j]] ; /* add j to list of its parent */
+ head [parent [j]] = j ;
+ }
+ for (j = 0 ; j < n ; j++)
+ {
+ if (parent [j] != -1) continue ; /* skip j if it is not a root */
+ k = cs_tdfs (j, k, head, next, post, stack) ;
+ }
+ return (cs_idone (post, NULL, w, 1)) ; /* success; free w, return post */
+}
+
+/* print a sparse matrix */
+int cs_print (const cs *A, int brief)
+{
+ int p, j, m, n, nzmax, nz, *Ap, *Ai ;
+ double *Ax ;
+ if (!A) { printf ("(null)\n") ; return (0) ; }
+ m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
+ nzmax = A->nzmax ; nz = A->nz ;
+ printf ("CSparse Version %d.%d.%d, %s. %s\n", CS_VER, CS_SUBVER,
+ CS_SUBSUB, CS_DATE, CS_COPYRIGHT) ;
+ if (nz < 0)
+ {
+ printf ("%d-by-%d, nzmax: %d nnz: %d, 1-norm: %g\n", m, n, nzmax,
+ Ap [n], cs_norm (A)) ;
+ for (j = 0 ; j < n ; j++)
+ {
+ printf (" col %d : locations %d to %d\n", j, Ap [j], Ap [j+1]-1);
+ for (p = Ap [j] ; p < Ap [j+1] ; p++)
+ {
+ printf (" %d : %g\n", Ai [p], Ax ? Ax [p] : 1) ;
+ if (brief && p > 20) { printf (" ...\n") ; return (1) ; }
+ }
+ }
+ }
+ else
+ {
+ printf ("triplet: %d-by-%d, nzmax: %d nnz: %d\n", m, n, nzmax, nz) ;
+ for (p = 0 ; p < nz ; p++)
+ {
+ printf (" %d %d : %g\n", Ai [p], Ap [p], Ax ? Ax [p] : 1) ;
+ if (brief && p > 20) { printf (" ...\n") ; return (1) ; }
+ }
+ }
+ return (1) ;
+}
+
+/* x = b(P), for dense vectors x and b; P=NULL denotes identity */
+int cs_pvec (int n, const int *P, const double *b, double *x)
+{
+ int k ;
+ if (!x || !b) return (0) ; /* check inputs */
+ for (k = 0 ; k < n ; k++) x [k] = b [P ? P [k] : k] ;
+ return (1) ;
+}
+
+/* sparse QR factorization [V,beta,p,R] = qr (A) */
+csn *cs_qr (const cs *A, const css *S)
+{
+ double *Rx, *Vx, *Ax, *Beta, *x ;
+ int i, k, p, m, n, vnz, p1, top, m2, len, col, rnz, *s, *leftmost, *Ap,
+ *Ai, *parent, *Rp, *Ri, *Vp, *Vi, *w, *Pinv, *Q ;
+ cs *R, *V ;
+ csn *N ;
+ if (!A || !S || !S->parent || !S->Pinv) return (NULL) ; /* check inputs */
+ m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
+ Q = S->Q ; parent = S->parent ; Pinv = S->Pinv ; m2 = S->m2 ;
+ vnz = S->lnz ; rnz = S->unz ;
+ leftmost = Pinv + m + n ;
+ w = cs_malloc (m2+n, sizeof (int)) ;
+ x = cs_malloc (m2, sizeof (double)) ;
+ N = cs_calloc (1, sizeof (csn)) ;
+ if (!w || !x || !N) return (cs_ndone (N, NULL, w, x, 0)) ;
+ s = w + m2 ; /* size n */
+ for (k = 0 ; k < m2 ; k++) x [k] = 0 ; /* clear workspace x */
+ N->L = V = cs_spalloc (m2, n, vnz, 1, 0) ; /* allocate V */
+ N->U = R = cs_spalloc (m2, n, rnz, 1, 0) ; /* allocate R, m2-by-n */
+ N->B = Beta = cs_malloc (n, sizeof (double)) ;
+ if (!R || !V || !Beta) return (cs_ndone (N, NULL, w, x, 0)) ;
+ Rp = R->p ; Ri = R->i ; Rx = R->x ;
+ Vp = V->p ; Vi = V->i ; Vx = V->x ;
+ for (i = 0 ; i < m2 ; i++) w [i] = -1 ; /* clear w, to mark nodes */
+ rnz = 0 ; vnz = 0 ;
+ for (k = 0 ; k < n ; k++) /* compute V and R */
+ {
+ Rp [k] = rnz ; /* R(:,k) starts here */
+ Vp [k] = p1 = vnz ; /* V(:,k) starts here */
+ w [k] = k ; /* add V(k,k) to pattern of V */
+ Vi [vnz++] = k ;
+ top = n ;
+ col = Q ? Q [k] : k ;
+ for (p = Ap [col] ; p < Ap [col+1] ; p++) /* find R(:,k) pattern */
+ {
+ i = leftmost [Ai [p]] ; /* i = min(find(A(i,Q))) */
+ for (len = 0 ; w [i] != k ; i = parent [i]) /* traverse up to k */
+ {
+ s [len++] = i ;
+ w [i] = k ;
+ }
+ while (len > 0) s [--top] = s [--len] ; /* push path on stack */
+ i = Pinv [Ai [p]] ; /* i = permuted row of A(:,col) */
+ x [i] = Ax [p] ; /* x (i) = A(.,col) */
+ if (i > k && w [i] < k) /* pattern of V(:,k) = x (k+1:m) */
+ {
+ Vi [vnz++] = i ; /* add i to pattern of V(:,k) */
+ w [i] = k ;
+ }
+ }
+ for (p = top ; p < n ; p++) /* for each i in pattern of R(:,k) */
+ {
+ i = s [p] ; /* R(i,k) is nonzero */
+ cs_happly (V, i, Beta [i], x) ; /* apply (V(i),Beta(i)) to x */
+ Ri [rnz] = i ; /* R(i,k) = x(i) */
+ Rx [rnz++] = x [i] ;
+ x [i] = 0 ;
+ if (parent [i] == k) vnz = cs_scatter (V, i, 0, w, NULL, k, V, vnz);
+ }
+ for (p = p1 ; p < vnz ; p++) /* gather V(:,k) = x */
+ {
+ Vx [p] = x [Vi [p]] ;
+ x [Vi [p]] = 0 ;
+ }
+ Ri [rnz] = k ; /* R(k,k) = norm (x) */
+ Rx [rnz++] = cs_house (Vx+p1, Beta+k, vnz-p1) ; /* [v,beta]=house(x) */
+ }
+ Rp [n] = rnz ; /* finalize R */
+ Vp [n] = vnz ; /* finalize V */
+ return (cs_ndone (N, NULL, w, x, 1)) ; /* success */
+}
+
+/* x=A\b where A can be rectangular; b overwritten with solution */
+int cs_qrsol (const cs *A, double *b, int order)
+{
+ double *x ;
+ css *S ;
+ csn *N ;
+ cs *AT = NULL ;
+ int k, m, n, ok ;
+ if (!A || !b) return (0) ; /* check inputs */
+ n = A->n ;
+ m = A->m ;
+ if (m >= n)
+ {
+ S = cs_sqr (A, order, 1) ; /* ordering and symbolic analysis */
+ N = cs_qr (A, S) ; /* numeric QR factorization */
+ x = cs_calloc (S ? S->m2 : 1, sizeof (double)) ;
+ ok = (S && N && x) ;
+ if (ok)
+ {
+ cs_ipvec (m, S->Pinv, b, x) ; /* x(0:m-1) = P*b(0:m-1) */
+ for (k = 0 ; k < n ; k++) /* apply Householder refl. to x */
+ {
+ cs_happly (N->L, k, N->B [k], x) ;
+ }
+ cs_usolve (N->U, x) ; /* x = R\x */
+ cs_ipvec (n, S->Q, x, b) ; /* b(0:n-1) = Q*x (permutation) */
+ }
+ }
+ else
+ {
+ AT = cs_transpose (A, 1) ; /* Ax=b is underdetermined */
+ S = cs_sqr (AT, order, 1) ; /* ordering and symbolic analysis */
+ N = cs_qr (AT, S) ; /* numeric QR factorization of A' */
+ x = cs_calloc (S ? S->m2 : 1, sizeof (double)) ;
+ ok = (AT && S && N && x) ;
+ if (ok)
+ {
+ cs_pvec (m, S->Q, b, x) ; /* x(0:m-1) = Q'*b (permutation) */
+ cs_utsolve (N->U, x) ; /* x = R'\x */
+ for (k = m-1 ; k >= 0 ; k--) /* apply Householder refl. to x */
+ {
+ cs_happly (N->L, k, N->B [k], x) ;
+ }
+ cs_pvec (n, S->Pinv, x, b) ; /* b (0:n-1) = P'*x */
+ }
+ }
+ cs_free (x) ;
+ cs_sfree (S) ;
+ cs_nfree (N) ;
+ cs_spfree (AT) ;
+ return (ok) ;
+}
+
+/* xi [top...n-1] = nodes reachable from graph of L*P' via nodes in B(:,k).
+ * xi [n...2n-1] used as workspace */
+int cs_reach (cs *L, const cs *B, int k, int *xi, const int *Pinv)
+{
+ int p, n, top, *Bp, *Bi, *Lp ;
+ if (!L || !B || !xi) return (-1) ;
+ n = L->n ; Bp = B->p ; Bi = B->i ; Lp = L->p ;
+ top = n ;
+ for (p = Bp [k] ; p < Bp [k+1] ; p++)
+ {
+ if (!CS_MARKED (Lp, Bi [p])) /* start a dfs at unmarked node i */
+ {
+ top = cs_dfs (Bi [p], L, top, xi, xi+n, Pinv) ;
+ }
+ }
+ for (p = top ; p < n ; p++) CS_MARK (Lp, xi [p]) ; /* restore L */
+ return (top) ;
+}
+
+/* x = x + beta * A(:,j), where x is a dense vector and A(:,j) is sparse */
+int cs_scatter (const cs *A, int j, double beta, int *w, double *x, int mark,
+ cs *C, int nz)
+{
+ int i, p, *Ap, *Ai, *Ci ;
+ double *Ax ;
+ if (!A || !w || !C) return (-1) ; /* ensure inputs are valid */
+ Ap = A->p ; Ai = A->i ; Ax = A->x ; Ci = C->i ;
+ for (p = Ap [j] ; p < Ap [j+1] ; p++)
+ {
+ i = Ai [p] ; /* A(i,j) is nonzero */
+ if (w [i] < mark)
+ {
+ w [i] = mark ; /* i is new entry in column j */
+ Ci [nz++] = i ; /* add i to pattern of C(:,j) */
+ if (x) x [i] = beta * Ax [p] ; /* x(i) = beta*A(i,j) */
+ }
+ else if (x) x [i] += beta * Ax [p] ; /* i exists in C(:,j) already */
+ }
+ return (nz) ;
+}
+
+/* find the strongly connected components of a square matrix */
+csd *cs_scc (cs *A) /* matrix A temporarily modified, then restored */
+{
+ int n, i, k, b = 0, top, *xi, *pstack, *P, *R, *Ap, *ATp ;
+ cs *AT ;
+ csd *D ;
+ if (!A) return (NULL) ;
+ n = A->n ; Ap = A->p ;
+ D = cs_dalloc (n, 0) ;
+ AT = cs_transpose (A, 0) ; /* AT = A' */
+ xi = cs_malloc (2*n, sizeof (int)) ; /* allocate workspace */
+ pstack = xi + n ;
+ if (!D || !AT || !xi) return (cs_ddone (D, AT, xi, 0)) ;
+ P = D->P ; R = D->R ; ATp = AT->p ;
+ top = n ;
+ for (i = 0 ; i < n ; i++) /* first dfs(A) to find finish times (xi) */
+ {
+ if (!CS_MARKED (Ap,i)) top = cs_dfs (i, A, top, xi, pstack, NULL) ;
+ }
+ for (i = 0 ; i < n ; i++) CS_MARK (Ap, i) ; /* restore A; unmark all nodes*/
+ top = n ;
+ b = n ;
+ for (k = 0 ; k < n ; k++) /* dfs(A') to find strongly connnected comp. */
+ {
+ i = xi [k] ; /* get i in reverse order of finish times */
+ if (CS_MARKED (ATp,i)) continue ; /* skip node i if already ordered */
+ R [b--] = top ; /* node i is the start of a component in P */
+ top = cs_dfs (i, AT, top, P, pstack, NULL) ;
+ }
+ R [b] = 0 ; /* first block starts at zero; shift R up */
+ for (k = b ; k <= n ; k++) R [k-b] = R [k] ;
+ D->nb = R [n+1] = b = n-b ; /* b = # of strongly connected components */
+ return (cs_ddone (D, AT, xi, 1)) ;
+}
+
+/* ordering and symbolic analysis for a Cholesky factorization */
+css *cs_schol (const cs *A, int order)
+{
+ int n, *c, *post, *P ;
+ cs *C ;
+ css *S ;
+ if (!A) return (NULL) ; /* check inputs */
+ n = A->n ;
+ S = cs_calloc (1, sizeof (css)) ; /* allocate symbolic analysis */
+ if (!S) return (NULL) ; /* out of memory */
+ P = cs_amd (A, order) ; /* P = amd(A+A'), or natural */
+ S->Pinv = cs_pinv (P, n) ; /* find inverse permutation */
+ cs_free (P) ;
+ if (order >= 0 && !S->Pinv) return (cs_sfree (S)) ;
+ C = cs_symperm (A, S->Pinv, 0) ; /* C = spones(triu(A(P,P))) */
+ S->parent = cs_etree (C, 0) ; /* find etree of C */
+ post = cs_post (n, S->parent) ; /* postorder the etree */
+ c = cs_counts (C, S->parent, post, 0) ; /* find column counts of chol(C) */
+ cs_free (post) ;
+ cs_spfree (C) ;
+ S->cp = cs_malloc (n+1, sizeof (int)) ; /* find column pointers for L */
+ S->unz = S->lnz = cs_cumsum (S->cp, c, n) ;
+ cs_free (c) ;
+ return ((S->lnz >= 0) ? S : cs_sfree (S)) ;
+}
+
+/* solve Lx=b(:,k), leaving pattern in xi[top..n-1], values scattered in x. */
+int cs_splsolve (cs *L, const cs *B, int k, int *xi, double *x, const int *Pinv)
+{
+ int j, jnew, p, px, top, n, *Lp, *Li, *Bp, *Bi ;
+ double *Lx, *Bx ;
+ if (!L || !B || !xi || !x) return (-1) ;
+ Lp = L->p ; Li = L->i ; Lx = L->x ; n = L->n ;
+ Bp = B->p ; Bi = B->i ; Bx = B->x ;
+ top = cs_reach (L, B, k, xi, Pinv) ; /* xi[top..n-1]=Reach(B(:,k)) */
+ for (p = top ; p < n ; p++) x [xi [p]] = 0 ;/* clear x */
+ for (p = Bp [k] ; p < Bp [k+1] ; p++) x [Bi [p]] = Bx [p] ; /* scatter B */
+ for (px = top ; px < n ; px++)
+ {
+ j = xi [px] ; /* x(j) is nonzero */
+ jnew = Pinv ? (Pinv [j]) : j ; /* j is column jnew of L */
+ if (jnew < 0) continue ; /* column jnew is empty */
+ for (p = Lp [jnew]+1 ; p < Lp [jnew+1] ; p++)
+ {
+ x [Li [p]] -= Lx [p] * x [j] ; /* x(i) -= L(i,j) * x(j) */
+ }
+ }
+ return (top) ; /* return top of stack */
+}
+
+/* compute vnz, Pinv, leftmost, m2 from A and parent */
+static int *cs_vcount (const cs *A, const int *parent, int *m2, int *vnz)
+{
+ int i, k, p, pa, n = A->n, m = A->m, *Ap = A->p, *Ai = A->i ;
+ int *Pinv = cs_malloc (2*m+n, sizeof (int)), *leftmost = Pinv + m + n ;
+ int *w = cs_malloc (m+3*n, sizeof (int)) ;
+ int *next = w, *head = w + m, *tail = w + m + n, *nque = w + m + 2*n ;
+ if (!Pinv || !w) return (cs_idone (Pinv, NULL, w, 0)) ;
+ for (k = 0 ; k < n ; k++) head [k] = -1 ; /* queue k is empty */
+ for (k = 0 ; k < n ; k++) tail [k] = -1 ;
+ for (k = 0 ; k < n ; k++) nque [k] = 0 ;
+ for (i = 0 ; i < m ; i++) leftmost [i] = -1 ;
+ for (k = n-1 ; k >= 0 ; k--)
+ {
+ for (p = Ap [k] ; p < Ap [k+1] ; p++)
+ {
+ leftmost [Ai [p]] = k ; /* leftmost[i] = min(find(A(i,:)))*/
+ }
+ }
+ for (i = m-1 ; i >= 0 ; i--) /* scan rows in reverse order */
+ {
+ Pinv [i] = -1 ; /* row i is not yet ordered */
+ k = leftmost [i] ;
+ if (k == -1) continue ; /* row i is empty */
+ if (nque [k]++ == 0) tail [k] = i ; /* first row in queue k */
+ next [i] = head [k] ; /* put i at head of queue k */
+ head [k] = i ;
+ }
+ (*vnz) = 0 ;
+ (*m2) = m ;
+ for (k = 0 ; k < n ; k++) /* find row permutation and nnz(V)*/
+ {
+ i = head [k] ; /* remove row i from queue k */
+ (*vnz)++ ; /* count V(k,k) as nonzero */
+ if (i < 0) i = (*m2)++ ; /* add a fictitious row */
+ Pinv [i] = k ; /* associate row i with V(:,k) */
+ if (--nque [k] <= 0) continue ; /* skip if V(k+1:m,k) is empty */
+ (*vnz) += nque [k] ; /* nque [k] = nnz (V(k+1:m,k)) */
+ if ((pa = parent [k]) != -1) /* move all rows to parent of k */
+ {
+ if (nque [pa] == 0) tail [pa] = tail [k] ;
+ next [tail [k]] = head [pa] ;
+ head [pa] = next [i] ;
+ nque [pa] += nque [k] ;
+ }
+ }
+ for (i = 0 ; i < m ; i++) if (Pinv [i] < 0) Pinv [i] = k++ ;
+ return (cs_idone (Pinv, NULL, w, 1)) ;
+}
+
+/* symbolic analysis for QR or LU */
+css *cs_sqr (const cs *A, int order, int qr)
+{
+ int n, k, ok = 1, *post ;
+ css *S ;
+ if (!A) return (NULL) ; /* check inputs */
+ n = A->n ;
+ S = cs_calloc (1, sizeof (css)) ; /* allocate symbolic analysis */
+ if (!S) return (NULL) ; /* out of memory */
+ S->Q = cs_amd (A, order) ; /* fill-reducing ordering */
+ if (order >= 0 && !S->Q) return (cs_sfree (S)) ;
+ if (qr) /* QR symbolic analysis */
+ {
+ cs *C = (order >= 0) ? cs_permute (A, NULL, S->Q, 0) : ((cs *) A) ;
+ S->parent = cs_etree (C, 1) ; /* etree of C'*C, where C=A(:,Q) */
+ post = cs_post (n, S->parent) ;
+ S->cp = cs_counts (C, S->parent, post, 1) ; /* col counts chol(C'*C) */
+ cs_free (post) ;
+ ok = C && S->parent && S->cp ;
+ if (ok) S->Pinv = cs_vcount (C, S->parent, &(S->m2), &(S->lnz)) ;
+ ok = ok && S->Pinv ;
+ if (ok) for (S->unz = 0, k = 0 ; k < n ; k++) S->unz += S->cp [k] ;
+ if (order >= 0) cs_spfree (C) ;
+ }
+ else
+ {
+ S->unz = 4*(A->p [n]) + n ; /* for LU factorization only, */
+ S->lnz = S->unz ; /* guess nnz(L) and nnz(U) */
+ }
+ return (ok ? S : cs_sfree (S)) ;
+}
+
+/* C = A(p,p) where A and C are symmetric the upper part stored, Pinv not P */
+cs *cs_symperm (const cs *A, const int *Pinv, int values)
+{
+ int i, j, p, q, i2, j2, n, *Ap, *Ai, *Cp, *Ci, *w ;
+ double *Cx, *Ax ;
+ cs *C ;
+ if (!A) return (NULL) ;
+ n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
+ C = cs_spalloc (n, n, Ap [n], values && (Ax != NULL), 0) ;
+ w = cs_calloc (n, sizeof (int)) ;
+ if (!C || !w) return (cs_done (C, w, NULL, 0)) ; /* out of memory */
+ Cp = C->p ; Ci = C->i ; Cx = C->x ;
+ for (j = 0 ; j < n ; j++) /* count entries in each column of C */
+ {
+ j2 = Pinv ? Pinv [j] : j ; /* column j of A is column j2 of C */
+ for (p = Ap [j] ; p < Ap [j+1] ; p++)
+ {
+ i = Ai [p] ;
+ if (i > j) continue ; /* skip lower triangular part of A */
+ i2 = Pinv ? Pinv [i] : i ; /* row i of A is row i2 of C */
+ w [CS_MAX (i2, j2)]++ ; /* column count of C */
+ }
+ }
+ cs_cumsum (Cp, w, n) ; /* compute column pointers of C */
+ for (j = 0 ; j < n ; j++)
+ {
+ j2 = Pinv ? Pinv [j] : j ; /* column j of A is column j2 of C */
+ for (p = Ap [j] ; p < Ap [j+1] ; p++)
+ {
+ i = Ai [p] ;
+ if (i > j) continue ; /* skip lower triangular part of A*/
+ i2 = Pinv ? Pinv [i] : i ; /* row i of A is row i2 of C */
+ Ci [q = w [CS_MAX (i2, j2)]++] = CS_MIN (i2, j2) ;
+ if (Cx) Cx [q] = Ax [p] ;
+ }
+ }
+ return (cs_done (C, w, NULL, 1)) ; /* success; free workspace, return C */
+}
+
+/* depth-first search and postorder of a tree rooted at node j */
+int cs_tdfs (int j, int k, int *head, const int *next, int *post, int *stack)
+{
+ int i, p, top = 0 ;
+ if (!head || !next || !post || !stack) return (-1) ; /* check inputs */
+ stack [0] = j ; /* place j on the stack */
+ while (top >= 0) /* while (stack is not empty) */
+ {
+ p = stack [top] ; /* p = top of stack */
+ i = head [p] ; /* i = youngest child of p */
+ if (i == -1)
+ {
+ top-- ; /* p has no unordered children left */
+ post [k++] = p ; /* node p is the kth postordered node */
+ }
+ else
+ {
+ head [p] = next [i] ; /* remove i from children of p */
+ stack [++top] = i ; /* start dfs on child node i */
+ }
+ }
+ return (k) ;
+}
+
+/* C = A' */
+cs *cs_transpose (const cs *A, int values)
+{
+ int p, q, j, *Cp, *Ci, n, m, *Ap, *Ai, *w ;
+ double *Cx, *Ax ;
+ cs *C ;
+ if (!A) return (NULL) ;
+ m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
+ C = cs_spalloc (n, m, Ap [n], values && Ax, 0) ; /* allocate result */
+ w = cs_calloc (m, sizeof (int)) ;
+ if (!C || !w) return (cs_done (C, w, NULL, 0)) ; /* out of memory */
+ Cp = C->p ; Ci = C->i ; Cx = C->x ;
+ for (p = 0 ; p < Ap [n] ; p++) w [Ai [p]]++ ; /* row counts */
+ cs_cumsum (Cp, w, m) ; /* row pointers */
+ for (j = 0 ; j < n ; j++)
+ {
+ for (p = Ap [j] ; p < Ap [j+1] ; p++)
+ {
+ Ci [q = w [Ai [p]]++] = j ; /* place A(i,j) as entry C(j,i) */
+ if (Cx) Cx [q] = Ax [p] ;
+ }
+ }
+ return (cs_done (C, w, NULL, 1)) ; /* success; free w and return C */
+}
+
+/* C = compressed-column form of a triplet matrix T */
+cs *cs_triplet (const cs *T)
+{
+ int m, n, nz, p, k, *Cp, *Ci, *w, *Ti, *Tj ;
+ double *Cx, *Tx ;
+ cs *C ;
+ if (!T) return (NULL) ; /* check inputs */
+ m = T->m ; n = T->n ; Ti = T->i ; Tj = T->p ; Tx = T->x ; nz = T->nz ;
+ C = cs_spalloc (m, n, nz, Tx != NULL, 0) ; /* allocate result */
+ w = cs_calloc (n, sizeof (int)) ; /* get workspace */
+ if (!C || !w) return (cs_done (C, w, NULL, 0)) ; /* out of memory */
+ Cp = C->p ; Ci = C->i ; Cx = C->x ;
+ for (k = 0 ; k < nz ; k++) w [Tj [k]]++ ; /* column counts */
+ cs_cumsum (Cp, w, n) ; /* column pointers */
+ for (k = 0 ; k < nz ; k++)
+ {
+ Ci [p = w [Tj [k]]++] = Ti [k] ; /* A(i,j) is the pth entry in C */
+ if (Cx) Cx [p] = Tx [k] ;
+ }
+ return (cs_done (C, w, NULL, 1)) ; /* success; free w and return C */
+}
+
+/* sparse Cholesky update/downdate, L*L' + sigma*w*w' (sigma = +1 or -1) */
+int cs_updown (cs *L, int sigma, const cs *C, const int *parent)
+{
+ int p, f, j, n, *Lp, *Li, *Cp, *Ci ;
+ double *Lx, *Cx, alpha, beta = 1, delta, gamma, w1, w2, *w, beta2 = 1 ;
+ if (!L || !C || !parent) return (0) ;
+ Lp = L->p ; Li = L->i ; Lx = L->x ; n = L->n ;
+ Cp = C->p ; Ci = C->i ; Cx = C->x ;
+ if ((p = Cp [0]) >= Cp [1]) return (1) ; /* return if C empty */
+ w = cs_malloc (n, sizeof (double)) ;
+ if (!w) return (0) ;
+ f = Ci [p] ;
+ for ( ; p < Cp [1] ; p++) f = CS_MIN (f, Ci [p]) ; /* f = min (find (C)) */
+ for (j = f ; j != -1 ; j = parent [j]) w [j] = 0 ; /* clear workspace w */
+ for (p = Cp [0] ; p < Cp [1] ; p++) w [Ci [p]] = Cx [p] ; /* w = C */
+ for (j = f ; j != -1 ; j = parent [j]) /* walk path f up to root */
+ {
+ p = Lp [j] ;
+ alpha = w [j] / Lx [p] ; /* alpha = w(j) / L(j,j) */
+ beta2 = beta*beta + sigma*alpha*alpha ;
+ if (beta2 <= 0) break ; /* not positive definite */
+ beta2 = sqrt (beta2) ;
+ delta = (sigma > 0) ? (beta / beta2) : (beta2 / beta) ;
+ gamma = sigma * alpha / (beta2 * beta) ;
+ Lx [p] = delta * Lx [p] + ((sigma > 0) ? (gamma * w [j]) : 0) ;
+ beta = beta2 ;
+ for (p++ ; p < Lp [j+1] ; p++)
+ {
+ w1 = w [Li [p]] ;
+ w [Li [p]] = w2 = w1 - alpha * Lx [p] ;
+ Lx [p] = delta * Lx [p] + gamma * ((sigma > 0) ? w1 : w2) ;
+ }
+ }
+ cs_free (w) ;
+ return (beta2 > 0) ;
+}
+
+/* solve Ux=b where x and b are dense. x=b on input, solution on output. */
+int cs_usolve (const cs *U, double *x)
+{
+ int p, j, n, *Up, *Ui ;
+ double *Ux ;
+ if (!U || !x) return (0) ; /* check inputs */
+ n = U->n ; Up = U->p ; Ui = U->i ; Ux = U->x ;
+ for (j = n-1 ; j >= 0 ; j--)
+ {
+ x [j] /= Ux [Up [j+1]-1] ;
+ for (p = Up [j] ; p < Up [j+1]-1 ; p++)
+ {
+ x [Ui [p]] -= Ux [p] * x [j] ;
+ }
+ }
+ return (1) ;
+}
+
+/* allocate a sparse matrix (triplet form or compressed-column form) */
+cs *cs_spalloc (int m, int n, int nzmax, int values, int triplet)
+{
+ cs *A = cs_calloc (1, sizeof (cs)) ; /* allocate the cs struct */
+ if (!A) return (NULL) ; /* out of memory */
+ A->m = m ; /* define dimensions and nzmax */
+ A->n = n ;
+ A->nzmax = nzmax = CS_MAX (nzmax, 1) ;
+ A->nz = triplet ? 0 : -1 ; /* allocate triplet or comp.col */
+ A->p = cs_malloc (triplet ? nzmax : n+1, sizeof (int)) ;
+ A->i = cs_malloc (nzmax, sizeof (int)) ;
+ A->x = values ? cs_malloc (nzmax, sizeof (double)) : NULL ;
+ return ((!A->p || !A->i || (values && !A->x)) ? cs_spfree (A) : A) ;
+}
+
+/* change the max # of entries sparse matrix */
+int cs_sprealloc (cs *A, int nzmax)
+{
+ int ok, oki, okj = 1, okx = 1 ;
+ if (!A) return (0) ;
+ nzmax = (nzmax <= 0) ? (A->p [A->n]) : nzmax ;
+ A->i = cs_realloc (A->i, nzmax, sizeof (int), &oki) ;
+ if (A->nz >= 0) A->p = cs_realloc (A->p, nzmax, sizeof (int), &okj) ;
+ if (A->x) A->x = cs_realloc (A->x, nzmax, sizeof (double), &okx) ;
+ ok = (oki && okj && okx) ;
+ if (ok) A->nzmax = nzmax ;
+ return (ok) ;
+}
+
+/* free a sparse matrix */
+cs *cs_spfree (cs *A)
+{
+ if (!A) return (NULL) ; /* do nothing if A already NULL */
+ cs_free (A->p) ;
+ cs_free (A->i) ;
+ cs_free (A->x) ;
+ return (cs_free (A)) ; /* free the cs struct and return NULL */
+}
+
+/* free a numeric factorization */
+csn *cs_nfree (csn *N)
+{
+ if (!N) return (NULL) ; /* do nothing if N already NULL */
+ cs_spfree (N->L) ;
+ cs_spfree (N->U) ;
+ cs_free (N->Pinv) ;
+ cs_free (N->B) ;
+ return (cs_free (N)) ; /* free the csn struct and return NULL */
+}
+
+/* free a symbolic factorization */
+css *cs_sfree (css *S)
+{
+ if (!S) return (NULL) ; /* do nothing if S already NULL */
+ cs_free (S->Pinv) ;
+ cs_free (S->Q) ;
+ cs_free (S->parent) ;
+ cs_free (S->cp) ;
+ return (cs_free (S)) ; /* free the css struct and return NULL */
+}
+
+/* allocate a cs_dmperm or cs_scc result */
+csd *cs_dalloc (int m, int n)
+{
+ csd *D ;
+ D = cs_calloc (1, sizeof (csd)) ;
+ if (!D) return (NULL) ;
+ D->P = cs_malloc (m, sizeof (int)) ;
+ D->R = cs_malloc (m+6, sizeof (int)) ;
+ D->Q = cs_malloc (n, sizeof (int)) ;
+ D->S = cs_malloc (n+6, sizeof (int)) ;
+ return ((!D->P || !D->R || !D->Q || !D->S) ? cs_dfree (D) : D) ;
+}
+
+/* free a cs_dmperm or cs_scc result */
+csd *cs_dfree (csd *D)
+{
+ if (!D) return (NULL) ; /* do nothing if D already NULL */
+ cs_free (D->P) ;
+ cs_free (D->Q) ;
+ cs_free (D->R) ;
+ cs_free (D->S) ;
+ return (cs_free (D)) ;
+}
+
+/* free workspace and return a sparse matrix result */
+cs *cs_done (cs *C, void *w, void *x, int ok)
+{
+ cs_free (w) ; /* free workspace */
+ cs_free (x) ;
+ return (ok ? C : cs_spfree (C)) ; /* return result if OK, else free it */
+}
+
+/* free workspace and return int array result */
+int *cs_idone (int *p, cs *C, void *w, int ok)
+{
+ cs_spfree (C) ; /* free temporary matrix */
+ cs_free (w) ; /* free workspace */
+ return (ok ? p : cs_free (p)) ; /* return result if OK, else free it */
+}
+
+/* free workspace and return a numeric factorization (Cholesky, LU, or QR) */
+csn *cs_ndone (csn *N, cs *C, void *w, void *x, int ok)
+{
+ cs_spfree (C) ; /* free temporary matrix */
+ cs_free (w) ; /* free workspace */
+ cs_free (x) ;
+ return (ok ? N : cs_nfree (N)) ; /* return result if OK, else free it */
+}
+
+/* free workspace and return a csd result */
+csd *cs_ddone (csd *D, cs *C, void *w, int ok)
+{
+ cs_spfree (C) ; /* free temporary matrix */
+ cs_free (w) ; /* free workspace */
+ return (ok ? D : cs_dfree (D)) ; /* return result if OK, else free it */
+}
+
+/* solve U'x=b where x and b are dense. x=b on input, solution on output. */
+int cs_utsolve (const cs *U, double *x)
+{
+ int p, j, n, *Up, *Ui ;
+ double *Ux ;
+ if (!U || !x) return (0) ; /* check inputs */
+ n = U->n ; Up = U->p ; Ui = U->i ; Ux = U->x ;
+ for (j = 0 ; j < n ; j++)
+ {
+ for (p = Up [j] ; p < Up [j+1]-1 ; p++)
+ {
+ x [j] -= Ux [p] * x [Ui [p]] ;
+ }
+ x [j] /= Ux [p] ;
+ }
+ return (1) ;
+}
=====================================
extlibs/csparse/csparse.h
=====================================
@@ -0,0 +1,144 @@
+#ifndef _CS_H
+#define _CS_H
+
+#ifdef MATLAB_MEX_FILE
+#include "mex.h"
+#endif
+#define CS_VER 1 /* CSparse Version 1.2.0 */
+#define CS_SUBVER 2
+#define CS_SUBSUB 0
+#define CS_DATE "Mar 6, 2006" /* CSparse release date */
+#define CS_COPYRIGHT "Copyright (c) Timothy A. Davis, 2006"
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+/* --- primary CSparse routines and data structures ------------------------- */
+typedef struct cs_sparse /* matrix in compressed-column or triplet form */
+{
+ int nzmax ; /* maximum number of entries */
+ int m ; /* number of rows */
+ int n ; /* number of columns */
+ int *p ; /* column pointers (size n+1) or col indices (size nzmax) */
+ int *i ; /* row indices, size nzmax */
+ double *x ; /* numerical values, size nzmax */
+ int nz ; /* # of entries in triplet matrix, -1 for compressed-col */
+} cs ;
+
+cs *cs_add (const cs *A, const cs *B, double alpha, double beta) ;
+int cs_cholsol (const cs *A, double *b, int order) ;
+int cs_dupl (cs *A) ;
+int cs_entry (cs *T, int i, int j, double x) ;
+int cs_lusol (const cs *A, double *b, int order, double tol) ;
+int cs_gaxpy (const cs *A, const double *x, double *y) ;
+cs *cs_multiply (const cs *A, const cs *B) ;
+int cs_qrsol (const cs *A, double *b, int order) ;
+cs *cs_transpose (const cs *A, int values) ;
+cs *cs_triplet (const cs *T) ;
+double cs_norm (const cs *A) ;
+int cs_print (const cs *A, int brief) ;
+cs *cs_load (FILE *f) ;
+/* utilities */
+void *cs_calloc (int n, size_t size) ;
+void *cs_free (void *p) ;
+void *cs_realloc (void *p, int n, size_t size, int *ok) ;
+cs *cs_spalloc (int m, int n, int nzmax, int values, int triplet) ;
+cs *cs_spfree (cs *A) ;
+int cs_sprealloc (cs *A, int nzmax) ;
+void *cs_malloc (int n, size_t size) ;
+
+/* --- secondary CSparse routines and data structures ----------------------- */
+typedef struct cs_symbolic /* symbolic Cholesky, LU, or QR analysis */
+{
+ int *Pinv ; /* inverse row perm. for QR, fill red. perm for Chol */
+ int *Q ; /* fill-reducing column permutation for LU and QR */
+ int *parent ; /* elimination tree for Cholesky and QR */
+ int *cp ; /* column pointers for Cholesky, row counts for QR */
+ int m2 ; /* # of rows for QR, after adding fictitious rows */
+ int lnz ; /* # entries in L for LU or Cholesky; in V for QR */
+ int unz ; /* # entries in U for LU; in R for QR */
+} css ;
+
+typedef struct cs_numeric /* numeric Cholesky, LU, or QR factorization */
+{
+ cs *L ; /* L for LU and Cholesky, V for QR */
+ cs *U ; /* U for LU, R for QR, not used for Cholesky */
+ int *Pinv ; /* partial pivoting for LU */
+ double *B ; /* beta [0..n-1] for QR */
+} csn ;
+
+typedef struct cs_dmperm_results /* cs_dmperm or cs_scc output */
+{
+ int *P ; /* size m, row permutation */
+ int *Q ; /* size n, column permutation */
+ int *R ; /* size nb+1, block k is rows R[k] to R[k+1]-1 in A(P,Q) */
+ int *S ; /* size nb+1, block k is cols S[k] to S[k+1]-1 in A(P,Q) */
+ int nb ; /* # of blocks in fine dmperm decomposition */
+ int rr [5] ; /* coarse row decomposition */
+ int cc [5] ; /* coarse column decomposition */
+} csd ;
+
+int *cs_amd (const cs *A, int order) ;
+csn *cs_chol (const cs *A, const css *S) ;
+csd *cs_dmperm (const cs *A) ;
+int cs_droptol (cs *A, double tol) ;
+int cs_dropzeros (cs *A) ;
+int cs_happly (const cs *V, int i, double beta, double *x) ;
+int cs_ipvec (int n, const int *P, const double *b, double *x) ;
+int cs_lsolve (const cs *L, double *x) ;
+int cs_ltsolve (const cs *L, double *x) ;
+csn *cs_lu (const cs *A, const css *S, double tol) ;
+cs *cs_permute (const cs *A, const int *P, const int *Q, int values) ;
+int *cs_pinv (const int *P, int n) ;
+int cs_pvec (int n, const int *P, const double *b, double *x) ;
+csn *cs_qr (const cs *A, const css *S) ;
+css *cs_schol (const cs *A, int order) ;
+css *cs_sqr (const cs *A, int order, int qr) ;
+cs *cs_symperm (const cs *A, const int *Pinv, int values) ;
+int cs_usolve (const cs *U, double *x) ;
+int cs_utsolve (const cs *U, double *x) ;
+int cs_updown (cs *L, int sigma, const cs *C, const int *parent) ;
+/* utilities */
+css *cs_sfree (css *S) ;
+csn *cs_nfree (csn *N) ;
+csd *cs_dfree (csd *D) ;
+
+/* --- tertiary CSparse routines -------------------------------------------- */
+int *cs_counts (const cs *A, const int *parent, const int *post, int ata) ;
+int cs_cumsum (int *p, int *c, int n) ;
+int cs_dfs (int j, cs *L, int top, int *xi, int *pstack, const int *Pinv) ;
+int *cs_etree (const cs *A, int ata) ;
+int cs_fkeep (cs *A, int (*fkeep) (int, int, double, void *), void *other) ;
+double cs_house (double *x, double *beta, int n) ;
+int *cs_maxtrans (const cs *A) ;
+int *cs_post (int n, const int *parent) ;
+int cs_reach (cs *L, const cs *B, int k, int *xi, const int *Pinv) ;
+csd *cs_scc (cs *A) ;
+int cs_scatter (const cs *A, int j, double beta, int *w, double *x, int mark,
+ cs *C, int nz) ;
+int cs_splsolve (cs *L, const cs *B, int k, int *xi, double *x,
+ const int *Pinv) ;
+int cs_tdfs (int j, int k, int *head, const int *next, int *post,
+ int *stack) ;
+/* utilities */
+csd *cs_dalloc (int m, int n) ;
+cs *cs_done (cs *C, void *w, void *x, int ok) ;
+int *cs_idone (int *p, cs *C, void *w, int ok) ;
+csn *cs_ndone (csn *N, cs *C, void *w, void *x, int ok) ;
+csd *cs_ddone (csd *D, cs *C, void *w, int ok) ;
+
+#define CS_MAX(a,b) (((a) > (b)) ? (a) : (b))
+#define CS_MIN(a,b) (((a) < (b)) ? (a) : (b))
+#define CS_FLIP(i) (-(i)-2)
+#define CS_UNFLIP(i) (((i) < 0) ? CS_FLIP(i) : (i))
+#define CS_MARKED(Ap,j) (Ap [j] < 0)
+#define CS_MARK(Ap,j) { Ap [j] = CS_FLIP (Ap [j]) ; }
+#define CS_OVERFLOW(n,size) (n > INT_MAX / (int) size)
+
+
+#if defined(__cplusplus)
+} //extern "C"
+#endif
+
+#endif
=====================================
extlibs/csparse/ldl.c
=====================================
@@ -0,0 +1,507 @@
+/* ========================================================================== */
+/* === ldl.c: sparse LDL' factorization and solve package =================== */
+/* ========================================================================== */
+
+/* LDL: a simple set of routines for sparse LDL' factorization. These routines
+ * are not terrifically fast (they do not use dense matrix kernels), but the
+ * code is very short. The purpose is to illustrate the algorithms in a very
+ * concise manner, primarily for educational purposes. Although the code is
+ * very concise, this package is slightly faster than the built-in sparse
+ * Cholesky factorization in MATLAB 7.0 (chol), when using the same input
+ * permutation.
+ *
+ * The routines compute the LDL' factorization of a real sparse symmetric
+ * matrix A (or PAP' if a permutation P is supplied), and solve upper
+ * and lower triangular systems with the resulting L and D factors. If A is
+ * positive definite then the factorization will be accurate. A can be
+ * indefinite (with negative values on the diagonal D), but in this case no
+ * guarantee of accuracy is provided, since no numeric pivoting is performed.
+ *
+ * The n-by-n sparse matrix A is in compressed-column form. The nonzero values
+ * in column j are stored in Ax [Ap [j] ... Ap [j+1]-1], with corresponding row
+ * indices in Ai [Ap [j] ... Ap [j+1]-1]. Ap [0] = 0 is required, and thus
+ * nz = Ap [n] is the number of nonzeros in A. Ap is an int array of size n+1.
+ * The int array Ai and the double array Ax are of size nz. This data structure
+ * is identical to the one used by MATLAB, except for the following
+ * generalizations. The row indices in each column of A need not be in any
+ * particular order, although they must be in the range 0 to n-1. Duplicate
+ * entries can be present; any duplicates are summed. That is, if row index i
+ * appears twice in a column j, then the value of A (i,j) is the sum of the two
+ * entries. The data structure used here for the input matrix A is more
+ * flexible than MATLAB's, which requires sorted columns with no duplicate
+ * entries.
+ *
+ * Only the diagonal and upper triangular part of A (or PAP' if a permutation
+ * P is provided) is accessed. The lower triangular parts of the matrix A or
+ * PAP' can be present, but they are ignored.
+ *
+ * The optional input permutation is provided as an array P of length n. If
+ * P [k] = j, the row and column j of A is the kth row and column of PAP'.
+ * If P is present then the factorization is LDL' = PAP' or L*D*L' = A(P,P) in
+ * 0-based MATLAB notation. If P is not present (a null pointer) then no
+ * permutation is performed, and the factorization is LDL' = A.
+ *
+ * The lower triangular matrix L is stored in the same compressed-column
+ * form (an int Lp array of size n+1, an int Li array of size Lp [n], and a
+ * double array Lx of the same size as Li). It has a unit diagonal, which is
+ * not stored. The row indices in each column of L are always returned in
+ * ascending order, with no duplicate entries. This format is compatible with
+ * MATLAB, except that it would be more convenient for MATLAB to include the
+ * unit diagonal of L. Doing so here would add additional complexity to the
+ * code, and is thus omitted in the interest of keeping this code short and
+ * readable.
+ *
+ * The elimination tree is held in the Parent [0..n-1] array. It is normally
+ * not required by the user, but it is required by ldl_numeric. The diagonal
+ * matrix D is held as an array D [0..n-1] of size n.
+ *
+ * --------------------
+ * C-callable routines:
+ * --------------------
+ *
+ * ldl_symbolic: Given the pattern of A, computes the Lp and Parent arrays
+ * required by ldl_numeric. Takes time proportional to the number of
+ * nonzeros in L. Computes the inverse Pinv of P if P is provided.
+ * Also returns Lnz, the count of nonzeros in each column of L below
+ * the diagonal (this is not required by ldl_numeric).
+ * ldl_numeric: Given the pattern and numerical values of A, the Lp array,
+ * the Parent array, and P and Pinv if applicable, computes the
+ * pattern and numerical values of L and D.
+ * ldl_lsolve: Solves Lx=b for a dense vector b.
+ * ldl_dsolve: Solves Dx=b for a dense vector b.
+ * ldl_ltsolve: Solves L'x=b for a dense vector b.
+ * ldl_perm: Computes x=Pb for a dense vector b.
+ * ldl_permt: Computes x=P'b for a dense vector b.
+ * ldl_valid_perm: checks the validity of a permutation vector
+ * ldl_valid_matrix: checks the validity of the sparse matrix A
+ *
+ * ----------------------------
+ * Limitations of this package:
+ * ----------------------------
+ *
+ * In the interest of keeping this code simple and readable, ldl_symbolic and
+ * ldl_numeric assume their inputs are valid. You can check your own inputs
+ * prior to calling these routines with the ldl_valid_perm and ldl_valid_matrix
+ * routines. Except for the two ldl_valid_* routines, no routine checks to see
+ * if the array arguments are present (non-NULL). Like all C routines, no
+ * routine can determine if the arrays are long enough and don't overlap.
+ *
+ * The ldl_numeric does check the numerical factorization, however. It returns
+ * n if the factorization is successful. If D (k,k) is zero, then k is
+ * returned, and L is only partially computed.
+ *
+ * No pivoting to control fill-in is performed, which is often critical for
+ * obtaining good performance. I recommend that you compute the permutation P
+ * using AMD or SYMAMD (approximate minimum degree ordering routines), or an
+ * appropriate graph-partitioning based ordering. See the ldldemo.m routine for
+ * an example in MATLAB, and the ldlmain.c stand-alone C program for examples of
+ * how to find P. Routines for manipulating compressed-column matrices are
+ * available in UMFPACK. AMD, SYMAMD, UMFPACK, and this LDL package are all
+ * available at http://www.cise.ufl.edu/research/sparse.
+ *
+ * -------------------------
+ * Possible simplifications:
+ * -------------------------
+ *
+ * These routines could be made even simpler with a few additional assumptions.
+ * If no input permutation were performed, the caller would have to permute the
+ * matrix first, but the computation of Pinv, and the use of P and Pinv could be
+ * removed. If only the diagonal and upper triangular part of A or PAP' are
+ * present, then the tests in the "if (i < k)" statement in ldl_symbolic and
+ * "if (i <= k)" in ldl_numeric, are always true, and could be removed (i can
+ * equal k in ldl_symbolic, but then the body of the if statement would
+ * correctly do no work since Flag [k] == k). If we could assume that no
+ * duplicate entries are present, then the statement Y [i] += Ax [p] could be
+ * replaced with Y [i] = Ax [p] in ldl_numeric.
+ *
+ * --------------------------
+ * Description of the method:
+ * --------------------------
+ *
+ * LDL computes the symbolic factorization by finding the pattern of L one row
+ * at a time. It does this based on the following theory. Consider a sparse
+ * system Lx=b, where L, x, and b, are all sparse, and where L comes from a
+ * Cholesky (or LDL') factorization. The elimination tree (etree) of L is
+ * defined as follows. The parent of node j is the smallest k > j such that
+ * L (k,j) is nonzero. Node j has no parent if column j of L is completely zero
+ * below the diagonal (j is a root of the etree in this case). The nonzero
+ * pattern of x is the union of the paths from each node i to the root, for
+ * each nonzero b (i). To compute the numerical solution to Lx=b, we can
+ * traverse the columns of L corresponding to nonzero values of x. This
+ * traversal does not need to be done in the order 0 to n-1. It can be done in
+ * any "topological" order, such that x (i) is computed before x (j) if i is a
+ * descendant of j in the elimination tree.
+ *
+ * The row-form of the LDL' factorization is shown in the MATLAB function
+ * ldlrow.m in this LDL package. Note that row k of L is found via a sparse
+ * triangular solve of L (1:k-1, 1:k-1) \ A (1:k-1, k), to use 1-based MATLAB
+ * notation. Thus, we can start with the nonzero pattern of the kth column of
+ * A (above the diagonal), follow the paths up to the root of the etree of the
+ * (k-1)-by-(k-1) leading submatrix of L, and obtain the pattern of the kth row
+ * of L. Note that we only need the leading (k-1)-by-(k-1) submatrix of L to
+ * do this. The elimination tree can be constructed as we go.
+ *
+ * The symbolic factorization does the same thing, except that it discards the
+ * pattern of L as it is computed. It simply counts the number of nonzeros in
+ * each column of L and then constructs the Lp index array when it's done. The
+ * symbolic factorization does not need to do this in topological order.
+ * Compare ldl_symbolic with the first part of ldl_numeric, and note that the
+ * while (len > 0) loop is not present in ldl_symbolic.
+ *
+ * LDL Version 1.3, Copyright (c) 2006 by Timothy A Davis,
+ * University of Florida. All Rights Reserved. Developed while on sabbatical
+ * at Stanford University and Lawrence Berkeley National Laboratory. Refer to
+ * the README file for the License. Available at
+ * http://www.cise.ufl.edu/research/sparse.
+ */
+
+#include "ldl.h"
+
+/* ========================================================================== */
+/* === ldl_symbolic ========================================================= */
+/* ========================================================================== */
+
+/* The input to this routine is a sparse matrix A, stored in column form, and
+ * an optional permutation P. The output is the elimination tree
+ * and the number of nonzeros in each column of L. Parent [i] = k if k is the
+ * parent of i in the tree. The Parent array is required by ldl_numeric.
+ * Lnz [k] gives the number of nonzeros in the kth column of L, excluding the
+ * diagonal.
+ *
+ * One workspace vector (Flag) of size n is required.
+ *
+ * If P is NULL, then it is ignored. The factorization will be LDL' = A.
+ * Pinv is not computed. In this case, neither P nor Pinv are required by
+ * ldl_numeric.
+ *
+ * If P is not NULL, then it is assumed to be a valid permutation. If
+ * row and column j of A is the kth pivot, the P [k] = j. The factorization
+ * will be LDL' = PAP', or A (p,p) in MATLAB notation. The inverse permutation
+ * Pinv is computed, where Pinv [j] = k if P [k] = j. In this case, both P
+ * and Pinv are required as inputs to ldl_numeric.
+ *
+ * The floating-point operation count of the subsequent call to ldl_numeric
+ * is not returned, but could be computed after ldl_symbolic is done. It is
+ * the sum of (Lnz [k]) * (Lnz [k] + 2) for k = 0 to n-1.
+ */
+
+void LDL_symbolic
+(
+ LDL_int n, /* A and L are n-by-n, where n >= 0 */
+ LDL_int Ap [ ], /* input of size n+1, not modified */
+ LDL_int Ai [ ], /* input of size nz=Ap[n], not modified */
+ LDL_int Lp [ ], /* output of size n+1, not defined on input */
+ LDL_int Parent [ ], /* output of size n, not defined on input */
+ LDL_int Lnz [ ], /* output of size n, not defined on input */
+ LDL_int Flag [ ], /* workspace of size n, not defn. on input or output */
+ LDL_int P [ ], /* optional input of size n */
+ LDL_int Pinv [ ] /* optional output of size n (used if P is not NULL) */
+)
+{
+ LDL_int i, k, p, kk, p2 ;
+ if (P)
+ {
+ /* If P is present then compute Pinv, the inverse of P */
+ for (k = 0 ; k < n ; k++)
+ {
+ Pinv [P [k]] = k ;
+ }
+ }
+ for (k = 0 ; k < n ; k++)
+ {
+ /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
+ Parent [k] = -1 ; /* parent of k is not yet known */
+ Flag [k] = k ; /* mark node k as visited */
+ Lnz [k] = 0 ; /* count of nonzeros in column k of L */
+ kk = (P) ? (P [k]) : (k) ; /* kth original, or permuted, column */
+ p2 = Ap [kk+1] ;
+ for (p = Ap [kk] ; p < p2 ; p++)
+ {
+ /* A (i,k) is nonzero (original or permuted A) */
+ i = (Pinv) ? (Pinv [Ai [p]]) : (Ai [p]) ;
+ if (i < k)
+ {
+ /* follow path from i to root of etree, stop at flagged node */
+ for ( ; Flag [i] != k ; i = Parent [i])
+ {
+ /* find parent of i if not yet determined */
+ if (Parent [i] == -1) Parent [i] = k ;
+ Lnz [i]++ ; /* L (k,i) is nonzero */
+ Flag [i] = k ; /* mark i as visited */
+ }
+ }
+ }
+ }
+ /* construct Lp index array from Lnz column counts */
+ Lp [0] = 0 ;
+ for (k = 0 ; k < n ; k++)
+ {
+ Lp [k+1] = Lp [k] + Lnz [k] ;
+ }
+}
+
+
+/* ========================================================================== */
+/* === ldl_numeric ========================================================== */
+/* ========================================================================== */
+
+/* Given a sparse matrix A (the arguments n, Ap, Ai, and Ax) and its symbolic
+ * analysis (Lp and Parent, and optionally P and Pinv), compute the numeric LDL'
+ * factorization of A or PAP'. The outputs of this routine are arguments Li,
+ * Lx, and D. It also requires three size-n workspaces (Y, Pattern, and Flag).
+ */
+
+LDL_int LDL_numeric /* returns n if successful, k if D (k,k) is zero */
+(
+ LDL_int n, /* A and L are n-by-n, where n >= 0 */
+ LDL_int Ap [ ], /* input of size n+1, not modified */
+ LDL_int Ai [ ], /* input of size nz=Ap[n], not modified */
+ double Ax [ ], /* input of size nz=Ap[n], not modified */
+ LDL_int Lp [ ], /* input of size n+1, not modified */
+ LDL_int Parent [ ], /* input of size n, not modified */
+ LDL_int Lnz [ ], /* output of size n, not defn. on input */
+ LDL_int Li [ ], /* output of size lnz=Lp[n], not defined on input */
+ double Lx [ ], /* output of size lnz=Lp[n], not defined on input */
+ double D [ ], /* output of size n, not defined on input */
+ double Y [ ], /* workspace of size n, not defn. on input or output */
+ LDL_int Pattern [ ],/* workspace of size n, not defn. on input or output */
+ LDL_int Flag [ ], /* workspace of size n, not defn. on input or output */
+ LDL_int P [ ], /* optional input of size n */
+ LDL_int Pinv [ ] /* optional input of size n */
+)
+{
+ double yi, l_ki ;
+ LDL_int i, k, p, kk, p2, len, top ;
+ for (k = 0 ; k < n ; k++)
+ {
+ /* compute nonzero Pattern of kth row of L, in topological order */
+ Y [k] = 0.0 ; /* Y(0:k) is now all zero */
+ top = n ; /* stack for pattern is empty */
+ Flag [k] = k ; /* mark node k as visited */
+ Lnz [k] = 0 ; /* count of nonzeros in column k of L */
+ kk = (P) ? (P [k]) : (k) ; /* kth original, or permuted, column */
+ p2 = Ap [kk+1] ;
+ for (p = Ap [kk] ; p < p2 ; p++)
+ {
+ i = (Pinv) ? (Pinv [Ai [p]]) : (Ai [p]) ; /* get A(i,k) */
+ if (i <= k)
+ {
+ Y [i] += Ax [p] ; /* scatter A(i,k) into Y (sum duplicates) */
+ for (len = 0 ; Flag [i] != k ; i = Parent [i])
+ {
+ Pattern [len++] = i ; /* L(k,i) is nonzero */
+ Flag [i] = k ; /* mark i as visited */
+ }
+ while (len > 0) Pattern [--top] = Pattern [--len] ;
+ }
+ }
+ /* compute numerical values kth row of L (a sparse triangular solve) */
+ D [k] = Y [k] ; /* get D(k,k) and clear Y(k) */
+ Y [k] = 0.0 ;
+ for ( ; top < n ; top++)
+ {
+ i = Pattern [top] ; /* Pattern [top:n-1] is pattern of L(:,k) */
+ yi = Y [i] ; /* get and clear Y(i) */
+ Y [i] = 0.0 ;
+ p2 = Lp [i] + Lnz [i] ;
+ for (p = Lp [i] ; p < p2 ; p++)
+ {
+ Y [Li [p]] -= Lx [p] * yi ;
+ }
+ l_ki = yi / D [i] ; /* the nonzero entry L(k,i) */
+ D [k] -= l_ki * yi ;
+ Li [p] = k ; /* store L(k,i) in column form of L */
+ Lx [p] = l_ki ;
+ Lnz [i]++ ; /* increment count of nonzeros in col i */
+ }
+ if (D [k] == 0.0) return (k) ; /* failure, D(k,k) is zero */
+ }
+ return (n) ; /* success, diagonal of D is all nonzero */
+}
+
+
+/* ========================================================================== */
+/* === ldl_lsolve: solve Lx=b ============================================== */
+/* ========================================================================== */
+
+void LDL_lsolve
+(
+ LDL_int n, /* L is n-by-n, where n >= 0 */
+ double X [ ], /* size n. right-hand-side on input, soln. on output */
+ LDL_int Lp [ ], /* input of size n+1, not modified */
+ LDL_int Li [ ], /* input of size lnz=Lp[n], not modified */
+ double Lx [ ] /* input of size lnz=Lp[n], not modified */
+)
+{
+ LDL_int j, p, p2 ;
+ for (j = 0 ; j < n ; j++)
+ {
+ p2 = Lp [j+1] ;
+ for (p = Lp [j] ; p < p2 ; p++)
+ {
+ X [Li [p]] -= Lx [p] * X [j] ;
+ }
+ }
+}
+
+
+/* ========================================================================== */
+/* === ldl_dsolve: solve Dx=b ============================================== */
+/* ========================================================================== */
+
+void LDL_dsolve
+(
+ LDL_int n, /* D is n-by-n, where n >= 0 */
+ double X [ ], /* size n. right-hand-side on input, soln. on output */
+ double D [ ] /* input of size n, not modified */
+)
+{
+ LDL_int j ;
+ for (j = 0 ; j < n ; j++)
+ {
+ X [j] /= D [j] ;
+ }
+}
+
+
+/* ========================================================================== */
+/* === ldl_ltsolve: solve L'x=b ============================================ */
+/* ========================================================================== */
+
+void LDL_ltsolve
+(
+ LDL_int n, /* L is n-by-n, where n >= 0 */
+ double X [ ], /* size n. right-hand-side on input, soln. on output */
+ LDL_int Lp [ ], /* input of size n+1, not modified */
+ LDL_int Li [ ], /* input of size lnz=Lp[n], not modified */
+ double Lx [ ] /* input of size lnz=Lp[n], not modified */
+)
+{
+ int j, p, p2 ;
+ for (j = n-1 ; j >= 0 ; j--)
+ {
+ p2 = Lp [j+1] ;
+ for (p = Lp [j] ; p < p2 ; p++)
+ {
+ X [j] -= Lx [p] * X [Li [p]] ;
+ }
+ }
+}
+
+
+/* ========================================================================== */
+/* === ldl_perm: permute a vector, x=Pb ===================================== */
+/* ========================================================================== */
+
+void LDL_perm
+(
+ LDL_int n, /* size of X, B, and P */
+ double X [ ], /* output of size n. */
+ double B [ ], /* input of size n. */
+ LDL_int P [ ] /* input permutation array of size n. */
+)
+{
+ LDL_int j ;
+ for (j = 0 ; j < n ; j++)
+ {
+ X [j] = B [P [j]] ;
+ }
+}
+
+
+/* ========================================================================== */
+/* === ldl_permt: permute a vector, x=P'b =================================== */
+/* ========================================================================== */
+
+void LDL_permt
+(
+ LDL_int n, /* size of X, B, and P */
+ double X [ ], /* output of size n. */
+ double B [ ], /* input of size n. */
+ LDL_int P [ ] /* input permutation array of size n. */
+)
+{
+ LDL_int j ;
+ for (j = 0 ; j < n ; j++)
+ {
+ X [P [j]] = B [j] ;
+ }
+}
+
+
+/* ========================================================================== */
+/* === ldl_valid_perm: check if a permutation vector is valid =============== */
+/* ========================================================================== */
+
+LDL_int LDL_valid_perm /* returns 1 if valid, 0 otherwise */
+(
+ LDL_int n,
+ LDL_int P [ ], /* input of size n, a permutation of 0:n-1 */
+ LDL_int Flag [ ] /* workspace of size n */
+)
+{
+ LDL_int j, k ;
+ if (n < 0 || !Flag)
+ {
+ return (0) ; /* n must be >= 0, and Flag must be present */
+ }
+ if (!P)
+ {
+ return (1) ; /* If NULL, P is assumed to be the identity perm. */
+ }
+ for (j = 0 ; j < n ; j++)
+ {
+ Flag [j] = 0 ; /* clear the Flag array */
+ }
+ for (k = 0 ; k < n ; k++)
+ {
+ j = P [k] ;
+ if (j < 0 || j >= n || Flag [j] != 0)
+ {
+ return (0) ; /* P is not valid */
+ }
+ Flag [j] = 1 ;
+ }
+ return (1) ; /* P is valid */
+}
+
+
+/* ========================================================================== */
+/* === ldl_valid_matrix: check if a sparse matrix is valid ================== */
+/* ========================================================================== */
+
+/* This routine checks to see if a sparse matrix A is valid for input to
+ * ldl_symbolic and ldl_numeric. It returns 1 if the matrix is valid, 0
+ * otherwise. A is in sparse column form. The numerical values in column j
+ * are stored in Ax [Ap [j] ... Ap [j+1]-1], with row indices in
+ * Ai [Ap [j] ... Ap [j+1]-1]. The Ax array is not checked.
+ */
+
+LDL_int LDL_valid_matrix
+(
+ LDL_int n,
+ LDL_int Ap [ ],
+ LDL_int Ai [ ]
+)
+{
+ LDL_int j, p ;
+ if (n < 0 || !Ap || !Ai || Ap [0] != 0)
+ {
+ return (0) ; /* n must be >= 0, and Ap and Ai must be present */
+ }
+ for (j = 0 ; j < n ; j++)
+ {
+ if (Ap [j] > Ap [j+1])
+ {
+ return (0) ; /* Ap must be monotonically nondecreasing */
+ }
+ }
+ for (p = 0 ; p < Ap [n] ; p++)
+ {
+ if (Ai [p] < 0 || Ai [p] >= n)
+ {
+ return (0) ; /* row indices must be in the range 0 to n-1 */
+ }
+ }
+ return (1) ; /* matrix is valid */
+}
=====================================
extlibs/csparse/ldl.h
=====================================
@@ -0,0 +1,111 @@
+/* ========================================================================== */
+/* === ldl.h: include file for the LDL package ============================= */
+/* ========================================================================== */
+
+/* LDL Copyright (c) Timothy A Davis,
+ * University of Florida. All Rights Reserved. See README for the License.
+ */
+
+#include "UFconfig.h"
+
+#ifdef LDL_LONG
+#define LDL_int UF_long
+#define LDL_ID UF_long_id
+
+#define LDL_symbolic ldl_l_symbolic
+#define LDL_numeric ldl_l_numeric
+#define LDL_lsolve ldl_l_lsolve
+#define LDL_dsolve ldl_l_dsolve
+#define LDL_ltsolve ldl_l_ltsolve
+#define LDL_perm ldl_l_perm
+#define LDL_permt ldl_l_permt
+#define LDL_valid_perm ldl_l_valid_perm
+#define LDL_valid_matrix ldl_l_valid_matrix
+
+#else
+#define LDL_int int
+#define LDL_ID "%d"
+
+#define LDL_symbolic ldl_symbolic
+#define LDL_numeric ldl_numeric
+#define LDL_lsolve ldl_lsolve
+#define LDL_dsolve ldl_dsolve
+#define LDL_ltsolve ldl_ltsolve
+#define LDL_perm ldl_perm
+#define LDL_permt ldl_permt
+#define LDL_valid_perm ldl_valid_perm
+#define LDL_valid_matrix ldl_valid_matrix
+
+#endif
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+/* ========================================================================== */
+/* === int version ========================================================== */
+/* ========================================================================== */
+
+void ldl_symbolic (int n, int Ap [ ], int Ai [ ], int Lp [ ],
+ int Parent [ ], int Lnz [ ], int Flag [ ], int P [ ], int Pinv [ ]) ;
+
+int ldl_numeric (int n, int Ap [ ], int Ai [ ], double Ax [ ],
+ int Lp [ ], int Parent [ ], int Lnz [ ], int Li [ ], double Lx [ ],
+ double D [ ], double Y [ ], int Pattern [ ], int Flag [ ],
+ int P [ ], int Pinv [ ]) ;
+
+void ldl_lsolve (int n, double X [ ], int Lp [ ], int Li [ ],
+ double Lx [ ]) ;
+
+void ldl_dsolve (int n, double X [ ], double D [ ]) ;
+
+void ldl_ltsolve (int n, double X [ ], int Lp [ ], int Li [ ],
+ double Lx [ ]) ;
+
+void ldl_perm (int n, double X [ ], double B [ ], int P [ ]) ;
+void ldl_permt (int n, double X [ ], double B [ ], int P [ ]) ;
+
+int ldl_valid_perm (int n, int P [ ], int Flag [ ]) ;
+int ldl_valid_matrix ( int n, int Ap [ ], int Ai [ ]) ;
+
+/* ========================================================================== */
+/* === long version ========================================================= */
+/* ========================================================================== */
+
+void ldl_l_symbolic (UF_long n, UF_long Ap [ ], UF_long Ai [ ], UF_long Lp [ ],
+ UF_long Parent [ ], UF_long Lnz [ ], UF_long Flag [ ], UF_long P [ ],
+ UF_long Pinv [ ]) ;
+
+UF_long ldl_l_numeric (UF_long n, UF_long Ap [ ], UF_long Ai [ ], double Ax [ ],
+ UF_long Lp [ ], UF_long Parent [ ], UF_long Lnz [ ], UF_long Li [ ],
+ double Lx [ ], double D [ ], double Y [ ], UF_long Pattern [ ],
+ UF_long Flag [ ], UF_long P [ ], UF_long Pinv [ ]) ;
+
+void ldl_l_lsolve (UF_long n, double X [ ], UF_long Lp [ ], UF_long Li [ ],
+ double Lx [ ]) ;
+
+void ldl_l_dsolve (UF_long n, double X [ ], double D [ ]) ;
+
+void ldl_l_ltsolve (UF_long n, double X [ ], UF_long Lp [ ], UF_long Li [ ],
+ double Lx [ ]) ;
+
+void ldl_l_perm (UF_long n, double X [ ], double B [ ], UF_long P [ ]) ;
+void ldl_l_permt (UF_long n, double X [ ], double B [ ], UF_long P [ ]) ;
+
+UF_long ldl_l_valid_perm (UF_long n, UF_long P [ ], UF_long Flag [ ]) ;
+UF_long ldl_l_valid_matrix ( UF_long n, UF_long Ap [ ], UF_long Ai [ ]) ;
+
+/* ========================================================================== */
+/* === LDL version ========================================================== */
+/* ========================================================================== */
+
+#define LDL_DATE "Nov 1, 2007"
+#define LDL_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
+#define LDL_MAIN_VERSION 2
+#define LDL_SUB_VERSION 0
+#define LDL_SUBSUB_VERSION 1
+#define LDL_VERSION LDL_VERSION_CODE(LDL_MAIN_VERSION,LDL_SUB_VERSION)
+
+#if defined(__cplusplus)
+} //extern "C"
+#endif
View it on GitLab: https://salsa.debian.org/med-team/sofa-framework/compare/e6e310e4e88b43cd9c764f54afa64e3bf18fcc0d...bbf9b000437bde7f888bea6c689fdb783d88747e
--
View it on GitLab: https://salsa.debian.org/med-team/sofa-framework/compare/e6e310e4e88b43cd9c764f54afa64e3bf18fcc0d...bbf9b000437bde7f888bea6c689fdb783d88747e
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20190919/6f0f5084/attachment-0001.html>
More information about the debian-med-commit
mailing list