Eigen: fold remaining OpenNL code into intern/eigen.
authorBrecht Van Lommel <brechtvanlommel@gmail.com>
Tue, 24 Nov 2015 19:42:10 +0000 (20:42 +0100)
committerBrecht Van Lommel <brechtvanlommel@gmail.com>
Thu, 10 Dec 2015 00:58:10 +0000 (01:58 +0100)
Differential Revision: https://developer.blender.org/D1662

56 files changed:
CMakeLists.txt
SConstruct
build_files/cmake/cmake_static_check_clang_array.py
build_files/cmake/cmake_static_check_cppcheck.py
build_files/cmake/cmake_static_check_smatch.py
build_files/cmake/cmake_static_check_sparse.py
build_files/cmake/cmake_static_check_splint.py
build_files/cmake/config/blender_full.cmake
build_files/cmake/config/blender_lite.cmake
build_files/cmake/macros.cmake
doc/doxygen/doxygen.intern.h
extern/CMakeLists.txt
extern/SConscript
extern/colamd/CMakeLists.txt [deleted file]
extern/colamd/Doc/ChangeLog [deleted file]
extern/colamd/Doc/lesser.txt [deleted file]
extern/colamd/Include/UFconfig.h [deleted file]
extern/colamd/Include/colamd.h [deleted file]
extern/colamd/README.txt [deleted file]
extern/colamd/SConscript [deleted file]
extern/colamd/Source/colamd.c [deleted file]
extern/colamd/Source/colamd_global.c [deleted file]
intern/CMakeLists.txt
intern/SConscript
intern/eigen/CMakeLists.txt
intern/eigen/eigen_capi.h
intern/eigen/intern/eigenvalues.cc
intern/eigen/intern/eigenvalues.h
intern/eigen/intern/linear_solver.cc [new file with mode: 0644]
intern/eigen/intern/linear_solver.h [new file with mode: 0644]
intern/eigen/intern/svd.cc
intern/eigen/intern/svd.h
intern/opennl/CMakeLists.txt [deleted file]
intern/opennl/SConscript [deleted file]
intern/opennl/extern/ONL_opennl.h [deleted file]
intern/opennl/intern/opennl.cpp [deleted file]
source/blender/blenkernel/intern/softbody.c
source/blender/blenlib/intern/math_solvers.c
source/blender/bmesh/CMakeLists.txt
source/blender/bmesh/SConscript
source/blender/bmesh/operators/bmo_smooth_laplacian.c
source/blender/editors/armature/CMakeLists.txt
source/blender/editors/armature/SConscript
source/blender/editors/armature/armature_skinning.c
source/blender/editors/armature/meshlaplacian.c
source/blender/editors/armature/reeb.c
source/blender/editors/uvedit/CMakeLists.txt
source/blender/editors/uvedit/SConscript
source/blender/editors/uvedit/uvedit_parametrizer.c
source/blender/modifiers/CMakeLists.txt
source/blender/modifiers/SConscript
source/blender/modifiers/intern/MOD_laplaciandeform.c
source/blender/modifiers/intern/MOD_laplaciansmooth.c
source/blender/modifiers/intern/MOD_util.c
source/blender/modifiers/intern/MOD_util.h
source/blenderplayer/CMakeLists.txt

index 26e0a3f13799a823a522387e2fea691e31e00341..1cf5d1c4620dff4e6a41b0b5ed37b594c8e304b1 100644 (file)
@@ -362,7 +362,6 @@ if(WIN32)
 endif()
 option(WITH_INPUT_NDOF "Enable NDOF input devices (SpaceNavigator and friends)" ${_init_INPUT_NDOF})
 option(WITH_RAYOPTIMIZATION    "Enable use of SIMD (SSE) optimizations for the raytracer" ON)
-option(WITH_OPENNL        "Enable use of Open Numerical Library" ON)
 if(UNIX AND NOT APPLE)
        option(WITH_INSTALL_PORTABLE "Install redistributeable runtime, otherwise install into CMAKE_INSTALL_PREFIX" ON)
        option(WITH_STATIC_LIBS "Try to link with static libraries, as much as possible, to make blender more portable across distributions" OFF)
@@ -2983,9 +2982,6 @@ if(FIRST_RUN)
                info_cfg_option(WITH_GL_ANGLE)
        endif()
 
-       info_cfg_text("Other:")
-       info_cfg_option(WITH_OPENNL)
-
        # debug
        message(STATUS "HAVE_STDBOOL_H = ${HAVE_STDBOOL_H}")
 
index a5efed552fe598a19c088b76f1aa345b6281dc32..de265df5b0232888bd873190e446016e0ce0e39a 100644 (file)
@@ -555,7 +555,6 @@ else:
 
 # TODO, make optional (as with CMake)
 env['CPPFLAGS'].append('-DWITH_AVI')
-env['CPPFLAGS'].append('-DWITH_OPENNL')
 
 if env['OURPLATFORM'] not in ('win32-vc', 'win64-vc'):
     env['CPPFLAGS'].append('-DHAVE_STDBOOL_H')
index 45b262a13ce289df3599349230157dfe21e9fda0..597d1d2b98002a12fba545a37ad717929e5df953 100644 (file)
@@ -32,7 +32,6 @@ USE_QUIET = (os.environ.get("QUIET", None) is not None)
 CHECKER_IGNORE_PREFIX = [
     "extern",
     "intern/moto",
-    "blender/intern/opennl",
     ]
 
 CHECKER_BIN = "python2"
index 9f012cb7f6de319302de43f64ce723791b97d56b..3504b005adcf7e9ea554821401f4be172a3c1b2c 100644 (file)
@@ -32,7 +32,6 @@ USE_QUIET = (os.environ.get("QUIET", None) is not None)
 CHECKER_IGNORE_PREFIX = [
     "extern",
     "intern/moto",
-    "blender/intern/opennl",
     ]
 
 CHECKER_BIN = "cppcheck"
index de13d141276f309a9ec0cc3b8d889efbe09c3dd1..f525187e22cd2db53cdbf01e699a7e1648b735ff 100644 (file)
@@ -25,7 +25,6 @@
 CHECKER_IGNORE_PREFIX = [
     "extern",
     "intern/moto",
-    "blender/intern/opennl",
     ]
 
 CHECKER_BIN = "smatch"
index 2ee99925adb68d2358168bc29177e420e562e922..6d1c4e3d6a27d8ab13d91f071ac87800011f16f2 100644 (file)
@@ -25,7 +25,6 @@
 CHECKER_IGNORE_PREFIX = [
     "extern",
     "intern/moto",
-    "blender/intern/opennl",
     ]
 
 CHECKER_BIN = "sparse"
index 5a967ecc7a3f026786f7a04af48c5b0d5e989c47..46d8808e89d1938adbebfd33e49cd957cec248bd 100644 (file)
@@ -25,7 +25,6 @@
 CHECKER_IGNORE_PREFIX = [
     "extern",
     "intern/moto",
-    "blender/intern/opennl",
     ]
 
 CHECKER_BIN = "splint"
index 0ca980050db6623c45dad173f0cee821e6240e29..ad8a681567526a1f4be71cf8aa9f91ca9394a910 100644 (file)
@@ -43,7 +43,6 @@ set(WITH_OPENAL              ON  CACHE BOOL "" FORCE)
 set(WITH_OPENCOLLADA         ON  CACHE BOOL "" FORCE)
 set(WITH_OPENCOLORIO         ON  CACHE BOOL "" FORCE)
 set(WITH_OPENMP              ON  CACHE BOOL "" FORCE)
-set(WITH_OPENNL              ON  CACHE BOOL "" FORCE)
 set(WITH_PYTHON_INSTALL      ON  CACHE BOOL "" FORCE)
 set(WITH_RAYOPTIMIZATION     ON  CACHE BOOL "" FORCE)
 set(WITH_SDL                 ON  CACHE BOOL "" FORCE)
index 9e7b4dcf39bb6545c0dd022a2a7498379e39b279..99e90caf7938056970f33723d3eb51059df0fd8a 100644 (file)
@@ -47,7 +47,6 @@ set(WITH_OPENCOLLADA         OFF CACHE BOOL "" FORCE)
 set(WITH_OPENCOLORIO         OFF CACHE BOOL "" FORCE)
 set(WITH_OPENIMAGEIO         OFF CACHE BOOL "" FORCE)
 set(WITH_OPENMP              OFF CACHE BOOL "" FORCE)
-set(WITH_OPENNL              OFF CACHE BOOL "" FORCE)
 set(WITH_RAYOPTIMIZATION     OFF CACHE BOOL "" FORCE)
 set(WITH_SDL                 OFF CACHE BOOL "" FORCE)
 set(WITH_X11_XINPUT          OFF CACHE BOOL "" FORCE)
index 39c12f2c4d922adc5fdf0e986ade00d2b0ba0bd6..4ba15c726778d5e6e861b3e367e5ad79c67b402b 100644 (file)
@@ -596,7 +596,6 @@ function(SETUP_BLENDER_SORTED_LIBS)
                ge_phys_bullet
                bf_intern_smoke
                extern_lzma
-               extern_colamd
                ge_logic_ketsji
                extern_recastnavigation
                ge_logic
@@ -698,10 +697,6 @@ function(SETUP_BLENDER_SORTED_LIBS)
                list(APPEND BLENDER_SORTED_LIBS bf_intern_locale)
        endif()
 
-       if(WITH_OPENNL)
-               list_insert_after(BLENDER_SORTED_LIBS "bf_render" "bf_intern_opennl")
-       endif()
-
        if(WITH_BULLET)
                list_insert_after(BLENDER_SORTED_LIBS "bf_blenkernel" "bf_intern_rigidbody")
        endif()
index 2c8ecae0ead5f1904a730dbf38888fb35e3648ba..e3cc11b3404d335ffa7cb3f264daef9df3e94ac0 100644 (file)
@@ -38,7 +38,7 @@
  *  \ingroup intern
  */
 
-/** \defgroup opennl opennl
+/** \defgroup eigen eigen
  *  \ingroup intern
  */
 
index d0c587b80e41303ae7f715325098d30dd5df2bbc..640de9d80e7161ba89e54731e873b8d5b9b80c3e 100644 (file)
@@ -30,10 +30,6 @@ add_subdirectory(rangetree)
 add_subdirectory(wcwidth)
 add_subdirectory(libmv)
 
-if(WITH_OPENNL)
-       add_subdirectory(colamd)
-endif()
-
 if(WITH_BULLET)
        if(NOT WITH_SYSTEM_BULLET)
                add_subdirectory(bullet2)
index 46c177c5bcba085bbdd3f2f6f4859f6150070f96..a5d8c1f078e1431410e2d898d4a0bc2dbe03f8c6 100644 (file)
@@ -7,7 +7,6 @@ if env['WITH_BF_GLEW_ES']:
 else:
     SConscript(['glew/SConscript'])
 
-SConscript(['colamd/SConscript'])
 SConscript(['rangetree/SConscript'])
 SConscript(['wcwidth/SConscript'])
 SConscript(['libmv/SConscript'])
