diff --git a/Sofa/Component/Compat/src/SofaSparseSolver/SparseCholeskySolver.h b/Sofa/Component/Compat/src/SofaSparseSolver/SparseCholeskySolver.h index c3829690869..e08378ac130 100644 --- a/Sofa/Component/Compat/src/SofaSparseSolver/SparseCholeskySolver.h +++ b/Sofa/Component/Compat/src/SofaSparseSolver/SparseCholeskySolver.h @@ -23,4 +23,6 @@ #include -SOFA_DISABLED_HEADER("v22.06", "v23.06", "sofa/component/linearsolver/direct/SparseCholeskySolver.h") +SOFA_PRAGMA_ERROR( \ + "This header has been DISABLED since v23.12. " \ + "To fix this error you must use the CSparseSolvers plugins. " ) diff --git a/Sofa/Component/Compat/src/SofaSparseSolver/SparseLUSolver.h b/Sofa/Component/Compat/src/SofaSparseSolver/SparseLUSolver.h index 406f08044d2..4e21871483e 100644 --- a/Sofa/Component/Compat/src/SofaSparseSolver/SparseLUSolver.h +++ b/Sofa/Component/Compat/src/SofaSparseSolver/SparseLUSolver.h @@ -23,4 +23,6 @@ #include -SOFA_DISABLED_HEADER("v22.06", "v23.06", "sofa/component/linearsolver/direct/SparseLUSolver.h") +SOFA_PRAGMA_ERROR( \ + "This header has been DISABLED since v23.12. " \ + "To fix this error you must use the CSparseSolvers plugins. " ) diff --git a/Sofa/Component/Compat/src/SofaSparseSolver/SparseLUSolver.inl b/Sofa/Component/Compat/src/SofaSparseSolver/SparseLUSolver.inl index 4610b2cc02d..a50573db9b0 100644 --- a/Sofa/Component/Compat/src/SofaSparseSolver/SparseLUSolver.inl +++ b/Sofa/Component/Compat/src/SofaSparseSolver/SparseLUSolver.inl @@ -23,4 +23,6 @@ #include -SOFA_DISABLED_HEADER("v22.06", "v23.06", "sofa/component/linearsolver/direct/SparseLUSolver.inl") +SOFA_PRAGMA_ERROR( \ + "This header has been DISABLED since v23.12. " \ + "To fix this error you must use the CSparseSolvers plugins. " ) diff --git a/Sofa/Component/LinearSolver/Direct/CMakeLists.txt b/Sofa/Component/LinearSolver/Direct/CMakeLists.txt index 898a186778f..51fa093de06 100644 --- a/Sofa/Component/LinearSolver/Direct/CMakeLists.txt +++ b/Sofa/Component/LinearSolver/Direct/CMakeLists.txt @@ -50,14 +50,10 @@ set(SOURCE_FILES ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/MatrixLinearSystem[BTDMatrix].cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/PrecomputedLinearSolver.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/SVDLinearSolver.cpp - ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/SparseCholeskySolver.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/SparseCommon.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/SparseLDLSolver.cpp - ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/SparseLUSolver.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/TypedMatrixLinearSystem[BTDMatrix].cpp ) -add_subdirectory(extlibs/csparse) -sofa_set_01(SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE VALUE TRUE) sofa_find_package(metis QUIET) # Unix users can have an installed version of metis if(NOT metis_FOUND) @@ -75,7 +71,7 @@ sofa_find_package(Sofa.Component.LinearSolver.Iterative REQUIRED) # Only for Mat add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES} ${WRAPPER_FILES}) target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Simulation.Core Sofa.Component.LinearSolver.Iterative) -target_link_libraries(${PROJECT_NAME} PUBLIC metis csparse) +target_link_libraries(${PROJECT_NAME} PUBLIC metis) target_link_libraries(${PROJECT_NAME} PUBLIC Threads::Threads) sofa_create_package_with_targets( diff --git a/Sofa/Component/LinearSolver/Direct/Sofa.Component.LinearSolver.DirectConfig.cmake.in b/Sofa/Component/LinearSolver/Direct/Sofa.Component.LinearSolver.DirectConfig.cmake.in index 0b4b584fb08..4025df0d43e 100644 --- a/Sofa/Component/LinearSolver/Direct/Sofa.Component.LinearSolver.DirectConfig.cmake.in +++ b/Sofa/Component/LinearSolver/Direct/Sofa.Component.LinearSolver.DirectConfig.cmake.in @@ -8,14 +8,10 @@ find_package(Sofa.Component.LinearSolver.Iterative QUIET REQUIRED) find_package(Threads QUIET REQUIRED) set(SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_METIS @SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_METIS@) -set(SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE @SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE@) if(SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_METIS) find_package(Metis QUIET REQUIRED HINTS "${CMAKE_CURRENT_LIST_DIR}/..") endif() -if(SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE) - find_package(CSparse QUIET REQUIRED HINTS "${CMAKE_CURRENT_LIST_DIR}/..") -endif() if(NOT TARGET @PROJECT_NAME@) include("${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@Targets.cmake") diff --git a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/CMakeLists.txt b/Sofa/Component/LinearSolver/Direct/extlibs/csparse/CMakeLists.txt deleted file mode 100644 index 266aa286dd4..00000000000 --- a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/CMakeLists.txt +++ /dev/null @@ -1,32 +0,0 @@ -cmake_minimum_required(VERSION 3.12) -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} SYSTEM PUBLIC "$") -target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC "$") - -if(WIN32) - # remove warnings about deprecation (CRT,etc) - target_compile_options(${PROJECT_NAME} PRIVATE "/wd4996") -endif() - -sofa_create_package_with_targets( - PACKAGE_NAME CSparse - PACKAGE_VERSION ${PROJECT_VERSION} - TARGETS ${PROJECT_NAME} AUTO_SET_TARGET_PROPERTIES - INCLUDE_INSTALL_DIR "extlibs/CSparse" - ) diff --git a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/COPYING.txt b/Sofa/Component/LinearSolver/Direct/extlibs/csparse/COPYING.txt deleted file mode 100644 index 49ba8beccea..00000000000 --- a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/COPYING.txt +++ /dev/null @@ -1,3 +0,0 @@ -CSPARSE: a Concise Sparse matrix package. -Copyright (c) 2006, Timothy A. Davis. -http://www.cise.ufl.edu/research/sparse/CSparse diff --git a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/CSparseConfig.cmake.in b/Sofa/Component/LinearSolver/Direct/extlibs/csparse/CSparseConfig.cmake.in deleted file mode 100644 index af532008d06..00000000000 --- a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/CSparseConfig.cmake.in +++ /dev/null @@ -1,12 +0,0 @@ -# CMake package configuration file for the csparse library - -@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@) diff --git a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/UFconfig.h b/Sofa/Component/LinearSolver/Direct/extlibs/csparse/UFconfig.h deleted file mode 100644 index 54208d58b79..00000000000 --- a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/UFconfig.h +++ /dev/null @@ -1,118 +0,0 @@ -/* ========================================================================== */ -/* === 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 - -/* ========================================================================== */ -/* === 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 diff --git a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/csparse.c b/Sofa/Component/LinearSolver/Direct/extlibs/csparse/csparse.c deleted file mode 100644 index 14c7b342622..00000000000 --- a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/csparse.c +++ /dev/null @@ -1,2122 +0,0 @@ -# include -# include -# include -# include - -# 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) ; -} diff --git a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/csparse.h b/Sofa/Component/LinearSolver/Direct/extlibs/csparse/csparse.h deleted file mode 100644 index 1c0b8e514a3..00000000000 --- a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/csparse.h +++ /dev/null @@ -1,144 +0,0 @@ -#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 diff --git a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/ldl.c b/Sofa/Component/LinearSolver/Direct/extlibs/csparse/ldl.c deleted file mode 100644 index a9b35c846ef..00000000000 --- a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/ldl.c +++ /dev/null @@ -1,507 +0,0 @@ -/* ========================================================================== */ -/* === 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 */ -} diff --git a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/ldl.h b/Sofa/Component/LinearSolver/Direct/extlibs/csparse/ldl.h deleted file mode 100644 index 64d3dbb2d2e..00000000000 --- a/Sofa/Component/LinearSolver/Direct/extlibs/csparse/ldl.h +++ /dev/null @@ -1,111 +0,0 @@ -/* ========================================================================== */ -/* === 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 diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenDirectSparseSolver.inl b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenDirectSparseSolver.inl index 6795096feaf..8416631c581 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenDirectSparseSolver.inl +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenDirectSparseSolver.inl @@ -66,11 +66,10 @@ void EigenDirectSparseSolver Mfiltered.compress(); } - if (!m_map) - { - m_map = std::make_unique(Mfiltered.rows(), Mfiltered.cols(), Mfiltered.getColsValue().size(), - (typename EigenSparseMatrixMap::StorageIndex*)Mfiltered.rowBegin.data(), (typename EigenSparseMatrixMap::StorageIndex*)Mfiltered.colsIndex.data(), Mfiltered.colsValue.data()); - } + m_map = std::make_unique(Mfiltered.rows(), Mfiltered.cols(), Mfiltered.getColsValue().size(), + (typename EigenSparseMatrixMap::StorageIndex*)Mfiltered.rowBegin.data(), + (typename EigenSparseMatrixMap::StorageIndex*)Mfiltered.colsIndex.data(), + Mfiltered.colsValue.data()); const bool analyzePattern = (MfilteredrowBegin != Mfiltered.rowBegin) || (MfilteredcolsIndex != Mfiltered.colsIndex); diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.h index 1b5f0363237..c813cee4291 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.h +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.h @@ -92,9 +92,7 @@ class PrecomputedLinearSolver : public sofa::component::linearsolver::MatrixLine void invert(TMatrix& M) override; void setSystemMBKMatrix(const core::MechanicalParams* mparams) override; void loadMatrix(TMatrix& M); -#if SOFASPARSESOLVER_HAVE_CSPARSE - void loadMatrixWithCSparse(TMatrix& M); -#endif + void loadMatrixWithCholeskyDecomposition(TMatrix& M); bool addJMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) override; /// Returns the sofa template name. By default the name of the c++ class signature is exposed... diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.inl b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.inl index 0eb64a99dfb..03b9d5538b6 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.inl +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.inl @@ -29,6 +29,8 @@ #include #include #include +#include +#include #include #include #include @@ -37,10 +39,6 @@ #include -#if SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE && !defined(SOFA_FLOAT) -#include -#endif - #include namespace sofa::component::linearsolver::direct @@ -89,13 +87,8 @@ void PrecomputedLinearSolver::loadMatrix(TMatrix& M) ss << this->getContext()->getName() << "-" << systemSize << "-" << dt << ".comp"; if(! use_file.getValue() || ! internalData.readFile(ss.str().c_str(),systemSize) ) { -#if SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE && !defined(SOFA_FLOAT) - loadMatrixWithCSparse(M); + loadMatrixWithCholeskyDecomposition(M); if (use_file.getValue()) internalData.writeFile(ss.str().c_str(),systemSize); -#else - SOFA_UNUSED(M); - msg_error()<< "CSPARSE support is required to invert the matrix"; -#endif } for (unsigned int j=0; j::loadMatrix(TMatrix& M) } } -#if SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE && !defined(SOFA_FLOAT) template -void PrecomputedLinearSolver::loadMatrixWithCSparse(TMatrix& M) +void PrecomputedLinearSolver::loadMatrixWithCholeskyDecomposition(TMatrix& M) { using namespace sofa::linearalgebra; msg_info() << "Compute the initial invert matrix with CS_PARSE" ; @@ -123,7 +115,7 @@ void PrecomputedLinearSolver::loadMatrixWithCSparse(TMatrix& M) matSolv.resize(systemSize,systemSize); r.resize(systemSize); b.resize(systemSize); - SparseCholeskySolver, FullVector > solver; + EigenSimplicialLLT solver; for (unsigned int j=0; j::loadMatrixWithCSparse(TMatrix& M) msg_info() << "Precomputing constraint correction : " << std::fixed << 100.0f << " % " << '\xd'; } -#endif // SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE && !defined(SOFA_FLOAT) template void PrecomputedLinearSolver::invert(TMatrix& /*M*/) {} diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.cpp b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.cpp deleted file mode 100644 index d283bea94ff..00000000000 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.cpp +++ /dev/null @@ -1,44 +0,0 @@ -/****************************************************************************** -* SOFA, Simulation Open-Framework Architecture * -* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * -* * -* This program is free software; you can redistribute it and/or modify it * -* under the terms of the GNU Lesser General Public License as published by * -* the Free Software Foundation; either version 2.1 of the License, or (at * -* your option) any later version. * -* * -* This program is distributed in the hope that it will be useful, but WITHOUT * -* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * -* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * -* for more details. * -* * -* You should have received a copy of the GNU Lesser General Public License * -* along with this program. If not, see . * -******************************************************************************* -* Authors: The SOFA Team and external contributors (see Authors.txt) * -* * -* Contact information: contact@sofa-framework.org * -******************************************************************************/ -#define SOFA_COMPONENT_LINEARSOLVER_SPARSECHOLESKYSOLVER_CPP -#include -#include - -namespace sofa::component::linearsolver::direct -{ - -using namespace sofa::linearalgebra; - -#ifdef SOFA_FLOAT -SOFA_PRAGMA_WARNING("SparseCholeskySolver does not support float as scalar.") -#else // SOFA_DOUBLE -int SparseCholeskySolverClass = - core::RegisterObject( - "Direct linear solver based on Sparse Cholesky factorization, implemented with the " - "CSPARSE library") - .add, FullVector > >(); - -template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API - SparseCholeskySolver, FullVector >; -#endif // SOFA_FLOAT - -} // namespace sofa::component::linearsolver::direct diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.h index 9fd72c37db4..e4bf3e0b8cf 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.h +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.h @@ -20,61 +20,8 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #pragma once +#include -#include - -#include -#include -#include -#include -#include - -namespace sofa::component::linearsolver::direct -{ - -// Direct linear solver based on Sparse Cholesky factorization, implemented with the CSPARSE library -template -class SparseCholeskySolver : public sofa::component::linearsolver::MatrixLinearSolver -{ -public: - SOFA_CLASS(SOFA_TEMPLATE2(SparseCholeskySolver,TMatrix,TVector),SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver,TMatrix,TVector)); - - typedef TMatrix Matrix; - typedef TVector Vector; - - SparseCholeskySolver(); - ~SparseCholeskySolver() override; - - void solve (Matrix& M, Vector& x, Vector& b) override; - void invert(Matrix& M) override; - - void parse(core::objectmodel::BaseObjectDescription *arg) override; - -protected: - - SOFA_ATTRIBUTE_DEPRECATED__SOLVER_DIRECT_VERBOSEDATA() - Data f_verbose; ///< Dump system state at each iteration - - cs A; - cs* permuted_A; - css *S; - csn *N; - int * A_i; ///< row indices, size nzmax - int * A_p; ///< column pointers (size n+1) or col indices (size nzmax) - type::vector Previous_colptr,Previous_rowind; ///< shape of the matrix at the previous step - type::vector perm,iperm; ///< fill reducing permutation - type::vector A_x,z_tmp,r_tmp,tmp; - bool notSameShape; - - Data d_typePermutation; - - void suiteSparseFactorization(bool applyPermutation); - - css* symbolic_Chol(cs *A); -}; - -#if !defined(SOFA_COMPONENT_LINEARSOLVER_SPARSECHOLESKYSOLVER_CPP) -extern template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API SparseCholeskySolver< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector >; -#endif - -} // namespace sofa::component::linearsolver::direct +SOFA_PRAGMA_ERROR( \ + "This header has been DISABLED since v23.12. " \ + "To fix this error you must use the CSparseSolvers plugins. " ) diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.inl b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.inl index 70190cd92c5..e4bf3e0b8cf 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.inl +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCholeskySolver.inl @@ -20,206 +20,8 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #pragma once +#include -#include -#include - -namespace sofa::component::linearsolver::direct -{ - -template -SparseCholeskySolver::SparseCholeskySolver() - : S(nullptr), N(nullptr) - , d_typePermutation(initData(&d_typePermutation, {"None", "SuiteSparse", "METIS"}, - "permutation", "Type of fill reducing permutation")) -{} - -template -SparseCholeskySolver::~SparseCholeskySolver() -{ - if (S) cs_sfree (S); - if (N) cs_nfree (N); -} - -template -void SparseCholeskySolver::solve (Matrix& /*M*/, Vector& x, Vector& b) -{ - const int n = A.n; - - SCOPED_TIMER_VARNAME(solveTimer, "solve"); - - switch( d_typePermutation.getValue().getSelectedId() ) - { - case 0://None->identity - case 1://SuiteSparse - if(N) - { - cs_ipvec (n, S->Pinv, (double*)b.ptr() , tmp.data() ); //x = P*b , permutation on rows - cs_lsolve (N->L, tmp.data() ); //x = L\x - cs_ltsolve (N->L, tmp.data() ); //x = L'\x/ - cs_pvec (n, S->Pinv, tmp.data() , (double*)x.ptr() ); //x = P'*x , permutation on columns - } - else - { - msg_error() << "Cannot solve system due to invalid factorization"; - } - break; - - case 2://METIS - if(N) - { - cs_ipvec (n, perm.data(), (double*)b.ptr() , tmp.data() ); //x = P*b , permutation on rows - cs_lsolve (N->L, tmp.data() ); //x = L\x - cs_ltsolve (N->L, tmp.data() ); //x = L'\x/ - cs_pvec (n, perm.data() , tmp.data() , (double*)x.ptr() ); //x = P'*x , permutation on columns - } - else - { - msg_error() << "Cannot solve system due to invalid factorization"; - } - break; - - default: - break; - - } - -} - -template -void SparseCholeskySolver::invert(Matrix& M) -{ - if (N) cs_nfree(N); - M.compress(); - - A.nzmax = M.getColsValue().size(); // maximum number of entries - A_p = (int *) &(M.getRowBegin()[0]); - A_i = (int *) &(M.getColsIndex()[0]); - A_x.resize(A.nzmax); - for (int i=0; iidentity - suiteSparseFactorization(false); - break; - - case 1:// SuiteSparse - suiteSparseFactorization(true); - break; - - case 2:// METIS - if( notSameShape ) - { - perm.resize(A.n); - iperm.resize(A.n); - - fillReducingPermutation( A , iperm.data(), perm.data() ); // compute the fill reducing permutation - } - - permuted_A = cs_permute( &A , perm.data() , iperm.data() , 1); - - if ( notSameShape ) - { - if (S) cs_sfree(S); - S = symbolic_Chol( permuted_A ); - } // symbolic analysis - - N = cs_chol (permuted_A, S) ; // numeric Cholesky factorization - assert(N); - - cs_free(permuted_A); - break; - } - } - - // store the shape of the matrix - if ( notSameShape ) - { - Previous_rowind.clear(); - Previous_colptr.resize(A.n +1); - for(int i=0 ; i -void SparseCholeskySolver::parse(core::objectmodel::BaseObjectDescription* arg) -{ - if (arg->getAttribute("verbose")) - { - msg_warning() << "Attribute 'verbose' has no use in this component. " - "To disable this warning, remove the attribute from the scene."; - } - - Inherit1::parse(arg); -} - -template -void SparseCholeskySolver::suiteSparseFactorization(bool applyPermutation) -{ - if( notSameShape ) - { - if (S) - { - cs_sfree(S); - } - const auto order = applyPermutation ? 0 : -1; - S = cs_schol (&A, order); - } - assert(S); - assert(S->cp); - assert(S->parent); - N = cs_chol (&A, S) ; // numeric Cholesky factorization - msg_error_when(!N) << "Matrix could not be factorized: possibly not positive-definite"; -} - -template -css* SparseCholeskySolver::symbolic_Chol(cs *A) -{ //based on cs_schol - int n, *c, *post; - cs *C ; - css *S ; - if (!A) return (NULL) ; // check inputs - n = A->n ; - S = (css*)cs_calloc (1, sizeof (css)) ; // allocate symbolic analysis - if (!S) return (NULL) ; // out of memory - 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 = (int*)cs_malloc (n+1, sizeof (int)) ; // find column pointers for L - S->unz = S->lnz = cs_cumsum (S->cp, c, n) ; - // we do not use the permutation of SuiteSparse - S->Q = nullptr ; // permutation on columns set to identity - S->Pinv = nullptr; // permutation on rows set to identity - cs_free (c) ; - return ((S->lnz >= 0) ? S : cs_sfree (S)) ; -} - -} // namespace sofa::component::linearsolver::direct +SOFA_PRAGMA_ERROR( \ + "This header has been DISABLED since v23.12. " \ + "To fix this error you must use the CSparseSolvers plugins. " ) diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCommon.cpp b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCommon.cpp index 91439016f57..af4764e0ac3 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCommon.cpp +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCommon.cpp @@ -19,10 +19,12 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ - - #include +extern "C" { +#include +} + namespace sofa::component::linearsolver::direct { void csrToAdj(int n, int * M_colptr, int * M_rowind, type::vector& adj, type::vector& xadj , type::vector& t_adj , type::vector& t_xadj, type::vector& tran_countvec) @@ -86,14 +88,13 @@ void csrToAdj(int n, int * M_colptr, int * M_rowind, type::vector& adj, typ } } - -void fillReducingPermutation(const cs &A,int * perm,int * invperm) +void fillReducingPermutation(int nbColumns, int *columns, int* rowIndices, + int * perm,int * invperm) { - int n = A.n; sofa::type::vector adj, xadj, t_adj, t_xadj, tran_countvec; - csrToAdj( A.n, A.p , A.i , adj, xadj, t_adj, t_xadj, tran_countvec ); - METIS_NodeND(&n, xadj.data(), adj.data(), nullptr, nullptr, perm, invperm); - + csrToAdj( nbColumns, columns , rowIndices , adj, xadj, t_adj, t_xadj, tran_countvec ); + METIS_NodeND(&nbColumns, xadj.data(), adj.data(), nullptr, nullptr, perm, invperm); } + } diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCommon.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCommon.h index a2334a40010..37cf3400d52 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCommon.h +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCommon.h @@ -24,11 +24,6 @@ #include #include -#include - -extern "C" { -#include -} namespace sofa::component::linearsolver::direct { @@ -46,7 +41,9 @@ SOFA_COMPONENT_LINEARSOLVER_DIRECT_API void csrToAdj(int n, int * M_colptr, int * M_rowind, type::vector& adj, type::vector& xadj, type::vector& t_adj, type::vector& t_xadj, type::vector& tran_countvec ); // compute the fill reducing permutation via METIS -void fillReducingPermutation(const cs &A,int * perm,int * invperm); +SOFA_COMPONENT_LINEARSOLVER_DIRECT_API +void fillReducingPermutation(int nbColumns, int *columns, int* rowIndices, + int * perm,int * invperm); // compare the shape of two matrix given in CSR format, return false if the matrices have the same shape and return true if their shapes are different inline bool compareMatrixShape(int s_M, int * M_colptr,int * M_rowind, int s_P, int * P_colptr,int * P_rowind) { diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLDLSolverImpl.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLDLSolverImpl.h index cf1d51c1f4a..2baf93b64ee 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLDLSolverImpl.h +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLDLSolverImpl.h @@ -26,7 +26,6 @@ #include #include #include -#include #include #include diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.cpp b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.cpp deleted file mode 100644 index dd4faff9814..00000000000 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.cpp +++ /dev/null @@ -1,46 +0,0 @@ -/****************************************************************************** -* SOFA, Simulation Open-Framework Architecture * -* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * -* * -* This program is free software; you can redistribute it and/or modify it * -* under the terms of the GNU Lesser General Public License as published by * -* the Free Software Foundation; either version 2.1 of the License, or (at * -* your option) any later version. * -* * -* This program is distributed in the hope that it will be useful, but WITHOUT * -* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * -* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * -* for more details. * -* * -* You should have received a copy of the GNU Lesser General Public License * -* along with this program. If not, see . * -******************************************************************************* -* Authors: The SOFA Team and external contributors (see Authors.txt) * -* * -* Contact information: contact@sofa-framework.org * -******************************************************************************/ -#define SOFA_COMPONENT_LINEARSOLVER_SPARSELUSOLVER_CPP -#include -#include -#include - -namespace sofa::component::linearsolver::direct -{ - -using namespace sofa::linearalgebra; - -// dont try to compile if floating point is float -#ifdef SOFA_FLOAT -SOFA_PRAGMA_WARNING("SparseLUSolver does not support float as scalar.") -#else // SOFA_DOUBLE -int SparseLUSolverClass = core::RegisterObject("Direct linear solver based on Sparse LU factorization, implemented with the CSPARSE library") - .add< SparseLUSolver< CompressedRowSparseMatrix, FullVector > >() - .add< SparseLUSolver< CompressedRowSparseMatrix >,FullVector > >() - ; - -template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API SparseLUSolver< CompressedRowSparseMatrix, FullVector >; -template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API SparseLUSolver< CompressedRowSparseMatrix >, FullVector >; -#endif // SOFA_FLOAT - -} // namespace sofa::component::linearsolver::direct - diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.h index b6ebba1a2ad..e4bf3e0b8cf 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.h +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.h @@ -20,89 +20,8 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #pragma once -#include +#include -#include -#include -#include -#include - -namespace sofa::component::linearsolver::direct -{ - -//defaut structure for a LU factorization -template -class SparseLUInvertData : public MatrixInvertData { -public : - - css *S; ///< store the permutations and the number of non null values by rows and by lines of the LU factorization - csn *N; ///< store the partial pivot and the LU factorization - cs A; - cs* permuted_A; - type::vector perm,iperm; ///< fill reducing permutation - type::vector Previous_colptr,Previous_rowind; ///< shape of the matrix at the previous step - type::vector A_i, A_p; - type::vector A_x; - Real * tmp; - bool notSameShape; - SparseLUInvertData() - { - S=nullptr; N=nullptr; tmp=nullptr; - } - - ~SparseLUInvertData() - { - if (S) cs_sfree (S); - if (N) cs_nfree (N); - if (tmp) cs_free (tmp); - } -}; - -// Direct linear solver based on Sparse LU factorization, implemented with the CSPARSE library -template -class SparseLUSolver : public sofa::component::linearsolver::MatrixLinearSolver -{ -public: - SOFA_CLASS(SOFA_TEMPLATE3(SparseLUSolver,TMatrix,TVector,TThreadManager),SOFA_TEMPLATE3(sofa::component::linearsolver::MatrixLinearSolver,TMatrix,TVector,TThreadManager)); - - typedef TMatrix Matrix; - typedef TVector Vector; - typedef typename Matrix::Real Real; - - typedef sofa::component::linearsolver::MatrixLinearSolver Inherit; - - SOFA_ATTRIBUTE_DEPRECATED__SOLVER_DIRECT_VERBOSEDATA() - Data f_verbose; ///< Dump system state at each iteration - - Data f_tol; ///< tolerance of factorization - - void solve (Matrix& M, Vector& x, Vector& b) override; - void invert(Matrix& M) override; - - SparseLUSolver(); - - bool supportNonSymmetricSystem() const override { return true; } - - void parse(core::objectmodel::BaseObjectDescription *arg) override; - -protected : - - Data d_typePermutation; - Data d_L_nnz; ///< Number of non-zero values in the lower triangular matrix of the factorization. The lower, the faster the system is solved. - - css* symbolic_LU(cs *A); - - MatrixInvertData * createInvertData() override { - return new SparseLUInvertData(); - } - - std::unique_ptr > Mfiltered; - -}; - - -#if !defined(SOFA_COMPONENT_LINEARSOLVER_SPARSELUSOLVER_CPP) -extern template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API SparseLUSolver< sofa::linearalgebra::CompressedRowSparseMatrix< SReal>, sofa::linearalgebra::FullVector >; -#endif - -} // namespace sofa::component::linearsolver::direct +SOFA_PRAGMA_ERROR( \ + "This header has been DISABLED since v23.12. " \ + "To fix this error you must use the CSparseSolvers plugins. " ) diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.inl b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.inl index 7c0339024f2..e4bf3e0b8cf 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.inl +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLUSolver.inl @@ -20,234 +20,8 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #pragma once +#include -#include -#include - -namespace sofa::component::linearsolver::direct -{ - -using namespace sofa::defaulttype; -using namespace sofa::core::behavior; -using namespace sofa::simulation; -using namespace sofa::core::objectmodel; -using sofa::helper::system::thread::CTime; -using sofa::helper::system::thread::ctime_t; -using std::cerr; -using std::endl; - -template -SparseLUSolver::SparseLUSolver() - : f_tol( initData(&f_tol,0.001,"tolerance","tolerance of factorization") ) - , d_typePermutation(initData(&d_typePermutation , "permutation", "Type of fill reducing permutation")) - , d_L_nnz(initData(&d_L_nnz, 0, "L_nnz", "Number of non-zero values in the lower triangular matrix of the factorization. The lower, the faster the system is solved.", true, true)) -{ - sofa::helper::OptionsGroup d_typePermutationOptions{"None", "SuiteSparse", "METIS"}; - d_typePermutationOptions.setSelectedItem(0); // default None - d_typePermutation.setValue(d_typePermutationOptions); -} - -template -void SparseLUSolver::parse( - core::objectmodel::BaseObjectDescription* arg) -{ - if (arg->getAttribute("verbose")) - { - msg_warning() << "Attribute 'verbose' has no use in this component. " - "To disable this warning, remove the attribute from the scene."; - } - - Inherit::parse(arg); -} - - -template -void SparseLUSolver::solve (Matrix& M, Vector& x, Vector& b) -{ - SparseLUInvertData * invertData = (SparseLUInvertData*) this->getMatrixInvertData(&M); - const int n = invertData->A.n; - - { - SCOPED_TIMER_VARNAME(solveTimer, "solve"); - switch( d_typePermutation.getValue().getSelectedId() ) - { - - case 0://None->Identity - - { - cs_ipvec (n, nullptr, b.ptr(), x.ptr()) ; // copy - cs_lsolve (invertData->N->L, x.ptr() ) ; // x = L\x - cs_usolve (invertData->N->U, x.ptr() ) ; // x = U\x - break; - } - - case 1://SuiteSparse - { - cs_ipvec (n, invertData->N->Pinv, b.ptr(), invertData->tmp) ; // x = P*b , partial pivot - cs_lsolve (invertData->N->L, invertData->tmp) ; // x = L\x - cs_usolve (invertData->N->U, invertData->tmp) ; // x = U\x - cs_ipvec (n, invertData->S->Q, invertData->tmp, x.ptr()) ; // x = Q*x fill reducing permutation on columns only - break; - } - - case 2://METIS - { - // no partial pivoting - cs_pvec (n, invertData->perm.data() , b.ptr(), invertData->tmp) ; // x = P*b permutation on rows - cs_lsolve (invertData->N->L, invertData->tmp) ; // x = L\x - cs_usolve (invertData->N->U, invertData->tmp) ; // x = U\x - cs_pvec (n, invertData->iperm.data() , invertData->tmp , x.ptr()) ; // x = Q*x permutation on columns - break; - } - - default: - break; - } - } -} - -template -void SparseLUSolver::invert(Matrix& M) -{ - SparseLUInvertData * invertData = (SparseLUInvertData*) this->getMatrixInvertData(&M); - - sofa::linearalgebra::CompressedRowSparseMatrix* matrix; - - if constexpr (!std::is_same_v) - { - if (!Mfiltered) - { - Mfiltered = std::make_unique >(); - } - Mfiltered->copyNonZeros(M); - Mfiltered->compress(); - - matrix = Mfiltered.get(); - } - else - { - M.compress(); - matrix = &M; - } - - if (invertData->N) cs_nfree(invertData->N); - if (invertData->tmp) cs_free(invertData->tmp); - - //build A with M - invertData->A.nzmax = matrix->getColsValue().size(); // maximum number of entries - invertData->A.m = matrix->rowSize(); // number of rows - invertData->A.n = matrix->colSize(); // number of columns - invertData->A_p = matrix->getRowBegin(); - invertData->A.p = (int *) &(invertData->A_p[0]); // column pointers (size n+1) or col indices (size nzmax) - invertData->A_i = matrix->getColsIndex(); - invertData->A.i = (int *) &(invertData->A_i[0]); // row indices, size nzmax - invertData->A_x = matrix->getColsValue(); - invertData->A.x = (Real *) &(invertData->A_x[0]); // numerical values, size nzmax - invertData->A.nz = -1; // # of entries in triplet matrix, -1 for compressed-col - cs_dropzeros( &invertData->A ); - - invertData->notSameShape = compareMatrixShape(invertData->A.n , (int*) invertData->A_p.data() , (int*) invertData->A_i.data(), (invertData->Previous_colptr.size())-1 ,invertData->Previous_colptr.data() ,invertData->Previous_rowind.data() ); - - invertData->tmp = (Real *) cs_malloc (invertData->A.n, sizeof (Real)) ; - - { - SCOPED_TIMER_VARNAME(factorizationTimer, "factorization"); - switch( d_typePermutation.getValue().getSelectedId() ) - { - case 0://None->Identity - { - if (invertData->S) - { - cs_sfree(invertData->S); - } - invertData->S = symbolic_LU( &(invertData->A) ) ; // ordering and symbolic analysis - invertData->N = cs_lu (&invertData->A, invertData->S, f_tol.getValue()) ; // numeric LU factorization - break; - } - - case 1://SuiteSparse - { - if (invertData->notSameShape ) - { - if (invertData->S) - { - cs_sfree(invertData->S); - } - invertData->S = cs_sqr (&invertData->A, 1 , 0) ; // Suitsparse compute fill reducing permutation for LU decomposition - } - invertData->N = cs_lu (&invertData->A, invertData->S, f_tol.getValue()) ; // numeric LU factorization - break; - } - - case 2://Metis - { - if (invertData->notSameShape ) - { - invertData->perm.resize(invertData->A.n); - invertData->iperm.resize(invertData->A.n); - - fillReducingPermutation(invertData->A, invertData->perm.data(), invertData->iperm.data() ); //METIS compute the fill reducing permutation - } - - invertData->permuted_A = cs_permute(&(invertData->A), invertData->iperm.data(), invertData->perm.data(), 1); - - if (invertData->notSameShape ) - { - if (invertData->S) - { - cs_sfree(invertData->S); - } - invertData->S = symbolic_LU( invertData->permuted_A ); - } - - invertData->N = cs_lu ( invertData->permuted_A, invertData->S, f_tol.getValue()) ; // numeric LU factorization - - cs_free(invertData->permuted_A); // prevent memory leak - - break; - } - - default: - break; - } - } - - d_L_nnz.setValue(invertData->N->L->nzmax); - - // store the shape of the matrix - if ( invertData->notSameShape ) - { - invertData->Previous_rowind.clear(); - invertData->Previous_colptr.resize( (invertData->A.n) +1); - - for (int i=0 ; iA.n ; i++) - { - invertData->Previous_colptr[i+1] = invertData->A_p[i+1]; - - for ( int j = (int) invertData->A_p[i] ; j < (int)invertData->A_p[i+1] ; j++) - { - invertData->Previous_rowind.push_back( invertData->A_i[j]); - } - } - - } - -} - -template -css* SparseLUSolver::symbolic_LU(cs *A) -{// based on cs_sqr - - int n; - css *S ; - if (!A) return (NULL) ; /* check inputs */ - n = A->n ; - S = (css*)cs_calloc (1, sizeof (css)) ; /* allocate symbolic analysis */ - if (!S) return (NULL) ; /* out of memory */ - S->unz = 4*(A->p [n]) + n ; /* for LU factorization only, */ - S->lnz = S->unz ; /* guess nnz(L) and nnz(U) */ - S->Q = nullptr; // should have been the fill permutation computed by SuiteSparse, not used here - return S ; -} - -} // namespace sofa::component::linearsolver::direct +SOFA_PRAGMA_ERROR( \ + "This header has been DISABLED since v23.12. " \ + "To fix this error you must use the CSparseSolvers plugins. " ) diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.h b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.h index 9c71ecd0aa7..0ab268218b2 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.h +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.h @@ -176,9 +176,7 @@ protected : PrecomputedWarpPreconditionerInternalData internalData; void rotateConstraints(); -#if SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE && !defined(SOFA_FLOAT) - void loadMatrixWithCSparse(TMatrix& M); -#endif + void loadMatrixWithCholeskyDecomposition(TMatrix& M); void loadMatrixWithSolver(); template diff --git a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.inl b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.inl index 5155ed3691f..01627d1ace6 100644 --- a/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.inl +++ b/Sofa/Component/LinearSolver/Preconditioner/src/sofa/component/linearsolver/preconditioner/PrecomputedWarpPreconditioner.inl @@ -39,11 +39,8 @@ #include #include - -#if SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE && !defined(SOFA_FLOAT) -#include -#endif -#include +#include +#include #include @@ -184,13 +181,7 @@ void PrecomputedWarpPreconditioner::loadMatrix(TMatrix& M) msg_info() << "Precompute : " << fname << " compliance."; if (l_linearSolver.empty()) { -#if SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE && !defined(SOFA_FLOAT) - loadMatrixWithCSparse(M); -#else - msg_error() << "Link \"linearSolver\" is empty, but it is required to load matrix."; - this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); - return; -#endif // SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE && !defined(SOFA_FLOAT) + loadMatrixWithCholeskyDecomposition(M); } else { @@ -222,9 +213,8 @@ void PrecomputedWarpPreconditioner::loadMatrix(TMatrix& M) } } -#if SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE && !defined(SOFA_FLOAT) template -void PrecomputedWarpPreconditioner::loadMatrixWithCSparse(TMatrix& M) +void PrecomputedWarpPreconditioner::loadMatrixWithCholeskyDecomposition(TMatrix& M) { msg_info() << "Compute the initial invert matrix with CS_PARSE" ; @@ -237,7 +227,7 @@ void PrecomputedWarpPreconditioner::loadMatrixWithCSparse(TMatrix& M b.resize(systemSize); for (unsigned int j=0; j, FullVector > solver; + direct::EigenSimplicialLLT solver; msg_info() << "Precomputing constraint correction LU decomposition " ; solver.invert(M); @@ -277,7 +267,6 @@ void PrecomputedWarpPreconditioner::loadMatrixWithCSparse(TMatrix& M tmpStr << "Precomputing constraint correction : " << std::fixed << 100.0f << " %" ; msg_info() << tmpStr.str(); } -#endif // SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CSPARSE && !defined(SOFA_FLOAT) template void PrecomputedWarpPreconditioner::loadMatrixWithSolver() diff --git a/Sofa/framework/Core/src/sofa/core/topology/TopologyData.cpp b/Sofa/framework/Core/src/sofa/core/topology/TopologyData.cpp index ac1a0847a0c..433fed3d6b8 100644 --- a/Sofa/framework/Core/src/sofa/core/topology/TopologyData.cpp +++ b/Sofa/framework/Core/src/sofa/core/topology/TopologyData.cpp @@ -26,7 +26,7 @@ namespace sofa::core::topology { - + template class SOFA_CORE_API sofa::core::topology::TopologyData < core::topology::BaseMeshTopology::Point, type::vector >; } //namespace sofa::core::topology diff --git a/Sofa/framework/Helper/src/sofa/helper/ComponentChange.cpp b/Sofa/framework/Helper/src/sofa/helper/ComponentChange.cpp index b7958a4da16..a3cb1196bed 100644 --- a/Sofa/framework/Helper/src/sofa/helper/ComponentChange.cpp +++ b/Sofa/framework/Helper/src/sofa/helper/ComponentChange.cpp @@ -412,9 +412,7 @@ const std::map > uncreatableComponents // SofaSparseSolver was deprecated in #2717 { "FillReducingOrdering", Moved("v22.06", "SofaGeneralLinearSolver", "Sofa.Component.LinearSolver.Direct") }, { "PrecomputedLinearSolver", Moved("v22.06", "SofaGeneralLinearSolver", "Sofa.Component.LinearSolver.Direct") }, - { "SparseCholeskySolver", Moved("v22.06", "SofaSparseSolver", "Sofa.Component.LinearSolver.Direct") }, { "SparseLDLSolver", Moved("v22.06", "SofaSparseSolver", "Sofa.Component.LinearSolver.Direct") }, - { "SparseLUSolver", Moved("v22.06", "SofaSparseSolver", "Sofa.Component.LinearSolver.Direct") }, // SofaDenseSolver was deprecated in #2717 { "SVDLinearSolver", Moved("v22.06", "SofaDenseSolver", "Sofa.Component.LinearSolver.Direct") }, @@ -720,6 +718,10 @@ const std::map > uncreatableComponents // Removed in #4040, deprecated in #2777 { "MechanicalMatrixMapper", Removed("v23.06", "v23.12") }, { "MappingGeometricStiffnessForceField", Removed("v23.06", "v23.12") }, + + // Moved to CSparseSolvers + { "SparseCholeskySolver", Moved("v23.12", "Sofa.Component.LinearSolver.Direct", "CSparseSolvers") }, + { "SparseLUSolver", Moved("v23.12", "Sofa.Component.LinearSolver.Direct", "CSparseSolvers") }, }; } // namespace sofa::helper::lifecycle diff --git a/applications/plugins/CMakeLists.txt b/applications/plugins/CMakeLists.txt index 6d05efc1520..1e72b41f483 100644 --- a/applications/plugins/CMakeLists.txt +++ b/applications/plugins/CMakeLists.txt @@ -51,6 +51,7 @@ sofa_add_subdirectory(plugin SoftRobots SoftRobots EXTERNAL) sofa_add_subdirectory(plugin CollisionAlgorithm CollisionAlgorithm EXTERNAL) sofa_add_subdirectory(plugin ConstraintGeometry ConstraintGeometry EXTERNAL) sofa_add_subdirectory(plugin ShapeMatchingPlugin ShapeMatchingPlugin EXTERNAL) +sofa_add_subdirectory(plugin CSparseSolvers CSparseSolvers EXTERNAL) sofa_add_subdirectory(plugin PSL PSL EXTERNAL) diff --git a/applications/plugins/CSparseSolvers/ExternalProjectConfig.cmake.in b/applications/plugins/CSparseSolvers/ExternalProjectConfig.cmake.in new file mode 100644 index 00000000000..989dbe70379 --- /dev/null +++ b/applications/plugins/CSparseSolvers/ExternalProjectConfig.cmake.in @@ -0,0 +1,14 @@ +cmake_minimum_required(VERSION 3.11) + +include(ExternalProject) +ExternalProject_Add(CSparseSolvers + GIT_REPOSITORY https://github.com/sofa-framework/CSparseSolvers + GIT_TAG origin/master + SOURCE_DIR "${CMAKE_SOURCE_DIR}/applications/plugins/CSparseSolvers" + BINARY_DIR "" + CONFIGURE_COMMAND "" + BUILD_COMMAND "" + INSTALL_COMMAND "" + TEST_COMMAND "" + GIT_CONFIG "remote.origin.fetch=+refs/pull/*:refs/remotes/origin/pr/*" +) diff --git a/applications/plugins/SofaCUDA/examples/beam10x10x46-warp-preconditioner-CUDA.scn b/applications/plugins/SofaCUDA/examples/beam10x10x46-warp-preconditioner-CUDA.scn index 169c19466eb..50031644065 100644 --- a/applications/plugins/SofaCUDA/examples/beam10x10x46-warp-preconditioner-CUDA.scn +++ b/applications/plugins/SofaCUDA/examples/beam10x10x46-warp-preconditioner-CUDA.scn @@ -32,7 +32,6 @@ - diff --git a/applications/plugins/SofaMatrix/examples/ComplianceMatrixExporter.scn b/applications/plugins/SofaMatrix/examples/ComplianceMatrixExporter.scn index d3eba819d21..4949de49328 100644 --- a/applications/plugins/SofaMatrix/examples/ComplianceMatrixExporter.scn +++ b/applications/plugins/SofaMatrix/examples/ComplianceMatrixExporter.scn @@ -6,7 +6,7 @@ - + @@ -24,7 +24,7 @@ - + diff --git a/applications/plugins/SofaMatrix/examples/ComplianceMatrixImage.scn b/applications/plugins/SofaMatrix/examples/ComplianceMatrixImage.scn index 98e93a04704..d443ea93089 100644 --- a/applications/plugins/SofaMatrix/examples/ComplianceMatrixImage.scn +++ b/applications/plugins/SofaMatrix/examples/ComplianceMatrixImage.scn @@ -6,7 +6,7 @@ - + @@ -24,7 +24,7 @@ - + diff --git a/applications/plugins/SofaMatrix/examples/FillReducingOrdering.scn b/applications/plugins/SofaMatrix/examples/FillReducingOrdering.scn index f1b20911242..6e5f3c40528 100644 --- a/applications/plugins/SofaMatrix/examples/FillReducingOrdering.scn +++ b/applications/plugins/SofaMatrix/examples/FillReducingOrdering.scn @@ -10,7 +10,7 @@ The scene compares two simulations in which only the vertices order differs: - + @@ -33,7 +33,7 @@ The scene compares two simulations in which only the vertices order differs: - + @@ -49,7 +49,7 @@ The scene compares two simulations in which only the vertices order differs: - + diff --git a/examples/Component/Constraint/Lagrangian/InextensiblePendulum.scn b/examples/Component/Constraint/Lagrangian/InextensiblePendulum.scn index 939ff0fecea..ab6bdf1cb57 100644 --- a/examples/Component/Constraint/Lagrangian/InextensiblePendulum.scn +++ b/examples/Component/Constraint/Lagrangian/InextensiblePendulum.scn @@ -10,7 +10,7 @@ - + @@ -35,7 +35,7 @@ - + diff --git a/examples/Component/LinearSolver/Direct/FEMBAR_SparseCholeskySolver.scn b/examples/Component/LinearSolver/Direct/FEMBAR_SparseCholeskySolver.scn deleted file mode 100644 index a9e826a2956..00000000000 --- a/examples/Component/LinearSolver/Direct/FEMBAR_SparseCholeskySolver.scn +++ /dev/null @@ -1,8 +0,0 @@ - - - - - - - - diff --git a/examples/Component/LinearSolver/Direct/FEMBAR_SparseLUSolver.scn b/examples/Component/LinearSolver/Direct/FEMBAR_SparseLUSolver.scn deleted file mode 100644 index 241b8fb05d3..00000000000 --- a/examples/Component/LinearSolver/Direct/FEMBAR_SparseLUSolver.scn +++ /dev/null @@ -1,8 +0,0 @@ - - - - - - - - diff --git a/examples/Component/Mapping/NonLinear/DistanceFromTargetMapping.scn b/examples/Component/Mapping/NonLinear/DistanceFromTargetMapping.scn index 0a2144bea18..1dd85728358 100644 --- a/examples/Component/Mapping/NonLinear/DistanceFromTargetMapping.scn +++ b/examples/Component/Mapping/NonLinear/DistanceFromTargetMapping.scn @@ -1,7 +1,7 @@ - + @@ -15,7 +15,7 @@ - + diff --git a/examples/Component/SolidMechanics/FEM/RestShapeSpringsForceField3.scn b/examples/Component/SolidMechanics/FEM/RestShapeSpringsForceField3.scn index b7680484750..6b034813594 100644 --- a/examples/Component/SolidMechanics/FEM/RestShapeSpringsForceField3.scn +++ b/examples/Component/SolidMechanics/FEM/RestShapeSpringsForceField3.scn @@ -5,7 +5,7 @@ - + @@ -21,7 +21,7 @@ - + - + @@ -23,7 +23,7 @@ - + diff --git a/examples/RegressionStateScenes.regression-tests b/examples/RegressionStateScenes.regression-tests index f430b475ae6..7ac5834edf8 100644 --- a/examples/RegressionStateScenes.regression-tests +++ b/examples/RegressionStateScenes.regression-tests @@ -33,9 +33,7 @@ Component/LinearSolver/Direct/FEMBAR_EigenSimplicialLLT.scn 100 1e-8 0 1 Component/LinearSolver/Direct/FEMBAR_EigenSparseLU.scn 100 1e-8 0 1 Component/LinearSolver/Direct/FEMBAR_EigenSparseQR.scn 100 1e-8 0 1 Component/LinearSolver/Direct/FEMBAR_PrecomputedLinearSolver.scn 100 1e-8 0 1 -Component/LinearSolver/Direct/FEMBAR_SparseCholeskySolver.scn 100 1e-8 0 1 Component/LinearSolver/Direct/FEMBAR_SparseLDLSolver.scn 100 1e-8 0 1 -Component/LinearSolver/Direct/FEMBAR_SparseLUSolver.scn 100 1e-8 0 1 Component/LinearSolver/Iterative/FEMBAR_MinResLinearSolver.scn 100 1e-8 0 1 Component/LinearSolver/Iterative/FEMBAR_ShewchukPCGLinearSolver.scn 1e-8 0 1 Component/LinearSolver/Iterative/FEMBAR_CGLinearSolver.scn 100 1e-8 0 1