===== scipy-1.18.1/LICENSES_bundled.txt ===== ---- The SciPy repository and source distributions bundle a number of libraries that are compatibly licensed. We list these here. Name: fast_matrix_market Files: scipy/io/_fast_matrix_market/* License: BSD-2-Clause For details, see scipy/io/_fast_matrix_market/LICENSE.txt Name: pystreambuf Files: scipy/io/_fast_matrix_market/src/pystreambuf.h License: BSD-3-Clause For details, see the header inside scipy/io/_fast_matrix_market/src/pystreambuf.h Name: fast_float Files: scipy/io/_fast_matrix_market/fast_matrix_market/dependencies/fast_float/* License: MIT For details, see scipy/io/_fast_matrix_market/fast_matrix_market/dependencies/fast_float/LICENSE-MIT Name: ryu Files: scipy/io/_fast_matrix_market/fast_matrix_market/dependencies/ryu/* License: BSL-1.0 For details, see scipy/io/_fast_matrix_market/fast_matrix_market/dependencies/ryu/LICENSE-Boost Name: L-BFGS-B Files: scipy/optimize/lbfgsb/[__lbfgsb.h,__lbfgsb.c] License: BSD-3-Clause For details, see scipy/optimize/__lbfgsb.c Name: LAPJVsp Files: scipy/sparse/csgraph/_matching.pyx License: BSD-3-Clause Copyright 1987-, A. Volgenant/Amsterdam School of Economics, University of Amsterdam Distributed under 3-clause BSD license with permission from University of Amsterdam. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. Name: SuperLU Files: scipy/sparse/linalg/_dsolve/SuperLU/* License: BSD-3-Clause For details, see scipy/sparse/linalg/_dsolve/SuperLU/License.txt Name: ARPACK Files: scipy/sparse/linalg/_eigen/arpack/ARPACK/* License: BSD-3-Clause For details, see scipy/sparse/linalg/_eigen/arpack/ARPACK/COPYING Name: Qhull Files: subprojects/qhull_r/libqhull_r/* License: Qhull For details, see subprojects/qhull_r/libqhull_r/COPYING.txt Name: xsf Files: subprojects/xsf/* License: BSD-3-Clause AND MIT AND BSD-3-Clause-LBNL AND Apache-2.0 WITH LLVM-exception For details, see subprojects/xsf/LICENSE and subprojects/xsf/LICENSES_bundled.txt Name: pyduccfft Files: scipy/fft/_duccfft/pyduccfft.cxx License: 3-Clause BSD For details, see scipy/fft/_duccfft/LICENSE.md Name: ducc0 Files: subprojects/duccfft/ducc0/* License: 3-Clause BSD For details, see subprojects/duccfft/LICENSE.md Name: uarray Files: scipy/_lib/_uarray/* License: BSD-3-Clause For details, see scipy/_lib/_uarray/LICENSE Name: ampgo Files: benchmarks/benchmarks/go_benchmark_functions/*.py License: MIT Functions for testing global optimizers, forked from the AMPGO project, https://code.google.com/archive/p/ampgo Name: pybind11 Files: no source files are included, however pybind11 binary artifacts are included with every binary build of SciPy. License: BSD-3-Clause Copyright (c) 2016 Wenzel Jakob , All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. Name: HiGHS Files: subprojects/highs/* License: MIT For details, see subprojects/highs/LICENSE.txt Name: Boost Files: subprojects/boost_math/math/* License: BSL-1.0 For details, see subprojects/boost_math/math/LICENSE Name: Biasedurn Files: scipy/stats/biasedurn/* License: BSD-3-Clause For details, see scipy/stats/biasedurn/license.txt Name: UNU.RAN Files: subprojects/unuran/unuran/* License: BSD-3-Clause For details, see subprojects/unuran/unuran/license.txt Name: NumPy Files: scipy/stats/libnpyrandom/[logfactorial.h,logfactorial.c] License: BSD-3-Clause For details, see header inside scipy/stats/libnpyrandom/logfactorial.h and scipy/stats/libnpyrandom/logfactorial.c Name: array-api-compat Files: subprojects/array_api_compat/array_api_compat/* License: MIT For details, see subprojects/array_api_compat/array_api_compat/LICENSE Name: array-api-extra Files: subprojects/array_api_extra/* License: MIT For details, see subprojects/array_api_extra/LICENSE Name: Tempita Files: scipy/_build_utils/tempita/* License: MIT For details, see scipy/_build_utils/tempita/LICENSE.txt Name: Chebfun Files: scipy/interpolate/[_bary_rational.py, tests/test_bary_rational.py] License: BSD-3-Clause For details, see scipy/interpolate/_bary_rational.py Name: getLebedevSphere Files: scipy/integrate/_lebedev.py License: BSD-2-Clause For details, see scipy/integrate/_lebedev.py Name: prima Files: subprojects/pyprima/pyprima License: BSD-3-Clause For details, see subprojects/pyprima/pyprima/LICENCE.txt Name: packaging.version Files: scipy/_external/packaging_version/* License: BSD-3-Clause For details, see scipy/_external/src/LICENSE.BSD ===== scipy-1.18.1/scipy/_build_utils/tempita/LICENSE.txt ===== Copyright (c) 2008 Ian Bicking and Contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/scipy/_external/packaging_version/src/LICENSE.BSD ===== Copyright (c) Donald Stufft and individual contributors. All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/_lib/_uarray/LICENSE ===== BSD 3-Clause License Copyright (c) 2018, Quansight-Labs All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. * Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/fft/_duccfft/LICENSE.md ===== Copyright (C) 2010-2019 Max-Planck-Society All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. * Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/io/_fast_matrix_market/LICENSE.txt ===== Copyright (C) 2020-2023 Adam Lugowski Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/io/_fast_matrix_market/fast_matrix_market/dependencies/fast_float/LICENSE-MIT ===== MIT License Copyright (c) 2021 The fast_float authors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/scipy/io/_fast_matrix_market/fast_matrix_market/dependencies/ryu/LICENSE-Boost ===== Boost Software License - Version 1.0 - August 17th, 2003 Permission is hereby granted, free of charge, to any person or organization obtaining a copy of the software and accompanying documentation covered by this license (the "Software") to use, reproduce, display, distribute, execute, and transmit the Software, and to prepare derivative works of the Software, and to permit third-parties to whom the Software is furnished to do so, all subject to the following: The copyright notices in the Software and this entire statement, including the above license grant, this restriction and the following disclaimer, must be included in all copies of the Software, in whole or in part, and all derivative works of the Software, unless such copies or derivative works are solely in the form of machine-executable object code generated by a source language processor. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/scipy/ndimage/LICENSE.txt ===== Copyright (C) 2003-2005 Peter J. Verveer Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. The name of the author may not be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/optimize/_direct/COPYING ===== This code is based on the DIRECT 2.0.4 Fortran code by Gablonsky et al. at http://www4.ncsu.edu/~ctk/SOFTWARE/DIRECTv204.tar.gz The C version was initially converted via f2c and then cleaned up and reorganized by Steven G. Johnson (stevenj@alum.mit.edu), August 2007. ******** Copyright and license for the original Fortran DIRECT code ******** Copyright (c) 1999, 2000, 2001 North Carolina State University This program is distributed under the MIT License (see http://www.opensource.org/licenses/mit-license.php): Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/scipy/optimize/tnc/LICENSE ===== Copyright (c) 2002-2005, Jean-Sebastien Roy (js@jeannot.org) Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/scipy/sparse/linalg/_dsolve/SuperLU/License.txt ===== Copyright (c) 2003, The Regents of the University of California, through Lawrence Berkeley National Laboratory (subject to receipt of any required approvals from U.S. Dept. of Energy) All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: (1) Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. (2) Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. (3) Neither the name of Lawrence Berkeley National Laboratory, U.S. Dept. of Energy nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/sparse/linalg/_eigen/arpack/arnaud/ARPACK_LICENSE.txt ===== This is the ARPACK package from http://www.caam.rice.edu/software/ARPACK/ Specifically the files are from http://www.caam.rice.edu/software/ARPACK/SRC/arpack96.tar.gz with the patch http://www.caam.rice.edu/software/ARPACK/SRC/patch.tar.gz The ARPACK README is at http://www.caam.rice.edu/software/ARPACK/SRC/readme.arpack --- ARPACK is a collection of Fortran77 subroutines designed to solve large scale eigenvalue problems. The package is designed to compute a few eigenvalues and corresponding eigenvectors of a general n by n matrix A. It is most appropriate for large sparse or structured matrices A where structured means that a matrix-vector product w <- Av requires order n rather than the usual order n**2 floating point operations. This software is based upon an algorithmic variant of the Arnoldi process called the Implicitly Restarted Arnoldi Method (IRAM). When the matrix A is symmetric it reduces to a variant of the Lanczos process called the Implicitly Restarted Lanczos Method (IRLM). These variants may be viewed as a synthesis of the Arnoldi/Lanczos process with the Implicitly Shifted QR technique that is suitable for large scale problems. For many standard problems, a matrix factorization is not required. Only the action of the matrix on a vector is needed. ARPACK software is capable of solving large scale symmetric, nonsymmetric, and generalized eigenproblems from significant application areas. The software is designed to compute a few (k) eigenvalues with user specified features such as those of largest real part or largest magnitude. Storage requirements are on the order of n*k locations. No auxiliary storage is required. A set of Schur basis vectors for the desired k-dimensional eigen-space is computed which is numerically orthogonal to working precision. Numerically accurate eigenvectors are available on request. Important Features: o Reverse Communication Interface. o Single and Double Precision Real Arithmetic Versions for Symmetric, Non-symmetric, Standard or Generalized Problems. o Single and Double Precision Complex Arithmetic Versions for Standard or Generalized Problems. o Routines for Banded Matrices - Standard or Generalized Problems. o Routines for The Singular Value Decomposition. o Example driver routines that may be used as templates to implement numerous Shift-Invert strategies for all problem types, data types and precision. --- The ARPACK license is the BSD 3-clause license ("New BSD License") http://www.caam.rice.edu/software/ARPACK/RiceBSD.txt --- BSD Software License Pertains to ARPACK and P_ARPACK Copyright (c) 1996-2008 Rice University. Developed by D.C. Sorensen, R.B. Lehoucq, C. Yang, and K. Maschhoff. All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: - Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. - Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer listed in this license in the documentation and/or other materials provided with the distribution. - Neither the name of the copyright holders nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/sparse/linalg/_eigen/arpack/arnaud/LICENSE ===== Author: Ilhan Polat Copyright (C) 2025 SciPy developers Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: a. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. b. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. c. Names of the SciPy Developers may not be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/sparse/linalg/_propack/PROPACK/LICENSE ===== Author: Ilhan Polat Copyright (C) 2025 SciPy developers Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: a. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. b. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. c. Names of the SciPy Developers may not be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/sparse/linalg/_propack/PROPACK/PROPACK_LICENSE.txt ===== Reproduced from http://soi.stanford.edu/~rmunk/PROPACK/license.txt BSD License for PROPACK numerical software http://www.opensource.org/licenses/bsd-license.php This notice applies to the software packages PROPACK77 (Fortran 77) and PROPACK (Matlab) made available at http://sun.stanford.edu/~rmunk/PROPACK by Solar Oscillations Investigation, W.W. Hansen Experimental Physics Laboratory (HEPL), 452 Lomita Mall, Stanford University, Stanford, CA 94305-4085, USA Copyright (c) 2005, Rasmus Munk Larsen, Stanford University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. * Neither the name of Stanford University nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/spatial/COPYING_QHULL.txt ===== Qhull, Copyright (c) 1993-2020 C.B. Barber Arlington, MA and The National Science and Technology Research Center for Computation and Visualization of Geometric Structures (The Geometry Center) University of Minnesota email: qhull@qhull.org This software includes Qhull from C.B. Barber and The Geometry Center. Files derived from Qhull 1.0 are copyrighted by the Geometry Center. The remaining files are copyrighted by C.B. Barber. Qhull is free software and may be obtained via http from www.qhull.org. It may be freely copied, modified, and redistributed under the following conditions: 1. All copyright notices must remain intact in all files. 2. A copy of this text file must be distributed along with any copies of Qhull that you redistribute; this includes copies that you have modified, or copies of programs or other software products that include Qhull. 3. If you modify Qhull, you must include a notice giving the name of the person performing the modification, the date of modification, and the reason for such modification. 4. When distributing modified versions of Qhull, or other software products that include Qhull, you must provide notice that the original source code may be obtained as noted above. 5. There is no warranty or other guarantee of fitness for Qhull, it is provided solely "as is". Bug reports or fixes may be sent to qhull_bug@qhull.org; the authors may or may not act on them as they desire. ===== scipy-1.18.1/scipy/stats/biasedurn/license.txt ===== ---------- Forwarded message --------- From: Agner Fog Date: Wed, Jul 22, 2020 at 8:54 PM Subject: Re: Add BiasedUrn to SciPy? To: Matt Haberland Cc: Warren Weckesser Dear Matt The C++ code for these distributions without binding to the R project are available here: https://www.agner.org/random/stocc.zip You have my permission to use it in SciPy under the new BSD license. Best regards, Agner Fog On 22/07/2020 17.05, Matt Haberland wrote: Hi Agner, I'm a maintainer of the Python scientific computing library SciPy, and I'm working on improving the capabilities of our statistics subpackage  scipy.stats. We'd like to add both Fisher's and Wallenius' noncentral hypergeometric distributions, but these look quite challenging to implement. We found your BiasedUrn package, and we'd love to wrap it, but we can't include GPL-3 licensed code under our BSD-new license. Would you consider re-licensing BiasedUrn under a more permissive license, such as BSD-new or MIT, so we can include it in SciPy? Thank you for considering it! Matt -- Matt Haberland Assistant Professor BioResource and Agricultural Engineering 08A-3K, Cal Poly ===== scipy-1.18.1/subprojects/array_api_compat/array_api_compat/LICENSE ===== MIT License Copyright (c) 2022 Consortium for Python Data API Standards Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/subprojects/array_api_extra/LICENSE ===== MIT License Copyright (c) 2024 Consortium for Python Data API Standards Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/subprojects/boost_math/math/LICENSE ===== Boost Software License - Version 1.0 - August 17th, 2003 Permission is hereby granted, free of charge, to any person or organization obtaining a copy of the software and accompanying documentation covered by this license (the "Software") to use, reproduce, display, distribute, execute, and transmit the Software, and to prepare derivative works of the Software, and to permit third-parties to whom the Software is furnished to do so, all subject to the following: The copyright notices in the Software and this entire statement, including the above license grant, this restriction and the following disclaimer, must be included in all copies of the Software, in whole or in part, and all derivative works of the Software, unless such copies or derivative works are solely in the form of machine-executable object code generated by a source language processor. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/subprojects/cobyqa/LICENSE ===== BSD 3-Clause License Copyright (c) 2021-2026, Tom M. Ragonneau and Zaikun Zhang Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/subprojects/duccfft/LICENSE.md ===== Copyright (C) 2010-2019 Max-Planck-Society All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. * Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/subprojects/highs/LICENSE.txt ===== MIT License Copyright (c) 2024 HiGHS Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/subprojects/highs/extern/LICENCE_1_0.txt ===== Boost Software License - Version 1.0 - August 17th, 2003 Permission is hereby granted, free of charge, to any person or organization obtaining a copy of the software and accompanying documentation covered by this license (the "Software") to use, reproduce, display, distribute, execute, and transmit the Software, and to prepare derivative works of the Software, and to permit third-parties to whom the Software is furnished to do so, all subject to the following: The copyright notices in the Software and this entire statement, including the above license grant, this restriction and the following disclaimer, must be included in all copies of the Software, in whole or in part, and all derivative works of the Software, unless such copies or derivative works are solely in the form of machine-executable object code generated by a source language processor. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/subprojects/highs/extern/filereaderlp/LICENSE ===== Copyright (c) 2020 Michael Feldmeier Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ===== scipy-1.18.1/subprojects/highs/extern/pdqsort/license.txt ===== Copyright (c) 2021 Orson Peters This software is provided 'as-is', without any express or implied warranty. In no event will the authors be held liable for any damages arising from the use of this software. Permission is granted to anyone to use this software for any purpose, including commercial applications, and to alter it and redistribute it freely, subject to the following restrictions: 1. The origin of this software must not be misrepresented; you must not claim that you wrote the original software. If you use this software in a product, an acknowledgment in the product documentation would be appreciated but is not required. 2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software. 3. This notice may not be removed or altered from any source distribution. ===== scipy-1.18.1/subprojects/pyprima/pyprima/LICENCE.txt ===== BSD 3-Clause License Copyright (c) 2020--2025, Zaikun ZHANG ( https://www.zhangzk.net ) Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/subprojects/qhull_r/libqhull_r/COPYING.txt ===== Qhull, Copyright (c) 1993-2020 C.B. Barber Arlington, MA and The National Science and Technology Research Center for Computation and Visualization of Geometric Structures (The Geometry Center) University of Minnesota email: qhull@qhull.org This software includes Qhull from C.B. Barber and The Geometry Center. Files derived from Qhull 1.0 are copyrighted by the Geometry Center. The remaining files are copyrighted by C.B. Barber. Qhull is free software and may be obtained via http from www.qhull.org. It may be freely copied, modified, and redistributed under the following conditions: 1. All copyright notices must remain intact in all files. 2. A copy of this text file must be distributed along with any copies of Qhull that you redistribute; this includes copies that you have modified, or copies of programs or other software products that include Qhull. 3. If you modify Qhull, you must include a notice giving the name of the person performing the modification, the date of modification, and the reason for such modification. 4. When distributing modified versions of Qhull, or other software products that include Qhull, you must provide notice that the original source code may be obtained as noted above. 5. There is no warranty or other guarantee of fitness for Qhull, it is provided solely "as is". Bug reports or fixes may be sent to qhull_bug@qhull.org; the authors may or may not act on them as they desire. ===== scipy-1.18.1/subprojects/unuran/license.txt ===== Permission was granted on 23 March 2021 by Josef Leydold (josef.leydold@wu.ac.at) and Wolfgang Hoermann (hormannw@bound.edu.tr) to use the UNU.RAN code excluding the uniform random number generator in src/uniform/mrg31k3p.c (copyright by Renee Touzin) under the BSD licence. The following files have been removed from the original library: unuran-1.8.1/COPYING unuran-1.8.1/INSTALL unuran-1.8.1/Makefile.am unuran-1.8.1/Makefile.in unuran-1.8.1/NEWS unuran-1.8.1/UPGRADE unuran-1.8.1/acinclude.m4 unuran-1.8.1/aclocal.m4 unuran-1.8.1/autoconf/config.guess unuran-1.8.1/autoconf/config.sub unuran-1.8.1/autoconf/depcomp unuran-1.8.1/autoconf/install-sh unuran-1.8.1/autoconf/ltmain.sh unuran-1.8.1/autoconf/mdate-sh unuran-1.8.1/autoconf/missing unuran-1.8.1/autogen.sh unuran-1.8.1/configure unuran-1.8.1/doc/Makefile.am unuran-1.8.1/doc/Makefile.in unuran-1.8.1/doc/figures/Makefile.am unuran-1.8.1/doc/figures/Makefile.in unuran-1.8.1/doc/figures/discrete_inversion.eps unuran-1.8.1/doc/figures/discrete_inversion.pdf unuran-1.8.1/doc/figures/discrete_inversion.png unuran-1.8.1/doc/figures/discrete_inversion.txt unuran-1.8.1/doc/figures/inversion.eps unuran-1.8.1/doc/figures/inversion.pdf unuran-1.8.1/doc/figures/inversion.png unuran-1.8.1/doc/figures/inversion.txt unuran-1.8.1/doc/figures/rejection_from_constant_hat.eps unuran-1.8.1/doc/figures/rejection_from_constant_hat.pdf unuran-1.8.1/doc/figures/rejection_from_constant_hat.png unuran-1.8.1/doc/figures/rejection_from_constant_hat.txt unuran-1.8.1/doc/figures/rou-cauchy.eps unuran-1.8.1/doc/figures/rou-cauchy.pdf unuran-1.8.1/doc/figures/rou-cauchy.png unuran-1.8.1/doc/figures/rou-cauchy.txt unuran-1.8.1/doc/figures/texpict2img.pl unuran-1.8.1/doc/src/Makefile.am unuran-1.8.1/doc/src/Makefile.in unuran-1.8.1/doc/src/README unuran-1.8.1/doc/src/arvag.dh unuran-1.8.1/doc/src/error_debug.dh unuran-1.8.1/doc/src/examples.dh unuran-1.8.1/doc/src/glossary.dh unuran-1.8.1/doc/src/index.dh unuran-1.8.1/doc/src/installation.dh unuran-1.8.1/doc/src/intro.dh unuran-1.8.1/doc/src/make_texi.pl unuran-1.8.1/doc/src/methods.dh unuran-1.8.1/doc/src/methods_cemp.texi unuran-1.8.1/doc/src/methods_cont.texi unuran-1.8.1/doc/src/methods_cvec.texi unuran-1.8.1/doc/src/methods_cvemp.texi unuran-1.8.1/doc/src/methods_discr.texi unuran-1.8.1/doc/src/methods_matr.texi unuran-1.8.1/doc/src/methods_mcmc.texi unuran-1.8.1/doc/src/misc.dh unuran-1.8.1/doc/src/ref_example0.texi unuran-1.8.1/doc/src/ref_example0_str.texi unuran-1.8.1/doc/src/ref_example1.texi unuran-1.8.1/doc/src/ref_example1_str.texi unuran-1.8.1/doc/src/ref_example2.texi unuran-1.8.1/doc/src/ref_example2_str.texi unuran-1.8.1/doc/src/ref_example3.texi unuran-1.8.1/doc/src/ref_example3_str.texi unuran-1.8.1/doc/src/ref_example_FuncStr.texi unuran-1.8.1/doc/src/ref_example_anti.texi unuran-1.8.1/doc/src/ref_example_anti_str.texi unuran-1.8.1/doc/src/ref_example_cext.texi unuran-1.8.1/doc/src/ref_example_cont.texi unuran-1.8.1/doc/src/ref_example_cont_str.texi unuran-1.8.1/doc/src/ref_example_dext.texi unuran-1.8.1/doc/src/ref_example_discr.texi unuran-1.8.1/doc/src/ref_example_discr_str.texi unuran-1.8.1/doc/src/ref_example_emp.texi unuran-1.8.1/doc/src/ref_example_emp_str.texi unuran-1.8.1/doc/src/ref_example_gsl.texi unuran-1.8.1/doc/src/ref_example_mixt.texi unuran-1.8.1/doc/src/ref_example_mixt_inv.texi unuran-1.8.1/doc/src/ref_example_rngstreams.texi unuran-1.8.1/doc/src/ref_example_vemp.texi unuran-1.8.1/doc/src/ref_example_win32_1.texi unuran-1.8.1/doc/src/ref_example_win32_2.texi unuran-1.8.1/doc/src/references.dh unuran-1.8.1/doc/src/stddist.dh unuran-1.8.1/doc/src/texinfo.tex unuran-1.8.1/doc/src/top.dh unuran-1.8.1/doc/src/unuran_src.texi unuran-1.8.1/doc/stamp-1 unuran-1.8.1/doc/stamp-vti unuran-1.8.1/doc/unuran.html unuran-1.8.1/doc/unuran.info unuran-1.8.1/doc/unuran.pdf unuran-1.8.1/doc/unuran.texi unuran-1.8.1/doc/unuran.txt unuran-1.8.1/doc/unuran_win32.html unuran-1.8.1/doc/unuran_win32.info unuran-1.8.1/doc/unuran_win32.pdf unuran-1.8.1/doc/unuran_win32.texi unuran-1.8.1/doc/unuran_win32_libname.texi unuran-1.8.1/doc/version.texi unuran-1.8.1/doc/version_win32.texi unuran-1.8.1/examples/Makefile.am unuran-1.8.1/examples/Makefile.in unuran-1.8.1/examples/example0.c unuran-1.8.1/examples/example0_str.c unuran-1.8.1/examples/example1.c unuran-1.8.1/examples/example1_str.c unuran-1.8.1/examples/example2.c unuran-1.8.1/examples/example2_str.c unuran-1.8.1/examples/example3.c unuran-1.8.1/examples/example3_str.c unuran-1.8.1/examples/example_FuncStr.c unuran-1.8.1/examples/example_anti.c unuran-1.8.1/examples/example_anti_str.c unuran-1.8.1/examples/example_cext.c unuran-1.8.1/examples/example_cont.c unuran-1.8.1/examples/example_cont_str.c unuran-1.8.1/examples/example_cpp.cpp unuran-1.8.1/examples/example_dext.c unuran-1.8.1/examples/example_discr.c unuran-1.8.1/examples/example_discr_str.c unuran-1.8.1/examples/example_emp.c unuran-1.8.1/examples/example_emp_str.c unuran-1.8.1/examples/example_errorhandler.c unuran-1.8.1/examples/example_gsl.c unuran-1.8.1/examples/example_mcorr.c unuran-1.8.1/examples/example_mixt.c unuran-1.8.1/examples/example_mixt_inv.c unuran-1.8.1/examples/example_prng.c unuran-1.8.1/examples/example_prng_str.c unuran-1.8.1/examples/example_rngstreams.c unuran-1.8.1/examples/example_vemp.c unuran-1.8.1/experiments/Makefile.am unuran-1.8.1/experiments/Makefile.in unuran-1.8.1/experiments/make_timing_table.c unuran-1.8.1/experiments/timings_chap_10.conf unuran-1.8.1/experiments/timings_chap_8.conf unuran-1.8.1/scripts/Makefile.am unuran-1.8.1/scripts/Makefile.in unuran-1.8.1/scripts/compute_machine_constants.c unuran-1.8.1/scripts/make_urng_header.pl unuran-1.8.1/scripts/merge_h.pl unuran-1.8.1/scripts/read_PDF.pl unuran-1.8.1/scripts/remove_comments.pl unuran-1.8.1/scripts/run_valgrind.sh unuran-1.8.1/scripts/win32/Makefile.win32 unuran-1.8.1/scripts/win32/build.sh unuran-1.8.1/scripts/win32/example1.c unuran-1.8.1/scripts/win32/example2.c unuran-1.8.1/scripts/win32/make_dll_def.pl unuran-1.8.1/src/Makefile.am unuran-1.8.1/src/Makefile.in unuran-1.8.1/src/README unuran-1.8.1/src/distr/Makefile.am unuran-1.8.1/src/distr/Makefile.in unuran-1.8.1/src/distributions/Makefile.am unuran-1.8.1/src/distributions/Makefile.in unuran-1.8.1/src/distributions/README unuran-1.8.1/src/methods/Makefile.am unuran-1.8.1/src/methods/Makefile.in unuran-1.8.1/src/parser/Makefile.am unuran-1.8.1/src/parser/Makefile.in unuran-1.8.1/src/parser/functparser_doc.dh unuran-1.8.1/src/parser/make_stringparser.pl unuran-1.8.1/src/parser/stringparser_doc.dh unuran-1.8.1/src/specfunct/Makefile.am unuran-1.8.1/src/specfunct/Makefile.in unuran-1.8.1/src/tests/Makefile.am unuran-1.8.1/src/tests/Makefile.in unuran-1.8.1/src/uniform/Makefile.am unuran-1.8.1/src/uniform/Makefile.in unuran-1.8.1/src/uniform/fish.c unuran-1.8.1/src/uniform/mrg31k3p.c unuran-1.8.1/src/uniform/mstd.c unuran-1.8.1/src/uniform/unuran_urng_gsl.h unuran-1.8.1/src/uniform/unuran_urng_prng.h unuran-1.8.1/src/uniform/unuran_urng_rngstreams.h unuran-1.8.1/src/uniform/urng_builtin.c unuran-1.8.1/src/uniform/urng_builtin.h unuran-1.8.1/src/uniform/urng_fvoid.c unuran-1.8.1/src/uniform/urng_fvoid.h unuran-1.8.1/src/uniform/urng_gsl.c unuran-1.8.1/src/uniform/urng_gsl.h unuran-1.8.1/src/uniform/urng_gslqrng.c unuran-1.8.1/src/uniform/urng_gslqrng.h unuran-1.8.1/src/uniform/urng_prng.c unuran-1.8.1/src/uniform/urng_prng.h unuran-1.8.1/src/uniform/urng_randomshift.c unuran-1.8.1/src/uniform/urng_randomshift.h unuran-1.8.1/src/uniform/urng_rngstreams.c unuran-1.8.1/src/uniform/urng_rngstreams.h unuran-1.8.1/src/unuran.h.in unuran-1.8.1/src/urng/Makefile.am unuran-1.8.1/src/urng/Makefile.in unuran-1.8.1/src/utils/Makefile.am unuran-1.8.1/src/utils/Makefile.in unuran-1.8.1/tests/Makefile.am unuran-1.8.1/tests/Makefile.in unuran-1.8.1/tests/README.conf unuran-1.8.1/tests/make_test_files.pl unuran-1.8.1/tests/t_arou.c unuran-1.8.1/tests/t_arou.conf unuran-1.8.1/tests/t_ars.c unuran-1.8.1/tests/t_ars.conf unuran-1.8.1/tests/t_auto.c unuran-1.8.1/tests/t_auto.conf unuran-1.8.1/tests/t_cext.c unuran-1.8.1/tests/t_cext.conf unuran-1.8.1/tests/t_cstd.c unuran-1.8.1/tests/t_cstd.conf unuran-1.8.1/tests/t_dari.c unuran-1.8.1/tests/t_dari.conf unuran-1.8.1/tests/t_dau.c unuran-1.8.1/tests/t_dau.conf unuran-1.8.1/tests/t_deprecated_tdrgw.c unuran-1.8.1/tests/t_deprecated_vmt.c unuran-1.8.1/tests/t_deprecated_vmt.conf unuran-1.8.1/tests/t_dext.c unuran-1.8.1/tests/t_dext.conf unuran-1.8.1/tests/t_dgt.c unuran-1.8.1/tests/t_dgt.conf unuran-1.8.1/tests/t_distr.c unuran-1.8.1/tests/t_distr.conf unuran-1.8.1/tests/t_distr_cemp.c unuran-1.8.1/tests/t_distr_cemp.conf unuran-1.8.1/tests/t_distr_condi.c unuran-1.8.1/tests/t_distr_condi.conf unuran-1.8.1/tests/t_distr_cont.c unuran-1.8.1/tests/t_distr_cont.conf unuran-1.8.1/tests/t_distr_corder.c unuran-1.8.1/tests/t_distr_corder.conf unuran-1.8.1/tests/t_distr_cvec.c unuran-1.8.1/tests/t_distr_cvec.conf unuran-1.8.1/tests/t_distr_cvemp.c unuran-1.8.1/tests/t_distr_cvemp.conf unuran-1.8.1/tests/t_distr_cxtrans.c unuran-1.8.1/tests/t_distr_cxtrans.conf unuran-1.8.1/tests/t_distr_discr.c unuran-1.8.1/tests/t_distr_discr.conf unuran-1.8.1/tests/t_distr_matr.c unuran-1.8.1/tests/t_distr_matr.conf unuran-1.8.1/tests/t_dsrou.c unuran-1.8.1/tests/t_dsrou.conf unuran-1.8.1/tests/t_dss.c unuran-1.8.1/tests/t_dss.conf unuran-1.8.1/tests/t_dstd.c unuran-1.8.1/tests/t_dstd.conf unuran-1.8.1/tests/t_empk.c unuran-1.8.1/tests/t_empk.conf unuran-1.8.1/tests/t_empl.c unuran-1.8.1/tests/t_empl.conf unuran-1.8.1/tests/t_gibbs.c unuran-1.8.1/tests/t_gibbs.conf unuran-1.8.1/tests/t_hinv.c unuran-1.8.1/tests/t_hinv.conf unuran-1.8.1/tests/t_hist.c unuran-1.8.1/tests/t_hist.conf unuran-1.8.1/tests/t_hitro.c unuran-1.8.1/tests/t_hrb.c unuran-1.8.1/tests/t_hrb.conf unuran-1.8.1/tests/t_hrd.c unuran-1.8.1/tests/t_hrd.conf unuran-1.8.1/tests/t_hri.c unuran-1.8.1/tests/t_hri.conf unuran-1.8.1/tests/t_info.c unuran-1.8.1/tests/t_info.conf unuran-1.8.1/tests/t_itdr.c unuran-1.8.1/tests/t_itdr.conf unuran-1.8.1/tests/t_mcorr.c unuran-1.8.1/tests/t_mcorr.conf unuran-1.8.1/tests/t_mixt.c unuran-1.8.1/tests/t_mvstd.c unuran-1.8.1/tests/t_mvstd.conf unuran-1.8.1/tests/t_mvtdr.c unuran-1.8.1/tests/t_ninv.c unuran-1.8.1/tests/t_ninv.conf unuran-1.8.1/tests/t_norta.c unuran-1.8.1/tests/t_norta.conf unuran-1.8.1/tests/t_nrou.c unuran-1.8.1/tests/t_nrou.conf unuran-1.8.1/tests/t_pinv unuran-1.8.1/tests/t_pinv.c unuran-1.8.1/tests/t_srou.c unuran-1.8.1/tests/t_srou.conf unuran-1.8.1/tests/t_ssr.c unuran-1.8.1/tests/t_ssr.conf unuran-1.8.1/tests/t_stringparser.c unuran-1.8.1/tests/t_stringparser.conf unuran-1.8.1/tests/t_tabl.c unuran-1.8.1/tests/t_tabl.conf unuran-1.8.1/tests/t_tdr.c unuran-1.8.1/tests/t_tdr.conf unuran-1.8.1/tests/t_tdr_gw.c unuran-1.8.1/tests/t_tdr_gw.conf unuran-1.8.1/tests/t_tdr_ia.c unuran-1.8.1/tests/t_tdr_ia.conf unuran-1.8.1/tests/t_tdr_ps.c unuran-1.8.1/tests/t_tdr_ps.conf unuran-1.8.1/tests/t_unif.c unuran-1.8.1/tests/t_unif.conf unuran-1.8.1/tests/t_unur_error.c unuran-1.8.1/tests/t_unur_error.conf unuran-1.8.1/tests/t_unur_tests.c unuran-1.8.1/tests/t_unur_tests.conf unuran-1.8.1/tests/t_urng.c unuran-1.8.1/tests/t_urng.conf unuran-1.8.1/tests/t_utdr.c unuran-1.8.1/tests/t_utdr.conf unuran-1.8.1/tests/t_util_matrix.c unuran-1.8.1/tests/t_util_matrix.conf unuran-1.8.1/tests/t_vempk.c unuran-1.8.1/tests/t_vempk.conf unuran-1.8.1/tests/t_vnrou.c unuran-1.8.1/tests/t_vnrou.conf unuran-1.8.1/tests/t_x_gen.c unuran-1.8.1/tests/t_x_gen.conf unuran-1.8.1/tests/test_StdDistr.c unuran-1.8.1/tests/test_StdDistr.m unuran-1.8.1/tests/test_functionparser.c unuran-1.8.1/tests/test_functionparser.m unuran-1.8.1/tests/testcounter.c unuran-1.8.1/tests/testdistributions/Makefile.am unuran-1.8.1/tests/testdistributions/Makefile.in unuran-1.8.1/tests/testdistributions/c_sawtooth_contpdf.c unuran-1.8.1/tests/testdistributions/c_sawtooth_discpdf.c unuran-1.8.1/tests/testdistributions/c_w_pdf_from_logpdf.c unuran-1.8.1/tests/testdistributions/c_with_multiple_of_pdf.c unuran-1.8.1/tests/testdistributions/c_wo_logpdf.c unuran-1.8.1/tests/testdistributions/testdistributions.h unuran-1.8.1/tests/testdistributions/vc_cauchy_RoU_ball.c unuran-1.8.1/tests/testdistributions/vc_special_correlations_ar1.c unuran-1.8.1/tests/testdistributions/vc_special_correlations_constantrho.c unuran-1.8.1/tests/testdistributions/vc_w_marginals.c unuran-1.8.1/tests/testdistributions/vc_wo_logpdf.c unuran-1.8.1/tests/testroutines.c unuran-1.8.1/tests/testunuran.h unuran-1.8.1/src/distr/deprecated_distr.c unuran-1.8.1/src/distr/deprecated_distr.h unuran-1.8.1/src/methods/deprecated_methods.c unuran-1.8.1/src/methods/deprecated_methods.h unuran-1.8.1/src/methods/deprecated_tdrgw.h unuran-1.8.1/src/methods/deprecated_vmt.c unuran-1.8.1/src/methods/deprecated_vmt.h unuran-1.8.1/src/methods/deprecated_vmt_struct.h unuran-1.8.1/src/parser/stringparser_lists.ch.in UNU.RAN relies on a modified version of the library Cephes, see ./src/src/specfunct/cephes_source.h ===== scipy-1.18.1/subprojects/xsf/LICENSE ===== BSD 3-Clause License Copyright (c) 2024, SciPy Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/subprojects/xsf/LICENSES_bundled.txt ===== xsf bundles several libraries that are compatibly licensed. We list these here. Name: Cephes Files: include/xsf/cephes/* License: BSD-3-Clause Distributed under 3-clause BSD license with permission from the author, see https://lists.debian.org/debian-legal/2004/12/msg00295.html Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1995, 2000 by Stephen L. Moshier This software is derived from the Cephes Math Library and is incorporated herein by permission of the author. All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. * Neither the name of the nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. Name: Faddeeva Files: include/xsf/faddeeva.h License: MIT Copyright (c) 2012 Massachusetts Institute of Technology Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. Name: qd Files: include/xsf/cephes/dd_real.h License: BSD-3-Clause-LBNL This work was supported by the Director, Office of Science, Division of Mathematical, Information, and Computational Sciences of the U.S. Department of Energy under contract numbers DE-AC03-76SF00098 and DE-AC02-05CH11231. Copyright (c) 2003-2009, The Regents of the University of California, through Lawrence Berkeley National Laboratory (subject to receipt of any required approvals from U.S. Dept. of Energy) All rights reserved. 1. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: (1) Redistributions of source code must retain the copyright notice, this list of conditions and the following disclaimer. (2) Redistributions in binary form must reproduce the copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory, U.S. Dept. of Energy nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. 2. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 3. You are under no obligation whatsoever to provide any bug fixes, patches, or upgrades to the features, functionality or performance of the source code ("Enhancements") to anyone; however, if you choose to make your Enhancements available either publicly, or directly to Lawrence Berkeley National Laboratory, without imposing a separate written license agreement for such Enhancements, then you hereby grant the following license: a non-exclusive, royalty-free perpetual license to install, use, modify, prepare derivative works, incorporate into other computer software, distribute, and sublicense such enhancements or derivative works thereof, in binary and source code form. Name: mdspan Files: include/xsf/third_party/kokkos/mdspan.hpp License: Apache-2.0 WITH LLVM-exception For details, see include/xsf/third_party/kokkos/mdspan.hpp ===== scipy-1.18.1/scipy/io/_fast_matrix_market/src/pystreambuf.h ===== // SPDX-License-Identifier: BSD-3-Clause /* Based on https://gist.github.com/asford/544323a5da7dddad2c9174490eb5ed06 Original license text --------------------- This component utilizes components derived from cctbx, available at http://cci.lbl.gov/cctbx_sources/boost_adaptbx/python_streambuf.h *** License agreement *** cctbx Copyright (c) 2006, The Regents of the University of California, through Lawrence Berkeley National Laboratory (subject to receipt of any required approvals from the U.S. Dept. of Energy). All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: (1) Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. (2) Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory, U.S. Dept. of Energy nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. You are under no obligation whatsoever to provide any bug fixes, patches, or upgrades to the features, functionality or performance of the source code ("Enhancements") to anyone; however, if you choose to make your Enhancements available either publicly, or directly to Lawrence Berkeley National Laboratory, without imposing a separate written license agreement for such Enhancements, then you hereby grant the following license: a non-exclusive, royalty-free perpetual license to install, use, modify, prepare derivative works, incorporate into other computer software, distribute, and sublicense such enhancements or derivative works thereof, in binary and source code form. */ ===== scipy-1.18.1/scipy/interpolate/_bary_rational.py ===== # Copyright (c) 2017, The Chancellor, Masters and Scholars of the University # of Oxford, and the Chebfun Developers. All rights reserved. # # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: # * Redistributions of source code must retain the above copyright # notice, this list of conditions and the following disclaimer. # * Redistributions in binary form must reproduce the above copyright # notice, this list of conditions and the following disclaimer in the # documentation and/or other materials provided with the distribution. # * Neither the name of the University of Oxford nor the names of its # contributors may be used to endorse or promote products derived from # this software without specific prior written permission. # # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED # WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE # DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR # ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES # (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND # ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ===== scipy-1.18.1/scipy/integrate/_lebedev.py ===== # getLebedevSphere # Copyright (c) 2010, Robert Parrish # All rights reserved. # # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are # met: # # * Redistributions of source code must retain the above copyright # notice, this list of conditions and the following disclaimer. # * Redistributions in binary form must reproduce the above copyright # notice, this list of conditions and the following disclaimer in # the documentation and/or other materials provided with the distribution # # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE # ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE # LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF # SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS # INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN # CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. # # Brainlessly translated to Python ===== scipy-1.18.1/scipy/stats/libnpyrandom/logfactorial.c ===== /* Copyright (c) 2005-2022, NumPy Developers. All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. * Neither the name of the NumPy Developers nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ ===== scipy-1.18.1/subprojects/xsf/include/xsf/third_party/kokkos/mdspan.hpp ===== #ifndef MDSPAN_SINGLE_HEADER_INCLUDE_GUARD_ #define MDSPAN_SINGLE_HEADER_INCLUDE_GUARD_ // BEGIN_FILE_INCLUDE: mdspan/include/mdspan/mdspan.hpp //@HEADER // ************************************************************************ // // Kokkos v. 4.0 // Copyright (2022) National Technology & Engineering // Solutions of Sandia, LLC (NTESS). // // Under the terms of Contract DE-NA0003525 with NTESS, // the U.S. Government retains certain rights in this software. // // Part of Kokkos, under the Apache License v2.0 with LLVM Exceptions. // See https://kokkos.org/LICENSE for license information. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // //@HEADER ===== scipy-1.18.1/subprojects/unuran/unuran/src/specfunct/cephes_source.h ===== /* This file was considerably changed for UNU.RAN. As a consequence it has been renamed from "mconf.h" to "cephes_source.h" to avoid confusion. We were only interested in the files: gamma.c igamma.c incbet.c ndtr.c ndtri.c to use these files we needed the auxiliary files: isnan.c mtherr.c polevl.c The aim was to enhance portability for ANSI C without expert knowledge of the floating point unit. The main changes are numeric constants (in cephes - const.c) are replaced by macros defined in this file (source_mconf.h). We changed these constants using the constants provided by ANSI C in the math.h file. We changed also the definitions of MAXGAM (in incbet.c) and of MAXSTIR (in gamma.c). Everything concerning NANS and INFINITY was moved into isnan.c. (Also the defines to turn it off or on.) We have only tested this version with NAN and INFINITY turned off, but we have not changed anything concerning that question. So turning this support on should (could) work. March 26th, 2001, Josef Leydold and Wolfgang Hoermann. */ /*---------------------------------------------------------------------------*/ /* mconf.h * * Common include file for math routines * * * * SYNOPSIS: * * #include "mconf.h" * * * * DESCRIPTION: * * This file contains definitions for error codes that are * passed to the common error handling routine mtherr() * (which see). * * The file also includes a conditional assembly definition * for the type of computer arithmetic (IEEE, DEC, Motorola * IEEE, or UNKnown). * * For Digital Equipment PDP-11 and VAX computers, certain * IBM systems, and others that use numbers with a 56-bit * significand, the symbol DEC should be defined. In this * mode, most floating point constants are given as arrays * of octal integers to eliminate decimal to binary conversion * errors that might be introduced by the compiler. * * For little-endian computers, such as IBM PC, that follow the * IEEE Standard for Binary Floating Point Arithmetic (ANSI/IEEE * Std 754-1985), the symbol IBMPC should be defined. These * numbers have 53-bit significands. In this mode, constants * are provided as arrays of hexadecimal 16 bit integers. * * Big-endian IEEE format is denoted MIEEE. On some RISC * systems such as Sun SPARC, double precision constants * must be stored on 8-byte address boundaries. Since integer * arrays may be aligned differently, the MIEEE configuration * may fail on such machines. * * To accommodate other types of computer arithmetic, all * constants are also provided in a normal decimal radix * which one can hope are correctly converted to a suitable * format by the available C language compiler. To invoke * this mode, define the symbol UNK. * * An important difference among these modes is a predefined * set of machine arithmetic constants for each. The numbers * MACHEP (the machine roundoff error), MAXNUM (largest number * represented), and several other parameters are preset by * the configuration symbol. Check the file const.c to * ensure that these values are correct for your computer. * * Configurations NANS, INFINITIES, MINUSZERO, and DENORMAL * may fail on many systems. Verify that they are supposed * to work on your computer. */ /* Cephes Math Library Release 2.3: June, 1995 Copyright 1984, 1987, 1989, 1995 by Stephen L. Moshier */ /*---------------------------------------------------------------------------*/ /* Include constants used in UNU.RAN. */ ===== scipy-1.18.1/scipy/optimize/src/lbfgsb.c ===== #include "lbfgsb.h" enum Status { START = 0, NEW_X = 1, RESTART = 2, // RESTART FROM LINESEARCH FG = 3, // REQUEST FG FOR THE CURRENT X CONVERGENCE = 4, STOP = 5, WARNING = 6, ERROR = 7, ABNORMAL = 8 // ABNORMAL TERMINATION IN LINE SEARCH (WHY NOT AN ERROR?) }; enum StatusMsg { NO_MSG = 0, FG_START = 301, FG_LNSRCH = 302, CONV_GRAD = 401, // CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL CONV_F = 402, // CONVERGENCE: PROGRESS <= FACTR*EPSMCH STOP_CPU = 501, // CPU EXCEEDING THE TIME LIMIT ?! STOP_ITER = 502, // TOTAL NO. OF F,G EVALUATIONS EXCEEDS LIMIT STOP_GRAD = 503, // PROJECTED GRADIENT IS SUFFICIENTLY SMALL STOP_ITERC = 504, // TOTAL NO. of ITERATIONS REACHED LIMIT STOP_CALLB = 505, // CALLBACK REQUESTED HALT WARN_ROUND = 601, // ROUNDING ERRORS PREVENT PROGRESS WARN_STPMAX = 602, // STP = STPMAX WARN_STPMIN = 603, // STP = STPMIN WARN_XTOL = 604, // XTOL TEST SATISFIED ?! ERROR_NOFEAS = 701, // NO FEASIBLE SOLUTION ERROR_FACTR = 702, // FACTR < 0 ERROR_FTOL = 703, // FTOL < 0 ERROR_GTOL = 704, // GTOL < 0 ERROR_XTOL = 705, // XTOL < 0 ERROR_STP1 = 706, // STP < STPMIN ERROR_STP2 = 707, // STP > STPMAX ERROR_STPMIN = 708, // STPMIN < 0 ERROR_STPMAX = 709, // STPMAX < STPMIN ERROR_INITG = 710, // INITIAL G >= 0 ERROR_M = 711, // M <= 0 ERROR_N = 712, // N <= 0 ERROR_NBD = 713 // INVALID NBD }; // Internal functions static void mainlb( CBLAS_INT n, CBLAS_INT m, double* x, double* l, double* u, CBLAS_INT* nbd, double* f, double* g, double factr, double pgtol, double* ws, double* wy, double* sy, double* ss, double* wt, double* wn, double* snd, double* z, double* r, double* d, double* t, double* xp, double* wa, CBLAS_INT* index, CBLAS_INT* iwhere, CBLAS_INT* indx2, CBLAS_INT* task, CBLAS_INT* task_msg, CBLAS_INT* lsave, CBLAS_INT* isave, double* dsave, CBLAS_INT maxls, CBLAS_INT* ln_task, CBLAS_INT* ln_taskmsg ); static void active( CBLAS_INT n, double* l, double* u, CBLAS_INT* nbd, double* x, CBLAS_INT* iwhere, CBLAS_INT* prjctd, CBLAS_INT* cnstnd, CBLAS_INT* boxed ); static void bmv( CBLAS_INT m, double* sy, double* wt, CBLAS_INT col, double* v, double* p, CBLAS_INT* info ); static void cauchy( CBLAS_INT n, double* x, double* l, double* u, CBLAS_INT* nbd, double *g, CBLAS_INT* iorder, CBLAS_INT* iwhere, double* t, double* d, double* xcp, CBLAS_INT m, double* wy, double* ws, double* sy, double* wt, double theta, CBLAS_INT col, CBLAS_INT head, double* p, double* c, double* wbp, double* v, CBLAS_INT* nseg, double sbgnrm, CBLAS_INT* info ); static void cmprlb( CBLAS_INT n, CBLAS_INT m, double* x, double* g, double* ws, double* wy, double* sy, double* wt, double* z, double* r, double* wa, CBLAS_INT* index, double theta, CBLAS_INT col, CBLAS_INT head, CBLAS_INT nfree, CBLAS_INT cnstnd, CBLAS_INT* info ); static void errclb( CBLAS_INT n, CBLAS_INT m, double factr, double* l, double* u, CBLAS_INT* nbd, CBLAS_INT* task, CBLAS_INT* task_msg, CBLAS_INT* info, CBLAS_INT* k ); static void formk( CBLAS_INT n, CBLAS_INT nsub, CBLAS_INT* ind, CBLAS_INT nenter, CBLAS_INT ileave, CBLAS_INT* indx2, CBLAS_INT iupdat, CBLAS_INT updatd, double* wn, double* wn1, CBLAS_INT m, double* ws, double* wy, double* sy, double theta, CBLAS_INT col, CBLAS_INT head, CBLAS_INT* info ); static void formt( CBLAS_INT m, double* wt, double* sy, double* ss, CBLAS_INT col, double theta, CBLAS_INT* info ); static void freev( CBLAS_INT n, CBLAS_INT* nfree, CBLAS_INT* idx, CBLAS_INT* nenter, CBLAS_INT* ileave, CBLAS_INT* idx2, CBLAS_INT* iwhere, CBLAS_INT* wrk, CBLAS_INT updatd, CBLAS_INT cnstnd, CBLAS_INT iter ); static void hpsolb(CBLAS_INT n, double* t, CBLAS_INT* iorder, CBLAS_INT iheap); static void lnsrlb( CBLAS_INT n, double* l, double* u, CBLAS_INT* nbd, double* x, double f, double* fold, double *gd, double *gdold, double* g, double* d, double* r, double* t, double* z, double* stp, double* dnorm, double* dtd, double* xstep, double* stpmx, CBLAS_INT iter, CBLAS_INT* ifun, CBLAS_INT* iback, CBLAS_INT* nfgv, CBLAS_INT* info, CBLAS_INT* task, CBLAS_INT* task_msg, CBLAS_INT boxed, CBLAS_INT cnstnd, CBLAS_INT* isave, double* dsave, CBLAS_INT* temp_task, CBLAS_INT* temp_task_msg ); static void matupd( CBLAS_INT n, CBLAS_INT m, double* ws, double *wy, double* sy, double* ss, double* d, double* r, CBLAS_INT* itail, CBLAS_INT iupdat, CBLAS_INT* col, CBLAS_INT* head, double* theta, double rr, double dr, double stp, double dtd ); static void projgr( CBLAS_INT n, double* l, double* u, CBLAS_INT* nbd, double* x, double* g, double* sbgnrm ); static void subsm( CBLAS_INT n, CBLAS_INT m, CBLAS_INT nsub, CBLAS_INT* ind, double* l, double* u, CBLAS_INT* nbd, double* x, double* d, double* xp, double* ws, double* wy, double theta, double* xx, double* gg, CBLAS_INT col, CBLAS_INT head, CBLAS_INT* iword, double* wv, double* wn, CBLAS_INT* info ); // Line search functions static void dcsrch( double f, double g, double* stp, double ftol, double gtol, double xtol, double stpmin, double stpmax, CBLAS_INT* task, CBLAS_INT* task_msg, CBLAS_INT* isave, double* dsave ); static void dcstep( double* stx, double* fx, double* dx, double* sty, double* fy, double* dy, double* stp, double fp, double dp, CBLAS_INT* brackt, double stpmin, double stpmax ); // Helper functions static inline void save_local_vars( CBLAS_INT prjctd, CBLAS_INT cnstnd, CBLAS_INT boxed, CBLAS_INT updatd, CBLAS_INT nintol, CBLAS_INT iback, CBLAS_INT nskip, CBLAS_INT head, CBLAS_INT col, CBLAS_INT itail, CBLAS_INT iter, CBLAS_INT iupdat, CBLAS_INT nseg, CBLAS_INT nfgv, CBLAS_INT info, CBLAS_INT ifun, CBLAS_INT iword, CBLAS_INT nfree, CBLAS_INT nact, CBLAS_INT ileave, CBLAS_INT nenter, double theta, double fold, double tol, double dnorm, double gd, double stpmx, double sbgnrm, double stp, double gdold, double dtd, CBLAS_INT* lsave, CBLAS_INT* isave, double* dsave ); static inline void save_vars( CBLAS_INT brackt, CBLAS_INT stage, double ginit, double gtest, double gx, double gy, double finit, double fx, double fy, double stx, double sty, double stmin, double stmax, double width, double width1, CBLAS_INT* isave, double* dsave ); static double epsmach = 2.220446049250313e-016; /* np.finfo(np.float64).eps */ void setulb(CBLAS_INT n, CBLAS_INT m, double* x, double* l, double* u, CBLAS_INT* nbd, double* f, double* g, double factr, double pgtol, double* wa, CBLAS_INT* iwa, CBLAS_INT* task, CBLAS_INT* lsave, CBLAS_INT* isave, double* dsave, CBLAS_INT maxls, CBLAS_INT* ln_task) { // ************ // // Subroutine setulb // // This subroutine partitions the working arrays wa and iwa, and // then uses the limited memory BFGS method to solve the bound // constrained optimization problem by calling mainlb. // (The direct method will be used in the subspace minimization.) // // n is an integer variable. // On entry n is the dimension of the problem. // On exit n is unchanged. // // m is an integer variable. // On entry m is the maximum number of variable metric corrections // used to define the limited memory matrix. // On exit m is unchanged. // // x is a double precision array of dimension n. // On entry x is an approximation to the solution. // On exit x is the current approximation. // // l is a double precision array of dimension n. // On entry l is the lower bound on x. // On exit l is unchanged. // // u is a double precision array of dimension n. // On entry u is the upper bound on x. // On exit u is unchanged. // // nbd is an integer array of dimension n. // On entry nbd represents the type of bounds imposed on the // variables, and must be specified as follows: // nbd(i)=0 if x(i) is unbounded, // 1 if x(i) has only a lower bound, // 2 if x(i) has both lower and upper bounds, and // 3 if x(i) has only an upper bound. // On exit nbd is unchanged. // // f is a double precision variable. // On first entry f is unspecified. // On final exit f is the value of the function at x. // // g is a double precision array of dimension n. // On first entry g is unspecified. // On final exit g is the value of the gradient at x. // // factr is a double precision variable. // On entry factr >= 0 is specified by the user. The iteration // will stop when // // (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch // // where epsmch is the machine precision, which is automatically // generated by the code. Typical values for factr: 1.d+12 for // low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely // high accuracy. // On exit factr is unchanged. // // pgtol is a double precision variable. // On entry pgtol >= 0 is specified by the user. The iteration // will stop when // // max{|proj g_i | i = 1, ..., n} <= pgtol // // where pg_i is the ith component of the projected gradient. // On exit pgtol is unchanged. // // wa is a double precision working array of length // (2mmax + 5)nmax + 12mmax^2 + 12mmax. // // iwa is an integer working array of length 3nmax. // // task is a working integer array of length 2. task[0] indicates // the current job when entering and quitting this subroutine. // task[1] indicates the status message. // // iprint is an integer variable that must be set by the user. // It controls the frequency and type of output generated: // iprint<0 no output is generated; // iprint=0 print only one line at the last iteration; // 0100 print details of every iteration including x and g; // // // csave is a working string of characters of length 60. // // lsave is a logical working array of dimension 4. // On exit with 'task' = NEW_X, the following information is // available: // If lsave(1) = .true. then the initial X has been replaced by // its projection in the feasible set; // If lsave(2) = .true. then the problem is constrained; // If lsave(3) = .true. then each variable has upper and lower // bounds; // // isave is an integer working array of dimension 44. // On exit with 'task' = NEW_X, the following information is // available: // isave(22) = the total number of intervals explored in the // search of Cauchy points; // isave(26) = the total number of skipped BFGS updates before // the current iteration; // isave(30) = the number of current iteration; // isave(31) = the total number of BFGS updates prior the current // iteration; // isave(33) = the number of intervals explored in the search of // Cauchy point in the current iteration; // isave(34) = the total number of function and gradient // evaluations; // isave(36) = the number of function value or gradient // evaluations in the current iteration; // if isave(37) = 0 then the subspace argmin is within the box; // if isave(37) = 1 then the subspace argmin is beyond the box; // isave(38) = the number of free variables in the current // iteration; // isave(39) = the number of active constraints in the current // iteration; // n + 1 - isave(40) = the number of variables leaving the set of // active constraints in the current iteration; // isave(41) = the number of variables entering the set of active // constraints in the current iteration. // // dsave is a double precision working array of dimension 29. // On exit with 'task' = NEW_X, the following information is // available: // dsave(1) = current 'theta' in the BFGS matrix; // dsave(2) = f(x) in the previous iteration; // dsave(3) = factr*epsmch; // dsave(4) = 2-norm of the line search direction vector; // dsave(5) = the machine precision epsmch generated by the code; // dsave(7) = the accumulated time spent on searching for // Cauchy points; // dsave(8) = the accumulated time spent on // subspace minimization; // dsave(9) = the accumulated time spent on line search; // dsave(11) = the slope of the line search function at // the current point of line search; // dsave(12) = the maximum relative step length imposed in // line search; // dsave(13) = the infinity norm of the projected gradient; // dsave(14) = the relative step length in the line search; // dsave(15) = the slope of the line search function at // the starting point of the line search; // dsave(16) = the square of the 2-norm of the line search // direction vector. // // Subprograms called: // // L-BFGS-B Library ... mainlb. // // // References: // // [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited // memory algorithm for bound constrained optimization'', // SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. // // [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a // limited memory FORTRAN code for solving bound constrained // optimization problems'', Tech. Report, NAM-11, EECS Department, // Northwestern University, 1994. // // (Postscript files of these papers are available via anonymous // ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ // CBLAS_INT lws, lr, lz, lt, ld, lxp, lwa, lwy, lsy, lss, lwt, lwn, lsnd; if (task[0] == START) { // Pointer olympics isave[0] = m*n; isave[1] = m*m; isave[2] = 4*m*m; isave[3] = 0; // ws m*n isave[4] = isave[3] + isave[0]; // wy m*n isave[5] = isave[4] + isave[0]; // wsy m**2 isave[6] = isave[5] + isave[1]; // wss m**2 isave[7] = isave[6] + isave[1]; // wt m**2 isave[8] = isave[7] + isave[1]; // wn 4*m**2 isave[9] = isave[8] + isave[2]; // wsnd 4*m**2 isave[10] = isave[9] + isave[2]; // wz n isave[11] = isave[10] + n; // wr n isave[12] = isave[11] + n; // wd n isave[13] = isave[12] + n; // wt n isave[14] = isave[13] + n; // wxp n isave[15] = isave[14] + n; // wa 8*m } lws = isave[3]; lwy = isave[4]; lsy = isave[5]; lss = isave[6]; lwt = isave[7]; lwn = isave[8]; lsnd = isave[9]; lz = isave[10]; lr = isave[11]; ld = isave[12]; lt = isave[13]; lxp = isave[14]; lwa = isave[15]; mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, &wa[lws], &wa[lwy], &wa[lsy], &wa[lss], &wa[lwt], &wa[lwn], &wa[lsnd], &wa[lz], &wa[lr], &wa[ld], &wa[lt], &wa[lxp], &wa[lwa], iwa, &iwa[n], &iwa[2*n], &task[0], &task[1], lsave, &isave[21], dsave, maxls, &ln_task[0], &ln_task[1]); return; } static inline void save_local_vars( CBLAS_INT prjctd, CBLAS_INT cnstnd, CBLAS_INT boxed, CBLAS_INT updatd, CBLAS_INT nintol, CBLAS_INT iback, CBLAS_INT nskip, CBLAS_INT head, CBLAS_INT col, CBLAS_INT itail, CBLAS_INT iter, CBLAS_INT iupdat, CBLAS_INT nseg, CBLAS_INT nfgv, CBLAS_INT info, CBLAS_INT ifun, CBLAS_INT iword, CBLAS_INT nfree, CBLAS_INT nact, CBLAS_INT ileave, CBLAS_INT nenter, double theta, double fold, double tol, double dnorm, double gd, double stpmx, double sbgnrm, double stp, double gdold, double dtd, CBLAS_INT* lsave, CBLAS_INT* isave, double* dsave ) { // Save local variables lsave[0] = prjctd; lsave[1] = cnstnd; lsave[2] = boxed; lsave[3] = updatd; isave[0] = nintol; isave[3] = iback; isave[4] = nskip; isave[5] = head; isave[6] = col; isave[7] = itail; isave[8] = iter; isave[9] = iupdat; isave[11] = nseg; isave[12] = nfgv; isave[13] = info; isave[14] = ifun; isave[15] = iword; isave[16] = nfree; isave[17] = nact; isave[18] = ileave; isave[19] = nenter; dsave[0] = theta; dsave[1] = fold; dsave[2] = tol; dsave[3] = dnorm; // dsave[4] = epsmch; // dsave[5] = cpu1; // dsave[6] = cachyt; // dsave[7] = sbtime; // dsave[8] = lnscht; // dsave[9] = time1; dsave[10] = gd; dsave[11] = stpmx; dsave[12] = sbgnrm; dsave[13] = stp; dsave[14] = gdold; dsave[15] = dtd; } void mainlb(CBLAS_INT n, CBLAS_INT m, double* x, double* l, double* u, CBLAS_INT* nbd, double* f, double* g, double factr, double pgtol, double* ws, double* wy, double* sy, double* ss, double* wt, double* wn, double* snd, double* z, double* r, double* d, double* t, double* xp, double* wa, CBLAS_INT* index, CBLAS_INT* iwhere, CBLAS_INT* indx2, CBLAS_INT* task, CBLAS_INT* task_msg, CBLAS_INT* lsave, CBLAS_INT* isave, double* dsave, CBLAS_INT maxls, CBLAS_INT* temp_task, CBLAS_INT* temp_taskmsg) { // ************ // // Subroutine mainlb // // This subroutine solves bound constrained optimization problems by // using the compact formula of the limited memory BFGS updates. // // n is an integer variable. // On entry n is the number of variables. // On exit n is unchanged. // // m is an integer variable. // On entry m is the maximum number of variable metric // corrections allowed in the limited memory matrix. // On exit m is unchanged. // // x is a double precision array of dimension n. // On entry x is an approximation to the solution. // On exit x is the current approximation. // // l is a double precision array of dimension n. // On entry l is the lower bound of x. // On exit l is unchanged. // // u is a double precision array of dimension n. // On entry u is the upper bound of x. // On exit u is unchanged. // // nbd is an integer array of dimension n. // On entry nbd represents the type of bounds imposed on the // variables, and must be specified as follows: // nbd(i)=0 if x(i) is unbounded, // 1 if x(i) has only a lower bound, // 2 if x(i) has both lower and upper bounds, // 3 if x(i) has only an upper bound. // On exit nbd is unchanged. // // f is a double precision variable. // On first entry f is unspecified. // On final exit f is the value of the function at x. // // g is a double precision array of dimension n. // On first entry g is unspecified. // On final exit g is the value of the gradient at x. // // factr is a double precision variable. // On entry factr >= 0 is specified by the user. The iteration // will stop when // // (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch // // where epsmch is the machine precision, which is automatically // generated by the code. // On exit factr is unchanged. // // pgtol is a double precision variable. // On entry pgtol >= 0 is specified by the user. The iteration // will stop when // // max{|proj g_i | i = 1, ..., n} <= pgtol // // where pg_i is the ith component of the projected gradient. // On exit pgtol is unchanged. // // ws, wy, sy, and wt are double precision working arrays used to // store the following information defining the limited memory // BFGS matrix: // ws, of dimension n x m, stores S, the matrix of s-vectors; // wy, of dimension n x m, stores Y, the matrix of y-vectors; // sy, of dimension m x m, stores S'Y; // ss, of dimension m x m, stores S'S; // yy, of dimension m x m, stores Y'Y; // wt, of dimension m x m, stores the Cholesky factorization // of (theta*S'S+LD^(-1)L'); see eq. // (2.26) in [3]. // // wn is a double precision working array of dimension 2m x 2m // used to store the LEL^T factorization of the indefinite matrix // K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] // [L_a -R_z theta*S'AA'S ] // // where E = [-I 0] // [ 0 I] // // snd is a double precision working array of dimension 2m x 2m // used to store the lower triangular part of // N = [Y' ZZ'Y L_a'+R_z'] // [L_a +R_z S'AA'S ] // // z(n),r(n),d(n),t(n), xp(n),wa(8*m) are double precision working arrays. // z is used at different times to store the Cauchy point and // the Newton point. // xp is used to safeguard the projected Newton direction // // sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays. // // index is an integer working array of dimension n. // In subroutine freev, index is used to store the free and fixed // variables at the Generalized Cauchy Point (GCP). // // iwhere is an integer working array of dimension n used to record // the status of the vector x for GCP computation. // iwhere(i)=0 or -3 if x(i) is free and has bounds, // 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) // 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) // 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) // -1 if x(i) is always free, i.e., no bounds on it. // // indx2 is an integer working array of dimension n. // Within subroutine cauchy, indx2 corresponds to the array iorder. // In subroutine freev, a list of variables entering and leaving // the free set is stored in indx2, and it is passed on to // subroutine formk with this information. // // task is a working integer indicating // the current job when entering and leaving this subroutine. // // task_msg is a working integer indicating the status message. // // iprint is an INTEGER variable that must be set by the user. // It controls the frequency and type of output generated: // iprint<0 no output is generated; // iprint=0 print only one line at the last iteration; // 0100 print details of every iteration including x and g; // // // csave is a working string of characters of length 60. // // lsave is a logical working array of dimension 4. // // isave is an integer working array of dimension 23. // // dsave is a double precision working array of dimension 29. // // // Subprograms called // // L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk, // // errclb, prn1lb, prn2lb, prn3lb, active, projgr, // // freev, cmprlb, matupd, formt. // // Minpack2 Library ... timer // // Linpack Library ... dcopy, ddot. // // // References: // // [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited // memory algorithm for bound constrained optimization'', // SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. // // [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN // Subroutines for Large Scale Bound Constrained Optimization'' // Tech. Report, NAM-11, EECS Department, Northwestern University, // 1994. // // [3] R. Byrd, J. Nocedal and R. Schnabel "Representations of // Quasi-Newton Matrices and their use in Limited Memory Methods'', // Mathematical Programming 63 (1994), no. 4, pp. 129-156. // // (Postscript files of these papers are available via anonymous // ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT prjctd, cnstnd, boxed, updatd, wrk; CBLAS_INT i, k, nintol, iback, nskip, head, col, iter, itail, iupdat; CBLAS_INT nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter; double theta, fold, dr, rr, tol, xstep, sbgnrm, stpmx, ddum, dnorm, dtd; double gd, gdold, stp; CBLAS_INT one_int = 1; stpmx = 0.0; if (*task == START) { // Initialize counters and scalars col = 0; // size : Actual number of variable metric corrections head = 0; // index : Location of the first s- or y-vector in S or Y theta = 1.0; iupdat = 0; // size : Total number of BFGS updates made so far updatd = 0; // bool : Whether L-BFGS matrix is updated iback = 0; itail = 0; // index iword = 0; nact = 0; ileave = 0; // index : idx2[ileave:n] are the variables leaving the free set nenter = 0; // size : No of variables entering the free set fold = 0.0; dnorm = 0.0; gd = 0.0; sbgnrm = 0.0; stp = 0.0; gdold = 0.0; dtd = 0.0; // For operation counts iter = 0; nfgv = 0; nseg = 0; nintol = 0; nskip = 0; nfree = n; ifun = 0; // For stopping tolerance tol = factr * epsmach; // 'info' records the termination information. info = 0; // Check the input arguments for errors errclb(n, m, factr, l, u, nbd, task, task_msg, &info, &k); if (*task == ERROR) { return; } // Initialize iwhere & project x onto the feasible set active(n, l, u, nbd, x, iwhere, &prjctd, &cnstnd, &boxed); // End of the initialization } else { // Restore local variables prjctd = lsave[0]; cnstnd = lsave[1]; boxed = lsave[2]; updatd = lsave[3]; nintol = isave[0]; // = isave[1]; // = isave[2]; iback = isave[3]; nskip = isave[4]; head = isave[5]; col = isave[6]; itail = isave[7]; iter = isave[8]; iupdat = isave[9]; // = isave[10]; nseg = isave[11]; nfgv = isave[12]; info = isave[13]; ifun = isave[14]; iword = isave[15]; nfree = isave[16]; nact = isave[17]; ileave = isave[18]; nenter = isave[19]; theta = dsave[0]; fold = dsave[1]; tol = dsave[2]; dnorm = dsave[3]; // = dsave[4]; // = dsave[5]; // = dsave[6]; // = dsave[7]; // = dsave[8]; // = dsave[9]; gd = dsave[10]; stpmx = dsave[11]; sbgnrm = dsave[12]; stp = dsave[13]; gdold = dsave[14]; dtd = dsave[15]; // After returning from the driver go to the point where execution is to resume. if ((*task == FG) && (*task_msg == FG_LNSRCH)) { goto LINE666; } if (*task == NEW_X) { goto LINE777; } if ((*task == FG) && (*task_msg == FG_START)) { goto LINE111; } if (*task == STOP) { if (*task_msg == STOP_CPU) { BLAS_FUNC(dcopy)(&n, t, &one_int, x, &one_int); BLAS_FUNC(dcopy)(&n, r, &one_int, g, &one_int); *f = fold; } save_local_vars( prjctd, cnstnd, boxed, updatd, nintol, iback, nskip, head, col, itail, iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, gd, stpmx, sbgnrm, stp, gdold, dtd, lsave, isave, dsave ); return; } } // Compute f0 and g0. *task = FG; *task_msg = FG_START; // Return to the driver to calculate f and g, then reenter at 111 save_local_vars( prjctd, cnstnd, boxed, updatd, nintol, iback, nskip, head, col, itail, iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, gd, stpmx, sbgnrm, stp, gdold, dtd, lsave, isave, dsave ); return; LINE111: nfgv = 1; // Compute the infinity norm of the (-) projected gradient. projgr(n, l, u, nbd, x, g, &sbgnrm); if (sbgnrm <= pgtol) { *task = CONVERGENCE; *task_msg = CONV_GRAD; save_local_vars( prjctd, cnstnd, boxed, updatd, nintol, iback, nskip, head, col, itail, iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, gd, stpmx, sbgnrm, stp, gdold, dtd, lsave, isave, dsave ); return; } // ----------------- the beginning of the loop -------------------------- while(1) { iword = -1; if ((cnstnd) || (col == 0)) { ///////////////////////////////////////////////////// // // Compute the Generalized Cauchy Point (GCP). // ///////////////////////////////////////////////////// cauchy(n, x, l, u, nbd, g, indx2, iwhere, t, d, z, m, wy, ws, sy, wt, theta, col, head, wa, &wa[2*m], &wa[4*m], &wa[6*m], &nseg, sbgnrm, &info); if (info != 0) { // Singular triangular system detected; refresh the lbfgs memory. info = 0; col = 0; head = 0; theta = 1.0; iupdat = 0; updatd = 0; continue; } nintol += nseg; // Count the entering and leaving variables for iter > 0; // find the index set of free and active variables at the GCP. freev(n, &nfree, index, &nenter, &ileave, indx2, iwhere, &wrk, updatd, cnstnd, iter); nact = n - nfree; } else { BLAS_FUNC(dcopy)(&n, x, &one_int, z, &one_int); wrk = updatd; nseg = 0; } // If there are no free variables or B=theta*I, then skip the subspace // minimization. if ((nfree != 0) && (col != 0)) { ///////////////////////////////////////////////////// // // Subspace minimization // ///////////////////////////////////////////////////// // Form the LEL^T factorization of the indefinite // matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] // [L_a -R_z theta*S'AA'S ] // where E = [-I 0] // [ 0 I] if (wrk) { formk(n, nfree, index, nenter, ileave, indx2, iupdat, updatd, wn, snd, m, ws, wy, sy, theta, col, head, &info); } if (info != 0) { // nonpositive definiteness in Cholesky factorization; // refresh the lbfgs memory and restart the iteration. info = 0; col = 0; head = 0; theta = 1.0; iupdat = 0; updatd = 0; continue; } // compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x) from 'cauchy'). cmprlb(n, m, x, g, ws, wy, sy, wt, z, r, wa, index, theta, col, head, nfree, cnstnd, &info); if (info == 0) { // Call the direct method. subsm(n, m, nfree, index, l, u, nbd, z, r, xp, ws, wy, theta, x, g, col, head, &iword, wa, wn, &info); } if (info != 0) { // singular triangular system detected; // refresh the lbfgs memory and restart the iteration. info = 0; col = 0; head = 0; theta = 1.0; iupdat = 0; updatd = 0; continue; } } ///////////////////////////////////////////////////// // // Line search and optimality tests. // ///////////////////////////////////////////////////// // Generate the search direction d := z - x. for (i = 0; i < n; i++) { d[i] = z[i] - x[i]; } LINE666: lnsrlb(n, l, u, nbd, x, *f, &fold, &gd, &gdold, g, d, r, t, z, &stp, &dnorm, &dtd, &xstep, &stpmx, iter, &ifun, &iback, &nfgv, &info, task, task_msg, boxed, cnstnd, &isave[21], &dsave[16], temp_task, temp_taskmsg); if ((info != 0) || (iback >= maxls)) { BLAS_FUNC(dcopy)(&n, t, &one_int, x, &one_int); BLAS_FUNC(dcopy)(&n, r, &one_int, g, &one_int); *f = fold; if (col == 0) { // Abnormal termination if (info == 0) { info = -9; // Restore the actual number of f and g evaluations etc. nfgv--; ifun--; iback--; } *task = ABNORMAL; *task_msg = NO_MSG; iter++; save_local_vars( prjctd, cnstnd, boxed, updatd, nintol, iback, nskip, head, col, itail, iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, gd, stpmx, sbgnrm, stp, gdold, dtd, lsave, isave, dsave ); return; } else { // Refresh the lbfgs memory and restart the iteration. if (info == 0) { nfgv--; } info = 0; col = 0; head = 0; theta = 1.0; iupdat = 0; updatd = 0; *task = RESTART; *task_msg = NO_MSG; continue; } } else if ((*task == FG) && (*task_msg == FG_LNSRCH)) { // Return to the driver for calculating f and g; renter at 666. save_local_vars( prjctd, cnstnd, boxed, updatd, nintol, iback, nskip, head, col, itail, iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, gd, stpmx, sbgnrm, stp, gdold, dtd, lsave, isave, dsave ); return; } else { // Calculate and print out the quantities related to the new X. iter++; // Compute the infinity norm of the projected (-)gradient. projgr(n, l, u, nbd, x, g, &sbgnrm); save_local_vars( prjctd, cnstnd, boxed, updatd, nintol, iback, nskip, head, col, itail, iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, gd, stpmx, sbgnrm, stp, gdold, dtd, lsave, isave, dsave ); return; } LINE777: // Test for termination. if (sbgnrm <= pgtol) { // Terminate the algorithm. *task = CONVERGENCE; *task_msg = CONV_GRAD; save_local_vars( prjctd, cnstnd, boxed, updatd, nintol, iback, nskip, head, col, itail, iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, gd, stpmx, sbgnrm, stp, gdold, dtd, lsave, isave, dsave ); return; } ddum = fmax(fmax(fabs(fold), fabs(*f)), 1.0); if ((fold - *f) <= tol*ddum) { // Terminate the algorithm. *task = CONVERGENCE; *task_msg = CONV_F; if (iback >= 10) { info = -5; } // i.e. to issue a warning if iback > 10 in the line search. save_local_vars( prjctd, cnstnd, boxed, updatd, nintol, iback, nskip, head, col, itail, iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, gd, stpmx, sbgnrm, stp, gdold, dtd, lsave, isave, dsave ); return; } // Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. for (i = 0; i < n; i++) { r[i] = g[i] - r[i]; } // 42 rr = pow(BLAS_FUNC(dnrm2)(&n, r, &one_int), 2.0); if (stp == 1.0) { dr = gd - gdold; ddum = -gdold; } else { dr = (gd - gdold)*stp; BLAS_FUNC(dscal)(&n, &stp, d, &one_int); ddum = -gdold*stp; } if (dr <= epsmach*ddum) { // Skip the L-BFGS update. nskip += 1; updatd = 0; continue; } ///////////////////////////////////////////////////// // // Update the L-BFGS matrix. // ///////////////////////////////////////////////////// updatd = 1; iupdat += 1; // Update matrices WS and WY and form the middle matrix in B. matupd(n, m, ws, wy, sy, ss, d, r, &itail, iupdat, &col, &head, &theta, rr, dr, stp, dtd); // Form the upper half of the pds T = theta*SS + L*D^(-1)*L'; // Store T in the upper triangular of the array wt; // Cholesky factorize T to J*J' with // J' stored in the upper triangular of wt. formt(m, wt, sy, ss, col, theta, &info); if (info != 0) { // Nonpositive definiteness in Cholesky factorization; // refresh the lbfgs memory and restart the iteration. info = 0; col = 0; head = 0; theta = 1.0; iupdat = 0; updatd = 0; continue; } // Now the inverse of the middle matrix in B is // [ D^(1/2) O ] [ -D^(1/2) D^(-1/2)*L' ] // [ -L*D^(-1/2) J ] [ 0 J' ] } // -------------------- the end of the loop ----------------------------- } void active(CBLAS_INT n, double* l, double* u, CBLAS_INT* nbd, double* x, CBLAS_INT* iwhere, CBLAS_INT* prjctd, CBLAS_INT* cnstnd, CBLAS_INT* boxed) { // ************ // // Subroutine active // // This subroutine initializes iwhere and projects the initial x to // the feasible set if necessary. // // iwhere is an integer array of dimension n. // On entry iwhere is unspecified. // On exit iwhere(i)=-1 if x(i) has no bounds // 3 if l(i)=u(i) // 0 otherwise. // In cauchy, iwhere is given finer gradations. // // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT i; // Initialize nbdd, prjctd, cnstnd, and boxed. *prjctd = 0; *cnstnd = 0; *boxed = 1; // Project the initial x to the feasible set if necessary. for (i = 0; i < n; i++) { if (nbd[i] > 0) { if ((nbd[i] <= 2) && (x[i] <= l[i])) { if (x[i] < l[i]) { *prjctd = 1; x[i] = l[i]; } } else if ((nbd[i] >= 2) && (x[i] >= u[i])) { if (x[i] > u[i]) { *prjctd = 1; x[i] = u[i]; } } } } // 10 // Initialize iwhere and assign values to cnstnd and boxed. for (i = 0; i < n; i++) { if (nbd[i] != 2) { *boxed = 0; } if (nbd[i] == 0) { // This variable is always free iwhere[i] = -1; } else { // Otherwise set x(i)=mid(x(i), u(i), l(i)). *cnstnd = 1; if ((nbd[i] == 2) && (u[i] - l[i] <= 0.0)) { // This variable is always fixed. iwhere[i] = 3; } else { iwhere[i] = 0; } } } // 20 return; } void bmv(CBLAS_INT m, double* sy, double* wt, CBLAS_INT col, double* v, double* p, CBLAS_INT* info) { // ************ // // Subroutine bmv // // This subroutine computes the product of the 2m x 2m middle matrix // in the compact L-BFGS formula of B and a 2m vector v; // it returns the product in p. // // m is an integer variable. // On entry m is the maximum number of variable metric corrections // used to define the limited memory matrix. // On exit m is unchanged. // // sy is a double precision array of dimension m x m. // On entry sy specifies the matrix S'Y. // On exit sy is unchanged. // // wt is a double precision array of dimension m x m. // On entry wt specifies the upper triangular matrix J' which is // the Cholesky factor of (thetaS'S+LD^(-1)L'). // On exit wt is unchanged. // // col is an integer variable. // On entry col specifies the number of s-vectors (or y-vectors) // stored in the compact L-BFGS formula. // On exit col is unchanged. // // v is a double precision array of dimension 2col. // On entry v specifies vector v. // On exit v is unchanged. // // p is a double precision array of dimension 2col. // On entry p is unspecified. // On exit p is the product Mv. // // info is an integer variable. // On entry info is unspecified. // On exit info = 0 for normal return, // = nonzero for abnormal return when the system // to be solved by dtrsl is singular. // // Subprograms called: // // Linpack ... dtrsl. // // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT i, k, i2, one_int = 1; double ssum; char *uplo = "U", *trans = "T", *diag = "N"; if (col == 0) { return; } // PART I: solve [ D^(1/2) O ] [ p1 ] = [ v1 ] // [ -L*D^(-1/2) J ] [ p2 ] [ v2 ]. // solve Jp2=v2+LD^(-1)v1. p[col] = v[col]; for (i = 1; i < col; i++) { i2 = col + i; ssum = 0.0; for (k = 0; k <= i-1; k++) { ssum = ssum + sy[i + m*k] * v[k] / sy[k + m*k]; } // 10 p[i2] = v[i2] + ssum; } // 20 // Solve the triangular system // dtrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info) BLAS_FUNC(dtrtrs)(uplo, trans, diag, &col, &one_int, wt, &m, &p[col], &col, info); if (*info != 0) { return; } // Solve D^(1/2)p1=v1. for (i = 0; i < col; i++) { p[i] = v[i] / sqrt(sy[i + m*i]); } // 30 // PART II: solve [ -D^(1/2) D^(-1/2)*L' ] [ p1 ] = [ p1 ] // [ 0 J' ] [ p2 ] [ p2 ]. // solve J^Tp2=p2. // dtrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info) trans = "N"; BLAS_FUNC(dtrtrs)(uplo, trans, diag, &col, &one_int, wt, &m, &p[col], &col, info); if (*info != 0) { return; } // compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2) // =-D^(-1/2)p1+D^(-1)L'p2. for (i = 0; i < col; i++) { p[i] = -p[i] / sqrt(sy[i + m*i]); } // 40 for (i = 0; i < col; i++) { ssum = 0.0; for (k = i+1; k < col; k++) { ssum = ssum + sy[k + m*i] * p[col + k] / sy[i + m*i]; } // 50 p[i] = p[i] + ssum; } // 60 return; } void cauchy(CBLAS_INT n, double* x, double* l, double* u, CBLAS_INT* nbd, double *g, CBLAS_INT* iorder, CBLAS_INT* iwhere, double* t, double* d, double* xcp, CBLAS_INT m, double* wy, double* ws, double* sy, double* wt, double theta, CBLAS_INT col, CBLAS_INT head, double* p, double* c, double* wbp, double* v, CBLAS_INT* nseg, double sbgnrm, CBLAS_INT* info) { // ************ // // Subroutine cauchy // // For given x, l, u, g (with sbgnrm > 0), and a limited memory // BFGS matrix B defined in terms of matrices WY, WS, WT, and // scalars head, col, and theta, this subroutine computes the // generalized Cauchy point (GCP), defined as the first local // minimizer of the quadratic // // Q(x + s) = g's + 1/2 s'Bs // // along the projected gradient direction P(x-tg,l,u). // The routine returns the GCP in xcp. // // n is an integer variable. // On entry n is the dimension of the problem. // On exit n is unchanged. // // x is a double precision array of dimension n. // On entry x is the starting point for the GCP computation. // On exit x is unchanged. // // l is a double precision array of dimension n. // On entry l is the lower bound of x. // On exit l is unchanged. // // u is a double precision array of dimension n. // On entry u is the upper bound of x. // On exit u is unchanged. // // nbd is an integer array of dimension n. // On entry nbd represents the type of bounds imposed on the // variables, and must be specified as follows: // nbd(i)=0 if x(i) is unbounded, // 1 if x(i) has only a lower bound, // 2 if x(i) has both lower and upper bounds, and // 3 if x(i) has only an upper bound. // On exit nbd is unchanged. // // g is a double precision array of dimension n. // On entry g is the gradient of f(x). g must be a nonzero vector. // On exit g is unchanged. // // iorder is an integer working array of dimension n. // iorder will be used to store the breakpoints in the piecewise // linear path and free variables encountered. On exit, // iorder(1),...,iorder(nleft) are indices of breakpoints // which have not been encountered; // iorder(nleft+1),...,iorder(nbreak) are indices of // encountered breakpoints; and // iorder(nfree),...,iorder(n) are indices of variables which // have no bound constraits along the search direction. // // iwhere is an integer array of dimension n. // On entry iwhere indicates only the permanently fixed (iwhere=3) // or free (iwhere= -1) components of x. // On exit iwhere records the status of the current x variables. // iwhere(i)=-3 if x(i) is free and has bounds, but is not moved // 0 if x(i) is free and has bounds, and is moved // 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) // 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) // 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) // -1 if x(i) is always free, i.e., it has no bounds. // // t is a double precision working array of dimension n. // t will be used to store the break points. // // d is a double precision array of dimension n used to store // the Cauchy direction P(x-tg)-x. // // xcp is a double precision array of dimension n used to return the // GCP on exit. // // m is an integer variable. // On entry m is the maximum number of variable metric corrections // used to define the limited memory matrix. // On exit m is unchanged. // // ws, wy, sy, and wt are double precision arrays. // On entry they store information that defines the // limited memory BFGS matrix: // ws(n,m) stores S, a set of s-vectors; // wy(n,m) stores Y, a set of y-vectors; // sy(m,m) stores S'Y; // wt(m,m) stores the // Cholesky factorization of (theta*S'S+LD^(-1)L'). // On exit these arrays are unchanged. // // theta is a double precision variable. // On entry theta is the scaling factor specifying B_0 = theta I. // On exit theta is unchanged. // // col is an integer variable. // On entry col is the actual number of variable metric // corrections stored so far. // On exit col is unchanged. // // head is an integer variable. // On entry head is the location of the first s-vector (or y-vector) // in S (or Y). // On exit col is unchanged. // // p is a double precision working array of dimension 2m. // p will be used to store the vector p = W^(T)d. // // c is a double precision working array of dimension 2m. // c will be used to store the vector c = W^(T)(xcp-x). // // wbp is a double precision working array of dimension 2m. // wbp will be used to store the row of W corresponding // to a breakpoint. // // v is a double precision working array of dimension 2m. // // nseg is an integer variable. // On exit nseg records the number of quadratic segments explored // in searching for the GCP. // // sg and yg are double precision arrays of dimension m. // On entry sg and yg store S'g and Y'g correspondingly. // On exit they are unchanged. // // iprint is an INTEGER variable that must be set by the user. // It controls the frequency and type of output generated: // iprint<0 no output is generated; // iprint=0 print only one line at the last iteration; // 0100 print details of every iteration including x and g; // // // sbgnrm is a double precision variable. // On entry sbgnrm is the norm of the projected gradient at x. // On exit sbgnrm is unchanged. // // info is an integer variable. // On entry info is 0. // On exit info = 0 for normal return, // = nonzero for abnormal return when the system // used in routine bmv is singular. // // Subprograms called: // // L-BFGS-B Library ... hpsolb, bmv. // // Linpack ... dscal dcopy, daxpy. // // // References: // // [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited // memory algorithm for bound constrained optimization'', // SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. // // [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN // Subroutines for Large Scale Bound Constrained Optimization'' // Tech. Report, NAM-11, EECS Department, Northwestern University, // 1994. // // (Postscript files of these papers are available via anonymous // ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT bnded, xlower, xupper; CBLAS_INT i, j, col2, nfree, nbreak, pointr, ibp, nleft, ibkmin, iter, one_int = 1; double f1, f2, dt, dtm, tsum, dibp, zibp, dibp2, bkmin, tu, tl, wmc, wmp; double wmw, tj, tj0, neggi, f2_org; f1 = 0.0; tl = 0.0; tu = 0.0; if (sbgnrm <= 0.0) { BLAS_FUNC(dcopy)(&n, x, &one_int, xcp, &one_int); return; } bnded = 1; nfree = n; nbreak = -1; ibkmin = 0; bkmin = 0.0; col2 = 2*col; // We set p to zero and build it up as we determine d. for (i = 0; i < col2; i++) { p[i] = 0.0; } // 20 // In the following loop we determine for each variable its boud status and // its breakpoint, and update p accordingly. Smallest breakpoint is identified. for (i = 0; i < n; i++) { neggi = -g[i]; if ((iwhere[i] != 3) && (iwhere[i] != -1)) { // If x[i] is not a constant and has bounds, compute the difference // between x[i] and its bounds. if (nbd[i] <= 2) { tl = x[i] - l[i]; } if (nbd[i] >= 2) { tu = u[i] - x[i]; } // If a varaible is close enough to a bound, we treat it as at bound. xlower = ((nbd[i] <= 2) && (tl <= 0.0)); xupper = ((nbd[i] >= 2) && (tu <= 0.0)); // Reset iwhere[i] iwhere[i] = 0; if (xlower) { if (neggi <= 0.0) { iwhere[i] = 1; } } else if (xupper) { if (neggi >= 0.0) { iwhere[i] = 2; } } else { if (fabs(neggi) <= 0.0) { iwhere[i] = -3; } } } pointr = head; if ((iwhere[i] != 0) && (iwhere[i] != -1)) { d[i] = 0.0; } else { d[i] = neggi; f1 = f1 - neggi * neggi; // Calculate p := p - W'e_i* (g_i). for (j = 0; j < col; j++) { p[j] = p[j] + wy[i + n*pointr]*neggi; p[col + j] = p[col + j] + ws[i + n*pointr]*neggi; pointr = (pointr + 1) % m; } // 40 if ((nbd[i] <= 2) && (nbd[i] != 0) && (neggi < 0.0)) { // x(i) + d(i) is bounded; compute t(i). nbreak += 1; iorder[nbreak] = i; t[nbreak] = tl / (-neggi); if ((nbreak == 0) || (t[nbreak] < bkmin)) { bkmin = t[nbreak]; ibkmin = nbreak; } } else if ((nbd[i] >= 2) && (neggi > 0.0)) { // x(i) + d(i) is bounded; compute t(i). nbreak += 1; iorder[nbreak] = i; t[nbreak] = tu / neggi; if ((nbreak == 0) || (t[nbreak] < bkmin)) { bkmin = t[nbreak]; ibkmin = nbreak; } } else { // x(i) + d(i) is not bounded. nfree -= 1; iorder[nfree] = i; if (fabs(neggi) > 0.0) { bnded = 0; } } } } // 50 // The indices of the nonzero components of d are now stored // in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n). // The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin. if (theta != 1.0) { // Complete the initialization of p for theta not= one. BLAS_FUNC(dscal)(&col, &theta, &p[col], &one_int); } // Initialize GCP xcp = x. BLAS_FUNC(dcopy)(&n, x, &one_int, xcp, &one_int); if ((nbreak == -1) && (nfree == n)) { // is a zero vector, return with the initial xcp as GCP. return; } // Initialize c = W'(xcp - x) = 0. for (j = 0; j < col2; j++) { c[j] = 0.0; } // 60 // Initialize derivative f2. f2 = -(theta) * f1; f2_org = f2; if (col > 0) { bmv(m, sy, wt, col, p, v, info); if (*info != 0) { return; } f2 = f2 - BLAS_FUNC(ddot)(&col2, v, &one_int, p, &one_int); } dtm = -f1 / f2; tsum = 0.0; *nseg = 1; nleft = nbreak; iter = 0; tj = 0.0; // ------------------- the beginning of the loop ------------------------- // If there are no breakpoints, we skip the loop, locate the GCP and return. while (nbreak > -1) { // Find the next smallest breakpoint; // compute dt = t(nleft) - t(nleft + 1). tj0 = tj; if (iter == 0) { // Since we already have the smallest breakpoint we need not do // heapsort yet. Often only one breakpoint is used and the // cost of heapsort is avoided. tj = bkmin; ibp = iorder[ibkmin]; } else { if (iter == 1) { // Replace the already used smallest breakpoint with the // breakpoint numbered nbreak > nlast, before heapsort call. if (ibkmin != nbreak) { t[ibkmin] = t[nbreak]; iorder[ibkmin] = iorder[nbreak]; } } // Update heap structure of breakpoints (if iter=2, initialize heap). hpsolb(nleft, t, iorder, iter - 1); tj = t[nleft]; ibp = iorder[nleft]; } dt = tj - tj0; // If a minimizer is within this interval, locate the GCP and return. if (dtm < dt) { if (dtm <= 0.0) { dtm = 0.0; } tsum = tsum + dtm; // Move free variables (i.e. the ones w/o breakpoints) and the variables // whose breakpoints haven't been reached. BLAS_FUNC(daxpy)(&n, &tsum, d, &one_int, xcp, &one_int); // Update c = c + dtm*p = W'(x^c - x) // which will be used in computing r = Z'(B(x^c - x) + g). if (col > 0) { BLAS_FUNC(daxpy)(&col2, &dtm, p, &one_int, c, &one_int); } return; } // Otherwise fix one variable and reset the corresponding component // of d to zero. tsum = tsum + dt; nleft -= 1; iter += 1; dibp = d[ibp]; d[ibp] = 0.0; if (dibp > 0.0) { zibp = u[ibp] - x[ibp]; xcp[ibp] = u[ibp]; iwhere[ibp] = 2; } else { zibp = l[ibp] - x[ibp]; xcp[ibp] = l[ibp]; iwhere[ibp] = 1; } if (nleft == -1 && nbreak == n) { // All n variables are fixed, return with xcp as GCP. dtm = dt; // Update c = c + dtm*p = W'(x^c - x) // which will be used in computing r = Z'(B(x^c - x) + g). if (col > 0) { BLAS_FUNC(daxpy)(&col2, &dtm, p, &one_int, c, &one_int); } return; } // Update the derivative information. *nseg += 1; dibp2 = pow(dibp, 2.0); // Update f1 and f2 // Temporarily set f1 and f2 for col=0. f1 = f1 + dt * f2 + dibp2 - theta*dibp*zibp; f2 = f2 - theta*dibp2; if (col > 0) { // Update c = c + dt*p. BLAS_FUNC(daxpy)(&col2, &dt, p, &one_int, c, &one_int); // Choose wbp, the row of W corresponding to the breakpoint encountered. pointr = head; for (j = 0; j < col; j++) { wbp[j] = wy[ibp + n*pointr]; wbp[col + j] = theta * ws[ibp + n*pointr]; pointr = (pointr + 1) % m; } // 70 // Compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'. bmv(m, sy, wt, col, wbp, v, info); if (*info != 0) { return; } wmc = BLAS_FUNC(ddot)(&col2, c, &one_int, v, &one_int); wmp = BLAS_FUNC(ddot)(&col2, p, &one_int, v, &one_int); wmw = BLAS_FUNC(ddot)(&col2, wbp, &one_int, v, &one_int); // Update p = p - dibp*wbp. dibp = -dibp; BLAS_FUNC(daxpy)(&col2, &dibp, wbp, &one_int, p, &one_int); dibp = -dibp; // Complete updating f1 and f2 while col > 0. f1 = f1 + dibp*wmc; f2 = f2 + dibp*2.0*wmp - dibp2*wmw; } f2 = fmax(epsmach*f2_org, f2); if (nleft >= 0) { // to repeat the loop for unsearched intervals. dtm = -f1 / f2; } else if (bnded) { dtm = 0.0; break; } else { dtm = -f1 / f2; break; } } // ------------------- the end of the loop ------------------------------- if (dtm <= 0.0) { dtm = 0.0; } tsum = tsum + dtm; // Move free variables (i.e. the ones w/o breakpoints) and the variables // whose breakpoints haven't been reached. BLAS_FUNC(daxpy)(&n, &tsum, d, &one_int, xcp, &one_int); // Update c = c + dtm*p = W'(x^c - x) // which will be used in computing r = Z'(B(x^c - x) + g). if (col > 0) { BLAS_FUNC(daxpy)(&col2, &dtm, p, &one_int, c, &one_int); } return; } void cmprlb(CBLAS_INT n, CBLAS_INT m, double* x, double* g, double* ws, double* wy, double* sy, double* wt, double* z, double* r, double* wa, CBLAS_INT* index, double theta, CBLAS_INT col, CBLAS_INT head, CBLAS_INT nfree, CBLAS_INT cnstnd, CBLAS_INT* info) { // ************ // // Subroutine cmprlb // // This subroutine computes r=-Z'B(xcp-xk)-Z'g by using // wa(2m+1)=W'(xcp-x) from subroutine cauchy. // // Subprograms called: // // L-BFGS-B Library ... bmv. // // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT i, j, k, pointr; double a1, a2; if ((!cnstnd) && (col > 0)) { for (i = 0; i < n; i++) { r[i] = -g[i]; } // 26 } else { for (i = 0; i < nfree; i++) { k = index[i]; r[i] = -theta*(z[k] - x[k]) - g[k]; } // 30 bmv(m, sy, wt, col, &wa[2*m], wa, info); if (*info != 0) { *info = -8; return; } pointr = head; for (j = 0; j < col; j++) { a1 = wa[j]; a2 = theta*wa[col + j]; for (i = 0; i < nfree; i++) { k = index[i]; r[i] = r[i] + wy[k + n*pointr]*a1 + ws[k + n*pointr]*a2; } // 32 pointr = (pointr + 1) % m; } // 34 } return; } void errclb(CBLAS_INT n, CBLAS_INT m, double factr, double* l, double* u, CBLAS_INT* nbd, CBLAS_INT* task, CBLAS_INT* task_msg, CBLAS_INT* info, CBLAS_INT* k) { // ************ // // Subroutine errclb // // This subroutine checks the validity of the input data. // // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT i; // Check the input arguments for errors. if (n <= 0) { *task = ERROR; *task_msg = ERROR_N; } if (m <= 0) { *task = ERROR; *task_msg = ERROR_M; } if (factr < 0) { *task = ERROR; *task_msg = ERROR_FACTR; } // Check the validity of the arrays nbd(i), u(i), and l(i). for (i = 0; i < n; i++) { if ((nbd[i] < 0) || (nbd[i] > 3)) { *task = ERROR; *task_msg = ERROR_NBD; *info = -6; *k = i; } if (nbd[i] == 2) { if (l[i] > u[i]) { // Return *task = ERROR; *task_msg = ERROR_NOFEAS; *info = -7; *k = i; } } } return; } void formk(CBLAS_INT n, CBLAS_INT nsub, CBLAS_INT* ind, CBLAS_INT nenter, CBLAS_INT ileave, CBLAS_INT* indx2, CBLAS_INT iupdat, CBLAS_INT updatd, double* wn, double* wn1, CBLAS_INT m, double* ws, double* wy, double* sy, double theta, CBLAS_INT col, CBLAS_INT head, CBLAS_INT* info) { // ************ // // Subroutine formk // // This subroutine forms the LEL^T factorization of the indefinite // // matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] // [L_a -R_z theta*S'AA'S ] // where E = [-I 0] // [ 0 I] // The matrix K can be shown to be equal to the matrix M^[-1]N // occurring in section 5.1 of [1], as well as to the matrix // Mbar^[-1] Nbar in section 5.3. // // n is an integer variable. // On entry n is the dimension of the problem. // On exit n is unchanged. // // nsub is an integer variable // On entry nsub is the number of subspace variables in free set. // On exit nsub is not changed. // // ind is an integer array of dimension nsub. // On entry ind specifies the indices of subspace variables. // On exit ind is unchanged. // // nenter is an integer variable. // On entry nenter is the number of variables entering the // free set. // On exit nenter is unchanged. // // ileave is an integer variable. // On entry indx2(ileave),...,indx2(n) are the variables leaving // the free set. // On exit ileave is unchanged. // // indx2 is an integer array of dimension n. // On entry indx2(1),...,indx2(nenter) are the variables entering // the free set, while indx2(ileave),...,indx2(n) are the // variables leaving the free set. // On exit indx2 is unchanged. // // iupdat is an integer variable. // On entry iupdat is the total number of BFGS updates made so far. // On exit iupdat is unchanged. // // updatd is a logical variable. // On entry 'updatd' is true if the L-BFGS matrix is updatd. // On exit 'updatd' is unchanged. // // wn is a double precision array of dimension 2m x 2m. // On entry wn is unspecified. // On exit the upper triangle of wn stores the LEL^T factorization // of the 2*col x 2*col indefinite matrix // [-D -Y'ZZ'Y/theta L_a'-R_z' ] // [L_a -R_z theta*S'AA'S ] // // wn1 is a double precision array of dimension 2m x 2m. // On entry wn1 stores the lower triangular part of // [Y' ZZ'Y L_a'+R_z'] // [L_a+R_z S'AA'S ] // in the previous iteration. // On exit wn1 stores the corresponding updated matrices. // The purpose of wn1 is just to store these inner products // so they can be easily updated and inserted into wn. // // m is an integer variable. // On entry m is the maximum number of variable metric corrections // used to define the limited memory matrix. // On exit m is unchanged. // // ws, wy, sy, and wtyy are double precision arrays; // theta is a double precision variable; // col is an integer variable; // head is an integer variable. // On entry they store the information defining the // limited memory BFGS matrix: // ws(n,m) stores S, a set of s-vectors; // wy(n,m) stores Y, a set of y-vectors; // sy(m,m) stores S'Y; // wtyy(m,m) stores the Cholesky factorization // of (theta*S'S+LD^(-1)L') // theta is the scaling factor specifying B_0 = theta I; // col is the number of variable metric corrections stored; // head is the location of the 1st s- (or y-) vector in S (or Y). // On exit they are unchanged. // // info is an integer variable. // On entry info is unspecified. // On exit info = 0 for normal return; // = -1 when the 1st Cholesky factorization failed; // = -2 when the 2st Cholesky factorization failed. // // Subprograms called: // // Linpack ... dcopy, dpofa, dtrsl. // // // References: // [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited // memory algorithm for bound constrained optimization'', // SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. // // [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a // limited memory FORTRAN code for solving bound constrained // optimization problems'', Tech. Report, NAM-11, EECS Department, // Northwestern University, 1994. // // (Postscript files of these papers are available via anonymous // ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT m2, ipntr, jpntr, iy, is, jy, js, is1, js1, k1, i, k, col2, pbegin, pend; CBLAS_INT dbegin, dend, upcl, one_int = 1, temp_int; double temp1, temp2, temp3, temp4; char *uplo = "U", *trans = "T", *diag = "N"; // Form the lower triangular part of // WN1 = [Y' ZZ'Y L_a'+R_z'] // [L_a+R_z S'AA'S ] // where L_a is the strictly lower triangular part of S'AA'Y // R_z is the upper triangular part of S'ZZ'Y. if (updatd) { if (iupdat > m) { // Shift old part of WN1. for (jy = 0; jy < m - 1; jy++) { js = m + jy; temp_int = m - (jy + 1); BLAS_FUNC(dcopy)(&temp_int, &wn1[(jy + 1) + 2*m*(jy + 1)], &one_int, &wn1[jy + 2*m*jy], &one_int); BLAS_FUNC(dcopy)(&temp_int, &wn1[(js + 1) + 2*m*(js + 1)], &one_int, &wn1[js + 2*m*js], &one_int); temp_int = m - 1; BLAS_FUNC(dcopy)(&temp_int, &wn1[(m + 1) + 2*m*(jy + 1)], &one_int, &wn1[m + 2*m*jy], &one_int); } // 10 } // Put new rows in blocks (1,1), (2,1) and (2,2). pbegin = 0; // Inclusive loop start pend = nsub; // Exclusive loop end dbegin = nsub; // Inclusive loop start dend = n; // Exclusive loop end iy = col - 1; is = m + col - 1; ipntr = head + col - 1; if (ipntr >= m) { ipntr -= m; } jpntr = head; for (jy = 0; jy < col; jy++) { js = m + jy; temp1 = 0.0; temp2 = 0.0; temp3 = 0.0; // Compute element jy of row 'col' of Y'ZZ'Y. for (k = pbegin; k < pend; k++) { k1 = ind[k]; temp1 = temp1 + wy[k1 + n*ipntr]*wy[k1 + n*jpntr]; } // 15 // Compute elements jy of row 'col' of L_a and S'AA'S. for (k = dbegin; k < dend; k++) { k1 = ind[k]; temp2 = temp2 + ws[k1 + n*ipntr]*ws[k1 + n*jpntr]; temp3 = temp3 + ws[k1 + n*ipntr]*wy[k1 + n*jpntr]; } // 16 wn1[iy + 2*m*jy] = temp1; wn1[is + 2*m*js] = temp2; wn1[is + 2*m*jy] = temp3; jpntr = (jpntr + 1) % m; } // 20 // Put new column in block (2,1) jy = col - 1; jpntr = head + col - 1; if (jpntr >= m) { jpntr -= m; } ipntr = head; for (i = 0; i < col; i++) { is = m + i; temp3 = 0.0; // Compute element i of column 'col' of R_z. for (k = pbegin; k < pend; k++) { k1 = ind[k]; temp3 = temp3 + ws[k1 + n*ipntr]*wy[k1 + n*jpntr]; } // 25 ipntr = (ipntr + 1) % m; wn1[is + 2*m*jy] = temp3; } // 30 upcl = col - 1; } else { upcl = col; } // modify the old parts in blocks (1,1) and (2,2) due to changes // in the set of free variables. ipntr = head; for (iy = 0; iy < upcl; iy++) { is = m + iy; jpntr = head; for (jy = 0; jy <= iy; jy++) { js = m + jy; temp1 = 0.0; temp2 = 0.0; temp3 = 0.0; temp4 = 0.0; for (k = 0; k < nenter; k++) { k1 = indx2[k]; temp1 = temp1 + wy[k1 + n*ipntr]*wy[k1 + n*jpntr]; temp2 = temp2 + ws[k1 + n*ipntr]*ws[k1 + n*jpntr]; } // 35 for (k = ileave; k < n; k++) { k1 = indx2[k]; temp3 = temp3 + wy[k1 + n*ipntr]*wy[k1 + n*jpntr]; temp4 = temp4 + ws[k1 + n*ipntr]*ws[k1 + n*jpntr]; } // 36 wn1[iy + 2*m*jy] = wn1[iy + 2*m*jy] + temp1 - temp3; wn1[is + 2*m*js] = wn1[is + 2*m*js] - temp2 + temp4; jpntr = (jpntr + 1) % m; } // 40 ipntr = (ipntr + 1) % m; } // 45 // Modify the old parts in block (2,1) ipntr = head; for (is = m; is < m + upcl; is++) { jpntr = head; for (jy = 0; jy < upcl; jy++) { temp1 = 0.0; temp3 = 0.0; for (k = 0; k < nenter; k++) { k1 = indx2[k]; temp1 = temp1 + ws[k1 + n*ipntr]*wy[k1 + n*jpntr]; } // 50 for (k = ileave; k < n; k++) { k1 = indx2[k]; temp3 = temp3 + ws[k1 + n*ipntr]*wy[k1 + n*jpntr]; } // 51 if (is <= jy + m) { wn1[is + 2*m*jy] = wn1[is + 2*m*jy] + temp1 - temp3; } else { wn1[is + 2*m*jy] = wn1[is + 2*m*jy] - temp1 + temp3; } jpntr = (jpntr + 1) % m; } // 55 ipntr = (ipntr + 1) % m; } // 60 // Form the upper triangle of WN = [D+Y' ZZ'Y/theta -L_a'+R_z' ] // [-L_a +R_z S'AA'S*theta] m2 = 2*m; for (iy = 0; iy < col; iy++) { is = col + iy; is1 = m + iy; for (jy = 0; jy <= iy; jy++) { js = col + jy; js1 = m + jy; wn[jy + m2*iy] = wn1[iy + m2*jy] / theta; wn[js + m2*is] = wn1[is1 + m2*js1] * theta; } // 65 for (jy = 0; jy < iy; jy++) { wn[jy + m2*is] = -wn1[is1 + m2*jy]; } // 66 for (jy = iy; jy < col; jy++) { wn[jy + m2*is] = wn1[is1 + m2*jy]; } // 67 wn[iy + m2*iy] = wn[iy + m2*iy] + sy[iy + m*iy]; } // 70 // Form the upper triangle of WN= [ LL' L^-1(-L_a'+R_z')] // [(-L_a +R_z)L'^-1 S'AA'S*theta ] // // first Cholesky factor (1,1) block of wn to get LL' // with L' stored in the upper triangle of wn. // dpotrf(uplo, n, a, lda, info) BLAS_FUNC(dpotrf)(uplo, &col, wn, &m2, info); if (*info != 0) { *info = -1; return; } // Then form L^-1(-L_a'+R_z') in the (1,2) block. col2 = 2*col; // dtrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info) BLAS_FUNC(dtrtrs)(uplo, trans, diag, &col, &col, wn, &m2, &wn[m2*col], &m2, info); // Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the upper // triangle of (2,2) block of wn. for (is = col; is < col2; is++) { for (js = is; js < col2; js++) { wn[is + m2*js] = wn[is + m2*js] + BLAS_FUNC(ddot)(&col, &wn[m2*is], &one_int, &wn[m2*js], &one_int); } // 74 } // 72 // dpotrf(uplo, n, a, lda, info) BLAS_FUNC(dpotrf)(uplo, &col, &wn[col + m2*col], &m2, info); if (*info != 0) { *info = -2; } return; } void formt(CBLAS_INT m, double* wt, double* sy, double* ss, CBLAS_INT col, double theta, CBLAS_INT* info) { // ************ // // Subroutine formt // // This subroutine forms the upper half of the pos. def. and symm. // T = theta*SS + L*D^(-1)*L', stores T in the upper triangle // of the array wt, and performs the Cholesky factorization of T // to produce J*J', with J' stored in the upper triangle of wt. // // Subprograms called: // // Linpack ... dpofa. // // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT i, j, k, k1; double ddum; char *uplo = "U"; // Form the upper half of T = theta*SS + L*D^(-1)*L', store T in the upper // triangle of the array wt. for (j = 0; j < col; j++) { wt[m*j] = theta*ss[m*j]; } // 52 for (i = 1; i < col; i++) { for (j = i; j < col; j++) { k1 = (i > j ? j : i); ddum = 0.0; for (k = 0; k < k1; k++) { ddum = ddum + ((sy[i + m*k] * sy[j + m*k]) / sy[k + m*k]); } // 53 wt[i + m*j] = ddum + theta*ss[i + m*j]; } // 54 } // 55 // Cholesky factorize T to J*J' with J' stored in the upper triangle of wt. // dpotrf(uplo, n, a, lda, info) BLAS_FUNC(dpotrf)(uplo, &col, wt, &m, info); if (*info != 0) { *info = -3; } return; } void freev(CBLAS_INT n, CBLAS_INT* nfree, CBLAS_INT* idx, CBLAS_INT* nenter, CBLAS_INT* ileave, CBLAS_INT* idx2, CBLAS_INT* iwhere, CBLAS_INT* wrk, CBLAS_INT updatd, CBLAS_INT cnstnd, CBLAS_INT iter) { // ************ // // Subroutine freev // // This subroutine counts the entering and leaving variables when // iter > 0, and finds the index set of free and active variables // at the GCP. // // cnstnd is a logical variable indicating whether bounds are present // // index is an integer array of dimension n // for i=1,...,nfree, index(i) are the indices of free variables // for i=nfree+1,...,n, index(i) are the indices of bound variables // On entry after the first iteration, index gives // the free variables at the previous iteration. // On exit it gives the free variables based on the determination // in cauchy using the array iwhere. // // indx2 is an integer array of dimension n // On entry indx2 is unspecified. // On exit with iter>0, indx2 indicates which variables // have changed status since the previous iteration. // For i= 1,...,nenter, indx2(i) have changed from bound to free. // For i= ileave+1,...,n, indx2(i) have changed from free to bound. // // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT i, iact, k; // in formk // nenter will be used as exclusive upper bound // ileave will be used as the inclusive lower bound *nenter = 0; *ileave = n; if ((iter > 0) && cnstnd) { // Count the entering and leaving variables. for (i = 0; i < *nfree; i++) { k = idx[i]; if (iwhere[k] > 0) { *ileave -= 1; idx2[*ileave] = k; } } // 20 for (i = *nfree; i < n; i++) { k = idx[i]; if (iwhere[k] <= 0) { idx2[*nenter] = k; *nenter += 1; } } // 22 } *wrk = ((*ileave < n) || (*nenter > 0) || updatd); // Find the index set of free and active variables at the GCP. *nfree = 0; iact = n; for (i = 0; i < n; i++) { if (iwhere[i] <= 0) { idx[*nfree] = i; *nfree += 1; } else { iact -= 1; idx[iact] = i; } } // 24 return; } void hpsolb(CBLAS_INT n, double* t, CBLAS_INT* iorder, CBLAS_INT iheap) { // ************ // // Subroutine hpsolb // // This subroutine sorts out the least element of t, and puts the // remaining elements of t in a heap. // // n is an integer variable. // On entry n is the dimension of the arrays t and iorder. // On exit n is unchanged. // // t is a double precision array of dimension n. // On entry t stores the elements to be sorted, // On exit t(n) stores the least elements of t, and t(1) to t(n-1) // stores the remaining elements in the form of a heap. // // iorder is an integer array of dimension n. // On entry iorder(i) is the index of t(i). // On exit iorder(i) is still the index of t(i), but iorder may be // permuted in accordance with t. // // iheap is an integer variable specifying the task. // On entry iheap should be set as follows: // iheap .eq. 0 if t(1) to t(n) is not in the form of a heap, // iheap .ne. 0 if otherwise. // On exit iheap is unchanged. // // // References: // Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT. // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // ************ CBLAS_INT i, j, k, indxin, indxout; double ddum, out; if (iheap == 0) { // Rearrange the elements t(1) to t(n) to form a heap. for (k = 2; k <= n + 1; k++) { ddum = t[k-1]; indxin = iorder[k-1]; // Add ddum to the heap i = k; while (1) { if (i > 1) { j = i / 2; if (ddum < t[j-1]) { t[i-1] = t[j-1]; iorder[i-1] = iorder[j-1]; i = j; } else { break; } } else { break; } } // 10 t[i-1] = ddum; iorder[i-1] = indxin; } // 20 } // Assign to 'out' the value of t(1), the least member of the heap, and // rearrange the remaining members to form a heap as elements 1 to n-1 of t. if (n > 0) { i = 1; out = t[0]; indxout = iorder[0]; ddum = t[n]; indxin = iorder[n]; // Restore the heap while (1) { j = i + i; if (j <= n) { if (t[j] < t[j-1]) { j++; } if (t[j-1] < ddum) { t[i-1] = t[j-1]; iorder[i-1] = iorder[j-1]; i = j; } else { break; } } else { break; } } t[i-1] = ddum; iorder[i-1] = indxin; // Put the least member in t(n). t[n] = out; iorder[n] = indxout; } return; } void lnsrlb(CBLAS_INT n, double* l, double* u, CBLAS_INT* nbd, double* x, double f, double* fold, double *gd, double *gdold, double* g, double* d, double* r, double* t, double* z, double* stp, double* dnorm, double* dtd, double* xstep, double* stpmx, CBLAS_INT iter, CBLAS_INT* ifun, CBLAS_INT* iback, CBLAS_INT* nfgv, CBLAS_INT* info, CBLAS_INT* task, CBLAS_INT* task_msg, CBLAS_INT boxed, CBLAS_INT cnstnd, CBLAS_INT* isave, double* dsave, CBLAS_INT* temp_task, CBLAS_INT* temp_taskmsg) { // ********** // // Subroutine lnsrlb // // This subroutine calls subroutine dcsrch from the Minpack2 library // to perform the line search. Subroutine dscrch is safeguarded so // that all trial points lie within the feasible region. // // Subprograms called: // // Minpack2 Library ... dcsrch. // // Linpack ... dtrsl, ddot. // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ********** CBLAS_INT i, one_int = 1; double a1, a2, ftol = 1e-3, gtol = 0.9, xtol = 0.1; if (*task_msg != FG_LNSRCH) { *dnorm = BLAS_FUNC(dnrm2)(&n, d, &one_int); *dtd = pow(*dnorm, 2.0); // Determine the maximum step length. *stpmx = 1.0e+10; if (cnstnd) { if (iter == 0) { *stpmx = 1.0; } else { for (i = 0; i < n; i++) { a1 = d[i]; if (nbd[i] != 0) { if ((a1 < 0.0) && (nbd[i] <= 2)) { a2 = l[i] - x[i]; if (a2 >= 0.0) { *stpmx = 0.0; } else if (a1*(*stpmx) < a2) { *stpmx = a2 / a1; } } else if ((a1 > 0.0) && (nbd[i] >= 2)) { a2 = u[i] - x[i]; if (a2 <= 0.0) { *stpmx = 0.0; } else if (a1*(*stpmx) > a2) { *stpmx = a2 / a1; } } } } // 43 } } if ((iter == 0) && (!(boxed))) { *stp = fmin(1.0 / (*dnorm), *stpmx); } else { *stp = 1.0; } for (i = 0; i < n; i++) { t[i] = x[i]; r[i] = g[i]; } *fold = f; // Later used in mainlb, see control flow after returning from this function. *ifun = 0; *iback = 0; *temp_task = START; *temp_taskmsg = NO_MSG; } *gd = BLAS_FUNC(ddot)(&n, g, &one_int, d, &one_int); if (*ifun == 0) { *gdold = *gd; if (*gd >= 0.0) { // The directional derivative >= 0, line search is impossible. *info = -4; return; } } dcsrch(f, *gd, stp, ftol, gtol, xtol, 0.0, *stpmx, temp_task, temp_taskmsg, isave, dsave); *xstep = (*stp) * (*dnorm); if ((*temp_task != CONVERGENCE) && (*temp_task != WARNING)) { *task = FG; *task_msg = FG_LNSRCH; *ifun += 1; *nfgv += 1; *iback = *ifun - 1; if (*stp == 1.0) { // for (i = 0; i < n; i++) { x[i] = z[i]; } BLAS_FUNC(dcopy)(&n, z, &one_int, x, &one_int); } else { for (i = 0; i < n; i++) { // Take step and prevent rounding error beyond bound. x[i] = (*stp) * d[i] + t[i]; if ((nbd[i] == 1) || (nbd[i] == 2)) { x[i] = fmax(x[i], l[i]); } if ((nbd[i] == 2) || (nbd[i] == 3)) { x[i] = fmin(x[i], u[i]); } } // 41 } } else { *task = NEW_X; *task_msg = NO_MSG; } return; } void matupd(CBLAS_INT n, CBLAS_INT m, double* ws, double *wy, double* sy, double* ss, double* d, double* r, CBLAS_INT* itail, CBLAS_INT iupdat, CBLAS_INT* col, CBLAS_INT* head, double* theta, double rr, double dr, double stp, double dtd) { // ************ // // Subroutine matupd // // This subroutine updates matrices WS and WY, and forms the // middle matrix in B. // // Subprograms called: // // Linpack ... dcopy, ddot. // // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT j, pointr, temp_int, one_int = 1; // Set pointers for matrices WS and WY. if (iupdat <= m) { *col = iupdat; *itail = (*head + iupdat - 1) % m; } else { *itail = (*itail + 1) % m; *head = (*head + 1) % m; } // Update matrices WS and WY. BLAS_FUNC(dcopy)(&n, d, &one_int, &ws[(*itail) * n], &one_int); BLAS_FUNC(dcopy)(&n, r, &one_int, &wy[(*itail) * n], &one_int); // Set theta=yy/ys. *theta = rr / dr; // Form the middle matrix in B. // update the upper triangle of SS (m, m), and the lower triangle of SY (m, m): if (iupdat > m) { for (j = 1; j < *col; j++) { BLAS_FUNC(dcopy)(&j, &ss[1 + m*j], &one_int, &ss[m*(j-1)], &one_int); temp_int = *col - j; BLAS_FUNC(dcopy)(&temp_int, &sy[j + m*j], &one_int, &sy[(j-1) + m*(j-1)], &one_int); } // 50 } // add new information: the last row of SY and the last column of SS pointr = *head; for (j = 0; j < *col - 1; j++) { sy[*col - 1 + m*j] = BLAS_FUNC(ddot)(&n, d, &one_int, &wy[pointr*n], &one_int); ss[j + m*(*col - 1)] = BLAS_FUNC(ddot)(&n, &ws[pointr*n], &one_int, &d[0], &one_int); pointr = (pointr + 1) % m; } // 51 ss[(*col - 1) + m*(*col - 1)] = (stp == 1.0 ? dtd : stp*stp*dtd ); sy[(*col - 1) + m*(*col - 1)] = dr; return; } void projgr(CBLAS_INT n, double* l, double* u, CBLAS_INT* nbd, double* x, double* g, double* sbgnrm) { // ************ // // Subroutine projgr // // This subroutine computes the infinity norm of the projected // gradient. // // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT i; double gi; *sbgnrm = 0.0; for (i = 0; i < n; i++) { gi = g[i]; if (gi != gi) { // NaN value in gradient: propagate it *sbgnrm = gi; return; } if (nbd[i] != 0) { if (gi < 0.0) { if (nbd[i] >= 2) { gi = fmax(x[i] - u[i], gi); } } else { if (nbd[i] <= 2) { gi = fmin(x[i] - l[i], gi); } } } *sbgnrm = fmax(*sbgnrm, fabs(gi)); } // 15 return; } void subsm(CBLAS_INT n, CBLAS_INT m, CBLAS_INT nsub, CBLAS_INT* ind, double* l, double* u, CBLAS_INT* nbd, double* x, double* d, double* xp, double* ws, double* wy, double theta, double* xx, double* gg, CBLAS_INT col, CBLAS_INT head, CBLAS_INT* iword, double* wv, double* wn, CBLAS_INT* info) { // ********************************************************************** // // This routine contains the major changes in the updated version. // The changes are described in the accompanying paper // // Jose Luis Morales, Jorge Nocedal // "Remark On Algorithm 788: L-BFGS-B: Fortran Subroutines for Large-Scale // Bound Constrained Optimization". Decemmber 27, 2010. // // J.L. Morales Departamento de Matematicas, // Instituto Tecnologico Autonomo de Mexico // Mexico D.F. // // J, Nocedal Department of Electrical Engineering and // Computer Science. // Northwestern University. Evanston, IL. USA // // January 17, 2011 // // ********************************************************************** // // // Subroutine subsm // // Given xcp, l, u, r, an index set that specifies // the active set at xcp, and an l-BFGS matrix B // (in terms of WY, WS, SY, WT, head, col, and theta), // this subroutine computes an approximate solution // of the subspace problem // // (P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp) // // subject to l<=x<=u // x_i=xcp_i for all i in A(xcp) // // along the subspace unconstrained Newton direction // // d = -(Z'BZ)^(-1) r. // // The formula for the Newton direction, given the L-BFGS matrix // and the Sherman-Morrison formula, is // // d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r. // // where // K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] // [L_a -R_z theta*S'AA'S ] // // Note that this procedure for computing d differs // from that described in [1]. One can show that the matrix K is // equal to the matrix M^[-1]N in that paper. // // n is an integer variable. // On entry n is the dimension of the problem. // On exit n is unchanged. // // m is an integer variable. // On entry m is the maximum number of variable metric corrections // used to define the limited memory matrix. // On exit m is unchanged. // // nsub is an integer variable. // On entry nsub is the number of free variables. // On exit nsub is unchanged. // // ind is an integer array of dimension nsub. // On entry ind specifies the coordinate indices of free variables. // On exit ind is unchanged. // // l is a double precision array of dimension n. // On entry l is the lower bound of x. // On exit l is unchanged. // // u is a double precision array of dimension n. // On entry u is the upper bound of x. // On exit u is unchanged. // // nbd is a integer array of dimension n. // On entry nbd represents the type of bounds imposed on the // variables, and must be specified as follows: // nbd(i)=0 if x(i) is unbounded, // 1 if x(i) has only a lower bound, // 2 if x(i) has both lower and upper bounds, and // 3 if x(i) has only an upper bound. // On exit nbd is unchanged. // // x is a double precision array of dimension n. // On entry x specifies the Cauchy point xcp. // On exit x(i) is the minimizer of Q over the subspace of // free variables. // // d is a double precision array of dimension n. // On entry d is the reduced gradient of Q at xcp. // On exit d is the Newton direction of Q. // // xp is a double precision array of dimension n. // used to safeguard the projected Newton direction // // xx is a double precision array of dimension n // On entry it holds the current iterate // On output it is unchanged // gg is a double precision array of dimension n // On entry it holds the gradient at the current iterate // On output it is unchanged // // ws and wy are double precision arrays; // theta is a double precision variable; // col is an integer variable; // head is an integer variable. // On entry they store the information defining the // limited memory BFGS matrix: // ws(n,m) stores S, a set of s-vectors; // wy(n,m) stores Y, a set of y-vectors; // theta is the scaling factor specifying B_0 = theta I; // col is the number of variable metric corrections stored; // head is the location of the 1st s- (or y-) vector in S (or Y). // On exit they are unchanged. // // iword is an integer variable. // On entry iword is unspecified. // On exit iword specifies the status of the subspace solution. // iword = 0 if the solution is in the box, // 1 if some bound is encountered. // // wv is a double precision working array of dimension 2m. // // wn is a double precision array of dimension 2m x 2m. // On entry the upper triangle of wn stores the LEL^T factorization // of the indefinite matrix // // K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] // [L_a -R_z theta*S'AA'S ] // where E = [-I 0] // [ 0 I] // On exit wn is unchanged. // // iprint is an INTEGER variable that must be set by the user. // It controls the frequency and type of output generated: // iprint<0 no output is generated; // iprint=0 print only one line at the last iteration; // 0100 print details of every iteration including x and g; // // // info is an integer variable. // On entry info is unspecified. // On exit info = 0 for normal return, // = nonzero for abnormal return // when the matrix K is ill-conditioned. // // Subprograms called: // // Linpack dtrsl. // // // References: // // [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited // memory algorithm for bound constrained optimization'', // SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. // // // // * * * // // NEOS, November 1994. (Latest revision June 1996.) // Optimization Technology Center. // Argonne National Laboratory and Northwestern University. // Written by // Ciyou Zhu // in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. // // // ************ CBLAS_INT pointr, m2, col2, ibd, jy, js, i, j, k, one_int = 1; double alpha, xk, dk, temp, temp1, temp2, dd_p; char *uplo = "U", *trans = "T", *diag = "N"; // n, m, col, nsub are size variables. // ind, head are index variables/arrays. if (nsub <= 0) { return; } // Compute wv = W'Zd. pointr = head; for (i = 0; i < col; i++) { temp1 = 0.0; temp2 = 0.0; for (j = 0; j < nsub; j++) { k = ind[j]; temp1 += wy[k + n*pointr]*d[j]; temp2 += ws[k + n*pointr]*d[j]; } wv[i] = temp1; wv[col + i] = theta*temp2; pointr = (pointr + 1) % m; } // Compute wv:=K^(-1)wv. m2 = 2*m; col2 = 2*col; // dtrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info) BLAS_FUNC(dtrtrs)(uplo, trans, diag, &col2, &one_int, wn, &m2, wv, &m2, info); if (*info != 0) { return; } for (i = 0; i < col; i++) { wv[i] = -wv[i]; } trans = "N"; BLAS_FUNC(dtrtrs)(uplo, trans, diag, &col2, &one_int, wn, &m2, wv, &m2, info); if (*info != 0) { return; } // Compute d = (1/theta)d + (1/theta**2)Z'W wv. pointr = head; for (jy = 0; jy < col; jy++) { js = col + jy; for (i = 0; i < nsub; i++) { k = ind[i]; d[i] = d[i] + (wy[k + n*pointr] * wv[jy] / theta) + (ws[k + n*pointr] * wv[js]); } // 30 pointr = (pointr + 1) % m; } // 40 temp = 1.0 / theta; // n, da, dx, incx BLAS_FUNC(dscal)(&nsub, &temp, d, &one_int); // ----------------------------------------------------- // Let us try the projection, d is the Newton direction. *iword = 0; // for (i = 0; i < n; i++) { xp[i] = x[i]; } BLAS_FUNC(dcopy)(&n, x, &one_int, xp, &one_int); for (i = 0; i < nsub; i++) { k = ind[i]; dk = d[i]; xk = x[k]; if (nbd[k] != 0) { if (nbd[k] == 1) // Lower bounds only { x[k] = fmax(l[k], xk + dk); if (x[k] == l[k]) { *iword = 1; } } else { if (nbd[k] == 2) // Upper and lower bounds { xk = fmax(l[k], xk + dk); x[k] = fmin(u[k], xk); if ((x[k] == l[k]) || (x[k] == u[k])) { *iword = 1; } } else { if (nbd[k] == 3) // Upper bounds only { x[k] = fmin(u[k], xk + dk); if (x[k] == u[k]) { *iword = 1; } } } } } else { // Free variables x[k] = xk + dk; } } // 50 if (*iword == 0) { return; } // Check sign of the directional derivative dd_p = 0.0; for (i = 0; i < n; i++) { dd_p = dd_p + (x[i] - xx[i])*gg[i]; } // 55 if (dd_p > 0.0) { // for (i = 0; i < n; i++) { x[i] = xp[i]; } BLAS_FUNC(dcopy)(&n, xp, &one_int, x, &one_int); } else { return; } // ------------------------------------------------------------- alpha = 1.0; temp1 = 1.0; ibd = 0; for (i = 0; i < nsub; i++) { k = ind[i]; dk = d[i]; if (nbd[k] != 0) { if ((dk < 0.0) && (nbd[k] <= 2)) { temp2 = l[k] - x[k]; if (temp2 >= 0.0) { temp1 = 0.0; } else if (dk * alpha < temp2) { temp1 = temp2 / dk; } } else if ((dk > 0.0) && (nbd[k] >= 2)) { temp2 = u[k] - x[k]; if (temp2 <= 0.0) { temp1 = 0.0; } else if (dk * alpha > temp2) { temp1 = temp2 / dk; } } if (temp1 < alpha) { alpha = temp1; ibd = i; } } } // 60 if (alpha < 1.0) { dk = d[ibd]; k = ind[ibd]; if (dk > 0.0) { x[k] = u[k]; d[ibd] = 0.0; } else if (dk < 0.0) { x[k] = l[k]; d[ibd] = 0.0; } } for (i = 0; i < nsub; i++) { k = ind[i]; x[k] = x[k] + alpha*d[i]; } // 70 return; } void dcsrch(double f, double g, double* stp, double ftol, double gtol, double xtol, double stpmin, double stpmax, CBLAS_INT* task, CBLAS_INT* task_msg, CBLAS_INT* isave, double* dsave) { // ********** // // Subroutine dcsrch // // This subroutine finds a step that satisfies a sufficient // decrease condition and a curvature condition. // // Each call of the subroutine updates an interval with // endpoints stx and sty. The interval is initially chosen // so that it contains a minimizer of the modified function // // psi(stp) = f(stp) - f(0) - ftol*stp*f'(0). // // If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the // interval is chosen so that it contains a minimizer of f. // // The algorithm is designed to find a step that satisfies // the sufficient decrease condition // // f(stp) <= f(0) + ftol*stp*f'(0), // // and the curvature condition // // abs(f'(stp)) <= gtol*abs(f'(0)). // // If ftol is less than gtol and if, for example, the function // is bounded below, then there is always a step which satisfies // both conditions. // // If no step can be found that satisfies both conditions, then // the algorithm stops with a warning. In this case stp only // satisfies the sufficient decrease condition. // // A typical invocation of dcsrch has the following outline: // // task = 'START' // 10 continue // call dcsrch( ... ) // if (task .eq. 'FG') then // Evaluate the function and the gradient at stp // goto 10 // end if // // NOTE: The user must no alter work arrays between calls. // // The subroutine statement is // // subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax, // task,isave,dsave) // where // // f is a double precision variable. // On initial entry f is the value of the function at 0. // On subsequent entries f is the value of the // function at stp. // On exit f is the value of the function at stp. // // g is a double precision variable. // On initial entry g is the derivative of the function at 0. // On subsequent entries g is the derivative of the // function at stp. // On exit g is the derivative of the function at stp. // // stp is a double precision variable. // On entry stp is the current estimate of a satisfactory // step. On initial entry, a positive initial estimate // must be provided. // On exit stp is the current estimate of a satisfactory step // if task = 'FG'. If task = 'CONV' then stp satisfies // the sufficient decrease and curvature condition. // // ftol is a double precision variable. // On entry ftol specifies a nonnegative tolerance for the // sufficient decrease condition. // On exit ftol is unchanged. // // gtol is a double precision variable. // On entry gtol specifies a nonnegative tolerance for the // curvature condition. // On exit gtol is unchanged. // // xtol is a double precision variable. // On entry xtol specifies a nonnegative relative tolerance // for an acceptable step. The subroutine exits with a // warning if the relative difference between sty and stx // is less than xtol. // On exit xtol is unchanged. // // stpmin is a double precision variable. // On entry stpmin is a nonnegative lower bound for the step. // On exit stpmin is unchanged. // // stpmax is a double precision variable. // On entry stpmax is a nonnegative upper bound for the step. // On exit stpmax is unchanged. // // task is an integer variable with values from enum Status. // On initial entry task must be set to 'START'. // On exit task indicates the required action: // // If task = 'FG' then evaluate the function and // derivative at stp and call dcsrch again. // // If task = 'CONV' then the search is successful. // // If task = 'WARN' then the subroutine is not able // to satisfy the convergence conditions. The exit value of // stp contains the best point found during the search. // // If task = 'ERROR' then there is an error in the // input arguments. // // On exit with convergence, a warning or an error, the // variable task_msg contains additional information. // // isave is an integer work array of dimension 2. // // dsave is a double precision work array of dimension 13. // // Subprograms called // // MINPACK-2 ... dcstep // // MINPACK-1 Project. June 1983. // Argonne National Laboratory. // Jorge J. More' and David J. Thuente. // // MINPACK-2 Project. October 1993. // Argonne National Laboratory and University of Minnesota. // Brett M. Averick, Richard G. Carter, and Jorge J. More'. // // ********** CBLAS_INT brackt, stage; double finit, ftest, fm, fx, fxm, fy, fym, ginit, gtest, gm, gx, gxm, gy; double gym, stx, sty, stmin, stmax, width, width1; double xtrapl = 1.1; double xtrapu = 4.0; // Initialization block. if (*task == START) { // Check the input arguments for errors. if (*stp < stpmin) { *task = ERROR; *task_msg = ERROR_STP1; } if (*stp > stpmax) { *task = ERROR; *task_msg = ERROR_STP2; } if (g >= 0.0) { *task = ERROR; *task_msg = ERROR_INITG; } if (ftol < 0.0) { *task = ERROR; *task_msg = ERROR_FTOL; } if (gtol < 0.0) { *task = ERROR; *task_msg = ERROR_GTOL; } if (xtol < 0.0) { *task = ERROR; *task_msg = ERROR_XTOL; } if (stpmin < 0.0) { *task = ERROR; *task_msg = ERROR_STPMIN; } if (stpmax < stpmin) { *task = ERROR; *task_msg = ERROR_STPMAX; } // Exit if there are errors on input. if (*task == ERROR) { return; } // Initialize local variables. brackt = 0; stage = 1; finit = f; ginit = g; gtest = ftol * ginit; width = stpmax - stpmin; width1 = width / 0.5; // The variables stx, fx, gx contain the values of the step, function, // and derivative at the best step. // The variables sty, fy, gy contain the value of the step, function, // and derivative at sty. // The variables stp, f, g contain the values of the step, function, // and derivative at stp. stx = 0.0; fx = finit; gx = ginit; sty = 0.0; fy = finit; gy = ginit; stmin = 0.0; stmax = *stp + xtrapu*(*stp); *task = FG; *task_msg = NO_MSG; save_vars(brackt, stage, ginit, gtest, gx, gy, finit, fx, fy, stx, sty, stmin, stmax, width, width1, isave, dsave); return; } else { // Restore local variables. if (isave[0] == 1) { brackt = 1; } else { brackt = 0; } stage = isave[1]; ginit = dsave[0]; gtest = dsave[1]; gx = dsave[2]; gy = dsave[3]; finit = dsave[4]; fx = dsave[5]; fy = dsave[6]; stx = dsave[7]; sty = dsave[8]; stmin = dsave[9]; stmax = dsave[10]; width = dsave[11]; width1 = dsave[12]; } // If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the algorithm // enters the second stage. ftest = finit + (*stp)*gtest; if ((stage == 1) && (f <= ftest) && (g >= 0.0)) { stage = 2; } // Test for warnings. if ((brackt) && ((*stp <= stmin) || (*stp >= stmax))) { *task = WARNING; *task_msg = WARN_ROUND; } if ((brackt) && ((stmax - stmin) <= xtol * stmax)) { *task = WARNING; *task_msg = WARN_XTOL; } if ((*stp == stpmax) && (f <= ftest) && (g <= gtest)) { *task = WARNING; *task_msg = WARN_STPMAX; } if ((*stp == stpmin) && ((f > ftest) || (g >= gtest))) { *task = WARNING; *task_msg = WARN_STPMIN; } // Test for convergence. if ((f <= ftest) && (fabs(g) <= gtol * (-ginit))) { *task = CONVERGENCE; } // Test for termination. if ((*task == WARNING) || (*task == CONVERGENCE)) { save_vars(brackt, stage, ginit, gtest, gx, gy, finit, fx, fy, stx, sty, stmin, stmax, width, width1, isave, dsave); return; } // A modified function is used to predict the step during the first stage if // a lower function value has been obtained but the decrease is not sufficient. if ((stage == 1) && (f <= fx) && (f > ftest)) { // Define the modified function and derivative values. fm = f - (*stp)*gtest; fxm = fx - stx*gtest; fym = fy - sty*gtest; gm = g - gtest; gxm = gx - gtest; gym = gy - gtest; // Call dcstep to update stx, sty, and to compute the new step. dcstep(&stx, &fxm, &gxm, &sty, &fym, &gym, stp, fm, gm, &brackt, stmin, stmax); // Reset the function and derivative values for f. fx = fxm + stx*gtest; fy = fym + sty*gtest; gx = gxm + gtest; gy = gym + gtest; } else { // Call dcstep to update stx, sty, and to compute the new step. dcstep(&stx, &fx, &gx, &sty, &fy, &gy, stp, f, g, &brackt, stmin, stmax); } // Decide if a bisection step is needed. if (brackt) { if (fabs(sty - stx) >= 0.66*width1) { *stp = stx + 0.5*(sty - stx); } width1 = width; width = fabs(sty - stx); } // Set the minimum and maximum steps allowed for stp. if (brackt) { stmin = fmin(stx,sty); stmax = fmax(stx,sty); } else { stmin = *stp + xtrapl*(*stp - stx); stmax = *stp + xtrapu*(*stp - stx); } // Force the step to be within the bounds stpmax and stpmin. *stp = fmax(*stp, stpmin); *stp = fmin(*stp, stpmax); // If further progress is not possible, let stp be the best point obtained // during the search. if (((brackt) && (*stp <= stmin || *stp >= stmax)) || ((brackt) && (stmax - stmin <= xtol * stmax))) { *stp = stx; } // Obtain another function and derivative. *task = FG; *task_msg = NO_MSG; save_vars(brackt, stage, ginit, gtest, gx, gy, finit, fx, fy, stx, sty, stmin, stmax, width, width1, isave, dsave); return; } static inline void save_vars( CBLAS_INT brackt, CBLAS_INT stage, double ginit, double gtest, double gx, double gy, double finit, double fx, double fy, double stx, double sty, double stmin, double stmax, double width, double width1, CBLAS_INT* isave, double* dsave) { if (brackt) { isave[0] = 1; } else { isave[0] = 0; } isave[1] = stage; dsave[0] = ginit; dsave[1] = gtest; dsave[2] = gx; dsave[3] = gy; dsave[4] = finit; dsave[5] = fx; dsave[6] = fy; dsave[7] = stx; dsave[8] = sty; dsave[9] = stmin; dsave[10] = stmax; dsave[11] = width; dsave[12] = width1; return; } void dcstep (double* stx, double* fx, double* dx, double* sty, double* fy, double* dy, double* stp, double fp, double dp, CBLAS_INT* brackt, double stpmin, double stpmax) { // ********** // // Subroutine dcstep // // This subroutine computes a safeguarded step for a search // procedure and updates an interval that contains a step that // satisfies a sufficient decrease and a curvature condition. // // The parameter stx contains the step with the least function // value. If brackt is set to .true. then a minimizer has // been bracketed in an interval with endpoints stx and sty. // The parameter stp contains the current step. // The subroutine assumes that if brackt is set to .true. then // // min(stx,sty) < stp < max(stx,sty), // // and that the derivative at stx is negative in the direction // of the step. // // The subroutine statement is // // subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt, // stpmin,stpmax) // // where // // stx is a double precision variable. // On entry stx is the best step obtained so far and is an // endpoint of the interval that contains the minimizer. // On exit stx is the updated best step. // // fx is a double precision variable. // On entry fx is the function at stx. // On exit fx is the function at stx. // // dx is a double precision variable. // On entry dx is the derivative of the function at // stx. The derivative must be negative in the direction of // the step, that is, dx and stp - stx must have opposite // signs. // On exit dx is the derivative of the function at stx. // // sty is a double precision variable. // On entry sty is the second endpoint of the interval that // contains the minimizer. // On exit sty is the updated endpoint of the interval that // contains the minimizer. // // fy is a double precision variable. // On entry fy is the function at sty. // On exit fy is the function at sty. // // dy is a double precision variable. // On entry dy is the derivative of the function at sty. // On exit dy is the derivative of the function at the exit sty. // // stp is a double precision variable. // On entry stp is the current step. If brackt is set to .true. // then on input stp must be between stx and sty. // On exit stp is a new trial step. // // fp is a double precision variable. // On entry fp is the function at stp // On exit fp is unchanged. // // dp is a double precision variable. // On entry dp is the derivative of the function at stp. // On exit dp is unchanged. // // brackt is a logical variable. // On entry brackt specifies if a minimizer has been bracketed. // Initially brackt must be set to .false. // On exit brackt specifies if a minimizer has been bracketed. // When a minimizer is bracketed brackt is set to .true. // // stpmin is a double precision variable. // On entry stpmin is a lower bound for the step. // On exit stpmin is unchanged. // // stpmax is a double precision variable. // On entry stpmax is an upper bound for the step. // On exit stpmax is unchanged. // // MINPACK-1 Project. June 1983 // Argonne National Laboratory. // Jorge J. More' and David J. Thuente. // // MINPACK-2 Project. October 1993. // Argonne National Laboratory and University of Minnesota. // Brett M. Averick and Jorge J. More'. // // ********** double gamma, p, q, r, s, sgnd, stpc, stpf, stpq, theta; sgnd = dp * (*dx / fabs(*dx)); // First case: A higher function value. The minimum is bracketed. // If the cubic step is closer to stx than the quadratic step, the // cubic step is taken, otherwise the average of the cubic and // quadratic steps is taken. if (fp > *fx) { theta = 3.0 * (*fx - fp) / (*stp - *stx) + *dx + dp; s = fmax(fabs(theta), fmax(fabs(*dx), fabs(dp))); gamma = s * sqrt(pow(theta / s, 2.0) - (*dx / s)*(dp / s)); if (*stp < *stx) { gamma = -gamma; } p = (gamma - *dx) + theta; q = ((gamma - *dx) + gamma) + dp; r = p / q; stpc = *stx + r * (*stp - *stx); stpq = *stx + ((*dx / ((*fx - fp) / (*stp - *stx) + *dx)) / 2.0) * (*stp - *stx); if (fabs(stpc - *stx) < fabs(stpq - *stx)) { stpf = stpc; } else { stpf = stpc + (stpq - stpc) / 2.0; } *brackt = 1; } else if (sgnd < 0.0) { // Second case: A lower function value and derivatives of opposite // sign. The minimum is bracketed. If the cubic step is farther from // stp than the secant step, the cubic step is taken, otherwise the // secant step is taken. theta = 3.0 * (*fx - fp) / (*stp - *stx) + *dx + dp; s = fmax(fabs(theta), fmax(fabs(*dx), fabs(dp))); gamma = s * sqrt(pow(theta / s, 2.0) - (*dx / s)*(dp / s)); if (*stp > *stx) { gamma = -gamma; } p = (gamma - dp) + theta; q = ((gamma - dp) + gamma) + *dx; r = p / q; stpc = *stp + r * (*stx - *stp); stpq = *stp + (dp / (dp - *dx)) * (*stx - *stp); if (fabs(stpc - *stp) > fabs(stpq - *stp)) { stpf = stpc; } else { stpf = stpq; } *brackt = 1; } else if (fabs(dp) < fabs(*dx)) { // Third case: A lower function value, derivatives of the same sign, // and the magnitude of the derivative decreases. // The cubic step is computed only if the cubic tends to infinity // in the direction of the step or if the minimum of the cubic // is beyond stp. Otherwise the cubic step is defined to be the // secant step. theta = 3.0 * (*fx - fp) / (*stp - *stx) + *dx + dp; s = fmax(fabs(theta),fmax(fabs(*dx), fabs(dp))); // The case gamma = 0 only arises if the cubic does not tend // to infinity in the direction of the step. gamma = s*sqrt(fmax(0.0, pow(theta/s, 2.0) - (*dx / s) * (dp / s))); if (*stp > *stx) { gamma = -gamma; } p = (gamma - dp) + theta; q = (gamma + (*dx - dp)) + gamma; r = p / q; if ((r < 0.0) && (gamma != 0.0)) { stpc = *stp + r*(*stx - *stp); } else if (*stp > *stx) { stpc = stpmax; } else { stpc = stpmin; } stpq = *stp + (dp / (dp - *dx)) * (*stx - *stp); if (*brackt) { // A minimizer has been bracketed. If the cubic step is // closer to stp than the secant step, the cubic step is // taken, otherwise the secant step is taken. if (fabs(stpc - *stp) < fabs(stpq - *stp)) { stpf = stpc; } else { stpf = stpq; } if (*stp > *stx) { stpf = fmin(*stp + 0.66*(*sty - *stp), stpf); } else { stpf = fmax(*stp + 0.66*(*sty - *stp), stpf); } } else { // A minimizer has not been bracketed. If the cubic step is // farther from stp than the secant step, the cubic step is // taken, otherwise the secant step is taken. if (fabs(stpc - *stp) > fabs(stpq - *stp)) { stpf = stpc; } else { stpf = stpq; } stpf = fmin(stpmax, stpf); stpf = fmax(stpmin, stpf); } } else { // Fourth case: A lower function value, derivatives of the same sign, // and the magnitude of the derivative does not decrease. If the // minimum is not bracketed, the step is either stpmin or stpmax, // otherwise the cubic step is taken. if (*brackt) { theta = 3.0 * (fp - *fy) / (*sty - *stp) + *dy + dp; s = fmax(fabs(theta), fmax(fabs(*dy), fabs(dp))); gamma = s*sqrt(pow(theta/s, 2.0) - (*dy / s) * (dp / s)); if (*stp > *sty) { gamma = -gamma; } p = (gamma - dp) + theta; q = ((gamma - dp) + gamma) + *dy; r = p / q; stpc = *stp + r*(*sty - *stp); stpf = stpc; } else if (*stp > *stx) { stpf = stpmax; } else { stpf = stpmin; } } // Update the interval which contains a minimizer. if (fp > *fx) { *sty = *stp; *fy = fp; *dy = dp; } else { if (sgnd < 0.0) { *sty = *stx; *fy = *fx; *dy = *dx; } *stx = *stp; *fx = fp; *dx = dp; } // Compute the new step. *stp = stpf; return; }