diff --git a/extern/colamd/CMakeLists.txt b/extern/colamd/CMakeLists.txt
deleted file mode 100644 (file)
index 3019ee5..0000000
+++ /dev/null
@@ -1,41 +0,0 @@
-# ***** BEGIN GPL LICENSE BLOCK *****
-#
-# This program is free software; you can redistribute it and/or
-# modify it under the terms of the GNU General Public License
-# as published by the Free Software Foundation; either version 2
-# 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 General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program; if not, write to the Free Software Foundation,
-# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-#
-# The Original Code is Copyright (C) 2011, Blender Foundation
-# All rights reserved.
-#
-# Contributor(s): Blender Foundation,
-#                 Sergey Sharybin
-#
-# ***** END GPL LICENSE BLOCK *****
-
-set(INC
-       Include
-)
-
-set(INC_SYS
-
-)
-
-set(SRC
-       Source/colamd.c
-       Source/colamd_global.c
-
-       Include/colamd.h
-       Include/UFconfig.h
-)
-
-blender_add_lib(extern_colamd "${SRC}" "${INC}" "${INC_SYS}")
diff --git a/extern/colamd/Doc/ChangeLog b/extern/colamd/Doc/ChangeLog
deleted file mode 100644 (file)
index 29308e9..0000000
+++ /dev/null
@@ -1,129 +0,0 @@
-May 31, 2007: version 2.7.0
-
-    * ported to 64-bit MATLAB
-
-    * subdirectories added (Source/, Include/, Lib/, Doc/, MATLAB/, Demo/)
-
-Dec 12, 2006, version 2.5.2
-
-    * minor MATLAB cleanup.  MATLAB functions renamed colamd2 and symamd2,
-       so that they do not conflict with the built-in versions.  Note that
-       the MATLAB built-in functions colamd and symamd are identical to
-       the colamd and symamd functions here.
-
-Aug 31, 2006: Version 2.5.1
-
-    * minor change to colamd.m and symamd.m, to use etree instead
-       of sparsfun.
-
-Apr. 30, 2006: Version 2.5
-
-    * colamd_recommended modified, to do more careful integer overflow
-       checking.  It now returns size_t, not int.  colamd_l_recommended
-       also returns size_t.  A zero is returned if an error occurs.  A
-       postive return value denotes success.  In v2.4 and earlier,
-       -1 was returned on error (an int or long).
-
-    * long replaced with UF_long integer, which is long except on WIN64.
-
-Nov 15, 2005:
-
-    * minor editting of comments; version number (2.4) unchanged.
-
-Changes from Version 2.3 to 2.4 (Aug 30, 2005)
-
-    * Makefile now relies on ../UFconfig/UFconfig.mk
-
-    * changed the dense row/col detection.  The meaning of the knobs
-       has thus changed.
-
-    * added an option to turn off aggressive absorption.  It was
-       always on in versions 2.3 and earlier.
-
-    * added a #define'd version number
-
-    * added a function pointer (colamd_printf) for COLAMD's printing.
-
-    * added a -DNPRINT option, to turn off printing at compile-time.
-
-    * added a check for integer overflow in colamd_recommended
-
-    * minor changes to allow for more simpler 100% test coverage
-
-    * bug fix.  If symamd v2.3 fails to allocate its copy of the input
-       matrix, then it erroneously frees a calloc'd workspace twice.
-       This bug has no effect on the MATLAB symamd mexFunction, since
-       mxCalloc terminates the mexFunction if it fails to allocate
-       memory.  Similarly, UMFPACK is not affected because it does not
-       use symamd.  The bug has no effect on the colamd ordering
-       routine in v2.3.
-
-Changes from Version 2.2 to 2.3 (Sept. 8, 2003)
-
-    * removed the call to the MATLAB spparms ('spumoni') function.
-       This can take a lot of time if you are ordering many small
-       matrices.  Only affects the MATLAB interface (colamdmex.c,
-       symamdmex.c, colamdtestmex.c, and symamdtestmex.c).  The
-       usage of the optional 2nd argument to the colamd and symamd
-       mexFunctions was changed accordingly.
-
-Changes from Version 2.1 to 2.2 (Sept. 23, 2002)
-
-    * extensive testing routines added (colamd_test.m, colamdtestmex.c,
-       and symamdtestmex.c), and the Makefile modified accordingly.
-
-    * a few typos in the comments corrected 
-
-    * use of the MATLAB "flops" command removed from colamd_demo, and an
-       m-file routine luflops.m added.
-
-    * an explicit typecast from unsigned to int added, for COLAMD_C and
-       COLAMD_R in colamd.h.
-
-    * #include <stdio.h> added to colamd_example.c
-
-
-Changes from Version 2.0 to 2.1 (May 4, 2001)
-
-    * TRUE and FALSE are predefined on some systems, so they are defined
-           here only if not already defined.
-    
-    * web site changed
-
-    * UNIX Makefile modified, to handle the case if "." is not in your path.
-
-
-Changes from Version 1.0 to 2.0 (January 31, 2000)
-
-    No bugs were found in version 1.1.  These changes merely add new
-    functionality.
-
-    * added the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro.
-
-    * moved the output statistics, from A, to a separate output argument.
-           The arguments changed for the C-callable routines.
-
-    * added colamd_report and symamd_report.
-
-    * added a C-callable symamd routine.  Formerly, symamd was only
-           available as a mexFunction from MATLAB.
-
-    * added error-checking to symamd.  Formerly, it assumed its input
-           was error-free.
-
-    * added the optional stats and knobs arguments to the symamd mexFunction
-
-    * deleted colamd_help.  A help message is still available from
-           "help colamd" and "help symamd" in MATLAB.
-
-    * deleted colamdtree.m and symamdtree.m.  Now, colamd.m and symamd.m
-           also do the elimination tree post-ordering.  The Version 1.1
-           colamd and symamd mexFunctions, which do not do the post-
-           ordering, are now visible as colamdmex and symamdmex from
-           MATLAB.  Essentialy, the post-ordering is now the default
-           behavior of colamd.m and symamd.m, to match the behavior of
-           colmmd and symmmd.  The post-ordering is only available in the
-           MATLAB interface, not the C-callable interface.
-
-    * made a slight change to the dense row/column detection in symamd,
-           to match the stated specifications.
diff --git a/extern/colamd/Doc/lesser.txt b/extern/colamd/Doc/lesser.txt
deleted file mode 100644 (file)
index 8add30a..0000000
+++ /dev/null
@@ -1,504 +0,0 @@
-                 GNU LESSER GENERAL PUBLIC LICENSE
-                      Version 2.1, February 1999
-
- Copyright (C) 1991, 1999 Free Software Foundation, Inc.
-     51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
- Everyone is permitted to copy and distribute verbatim copies
- of this license document, but changing it is not allowed.
-
-[This is the first released version of the Lesser GPL.  It also counts
- as the successor of the GNU Library Public License, version 2, hence
- the version number 2.1.]
-
-                           Preamble
-
-  The licenses for most software are designed to take away your
-freedom to share and change it.  By contrast, the GNU General Public
-Licenses are intended to guarantee your freedom to share and change
-free software--to make sure the software is free for all its users.
-
-  This license, the Lesser General Public License, applies to some
-specially designated software packages--typically libraries--of the
-Free Software Foundation and other authors who decide to use it.  You
-can use it too, but we suggest you first think carefully about whether
-this license or the ordinary General Public License is the better
-strategy to use in any particular case, based on the explanations below.
-
-  When we speak of free software, we are referring to freedom of use,
-not price.  Our General Public Licenses are designed to make sure that
-you have the freedom to distribute copies of free software (and charge
-for this service if you wish); that you receive source code or can get
-it if you want it; that you can change the software and use pieces of
-it in new free programs; and that you are informed that you can do
-these things.
-
-  To protect your rights, we need to make restrictions that forbid
-distributors to deny you these rights or to ask you to surrender these
-rights.  These restrictions translate to certain responsibilities for
-you if you distribute copies of the library or if you modify it.
-
-  For example, if you distribute copies of the library, whether gratis
-or for a fee, you must give the recipients all the rights that we gave
-you.  You must make sure that they, too, receive or can get the source
-code.  If you link other code with the library, you must provide
-complete object files to the recipients, so that they can relink them
-with the library after making changes to the library and recompiling
-it.  And you must show them these terms so they know their rights.
-
-  We protect your rights with a two-step method: (1) we copyright the
-library, and (2) we offer you this license, which gives you legal
-permission to copy, distribute and/or modify the library.
-
-  To protect each distributor, we want to make it very clear that
-there is no warranty for the free library.  Also, if the library is
-modified by someone else and passed on, the recipients should know
-that what they have is not the original version, so that the original
-author's reputation will not be affected by problems that might be
-introduced by others.
-\f
-  Finally, software patents pose a constant threat to the existence of
-any free program.  We wish to make sure that a company cannot
-effectively restrict the users of a free program by obtaining a
-restrictive license from a patent holder.  Therefore, we insist that
-any patent license obtained for a version of the library must be
-consistent with the full freedom of use specified in this license.
-
-  Most GNU software, including some libraries, is covered by the
-ordinary GNU General Public License.  This license, the GNU Lesser
-General Public License, applies to certain designated libraries, and
-is quite different from the ordinary General Public License.  We use
-this license for certain libraries in order to permit linking those
-libraries into non-free programs.
-
-  When a program is linked with a library, whether statically or using
-a shared library, the combination of the two is legally speaking a
-combined work, a derivative of the original library.  The ordinary
-General Public License therefore permits such linking only if the
-entire combination fits its criteria of freedom.  The Lesser General
-Public License permits more lax criteria for linking other code with
-the library.
-
-  We call this license the "Lesser" General Public License because it
-does Less to protect the user's freedom than the ordinary General
-Public License.  It also provides other free software developers Less
-of an advantage over competing non-free programs.  These disadvantages
-are the reason we use the ordinary General Public License for many
-libraries.  However, the Lesser license provides advantages in certain
-special circumstances.
-
-  For example, on rare occasions, there may be a special need to
-encourage the widest possible use of a certain library, so that it becomes
-a de-facto standard.  To achieve this, non-free programs must be
-allowed to use the library.  A more frequent case is that a free
-library does the same job as widely used non-free libraries.  In this
-case, there is little to gain by limiting the free library to free
-software only, so we use the Lesser General Public License.
-
-  In other cases, permission to use a particular library in non-free
-programs enables a greater number of people to use a large body of
-free software.  For example, permission to use the GNU C Library in
-non-free programs enables many more people to use the whole GNU
-operating system, as well as its variant, the GNU/Linux operating
-system.
-
-  Although the Lesser General Public License is Less protective of the
-users' freedom, it does ensure that the user of a program that is
-linked with the Library has the freedom and the wherewithal to run
-that program using a modified version of the Library.
-
-  The precise terms and conditions for copying, distribution and
-modification follow.  Pay close attention to the difference between a
-"work based on the library" and a "work that uses the library".  The
-former contains code derived from the library, whereas the latter must
-be combined with the library in order to run.
-\f
-                 GNU LESSER GENERAL PUBLIC LICENSE
-   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
-
-  0. This License Agreement applies to any software library or other
-program which contains a notice placed by the copyright holder or
-other authorized party saying it may be distributed under the terms of
-this Lesser General Public License (also called "this License").
-Each licensee is addressed as "you".
-
-  A "library" means a collection of software functions and/or data
-prepared so as to be conveniently linked with application programs
-(which use some of those functions and data) to form executables.
-
-  The "Library", below, refers to any such software library or work
-which has been distributed under these terms.  A "work based on the
-Library" means either the Library or any derivative work under
-copyright law: that is to say, a work containing the Library or a
-portion of it, either verbatim or with modifications and/or translated
-straightforwardly into another language.  (Hereinafter, translation is
-included without limitation in the term "modification".)
-
-  "Source code" for a work means the preferred form of the work for
-making modifications to it.  For a library, complete source code means
-all the source code for all modules it contains, plus any associated
-interface definition files, plus the scripts used to control compilation
-and installation of the library.
-
-  Activities other than copying, distribution and modification are not
-covered by this License; they are outside its scope.  The act of
-running a program using the Library is not restricted, and output from
-such a program is covered only if its contents constitute a work based
-on the Library (independent of the use of the Library in a tool for
-writing it).  Whether that is true depends on what the Library does
-and what the program that uses the Library does.
-  
-  1. You may copy and distribute verbatim copies of the Library's
-complete source code as you receive it, in any medium, provided that
-you conspicuously and appropriately publish on each copy an
-appropriate copyright notice and disclaimer of warranty; keep intact
-all the notices that refer to this License and to the absence of any
-warranty; and distribute a copy of this License along with the
-Library.
-
-  You may charge a fee for the physical act of transferring a copy,
-and you may at your option offer warranty protection in exchange for a
-fee.
-\f
-  2. You may modify your copy or copies of the Library or any portion
-of it, thus forming a work based on the Library, and copy and
-distribute such modifications or work under the terms of Section 1
-above, provided that you also meet all of these conditions:
-
-    a) The modified work must itself be a software library.
-
-    b) You must cause the files modified to carry prominent notices
-    stating that you changed the files and the date of any change.
-
-    c) You must cause the whole of the work to be licensed at no
-    charge to all third parties under the terms of this License.
-
-    d) If a facility in the modified Library refers to a function or a
-    table of data to be supplied by an application program that uses
-    the facility, other than as an argument passed when the facility
-    is invoked, then you must make a good faith effort to ensure that,
-    in the event an application does not supply such function or
-    table, the facility still operates, and performs whatever part of
-    its purpose remains meaningful.
-
-    (For example, a function in a library to compute square roots has
-    a purpose that is entirely well-defined independent of the
-    application.  Therefore, Subsection 2d requires that any
-    application-supplied function or table used by this function must
-    be optional: if the application does not supply it, the square
-    root function must still compute square roots.)
-
-These requirements apply to the modified work as a whole.  If
-identifiable sections of that work are not derived from the Library,
-and can be reasonably considered independent and separate works in
-themselves, then this License, and its terms, do not apply to those
-sections when you distribute them as separate works.  But when you
-distribute the same sections as part of a whole which is a work based
-on the Library, the distribution of the whole must be on the terms of
-this License, whose permissions for other licensees extend to the
-entire whole, and thus to each and every part regardless of who wrote
-it.
-
-Thus, it is not the intent of this section to claim rights or contest
-your rights to work written entirely by you; rather, the intent is to
-exercise the right to control the distribution of derivative or
-collective works based on the Library.
-
-In addition, mere aggregation of another work not based on the Library
-with the Library (or with a work based on the Library) on a volume of
-a storage or distribution medium does not bring the other work under
-the scope of this License.
-
-  3. You may opt to apply the terms of the ordinary GNU General Public
-License instead of this License to a given copy of the Library.  To do
-this, you must alter all the notices that refer to this License, so
-that they refer to the ordinary GNU General Public License, version 2,
-instead of to this License.  (If a newer version than version 2 of the
-ordinary GNU General Public License has appeared, then you can specify
-that version instead if you wish.)  Do not make any other change in
-these notices.
-\f
-  Once this change is made in a given copy, it is irreversible for
-that copy, so the ordinary GNU General Public License applies to all
-subsequent copies and derivative works made from that copy.
-
-  This option is useful when you wish to copy part of the code of
-the Library into a program that is not a library.
-
-  4. You may copy and distribute the Library (or a portion or
-derivative of it, under Section 2) in object code or executable form
-under the terms of Sections 1 and 2 above provided that you accompany
-it with the complete corresponding machine-readable source code, which
-must be distributed under the terms of Sections 1 and 2 above on a
-medium customarily used for software interchange.
-
-  If distribution of object code is made by offering access to copy
-from a designated place, then offering equivalent access to copy the
-source code from the same place satisfies the requirement to
-distribute the source code, even though third parties are not
-compelled to copy the source along with the object code.
-
-  5. A program that contains no derivative of any portion of the
-Library, but is designed to work with the Library by being compiled or
-linked with it, is called a "work that uses the Library".  Such a
-work, in isolation, is not a derivative work of the Library, and
-therefore falls outside the scope of this License.
-
-  However, linking a "work that uses the Library" with the Library
-creates an executable that is a derivative of the Library (because it
-contains portions of the Library), rather than a "work that uses the
-library".  The executable is therefore covered by this License.
-Section 6 states terms for distribution of such executables.
-
-  When a "work that uses the Library" uses material from a header file
-that is part of the Library, the object code for the work may be a
-derivative work of the Library even though the source code is not.
-Whether this is true is especially significant if the work can be
-linked without the Library, or if the work is itself a library.  The
-threshold for this to be true is not precisely defined by law.
-
-  If such an object file uses only numerical parameters, data
-structure layouts and accessors, and small macros and small inline
-functions (ten lines or less in length), then the use of the object
-file is unrestricted, regardless of whether it is legally a derivative
-work.  (Executables containing this object code plus portions of the
-Library will still fall under Section 6.)
-
-  Otherwise, if the work is a derivative of the Library, you may
-distribute the object code for the work under the terms of Section 6.
-Any executables containing that work also fall under Section 6,
-whether or not they are linked directly with the Library itself.
-\f
-  6. As an exception to the Sections above, you may also combine or
-link a "work that uses the Library" with the Library to produce a
-work containing portions of the Library, and distribute that work
-under terms of your choice, provided that the terms permit
-modification of the work for the customer's own use and reverse
-engineering for debugging such modifications.
-
-  You must give prominent notice with each copy of the work that the
-Library is used in it and that the Library and its use are covered by
-this License.  You must supply a copy of this License.  If the work
-during execution displays copyright notices, you must include the
-copyright notice for the Library among them, as well as a reference
-directing the user to the copy of this License.  Also, you must do one
-of these things:
-
-    a) Accompany the work with the complete corresponding
-    machine-readable source code for the Library including whatever
-    changes were used in the work (which must be distributed under
-    Sections 1 and 2 above); and, if the work is an executable linked
-    with the Library, with the complete machine-readable "work that
-    uses the Library", as object code and/or source code, so that the
-    user can modify the Library and then relink to produce a modified
-    executable containing the modified Library.  (It is understood
-    that the user who changes the contents of definitions files in the
-    Library will not necessarily be able to recompile the application
-    to use the modified definitions.)
-
-    b) Use a suitable shared library mechanism for linking with the
-    Library.  A suitable mechanism is one that (1) uses at run time a
-    copy of the library already present on the user's computer system,
-    rather than copying library functions into the executable, and (2)
-    will operate properly with a modified version of the library, if
-    the user installs one, as long as the modified version is
-    interface-compatible with the version that the work was made with.
-
-    c) Accompany the work with a written offer, valid for at
-    least three years, to give the same user the materials
-    specified in Subsection 6a, above, for a charge no more
-    than the cost of performing this distribution.
-
-    d) If distribution of the work is made by offering access to copy
-    from a designated place, offer equivalent access to copy the above
-    specified materials from the same place.
-
-    e) Verify that the user has already received a copy of these
-    materials or that you have already sent this user a copy.
-
-  For an executable, the required form of the "work that uses the
-Library" must include any data and utility programs needed for
-reproducing the executable from it.  However, as a special exception,
-the materials to be distributed need not include anything that is
-normally distributed (in either source or binary form) with the major
-components (compiler, kernel, and so on) of the operating system on
-which the executable runs, unless that component itself accompanies
-the executable.
-
-  It may happen that this requirement contradicts the license
-restrictions of other proprietary libraries that do not normally
-accompany the operating system.  Such a contradiction means you cannot
-use both them and the Library together in an executable that you
-distribute.
-\f
-  7. You may place library facilities that are a work based on the
-Library side-by-side in a single library together with other library
-facilities not covered by this License, and distribute such a combined
-library, provided that the separate distribution of the work based on
-the Library and of the other library facilities is otherwise
-permitted, and provided that you do these two things:
-
-    a) Accompany the combined library with a copy of the same work
-    based on the Library, uncombined with any other library
-    facilities.  This must be distributed under the terms of the
-    Sections above.
-
-    b) Give prominent notice with the combined library of the fact
-    that part of it is a work based on the Library, and explaining
-    where to find the accompanying uncombined form of the same work.
-
-  8. You may not copy, modify, sublicense, link with, or distribute
-the Library except as expressly provided under this License.  Any
-attempt otherwise to copy, modify, sublicense, link with, or
-distribute the Library is void, and will automatically terminate your
-rights under this License.  However, parties who have received copies,
-or rights, from you under this License will not have their licenses
-terminated so long as such parties remain in full compliance.
-
-  9. You are not required to accept this License, since you have not
-signed it.  However, nothing else grants you permission to modify or
-distribute the Library or its derivative works.  These actions are
-prohibited by law if you do not accept this License.  Therefore, by
-modifying or distributing the Library (or any work based on the
-Library), you indicate your acceptance of this License to do so, and
-all its terms and conditions for copying, distributing or modifying
-the Library or works based on it.
-
-  10. Each time you redistribute the Library (or any work based on the
-Library), the recipient automatically receives a license from the
-original licensor to copy, distribute, link with or modify the Library
-subject to these terms and conditions.  You may not impose any further
-restrictions on the recipients' exercise of the rights granted herein.
-You are not responsible for enforcing compliance by third parties with
-this License.
-\f
-  11. If, as a consequence of a court judgment or allegation of patent
-infringement or for any other reason (not limited to patent issues),
-conditions are imposed on you (whether by court order, agreement or
-otherwise) that contradict the conditions of this License, they do not
-excuse you from the conditions of this License.  If you cannot
-distribute so as to satisfy simultaneously your obligations under this
-License and any other pertinent obligations, then as a consequence you
-may not distribute the Library at all.  For example, if a patent
-license would not permit royalty-free redistribution of the Library by
-all those who receive copies directly or indirectly through you, then
-the only way you could satisfy both it and this License would be to
-refrain entirely from distribution of the Library.
-
-If any portion of this section is held invalid or unenforceable under any
-particular circumstance, the balance of the section is intended to apply,
-and the section as a whole is intended to apply in other circumstances.
-
-It is not the purpose of this section to induce you to infringe any
-patents or other property right claims or to contest validity of any
-such claims; this section has the sole purpose of protecting the
-integrity of the free software distribution system which is
-implemented by public license practices.  Many people have made
-generous contributions to the wide range of software distributed
-through that system in reliance on consistent application of that
-system; it is up to the author/donor to decide if he or she is willing
-to distribute software through any other system and a licensee cannot
-impose that choice.
-
-This section is intended to make thoroughly clear what is believed to
-be a consequence of the rest of this License.
-
-  12. If the distribution and/or use of the Library is restricted in
-certain countries either by patents or by copyrighted interfaces, the
-original copyright holder who places the Library under this License may add
-an explicit geographical distribution limitation excluding those countries,
-so that distribution is permitted only in or among countries not thus
-excluded.  In such case, this License incorporates the limitation as if
-written in the body of this License.
-
-  13. The Free Software Foundation may publish revised and/or new
-versions of the Lesser General Public License from time to time.
-Such new versions will be similar in spirit to the present version,
-but may differ in detail to address new problems or concerns.
-
-Each version is given a distinguishing version number.  If the Library
-specifies a version number of this License which applies to it and
-"any later version", you have the option of following the terms and
-conditions either of that version or of any later version published by
-the Free Software Foundation.  If the Library does not specify a
-license version number, you may choose any version ever published by
-the Free Software Foundation.
-\f
-  14. If you wish to incorporate parts of the Library into other free
-programs whose distribution conditions are incompatible with these,
-write to the author to ask for permission.  For software which is
-copyrighted by the Free Software Foundation, write to the Free
-Software Foundation; we sometimes make exceptions for this.  Our
-decision will be guided by the two goals of preserving the free status
-of all derivatives of our free software and of promoting the sharing
-and reuse of software generally.
-
-                           NO WARRANTY
-
-  15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO
-WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW.
-EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR
-OTHER PARTIES PROVIDE THE LIBRARY "AS IS" WITHOUT WARRANTY OF ANY
-KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE
-IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE
-LIBRARY IS WITH YOU.  SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME
-THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
-
-  16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN
-WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY
-AND/OR REDISTRIBUTE THE LIBRARY AS PERMITTED ABOVE, BE LIABLE TO YOU
-FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR
-CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE
-LIBRARY (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING
-RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A
-FAILURE OF THE LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF
-SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH
-DAMAGES.
-
-                    END OF TERMS AND CONDITIONS
-\f
-           How to Apply These Terms to Your New Libraries
-
-  If you develop a new library, and you want it to be of the greatest
-possible use to the public, we recommend making it free software that
-everyone can redistribute and change.  You can do so by permitting
-redistribution under these terms (or, alternatively, under the terms of the
-ordinary General Public License).
-
-  To apply these terms, attach the following notices to the library.  It is
-safest to attach them to the start of each source file to most effectively
-convey the exclusion of warranty; and each file should have at least the
-"copyright" line and a pointer to where the full notice is found.
-
-    <one line to give the library's name and a brief idea of what it does.>
-    Copyright (C) <year>  <name of author>
-
-    This library 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 library 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 library; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-
-Also add information on how to contact you by electronic and paper mail.
-
-You should also get your employer (if you work as a programmer) or your
-school, if any, to sign a "copyright disclaimer" for the library, if
-necessary.  Here is a sample; alter the names:
-
-  Yoyodyne, Inc., hereby disclaims all copyright interest in the
-  library `Frob' (a library for tweaking knobs) written by James Random Hacker.
-
-  <signature of Ty Coon>, 1 April 1990
-  Ty Coon, President of Vice
-
-That's all there is to it!
-
-
diff --git a/extern/colamd/Include/UFconfig.h b/extern/colamd/Include/UFconfig.h
deleted file mode 100644 (file)
index 7b5e79e..0000000
+++ /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 <limits.h>
-
-/* ========================================================================== */
-/* === UF_long ============================================================== */
-/* ========================================================================== */
-
-#ifndef UF_long
-
-#ifdef _WIN64
-
-#define UF_long __int64
-#define UF_long_max _I64_MAX
-#define UF_long_id "%I64d"
-
-#else
-
-#define UF_long long
-#define UF_long_max LONG_MAX
-#define UF_long_id "%ld"
-
-#endif
-#endif
-
-/* ========================================================================== */
-/* === SuiteSparse version ================================================== */
-/* ========================================================================== */
-
-/* SuiteSparse is not a package itself, but a collection of packages, some of
- * which must be used together (UMFPACK requires AMD, CHOLMOD requires AMD,
- * COLAMD, CAMD, and CCOLAMD, etc).  A version number is provided here for the
- * collection itself.  The versions of packages within each version of
- * SuiteSparse are meant to work together.  Combining one packge from one
- * version of SuiteSparse, with another package from another version of
- * SuiteSparse, may or may not work.
- *
- * SuiteSparse Version 3.4.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.1
- *  CSparse        version 2.2.3
- *  CXSparse       version 2.2.3
- *  KLU                    version 1.1.0
- *  BTF                    version 1.1.0
- *  LDL                    version 2.0.1
- *  UFconfig       version number is the same as SuiteSparse
- *  UMFPACK        version 5.4.0
- *  RBio           version 1.1.2
- *  UFcollection    version 1.2.0
- *  LINFACTOR       version 1.1.0
- *  MESHND          version 1.1.1
- *  SSMULT          version 2.0.0
- *  MATLAB_Tools    no specific version number
- *  SuiteSparseQR   version 1.1.2
- *
- * 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 "May 20, 2009"
-#define SUITESPARSE_VER_CODE(main,sub) ((main) * 1000 + (sub))
-#define SUITESPARSE_MAIN_VERSION 3
-#define SUITESPARSE_SUB_VERSION 4
-#define SUITESPARSE_SUBSUB_VERSION 0
-#define SUITESPARSE_VERSION \
-    SUITESPARSE_VER_CODE(SUITESPARSE_MAIN_VERSION,SUITESPARSE_SUB_VERSION)
-
-#ifdef __cplusplus
-}
-#endif
-#endif
diff --git a/extern/colamd/Include/colamd.h b/extern/colamd/Include/colamd.h
deleted file mode 100644 (file)
index 26372d8..0000000
+++ /dev/null
@@ -1,255 +0,0 @@
-/* ========================================================================== */
-/* === colamd/symamd prototypes and definitions ============================= */
-/* ========================================================================== */
-
-/* COLAMD / SYMAMD include file
-
-    You must include this file (colamd.h) in any routine that uses colamd,
-    symamd, or the related macros and definitions.
-
-    Authors:
-
-       The authors of the code itself are Stefan I. Larimore and Timothy A.
-       Davis (davis at cise.ufl.edu), University of Florida.  The algorithm was
-       developed in collaboration with John Gilbert, Xerox PARC, and Esmond
-       Ng, Oak Ridge National Laboratory.
-
-    Acknowledgements:
-
-       This work was supported by the National Science Foundation, under
-       grants DMS-9504974 and DMS-9803599.
-
-    Notice:
-
-       Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
-
-       THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
-       EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
-
-       Permission is hereby granted to use, copy, modify, and/or distribute
-       this program, provided that the Copyright, this License, and the
-       Availability of the original version is retained on all copies and made
-       accessible to the end-user of any code or package that includes COLAMD
-       or any modified version of COLAMD. 
-
-    Availability:
-
-       The colamd/symamd library is available at
-
-           http://www.cise.ufl.edu/research/sparse/colamd/
-
-       This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.h
-       file.  It is required by the colamd.c, colamdmex.c, and symamdmex.c
-       files, and by any C code that calls the routines whose prototypes are
-       listed below, or that uses the colamd/symamd definitions listed below.
-
-*/
-
-#ifndef COLAMD_H
-#define COLAMD_H
-
-/* make it easy for C++ programs to include COLAMD */
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-/* ========================================================================== */
-/* === Include files ======================================================== */
-/* ========================================================================== */
-
-#include <stdlib.h>
-
-/* ========================================================================== */
-/* === COLAMD version ======================================================= */
-/* ========================================================================== */
-
-/* COLAMD Version 2.4 and later will include the following definitions.
- * As an example, to test if the version you are using is 2.4 or later:
- *
- * #ifdef COLAMD_VERSION
- *     if (COLAMD_VERSION >= COLAMD_VERSION_CODE (2,4)) ...
- * #endif
- *
- * This also works during compile-time:
- *
- *  #if defined(COLAMD_VERSION) && (COLAMD_VERSION >= COLAMD_VERSION_CODE (2,4))
- *    printf ("This is version 2.4 or later\n") ;
- *  #else
- *    printf ("This is an early version\n") ;
- *  #endif
- *
- * Versions 2.3 and earlier of COLAMD do not include a #define'd version number.
- */
-
-#define COLAMD_DATE "Nov 1, 2007"
-#define COLAMD_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
-#define COLAMD_MAIN_VERSION 2
-#define COLAMD_SUB_VERSION 7
-#define COLAMD_SUBSUB_VERSION 1
-#define COLAMD_VERSION \
-       COLAMD_VERSION_CODE(COLAMD_MAIN_VERSION,COLAMD_SUB_VERSION)
-
-/* ========================================================================== */
-/* === Knob and statistics definitions ====================================== */
-/* ========================================================================== */
-
-/* size of the knobs [ ] array.  Only knobs [0..1] are currently used. */
-#define COLAMD_KNOBS 20
-
-/* number of output statistics.  Only stats [0..6] are currently used. */
-#define COLAMD_STATS 20
-
-/* knobs [0] and stats [0]: dense row knob and output statistic. */
-#define COLAMD_DENSE_ROW 0
-
-/* knobs [1] and stats [1]: dense column knob and output statistic. */
-#define COLAMD_DENSE_COL 1
-
-/* knobs [2]: aggressive absorption */
-#define COLAMD_AGGRESSIVE 2
-
-/* stats [2]: memory defragmentation count output statistic */
-#define COLAMD_DEFRAG_COUNT 2
-
-/* stats [3]: colamd status:  zero OK, > 0 warning or notice, < 0 error */
-#define COLAMD_STATUS 3
-
-/* stats [4..6]: error info, or info on jumbled columns */ 
-#define COLAMD_INFO1 4
-#define COLAMD_INFO2 5
-#define COLAMD_INFO3 6
-
-/* error codes returned in stats [3]: */
-#define COLAMD_OK                              (0)
-#define COLAMD_OK_BUT_JUMBLED                  (1)
-#define COLAMD_ERROR_A_not_present             (-1)
-#define COLAMD_ERROR_p_not_present             (-2)
-#define COLAMD_ERROR_nrow_negative             (-3)
-#define COLAMD_ERROR_ncol_negative             (-4)
-#define COLAMD_ERROR_nnz_negative              (-5)
-#define COLAMD_ERROR_p0_nonzero                        (-6)
-#define COLAMD_ERROR_A_too_small               (-7)
-#define COLAMD_ERROR_col_length_negative       (-8)
-#define COLAMD_ERROR_row_index_out_of_bounds   (-9)
-#define COLAMD_ERROR_out_of_memory             (-10)
-#define COLAMD_ERROR_internal_error            (-999)
-
-
-/* ========================================================================== */
-/* === Prototypes of user-callable routines ================================= */
-/* ========================================================================== */
-
-/* define UF_long */
-#include "UFconfig.h"
-
-size_t colamd_recommended      /* returns recommended value of Alen, */
-                               /* or 0 if input arguments are erroneous */
-(
-    int nnz,                   /* nonzeros in A */
-    int n_row,                 /* number of rows in A */
-    int n_col                  /* number of columns in A */
-) ;
-
-size_t colamd_l_recommended    /* returns recommended value of Alen, */
-                               /* or 0 if input arguments are erroneous */
-(
-    UF_long nnz,               /* nonzeros in A */
-    UF_long n_row,             /* number of rows in A */
-    UF_long n_col              /* number of columns in A */
-) ;
-
-void colamd_set_defaults       /* sets default parameters */
-(                              /* knobs argument is modified on output */
-    double knobs [COLAMD_KNOBS]        /* parameter settings for colamd */
-) ;
-
-void colamd_l_set_defaults     /* sets default parameters */
-(                              /* knobs argument is modified on output */
-    double knobs [COLAMD_KNOBS]        /* parameter settings for colamd */
-) ;
-
-int colamd                     /* returns (1) if successful, (0) otherwise*/
-(                              /* A and p arguments are modified on output */
-    int n_row,                 /* number of rows in A */
-    int n_col,                 /* number of columns in A */
-    int Alen,                  /* size of the array A */
-    int A [],                  /* row indices of A, of size Alen */
-    int p [],                  /* column pointers of A, of size n_col+1 */
-    double knobs [COLAMD_KNOBS],/* parameter settings for colamd */
-    int stats [COLAMD_STATS]   /* colamd output statistics and error codes */
-) ;
-
-UF_long colamd_l               /* returns (1) if successful, (0) otherwise*/
-(                              /* A and p arguments are modified on output */
-    UF_long n_row,             /* number of rows in A */
-    UF_long n_col,             /* number of columns in A */
-    UF_long Alen,              /* size of the array A */
-    UF_long A [],              /* row indices of A, of size Alen */
-    UF_long p [],              /* column pointers of A, of size n_col+1 */
-    double knobs [COLAMD_KNOBS],/* parameter settings for colamd */
-    UF_long stats [COLAMD_STATS]/* colamd output statistics and error codes */
-) ;
-
-int symamd                             /* return (1) if OK, (0) otherwise */
-(
-    int n,                             /* number of rows and columns of A */
-    int A [],                          /* row indices of A */
-    int p [],                          /* column pointers of A */
-    int perm [],                       /* output permutation, size n_col+1 */
-    double knobs [COLAMD_KNOBS],       /* parameters (uses defaults if NULL) */
-    int stats [COLAMD_STATS],          /* output statistics and error codes */
-    void * (*allocate) (size_t, size_t),
-                                       /* pointer to calloc (ANSI C) or */
-                                       /* mxCalloc (for MATLAB mexFunction) */
-    void (*release) (void *)
-                                       /* pointer to free (ANSI C) or */
-                                       /* mxFree (for MATLAB mexFunction) */
-) ;
-
-UF_long symamd_l                       /* return (1) if OK, (0) otherwise */
-(
-    UF_long n,                         /* number of rows and columns of A */
-    UF_long A [],                      /* row indices of A */
-    UF_long p [],                      /* column pointers of A */
-    UF_long perm [],                   /* output permutation, size n_col+1 */
-    double knobs [COLAMD_KNOBS],       /* parameters (uses defaults if NULL) */
-    UF_long stats [COLAMD_STATS],      /* output statistics and error codes */
-    void * (*allocate) (size_t, size_t),
-                                       /* pointer to calloc (ANSI C) or */
-                                       /* mxCalloc (for MATLAB mexFunction) */
-    void (*release) (void *)
-                                       /* pointer to free (ANSI C) or */
-                                       /* mxFree (for MATLAB mexFunction) */
-) ;
-
-void colamd_report
-(
-    int stats [COLAMD_STATS]
-) ;
-
-void colamd_l_report
-(
-    UF_long stats [COLAMD_STATS]
-) ;
-
-void symamd_report
-(
-    int stats [COLAMD_STATS]
-) ;
-
-void symamd_l_report
-(
-    UF_long stats [COLAMD_STATS]
-) ;
-
-#ifndef EXTERN
-#define EXTERN extern
-#endif
-
-EXTERN int (*colamd_printf) (const char *, ...) ;
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif /* COLAMD_H */
diff --git a/extern/colamd/README.txt b/extern/colamd/README.txt
deleted file mode 100644 (file)
index 5ed81c7..0000000
+++ /dev/null
@@ -1,127 +0,0 @@
-The COLAMD ordering method - Version 2.7
--------------------------------------------------------------------------------
-
-The COLAMD column approximate minimum degree ordering algorithm computes
-a permutation vector P such that the LU factorization of A (:,P)
-tends to be sparser than that of A.  The Cholesky factorization of
-(A (:,P))'*(A (:,P)) will also tend to be sparser than that of A'*A.
-SYMAMD is a symmetric minimum degree ordering method based on COLAMD,
-available as a MATLAB-callable function.  It constructs a matrix M such
-that M'*M has the same pattern as A, and then uses COLAMD to compute a column
-ordering of M.  Colamd and symamd tend to be faster and generate better
-orderings than their MATLAB counterparts, colmmd and symmmd.
-
-To compile and test the colamd m-files and mexFunctions, just unpack the
-COLAMD/ directory from the COLAMD.tar.gz file, and run MATLAB from
-within that directory.  Next, type colamd_test to compile and test colamd
-and symamd.  This will work on any computer with MATLAB (Unix, PC, or Mac).
-Alternatively, type "make" (in Unix) to compile and run a simple example C
-code, without using MATLAB.
-
-To compile and install the colamd m-files and mexFunctions, just cd to
-COLAMD/MATLAB and type colamd_install in the MATLAB command window.
-A short demo will run.  Optionally, type colamd_test to run an extensive tests.
-Type "make" in Unix in the COLAMD directory to compile the C-callable
-library and to run a short demo.
-
-If you have MATLAB 7.2 or earlier, you must first edit UFconfig/UFconfig.h to
-remove the "-largeArrayDims" option from the MEX command (or just use
-colamd_make.m inside MATLAB).
-
-Colamd is a built-in routine in MATLAB, available from The 
-Mathworks, Inc.  Under most cases, the compiled COLAMD from Versions 2.0 to the
-current version do not differ.  Colamd Versions 2.2 and 2.3 differ only in their
-mexFunction interaces to MATLAB.  v2.4 fixes a bug in the symamd routine in
-v2.3.  The bug (in v2.3 and earlier) has no effect on the MATLAB symamd
-mexFunction.  v2.5 adds additional checks for integer overflow, so that
-the "int" version can be safely used with 64-bit pointers.  Refer to the
-ChangeLog for more details.
-
-To use colamd and symamd within an application written in C, all you need are
-colamd.c, colamd_global.c, and colamd.h, which are the C-callable
-colamd/symamd codes.  See colamd.c for more information on how to call
-colamd from a C program.
-
-Requires UFconfig, in the ../UFconfig directory relative to this directory.
-
-       Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
-
-       See http://www.cise.ufl.edu/research/sparse/colamd (the colamd.c
-       file) for the License.
-
-
-Related papers:
-
-       T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, An approximate column
-       minimum degree ordering algorithm, ACM Transactions on Mathematical
-       Software, vol. 30, no. 3., pp. 353-376, 2004.
-
-       T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD,
-       an approximate column minimum degree ordering algorithm, ACM
-       Transactions on Mathematical Software, vol. 30, no. 3., pp. 377-380,
-       2004.
-
-       "An approximate minimum degree column ordering algorithm",
-       S. I. Larimore, MS Thesis, Dept. of Computer and Information
-       Science and Engineering, University of Florida, Gainesville, FL,
-       1998.  CISE Tech Report TR-98-016.  Available at 
-       ftp://ftp.cise.ufl.edu/cis/tech-reports/tr98/tr98-016.ps
-       via anonymous ftp.
-
-       Approximate Deficiency for Ordering the Columns of a Matrix,
-       J. L. Kern, Senior Thesis, Dept. of Computer and Information
-       Science and Engineering, University of Florida, Gainesville, FL,
-       1999.  Available at http://www.cise.ufl.edu/~davis/Kern/kern.ps 
-
-
-Authors:  Stefan I. Larimore and Timothy A. Davis, University of Florida,
-in collaboration with John Gilbert, Xerox PARC (now at UC Santa Barbara),
-and Esmong Ng, Lawrence Berkeley National Laboratory (much of this work
-he did while at Oak Ridge National Laboratory). 
-
-COLAMD files:
-
-    Demo           simple demo
-    Doc                    additional documentation (see colamd.c for more)
-    Include        include file
-    Lib                    compiled C-callable library
-    Makefile       primary Unix Makefile
-    MATLAB         MATLAB functions
-    README.txt     this file
-    Source         C source code
-
-    ./Demo:
-    colamd_example.c       simple example
-    colamd_example.out     output of colamd_example.c
-    colamd_l_example.c     simple example, long integers
-    colamd_l_example.out    output of colamd_l_example.c
-    Makefile               Makefile for C demos
-
-    ./Doc:
-    ChangeLog      change log
-    lesser.txt     license
-
-    ./Include:
-    colamd.h       include file
-
-    ./Lib:
-    Makefile       Makefile for C-callable library
-
-    ./MATLAB:
-    colamd2.m          MATLAB interface for colamd2
-    colamd_demo.m      simple demo
-    colamd_install.m   compile and install colamd2 and symamd2
-    colamd_make.m      compile colamd2 and symamd2
-    colamdmex.ca       MATLAB mexFunction for colamd2
-    colamd_test.m      extensive test
-    colamdtestmex.c    test function for colamd
-    Contents.m         contents of the MATLAB directory
-    luflops.m          test code
-    Makefile           Makefile for MATLAB functions
-    symamd2.m          MATLAB interface for symamd2
-    symamdmex.c                MATLAB mexFunction for symamd2
-    symamdtestmex.c    test function for symamd
-
-    ./Source:
-    colamd.c           primary source code
-    colamd_global.c    globally defined function pointers (malloc, free, ...)
diff --git a/extern/colamd/SConscript b/extern/colamd/SConscript
deleted file mode 100644 (file)
index 7930e3a..0000000
+++ /dev/null
@@ -1,14 +0,0 @@
-#!/usr/bin/python
-import sys
-import os
-
-Import('env')
-
-defs = ''
-cflags = []
-
-src = env.Glob('Source/*.c')
-
-incs = './Include'
-
-env.BlenderLib ( libname = 'extern_colamd', sources=src, includes=Split(incs), defines=Split(defs), libtype=['extern', 'player'], priority=[20,137], compileflags=cflags )
diff --git a/extern/colamd/Source/colamd.c b/extern/colamd/Source/colamd.c
deleted file mode 100644 (file)
index 5fe20d6..0000000
+++ /dev/null
@@ -1,3611 +0,0 @@
-/* ========================================================================== */
-/* === colamd/symamd - a sparse matrix column ordering algorithm ============ */
-/* ========================================================================== */
-
-/* COLAMD / SYMAMD
-
-    colamd:  an approximate minimum degree column ordering algorithm,
-       for LU factorization of symmetric or unsymmetric matrices,
-       QR factorization, least squares, interior point methods for
-       linear programming problems, and other related problems.
-
-    symamd:  an approximate minimum degree ordering algorithm for Cholesky
-       factorization of symmetric matrices.
-
-    Purpose:
-
-       Colamd computes a permutation Q such that the Cholesky factorization of
-       (AQ)'(AQ) has less fill-in and requires fewer floating point operations
-       than A'A.  This also provides a good ordering for sparse partial
-       pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
-       factorization, and P is computed during numerical factorization via
-       conventional partial pivoting with row interchanges.  Colamd is the
-       column ordering method used in SuperLU, part of the ScaLAPACK library.
-       It is also available as built-in function in MATLAB Version 6,
-       available from MathWorks, Inc. (http://www.mathworks.com).  This
-       routine can be used in place of colmmd in MATLAB.
-
-       Symamd computes a permutation P of a symmetric matrix A such that the
-       Cholesky factorization of PAP' has less fill-in and requires fewer
-       floating point operations than A.  Symamd constructs a matrix M such
-       that M'M has the same nonzero pattern of A, and then orders the columns
-       of M using colmmd.  The column ordering of M is then returned as the
-       row and column ordering P of A. 
-
-    Authors:
-
-       The authors of the code itself are Stefan I. Larimore and Timothy A.
-       Davis (davis at cise.ufl.edu), University of Florida.  The algorithm was
-       developed in collaboration with John Gilbert, Xerox PARC, and Esmond
-       Ng, Oak Ridge National Laboratory.
-
-    Acknowledgements:
-
-       This work was supported by the National Science Foundation, under
-       grants DMS-9504974 and DMS-9803599.
-
-    Copyright and License:
-
-       Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
-       COLAMD is also available under alternate licenses, contact T. Davis
-       for details.
-
-       This library 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 library 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 library; if not, write to the Free Software
-       Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
-       USA
-
-       Permission is hereby granted to use or copy this program under the
-       terms of the GNU LGPL, provided that the Copyright, this License,
-       and the Availability of the original version is retained on all copies.
-       User documentation of any code that uses this code or any modified
-       version of this code must cite the Copyright, this License, the
-       Availability note, and "Used by permission." Permission to modify
-       the code and to distribute modified code is granted, provided the
-       Copyright, this License, and the Availability note are retained,
-       and a notice that the code was modified is included.
-
-    Availability:
-
-       The colamd/symamd library is available at
-
-           http://www.cise.ufl.edu/research/sparse/colamd/
-
-       This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c
-       file.  It requires the colamd.h file.  It is required by the colamdmex.c
-       and symamdmex.c files, for the MATLAB interface to colamd and symamd.
-       Appears as ACM Algorithm 836.
-
-    See the ChangeLog file for changes since Version 1.0.
-
-    References:
-
-       T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, An approximate column
-       minimum degree ordering algorithm, ACM Transactions on Mathematical
-       Software, vol. 30, no. 3., pp. 353-376, 2004.
-
-       T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD,
-       an approximate column minimum degree ordering algorithm, ACM
-       Transactions on Mathematical Software, vol. 30, no. 3., pp. 377-380,
-       2004.
-
-*/
-
-/* ========================================================================== */
-/* === Description of user-callable routines ================================ */
-/* ========================================================================== */
-
-/* COLAMD includes both int and UF_long versions of all its routines.  The
- * description below is for the int version.  For UF_long, all int arguments
- * become UF_long.  UF_long is normally defined as long, except for WIN64.
-
-    ----------------------------------------------------------------------------
-    colamd_recommended:
-    ----------------------------------------------------------------------------
-
-       C syntax:
-
-           #include "colamd.h"
-           size_t colamd_recommended (int nnz, int n_row, int n_col) ;
-           size_t colamd_l_recommended (UF_long nnz, UF_long n_row,
-               UF_long n_col) ;
-
-       Purpose:
-
-           Returns recommended value of Alen for use by colamd.  Returns 0
-           if any input argument is negative.  The use of this routine
-           is optional.  Not needed for symamd, which dynamically allocates
-           its own memory.
-
-           Note that in v2.4 and earlier, these routines returned int or long.
-           They now return a value of type size_t.
-
-       Arguments (all input arguments):
-
-           int nnz ;           Number of nonzeros in the matrix A.  This must
-                               be the same value as p [n_col] in the call to
-                               colamd - otherwise you will get a wrong value
-                               of the recommended memory to use.
-
-           int n_row ;         Number of rows in the matrix A.
-
-           int n_col ;         Number of columns in the matrix A.
-
-    ----------------------------------------------------------------------------
-    colamd_set_defaults:
-    ----------------------------------------------------------------------------
-
-       C syntax:
-
-           #include "colamd.h"
-           colamd_set_defaults (double knobs [COLAMD_KNOBS]) ;
-           colamd_l_set_defaults (double knobs [COLAMD_KNOBS]) ;
-
-       Purpose:
-
-           Sets the default parameters.  The use of this routine is optional.
-
-       Arguments:
-
-           double knobs [COLAMD_KNOBS] ;       Output only.
-
-               NOTE: the meaning of the dense row/col knobs has changed in v2.4
-
-               knobs [0] and knobs [1] control dense row and col detection:
-
-               Colamd: rows with more than
-               max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n_col))
-               entries are removed prior to ordering.  Columns with more than
-               max (16, knobs [COLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col)))
-               entries are removed prior to
-               ordering, and placed last in the output column ordering. 
-
-               Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0].
-               Rows and columns with more than
-               max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n))
-               entries are removed prior to ordering, and placed last in the
-               output ordering.
-
-               COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
-               respectively, in colamd.h.  Default values of these two knobs
-               are both 10.  Currently, only knobs [0] and knobs [1] are
-               used, but future versions may use more knobs.  If so, they will
-               be properly set to their defaults by the future version of
-               colamd_set_defaults, so that the code that calls colamd will
-               not need to change, assuming that you either use
-               colamd_set_defaults, or pass a (double *) NULL pointer as the
-               knobs array to colamd or symamd.
-
-           knobs [2]: aggressive absorption
-
-               knobs [COLAMD_AGGRESSIVE] controls whether or not to do
-               aggressive absorption during the ordering.  Default is TRUE.
-
-
-    ----------------------------------------------------------------------------
-    colamd:
-    ----------------------------------------------------------------------------
-
-       C syntax:
-
-           #include "colamd.h"
-           int colamd (int n_row, int n_col, int Alen, int *A, int *p,
-               double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ;
-           UF_long colamd_l (UF_long n_row, UF_long n_col, UF_long Alen,
-               UF_long *A, UF_long *p, double knobs [COLAMD_KNOBS],
-               UF_long stats [COLAMD_STATS]) ;
-
-       Purpose:
-
-           Computes a column ordering (Q) of A such that P(AQ)=LU or
-           (AQ)'AQ=LL' have less fill-in and require fewer floating point
-           operations than factorizing the unpermuted matrix A or A'A,
-           respectively.
-           
-       Returns:
-
-           TRUE (1) if successful, FALSE (0) otherwise.
-
-       Arguments:
-
-           int n_row ;         Input argument.
-
-               Number of rows in the matrix A.
-               Restriction:  n_row >= 0.
-               Colamd returns FALSE if n_row is negative.
-
-           int n_col ;         Input argument.
-
-               Number of columns in the matrix A.
-               Restriction:  n_col >= 0.
-               Colamd returns FALSE if n_col is negative.
-
-           int Alen ;          Input argument.
-
-               Restriction (see note):
-               Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col
-               Colamd returns FALSE if these conditions are not met.
-
-               Note:  this restriction makes an modest assumption regarding
-               the size of the two typedef's structures in colamd.h.
-               We do, however, guarantee that
-
-                       Alen >= colamd_recommended (nnz, n_row, n_col)
-
-               will be sufficient.  Note: the macro version does not check
-               for integer overflow, and thus is not recommended.  Use
-               the colamd_recommended routine instead.
-
-           int A [Alen] ;      Input argument, undefined on output.
-
-               A is an integer array of size Alen.  Alen must be at least as
-               large as the bare minimum value given above, but this is very
-               low, and can result in excessive run time.  For best
-               performance, we recommend that Alen be greater than or equal to
-               colamd_recommended (nnz, n_row, n_col), which adds
-               nnz/5 to the bare minimum value given above.
-
-               On input, the row indices of the entries in column c of the
-               matrix are held in A [(p [c]) ... (p [c+1]-1)].  The row indices
-               in a given column c need not be in ascending order, and
-               duplicate row indices may be be present.  However, colamd will
-               work a little faster if both of these conditions are met
-               (Colamd puts the matrix into this format, if it finds that the
-               the conditions are not met).
-
-               The matrix is 0-based.  That is, rows are in the range 0 to
-               n_row-1, and columns are in the range 0 to n_col-1.  Colamd
-               returns FALSE if any row index is out of range.
-
-               The contents of A are modified during ordering, and are
-               undefined on output.
-
-           int p [n_col+1] ;   Both input and output argument.
-
-               p is an integer array of size n_col+1.  On input, it holds the
-               "pointers" for the column form of the matrix A.  Column c of
-               the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
-               entry, p [0], must be zero, and p [c] <= p [c+1] must hold
-               for all c in the range 0 to n_col-1.  The value p [n_col] is
-               thus the total number of entries in the pattern of the matrix A.
-               Colamd returns FALSE if these conditions are not met.
-
-               On output, if colamd returns TRUE, the array p holds the column
-               permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
-               the first column index in the new ordering, and p [n_col-1] is
-               the last.  That is, p [k] = j means that column j of A is the
-               kth pivot column, in AQ, where k is in the range 0 to n_col-1
-               (p [0] = j means that column j of A is the first column in AQ).
-
-               If colamd returns FALSE, then no permutation is returned, and
-               p is undefined on output.
-
-           double knobs [COLAMD_KNOBS] ;       Input argument.
-
-               See colamd_set_defaults for a description.
-
-           int stats [COLAMD_STATS] ;          Output argument.
-
-               Statistics on the ordering, and error status.
-               See colamd.h for related definitions.
-               Colamd returns FALSE if stats is not present.
-
-               stats [0]:  number of dense or empty rows ignored.
-
-               stats [1]:  number of dense or empty columns ignored (and
-                               ordered last in the output permutation p)
-                               Note that a row can become "empty" if it
-                               contains only "dense" and/or "empty" columns,
-                               and similarly a column can become "empty" if it
-                               only contains "dense" and/or "empty" rows.
-
-               stats [2]:  number of garbage collections performed.
-                               This can be excessively high if Alen is close
-                               to the minimum required value.
-
-               stats [3]:  status code.  < 0 is an error code.
-                           > 1 is a warning or notice.
-
-                       0       OK.  Each column of the input matrix contained
-                               row indices in increasing order, with no
-                               duplicates.
-
-                       1       OK, but columns of input matrix were jumbled
-                               (unsorted columns or duplicate entries).  Colamd
-                               had to do some extra work to sort the matrix
-                               first and remove duplicate entries, but it
-                               still was able to return a valid permutation
-                               (return value of colamd was TRUE).
-
-                                       stats [4]: highest numbered column that
-                                               is unsorted or has duplicate
-                                               entries.
-                                       stats [5]: last seen duplicate or
-                                               unsorted row index.
-                                       stats [6]: number of duplicate or
-                                               unsorted row indices.
-
-                       -1      A is a null pointer
-
-                       -2      p is a null pointer
-
-                       -3      n_row is negative
-
-                                       stats [4]: n_row
-
-                       -4      n_col is negative
-
-                                       stats [4]: n_col
-
-                       -5      number of nonzeros in matrix is negative
-
-                                       stats [4]: number of nonzeros, p [n_col]
-
-                       -6      p [0] is nonzero
-
-                                       stats [4]: p [0]
-
-                       -7      A is too small
-
-                                       stats [4]: required size
-                                       stats [5]: actual size (Alen)
-
-                       -8      a column has a negative number of entries
-
-                                       stats [4]: column with < 0 entries
-                                       stats [5]: number of entries in col
-
-                       -9      a row index is out of bounds
-
-                                       stats [4]: column with bad row index
-                                       stats [5]: bad row index
-                                       stats [6]: n_row, # of rows of matrx
-
-                       -10     (unused; see symamd.c)
-
-                       -999    (unused; see symamd.c)
-
-               Future versions may return more statistics in the stats array.
-
-       Example:
-       
-           See http://www.cise.ufl.edu/research/sparse/colamd/example.c
-           for a complete example.
-
-           To order the columns of a 5-by-4 matrix with 11 nonzero entries in
-           the following nonzero pattern
-
-               x 0 x 0
-               x 0 x x
-               0 x x 0
-               0 0 x x
-               x x 0 0
-
-           with default knobs and no output statistics, do the following:
-
-               #include "colamd.h"
-               #define ALEN 100
-               int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;
-               int p [ ] = {0, 3, 5, 9, 11} ;
-               int stats [COLAMD_STATS] ;
-               colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ;
-
-           The permutation is returned in the array p, and A is destroyed.
-
-    ----------------------------------------------------------------------------
-    symamd:
-    ----------------------------------------------------------------------------
-
-       C syntax:
-
-           #include "colamd.h"
-           int symamd (int n, int *A, int *p, int *perm,
-               double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS],
-               void (*allocate) (size_t, size_t), void (*release) (void *)) ;
-           UF_long symamd_l (UF_long n, UF_long *A, UF_long *p, UF_long *perm,
-               double knobs [COLAMD_KNOBS], UF_long stats [COLAMD_STATS],
-               void (*allocate) (size_t, size_t), void (*release) (void *)) ;
-
-       Purpose:
-
-           The symamd routine computes an ordering P of a symmetric sparse
-           matrix A such that the Cholesky factorization PAP' = LL' remains
-           sparse.  It is based on a column ordering of a matrix M constructed
-           so that the nonzero pattern of M'M is the same as A.  The matrix A
-           is assumed to be symmetric; only the strictly lower triangular part
-           is accessed.  You must pass your selected memory allocator (usually
-           calloc/free or mxCalloc/mxFree) to symamd, for it to allocate
-           memory for the temporary matrix M.
-
-       Returns:
-
-           TRUE (1) if successful, FALSE (0) otherwise.
-
-       Arguments:
-
-           int n ;             Input argument.
-
-               Number of rows and columns in the symmetrix matrix A.
-               Restriction:  n >= 0.
-               Symamd returns FALSE if n is negative.
-
-           int A [nnz] ;       Input argument.
-
-               A is an integer array of size nnz, where nnz = p [n].
-               
-               The row indices of the entries in column c of the matrix are
-               held in A [(p [c]) ... (p [c+1]-1)].  The row indices in a
-               given column c need not be in ascending order, and duplicate
-               row indices may be present.  However, symamd will run faster
-               if the columns are in sorted order with no duplicate entries. 
-
-               The matrix is 0-based.  That is, rows are in the range 0 to
-               n-1, and columns are in the range 0 to n-1.  Symamd
-               returns FALSE if any row index is out of range.
-
-               The contents of A are not modified.
-
-           int p [n+1] ;       Input argument.
-
-               p is an integer array of size n+1.  On input, it holds the
-               "pointers" for the column form of the matrix A.  Column c of
-               the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
-               entry, p [0], must be zero, and p [c] <= p [c+1] must hold
-               for all c in the range 0 to n-1.  The value p [n] is
-               thus the total number of entries in the pattern of the matrix A.
-               Symamd returns FALSE if these conditions are not met.
-
-               The contents of p are not modified.
-
-           int perm [n+1] ;    Output argument.
-
-               On output, if symamd returns TRUE, the array perm holds the
-               permutation P, where perm [0] is the first index in the new
-               ordering, and perm [n-1] is the last.  That is, perm [k] = j
-               means that row and column j of A is the kth column in PAP',
-               where k is in the range 0 to n-1 (perm [0] = j means
-               that row and column j of A are the first row and column in
-               PAP').  The array is used as a workspace during the ordering,
-               which is why it must be of length n+1, not just n.
-
-           double knobs [COLAMD_KNOBS] ;       Input argument.
-
-               See colamd_set_defaults for a description.
-
-           int stats [COLAMD_STATS] ;          Output argument.
-
-               Statistics on the ordering, and error status.
-               See colamd.h for related definitions.
-               Symamd returns FALSE if stats is not present.
-
-               stats [0]:  number of dense or empty row and columns ignored
-                               (and ordered last in the output permutation 
-                               perm).  Note that a row/column can become
-                               "empty" if it contains only "dense" and/or
-                               "empty" columns/rows.
-
-               stats [1]:  (same as stats [0])
-
-               stats [2]:  number of garbage collections performed.
-
-               stats [3]:  status code.  < 0 is an error code.
-                           > 1 is a warning or notice.
-
-                       0       OK.  Each column of the input matrix contained
-                               row indices in increasing order, with no
-                               duplicates.
-
-                       1       OK, but columns of input matrix were jumbled
-                               (unsorted columns or duplicate entries).  Symamd
-                               had to do some extra work to sort the matrix
-                               first and remove duplicate entries, but it
-                               still was able to return a valid permutation
-                               (return value of symamd was TRUE).
-
-                                       stats [4]: highest numbered column that
-                                               is unsorted or has duplicate
-                                               entries.
-                                       stats [5]: last seen duplicate or
-                                               unsorted row index.
-                                       stats [6]: number of duplicate or
-                                               unsorted row indices.
-
-                       -1      A is a null pointer
-
-                       -2      p is a null pointer
-
-                       -3      (unused, see colamd.c)
-
-                       -4      n is negative
-
-                                       stats [4]: n
-
-                       -5      number of nonzeros in matrix is negative
-
-                                       stats [4]: # of nonzeros (p [n]).
-
-                       -6      p [0] is nonzero
-
-                                       stats [4]: p [0]
-
-                       -7      (unused)
-
-                       -8      a column has a negative number of entries
-
-                                       stats [4]: column with < 0 entries
-                                       stats [5]: number of entries in col
-
-                       -9      a row index is out of bounds
-
-                                       stats [4]: column with bad row index
-                                       stats [5]: bad row index
-                                       stats [6]: n_row, # of rows of matrx
-
-                       -10     out of memory (unable to allocate temporary
-                               workspace for M or count arrays using the
-                               "allocate" routine passed into symamd).
-
-               Future versions may return more statistics in the stats array.
-
-           void * (*allocate) (size_t, size_t)
-
-               A pointer to a function providing memory allocation.  The
-               allocated memory must be returned initialized to zero.  For a
-               C application, this argument should normally be a pointer to
-               calloc.  For a MATLAB mexFunction, the routine mxCalloc is
-               passed instead.
-
-           void (*release) (size_t, size_t)
-
-               A pointer to a function that frees memory allocated by the
-               memory allocation routine above.  For a C application, this
-               argument should normally be a pointer to free.  For a MATLAB
-               mexFunction, the routine mxFree is passed instead.
-
-
-    ----------------------------------------------------------------------------
-    colamd_report:
-    ----------------------------------------------------------------------------
-
-       C syntax:
-
-           #include "colamd.h"
-           colamd_report (int stats [COLAMD_STATS]) ;
-           colamd_l_report (UF_long stats [COLAMD_STATS]) ;
-
-       Purpose:
-
-           Prints the error status and statistics recorded in the stats
-           array on the standard error output (for a standard C routine)
-           or on the MATLAB output (for a mexFunction).
-
-       Arguments:
-
-           int stats [COLAMD_STATS] ;  Input only.  Statistics from colamd.
-
-
-    ----------------------------------------------------------------------------
-    symamd_report:
-    ----------------------------------------------------------------------------
-
-       C syntax:
-
-           #include "colamd.h"
-           symamd_report (int stats [COLAMD_STATS]) ;
-           symamd_l_report (UF_long stats [COLAMD_STATS]) ;
-
-       Purpose:
-
-           Prints the error status and statistics recorded in the stats
-           array on the standard error output (for a standard C routine)
-           or on the MATLAB output (for a mexFunction).
-
-       Arguments:
-
-           int stats [COLAMD_STATS] ;  Input only.  Statistics from symamd.
-
-
-*/
-
-/* ========================================================================== */
-/* === Scaffolding code definitions  ======================================== */
-/* ========================================================================== */
-
-/* Ensure that debugging is turned off: */
-#ifndef NDEBUG
-#define NDEBUG
-#endif
-
-/* turn on debugging by uncommenting the following line
- #undef NDEBUG
-*/
-
-/*
-   Our "scaffolding code" philosophy:  In our opinion, well-written library
-   code should keep its "debugging" code, and just normally have it turned off
-   by the compiler so as not to interfere with performance.  This serves
-   several purposes:
-
-   (1) assertions act as comments to the reader, telling you what the code
-       expects at that point.  All assertions will always be true (unless
-       there really is a bug, of course).
-
-   (2) leaving in the scaffolding code assists anyone who would like to modify
-       the code, or understand the algorithm (by reading the debugging output,
-       one can get a glimpse into what the code is doing).
-
-   (3) (gasp!) for actually finding bugs.  This code has been heavily tested
-       and "should" be fully functional and bug-free ... but you never know...
-
-    The code will become outrageously slow when debugging is
-    enabled.  To control the level of debugging output, set an environment
-    variable D to 0 (little), 1 (some), 2, 3, or 4 (lots).  When debugging,
-    you should see the following message on the standard output:
-
-       colamd: debug version, D = 1 (THIS WILL BE SLOW!)
-
-    or a similar message for symamd.  If you don't, then debugging has not
-    been enabled.
-
-*/
-
-/* ========================================================================== */
-/* === Include files ======================================================== */
-/* ========================================================================== */
-
-#include "colamd.h"
-#include <limits.h>
-#include <math.h>
-
-#ifdef MATLAB_MEX_FILE
-#include "mex.h"
-#include "matrix.h"
-#endif /* MATLAB_MEX_FILE */
-
-#if !defined (NPRINT) || !defined (NDEBUG)
-#include <stdio.h>
-#endif
-
-#ifndef NULL
-#define NULL ((void *) 0)
-#endif
-
-/* ========================================================================== */
-/* === int or UF_long ======================================================= */
-/* ========================================================================== */
-
-/* define UF_long */
-#include "UFconfig.h"
-
-#ifdef DLONG
-
-#define Int UF_long
-#define ID  UF_long_id
-#define Int_MAX UF_long_max
-
-#define COLAMD_recommended colamd_l_recommended
-#define COLAMD_set_defaults colamd_l_set_defaults
-#define COLAMD_MAIN colamd_l
-#define SYMAMD_MAIN symamd_l
-#define COLAMD_report colamd_l_report
-#define SYMAMD_report symamd_l_report
-
-#else
-
-#define Int int
-#define ID "%d"
-#define Int_MAX INT_MAX
-
-#define COLAMD_recommended colamd_recommended
-#define COLAMD_set_defaults colamd_set_defaults
-#define COLAMD_MAIN colamd
-#define SYMAMD_MAIN symamd
-#define COLAMD_report colamd_report
-#define SYMAMD_report symamd_report
-
-#endif
-
-/* ========================================================================== */
-/* === Row and Column structures ============================================ */
-/* ========================================================================== */
-
-/* User code that makes use of the colamd/symamd routines need not directly */
-/* reference these structures.  They are used only for colamd_recommended. */
-
-typedef struct Colamd_Col_struct
-{
-    Int start ;                /* index for A of first row in this column, or DEAD */
-                       /* if column is dead */
-    Int length ;       /* number of rows in this column */
-    union
-    {
-       Int thickness ; /* number of original columns represented by this */
-                       /* col, if the column is alive */
-       Int parent ;    /* parent in parent tree super-column structure, if */
-                       /* the column is dead */
-    } shared1 ;
-    union
-    {
-       Int score ;     /* the score used to maintain heap, if col is alive */
-       Int order ;     /* pivot ordering of this column, if col is dead */
-    } shared2 ;
-    union
-    {
-       Int headhash ;  /* head of a hash bucket, if col is at the head of */
-                       /* a degree list */
-       Int hash ;      /* hash value, if col is not in a degree list */
-       Int prev ;      /* previous column in degree list, if col is in a */
-                       /* degree list (but not at the head of a degree list) */
-    } shared3 ;
-    union
-    {
-       Int degree_next ;       /* next column, if col is in a degree list */
-       Int hash_next ;         /* next column, if col is in a hash list */
-    } shared4 ;
-
-} Colamd_Col ;
-
-typedef struct Colamd_Row_struct
-{
-    Int start ;                /* index for A of first col in this row */
-    Int length ;       /* number of principal columns in this row */
-    union
-    {
-       Int degree ;    /* number of principal & non-principal columns in row */
-       Int p ;         /* used as a row pointer in init_rows_cols () */
-    } shared1 ;
-    union
-    {
-       Int mark ;      /* for computing set differences and marking dead rows*/
-       Int first_column ;/* first column in row (used in garbage collection) */
-    } shared2 ;
-
-} Colamd_Row ;
-
-/* ========================================================================== */
-/* === Definitions ========================================================== */
-/* ========================================================================== */
-
-/* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
-#define PUBLIC
-#define PRIVATE static
-
-#define DENSE_DEGREE(alpha,n) \
-    ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))
-
-#define MAX(a,b) (((a) > (b)) ? (a) : (b))
-#define MIN(a,b) (((a) < (b)) ? (a) : (b))
-
-#define ONES_COMPLEMENT(r) (-(r)-1)
-
-/* -------------------------------------------------------------------------- */
-/* Change for version 2.1:  define TRUE and FALSE only if not yet defined */  
-/* -------------------------------------------------------------------------- */
-
-#ifndef TRUE
-#define TRUE (1)
-#endif
-
-#ifndef FALSE
-#define FALSE (0)
-#endif
-
-/* -------------------------------------------------------------------------- */
-
-#define EMPTY  (-1)
-
-/* Row and column status */
-#define ALIVE  (0)
-#define DEAD   (-1)
-
-/* Column status */
-#define DEAD_PRINCIPAL         (-1)
-#define DEAD_NON_PRINCIPAL     (-2)
-
-/* Macros for row and column status update and checking. */
-#define ROW_IS_DEAD(r)                 ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
-#define ROW_IS_MARKED_DEAD(row_mark)   (row_mark < ALIVE)
-#define ROW_IS_ALIVE(r)                        (Row [r].shared2.mark >= ALIVE)
-#define COL_IS_DEAD(c)                 (Col [c].start < ALIVE)
-#define COL_IS_ALIVE(c)                        (Col [c].start >= ALIVE)
-#define COL_IS_DEAD_PRINCIPAL(c)       (Col [c].start == DEAD_PRINCIPAL)
-#define KILL_ROW(r)                    { Row [r].shared2.mark = DEAD ; }
-#define KILL_PRINCIPAL_COL(c)          { Col [c].start = DEAD_PRINCIPAL ; }
-#define KILL_NON_PRINCIPAL_COL(c)      { Col [c].start = DEAD_NON_PRINCIPAL ; }
-
-/* ========================================================================== */
-/* === Colamd reporting mechanism =========================================== */
-/* ========================================================================== */
-
-#if defined (MATLAB_MEX_FILE) || defined (MATHWORKS)
-/* In MATLAB, matrices are 1-based to the user, but 0-based internally */
-#define INDEX(i) ((i)+1)
-#else
-/* In C, matrices are 0-based and indices are reported as such in *_report */
-#define INDEX(i) (i)
-#endif
-
-/* All output goes through the PRINTF macro.  */
-#define PRINTF(params) { if (colamd_printf != NULL) (void) colamd_printf params ; }
-
-/* ========================================================================== */
-/* === Prototypes of PRIVATE routines ======================================= */
-/* ========================================================================== */
-
-PRIVATE Int init_rows_cols
-(
-    Int n_row,
-    Int n_col,
-    Colamd_Row Row [],
-    Colamd_Col Col [],
-    Int A [],
-    Int p [],
-    Int stats [COLAMD_STATS]
-) ;
-
-PRIVATE void init_scoring
-(
-    Int n_row,
-    Int n_col,
-    Colamd_Row Row [],
-    Colamd_Col Col [],
-    Int A [],
-    Int head [],
-    double knobs [COLAMD_KNOBS],
-    Int *p_n_row2,
-    Int *p_n_col2,
-    Int *p_max_deg
-) ;
-
-PRIVATE Int find_ordering
-(
-    Int n_row,
-    Int n_col,
-    Int Alen,
-    Colamd_Row Row [],
-    Colamd_Col Col [],
-    Int A [],
-    Int head [],
-    Int n_col2,
-    Int max_deg,
-    Int pfree,
-    Int aggressive
-) ;
-
-PRIVATE void order_children
-(
-    Int n_col,
-    Colamd_Col Col [],
-    Int p []
-) ;
-
-PRIVATE void detect_super_cols
-(
-
-#ifndef NDEBUG
-    Int n_col,
-    Colamd_Row Row [],
-#endif /* NDEBUG */
-
-    Colamd_Col Col [],
-    Int A [],
-    Int head [],
-    Int row_start,
-    Int row_length
-) ;
-
-PRIVATE Int garbage_collection
-(
-    Int n_row,
-    Int n_col,
-    Colamd_Row Row [],
-    Colamd_Col Col [],
-    Int A [],
-    Int *pfree
-) ;
-
-PRIVATE Int clear_mark
-(
-    Int tag_mark,
-    Int max_mark,
-    Int n_row,
-    Colamd_Row Row []
-) ;
-
-PRIVATE void print_report
-(
-    char *method,
-    Int stats [COLAMD_STATS]
-) ;
-
-/* ========================================================================== */
-/* === Debugging prototypes and definitions ================================= */
-/* ========================================================================== */
-
-#ifndef NDEBUG
-
-#include <assert.h>
-
-/* colamd_debug is the *ONLY* global variable, and is only */
-/* present when debugging */
-
-PRIVATE Int colamd_debug = 0 ; /* debug print level */
-
-#define DEBUG0(params) { PRINTF (params) ; }
-#define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; }
-#define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; }
-#define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; }
-#define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; }
-
-#ifdef MATLAB_MEX_FILE
-#define ASSERT(expression) (mxAssert ((expression), ""))
-#else
-#define ASSERT(expression) (assert (expression))
-#endif /* MATLAB_MEX_FILE */
-
-PRIVATE void colamd_get_debug  /* gets the debug print level from getenv */
-(
-    char *method
-) ;
-
-PRIVATE void debug_deg_lists
-(
-    Int n_row,
-    Int n_col,
-    Colamd_Row Row [],
-    Colamd_Col Col [],
-    Int head [],
-    Int min_score,
-    Int should,
-    Int max_deg
-) ;
-
-PRIVATE void debug_mark
-(
-    Int n_row,
-    Colamd_Row Row [],
-    Int tag_mark,
-    Int max_mark
-) ;
-
-PRIVATE void debug_matrix
-(
-    Int n_row,
-    Int n_col,
-    Colamd_Row Row [],
-    Colamd_Col Col [],
-    Int A []
-) ;
-
-PRIVATE void debug_structures
-(
-    Int n_row,
-    Int n_col,
-    Colamd_Row Row [],
-    Colamd_Col Col [],
-    Int A [],
-    Int n_col2
-) ;
-
-#else /* NDEBUG */
-
-/* === No debugging ========================================================= */
-
-#define DEBUG0(params) ;
-#define DEBUG1(params) ;
-#define DEBUG2(params) ;
-#define DEBUG3(params) ;
-#define DEBUG4(params) ;
-
-#define ASSERT(expression)
-
-#endif /* NDEBUG */
-
-/* ========================================================================== */
-/* === USER-CALLABLE ROUTINES: ============================================== */
-/* ========================================================================== */
-
-/* ========================================================================== */
-/* === colamd_recommended =================================================== */
-/* ========================================================================== */
-
-/*
-    The colamd_recommended routine returns the suggested size for Alen.  This
-    value has been determined to provide good balance between the number of
-    garbage collections and the memory requirements for colamd.  If any
-    argument is negative, or if integer overflow occurs, a 0 is returned as an
-    error condition.  2*nnz space is required for the row and column
-    indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is
-    required for the Col and Row arrays, respectively, which are internal to
-    colamd (roughly 6*n_col + 4*n_row).  An additional n_col space is the
-    minimal amount of "elbow room", and nnz/5 more space is recommended for
-    run time efficiency.
-
-    Alen is approximately 2.2*nnz + 7*n_col + 4*n_row + 10.
-
-    This function is not needed when using symamd.
-*/
-
-/* add two values of type size_t, and check for integer overflow */
-static size_t t_add (size_t a, size_t b, int *ok)
-{
-    (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ;
-    return ((*ok) ? (a + b) : 0) ;
-}
-
-/* compute a*k where k is a small integer, and check for integer overflow */
-static size_t t_mult (size_t a, size_t k, int *ok)
-{
-    size_t i, s = 0 ;
-    for (i = 0 ; i < k ; i++)
-    {
-       s = t_add (s, a, ok) ;
-    }
-    return (s) ;
-}
-
-/* size of the Col and Row structures */
-#define COLAMD_C(n_col,ok) \
-    ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int)))
-
-#define COLAMD_R(n_row,ok) \
-    ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int)))
-
-
-PUBLIC size_t COLAMD_recommended       /* returns recommended value of Alen. */
-(
-    /* === Parameters ======================================================= */
-
-    Int nnz,                   /* number of nonzeros in A */
-    Int n_row,                 /* number of rows in A */
-    Int n_col                  /* number of columns in A */
-)
-{
-    size_t s, c, r ;
-    int ok = TRUE ;
-    if (nnz < 0 || n_row < 0 || n_col < 0)
-    {
-       return (0) ;
-    }
-    s = t_mult (nnz, 2, &ok) ;     /* 2*nnz */
-    c = COLAMD_C (n_col, &ok) ;            /* size of column structures */
-    r = COLAMD_R (n_row, &ok) ;            /* size of row structures */
-    s = t_add (s, c, &ok) ;
-    s = t_add (s, r, &ok) ;
-    s = t_add (s, n_col, &ok) ;            /* elbow room */
-    s = t_add (s, nnz/5, &ok) ;            /* elbow room */
-    ok = ok && (s < Int_MAX) ;
-    return (ok ? s : 0) ;
-}
-
-
-/* ========================================================================== */
-/* === colamd_set_defaults ================================================== */
-/* ========================================================================== */
-
-/*
-    The colamd_set_defaults routine sets the default values of the user-
-    controllable parameters for colamd and symamd:
-
-       Colamd: rows with more than max (16, knobs [0] * sqrt (n_col))
-       entries are removed prior to ordering.  Columns with more than
-       max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed
-       prior to ordering, and placed last in the output column ordering. 
-
-       Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n))
-       entries are removed prior to ordering, and placed last in the
-       output ordering.
-
-       knobs [0]       dense row control
-
-       knobs [1]       dense column control
-
-       knobs [2]       if nonzero, do aggresive absorption
-
-       knobs [3..19]   unused, but future versions might use this
-
-*/
-
-PUBLIC void COLAMD_set_defaults
-(
-    /* === Parameters ======================================================= */
-
-    double knobs [COLAMD_KNOBS]                /* knob array */
-)
-{
-    /* === Local variables ================================================== */
-
-    Int i ;
-
-    if (!knobs)
-    {
-       return ;                        /* no knobs to initialize */
-    }
-    for (i = 0 ; i < COLAMD_KNOBS ; i++)
-    {
-       knobs [i] = 0 ;
-    }
-    knobs [COLAMD_DENSE_ROW] = 10 ;
-    knobs [COLAMD_DENSE_COL] = 10 ;
-    knobs [COLAMD_AGGRESSIVE] = TRUE ; /* default: do aggressive absorption*/
-}
-
-
-/* ========================================================================== */
-/* === symamd =============================================================== */
-/* ========================================================================== */
-
-PUBLIC Int SYMAMD_MAIN                 /* return TRUE if OK, FALSE otherwise */
-(
-    /* === Parameters ======================================================= */
-
-    Int n,                             /* number of rows and columns of A */
-    Int A [],                          /* row indices of A */
-    Int p [],                          /* column pointers of A */
-    Int perm [],                       /* output permutation, size n+1 */
-    double knobs [COLAMD_KNOBS],       /* parameters (uses defaults if NULL) */
-    Int stats [COLAMD_STATS],          /* output statistics and error codes */
-    void * (*allocate) (size_t, size_t),
-                                       /* pointer to calloc (ANSI C) or */
-                                       /* mxCalloc (for MATLAB mexFunction) */
-    void (*release) (void *)
-                                       /* pointer to free (ANSI C) or */
-                                       /* mxFree (for MATLAB mexFunction) */
-)
-{
-    /* === Local variables ================================================== */
-
-    Int *count ;               /* length of each column of M, and col pointer*/
-    Int *mark ;                        /* mark array for finding duplicate entries */
-    Int *M ;                   /* row indices of matrix M */
-    size_t Mlen ;              /* length of M */
-    Int n_row ;                        /* number of rows in M */
-    Int nnz ;                  /* number of entries in A */
-    Int i ;                    /* row index of A */
-    Int j ;                    /* column index of A */
-    Int k ;                    /* row index of M */ 
-    Int mnz ;                  /* number of nonzeros in M */
-    Int pp ;                   /* index into a column of A */
-    Int last_row ;             /* last row seen in the current column */
-    Int length ;               /* number of nonzeros in a column */
-
-    double cknobs [COLAMD_KNOBS] ;             /* knobs for colamd */
-    double default_knobs [COLAMD_KNOBS] ;      /* default knobs for colamd */
-
-#ifndef NDEBUG
-    colamd_get_debug ("symamd") ;
-#endif /* NDEBUG */
-
-    /* === Check the input arguments ======================================== */
-
-    if (!stats)
-    {
-       DEBUG0 (("symamd: stats not present\n")) ;
-       return (FALSE) ;
-    }
-    for (i = 0 ; i < COLAMD_STATS ; i++)
-    {
-       stats [i] = 0 ;
-    }
-    stats [COLAMD_STATUS] = COLAMD_OK ;
-    stats [COLAMD_INFO1] = -1 ;
-    stats [COLAMD_INFO2] = -1 ;
-
-    if (!A)
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
-       DEBUG0 (("symamd: A not present\n")) ;
-       return (FALSE) ;
-    }
-
-    if (!p)            /* p is not present */
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
-       DEBUG0 (("symamd: p not present\n")) ;
-       return (FALSE) ;
-    }
-
-    if (n < 0)         /* n must be >= 0 */
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
-       stats [COLAMD_INFO1] = n ;
-       DEBUG0 (("symamd: n negative %d\n", n)) ;
-       return (FALSE) ;
-    }
-
-    nnz = p [n] ;
-    if (nnz < 0)       /* nnz must be >= 0 */
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
-       stats [COLAMD_INFO1] = nnz ;
-       DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ;
-       return (FALSE) ;
-    }
-
-    if (p [0] != 0)
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
-       stats [COLAMD_INFO1] = p [0] ;
-       DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ;
-       return (FALSE) ;
-    }
-
-    /* === If no knobs, set default knobs =================================== */
-
-    if (!knobs)
-    {
-       COLAMD_set_defaults (default_knobs) ;
-       knobs = default_knobs ;
-    }
-
-    /* === Allocate count and mark ========================================== */
-
-    count = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
-    if (!count)
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
-       DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ;
-       return (FALSE) ;
-    }
-
-    mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
-    if (!mark)
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
-       (*release) ((void *) count) ;
-       DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ;
-       return (FALSE) ;
-    }
-
-    /* === Compute column counts of M, check if A is valid ================== */
-
-    stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
-
-    for (i = 0 ; i < n ; i++)
-    {
-       mark [i] = -1 ;
-    }
-
-    for (j = 0 ; j < n ; j++)
-    {
-       last_row = -1 ;
-
-       length = p [j+1] - p [j] ;
-       if (length < 0)
-       {
-           /* column pointers must be non-decreasing */
-           stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
-           stats [COLAMD_INFO1] = j ;
-           stats [COLAMD_INFO2] = length ;
-           (*release) ((void *) count) ;
-           (*release) ((void *) mark) ;
-           DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ;
-           return (FALSE) ;
-       }
-
-       for (pp = p [j] ; pp < p [j+1] ; pp++)
-       {
-           i = A [pp] ;
-           if (i < 0 || i >= n)
-           {
-               /* row index i, in column j, is out of bounds */
-               stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
-               stats [COLAMD_INFO1] = j ;
-               stats [COLAMD_INFO2] = i ;
-               stats [COLAMD_INFO3] = n ;
-               (*release) ((void *) count) ;
-               (*release) ((void *) mark) ;
-               DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ;
-               return (FALSE) ;
-           }
-
-           if (i <= last_row || mark [i] == j)
-           {
-               /* row index is unsorted or repeated (or both), thus col */
-               /* is jumbled.  This is a notice, not an error condition. */
-               stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
-               stats [COLAMD_INFO1] = j ;
-               stats [COLAMD_INFO2] = i ;
-               (stats [COLAMD_INFO3]) ++ ;
-               DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ;
-           }
-
-           if (i > j && mark [i] != j)
-           {
-               /* row k of M will contain column indices i and j */
-               count [i]++ ;
-               count [j]++ ;
-           }
-
-           /* mark the row as having been seen in this column */
-           mark [i] = j ;
-
-           last_row = i ;
-       }
-    }
-
-    /* v2.4: removed free(mark) */
-
-    /* === Compute column pointers of M ===================================== */
-
-    /* use output permutation, perm, for column pointers of M */
-    perm [0] = 0 ;
-    for (j = 1 ; j <= n ; j++)
-    {
-       perm [j] = perm [j-1] + count [j-1] ;
-    }
-    for (j = 0 ; j < n ; j++)
-    {
-       count [j] = perm [j] ;
-    }
-
-    /* === Construct M ====================================================== */
-
-    mnz = perm [n] ;
-    n_row = mnz / 2 ;
-    Mlen = COLAMD_recommended (mnz, n_row, n) ;
-    M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ;
-    DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %g\n",
-       n_row, n, mnz, (double) Mlen)) ;
-
-    if (!M)
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
-       (*release) ((void *) count) ;
-       (*release) ((void *) mark) ;
-       DEBUG0 (("symamd: allocate M (size %g) failed\n", (double) Mlen)) ;
-       return (FALSE) ;
-    }
-
-    k = 0 ;
-
-    if (stats [COLAMD_STATUS] == COLAMD_OK)
-    {
-       /* Matrix is OK */
-       for (j = 0 ; j < n ; j++)
-       {
-           ASSERT (p [j+1] - p [j] >= 0) ;
-           for (pp = p [j] ; pp < p [j+1] ; pp++)
-           {
-               i = A [pp] ;
-               ASSERT (i >= 0 && i < n) ;
-               if (i > j)
-               {
-                   /* row k of M contains column indices i and j */
-                   M [count [i]++] = k ;
-                   M [count [j]++] = k ;
-                   k++ ;
-               }
-           }
-       }
-    }
-    else
-    {
-       /* Matrix is jumbled.  Do not add duplicates to M.  Unsorted cols OK. */
-       DEBUG0 (("symamd: Duplicates in A.\n")) ;
-       for (i = 0 ; i < n ; i++)
-       {
-           mark [i] = -1 ;
-       }
-       for (j = 0 ; j < n ; j++)
-       {
-           ASSERT (p [j+1] - p [j] >= 0) ;
-           for (pp = p [j] ; pp < p [j+1] ; pp++)
-           {
-               i = A [pp] ;
-               ASSERT (i >= 0 && i < n) ;
-               if (i > j && mark [i] != j)
-               {
-                   /* row k of M contains column indices i and j */
-                   M [count [i]++] = k ;
-                   M [count [j]++] = k ;
-                   k++ ;
-                   mark [i] = j ;
-               }
-           }
-       }
-       /* v2.4: free(mark) moved below */
-    }
-
-    /* count and mark no longer needed */
-    (*release) ((void *) count) ;
-    (*release) ((void *) mark) ;       /* v2.4: free (mark) moved here */
-    ASSERT (k == n_row) ;
-
-    /* === Adjust the knobs for M =========================================== */
-
-    for (i = 0 ; i < COLAMD_KNOBS ; i++)
-    {
-       cknobs [i] = knobs [i] ;
-    }
-
-    /* there are no dense rows in M */
-    cknobs [COLAMD_DENSE_ROW] = -1 ;
-    cknobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW] ;
-
-    /* === Order the columns of M =========================================== */
-
-    /* v2.4: colamd cannot fail here, so the error check is removed */
-    (void) COLAMD_MAIN (n_row, n, (Int) Mlen, M, perm, cknobs, stats) ;
-
-    /* Note that the output permutation is now in perm */
-
-    /* === get the statistics for symamd from colamd ======================== */
-
-    /* a dense column in colamd means a dense row and col in symamd */
-    stats [COLAMD_DENSE_ROW] = stats [COLAMD_DENSE_COL] ;
-
-    /* === Free M =========================================================== */
-
-    (*release) ((void *) M) ;
-    DEBUG0 (("symamd: done.\n")) ;
-    return (TRUE) ;
-
-}
-
-/* ========================================================================== */
-/* === colamd =============================================================== */
-/* ========================================================================== */
-
-/*
-    The colamd routine computes a column ordering Q of a sparse matrix
-    A such that the LU factorization P(AQ) = LU remains sparse, where P is
-    selected via partial pivoting.   The routine can also be viewed as
-    providing a permutation Q such that the Cholesky factorization
-    (AQ)'(AQ) = LL' remains sparse.
-*/
-
-PUBLIC Int COLAMD_MAIN         /* returns TRUE if successful, FALSE otherwise*/
-(
-    /* === Parameters ======================================================= */
-
-    Int n_row,                 /* number of rows in A */
-    Int n_col,                 /* number of columns in A */
-    Int Alen,                  /* length of A */
-    Int A [],                  /* row indices of A */
-    Int p [],                  /* pointers to columns in A */
-    double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */
-    Int stats [COLAMD_STATS]   /* output statistics and error codes */
-)
-{
-    /* === Local variables ================================================== */
-
-    Int i ;                    /* loop index */
-    Int nnz ;                  /* nonzeros in A */
-    size_t Row_size ;          /* size of Row [], in integers */
-    size_t Col_size ;          /* size of Col [], in integers */
-    size_t need ;              /* minimum required length of A */
-    Colamd_Row *Row ;          /* pointer into A of Row [0..n_row] array */
-    Colamd_Col *Col ;          /* pointer into A of Col [0..n_col] array */
-    Int n_col2 ;               /* number of non-dense, non-empty columns */
-    Int n_row2 ;               /* number of non-dense, non-empty rows */
-    Int ngarbage ;             /* number of garbage collections performed */
-    Int max_deg ;              /* maximum row degree */
-    double default_knobs [COLAMD_KNOBS] ;      /* default knobs array */
-    Int aggressive ;           /* do aggressive absorption */
-    int ok ;
-
-#ifndef NDEBUG
-    colamd_get_debug ("colamd") ;
-#endif /* NDEBUG */
-
-    /* === Check the input arguments ======================================== */
-
-    if (!stats)
-    {
-       DEBUG0 (("colamd: stats not present\n")) ;
-       return (FALSE) ;
-    }
-    for (i = 0 ; i < COLAMD_STATS ; i++)
-    {
-       stats [i] = 0 ;
-    }
-    stats [COLAMD_STATUS] = COLAMD_OK ;
-    stats [COLAMD_INFO1] = -1 ;
-    stats [COLAMD_INFO2] = -1 ;
-
-    if (!A)            /* A is not present */
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
-       DEBUG0 (("colamd: A not present\n")) ;
-       return (FALSE) ;
-    }
-
-    if (!p)            /* p is not present */
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
-       DEBUG0 (("colamd: p not present\n")) ;
-       return (FALSE) ;
-    }
-
-    if (n_row < 0)     /* n_row must be >= 0 */
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
-       stats [COLAMD_INFO1] = n_row ;
-       DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
-       return (FALSE) ;
-    }
-
-    if (n_col < 0)     /* n_col must be >= 0 */
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
-       stats [COLAMD_INFO1] = n_col ;
-       DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
-       return (FALSE) ;
-    }
-
-    nnz = p [n_col] ;
-    if (nnz < 0)       /* nnz must be >= 0 */
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
-       stats [COLAMD_INFO1] = nnz ;
-       DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
-       return (FALSE) ;
-    }
-
-    if (p [0] != 0)
-    {
-       stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
-       stats [COLAMD_INFO1] = p [0] ;
-       DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
-       return (FALSE) ;
-    }
-
-    /* === If no knobs, set default knobs =================================== */
-
-    if (!knobs)
-    {
-       COLAMD_set_defaults (default_knobs) ;
-       knobs = default_knobs ;
-    }
-
-    aggressive = (knobs [COLAMD_AGGRESSIVE] != FALSE) ;
-
-    /* === Allocate the Row and Col arrays from array A ===================== */
-
-    ok = TRUE ;
-    Col_size = COLAMD_C (n_col, &ok) ;     /* size of Col array of structs */
-    Row_size = COLAMD_R (n_row, &ok) ;     /* size of Row array of structs */
-
-    /* need = 2*nnz + n_col + Col_size + Row_size ; */
-    need = t_mult (nnz, 2, &ok) ;
-    need = t_add (need, n_col, &ok) ;
-    need = t_add (need, Col_size, &ok) ;
-    need = t_add (need, Row_size, &ok) ;
-
-    if (!ok || need > (size_t) Alen || need > Int_MAX)
-    {
-       /* not enough space in array A to perform the ordering */
-       stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
-       stats [COLAMD_INFO1] = need ;
-       stats [COLAMD_INFO2] = Alen ;
-       DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
-       return (FALSE) ;
-    }
-
-    Alen -= Col_size + Row_size ;
-    Col = (Colamd_Col *) &A [Alen] ;
-    Row = (Colamd_Row *) &A [Alen + Col_size] ;
-
-    /* === Construct the row and column data structures ===================== */
-
-    if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
-    {
-       /* input matrix is invalid */
-       DEBUG0 (("colamd: Matrix invalid\n")) ;
-       return (FALSE) ;
-    }
-
-    /* === Initialize scores, kill dense rows/columns ======================= */
-
-    init_scoring (n_row, n_col, Row, Col, A, p, knobs,
-       &n_row2, &n_col2, &max_deg) ;
-
-    /* === Order the supercolumns =========================================== */
-
-    ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
-       n_col2, max_deg, 2*nnz, aggressive) ;
-
-    /* === Order the non-principal columns ================================== */
-
-    order_children (n_col, Col, p) ;
-
-    /* === Return statistics in stats ======================================= */
-
-    stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
-    stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
-    stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
-    DEBUG0 (("colamd: done.\n")) ; 
-    return (TRUE) ;
-}
-
-
-/* ========================================================================== */
-/* === colamd_report ======================================================== */
-/* ========================================================================== */
-
-PUBLIC void COLAMD_report
-(
-    Int stats [COLAMD_STATS]
-)
-{
-    print_report ("colamd", stats) ;
-}
-
-
-/* ========================================================================== */
-/* === symamd_report ======================================================== */
-/* ========================================================================== */
-
-PUBLIC void SYMAMD_report
-(
-    Int stats [COLAMD_STATS]
-)
-{
-    print_report ("symamd", stats) ;
-}
-
-
-
-/* ========================================================================== */
-/* === NON-USER-CALLABLE ROUTINES: ========================================== */
-/* ========================================================================== */
-
-/* There are no user-callable routines beyond this point in the file */
-
-
-/* ========================================================================== */
-/* === init_rows_cols ======================================================= */
-/* ========================================================================== */
-
-/*
-    Takes the column form of the matrix in A and creates the row form of the
-    matrix.  Also, row and column attributes are stored in the Col and Row
-    structs.  If the columns are un-sorted or contain duplicate row indices,
-    this routine will also sort and remove duplicate row indices from the
-    column form of the matrix.  Returns FALSE if the matrix is invalid,
-    TRUE otherwise.  Not user-callable.
-*/
-
-PRIVATE Int init_rows_cols     /* returns TRUE if OK, or FALSE otherwise */
-(
-    /* === Parameters ======================================================= */
-
-    Int n_row,                 /* number of rows of A */
-    Int n_col,                 /* number of columns of A */
-    Colamd_Row Row [],         /* of size n_row+1 */
-    Colamd_Col Col [],         /* of size n_col+1 */
-    Int A [],                  /* row indices of A, of size Alen */
-    Int p [],                  /* pointers to columns in A, of size n_col+1 */
-    Int stats [COLAMD_STATS]   /* colamd statistics */ 
-)
-{
-    /* === Local variables ================================================== */
-
-    Int col ;                  /* a column index */
-    Int row ;                  /* a row index */
-    Int *cp ;                  /* a column pointer */
-    Int *cp_end ;              /* a pointer to the end of a column */
-    Int *rp ;                  /* a row pointer */
-    Int *rp_end ;              /* a pointer to the end of a row */
-    Int last_row ;             /* previous row */
-
-    /* === Initialize columns, and check column pointers ==================== */
-
-    for (col = 0 ; col < n_col ; col++)
-    {
-       Col [col].start = p [col] ;
-       Col [col].length = p [col+1] - p [col] ;
-
-       if (Col [col].length < 0)
-       {
-           /* column pointers must be non-decreasing */
-           stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
-           stats [COLAMD_INFO1] = col ;
-           stats [COLAMD_INFO2] = Col [col].length ;
-           DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
-           return (FALSE) ;
-       }
-
-       Col [col].shared1.thickness = 1 ;
-       Col [col].shared2.score = 0 ;
-       Col [col].shared3.prev = EMPTY ;
-       Col [col].shared4.degree_next = EMPTY ;
-    }
-
-    /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
-
-    /* === Scan columns, compute row degrees, and check row indices ========= */
-
-    stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
-
-    for (row = 0 ; row < n_row ; row++)
-    {
-       Row [row].length = 0 ;
-       Row [row].shared2.mark = -1 ;
-    }
-
-    for (col = 0 ; col < n_col ; col++)
-    {
-       last_row = -1 ;
-
-       cp = &A [p [col]] ;
-       cp_end = &A [p [col+1]] ;
-
-       while (cp < cp_end)
-       {
-           row = *cp++ ;
-
-           /* make sure row indices within range */
-           if (row < 0 || row >= n_row)
-           {
-               stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
-               stats [COLAMD_INFO1] = col ;
-               stats [COLAMD_INFO2] = row ;
-               stats [COLAMD_INFO3] = n_row ;
-               DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
-               return (FALSE) ;
-           }
-
-           if (row <= last_row || Row [row].shared2.mark == col)
-           {
-               /* row index are unsorted or repeated (or both), thus col */
-               /* is jumbled.  This is a notice, not an error condition. */
-               stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
-               stats [COLAMD_INFO1] = col ;
-               stats [COLAMD_INFO2] = row ;
-               (stats [COLAMD_INFO3]) ++ ;
-               DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
-           }
-
-           if (Row [row].shared2.mark != col)
-           {
-               Row [row].length++ ;
-           }
-           else
-           {
-               /* this is a repeated entry in the column, */
-               /* it will be removed */
-               Col [col].length-- ;
-           }
-
-           /* mark the row as having been seen in this column */
-           Row [row].shared2.mark = col ;
-
-           last_row = row ;
-       }
-    }
-
-    /* === Compute row pointers ============================================= */
-
-    /* row form of the matrix starts directly after the column */
-    /* form of matrix in A */
-    Row [0].start = p [n_col] ;
-    Row [0].shared1.p = Row [0].start ;
-    Row [0].shared2.mark = -1 ;
-    for (row = 1 ; row < n_row ; row++)
-    {
-       Row [row].start = Row [row-1].start + Row [row-1].length ;
-       Row [row].shared1.p = Row [row].start ;
-       Row [row].shared2.mark = -1 ;
-    }
-
-    /* === Create row form ================================================== */
-
-    if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
-    {
-       /* if cols jumbled, watch for repeated row indices */
-       for (col = 0 ; col < n_col ; col++)
-       {
-           cp = &A [p [col]] ;
-           cp_end = &A [p [col+1]] ;
-           while (cp < cp_end)
-           {
-               row = *cp++ ;
-               if (Row [row].shared2.mark != col)
-               {
-                   A [(Row [row].shared1.p)++] = col ;
-                   Row [row].shared2.mark = col ;
-               }
-           }
-       }
-    }
-    else
-    {
-       /* if cols not jumbled, we don't need the mark (this is faster) */
-       for (col = 0 ; col < n_col ; col++)
-       {
-           cp = &A [p [col]] ;
-           cp_end = &A [p [col+1]] ;
-           while (cp < cp_end)
-           {
-               A [(Row [*cp++].shared1.p)++] = col ;
-           }
-       }
-    }
-
-    /* === Clear the row marks and set row degrees ========================== */
-
-    for (row = 0 ; row < n_row ; row++)
-    {
-       Row [row].shared2.mark = 0 ;
-       Row [row].shared1.degree = Row [row].length ;
-    }
-
-    /* === See if we need to re-create columns ============================== */
-
-    if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
-    {
-       DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
-
-#ifndef NDEBUG
-       /* make sure column lengths are correct */
-       for (col = 0 ; col < n_col ; col++)
-       {
-           p [col] = Col [col].length ;
-       }
-       for (row = 0 ; row < n_row ; row++)
-       {
-           rp = &A [Row [row].start] ;
-           rp_end = rp + Row [row].length ;
-           while (rp < rp_end)
-           {
-               p [*rp++]-- ;
-           }
-       }
-       for (col = 0 ; col < n_col ; col++)
-       {
-           ASSERT (p [col] == 0) ;
-       }
-       /* now p is all zero (different than when debugging is turned off) */
-#endif /* NDEBUG */
-
-       /* === Compute col pointers ========================================= */
-
-       /* col form of the matrix starts at A [0]. */
-       /* Note, we may have a gap between the col form and the row */
-       /* form if there were duplicate entries, if so, it will be */
-       /* removed upon the first garbage collection */
-       Col [0].start = 0 ;
-       p [0] = Col [0].start ;
-       for (col = 1 ; col < n_col ; col++)
-       {
-           /* note that the lengths here are for pruned columns, i.e. */
-           /* no duplicate row indices will exist for these columns */
-           Col [col].start = Col [col-1].start + Col [col-1].length ;
-           p [col] = Col [col].start ;
-       }
-
-       /* === Re-create col form =========================================== */
-
-       for (row = 0 ; row < n_row ; row++)
-       {
-           rp = &A [Row [row].start] ;
-           rp_end = rp + Row [row].length ;
-           while (rp < rp_end)
-           {
-               A [(p [*rp++])++] = row ;
-           }
-       }
-    }
-
-    /* === Done.  Matrix is not (or no longer) jumbled ====================== */
-
-    return (TRUE) ;
-}
-
-
-/* ========================================================================== */
-/* === init_scoring ========================================================= */
-/* ========================================================================== */
-
-/*
-    Kills dense or empty columns and rows, calculates an initial score for
-    each column, and places all columns in the degree lists.  Not user-callable.
-*/
-
-PRIVATE void init_scoring
-(
-    /* === Parameters ======================================================= */
-
-    Int n_row,                 /* number of rows of A */
-    Int n_col,                 /* number of columns of A */
-    Colamd_Row Row [],         /* of size n_row+1 */
-    Colamd_Col Col [],         /* of size n_col+1 */
-    Int A [],                  /* column form and row form of A */
-    Int head [],               /* of size n_col+1 */
-    double knobs [COLAMD_KNOBS],/* parameters */
-    Int *p_n_row2,             /* number of non-dense, non-empty rows */
-    Int *p_n_col2,             /* number of non-dense, non-empty columns */
-    Int *p_max_deg             /* maximum row degree */
-)
-{
-    /* === Local variables ================================================== */
-
-    Int c ;                    /* a column index */
-    Int r, row ;               /* a row index */
-    Int *cp ;                  /* a column pointer */
-    Int deg ;                  /* degree of a row or column */
-    Int *cp_end ;              /* a pointer to the end of a column */
-    Int *new_cp ;              /* new column pointer */
-    Int col_length ;           /* length of pruned column */
-    Int score ;                        /* current column score */
-    Int n_col2 ;               /* number of non-dense, non-empty columns */
-    Int n_row2 ;               /* number of non-dense, non-empty rows */
-    Int dense_row_count ;      /* remove rows with more entries than this */
-    Int dense_col_count ;      /* remove cols with more entries than this */
-    Int min_score ;            /* smallest column score */
-    Int max_deg ;              /* maximum row degree */
-    Int next_col ;             /* Used to add to degree list.*/
-
-#ifndef NDEBUG
-    Int debug_count ;          /* debug only. */
-#endif /* NDEBUG */
-
-    /* === Extract knobs ==================================================== */
-
-    /* Note: if knobs contains a NaN, this is undefined: */
-    if (knobs [COLAMD_DENSE_ROW] < 0)
-    {
-       /* only remove completely dense rows */
-       dense_row_count = n_col-1 ;
-    }
-    else
-    {
-       dense_row_count = DENSE_DEGREE (knobs [COLAMD_DENSE_ROW], n_col) ;
-    }
-    if (knobs [COLAMD_DENSE_COL] < 0)
-    {
-       /* only remove completely dense columns */
-       dense_col_count = n_row-1 ;
-    }
-    else
-    {
-       dense_col_count =
-           DENSE_DEGREE (knobs [COLAMD_DENSE_COL], MIN (n_row, n_col)) ;
-    }
-
-    DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
-    max_deg = 0 ;
-    n_col2 = n_col ;
-    n_row2 = n_row ;
-
-    /* === Kill empty columns =============================================== */
-
-    /* Put the empty columns at the end in their natural order, so that LU */
-    /* factorization can proceed as far as possible. */
-    for (c = n_col-1 ; c >= 0 ; c--)
-    {
-       deg = Col [c].length ;
-       if (deg == 0)
-       {
-           /* this is a empty column, kill and order it last */
-           Col [c].shared2.order = --n_col2 ;
-           KILL_PRINCIPAL_COL (c) ;
-       }
-    }
-    DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
-
-    /* === Kill dense columns =============================================== */
-
-    /* Put the dense columns at the end, in their natural order */
-    for (c = n_col-1 ; c >= 0 ; c--)
-    {
-       /* skip any dead columns */
-       if (COL_IS_DEAD (c))
-       {
-           continue ;
-       }
-       deg = Col [c].length ;
-       if (deg > dense_col_count)
-       {
-           /* this is a dense column, kill and order it last */
-           Col [c].shared2.order = --n_col2 ;
-           /* decrement the row degrees */
-           cp = &A [Col [c].start] ;
-           cp_end = cp + Col [c].length ;
-           while (cp < cp_end)
-           {
-               Row [*cp++].shared1.degree-- ;
-           }
-           KILL_PRINCIPAL_COL (c) ;
-       }
-    }
-    DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
-
-    /* === Kill dense and empty rows ======================================== */
-
-    for (r = 0 ; r < n_row ; r++)
-    {
-       deg = Row [r].shared1.degree ;
-       ASSERT (deg >= 0 && deg <= n_col) ;
-       if (deg > dense_row_count || deg == 0)
-       {
-           /* kill a dense or empty row */
-           KILL_ROW (r) ;
-           --n_row2 ;
-       }
-       else
-       {
-           /* keep track of max degree of remaining rows */
-           max_deg = MAX (max_deg, deg) ;
-       }
-    }
-    DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
-
-    /* === Compute initial column scores ==================================== */
-
-    /* At this point the row degrees are accurate.  They reflect the number */
-    /* of "live" (non-dense) columns in each row.  No empty rows exist. */
-    /* Some "live" columns may contain only dead rows, however.  These are */
-    /* pruned in the code below. */
-
-    /* now find the initial matlab score for each column */
-    for (c = n_col-1 ; c >= 0 ; c--)
-    {
-       /* skip dead column */
-       if (COL_IS_DEAD (c))
-       {
-           continue ;
-       }
-       score = 0 ;
-       cp = &A [Col [c].start] ;
-       new_cp = cp ;
-       cp_end = cp + Col [c].length ;
-       while (cp < cp_end)
-       {
-           /* get a row */
-           row = *cp++ ;
-           /* skip if dead */
-           if (ROW_IS_DEAD (row))
-           {
-               continue ;
-           }
-           /* compact the column */
-           *new_cp++ = row ;
-           /* add row's external degree */
-           score += Row [row].shared1.degree - 1 ;
-           /* guard against integer overflow */
-           score = MIN (score, n_col) ;
-       }
-       /* determine pruned column length */
-       col_length = (Int) (new_cp - &A [Col [c].start]) ;
-       if (col_length == 0)
-       {
-           /* a newly-made null column (all rows in this col are "dense" */
-           /* and have already been killed) */
-           DEBUG2 (("Newly null killed: %d\n", c)) ;
-           Col [c].shared2.order = --n_col2 ;
-           KILL_PRINCIPAL_COL (c) ;
-       }
-       else
-       {
-           /* set column length and set score */
-           ASSERT (score >= 0) ;
-           ASSERT (score <= n_col) ;
-           Col [c].length = col_length ;
-           Col [c].shared2.score = score ;
-       }
-    }
-    DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
-       n_col-n_col2)) ;
-
-    /* At this point, all empty rows and columns are dead.  All live columns */
-    /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
-    /* yet).  Rows may contain dead columns, but all live rows contain at */
-    /* least one live column. */
-
-#ifndef NDEBUG
-    debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
-#endif /* NDEBUG */
-
-    /* === Initialize degree lists ========================================== */
-
-#ifndef NDEBUG
-    debug_count = 0 ;
-#endif /* NDEBUG */
-
-    /* clear the hash buckets */
-    for (c = 0 ; c <= n_col ; c++)
-    {
-       head [c] = EMPTY ;
-    }
-    min_score = n_col ;
-    /* place in reverse order, so low column indices are at the front */
-    /* of the lists.  This is to encourage natural tie-breaking */
-    for (c = n_col-1 ; c >= 0 ; c--)
-    {
-       /* only add principal columns to degree lists */
-       if (COL_IS_ALIVE (c))
-       {
-           DEBUG4 (("place %d score %d minscore %d ncol %d\n",
-               c, Col [c].shared2.score, min_score, n_col)) ;
-
-           /* === Add columns score to DList =============================== */
-
-           score = Col [c].shared2.score ;
-
-           ASSERT (min_score >= 0) ;
-           ASSERT (min_score <= n_col) ;
-           ASSERT (score >= 0) ;
-           ASSERT (score <= n_col) ;
-           ASSERT (head [score] >= EMPTY) ;
-
-           /* now add this column to dList at proper score location */
-           next_col = head [score] ;
-           Col [c].shared3.prev = EMPTY ;
-           Col [c].shared4.degree_next = next_col ;
-
-           /* if there already was a column with the same score, set its */
-           /* previous pointer to this new column */
-           if (next_col != EMPTY)
-           {
-               Col [next_col].shared3.prev = c ;
-           }
-           head [score] = c ;
-
-           /* see if this score is less than current min */
-           min_score = MIN (min_score, score) ;
-
-#ifndef NDEBUG
-           debug_count++ ;
-#endif /* NDEBUG */
-
-       }
-    }
-
-#ifndef NDEBUG
-    DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n",
-       debug_count, n_col, n_col-debug_count)) ;
-    ASSERT (debug_count == n_col2) ;
-    debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
-#endif /* NDEBUG */
-
-    /* === Return number of remaining columns, and max row degree =========== */
-
-    *p_n_col2 = n_col2 ;
-    *p_n_row2 = n_row2 ;
-    *p_max_deg = max_deg ;
-}
-
-
-/* ========================================================================== */
-/* === find_ordering ======================================================== */
-/* ========================================================================== */
-
-/*
-    Order the principal columns of the supercolumn form of the matrix
-    (no supercolumns on input).  Uses a minimum approximate column minimum
-    degree ordering method.  Not user-callable.
-*/
-
-PRIVATE Int find_ordering      /* return the number of garbage collections */
-(
-    /* === Parameters ======================================================= */
-
-    Int n_row,                 /* number of rows of A */
-    Int n_col,                 /* number of columns of A */
-    Int Alen,                  /* size of A, 2*nnz + n_col or larger */
-    Colamd_Row Row [],         /* of size n_row+1 */
-    Colamd_Col Col [],         /* of size n_col+1 */
-    Int A [],                  /* column form and row form of A */
-    Int head [],               /* of size n_col+1 */
-    Int n_col2,                        /* Remaining columns to order */
-    Int max_deg,               /* Maximum row degree */
-    Int pfree,                 /* index of first free slot (2*nnz on entry) */
-    Int aggressive
-)
-{
-    /* === Local variables ================================================== */
-
-    Int k ;                    /* current pivot ordering step */
-    Int pivot_col ;            /* current pivot column */
-    Int *cp ;                  /* a column pointer */
-    Int *rp ;                  /* a row pointer */
-    Int pivot_row ;            /* current pivot row */
-    Int *new_cp ;              /* modified column pointer */
-    Int *new_rp ;              /* modified row pointer */
-    Int pivot_row_start ;      /* pointer to start of pivot row */
-    Int pivot_row_degree ;     /* number of columns in pivot row */
-    Int pivot_row_length ;     /* number of supercolumns in pivot row */
-    Int pivot_col_score ;      /* score of pivot column */
-    Int needed_memory ;                /* free space needed for pivot row */
-    Int *cp_end ;              /* pointer to the end of a column */
-    Int *rp_end ;              /* pointer to the end of a row */
-    Int row ;                  /* a row index */
-    Int col ;                  /* a column index */
-    Int max_score ;            /* maximum possible score */
-    Int cur_score ;            /* score of current column */
-    unsigned Int hash ;                /* hash value for supernode detection */
-    Int head_column ;          /* head of hash bucket */
-    Int first_col ;            /* first column in hash bucket */
-    Int tag_mark ;             /* marker value for mark array */
-    Int row_mark ;             /* Row [row].shared2.mark */
-    Int set_difference ;       /* set difference size of row with pivot row */
-    Int min_score ;            /* smallest column score */
-    Int col_thickness ;                /* "thickness" (no. of columns in a supercol) */
-    Int max_mark ;             /* maximum value of tag_mark */
-    Int pivot_col_thickness ;  /* number of columns represented by pivot col */
-    Int prev_col ;             /* Used by Dlist operations. */
-    Int next_col ;             /* Used by Dlist operations. */
-    Int ngarbage ;             /* number of garbage collections performed */
-
-#ifndef NDEBUG
-    Int debug_d ;              /* debug loop counter */
-    Int debug_step = 0 ;       /* debug loop counter */
-#endif /* NDEBUG */
-
-    /* === Initialization and clear mark ==================================== */
-
-    max_mark = INT_MAX - n_col ;       /* INT_MAX defined in <limits.h> */
-    tag_mark = clear_mark (0, max_mark, n_row, Row) ;
-    min_score = 0 ;
-    ngarbage = 0 ;
-    DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
-
-    /* === Order the columns ================================================ */
-
-    for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
-    {
-
-#ifndef NDEBUG
-       if (debug_step % 100 == 0)
-       {
-           DEBUG2 (("\n...       Step k: %d out of n_col2: %d\n", k, n_col2)) ;
-       }
-       else
-       {
-           DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
-       }
-       debug_step++ ;
-       debug_deg_lists (n_row, n_col, Row, Col, head,
-               min_score, n_col2-k, max_deg) ;
-       debug_matrix (n_row, n_col, Row, Col, A) ;
-#endif /* NDEBUG */
-
-       /* === Select pivot column, and order it ============================ */
-
-       /* make sure degree list isn't empty */
-       ASSERT (min_score >= 0) ;
-       ASSERT (min_score <= n_col) ;
-       ASSERT (head [min_score] >= EMPTY) ;
-
-#ifndef NDEBUG
-       for (debug_d = 0 ; debug_d < min_score ; debug_d++)
-       {
-           ASSERT (head [debug_d] == EMPTY) ;
-       }
-#endif /* NDEBUG */
-
-       /* get pivot column from head of minimum degree list */
-       while (head [min_score] == EMPTY && min_score < n_col)
-       {
-           min_score++ ;
-       }
-       pivot_col = head [min_score] ;
-       ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
-       next_col = Col [pivot_col].shared4.degree_next ;
-       head [min_score] = next_col ;
-       if (next_col != EMPTY)
-       {
-           Col [next_col].shared3.prev = EMPTY ;
-       }
-
-       ASSERT (COL_IS_ALIVE (pivot_col)) ;
-
-       /* remember score for defrag check */
-       pivot_col_score = Col [pivot_col].shared2.score ;
-
-       /* the pivot column is the kth column in the pivot order */
-       Col [pivot_col].shared2.order = k ;
-
-       /* increment order count by column thickness */
-       pivot_col_thickness = Col [pivot_col].shared1.thickness ;
-       k += pivot_col_thickness ;
-       ASSERT (pivot_col_thickness > 0) ;
-       DEBUG3 (("Pivot col: %d thick %d\n", pivot_col, pivot_col_thickness)) ;
-
-       /* === Garbage_collection, if necessary ============================= */
-
-       needed_memory = MIN (pivot_col_score, n_col - k) ;
-       if (pfree + needed_memory >= Alen)
-       {
-           pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
-           ngarbage++ ;
-           /* after garbage collection we will have enough */
-           ASSERT (pfree + needed_memory < Alen) ;
-           /* garbage collection has wiped out the Row[].shared2.mark array */
-           tag_mark = clear_mark (0, max_mark, n_row, Row) ;
-
-#ifndef NDEBUG
-           debug_matrix (n_row, n_col, Row, Col, A) ;
-#endif /* NDEBUG */
-       }
-
-       /* === Compute pivot row pattern ==================================== */
-
-       /* get starting location for this new merged row */
-       pivot_row_start = pfree ;
-
-       /* initialize new row counts to zero */
-       pivot_row_degree = 0 ;
-
-       /* tag pivot column as having been visited so it isn't included */
-       /* in merged pivot row */
-       Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
-
-       /* pivot row is the union of all rows in the pivot column pattern */
-       cp = &A [Col [pivot_col].start] ;
-       cp_end = cp + Col [pivot_col].length ;
-       while (cp < cp_end)
-       {
-           /* get a row */
-           row = *cp++ ;
-           DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
-           /* skip if row is dead */
-           if (ROW_IS_ALIVE (row))
-           {
-               rp = &A [Row [row].start] ;
-               rp_end = rp + Row [row].length ;
-               while (rp < rp_end)
-               {
-                   /* get a column */
-                   col = *rp++ ;
-                   /* add the column, if alive and untagged */
-                   col_thickness = Col [col].shared1.thickness ;
-                   if (col_thickness > 0 && COL_IS_ALIVE (col))
-                   {
-                       /* tag column in pivot row */
-                       Col [col].shared1.thickness = -col_thickness ;
-                       ASSERT (pfree < Alen) ;
-                       /* place column in pivot row */
-                       A [pfree++] = col ;
-                       pivot_row_degree += col_thickness ;
-                   }
-               }
-           }
-       }
-
-       /* clear tag on pivot column */
-       Col [pivot_col].shared1.thickness = pivot_col_thickness ;
-       max_deg = MAX (max_deg, pivot_row_degree) ;
-
-#ifndef NDEBUG
-       DEBUG3 (("check2\n")) ;
-       debug_mark (n_row, Row, tag_mark, max_mark) ;
-#endif /* NDEBUG */
-
-       /* === Kill all rows used to construct pivot row ==================== */
-
-       /* also kill pivot row, temporarily */
-       cp = &A [Col [pivot_col].start] ;
-       cp_end = cp + Col [pivot_col].length ;
-       while (cp < cp_end)
-       {
-           /* may be killing an already dead row */
-           row = *cp++ ;
-           DEBUG3 (("Kill row in pivot col: %d\n", row)) ;
-           KILL_ROW (row) ;
-       }
-
-       /* === Select a row index to use as the new pivot row =============== */
-
-       pivot_row_length = pfree - pivot_row_start ;
-       if (pivot_row_length > 0)
-       {
-           /* pick the "pivot" row arbitrarily (first row in col) */
-           pivot_row = A [Col [pivot_col].start] ;
-           DEBUG3 (("Pivotal row is %d\n", pivot_row)) ;
-       }
-       else
-       {
-           /* there is no pivot row, since it is of zero length */
-           pivot_row = EMPTY ;
-           ASSERT (pivot_row_length == 0) ;
-       }
-       ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
-
-       /* === Approximate degree computation =============================== */
-
-       /* Here begins the computation of the approximate degree.  The column */
-       /* score is the sum of the pivot row "length", plus the size of the */
-       /* set differences of each row in the column minus the pattern of the */
-       /* pivot row itself.  The column ("thickness") itself is also */
-       /* excluded from the column score (we thus use an approximate */
-       /* external degree). */
-
-       /* The time taken by the following code (compute set differences, and */
-       /* add them up) is proportional to the size of the data structure */
-       /* being scanned - that is, the sum of the sizes of each column in */
-       /* the pivot row.  Thus, the amortized time to compute a column score */
-       /* is proportional to the size of that column (where size, in this */
-       /* context, is the column "length", or the number of row indices */
-       /* in that column).  The number of row indices in a column is */
-       /* monotonically non-decreasing, from the length of the original */
-       /* column on input to colamd. */
-
-       /* === Compute set differences ====================================== */
-
-       DEBUG3 (("** Computing set differences phase. **\n")) ;
-
-       /* pivot row is currently dead - it will be revived later. */
-
-       DEBUG3 (("Pivot row: ")) ;
-       /* for each column in pivot row */
-       rp = &A [pivot_row_start] ;
-       rp_end = rp + pivot_row_length ;
-       while (rp < rp_end)
-       {
-           col = *rp++ ;
-           ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
-           DEBUG3 (("Col: %d\n", col)) ;
-
-           /* clear tags used to construct pivot row pattern */
-           col_thickness = -Col [col].shared1.thickness ;
-           ASSERT (col_thickness > 0) ;
-           Col [col].shared1.thickness = col_thickness ;
-
-           /* === Remove column from degree list =========================== */
-
-           cur_score = Col [col].shared2.score ;
-           prev_col = Col [col].shared3.prev ;
-           next_col = Col [col].shared4.degree_next ;
-           ASSERT (cur_score >= 0) ;
-           ASSERT (cur_score <= n_col) ;
-           ASSERT (cur_score >= EMPTY) ;
-           if (prev_col == EMPTY)
-           {
-               head [cur_score] = next_col ;
-           }
-           else
-           {
-               Col [prev_col].shared4.degree_next = next_col ;
-           }
-           if (next_col != EMPTY)
-           {
-               Col [next_col].shared3.prev = prev_col ;
-           }
-
-           /* === Scan the column ========================================== */
-
-           cp = &A [Col [col].start] ;
-           cp_end = cp + Col [col].length ;
-           while (cp < cp_end)
-           {
-               /* get a row */
-               row = *cp++ ;
-               row_mark = Row [row].shared2.mark ;
-               /* skip if dead */
-               if (ROW_IS_MARKED_DEAD (row_mark))
-               {
-                   continue ;
-               }
-               ASSERT (row != pivot_row) ;
-               set_difference = row_mark - tag_mark ;
-               /* check if the row has been seen yet */
-               if (set_difference < 0)
-               {
-                   ASSERT (Row [row].shared1.degree <= max_deg) ;
-                   set_difference = Row [row].shared1.degree ;
-               }
-               /* subtract column thickness from this row's set difference */
-               set_difference -= col_thickness ;
-               ASSERT (set_difference >= 0) ;
-               /* absorb this row if the set difference becomes zero */
-               if (set_difference == 0 && aggressive)
-               {
-                   DEBUG3 (("aggressive absorption. Row: %d\n", row)) ;
-                   KILL_ROW (row) ;
-               }
-               else
-               {
-                   /* save the new mark */
-                   Row [row].shared2.mark = set_difference + tag_mark ;
-               }
-           }
-       }
-
-#ifndef NDEBUG
-       debug_deg_lists (n_row, n_col, Row, Col, head,
-               min_score, n_col2-k-pivot_row_degree, max_deg) ;
-#endif /* NDEBUG */
-
-       /* === Add up set differences for each column ======================= */
-
-       DEBUG3 (("** Adding set differences phase. **\n")) ;
-
-       /* for each column in pivot row */
-       rp = &A [pivot_row_start] ;
-       rp_end = rp + pivot_row_length ;
-       while (rp < rp_end)
-       {
-           /* get a column */
-           col = *rp++ ;
-           ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
-           hash = 0 ;
-           cur_score = 0 ;
-           cp = &A [Col [col].start] ;
-           /* compact the column */
-           new_cp = cp ;
-           cp_end = cp + Col [col].length ;
-
-           DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ;
-
-           while (cp < cp_end)
-           {
-               /* get a row */
-               row = *cp++ ;
-               ASSERT(row >= 0 && row < n_row) ;
-               row_mark = Row [row].shared2.mark ;
-               /* skip if dead */
-               if (ROW_IS_MARKED_DEAD (row_mark))
-               {
-                   DEBUG4 ((" Row %d, dead\n", row)) ;
-                   continue ;
-               }
-               DEBUG4 ((" Row %d, set diff %d\n", row, row_mark-tag_mark));
-               ASSERT (row_mark >= tag_mark) ;
-               /* compact the column */
-               *new_cp++ = row ;
-               /* compute hash function */
-               hash += row ;
-               /* add set difference */
-               cur_score += row_mark - tag_mark ;
-               /* integer overflow... */
-               cur_score = MIN (cur_score, n_col) ;
-           }
-
-           /* recompute the column's length */
-           Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
-
-           /* === Further mass elimination ================================= */
-
-           if (Col [col].length == 0)
-           {
-               DEBUG4 (("further mass elimination. Col: %d\n", col)) ;
-               /* nothing left but the pivot row in this column */
-               KILL_PRINCIPAL_COL (col) ;
-               pivot_row_degree -= Col [col].shared1.thickness ;
-               ASSERT (pivot_row_degree >= 0) ;
-               /* order it */
-               Col [col].shared2.order = k ;
-               /* increment order count by column thickness */
-               k += Col [col].shared1.thickness ;
-           }
-           else
-           {
-               /* === Prepare for supercolumn detection ==================== */
-
-               DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ;
-
-               /* save score so far */
-               Col [col].shared2.score = cur_score ;
-
-               /* add column to hash table, for supercolumn detection */
-               hash %= n_col + 1 ;
-
-               DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ;
-               ASSERT (((Int) hash) <= n_col) ;
-
-               head_column = head [hash] ;
-               if (head_column > EMPTY)
-               {
-                   /* degree list "hash" is non-empty, use prev (shared3) of */
-                   /* first column in degree list as head of hash bucket */
-                   first_col = Col [head_column].shared3.headhash ;
-                   Col [head_column].shared3.headhash = col ;
-               }
-               else
-               {
-                   /* degree list "hash" is empty, use head as hash bucket */
-                   first_col = - (head_column + 2) ;
-                   head [hash] = - (col + 2) ;
-               }
-               Col [col].shared4.hash_next = first_col ;
-
-               /* save hash function in Col [col].shared3.hash */
-               Col [col].shared3.hash = (Int) hash ;
-               ASSERT (COL_IS_ALIVE (col)) ;
-           }
-       }
-
-       /* The approximate external column degree is now computed.  */
-
-       /* === Supercolumn detection ======================================== */
-
-       DEBUG3 (("** Supercolumn detection phase. **\n")) ;
-
-       detect_super_cols (
-
-#ifndef NDEBUG
-               n_col, Row,
-#endif /* NDEBUG */
-
-               Col, A, head, pivot_row_start, pivot_row_length) ;
-
-       /* === Kill the pivotal column ====================================== */
-
-       KILL_PRINCIPAL_COL (pivot_col) ;
-
-       /* === Clear mark =================================================== */
-
-       tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ;
-
-#ifndef NDEBUG
-       DEBUG3 (("check3\n")) ;
-       debug_mark (n_row, Row, tag_mark, max_mark) ;
-#endif /* NDEBUG */
-
-       /* === Finalize the new pivot row, and column scores ================ */
-
-       DEBUG3 (("** Finalize scores phase. **\n")) ;
-
-       /* for each column in pivot row */
-       rp = &A [pivot_row_start] ;
-       /* compact the pivot row */
-       new_rp = rp ;
-       rp_end = rp + pivot_row_length ;
-       while (rp < rp_end)
-       {
-           col = *rp++ ;
-           /* skip dead columns */
-           if (COL_IS_DEAD (col))
-           {
-               continue ;
-           }
-           *new_rp++ = col ;
-           /* add new pivot row to column */
-           A [Col [col].start + (Col [col].length++)] = pivot_row ;
-
-           /* retrieve score so far and add on pivot row's degree. */
-           /* (we wait until here for this in case the pivot */
-           /* row's degree was reduced due to mass elimination). */
-           cur_score = Col [col].shared2.score + pivot_row_degree ;
-
-           /* calculate the max possible score as the number of */
-           /* external columns minus the 'k' value minus the */
-           /* columns thickness */
-           max_score = n_col - k - Col [col].shared1.thickness ;
-
-           /* make the score the external degree of the union-of-rows */
-           cur_score -= Col [col].shared1.thickness ;
-
-           /* make sure score is less or equal than the max score */
-           cur_score = MIN (cur_score, max_score) ;
-           ASSERT (cur_score >= 0) ;
-
-           /* store updated score */
-           Col [col].shared2.score = cur_score ;
-
-           /* === Place column back in degree list ========================= */
-
-           ASSERT (min_score >= 0) ;
-           ASSERT (min_score <= n_col) ;
-           ASSERT (cur_score >= 0) ;
-           ASSERT (cur_score <= n_col) ;
-           ASSERT (head [cur_score] >= EMPTY) ;
-           next_col = head [cur_score] ;
-           Col [col].shared4.degree_next = next_col ;
-           Col [col].shared3.prev = EMPTY ;
-           if (next_col != EMPTY)
-           {
-               Col [next_col].shared3.prev = col ;
-           }
-           head [cur_score] = col ;
-
-           /* see if this score is less than current min */
-           min_score = MIN (min_score, cur_score) ;
-
-       }
-
-#ifndef NDEBUG
-       debug_deg_lists (n_row, n_col, Row, Col, head,
-               min_score, n_col2-k, max_deg) ;
-#endif /* NDEBUG */
-
-       /* === Resurrect the new pivot row ================================== */
-
-       if (pivot_row_degree > 0)
-       {
-           /* update pivot row length to reflect any cols that were killed */
-           /* during super-col detection and mass elimination */
-           Row [pivot_row].start  = pivot_row_start ;
-           Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
-           ASSERT (Row [pivot_row].length > 0) ;
-           Row [pivot_row].shared1.degree = pivot_row_degree ;
-           Row [pivot_row].shared2.mark = 0 ;
-           /* pivot row is no longer dead */
-
-           DEBUG1 (("Resurrect Pivot_row %d deg: %d\n",
-                       pivot_row, pivot_row_degree)) ;
-       }
-    }
-
-    /* === All principal columns have now been ordered ====================== */
-
-    return (ngarbage) ;
-}
-
-
-/* ========================================================================== */
-/* === order_children ======================================================= */
-/* ========================================================================== */
-
-/*
-    The find_ordering routine has ordered all of the principal columns (the
-    representatives of the supercolumns).  The non-principal columns have not
-    yet been ordered.  This routine orders those columns by walking up the
-    parent tree (a column is a child of the column which absorbed it).  The
-    final permutation vector is then placed in p [0 ... n_col-1], with p [0]
-    being the first column, and p [n_col-1] being the last.  It doesn't look
-    like it at first glance, but be assured that this routine takes time linear
-    in the number of columns.  Although not immediately obvious, the time
-    taken by this routine is O (n_col), that is, linear in the number of
-    columns.  Not user-callable.
-*/
-
-PRIVATE void order_children
-(
-    /* === Parameters ======================================================= */
-
-    Int n_col,                 /* number of columns of A */
-    Colamd_Col Col [],         /* of size n_col+1 */
-    Int p []                   /* p [0 ... n_col-1] is the column permutation*/
-)
-{
-    /* === Local variables ================================================== */
-
-    Int i ;                    /* loop counter for all columns */
-    Int c ;                    /* column index */
-    Int parent ;               /* index of column's parent */
-    Int order ;                        /* column's order */
-
-    /* === Order each non-principal column ================================== */
-
-    for (i = 0 ; i < n_col ; i++)
-    {
-       /* find an un-ordered non-principal column */
-       ASSERT (COL_IS_DEAD (i)) ;
-       if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY)
-       {
-           parent = i ;
-           /* once found, find its principal parent */
-           do
-           {
-               parent = Col [parent].shared1.parent ;
-           } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
-
-           /* now, order all un-ordered non-principal columns along path */
-           /* to this parent.  collapse tree at the same time */
-           c = i ;
-           /* get order of parent */
-           order = Col [parent].shared2.order ;
-
-           do
-           {
-               ASSERT (Col [c].shared2.order == EMPTY) ;
-
-               /* order this column */
-               Col [c].shared2.order = order++ ;
-               /* collaps tree */
-               Col [c].shared1.parent = parent ;
-
-               /* get immediate parent of this column */
-               c = Col [c].shared1.parent ;
-
-               /* continue until we hit an ordered column.  There are */
-               /* guarranteed not to be anymore unordered columns */
-               /* above an ordered column */
-           } while (Col [c].shared2.order == EMPTY) ;
-
-           /* re-order the super_col parent to largest order for this group */
-           Col [parent].shared2.order = order ;
-       }
-    }
-
-    /* === Generate the permutation ========================================= */
-
-    for (c = 0 ; c < n_col ; c++)
-    {
-       p [Col [c].shared2.order] = c ;
-    }
-}
-
-
-/* ========================================================================== */
-/* === detect_super_cols ==================================================== */
-/* ========================================================================== */
-
-/*
-    Detects supercolumns by finding matches between columns in the hash buckets.
-    Check amongst columns in the set A [row_start ... row_start + row_length-1].
-    The columns under consideration are currently *not* in the degree lists,
-    and have already been placed in the hash buckets.
-
-    The hash bucket for columns whose hash function is equal to h is stored
-    as follows:
-
-       if head [h] is >= 0, then head [h] contains a degree list, so:
-
-               head [h] is the first column in degree bucket h.
-               Col [head [h]].headhash gives the first column in hash bucket h.
-
-       otherwise, the degree list is empty, and:
-
-               -(head [h] + 2) is the first column in hash bucket h.
-
-    For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
-    column" pointer.  Col [c].shared3.hash is used instead as the hash number
-    for that column.  The value of Col [c].shared4.hash_next is the next column
-    in the same hash bucket.
-
-    Assuming no, or "few" hash collisions, the time taken by this routine is
-    linear in the sum of the sizes (lengths) of each column whose score has
-    just been computed in the approximate degree computation.
-    Not user-callable.
-*/
-
-PRIVATE void detect_super_cols
-(
-    /* === Parameters ======================================================= */
-
-#ifndef NDEBUG
-    /* these two parameters are only needed when debugging is enabled: */
-    Int n_col,                 /* number of columns of A */
-    Colamd_Row Row [],         /* of size n_row+1 */
-#endif /* NDEBUG */
-
-    Colamd_Col Col [],         /* of size n_col+1 */
-    Int A [],                  /* row indices of A */
-    Int head [],               /* head of degree lists and hash buckets */
-    Int row_start,             /* pointer to set of columns to check */
-    Int row_length             /* number of columns to check */
-)
-{
-    /* === Local variables ================================================== */
-
-    Int hash ;                 /* hash value for a column */
-    Int *rp ;                  /* pointer to a row */
-    Int c ;                    /* a column index */
-    Int super_c ;              /* column index of the column to absorb into */
-    Int *cp1 ;                 /* column pointer for column super_c */
-    Int *cp2 ;                 /* column pointer for column c */
-    Int length ;               /* length of column super_c */
-    Int prev_c ;               /* column preceding c in hash bucket */
-    Int i ;                    /* loop counter */
-    Int *rp_end ;              /* pointer to the end of the row */
-    Int col ;                  /* a column index in the row to check */
-    Int head_column ;          /* first column in hash bucket or degree list */
-    Int first_col ;            /* first column in hash bucket */
-
-    /* === Consider each column in the row ================================== */
-
-    rp = &A [row_start] ;
-    rp_end = rp + row_length ;
-    while (rp < rp_end)
-    {
-       col = *rp++ ;
-       if (COL_IS_DEAD (col))
-       {
-           continue ;
-       }
-
-       /* get hash number for this column */
-       hash = Col [col].shared3.hash ;
-       ASSERT (hash <= n_col) ;
-
-       /* === Get the first column in this hash bucket ===================== */
-
-       head_column = head [hash] ;
-       if (head_column > EMPTY)
-       {
-           first_col = Col [head_column].shared3.headhash ;
-       }
-       else
-       {
-           first_col = - (head_column + 2) ;
-       }
-
-       /* === Consider each column in the hash bucket ====================== */
-
-       for (super_c = first_col ; super_c != EMPTY ;
-           super_c = Col [super_c].shared4.hash_next)
-       {
-           ASSERT (COL_IS_ALIVE (super_c)) ;
-           ASSERT (Col [super_c].shared3.hash == hash) ;
-           length = Col [super_c].length ;
-
-           /* prev_c is the column preceding column c in the hash bucket */
-           prev_c = super_c ;
-
-           /* === Compare super_c with all columns after it ================ */
-
-           for (c = Col [super_c].shared4.hash_next ;
-                c != EMPTY ; c = Col [c].shared4.hash_next)
-           {
-               ASSERT (c != super_c) ;
-               ASSERT (COL_IS_ALIVE (c)) ;
-               ASSERT (Col [c].shared3.hash == hash) ;
-
-               /* not identical if lengths or scores are different */
-               if (Col [c].length != length ||
-                   Col [c].shared2.score != Col [super_c].shared2.score)
-               {
-                   prev_c = c ;
-                   continue ;
-               }
-
-               /* compare the two columns */
-               cp1 = &A [Col [super_c].start] ;
-               cp2 = &A [Col [c].start] ;
-
-               for (i = 0 ; i < length ; i++)
-               {
-                   /* the columns are "clean" (no dead rows) */
-                   ASSERT (ROW_IS_ALIVE (*cp1))  ;
-                   ASSERT (ROW_IS_ALIVE (*cp2))  ;
-                   /* row indices will same order for both supercols, */
-                   /* no gather scatter nessasary */
-                   if (*cp1++ != *cp2++)
-                   {
-                       break ;
-                   }
-               }
-
-               /* the two columns are different if the for-loop "broke" */
-               if (i != length)
-               {
-                   prev_c = c ;
-                   continue ;
-               }
-
-               /* === Got it!  two columns are identical =================== */
-
-               ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
-
-               Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
-               Col [c].shared1.parent = super_c ;
-               KILL_NON_PRINCIPAL_COL (c) ;
-               /* order c later, in order_children() */
-               Col [c].shared2.order = EMPTY ;
-               /* remove c from hash bucket */
-               Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
-           }
-       }
-
-       /* === Empty this hash bucket ======================================= */
-
-       if (head_column > EMPTY)
-       {
-           /* corresponding degree list "hash" is not empty */
-           Col [head_column].shared3.headhash = EMPTY ;
-       }
-       else
-       {
-           /* corresponding degree list "hash" is empty */
-           head [hash] = EMPTY ;
-       }
-    }
-}
-
-
-/* ========================================================================== */
-/* === garbage_collection =================================================== */
-/* ========================================================================== */
-
-/*
-    Defragments and compacts columns and rows in the workspace A.  Used when
-    all avaliable memory has been used while performing row merging.  Returns
-    the index of the first free position in A, after garbage collection.  The
-    time taken by this routine is linear is the size of the array A, which is
-    itself linear in the number of nonzeros in the input matrix.
-    Not user-callable.
-*/
-
-PRIVATE Int garbage_collection  /* returns the new value of pfree */
-(
-    /* === Parameters ======================================================= */
-
-    Int n_row,                 /* number of rows */
-    Int n_col,                 /* number of columns */
-    Colamd_Row Row [],         /* row info */
-    Colamd_Col Col [],         /* column info */
-    Int A [],                  /* A [0 ... Alen-1] holds the matrix */
-    Int *pfree                 /* &A [0] ... pfree is in use */
-)
-{
-    /* === Local variables ================================================== */
-
-    Int *psrc ;                        /* source pointer */
-    Int *pdest ;               /* destination pointer */
-    Int j ;                    /* counter */
-    Int r ;                    /* a row index */
-    Int c ;                    /* a column index */
-    Int length ;               /* length of a row or column */
-
-#ifndef NDEBUG
-    Int debug_rows ;
-    DEBUG2 (("Defrag..\n")) ;
-    for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
-    debug_rows = 0 ;
-#endif /* NDEBUG */
-
-    /* === Defragment the columns =========================================== */
-
-    pdest = &A[0] ;
-    for (c = 0 ; c < n_col ; c++)
-    {
-       if (COL_IS_ALIVE (c))
-       {
-           psrc = &A [Col [c].start] ;
-
-           /* move and compact the column */
-           ASSERT (pdest <= psrc) ;
-           Col [c].start = (Int) (pdest - &A [0]) ;
-           length = Col [c].length ;
-           for (j = 0 ; j < length ; j++)
-           {
-               r = *psrc++ ;
-               if (ROW_IS_ALIVE (r))
-               {
-                   *pdest++ = r ;
-               }
-           }
-           Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
-       }
-    }
-
-    /* === Prepare to defragment the rows =================================== */
-
-    for (r = 0 ; r < n_row ; r++)
-    {
-       if (ROW_IS_DEAD (r) || (Row [r].length == 0))
-       {
-           /* This row is already dead, or is of zero length.  Cannot compact
-            * a row of zero length, so kill it.  NOTE: in the current version,
-            * there are no zero-length live rows.  Kill the row (for the first
-            * time, or again) just to be safe. */
-           KILL_ROW (r) ;
-       }
-       else
-       {
-           /* save first column index in Row [r].shared2.first_column */
-           psrc = &A [Row [r].start] ;
-           Row [r].shared2.first_column = *psrc ;
-           ASSERT (ROW_IS_ALIVE (r)) ;
-           /* flag the start of the row with the one's complement of row */
-           *psrc = ONES_COMPLEMENT (r) ;
-#ifndef NDEBUG
-           debug_rows++ ;
-#endif /* NDEBUG */
-       }
-    }
-
-    /* === Defragment the rows ============================================== */
-
-    psrc = pdest ;
-    while (psrc < pfree)
-    {
-       /* find a negative number ... the start of a row */
-       if (*psrc++ < 0)
-       {
-           psrc-- ;
-           /* get the row index */
-           r = ONES_COMPLEMENT (*psrc) ;
-           ASSERT (r >= 0 && r < n_row) ;
-           /* restore first column index */
-           *psrc = Row [r].shared2.first_column ;
-           ASSERT (ROW_IS_ALIVE (r)) ;
-           ASSERT (Row [r].length > 0) ;
-           /* move and compact the row */
-           ASSERT (pdest <= psrc) ;
-           Row [r].start = (Int) (pdest - &A [0]) ;
-           length = Row [r].length ;
-           for (j = 0 ; j < length ; j++)
-           {
-               c = *psrc++ ;
-               if (COL_IS_ALIVE (c))
-               {
-                   *pdest++ = c ;
-               }
-           }
-           Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
-           ASSERT (Row [r].length > 0) ;
-#ifndef NDEBUG
-           debug_rows-- ;
-#endif /* NDEBUG */
-       }
-    }
-    /* ensure we found all the rows */
-    ASSERT (debug_rows == 0) ;
-
-    /* === Return the new value of pfree ==================================== */
-
-    return ((Int) (pdest - &A [0])) ;
-}
-
-
-/* ========================================================================== */
-/* === clear_mark =========================================================== */
-/* ========================================================================== */
-
-/*
-    Clears the Row [].shared2.mark array, and returns the new tag_mark.
-    Return value is the new tag_mark.  Not user-callable.
-*/
-
-PRIVATE Int clear_mark /* return the new value for tag_mark */
-(
-    /* === Parameters ======================================================= */
-
-    Int tag_mark,      /* new value of tag_mark */
-    Int max_mark,      /* max allowed value of tag_mark */
-
-    Int n_row,         /* number of rows in A */
-    Colamd_Row Row []  /* Row [0 ... n_row-1].shared2.mark is set to zero */
-)
-{
-    /* === Local variables ================================================== */
-
-    Int r ;
-
-    if (tag_mark <= 0 || tag_mark >= max_mark)
-    {
-       for (r = 0 ; r < n_row ; r++)
-       {
-           if (ROW_IS_ALIVE (r))
-           {
-               Row [r].shared2.mark = 0 ;
-           }
-       }
-       tag_mark = 1 ;
-    }
-
-    return (tag_mark) ;
-}
-
-
-/* ========================================================================== */
-/* === print_report ========================================================= */
-/* ========================================================================== */
-
-PRIVATE void print_report
-(
-    char *method,
-    Int stats [COLAMD_STATS]
-)
-{
-
-    Int i1, i2, i3 ;
-
-    PRINTF (("\n%s version %d.%d, %s: ", method,
-           COLAMD_MAIN_VERSION, COLAMD_SUB_VERSION, COLAMD_DATE)) ;
-
-    if (!stats)
-    {
-       PRINTF (("No statistics available.\n")) ;
-       return ;
-    }
-
-    i1 = stats [COLAMD_INFO1] ;
-    i2 = stats [COLAMD_INFO2] ;
-    i3 = stats [COLAMD_INFO3] ;
-
-    if (stats [COLAMD_STATUS] >= 0)
-    {
-       PRINTF (("OK.  ")) ;
-    }
-    else
-    {
-       PRINTF (("ERROR.  ")) ;
-    }
-
-    switch (stats [COLAMD_STATUS])
-    {
-
-       case COLAMD_OK_BUT_JUMBLED:
-
-           PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ;
-
-           PRINTF(("%s: number of duplicate or out-of-order row indices: %d\n",
-           method, i3)) ;
-
-           PRINTF(("%s: last seen duplicate or out-of-order row index:   %d\n",
-           method, INDEX (i2))) ;
-
-           PRINTF(("%s: last seen in column:                             %d",
-           method, INDEX (i1))) ;
-
-           /* no break - fall through to next case instead */
-
-       case COLAMD_OK:
-
-           PRINTF(("\n")) ;
-
-           PRINTF(("%s: number of dense or empty rows ignored:           %d\n",
-           method, stats [COLAMD_DENSE_ROW])) ;
-
-           PRINTF(("%s: number of dense or empty columns ignored:        %d\n",
-           method, stats [COLAMD_DENSE_COL])) ;
-
-           PRINTF(("%s: number of garbage collections performed:         %d\n",
-           method, stats [COLAMD_DEFRAG_COUNT])) ;
-           break ;
-
-       case COLAMD_ERROR_A_not_present:
-
-           PRINTF(("Array A (row indices of matrix) not present.\n")) ;
-           break ;
-
-       case COLAMD_ERROR_p_not_present:
-
-           PRINTF(("Array p (column pointers for matrix) not present.\n")) ;
-           break ;
-
-       case COLAMD_ERROR_nrow_negative:
-
-           PRINTF(("Invalid number of rows (%d).\n", i1)) ;
-           break ;
-
-       case COLAMD_ERROR_ncol_negative:
-
-           PRINTF(("Invalid number of columns (%d).\n", i1)) ;
-           break ;
-
-       case COLAMD_ERROR_nnz_negative:
-
-           PRINTF(("Invalid number of nonzero entries (%d).\n", i1)) ;
-           break ;
-
-       case COLAMD_ERROR_p0_nonzero:
-
-           PRINTF(("Invalid column pointer, p [0] = %d, must be zero.\n", i1));
-           break ;
-
-       case COLAMD_ERROR_A_too_small:
-
-           PRINTF(("Array A too small.\n")) ;
-           PRINTF(("        Need Alen >= %d, but given only Alen = %d.\n",
-           i1, i2)) ;
-           break ;
-
-       case COLAMD_ERROR_col_length_negative:
-
-           PRINTF
-           (("Column %d has a negative number of nonzero entries (%d).\n",
-           INDEX (i1), i2)) ;
-           break ;
-
-       case COLAMD_ERROR_row_index_out_of_bounds:
-
-           PRINTF
-           (("Row index (row %d) out of bounds (%d to %d) in column %d.\n",
-           INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1))) ;
-           break ;
-
-       case COLAMD_ERROR_out_of_memory:
-
-           PRINTF(("Out of memory.\n")) ;
-           break ;
-
-       /* v2.4: internal-error case deleted */
-    }
-}
-
-
-
-
-/* ========================================================================== */
-/* === colamd debugging routines ============================================ */
-/* ========================================================================== */
-
-/* When debugging is disabled, the remainder of this file is ignored. */
-
-#ifndef NDEBUG
-
-
-/* ========================================================================== */
-/* === debug_structures ===================================================== */
-/* ========================================================================== */
-
-/*
-    At this point, all empty rows and columns are dead.  All live columns
-    are "clean" (containing no dead rows) and simplicial (no supercolumns
-    yet).  Rows may contain dead columns, but all live rows contain at
-    least one live column.
-*/
-
-PRIVATE void debug_structures
-(
-    /* === Parameters ======================================================= */
-
-    Int n_row,
-    Int n_col,
-    Colamd_Row Row [],
-    Colamd_Col Col [],
-    Int A [],
-    Int n_col2
-)
-{
-    /* === Local variables ================================================== */
-
-    Int i ;
-    Int c ;
-    Int *cp ;
-    Int *cp_end ;
-    Int len ;
-    Int score ;
-    Int r ;
-    Int *rp ;
-    Int *rp_end ;
-    Int deg ;
-
-    /* === Check A, Row, and Col ============================================ */
-
-    for (c = 0 ; c < n_col ; c++)
-    {
-       if (COL_IS_ALIVE (c))
-       {
-           len = Col [c].length ;
-           score = Col [c].shared2.score ;
-           DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ;
-           ASSERT (len > 0) ;
-           ASSERT (score >= 0) ;
-           ASSERT (Col [c].shared1.thickness == 1) ;
-           cp = &A [Col [c].start] ;
-           cp_end = cp + len ;
-           while (cp < cp_end)
-           {
-               r = *cp++ ;
-               ASSERT (ROW_IS_ALIVE (r)) ;
-           }
-       }
-       else
-       {
-           i = Col [c].shared2.order ;
-           ASSERT (i >= n_col2 && i < n_col) ;
-       }
-    }
-
-    for (r = 0 ; r < n_row ; r++)
-    {
-       if (ROW_IS_ALIVE (r))
-       {
-           i = 0 ;
-           len = Row [r].length ;
-           deg = Row [r].shared1.degree ;
-           ASSERT (len > 0) ;
-           ASSERT (deg > 0) ;
-           rp = &A [Row [r].start] ;
-           rp_end = rp + len ;
-           while (rp < rp_end)
-           {
-               c = *rp++ ;
-               if (COL_IS_ALIVE (c))
-               {
-                   i++ ;
-               }
-           }
-           ASSERT (i > 0) ;
-       }
-    }
-}
-
-
-/* ========================================================================== */
-/* === debug_deg_lists ====================================================== */
-/* ========================================================================== */
-
-/*
-    Prints the contents of the degree lists.  Counts the number of columns
-    in the degree list and compares it to the total it should have.  Also
-    checks the row degrees.
-*/
-
-PRIVATE void debug_deg_lists
-(
-    /* === Parameters ======================================================= */
-
-    Int n_row,
-    Int n_col,
-    Colamd_Row Row [],
-    Colamd_Col Col [],
-    Int head [],
-    Int min_score,
-    Int should,
-    Int max_deg
-)
-{
-    /* === Local variables ================================================== */
-
-    Int deg ;
-    Int col ;
-    Int have ;
-    Int row ;
-
-    /* === Check the degree lists =========================================== */
-
-    if (n_col > 10000 && colamd_debug <= 0)
-    {
-       return ;
-    }
-    have = 0 ;
-    DEBUG4 (("Degree lists: %d\n", min_score)) ;
-    for (deg = 0 ; deg <= n_col ; deg++)
-    {
-       col = head [deg] ;
-       if (col == EMPTY)
-       {
-           continue ;
-       }
-       DEBUG4 (("%d:", deg)) ;
-       while (col != EMPTY)
-       {
-           DEBUG4 ((" %d", col)) ;
-           have += Col [col].shared1.thickness ;
-           ASSERT (COL_IS_ALIVE (col)) ;
-           col = Col [col].shared4.degree_next ;
-       }
-       DEBUG4 (("\n")) ;
-    }
-    DEBUG4 (("should %d have %d\n", should, have)) ;
-    ASSERT (should == have) ;
-
-    /* === Check the row degrees ============================================ */
-
-    if (n_row > 10000 && colamd_debug <= 0)
-    {
-       return ;
-    }
-    for (row = 0 ; row < n_row ; row++)
-    {
-       if (ROW_IS_ALIVE (row))
-       {
-           ASSERT (Row [row].shared1.degree <= max_deg) ;
-       }
-    }
-}
-
-
-/* ========================================================================== */
-/* === debug_mark =========================================================== */
-/* ========================================================================== */
-
-/*
-    Ensures that the tag_mark is less that the maximum and also ensures that
-    each entry in the mark array is less than the tag mark.
-*/
-
-PRIVATE void debug_mark
-(
-    /* === Parameters ======================================================= */
-
-    Int n_row,
-    Colamd_Row Row [],
-    Int tag_mark,
-    Int max_mark
-)
-{
-    /* === Local variables ================================================== */
-
-    Int r ;
-
-    /* === Check the Row marks ============================================== */
-
-    ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
-    if (n_row > 10000 && colamd_debug <= 0)
-    {
-       return ;
-    }
-    for (r = 0 ; r < n_row ; r++)
-    {
-       ASSERT (Row [r].shared2.mark < tag_mark) ;
-    }
-}
-
-
-/* ========================================================================== */
-/* === debug_matrix ========================================================= */
-/* ========================================================================== */
-
-/*
-    Prints out the contents of the columns and the rows.
-*/
-
-PRIVATE void debug_matrix
-(
-    /* === Parameters ======================================================= */
-
-    Int n_row,
-    Int n_col,
-    Colamd_Row Row [],
-    Colamd_Col Col [],
-    Int A []
-)
-{
-    /* === Local variables ================================================== */
-
-    Int r ;
-    Int c ;
-    Int *rp ;
-    Int *rp_end ;
-    Int *cp ;
-    Int *cp_end ;
-
-    /* === Dump the rows and columns of the matrix ========================== */
-
-    if (colamd_debug < 3)
-    {
-       return ;
-    }
-    DEBUG3 (("DUMP MATRIX:\n")) ;
-    for (r = 0 ; r < n_row ; r++)
-    {
-       DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ;
-       if (ROW_IS_DEAD (r))
-       {
-           continue ;
-       }
-       DEBUG3 (("start %d length %d degree %d\n",
-               Row [r].start, Row [r].length, Row [r].shared1.degree)) ;
-       rp = &A [Row [r].start] ;
-       rp_end = rp + Row [r].length ;
-       while (rp < rp_end)
-       {
-           c = *rp++ ;
-           DEBUG4 (("  %d col %d\n", COL_IS_ALIVE (c), c)) ;
-       }
-    }
-
-    for (c = 0 ; c < n_col ; c++)
-    {
-       DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ;
-       if (COL_IS_DEAD (c))
-       {
-           continue ;
-       }
-       DEBUG3 (("start %d length %d shared1 %d shared2 %d\n",
-               Col [c].start, Col [c].length,
-               Col [c].shared1.thickness, Col [c].shared2.score)) ;
-       cp = &A [Col [c].start] ;
-       cp_end = cp + Col [c].length ;
-       while (cp < cp_end)
-       {
-           r = *cp++ ;
-           DEBUG4 (("  %d row %d\n", ROW_IS_ALIVE (r), r)) ;
-       }
-    }
-}
-
-PRIVATE void colamd_get_debug
-(
-    char *method
-)
-{
-    FILE *f ;
-    colamd_debug = 0 ;         /* no debug printing */
-    f = fopen ("debug", "r") ;
-    if (f == (FILE *) NULL)
-    {
-       colamd_debug = 0 ;
-    }
-    else
-    {
-       fscanf (f, "%d", &colamd_debug) ;
-       fclose (f) ;
-    }
-    DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n",
-       method, colamd_debug)) ;
-}
-
-#endif /* NDEBUG */
diff --git a/extern/colamd/Source/colamd_global.c b/extern/colamd/Source/colamd_global.c
deleted file mode 100644 (file)
index 4d1ae22..0000000
+++ /dev/null
@@ -1,24 +0,0 @@
-/* ========================================================================== */
-/* === colamd_global.c ====================================================== */
-/* ========================================================================== */
-
-/* ----------------------------------------------------------------------------
- * COLAMD, Copyright (C) 2007, Timothy A. Davis.
- * See License.txt for the Version 2.1 of the GNU Lesser General Public License
- * http://www.cise.ufl.edu/research/sparse
- * -------------------------------------------------------------------------- */
-
-/* Global variables for COLAMD */
-
-#ifndef NPRINT
-#ifdef MATLAB_MEX_FILE
-#include "mex.h"
-int (*colamd_printf) (const char *, ...) = mexPrintf ;
-#else
-#include <stdio.h>
-int (*colamd_printf) (const char *, ...) = printf ;
-#endif
-#else
-int (*colamd_printf) (const char *, ...) = ((void *) 0) ;
-#endif
-
index 9b0f2de22f2c5565fb48edb8936bab04de057410..275524f788fcd7175265db1ee21b6027e7e74d80 100644 (file)
@@ -74,10 +74,6 @@ if(WITH_BULLET)
        add_subdirectory(rigidbody)
 endif()
 
-if(WITH_OPENNL)
-       add_subdirectory(opennl)
-endif()
-
 if(WITH_OPENSUBDIV)
        add_subdirectory(opensubdiv)
 endif()
index 3b855d60ac883aa68324b870d4d7db86c627a5e3..736242102cc3a4b6baa10022f6a4eb865d1c27a9 100644 (file)
@@ -37,7 +37,6 @@ SConscript(['string/SConscript',
             'itasc/SConscript',
             'eigen/SConscript',
             'opencolorio/SConscript',
-            'opennl/SConscript',
             'mikktspace/SConscript',
             'smoke/SConscript',
             'raskter/SConscript'])
index 58964e432240143fe3b776abe9a2ab47c14326a9..5811b71de94d45db61c94d5780f9b02a6551b56a 100644 (file)
@@ -35,9 +35,11 @@ set(SRC
        eigen_capi.h
 
        intern/eigenvalues.cc
+       intern/linear_solver.cc
        intern/svd.cc
 
        intern/eigenvalues.h
+       intern/linear_solver.h
        intern/svd.h
 )
 
index 45ee1c015ece5e559aadeaca6ac4695e56f91508..be42e3402744d6174d60a0aebebef66a9476bb41 100644 (file)
@@ -28,6 +28,7 @@
 #define __EIGEN_C_API_H__
 
 #include "intern/eigenvalues.h"
+#include "intern/linear_solver.h"
 #include "intern/svd.h"
 
 #endif  /* __EIGEN_C_API_H__ */
index dcaaee8e9c2c850aa1abf96f8b0d6c11f2061a09..57942a4dc55014fcb81285b58a4d7e921a8a0a84 100644 (file)
@@ -45,7 +45,7 @@ using Eigen::Map;
 
 using Eigen::Success;
 
-bool EG3_self_adjoint_eigen_solve(const int size, const float *matrix, float *r_eigen_values, float *r_eigen_vectors)
+bool EIG_self_adjoint_eigen_solve(const int size, const float *matrix, float *r_eigen_values, float *r_eigen_vectors)
 {
        SelfAdjointEigenSolver<MatrixXf> eigen_solver;
 
index 93fc06c233906d62137d13c2d2e6209cc1ba1004..5c08ab5be39ec702b1a98f264be91693ae0d00d4 100644 (file)
@@ -31,7 +31,7 @@
 extern "C" {
 #endif
 
-bool EG3_self_adjoint_eigen_solve(const int size, const float *matrix, float *r_eigen_values, float *r_eigen_vectors);
+bool EIG_self_adjoint_eigen_solve(const int size, const float *matrix, float *r_eigen_values, float *r_eigen_vectors);
 
 #ifdef __cplusplus
 }
diff --git a/intern/eigen/intern/linear_solver.cc b/intern/eigen/intern/linear_solver.cc
new file mode 100644 (file)
index 0000000..181b278
--- /dev/null
@@ -0,0 +1,354 @@
+/*
+ *  Sparse linear solver.
+ *  Copyright (C) 2004 Bruno Levy
+ *  Copyright (C) 2005-2015 Blender Foundation
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 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 General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ *
+ *  If you modify this software, you should include a notice giving the
+ *  name of the person performing the modification, the date of modification,
+ *  and the reason for such modification.
+ */
+
+#include "linear_solver.h"
+
+#include <Eigen/Sparse>
+
+#include <algorithm>
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+
+/* Eigen data structures */
+
+typedef Eigen::SparseMatrix<double, Eigen::ColMajor> EigenSparseMatrix;
+typedef Eigen::SparseLU<EigenSparseMatrix> EigenSparseLU;
+typedef Eigen::VectorXd EigenVectorX;
+typedef Eigen::Triplet<double> EigenTriplet;
+
+/* Linear Solver data structure */
+
+struct LinearSolver
+{
+       struct Coeff
+       {
+               Coeff()
+               {
+                       index = 0;
+                       value = 0.0;
+               }
+
+               int index;
+               double value;
+       };
+
+       struct Variable
+       {
+               Variable()
+               {
+                       memset(value, 0, sizeof(value));
+                       locked = false;
+                       index = 0;
+               }
+
+               double value[4];
+               bool locked;
+               int index;
+               std::vector<Coeff> a;
+       };
+
+       enum State
+       {
+               STATE_VARIABLES_CONSTRUCT,
+               STATE_MATRIX_CONSTRUCT,
+               STATE_MATRIX_SOLVED
+       };
+
+       LinearSolver(int num_rows_, int num_variables_, int num_rhs_, bool lsq_)
+       {
+               assert(num_variables_ > 0);
+               assert(num_rhs_ <= 4);
+
+               state = STATE_VARIABLES_CONSTRUCT;
+               m = 0;
+               n = 0;
+               sparseLU = NULL;
+               num_variables = num_variables_;
+               num_rhs = num_rhs_;
+               num_rows = num_rows_;
+               least_squares = lsq_;
+
+               variable.resize(num_variables);
+       }
+
+       ~LinearSolver()
+       {
+               delete sparseLU;
+       }
+
+       State state;
+
+       int n;
+       int m;
+
+       std::vector<EigenTriplet> Mtriplets;
+       EigenSparseMatrix M;
+       EigenSparseMatrix MtM;
+       std::vector<EigenVectorX> b;
+       std::vector<EigenVectorX> x;
+
+       EigenSparseLU *sparseLU;
+
+       int num_variables;
+       std::vector<Variable> variable;
+
+       int num_rows;
+       int num_rhs;
+
+       bool least_squares;
+};
+
+LinearSolver *EIG_linear_solver_new(int num_rows, int num_columns, int num_rhs)
+{
+       return new LinearSolver(num_rows, num_columns, num_rhs, false);
+}
+
+LinearSolver *EIG_linear_least_squares_solver_new(int num_rows, int num_columns, int num_rhs)
+{
+       return new LinearSolver(num_rows, num_columns, num_rhs, true);
+}
+
+void EIG_linear_solver_delete(LinearSolver *solver)
+{
+       delete solver;
+}
+
+/* Variables */
+
+void EIG_linear_solver_variable_set(LinearSolver *solver, int rhs, int index, double value)
+{
+       solver->variable[index].value[rhs] = value;
+}
+
+double EIG_linear_solver_variable_get(LinearSolver *solver, int rhs, int index)
+{
+       return solver->variable[index].value[rhs];
+}
+
+void EIG_linear_solver_variable_lock(LinearSolver *solver, int index)
+{
+       if (!solver->variable[index].locked) {
+               assert(solver->state == LinearSolver::STATE_VARIABLES_CONSTRUCT);
+               solver->variable[index].locked = true;
+       }
+}
+
+static void linear_solver_variables_to_vector(LinearSolver *solver)
+{
+       int num_rhs = solver->num_rhs;
+
+       for (int i = 0; i < solver->num_variables; i++) {
+               LinearSolver::Variable* v = &solver->variable[i];
+               if (!v->locked) {
+                       for (int j = 0; j < num_rhs; j++)
+                               solver->x[j][v->index] = v->value[j];
+               }
+       }
+}
+
+static void linear_solver_vector_to_variables(LinearSolver *solver)
+{
+       int num_rhs = solver->num_rhs;
+
+       for (int i = 0; i < solver->num_variables; i++) {
+               LinearSolver::Variable* v = &solver->variable[i];
+               if (!v->locked) {
+                       for (int j = 0; j < num_rhs; j++)
+                               v->value[j] = solver->x[j][v->index];
+               }
+       }
+}
+
+/* Matrix */
+
+static void linear_solver_ensure_matrix_construct(LinearSolver *solver)
+{
+       /* transition to matrix construction if necessary */
+       if (solver->state == LinearSolver::STATE_VARIABLES_CONSTRUCT) {
+               int n = 0;
+
+               for (int i = 0; i < solver->num_variables; i++) {
+                       if (solver->variable[i].locked)
+                               solver->variable[i].index = ~0;
+                       else
+                               solver->variable[i].index = n++;
+               }
+
+               int m = (solver->num_rows == 0)? n: solver->num_rows;
+
+               solver->m = m;
+               solver->n = n;
+
+               assert(solver->least_squares || m == n);
+
+               /* reserve reasonable estimate */
+               solver->Mtriplets.clear();
+               solver->Mtriplets.reserve(std::max(m, n)*3);
+
+               solver->b.resize(solver->num_rhs);
+               solver->x.resize(solver->num_rhs);
+
+               for (int i = 0; i < solver->num_rhs; i++) {
+                       solver->b[i].setZero(m);
+                       solver->x[i].setZero(n);
+               }
+
+               linear_solver_variables_to_vector(solver);
+
+               solver->state = LinearSolver::STATE_MATRIX_CONSTRUCT;
+       }
+}
+
+void EIG_linear_solver_matrix_add(LinearSolver *solver, int row, int col, double value)
+{
+       if (solver->state == LinearSolver::STATE_MATRIX_SOLVED)
+               return;
+
+       linear_solver_ensure_matrix_construct(solver);
+
+       if (!solver->least_squares && solver->variable[row].locked);
+       else if (solver->variable[col].locked) {
+               if (!solver->least_squares)
+                       row = solver->variable[row].index;
+
+               LinearSolver::Coeff coeff;
+               coeff.index = row;
+               coeff.value = value;
+               solver->variable[col].a.push_back(coeff);
+       }
+       else {
+               if (!solver->least_squares)
+                       row = solver->variable[row].index;
+               col = solver->variable[col].index;
+
+               /* direct insert into matrix is too slow, so use triplets */
+               EigenTriplet triplet(row, col, value);
+               solver->Mtriplets.push_back(triplet);
+       }
+}
+
+/* Right hand side */
+
+void EIG_linear_solver_right_hand_side_add(LinearSolver *solver, int rhs, int index, double value)
+{
+       linear_solver_ensure_matrix_construct(solver);
+
+       if (solver->least_squares) {
+               solver->b[rhs][index] += value;
+       }
+       else if (!solver->variable[index].locked) {
+               index = solver->variable[index].index;
+               solver->b[rhs][index] += value;
+       }
+}
+
+/* Solve */
+
+bool EIG_linear_solver_solve(LinearSolver *solver)
+{
+       bool result = true;
+
+       assert(solver->state != LinearSolver::STATE_VARIABLES_CONSTRUCT);
+
+       if (solver->state == LinearSolver::STATE_MATRIX_CONSTRUCT) {
+               /* create matrix from triplets */
+               solver->M.resize(solver->m, solver->n);
+               solver->M.setFromTriplets(solver->Mtriplets.begin(), solver->Mtriplets.end());
+               solver->Mtriplets.clear();
+
+               /* create least squares matrix */
+               if (solver->least_squares)
+                       solver->MtM = solver->M.transpose() * solver->M;
+
+               /* convert M to compressed column format */
+               EigenSparseMatrix& M = (solver->least_squares)? solver->MtM: solver->M;
+               M.makeCompressed();
+
+               /* perform sparse LU factorization */
+               EigenSparseLU *sparseLU = new EigenSparseLU();
+               solver->sparseLU = sparseLU;
+
+               sparseLU->compute(M);
+               result = (sparseLU->info() == Eigen::Success);
+
+               solver->state = LinearSolver::STATE_MATRIX_SOLVED;
+       }
+
+       if (result) {
+               /* solve for each right hand side */
+               for (int rhs = 0; rhs < solver->num_rhs; rhs++) {
+                       /* modify for locked variables */
+                       EigenVectorX& b = solver->b[rhs];
+
+                       for (int i = 0; i < solver->num_variables; i++) {
+                               LinearSolver::Variable *variable = &solver->variable[i];
+
+                               if (variable->locked) {
+                                       std::vector<LinearSolver::Coeff>& a = variable->a;
+
+                                       for (int j = 0; j < a.size(); j++)
+                                               b[a[j].index] -= a[j].value*variable->value[rhs];
+                               }
+                       }
+
+                       /* solve */
+                       if (solver->least_squares) {
+                               EigenVectorX Mtb = solver->M.transpose() * b;
+                               solver->x[rhs] = solver->sparseLU->solve(Mtb);
+                       }
+                       else {
+                               EigenVectorX& b = solver->b[rhs];
+                               solver->x[rhs] = solver->sparseLU->solve(b);
+                       }
+
+                       if (solver->sparseLU->info() != Eigen::Success)
+                               result = false;
+               }
+
+               if (result)
+                       linear_solver_vector_to_variables(solver);
+       }
+
+       /* clear for next solve */
+       for (int rhs = 0; rhs < solver->num_rhs; rhs++)
+               solver->b[rhs].setZero(solver->m);
+
+       return result;
+}
+
+/* Debugging */
+
+void EIG_linear_solver_print_matrix(LinearSolver *solver)
+{
+       std::cout << "A:" << solver->M << std::endl;
+
+       for (int rhs = 0; rhs < solver->num_rhs; rhs++)
+               std::cout << "b " << rhs << ":" << solver->b[rhs] << std::endl;
+
+       if (solver->MtM.rows() && solver->MtM.cols())
+               std::cout << "AtA:" << solver->MtM << std::endl;
+}
+
diff --git a/intern/eigen/intern/linear_solver.h b/intern/eigen/intern/linear_solver.h
new file mode 100644 (file)
index 0000000..2dbea4d
--- /dev/null
@@ -0,0 +1,71 @@
+/*
+ *  Sparse linear solver.
+ *  Copyright (C) 2004 Bruno Levy
+ *  Copyright (C) 2005-2015 Blender Foundation
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 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 General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ *
+ *  If you modify this software, you should include a notice giving the
+ *  name of the person performing the modification, the date of modification,
+ *  and the reason for such modification.
+ */
+
+#pragma once
+
+#include <stdbool.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Solvers for Ax = b and AtAx = Atb */
+
+typedef struct LinearSolver LinearSolver;
+
+LinearSolver *EIG_linear_solver_new(
+       int num_rows,
+       int num_columns,
+       int num_right_hand_sides);
+
+LinearSolver *EIG_linear_least_squares_solver_new(
+       int num_rows,
+       int num_columns,
+       int num_right_hand_sides);
+
+void EIG_linear_solver_delete(LinearSolver *solver);
+
+/* Variables (x). Any locking must be done before matrix construction. */
+
+void EIG_linear_solver_variable_set(LinearSolver *solver, int rhs, int index, double value);
+double EIG_linear_solver_variable_get(LinearSolver *solver, int rhs, int index);
+void EIG_linear_solver_variable_lock(LinearSolver *solver, int index);
+
+/* Matrix (A) and right hand side (b) */
+
+void EIG_linear_solver_matrix_add(LinearSolver *solver, int row, int col, double value);
+void EIG_linear_solver_right_hand_side_add(LinearSolver *solver, int rhs, int index, double value);
+
+/* Solve. Repeated solves are supported, by changing b between solves. */
+
+bool EIG_linear_solver_solve(LinearSolver *solver);
+
+/* Debugging */
+
+void EIG_linear_solver_print_matrix(LinearSolver *solver);
+
+#ifdef __cplusplus
+}
+#endif
+
index e39a8261edb89990a1ecf57754fad45f71ccda57..4ed65af2a8f9c9c52a168bdf0207f6da3781fe5e 100644 (file)
@@ -48,7 +48,7 @@ using Eigen::MatrixXf;
 using Eigen::VectorXf;
 using Eigen::Map;
 
-void EG3_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V)
+void EIG_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V)
 {
        /* Since our matrix is squared, we can use thinU/V. */
        unsigned int flags = (r_U ? ComputeThinU : 0) | (r_V ? ComputeThinV : 0);
index 0ac5110897727e166f3b97f3c4f6d2f4a27dd24f..feadcc3520a6dac31857f7f3eca06cdbf75a1ba5 100644 (file)
@@ -31,7 +31,7 @@
 extern "C" {
 #endif
 
-void EG3_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V);
+void EIG_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V);
 
 #ifdef __cplusplus
 }
diff --git a/intern/opennl/CMakeLists.txt b/intern/opennl/CMakeLists.txt
deleted file mode 100644 (file)
index 9416bc2..0000000
+++ /dev/null
@@ -1,58 +0,0 @@
-# ***** BEGIN GPL LICENSE BLOCK *****
-#
-# This program is free software; you can redistribute it and/or
-# modify it under the terms of the GNU General Public License
-# as published by the Free Software Foundation; either version 2
-# 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 General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program; if not, write to the Free Software Foundation,
-# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-#
-# The Original Code is Copyright (C) 2006, Blender Foundation
-# All rights reserved.
-#
-# The Original Code is: all of this file.
-#
-# Contributor(s): Jacques Beaurain.
-#
-# ***** END GPL LICENSE BLOCK *****
-
-# External project, better not fix warnings.
-remove_strict_flags()
-
-# remove debug flag here since this is not a blender maintained library
-# and debug gives a lot of prints on UV unwrapping. developers can enable if they need to.
-if(MSVC)
-       remove_definitions(-DDEBUG)
-else()
-       add_definitions(-UDEBUG)
-endif()
-
-
-# quiet compiler warnings about undefined defines
-add_definitions(
-       -DDEBUGlevel=0
-       -DPRNTlevel=0
-)
-
-set(INC
-       extern
-)
-
-set(INC_SYS
-       ../../extern/colamd/Include
-       ../../extern/Eigen3
-)
-
-set(SRC
-       intern/opennl.cpp
-       extern/ONL_opennl.h
-)
-
-blender_add_lib(bf_intern_opennl "${SRC}" "${INC}" "${INC_SYS}")
diff --git a/intern/opennl/SConscript b/intern/opennl/SConscript
deleted file mode 100644 (file)
index 99df29b..0000000
+++ /dev/null
@@ -1,35 +0,0 @@
-#!/usr/bin/env python
-#
-# ***** BEGIN GPL LICENSE BLOCK *****
-#
-# This program is free software; you can redistribute it and/or
-# modify it under the terms of the GNU General Public License
-# as published by the Free Software Foundation; either version 2
-# 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 General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program; if not, write to the Free Software Foundation,
-# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-#
-# The Original Code is Copyright (C) 2006, Blender Foundation
-# All rights reserved.
-#
-# The Original Code is: all of this file.
-#
-# Contributor(s): Nathan Letwory.
-#
-# ***** END GPL LICENSE BLOCK *****
-
-Import ('env')
-
-sources = env.Glob('intern/*.cpp')
-
-incs = 'extern ../../extern/colamd/Include ../../extern/Eigen3'
-
-env.BlenderLib ('bf_intern_opennl', sources, Split(incs), [], libtype=['intern','player'], priority=[100,90] )
-
diff --git a/intern/opennl/extern/ONL_opennl.h b/intern/opennl/extern/ONL_opennl.h
deleted file mode 100644 (file)
index a414f40..0000000
+++ /dev/null
@@ -1,113 +0,0 @@
-/** \file opennl/extern/ONL_opennl.h
- *  \ingroup opennlextern
- */
-/*
- *  OpenNL: Numerical Library
- *  Copyright (C) 2004 Bruno Levy
- *
- *  This program is free software; you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation; either version 2 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 General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
- *
- *  If you modify this software, you should include a notice giving the
- *  name of the person performing the modification, the date of modification,
- *  and the reason for such modification.
- *
- *  Contact: Bruno Levy
- *
- *     levy@loria.fr
- *
- *     ISA Project
- *     LORIA, INRIA Lorraine, 
- *     Campus Scientifique, BP 239
- *     54506 VANDOEUVRE LES NANCY CEDEX 
- *     FRANCE
- *
- *  Note that the GNU General Public License does not permit incorporating
- *  the Software into proprietary programs. 
- */
-
-#ifndef nlOPENNL_H
-#define nlOPENNL_H
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-/* Datatypes */
-
-typedef unsigned int   NLenum;
-typedef unsigned char  NLboolean;
-typedef int                            NLint;          /* 4-byte signed */
-typedef unsigned int   NLuint;         /* 4-byte unsigned */
-typedef double                 NLdouble;       /* double precision float */
-
-typedef struct NLContext NLContext;
-
-/* Constants */
-
-#define NL_FALSE   0x0
-#define NL_TRUE    0x1
-
-/* Primitives */
-
-#define NL_SYSTEM  0x0
-#define NL_MATRIX  0x1
-
-/* Solver Parameters */
-
-#define NL_SOLVER              0x100
-#define NL_NB_VARIABLES        0x101
-#define NL_LEAST_SQUARES       0x102
-#define NL_ERROR               0x108
-#define NL_NB_ROWS             0x110
-#define NL_NB_RIGHT_HAND_SIDES 0x112 /* 4 max */
-
-/* Contexts */
-
-NLContext *nlNewContext(void);
-void nlDeleteContext(NLContext *context);
-
-/* State get/set */
-
-void nlSolverParameteri(NLContext *context, NLenum pname, NLint param);
-
-/* Variables */
-
-void nlSetVariable(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value);
-NLdouble nlGetVariable(NLContext *context, NLuint rhsindex, NLuint index);
-void nlLockVariable(NLContext *context, NLuint index);
-void nlUnlockVariable(NLContext *context, NLuint index);
-
-/* Begin/End */
-
-void nlBegin(NLContext *context, NLenum primitive);
-void nlEnd(NLContext *context, NLenum primitive);
-
-/* Setting elements in matrix/vector */
-
-void nlMatrixAdd(NLContext *context, NLuint row, NLuint col, NLdouble value);
-void nlRightHandSideAdd(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value);
-void nlRightHandSideSet(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value);
-
-/* Solve */
-
-void nlPrintMatrix(NLContext *context);
-NLboolean nlSolve(NLContext *context, NLboolean solveAgain);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
-
diff --git a/intern/opennl/intern/opennl.cpp b/intern/opennl/intern/opennl.cpp
deleted file mode 100644 (file)
index 331291e..0000000
+++ /dev/null
@@ -1,492 +0,0 @@
-/** \file opennl/intern/opennl.c
- *  \ingroup opennlintern
- */
-/*
- *
- *  OpenNL: Numerical Library
- *  Copyright (C) 2004 Bruno Levy
- *
- *  This program is free software; you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation; either version 2 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 General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
- *
- *  If you modify this software, you should include a notice giving the
- *  name of the person performing the modification, the date of modification,
- *  and the reason for such modification.
- *
- *  Contact: Bruno Levy
- *
- *      levy@loria.fr
- *
- *      ISA Project
- *      LORIA, INRIA Lorraine,
- *      Campus Scientifique, BP 239
- *      54506 VANDOEUVRE LES NANCY CEDEX
- *      FRANCE
- *
- *  Note that the GNU General Public License does not permit incorporating
- *  the Software into proprietary programs.
- */
-
-#include "ONL_opennl.h"
-
-#include <Eigen/Sparse>
-
-#include <algorithm>
-#include <cassert>
-#include <cstdlib>
-#include <iostream>
-#include <vector>
-
-/* Eigen data structures */
-
-typedef Eigen::SparseMatrix<double, Eigen::ColMajor> EigenSparseMatrix;
-typedef Eigen::SparseLU<EigenSparseMatrix> EigenSparseSolver;
-typedef Eigen::VectorXd EigenVectorX;
-typedef Eigen::Triplet<double> EigenTriplet;
-
-/* NLContext data structure */
-
-struct NLCoeff {
-       NLCoeff()
-       {
-               index = 0;
-               value = 0.0;
-       }
-       NLuint index;
-       NLdouble value;
-};
-
-struct NLVariable {
-       NLVariable()
-       {
-               memset(value, 0, sizeof(value));
-               locked = false;
-               index = 0;
-       }
-       NLdouble value[4];
-       NLboolean locked;
-       NLuint index;
-       std::vector<NLCoeff> a;
-};
-
-#define NL_STATE_INITIAL            0
-#define NL_STATE_SYSTEM             1
-#define NL_STATE_MATRIX             2
-#define NL_STATE_MATRIX_CONSTRUCTED 3
-#define NL_STATE_SYSTEM_CONSTRUCTED 4
-#define NL_STATE_SYSTEM_SOLVED      5
-
-struct NLContext {
-   NLContext()
-   {
-      state = NL_STATE_INITIAL;
-      n = 0;
-      m = 0;
-      sparse_solver = NULL;
-      nb_variables = 0;
-      nb_rhs = 1;
-      nb_rows = 0;
-      least_squares = false;
-      solve_again = false;
-   }
-
-   ~NLContext()
-   {
-      delete sparse_solver;
-   }
-
-       NLenum state;
-
-       NLuint n;
-       NLuint m;
-
-       std::vector<EigenTriplet> Mtriplets;
-       EigenSparseMatrix M;
-       EigenSparseMatrix MtM;
-       std::vector<EigenVectorX> b;
-       std::vector<EigenVectorX> Mtb;
-       std::vector<EigenVectorX> x;
-
-       EigenSparseSolver *sparse_solver;
-
-       NLuint nb_variables;
-       std::vector<NLVariable> variable;
-
-       NLuint nb_rows;
-       NLuint nb_rhs;
-
-       NLboolean least_squares;
-       NLboolean solve_again;
-};
-
-NLContext *nlNewContext(void)
-{
-       return new NLContext();
-}
-
-void nlDeleteContext(NLContext *context)
-{
-       delete context;
-}
-
-static void __nlCheckState(NLContext *context, NLenum state)
-{
-       assert(context->state == state);
-}
-
-static void __nlTransition(NLContext *context, NLenum from_state, NLenum to_state)
-{
-       __nlCheckState(context, from_state);
-       context->state = to_state;
-}
-
-/* Get/Set parameters */
-
-void nlSolverParameteri(NLContext *context, NLenum pname, NLint param)
-{
-       __nlCheckState(context, NL_STATE_INITIAL);
-       switch(pname) {
-       case NL_NB_VARIABLES: {
-               assert(param > 0);
-               context->nb_variables = (NLuint)param;
-       } break;
-       case NL_NB_ROWS: {
-               assert(param > 0);
-               context->nb_rows = (NLuint)param;
-       } break;
-       case NL_LEAST_SQUARES: {
-               context->least_squares = (NLboolean)param;
-       } break;
-       case NL_NB_RIGHT_HAND_SIDES: {
-               context->nb_rhs = (NLuint)param;
-       } break;
-       default: {
-               assert(0);
-       } break;
-       }
-}
-
-/* Get/Set Lock/Unlock variables */
-
-void nlSetVariable(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value)
-{
-       __nlCheckState(context, NL_STATE_SYSTEM);
-       context->variable[index].value[rhsindex] = value;
-}
-
-NLdouble nlGetVariable(NLContext *context, NLuint rhsindex, NLuint index)
-{
-       assert(context->state != NL_STATE_INITIAL);
-       return context->variable[index].value[rhsindex];
-}
-
-void nlLockVariable(NLContext *context, NLuint index)
-{
-       __nlCheckState(context, NL_STATE_SYSTEM);
-       context->variable[index].locked = true;
-}
-
-void nlUnlockVariable(NLContext *context, NLuint index)
-{
-       __nlCheckState(context, NL_STATE_SYSTEM);
-       context->variable[index].locked = false;
-}
-
-/* System construction */
-
-static void __nlVariablesToVector(NLContext *context)
-{
-       NLuint i, j, nb_rhs;
-
-       nb_rhs= context->nb_rhs;
-
-       for(i=0; i<context->nb_variables; i++) {
-               NLVariable* v = &(context->variable[i]);
-               if(!v->locked) {
-                       for(j=0; j<nb_rhs; j++)
-                               context->x[j][v->index] = v->value[j];
-               }
-       }
-}
-
-static void __nlVectorToVariables(NLContext *context)
-{
-       NLuint i, j, nb_rhs;
-
-       nb_rhs= context->nb_rhs;
-
-       for(i=0; i<context->nb_variables; i++) {
-               NLVariable* v = &(context->variable[i]);
-               if(!v->locked) {
-                       for(j=0; j<nb_rhs; j++)
-                               v->value[j] = context->x[j][v->index];
-               }
-       }
-}
-
-static void __nlBeginSystem(NLContext *context)
-{
-       assert(context->nb_variables > 0);
-
-       if (context->solve_again)
-               __nlTransition(context, NL_STATE_SYSTEM_SOLVED, NL_STATE_SYSTEM);
-       else {
-               __nlTransition(context, NL_STATE_INITIAL, NL_STATE_SYSTEM);
-
-               context->variable.resize(context->nb_variables);
-       }
-}
-
-static void __nlEndSystem(NLContext *context)
-{
-       __nlTransition(context, NL_STATE_MATRIX_CONSTRUCTED, NL_STATE_SYSTEM_CONSTRUCTED);
-}
-
-static void __nlBeginMatrix(NLContext *context)
-{
-       NLuint i;
-       NLuint m = 0, n = 0;
-
-       __nlTransition(context, NL_STATE_SYSTEM, NL_STATE_MATRIX);
-
-       if (!context->solve_again) {
-               for(i=0; i<context->nb_variables; i++) {
-                       if(context->variable[i].locked)
-                               context->variable[i].index = ~0;
-                       else
-                               context->variable[i].index = n++;
-               }
-
-               m = (context->nb_rows == 0)? n: context->nb_rows;
-
-               context->m = m;
-               context->n = n;
-
-               /* reserve reasonable estimate */
-               context->Mtriplets.clear();
-               context->Mtriplets.reserve(std::max(m, n)*3);
-
-               context->b.resize(context->nb_rhs);
-               context->x.resize(context->nb_rhs);
-
-               for (i=0; i<context->nb_rhs; i++) {
-                       context->b[i].setZero(m);
-                       context->x[i].setZero(n);
-               }
-       }
-       else {
-               /* need to recompute b only, A is not constructed anymore */
-               for (i=0; i<context->nb_rhs; i++)
-                       context->b[i].setZero(context->m);
-       }
-
-       __nlVariablesToVector(context);
-}
-
-static void __nlEndMatrixRHS(NLContext *context, NLuint rhs)
-{
-       NLVariable *variable;
-       NLuint i, j;
-
-       EigenVectorX& b = context->b[rhs];
-
-       for(i=0; i<context->nb_variables; i++) {
-               variable = &(context->variable[i]);
-
-               if(variable->locked) {
-                       std::vector<NLCoeff>& a = variable->a;
-
-                       for(j=0; j<a.size(); j++) {
-                               b[a[j].index] -= a[j].value*variable->value[rhs];
-                       }
-               }
-       }
-
-       if(context->least_squares)
-               context->Mtb[rhs] = context->M.transpose() * b;
-}
-
-static void __nlEndMatrix(NLContext *context)
-{
-       __nlTransition(context, NL_STATE_MATRIX, NL_STATE_MATRIX_CONSTRUCTED);
-
-       if(!context->solve_again) {
-               context->M.resize(context->m, context->n);
-               context->M.setFromTriplets(context->Mtriplets.begin(), context->Mtriplets.end());
-               context->Mtriplets.clear();
-
-               if(context->least_squares) {
-                       context->MtM = context->M.transpose() * context->M;
-
-                       context->Mtb.resize(context->nb_rhs);
-                       for (NLuint rhs=0; rhs<context->nb_rhs; rhs++)
-                               context->Mtb[rhs].setZero(context->n);
-               }
-       }
-
-       for (NLuint rhs=0; rhs<context->nb_rhs; rhs++)
-               __nlEndMatrixRHS(context, rhs);
-}
-
-void nlMatrixAdd(NLContext *context, NLuint row, NLuint col, NLdouble value)
-{
-       __nlCheckState(context, NL_STATE_MATRIX);
-
-       if(context->solve_again)
-               return;
-
-       if (!context->least_squares && context->variable[row].locked);
-       else if (context->variable[col].locked) {
-               if(!context->least_squares)
-                       row = context->variable[row].index;
-
-               NLCoeff coeff;
-               coeff.index = row;
-               coeff.value = value;
-               context->variable[col].a.push_back(coeff);
-       }
-       else {
-               if(!context->least_squares)
-                       row = context->variable[row].index;
-               col = context->variable[col].index;
-
-               // direct insert into matrix is too slow, so use triplets
-               EigenTriplet triplet(row, col, value);
-               context->Mtriplets.push_back(triplet);
-       }
-}
-
-void nlRightHandSideAdd(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value)
-{
-       __nlCheckState(context, NL_STATE_MATRIX);
-
-       if(context->least_squares) {
-               context->b[rhsindex][index] += value;
-       }
-       else {
-               if(!context->variable[index].locked) {
-                       index = context->variable[index].index;
-                       context->b[rhsindex][index] += value;
-               }
-       }
-}
-
-void nlRightHandSideSet(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value)
-{
-       __nlCheckState(context, NL_STATE_MATRIX);
-
-       if(context->least_squares) {
-               context->b[rhsindex][index] = value;
-       }
-       else {
-               if(!context->variable[index].locked) {
-                       index = context->variable[index].index;
-                       context->b[rhsindex][index] = value;
-               }
-       }
-}
-
-void nlBegin(NLContext *context, NLenum prim)
-{
-       switch(prim) {
-       case NL_SYSTEM: {
-               __nlBeginSystem(context);
-       } break;
-       case NL_MATRIX: {
-               __nlBeginMatrix(context);
-       } break;
-       default: {
-               assert(0);
-       }
-       }
-}
-
-void nlEnd(NLContext *context, NLenum prim)
-{
-       switch(prim) {
-       case NL_SYSTEM: {
-               __nlEndSystem(context);
-       } break;
-       case NL_MATRIX: {
-               __nlEndMatrix(context);
-       } break;
-       default: {
-               assert(0);
-       }
-       }
-}
-
-void nlPrintMatrix(NLContext *context)
-{
-       std::cout << "A:" << context->M << std::endl;
-
-       for(NLuint rhs=0; rhs<context->nb_rhs; rhs++)
-               std::cout << "b " << rhs << ":" << context->b[rhs] << std::endl;
-
-       if (context->MtM.rows() && context->MtM.cols())
-               std::cout << "AtA:" << context->MtM << std::endl;
-}
-
-/* Solving */
-
-NLboolean nlSolve(NLContext *context, NLboolean solveAgain)
-{
-       NLboolean result = true;
-
-       __nlCheckState(context, NL_STATE_SYSTEM_CONSTRUCTED);
-
-       if (!context->solve_again) {
-               EigenSparseMatrix& M = (context->least_squares)? context->MtM: context->M;
-
-               assert(M.rows() == M.cols());
-
-               /* Convert M to compressed column format */
-               M.makeCompressed();
-
-               /* Perform sparse LU factorization */
-               EigenSparseSolver *sparse_solver = new EigenSparseSolver();
-               context->sparse_solver = sparse_solver;
-
-               sparse_solver->analyzePattern(M);
-               sparse_solver->factorize(M);
-
-               result = (sparse_solver->info() == Eigen::Success);
-
-               /* Free M, don't need it anymore at this point */
-               M.resize(0, 0);
-       }
-
-       if (result) {
-               /* Solve each right hand side */
-               for(NLuint rhs=0; rhs<context->nb_rhs; rhs++) {
-                       EigenVectorX& b = (context->least_squares)? context->Mtb[rhs]: context->b[rhs];
-                       context->x[rhs] = context->sparse_solver->solve(b);
-
-                       if (context->sparse_solver->info() != Eigen::Success)
-                               result = false;
-               }
-
-               if (result) {
-                       __nlVectorToVariables(context);
-
-                       if (solveAgain)
-                               context->solve_again = true;
-
-                       __nlTransition(context, NL_STATE_SYSTEM_CONSTRUCTED, NL_STATE_SYSTEM_SOLVED);
-               }
-       }
-
-       return result;
-}
-
index 86ffb541bf8e87266e76de33897bf9d05e2b1a60..1c3ff3d3e0262c52ecd7d4824abcedc042910c3b 100644 (file)
@@ -81,7 +81,6 @@ variables on the UI for now
 #include "BKE_scene.h"
 
 #include  "PIL_time.h"
-// #include  "ONL_opennl.h" remove linking to ONL for now
 
 /* callbacks for errors and interrupts and some goo */
 static int (*SB_localInterruptCallBack)(void) = NULL;
@@ -1811,14 +1810,14 @@ static void dfdx_spring(int ia, int ic, int op, float dir[3], float L, float len
                        for (j=0;j<3;j++) {
                                delta_ij = (i==j ? (1.0f): (0.0f));
                                m=factor*(dir[i]*dir[j] + (1-L/len)*(delta_ij - dir[i]*dir[j]));
-                               nlMatrixAdd(ia+i, op+ic+j, m);
+                               EIG_linear_solver_matrix_add(ia+i, op+ic+j, m);
                        }
        }
        else {
                for (i=0;i<3;i++)
                        for (j=0;j<3;j++) {
                                m=factor*dir[i]*dir[j];
-                               nlMatrixAdd(ia+i, op+ic+j, m);
+                               EIG_linear_solver_matrix_add(ia+i, op+ic+j, m);
                        }
        }
 }
@@ -1827,13 +1826,13 @@ static void dfdx_spring(int ia, int ic, int op, float dir[3], float L, float len
 static void dfdx_goal(int ia, int ic, int op, float factor)
 {
        int i;
-       for (i=0;i<3;i++) nlMatrixAdd(ia+i, op+ic+i, factor);
+       for (i=0;i<3;i++) EIG_linear_solver_matrix_add(ia+i, op+ic+i, factor);
 }
 
 static void dfdv_goal(int ia, int ic, float factor)
 {
        int i;
-       for (i=0;i<3;i++) nlMatrixAdd(ia+i, ic+i, factor);
+       for (i=0;i<3;i++) EIG_linear_solver_matrix_add(ia+i, ic+i, factor);
 }
 */
 static void sb_spring_force(Object *ob, int bpi, BodySpring *bs, float iks, float UNUSED(forcetime))
index c27c6bea160233fdab2b73c818db675c62c113bb..641e50c9bdecf0c10619ed2856c6efc131b136b3 100644 (file)
@@ -57,7 +57,7 @@ bool BLI_eigen_solve_selfadjoint_m3(const float m3[3][3], float r_eigen_values[3
        }
 #endif
 
-       return EG3_self_adjoint_eigen_solve(3, (const float *)m3, r_eigen_values, (float *)r_eigen_vectors);
+       return EIG_self_adjoint_eigen_solve(3, (const float *)m3, r_eigen_values, (float *)r_eigen_vectors);
 }
 
 /**
@@ -70,5 +70,5 @@ bool BLI_eigen_solve_selfadjoint_m3(const float m3[3][3], float r_eigen_values[3
  */
 void BLI_svd_m3(const float m3[3][3], float r_U[3][3], float r_S[3], float r_V[3][3])
 {
-       EG3_svd_square_matrix(3, (const float *)m3, (float *)r_U, (float *)r_S, (float *)r_V);
+       EIG_svd_square_matrix(3, (const float *)m3, (float *)r_U, (float *)r_S, (float *)r_V);
 }
index 686b8507291732769e2f10a7fb52c66e16858d52..d13d3e1cd0fbab261fef680d70c06c166f944aff 100644 (file)
@@ -30,6 +30,7 @@ set(INC
        ../blentranslation
        ../makesdna
        ../../../intern/guardedalloc
+       ../../../intern/eigen
        ../../../extern/rangetree
 )
 
@@ -176,13 +177,6 @@ if(WITH_INTERNATIONAL)
        add_definitions(-DWITH_INTERNATIONAL)
 endif()
 
-if(WITH_OPENNL)
-       add_definitions(-DWITH_OPENNL)
-       list(APPEND INC_SYS
-               ../../../intern/opennl/extern
-       )
-endif()
-
 if(WITH_FREESTYLE)
        add_definitions(-DWITH_FREESTYLE)
 endif()
index c53974be1a77686f254f8131826da24f52d8ca78..6bf4fdf7c7872dd170519083de61cae0d8ec2d55 100644 (file)
@@ -40,9 +40,9 @@ incs = [
     '../makesdna',
     '../blenkernel',
     '#/intern/guardedalloc',
+    '#/intern/eigen',
     '#/extern/bullet2/src',
     '#/extern/rangetree',
-    '#/intern/opennl/extern'
     ]
 
 defs = []
index 624fbd938183ba7225044b1cb0663c1e7f7bde56..95ff3eb48588d740f5c647904cff3e45014148cb 100644 (file)
 
 #include "MEM_guardedalloc.h"
 
-
 #include "BLI_math.h"
 
+#include "eigen_capi.h"
+
 
 #include "bmesh.h"
 
 #include "intern/bmesh_operators_private.h" /* own include */
 
-#ifdef WITH_OPENNL
-
-#include "ONL_opennl.h"
-
 // #define SMOOTH_LAPLACIAN_AREA_FACTOR 4.0f  /* UNUSED */
 // #define SMOOTH_LAPLACIAN_EDGE_FACTOR 2.0f  /* UNUSED */
 #define SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE 1.8f
@@ -59,7 +56,7 @@ struct BLaplacianSystem {
        /* Pointers to data*/
        BMesh *bm;
        BMOperator *op;
-       NLContext *context;
+       LinearSolver *context;
 
        /*Data*/
        float min_area;
@@ -92,7 +89,7 @@ static void delete_laplacian_system(LaplacianSystem *sys)
        delete_void_pointer(sys->vweights);
        delete_void_pointer(sys->zerola);
        if (sys->context) {
-               nlDeleteContext(sys->context);
+               EIG_linear_solver_delete(sys->context);
        }
        sys->bm = NULL;
        sys->op = NULL;
@@ -333,9 +330,9 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
                                        w4 = w4 / 4.0f;
 
                                        if (!vert_is_boundary(vf[j]) && sys->zerola[idv1] == 0) {
-                                               nlMatrixAdd(sys->context, idv1, idv2, w2 * sys->vweights[idv1]);
-                                               nlMatrixAdd(sys->context, idv1, idv3, w3 * sys->vweights[idv1]);
-                                               nlMatrixAdd(sys->context, idv1, idv4, w4 * sys->vweights[idv1]);
+                                               EIG_linear_solver_matrix_add(sys->context, idv1, idv2, w2 * sys->vweights[idv1]);
+                                               EIG_linear_solver_matrix_add(sys->context, idv1, idv3, w3 * sys->vweights[idv1]);
+                                               EIG_linear_solver_matrix_add(sys->context, idv1, idv4, w4 * sys->vweights[idv1]);
                                        }
                                }
                        }
@@ -346,16 +343,16 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
                                /* Is ring if number of faces == number of edges around vertice*/
                                i = BM_elem_index_get(f);
                                if (!vert_is_boundary(vf[0]) && sys->zerola[idv1] == 0) {
-                                       nlMatrixAdd(sys->context, idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
-                                       nlMatrixAdd(sys->context, idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
+                                       EIG_linear_solver_matrix_add(sys->context, idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
+                                       EIG_linear_solver_matrix_add(sys->context, idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
                                }
                                if (!vert_is_boundary(vf[1]) && sys->zerola[idv2] == 0) {
-                                       nlMatrixAdd(sys->context, idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
-                                       nlMatrixAdd(sys->context, idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
+                                       EIG_linear_solver_matrix_add(sys->context, idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
+                                       EIG_linear_solver_matrix_add(sys->context, idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
                                }
                                if (!vert_is_boundary(vf[2]) && sys->zerola[idv3] == 0) {
-                                       nlMatrixAdd(sys->context, idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
-                                       nlMatrixAdd(sys->context, idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
+                                       EIG_linear_solver_matrix_add(sys->context, idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
+                                       EIG_linear_solver_matrix_add(sys->context, idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
                                }
                        }
                }
@@ -368,8 +365,8 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
                        idv2 = BM_elem_index_get(e->v2);
                        if (sys->zerola[idv1] == 0 && sys->zerola[idv2] == 0) {
                                i = BM_elem_index_get(e);
-                               nlMatrixAdd(sys->context, idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
-                               nlMatrixAdd(sys->context, idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
+                               EIG_linear_solver_matrix_add(sys->context, idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
+                               EIG_linear_solver_matrix_add(sys->context, idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
                        }
                }
        }
@@ -434,12 +431,12 @@ static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez
                idv2 = BM_elem_index_get(e->v2);
                vi1 = e->v1->co;
                vi2 =  e->v2->co;
-               ve1[0] = nlGetVariable(sys->context, 0, idv1);
-               ve1[1] = nlGetVariable(sys->context, 1, idv1);
-               ve1[2] = nlGetVariable(sys->context, 2, idv1);
-               ve2[0] = nlGetVariable(sys->context, 0, idv2);
-               ve2[1] = nlGetVariable(sys->context, 1, idv2);
-               ve2[2] = nlGetVariable(sys->context, 2, idv2);
+               ve1[0] = EIG_linear_solver_variable_get(sys->context, 0, idv1);
+               ve1[1] = EIG_linear_solver_variable_get(sys->context, 1, idv1);
+               ve1[2] = EIG_linear_solver_variable_get(sys->context, 2, idv1);
+               ve2[0] = EIG_linear_solver_variable_get(sys->context, 0, idv2);
+               ve2[1] = EIG_linear_solver_variable_get(sys->context, 1, idv2);
+               ve2[2] = EIG_linear_solver_variable_get(sys->context, 2, idv2);
                leni = len_v3v3(vi1, vi2);
                lene = len_v3v3(ve1, ve2);
                if (lene > leni * SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE || lene < leni * SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE) {
@@ -455,13 +452,13 @@ static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez
                m_vertex_id = BM_elem_index_get(v);
                if (sys->zerola[m_vertex_id] == 0) {
                        if (usex) {
-                               v->co[0] =  nlGetVariable(sys->context, 0, m_vertex_id);
+                               v->co[0] =  EIG_linear_solver_variable_get(sys->context, 0, m_vertex_id);
                        }
                        if (usey) {
-                               v->co[1] =  nlGetVariable(sys->context, 1, m_vertex_id);
+                               v->co[1] =  EIG_linear_solver_variable_get(sys->context, 1, m_vertex_id);
                        }
                        if (usez) {
-                               v->co[2] =  nlGetVariable(sys->context, 2, m_vertex_id);
+                               v->co[2] =  EIG_linear_solver_variable_get(sys->context, 2, m_vertex_id);
                        }
                }
        }
@@ -501,32 +498,24 @@ void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op)
        preserve_volume = BMO_slot_bool_get(op->slots_in, "preserve_volume");
 
 
-       sys->context = nlNewContext();
-
-       nlSolverParameteri(sys->context, NL_NB_VARIABLES, bm->totvert);
-       nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
-       nlSolverParameteri(sys->context, NL_NB_ROWS, bm->totvert);
-       nlSolverParameteri(sys->context, NL_NB_RIGHT_HAND_SIDES, 3);
+       sys->context = EIG_linear_least_squares_solver_new(bm->totvert, bm->totvert, 3);
 
-       nlBegin(sys->context, NL_SYSTEM);
        for (i = 0; i < bm->totvert; i++) {
-               nlLockVariable(sys->context, i);
+               EIG_linear_solver_variable_lock(sys->context, i);
        }
        BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
                m_vertex_id = BM_elem_index_get(v);
-               nlUnlockVariable(sys->context, m_vertex_id);
-               nlSetVariable(sys->context, 0, m_vertex_id, v->co[0]);
-               nlSetVariable(sys->context, 1, m_vertex_id, v->co[1]);
-               nlSetVariable(sys->context, 2, m_vertex_id, v->co[2]);
+               EIG_linear_solver_variable_set(sys->context, 0, m_vertex_id, v->co[0]);
+               EIG_linear_solver_variable_set(sys->context, 1, m_vertex_id, v->co[1]);
+               EIG_linear_solver_variable_set(sys->context, 2, m_vertex_id, v->co[2]);
        }
 
-       nlBegin(sys->context, NL_MATRIX);
        init_laplacian_matrix(sys);
        BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
                m_vertex_id = BM_elem_index_get(v);
-               nlRightHandSideAdd(sys->context, 0, m_vertex_id, v->co[0]);
-               nlRightHandSideAdd(sys->context, 1, m_vertex_id, v->co[1]);
-               nlRightHandSideAdd(sys->context, 2, m_vertex_id, v->co[2]);
+               EIG_linear_solver_right_hand_side_add(sys->context, 0, m_vertex_id, v->co[0]);
+               EIG_linear_solver_right_hand_side_add(sys->context, 1, m_vertex_id, v->co[1]);
+               EIG_linear_solver_right_hand_side_add(sys->context, 2, m_vertex_id, v->co[2]);
                i = m_vertex_id;
                if (sys->zerola[i] == 0) {
                        w = sys->vweights[i] * sys->ring_areas[i];
@@ -535,34 +524,22 @@ void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op)
                        sys->vlengths[i] = (w == 0.0f) ? 0.0f : -lambda_border  * 2.0f / w;
 
                        if (!vert_is_boundary(v)) {
-                               nlMatrixAdd(sys->context, i, i,  1.0f + lambda_factor / (4.0f * sys->ring_areas[i]));
+                               EIG_linear_solver_matrix_add(sys->context, i, i,  1.0f + lambda_factor / (4.0f * sys->ring_areas[i]));
                        }
                        else {
-                               nlMatrixAdd(sys->context, i, i,  1.0f + lambda_border * 2.0f);
+                               EIG_linear_solver_matrix_add(sys->context, i, i,  1.0f + lambda_border * 2.0f);
                        }
                }
                else {
-                       nlMatrixAdd(sys->context, i, i, 1.0f);
+                       EIG_linear_solver_matrix_add(sys->context, i, i, 1.0f);
                }
        }
        fill_laplacian_matrix(sys);
 
-       nlEnd(sys->context, NL_MATRIX);
-       nlEnd(sys->context, NL_SYSTEM);
-
-       if (nlSolve(sys->context, NL_TRUE) ) {
+       if (EIG_linear_solver_solve(sys->context) ) {
                validate_solution(sys, usex, usey, usez, preserve_volume);
        }
 
        delete_laplacian_system(sys);
 }
 
-#else  /* WITH_OPENNL */
-
-#ifdef __GNUC__
-#  pragma GCC diagnostic ignored "-Wunused-parameter"
-#endif
-
-void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op) {}
-
-#endif  /* WITH_OPENNL */
index 1ed70b3cd986a20023f13401c1ff5f037127849b..b213aca478f01bde550891241b08819d161bd7ba 100644 (file)
@@ -28,6 +28,7 @@ set(INC
        ../../makesrna
        ../../windowmanager
        ../../../../intern/guardedalloc
+       ../../../../intern/eigen
        ../../../../intern/glew-mx
 )
 
@@ -68,13 +69,6 @@ if(WITH_INTERNATIONAL)
        add_definitions(-DWITH_INTERNATIONAL)
 endif()
 
-if(WITH_OPENNL)
-       add_definitions(-DWITH_OPENNL)
-       list(APPEND INC_SYS
-               ../../../../intern/opennl/extern
-       )
-endif()
-
 add_definitions(${GL_DEFINITIONS})
 
 blender_add_lib(bf_editor_armature "${SRC}" "${INC}" "${INC_SYS}")
index 9c3959ecb5be4e0aa6b08e1ce5152cef958db20e..850834d2cd4e2b18ab21b0683e26d6fefc28b4a9 100644 (file)
@@ -31,9 +31,9 @@ sources = env.Glob('*.c')
 
 incs = [
     '#/intern/guardedalloc',
+    '#/intern/eigen',
     env['BF_GLEW_INC'],
     '#/intern/glew-mx',
-    '#/intern/opennl/extern',
     '../include',
     '../../blenkernel',
     '../../blenlib',
index ea1a94fbba6323b3382ebe80afcd0387e33590c6..c621d6f99c0051d46083a3b347e5214f847d96a0 100644 (file)
 #include "ED_armature.h"
 #include "ED_mesh.h"
 
+#include "eigen_capi.h"
 
 #include "armature_intern.h"
-
-#ifdef WITH_OPENNL
-#  include "meshlaplacian.h"
-#endif
+#include "meshlaplacian.h"
 
 #if 0
 #include "reeb.h"
@@ -401,12 +399,8 @@ static void add_verts_to_dgroups(ReportList *reports, Scene *scene, Object *ob,
        if (heat) {
                const char *error = NULL;
 
-#ifdef WITH_OPENNL
                heat_bone_weighting(ob, mesh, verts, numbones, dgrouplist, dgroupflip,
                                    root, tip, selected, &error);
-#else
-               error = "Built without OpenNL";
-#endif
                if (error) {
                        BKE_report(reports, RPT_WARNING, error);
                }
index 3f27a58930d204a6ffd8b9b8472ab1d98afe1a8a..6bf75413e81c69a4788b14351ef237c6bb093c08 100644 (file)
 #include "ED_mesh.h"
 #include "ED_armature.h"
 
-#include "meshlaplacian.h"
-
-#ifdef WITH_OPENNL
+#include "eigen_capi.h"
 
-#include "ONL_opennl.h"
+#include "meshlaplacian.h"
 
 /* ************* XXX *************** */
 static void waitcursor(int UNUSED(val)) {}
@@ -64,7 +62,7 @@ static void error(const char *str) { printf("error: %s\n", str); }
 /************************** Laplacian System *****************************/
 
 struct LaplacianSystem {
-       NLContext *context;  /* opennl context */
+       LinearSolver *context;  /* linear solver */
 
        int totvert, totface;
 
@@ -76,7 +74,7 @@ struct LaplacianSystem {
 
        int areaweights;        /* use area in cotangent weights? */
        int storeweights;       /* store cotangent weights in fweights */
-       int nlbegun;            /* nlBegin(NL_SYSTEM/NL_MATRIX) done */
+       bool variablesdone;     /* variables set in linear system */
 
        EdgeHash *edgehash;     /* edge hash for construction */
 
@@ -182,18 +180,18 @@ static void laplacian_triangle_weights(LaplacianSystem *sys, int f, int i1, int
        t2 = cotangent_tri_weight_v3(v2, v3, v1) / laplacian_edge_count(sys->edgehash, i3, i1);
        t3 = cotangent_tri_weight_v3(v3, v1, v2) / laplacian_edge_count(sys->edgehash, i1, i2);
 
-       nlMatrixAdd(sys->context, i1, i1, (t2 + t3) * varea[i1]);
-       nlMatrixAdd(sys->context, i2, i2, (t1 + t3) * varea[i2]);
-       nlMatrixAdd(sys->context, i3, i3, (t1 + t2) * varea[i3]);
+       EIG_linear_solver_matrix_add(sys->context, i1, i1, (t2 + t3) * varea[i1]);
+       EIG_linear_solver_matrix_add(sys->context, i2, i2, (t1 + t3) * varea[i2]);
+       EIG_linear_solver_matrix_add(sys->context, i3, i3, (t1 + t2) * varea[i3]);
 
-       nlMatrixAdd(sys->context, i1, i2, -t3 * varea[i1]);
-       nlMatrixAdd(sys->context, i2, i1, -t3 * varea[i2]);
+       EIG_linear_solver_matrix_add(sys->context, i1, i2, -t3 * varea[i1]);
+       EIG_linear_solver_matrix_add(sys->context, i2, i1, -t3 * varea[i2]);
 
-       nlMatrixAdd(sys->context, i2, i3, -t1 * varea[i2]);
-       nlMatrixAdd(sys->context, i3, i2, -t1 * varea[i3]);
+       EIG_linear_solver_matrix_add(sys->context, i2, i3, -t1 * varea[i2]);
+       EIG_linear_solver_matrix_add(sys->context, i3, i2, -t1 * varea[i3]);
 
-       nlMatrixAdd(sys->context, i3, i1, -t2 * varea[i3]);
-       nlMatrixAdd(sys->context, i1, i3, -t2 * varea[i1]);
+       EIG_linear_solver_matrix_add(sys->context, i3, i1, -t2 * varea[i3]);
+       EIG_linear_solver_matrix_add(sys->context, i1, i3, -t2 * varea[i1]);
 
        if (sys->storeweights) {
                sys->fweights[f][0] = t1 * varea[i1];
@@ -218,11 +216,11 @@ static LaplacianSystem *laplacian_system_construct_begin(int totvert, int totfac
        sys->areaweights = 1;
        sys->storeweights = 0;
 
-       /* create opennl context */
-       sys->context = nlNewContext();
-       nlSolverParameteri(sys->context, NL_NB_VARIABLES, totvert);
+       /* create linear solver */
        if (lsq)
-               nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
+               sys->context = EIG_linear_least_squares_solver_new(0, totvert, 1);
+       else
+               sys->context = EIG_linear_solver_new(0, totvert, 1);
 
        return sys;
 }
@@ -272,7 +270,7 @@ static void laplacian_system_construct_end(LaplacianSystem *sys)
 
                /* for heat weighting */
                if (sys->heat.H)
-                       nlMatrixAdd(sys->context, a, a, sys->heat.H[a]);
+                       EIG_linear_solver_matrix_add(sys->context, a, a, sys->heat.H[a]);
        }
 
        if (sys->storeweights)
@@ -301,7 +299,7 @@ static void laplacian_system_delete(LaplacianSystem *sys)
        if (sys->faces) MEM_freeN(sys->faces);
        if (sys->fweights) MEM_freeN(sys->fweights);
 
-       nlDeleteContext(sys->context);
+       EIG_linear_solver_delete(sys->context);
        MEM_freeN(sys);
 }
 
@@ -309,42 +307,37 @@ void laplacian_begin_solve(LaplacianSystem *sys, int index)
 {
        int a;
 
-       if (!sys->nlbegun) {
-               nlBegin(sys->context, NL_SYSTEM);
-
+       if (!sys->variablesdone) {
                if (index >= 0) {
                        for (a = 0; a < sys->totvert; a++) {
                                if (sys->vpinned[a]) {
-                                       nlSetVariable(sys->context, 0, a, sys->verts[a][index]);
-                                       nlLockVariable(sys->context, a);
+                                       EIG_linear_solver_variable_set(sys->context, 0, a, sys->verts[a][index]);
+                                       EIG_linear_solver_variable_lock(sys->context, a);
                                }
                        }
                }
 
-               nlBegin(sys->context, NL_MATRIX);
-               sys->nlbegun = 1;
+               sys->variablesdone = true;
        }
 }
 
 void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value)
 {
-       nlRightHandSideAdd(sys->context, 0, v, value);
+       EIG_linear_solver_right_hand_side_add(sys->context, 0, v, value);
 }
 
 int laplacian_system_solve(LaplacianSystem *sys)
 {
-       nlEnd(sys->context, NL_MATRIX);
-       nlEnd(sys->context, NL_SYSTEM);
-       sys->nlbegun = 0;
+       sys->variablesdone = false;
 
-       //nlPrintMatrix(sys->context, );
+       //EIG_linear_solver_print_matrix(sys->context, );
 
-       return nlSolve(sys->context, NL_TRUE);
+       return EIG_linear_solver_solve(sys->context);
 }
 
 float laplacian_system_get_solution(LaplacianSystem *sys, int v)
 {
-       return nlGetVariable(sys->context, 0, v);
+       return EIG_linear_solver_variable_get(sys->context, 0, v);
 }
 
 /************************* Heat Bone Weighting ******************************/
@@ -1284,7 +1277,7 @@ static float meshdeform_boundary_total_weight(MeshDeformBind *mdb, int x, int y,
        return totweight;
 }
 
-static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, NLContext *context, int x, int y, int z)
+static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, LinearSolver *context, int x, int y, int z)
 {
        MDefBoundIsect *isect;
        float weight, totweight;
@@ -1294,7 +1287,7 @@ static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, NLContext *context,
        if (mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
                return;
 
-       nlMatrixAdd(context, mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
+       EIG_linear_solver_matrix_add(context, mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
        
        totweight = meshdeform_boundary_total_weight(mdb, x, y, z);
        for (i = 1; i <= 6; i++) {
@@ -1305,12 +1298,12 @@ static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, NLContext *context,
                isect = mdb->boundisect[acenter][i - 1];
                if (!isect) {
                        weight = (1.0f / mdb->width[0]) / totweight;
-                       nlMatrixAdd(context, mdb->varidx[acenter], mdb->varidx[a], -weight);
+                       EIG_linear_solver_matrix_add(context, mdb->varidx[acenter], mdb->varidx[a], -weight);
                }
        }
 }
 
-static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, NLContext *context, int x, int y, int z, int cagevert)
+static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, LinearSolver *context, int x, int y, int z, int cagevert)
 {
        MDefBoundIsect *isect;
        float rhs, weight, totweight;
@@ -1331,7 +1324,7 @@ static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, NLContext *context, i
                if (isect) {
                        weight = (1.0f / isect->len) / totweight;
                        rhs = weight * meshdeform_boundary_phi(mdb, isect, cagevert);
-                       nlRightHandSideAdd(context, 0, mdb->varidx[acenter], rhs);
+                       EIG_linear_solver_right_hand_side_add(context, 0, mdb->varidx[acenter], rhs);
                }
        }
 }
@@ -1386,7 +1379,7 @@ static void meshdeform_matrix_add_exterior_phi(MeshDeformBind *mdb, int x, int y
 
 static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind *mdb)
 {
-       NLContext *context;
+       LinearSolver *context;
        float vec[3], gridvec[3];
        int a, b, x, y, z, totvar;
        char message[256];
@@ -1403,15 +1396,8 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind
 
        progress_bar(0, "Starting mesh deform solve");
 
-       /* setup opennl solver */
-       context = nlNewContext();
-
-       nlSolverParameteri(context, NL_NB_VARIABLES, totvar);
-       nlSolverParameteri(context, NL_NB_ROWS, totvar);
-       nlSolverParameteri(context, NL_NB_RIGHT_HAND_SIDES, 1);
-
-       nlBegin(context, NL_SYSTEM);
-       nlBegin(context, NL_MATRIX);
+       /* setup linear solver */
+       context = EIG_linear_solver_new(totvar, totvar, 1);
 
        /* build matrix */
        for (z = 0; z < mdb->size; z++)
@@ -1421,21 +1407,13 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind
 
        /* solve for each cage vert */
        for (a = 0; a < mdb->totcagevert; a++) {
-               if (a != 0) {
-                       nlBegin(context, NL_SYSTEM);
-                       nlBegin(context, NL_MATRIX);
-               }
-
                /* fill in right hand side and solve */
                for (z = 0; z < mdb->size; z++)
                        for (y = 0; y < mdb->size; y++)
                                for (x = 0; x < mdb->size; x++)
                                        meshdeform_matrix_add_rhs(mdb, context, x, y, z, a);
 
-               nlEnd(context, NL_MATRIX);
-               nlEnd(context, NL_SYSTEM);
-
-               if (nlSolve(context, NL_TRUE)) {
+               if (EIG_linear_solver_solve(context)) {
                        for (z = 0; z < mdb->size; z++)
                                for (y = 0; y < mdb->size; y++)
                                        for (x = 0; x < mdb->size; x++)
@@ -1448,7 +1426,7 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind
 
                        for (b = 0; b < mdb->size3; b++) {
                                if (mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
-                                       mdb->phi[b] = nlGetVariable(context, 0, mdb->varidx[b]);
+                                       mdb->phi[b] = EIG_linear_solver_variable_get(context, 0, mdb->varidx[b]);
                                mdb->totalphi[b] += mdb->phi[b];
                        }
 
@@ -1502,7 +1480,7 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind
        /* free */
        MEM_freeN(mdb->varidx);
 
-       nlDeleteContext(context);
+       EIG_linear_solver_delete(context);
 }
 
 static void harmonic_coordinates_bind(Scene *UNUSED(scene), MeshDeformModifierData *mmd, MeshDeformBind *mdb)
@@ -1705,13 +1683,3 @@ void mesh_deform_bind(Scene *scene, MeshDeformModifierData *mmd, float *vertexco
        waitcursor(0);
 }
 
-#else  /* WITH_OPENNL */
-
-#ifdef __GNUC__
-#  pragma GCC diagnostic ignored "-Wunused-parameter"
-#endif
-
-void mesh_deform_bind(Scene *scene, MeshDeformModifierData *mmd, float *vertexcos, int totvert, float cagemat[4][4]) {}
-void *modifier_mdef_compact_influences_link_kludge = modifier_mdef_compact_influences;
-
-#endif  /* WITH_OPENNL */
index d53e5b857a0718e08abd919b760aa17bde2f97da..a1c1f4f115ec9a92116f8146f9cc04d0ca1e2be0 100644 (file)
@@ -2497,7 +2497,7 @@ int weightFromLoc(EditMesh *em, int axis)
        return 1;
 }
 
-static void addTriangle(NLContext *context, EditVert *v1, EditVert *v2, EditVert *v3, int e1, int e2, int e3)
+static void addTriangle(LinearSolver *context, EditVert *v1, EditVert *v2, EditVert *v3, int e1, int e2, int e3)
 {
        /* Angle opposite e1 */
        float t1 = cotangent_tri_weight_v3(v1->co, v2->co, v3->co) / e2;
@@ -2512,23 +2512,23 @@ static void addTriangle(NLContext *context, EditVert *v1, EditVert *v2, EditVert
        int i2 = indexData(v2);
        int i3 = indexData(v3);
        
-       nlMatrixAdd(context, i1, i1, t2 + t3);
-       nlMatrixAdd(context, i2, i2, t1 + t3);
-       nlMatrixAdd(context, i3, i3, t1 + t2);
+       EIG_linear_solver_matrix_add(context, i1, i1, t2 + t3);
+       EIG_linear_solver_matrix_add(context, i2, i2, t1 + t3);
+       EIG_linear_solver_matrix_add(context, i3, i3, t1 + t2);
 
-       nlMatrixAdd(context, i1, i2, -t3);
-       nlMatrixAdd(context, i2, i1, -t3);
+       EIG_linear_solver_matrix_add(context, i1, i2, -t3);
+       EIG_linear_solver_matrix_add(context, i2, i1, -t3);
 
-       nlMatrixAdd(context, i2, i3, -t1);
-       nlMatrixAdd(context, i3, i2, -t1);
+       EIG_linear_solver_matrix_add(context, i2, i3, -t1);
+       EIG_linear_solver_matrix_add(context, i3, i2, -t1);
 
-       nlMatrixAdd(context, i3, i1, -t2);
-       nlMatrixAdd(context, i1, i3, -t2);
+       EIG_linear_solver_matrix_add(context, i3, i1, -t2);
+       EIG_linear_solver_matrix_add(context, i1, i3, -t2);
 }
 
 int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 {
-       NLContext *context;
+       LinearSolver *context;
        NLboolean success;
        EditVert *eve;
        EditEdge *eed;
@@ -2542,14 +2542,10 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
                totvert++;
        }
 
-       /* Solve with openNL */
+       /* Solve */
        
-       context = nlNewContext();
+       context = EIG_linear_solver_new(, 0, totvert, 1);
 
-       nlSolverParameteri(context, NL_NB_VARIABLES, totvert);
-
-       nlBegin(context, NL_SYSTEM);
-       
        /* Find local extrema */
        for (index = 0, eve = em->verts.first; eve; index++, eve = eve->next) {
                if (eve->h == 0) {
@@ -2583,8 +2579,8 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
                        if (maximum || minimum) {
                                float w = weightData(eve);
                                eve->f1 = 0;
-                               nlSetVariable(context, 0, index, w);
-                               nlLockVariable(context, index);
+                               EIG_linear_solver_variable_set(context, 0, index, w);
+                               EIG_linear_solver_variable_lock(context, index);
                        }
                        else {
                                eve->f1 = 1;
@@ -2592,8 +2588,6 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
                }
        }
        
-       nlBegin(context, NL_MATRIX);
-
        /* Zero edge weight */
        for (eed = em->edges.first; eed; eed = eed->next) {
                eed->tmp.l = 0;
@@ -2625,23 +2619,19 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
                }
        }
        
-       nlEnd(context, NL_MATRIX);
-
-       nlEnd(context, NL_SYSTEM);
-
-       success = nlSolve(context, NL_TRUE);
+       success = EIG_linear_solver_solve(context);
 
        if (success) {
                rval = 1;
                for (index = 0, eve = em->verts.first; eve; index++, eve = eve->next) {
-                       weightSetData(eve, nlGetVariable(context, 0, index));
+                       weightSetData(eve, EIG_linear_solver_variable_get(context, 0, index));
                }
        }
        else {
                rval = 0;
        }
 
-       nlDeleteContext(context);
+       EIG_linear_solver_delete(context);
 
        return rval;
 }
index a90763eed4eb82ea841cc88abfbe4dd2c2115217..543ef0e06632ac34bc6bbc7e6b8b395d614806fe 100644 (file)
@@ -29,6 +29,7 @@ set(INC
        ../../makesrna
        ../../windowmanager
        ../../../../intern/guardedalloc
+       ../../../../intern/eigen
        ../../../../intern/glew-mx
 )
 
@@ -52,13 +53,6 @@ if(WITH_INTERNATIONAL)
        add_definitions(-DWITH_INTERNATIONAL)
 endif()
 
-if(WITH_OPENNL)
-       add_definitions(-DWITH_OPENNL)
-       list(APPEND INC_SYS
-               ../../../../intern/opennl/extern
-       )
-endif()
-
 add_definitions(${GL_DEFINITIONS})
 
 blender_add_lib(bf_editor_uvedit "${SRC}" "${INC}" "${INC_SYS}")
index b5cccab400202ed8457bc4c3af821c012a710f25..85ea7f45536e0f32a00773f28807531ceb1bcb2f 100644 (file)
@@ -34,9 +34,9 @@ sources = env.Glob('*.c')
 
 incs = [
     '#/intern/guardedalloc',
+    '#/intern/eigen',
     env['BF_GLEW_INC'],
     '#/intern/glew-mx',
-    '#/intern/opennl/extern',
     '../include',
     '../../blenkernel',
     '../../blenlib',
index 311c152a239b357d4f039c0e4316381f57d0e8e5..e1495b617f85ab8797152cf157d8f7cfc0c0993f 100644 (file)
@@ -44,9 +44,7 @@
 
 #include "BLI_sys_types.h"  /* for intptr_t support */
 
-#ifdef WITH_OPENNL
-
-#include "ONL_opennl.h"
+#include "eigen_capi.h"
 
 /* Utils */
 
@@ -193,7 +191,7 @@ typedef struct PChart {
 
        union PChartUnion {
                struct PChartLscm {
-                       NLContext *context;
+                       LinearSolver *context;
                        float *abf_alpha;
                        PVert *pin1, *pin2;
                } lscm;
@@ -2471,17 +2469,12 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
        PEdge *e;
        int i, j, ninterior = sys->ninterior, nvar = 2 * sys->ninterior;
        PBool success;
-       NLContext *context;
-
-       context = nlNewContext();
-       nlSolverParameteri(context, NL_NB_VARIABLES, nvar);
+       LinearSolver *context;
 
-       nlBegin(context, NL_SYSTEM);
-
-       nlBegin(context, NL_MATRIX);
+       context = EIG_linear_solver_new(0, nvar, 1);
 
        for (i = 0; i < nvar; i++)
-               nlRightHandSideAdd(context, 0, i, sys->bInterior[i]);
+               EIG_linear_solver_right_hand_side_add(context, 0, i, sys->bInterior[i]);
 
        for (f = chart->faces; f; f = f->nextlink) {
                float wi1, wi2, wi3, b, si, beta[3], j2[3][3], W[3][3];
@@ -2527,8 +2520,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        sys->J2dt[e2->u.id][0] = j2[1][0] = p_abf_compute_sin_product(sys, v1, e2->u.id) * wi2;
                        sys->J2dt[e3->u.id][0] = j2[2][0] = p_abf_compute_sin_product(sys, v1, e3->u.id) * wi3;
 
-                       nlRightHandSideAdd(context, 0, v1->u.id, j2[0][0] * beta[0]);
-                       nlRightHandSideAdd(context, 0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
+                       EIG_linear_solver_right_hand_side_add(context, 0, v1->u.id, j2[0][0] * beta[0]);
+                       EIG_linear_solver_right_hand_side_add(context, 0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
 
                        row1[0] = j2[0][0] * W[0][0];
                        row2[0] = j2[0][0] * W[1][0];
@@ -2547,8 +2540,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        sys->J2dt[e2->u.id][1] = j2[1][1] = 1.0f * wi2;
                        sys->J2dt[e3->u.id][1] = j2[2][1] = p_abf_compute_sin_product(sys, v2, e3->u.id) * wi3;
 
-                       nlRightHandSideAdd(context, 0, v2->u.id, j2[1][1] * beta[1]);
-                       nlRightHandSideAdd(context, 0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
+                       EIG_linear_solver_right_hand_side_add(context, 0, v2->u.id, j2[1][1] * beta[1]);
+                       EIG_linear_solver_right_hand_side_add(context, 0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
 
                        row1[1] = j2[1][1] * W[0][1];
                        row2[1] = j2[1][1] * W[1][1];
@@ -2567,8 +2560,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        sys->J2dt[e2->u.id][2] = j2[1][2] = p_abf_compute_sin_product(sys, v3, e2->u.id) * wi2;
                        sys->J2dt[e3->u.id][2] = j2[2][2] = 1.0f * wi3;
 
-                       nlRightHandSideAdd(context, 0, v3->u.id, j2[2][2] * beta[2]);
-                       nlRightHandSideAdd(context, 0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
+                       EIG_linear_solver_right_hand_side_add(context, 0, v3->u.id, j2[2][2] * beta[2]);
+                       EIG_linear_solver_right_hand_side_add(context, 0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
 
                        row1[2] = j2[2][2] * W[0][2];
                        row2[2] = j2[2][2] * W[1][2];
@@ -2592,29 +2585,25 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                                        continue;
 
                                if (i == 0)
-                                       nlMatrixAdd(context, r, c, j2[0][i] * row1[j]);
+                                       EIG_linear_solver_matrix_add(context, r, c, j2[0][i] * row1[j]);
                                else
-                                       nlMatrixAdd(context, r + ninterior, c, j2[0][i] * row1[j]);
+                                       EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[0][i] * row1[j]);
 
                                if (i == 1)
-                                       nlMatrixAdd(context, r, c, j2[1][i] * row2[j]);
+                                       EIG_linear_solver_matrix_add(context, r, c, j2[1][i] * row2[j]);
                                else
-                                       nlMatrixAdd(context, r + ninterior, c, j2[1][i] * row2[j]);
+                                       EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[1][i] * row2[j]);
 
 
                                if (i == 2)
-                                       nlMatrixAdd(context, r, c, j2[2][i] * row3[j]);
+                                       EIG_linear_solver_matrix_add(context, r, c, j2[2][i] * row3[j]);
                                else
-                                       nlMatrixAdd(context, r + ninterior, c, j2[2][i] * row3[j]);
+                                       EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[2][i] * row3[j]);
                        }
                }
        }
 
-       nlEnd(context, NL_MATRIX);
-
-       nlEnd(context, NL_SYSTEM);
-
-       success = nlSolve(context, NL_FALSE);
+       success = EIG_linear_solver_solve(context);
 
        if (success) {
                for (f = chart->faces; f; f = f->nextlink) {
@@ -2625,24 +2614,24 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        pre[0] = pre[1] = pre[2] = 0.0;
 
                        if (v1->flag & PVERT_INTERIOR) {
-                               float x = nlGetVariable(context, 0, v1->u.id);
-                               float x2 = nlGetVariable(context, 0, ninterior + v1->u.id);
+                               float x = EIG_linear_solver_variable_get(context, 0, v1->u.id);
+                               float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v1->u.id);
                                pre[0] += sys->J2dt[e1->u.id][0] * x;
                                pre[1] += sys->J2dt[e2->u.id][0] * x2;
                                pre[2] += sys->J2dt[e3->u.id][0] * x2;
                        }
 
                        if (v2->flag & PVERT_INTERIOR) {
-                               float x = nlGetVariable(context, 0, v2->u.id);
-                               float x2 = nlGetVariable(context, 0, ninterior + v2->u.id);
+                               float x = EIG_linear_solver_variable_get(context, 0, v2->u.id);
+                               float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v2->u.id);
                                pre[0] += sys->J2dt[e1->u.id][1] * x2;
                                pre[1] += sys->J2dt[e2->u.id][1] * x;
                                pre[2] += sys->J2dt[e3->u.id][1] * x2;
                        }
 
                        if (v3->flag & PVERT_INTERIOR) {
-                               float x = nlGetVariable(context, 0, v3->u.id);
-                               float x2 = nlGetVariable(context, 0, ninterior + v3->u.id);
+                               float x = EIG_linear_solver_variable_get(context, 0, v3->u.id);
+                               float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v3->u.id);
                                pre[0] += sys->J2dt[e1->u.id][2] * x2;
                                pre[1] += sys->J2dt[e2->u.id][2] * x2;
                                pre[2] += sys->J2dt[e3->u.id][2] * x;
@@ -2673,12 +2662,12 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                }
 
                for (i = 0; i < ninterior; i++) {
-                       sys->lambdaPlanar[i] += (float)nlGetVariable(context, 0, i);
-                       sys->lambdaLength[i] += (float)nlGetVariable(context, 0, ninterior + i);
+                       sys->lambdaPlanar[i] += (float)EIG_linear_solver_variable_get(context, 0, i);
+                       sys->lambdaLength[i] += (float)EIG_linear_solver_variable_get(context, 0, ninterior + i);
                }
        }
 
-       nlDeleteContext(context);
+       EIG_linear_solver_delete(context);
 
        return success;
 }
@@ -3004,12 +2993,12 @@ static void p_chart_extrema_verts(PChart *chart, PVert **pin1, PVert **pin2)
 
 static void p_chart_lscm_load_solution(PChart *chart)
 {
-       NLContext *context = chart->u.lscm.context;
+       LinearSolver *context = chart->u.lscm.context;
        PVert *v;
 
        for (v = chart->verts; v; v = v->nextlink) {
-               v->uv[0] = nlGetVariable(context, 0, 2 * v->u.id);
-               v->uv[1] = nlGetVariable(context, 0, 2 * v->u.id + 1);
+               v->uv[0] = EIG_linear_solver_variable_get(context, 0, 2 * v->u.id);
+               v->uv[1] = EIG_linear_solver_variable_get(context, 0, 2 * v->u.id + 1);
        }
 }
 
@@ -3064,16 +3053,13 @@ static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
                for (v = chart->verts; v; v = v->nextlink)
                        v->u.id = id++;
 
-               chart->u.lscm.context = nlNewContext();
-               nlSolverParameteri(chart->u.lscm.context, NL_NB_VARIABLES, 2 * chart->nverts);
-               nlSolverParameteri(chart->u.lscm.context, NL_NB_ROWS, 2 * chart->nfaces);
-               nlSolverParameteri(chart->u.lscm.context, NL_LEAST_SQUARES, NL_TRUE);
+               chart->u.lscm.context = EIG_linear_least_squares_solver_new(2 * chart->nfaces, 2 * chart->nverts, 1);
        }
 }
 
 static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 {
-       NLContext *context = chart->u.lscm.context;
+       LinearSolver *context = chart->u.lscm.context;
        PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
        PFace *f;
        float *alpha = chart->u.lscm.abf_alpha;
@@ -3081,8 +3067,6 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
        bool flip_faces;
        int row;
 
-       nlBegin(context, NL_SYSTEM);
-
 #if 0
        /* TODO: make loading pins work for simplify/complexify. */
 #endif
@@ -3092,25 +3076,25 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
                        p_vert_load_pin_select_uvs(handle, v);  /* reload for live */
 
        if (chart->u.lscm.pin1) {
-               nlLockVariable(context, 2 * pin1->u.id);
-               nlLockVariable(context, 2 * pin1->u.id + 1);
-               nlLockVariable(context, 2 * pin2->u.id);
-               nlLockVariable(context, 2 * pin2->u.id + 1);
+               EIG_linear_solver_variable_lock(context, 2 * pin1->u.id);
+               EIG_linear_solver_variable_lock(context, 2 * pin1->u.id + 1);
+               EIG_linear_solver_variable_lock(context, 2 * pin2->u.id);
+               EIG_linear_solver_variable_lock(context, 2 * pin2->u.id + 1);
 
-               nlSetVariable(context, 0, 2 * pin1->u.id, pin1->uv[0]);
-               nlSetVariable(context, 0, 2 * pin1->u.id + 1, pin1->uv[1]);
-               nlSetVariable(context, 0, 2 * pin2->u.id, pin2->uv[0]);
-               nlSetVariable(context, 0, 2 * pin2->u.id + 1, pin2->uv[1]);
+               EIG_linear_solver_variable_set(context, 0, 2 * pin1->u.id, pin1->uv[0]);
+               EIG_linear_solver_variable_set(context, 0, 2 * pin1->u.id + 1, pin1->uv[1]);
+               EIG_linear_solver_variable_set(context, 0, 2 * pin2->u.id, pin2->uv[0]);
+               EIG_linear_solver_variable_set(context, 0, 2 * pin2->u.id + 1, pin2->uv[1]);
        }
        else {
                /* set and lock the pins */
                for (v = chart->verts; v; v = v->nextlink) {
                        if (v->flag & PVERT_PIN) {
-                               nlLockVariable(context, 2 * v->u.id);
-                               nlLockVariable(context, 2 * v->u.id + 1);
+                               EIG_linear_solver_variable_lock(context, 2 * v->u.id);
+                               EIG_linear_solver_variable_lock(context, 2 * v->u.id + 1);
 
-                               nlSetVariable(context, 0, 2 * v->u.id, v->uv[0]);
-                               nlSetVariable(context, 0, 2 * v->u.id + 1, v->uv[1]);
+                               EIG_linear_solver_variable_set(context, 0, 2 * v->u.id, v->uv[0]);
+                               EIG_linear_solver_variable_set(context, 0, 2 * v->u.id + 1, v->uv[1]);
                        }
                }
        }
@@ -3137,8 +3121,6 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 
        /* construct matrix */
 
-       nlBegin(context, NL_MATRIX);
-
        row = 0;
        for (f = chart->faces; f; f = f->nextlink) {
                PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
@@ -3185,26 +3167,22 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
                cosine = cosf(a1) * ratio;
                sine = sina1 * ratio;
 
-               nlMatrixAdd(context, row, 2 * v1->u.id,   cosine - 1.0f);
-               nlMatrixAdd(context, row, 2 * v1->u.id + 1, -sine);
-               nlMatrixAdd(context, row, 2 * v2->u.id,   -cosine);
-               nlMatrixAdd(context, row, 2 * v2->u.id + 1, sine);
-               nlMatrixAdd(context, row, 2 * v3->u.id,   1.0);
+               EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id,   cosine - 1.0f);
+               EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id + 1, -sine);
+               EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id,   -cosine);
+               EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id + 1, sine);
+               EIG_linear_solver_matrix_add(context, row, 2 * v3->u.id,   1.0);
                row++;
 
-               nlMatrixAdd(context, row, 2 * v1->u.id,   sine);
-               nlMatrixAdd(context, row, 2 * v1->u.id + 1, cosine - 1.0f);
-               nlMatrixAdd(context, row, 2 * v2->u.id,   -sine);
-               nlMatrixAdd(context, row, 2 * v2->u.id + 1, -cosine);
-               nlMatrixAdd(context, row, 2 * v3->u.id + 1, 1.0);
+               EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id,   sine);
+               EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id + 1, cosine - 1.0f);
+               EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id,   -sine);
+               EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id + 1, -cosine);
+               EIG_linear_solver_matrix_add(context, row, 2 * v3->u.id + 1, 1.0);
                row++;
        }
 
-       nlEnd(context, NL_MATRIX);
-
-       nlEnd(context, NL_SYSTEM);
-
-       if (nlSolve(context, NL_TRUE)) {
+       if (EIG_linear_solver_solve(context)) {
                p_chart_lscm_load_solution(chart);
                return P_TRUE;
        }
@@ -3221,7 +3199,7 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 static void p_chart_lscm_end(PChart *chart)
 {
        if (chart->u.lscm.context)
-               nlDeleteContext(chart->u.lscm.context);
+               EIG_linear_solver_delete(chart->u.lscm.context);
        
        if (chart->u.lscm.abf_alpha) {
                MEM_freeN(chart->u.lscm.abf_alpha);
@@ -4700,36 +4678,3 @@ void param_flush_restore(ParamHandle *handle)
        }
 }
 
-#else  /* WITH_OPENNL */
-
-#ifdef __GNUC__
-#  pragma GCC diagnostic ignored "-Wunused-parameter"
-#endif
-
-/* stubs */
-void param_face_add(ParamHandle *handle, ParamKey key, int nverts,
-                    ParamKey *vkeys, float **co, float **uv,
-                    ParamBool *pin, ParamBool *select, float normal[3]) {}
-void param_edge_set_seam(ParamHandle *handle,
-                         ParamKey *vkeys) {}
-void param_aspect_ratio(ParamHandle *handle, float aspx, float aspy) {}
-ParamHandle *param_construct_begin(void) { return NULL; }
-void param_construct_end(ParamHandle *handle, ParamBool fill, ParamBool impl) {}
-void param_delete(ParamHandle *handle) {}
-
-void param_stretch_begin(ParamHandle *handle) {}
-void param_stretch_blend(ParamHandle *handle, float blend) {}
-void param_stretch_iter(ParamHandle *handle) {}
-void param_stretch_end(ParamHandle *handle) {}
-
-void param_pack(ParamHandle *handle, float margin, bool do_rotate) {}
-void param_average(ParamHandle *handle) {}
-
-void param_flush(ParamHandle *handle) {}
-void param_flush_restore(ParamHandle *handle) {}
-
-void param_lscm_begin(ParamHandle *handle, ParamBool live, ParamBool abf) {}
-void param_lscm_solve(ParamHandle *handle) {}
-void param_lscm_end(ParamHandle *handle) {}
-
-#endif  /* WITH_OPENNL */
index ad230dede2402fe00d0e1090de2d9dc5af2b5def..0de7676e8f84eb571b0c0122c676df55b52f765e 100644 (file)
@@ -37,6 +37,7 @@ set(INC
        ../render/extern/include
        ../../../intern/elbeem/extern
        ../../../intern/guardedalloc
+       ../../../intern/eigen
 )
 
 set(INC_SYS
@@ -144,13 +145,6 @@ if(WITH_INTERNATIONAL)
        add_definitions(-DWITH_INTERNATIONAL)
 endif()
 
-if(WITH_OPENNL)
-       add_definitions(-DWITH_OPENNL)
-       list(APPEND INC_SYS
-               ../../../intern/opennl/extern
-       )
-endif()
-
 if(WITH_OPENSUBDIV)
        add_definitions(-DWITH_OPENSUBDIV)
 endif()
index 7be295aa1a02c1e77a4db9efefa57ca816660dab..300071185622e98fb10e24693776531186149f2b 100644 (file)
@@ -33,8 +33,8 @@ incs = [
     '.',
     './intern', 
     '#/intern/guardedalloc',
+    '#/intern/eigen',
     '#/intern/elbeem/extern',
-    '#/intern/opennl/extern',
     '../render/extern/include',
     '../bmesh',
     '../include',
index fdaacc7cd9e09cff5c98c11394242018faa71300..d4f02d923d31a8530b9803af8ac55181160cdb35 100644 (file)
@@ -42,6 +42,7 @@
 
 #include "MOD_util.h"
 
+#include "eigen_capi.h"
 
 enum {
        LAPDEFORM_SYSTEM_NOT_CHANGE = 0,
@@ -54,10 +55,6 @@ enum {
        LAPDEFORM_SYSTEM_CHANGE_NOT_VALID_GROUP,
 };
 
-#ifdef WITH_OPENNL
-
-#include "ONL_opennl.h"
-
 typedef struct LaplacianSystem {
        bool is_matrix_computed;
        bool has_solution;
@@ -75,7 +72,7 @@ typedef struct LaplacianSystem {
        int *unit_verts;                        /* Unit vectors of projected edges onto the plane orthogonal to n */
        int *ringf_indices;                     /* Indices of faces per vertex */
        int *ringv_indices;                     /* Indices of neighbors(vertex) per vertex */
-       NLContext *context;                     /* System for solve general implicit rotations */
+       LinearSolver *context;                  /* System for solve general implicit rotations */
        MeshElemMap *ringf_map;         /* Map of faces per vertex */
        MeshElemMap *ringv_map;         /* Map of vertex per vertex */
 } LaplacianSystem;
@@ -134,7 +131,7 @@ static void deleteLaplacianSystem(LaplacianSystem *sys)
        MEM_SAFE_FREE(sys->ringv_map);
 
        if (sys->context) {
-               nlDeleteContext(sys->context);
+               EIG_linear_solver_delete(sys->context);
        }
        MEM_SAFE_FREE(sys);
 }
@@ -283,9 +280,9 @@ static void initLaplacianMatrix(LaplacianSystem *sys)
                        sys->delta[idv[0]][1] -= v3[1] * w3;
                        sys->delta[idv[0]][2] -= v3[2] * w3;
 
-                       nlMatrixAdd(sys->context, idv[0], idv[1], -w2);
-                       nlMatrixAdd(sys->context, idv[0], idv[2], -w3);
-                       nlMatrixAdd(sys->context, idv[0], idv[0], w2 + w3);
+                       EIG_linear_solver_matrix_add(sys->context, idv[0], idv[1], -w2);
+                       EIG_linear_solver_matrix_add(sys->context, idv[0], idv[2], -w3);
+                       EIG_linear_solver_matrix_add(sys->context, idv[0], idv[0], w2 + w3);
                }
        }
 }
@@ -338,9 +335,9 @@ static void rotateDifferentialCoordinates(LaplacianSystem *sys)
                beta = dot_v3v3(uij, di);
                gamma = dot_v3v3(e2, di);
 
-               pi[0] = nlGetVariable(sys->context, 0, i);
-               pi[1] = nlGetVariable(sys->context, 1, i);
-               pi[2] = nlGetVariable(sys->context, 2, i);
+               pi[0] = EIG_linear_solver_variable_get(sys->context, 0, i);
+               pi[1] = EIG_linear_solver_variable_get(sys->context, 1, i);
+               pi[2] = EIG_linear_solver_variable_get(sys->context, 2, i);
                zero_v3(ni);
                num_fni = 0;
                num_fni = sys->ringf_map[i].count;
@@ -349,9 +346,9 @@ static void rotateDifferentialCoordinates(LaplacianSystem *sys)
                        fidn = sys->ringf_map[i].indices;
                        vin = sys->tris[fidn[fi]];
                        for (j = 0; j < 3; j++) {
-                               vn[j][0] = nlGetVariable(sys->context, 0, vin[j]);
-                               vn[j][1] = nlGetVariable(sys->context, 1, vin[j]);
-                               vn[j][2] = nlGetVariable(sys->context, 2, vin[j]);
+                               vn[j][0] = EIG_linear_solver_variable_get(sys->context, 0, vin[j]);
+                               vn[j][1] = EIG_linear_solver_variable_get(sys->context, 1, vin[j]);
+                               vn[j][2] = EIG_linear_solver_variable_get(sys->context, 2, vin[j]);
                                if (vin[j] == sys->unit_verts[i]) {
                                        copy_v3_v3(pj, vn[j]);
                                }
@@ -372,14 +369,14 @@ static void rotateDifferentialCoordinates(LaplacianSystem *sys)
                fni[2] = alpha * ni[2] + beta * uij[2] + gamma * e2[2];
 
                if (len_squared_v3(fni) > FLT_EPSILON) {
-                       nlRightHandSideSet(sys->context, 0, i, fni[0]);
-                       nlRightHandSideSet(sys->context, 1, i, fni[1]);
-                       nlRightHandSideSet(sys->context, 2, i, fni[2]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 0, i, fni[0]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 1, i, fni[1]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 2, i, fni[2]);
                }
                else {
-                       nlRightHandSideSet(sys->context, 0, i, sys->delta[i][0]);
-                       nlRightHandSideSet(sys->context, 1, i, sys->delta[i][1]);
-                       nlRightHandSideSet(sys->context, 2, i, sys->delta[i][2]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 0, i, sys->delta[i][0]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 1, i, sys->delta[i][1]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 2, i, sys->delta[i][2]);
                }
        }
 }
@@ -390,75 +387,59 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3])
        n = sys->total_verts;
        na = sys->total_anchors;
 
-#ifdef OPENNL_THREADING_HACK
-       modifier_opennl_lock();
-#endif
-
        if (!sys->is_matrix_computed) {
-               sys->context = nlNewContext();
+               sys->context = EIG_linear_least_squares_solver_new(n + na, n, 3);
 
-               nlSolverParameteri(sys->context, NL_NB_VARIABLES, n);
-               nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
-               nlSolverParameteri(sys->context, NL_NB_ROWS, n + na);
-               nlSolverParameteri(sys->context, NL_NB_RIGHT_HAND_SIDES, 3);
-               nlBegin(sys->context, NL_SYSTEM);
                for (i = 0; i < n; i++) {
-                       nlSetVariable(sys->context, 0, i, sys->co[i][0]);
-                       nlSetVariable(sys->context, 1, i, sys->co[i][1]);
-                       nlSetVariable(sys->context, 2, i, sys->co[i][2]);
+                       EIG_linear_solver_variable_set(sys->context, 0, i, sys->co[i][0]);
+                       EIG_linear_solver_variable_set(sys->context, 1, i, sys->co[i][1]);
+                       EIG_linear_solver_variable_set(sys->context, 2, i, sys->co[i][2]);
                }
                for (i = 0; i < na; i++) {
                        vid = sys->index_anchors[i];
-                       nlSetVariable(sys->context, 0, vid, vertexCos[vid][0]);
-                       nlSetVariable(sys->context, 1, vid, vertexCos[vid][1]);
-                       nlSetVariable(sys->context, 2, vid, vertexCos[vid][2]);
+                       EIG_linear_solver_variable_set(sys->context, 0, vid, vertexCos[vid][0]);
+                       EIG_linear_solver_variable_set(sys->context, 1, vid, vertexCos[vid][1]);
+                       EIG_linear_solver_variable_set(sys->context, 2, vid, vertexCos[vid][2]);
                }
-               nlBegin(sys->context, NL_MATRIX);
 
                initLaplacianMatrix(sys);
                computeImplictRotations(sys);
 
                for (i = 0; i < n; i++) {
-                       nlRightHandSideSet(sys->context, 0, i, sys->delta[i][0]);
-                       nlRightHandSideSet(sys->context, 1, i, sys->delta[i][1]);
-                       nlRightHandSideSet(sys->context, 2, i, sys->delta[i][2]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 0, i, sys->delta[i][0]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 1, i, sys->delta[i][1]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 2, i, sys->delta[i][2]);
                }
                for (i = 0; i < na; i++) {
                        vid = sys->index_anchors[i];
-                       nlRightHandSideSet(sys->context, 0, n + i, vertexCos[vid][0]);
-                       nlRightHandSideSet(sys->context, 1, n + i, vertexCos[vid][1]);
-                       nlRightHandSideSet(sys->context, 2, n + i, vertexCos[vid][2]);
-                       nlMatrixAdd(sys->context, n + i, vid, 1.0f);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 0, n + i, vertexCos[vid][0]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 1, n + i, vertexCos[vid][1]);
+                       EIG_linear_solver_right_hand_side_add(sys->context, 2, n + i, vertexCos[vid][2]);
+                       EIG_linear_solver_matrix_add(sys->context, n + i, vid, 1.0f);
                }
-               nlEnd(sys->context, NL_MATRIX);
-               nlEnd(sys->context, NL_SYSTEM);
-               if (nlSolve(sys->context, NL_TRUE)) {
+               if (EIG_linear_solver_solve(sys->context)) {
                        sys->has_solution = true;
 
                        for (j = 1; j <= sys->repeat; j++) {
-                               nlBegin(sys->context, NL_SYSTEM);
-                               nlBegin(sys->context, NL_MATRIX);
                                rotateDifferentialCoordinates(sys);
 
       &n