diff --git a/contrib/lbfgs/solvers.h b/contrib/lbfgs/solvers.h
deleted file mode 100755
index 2f7e1b5ae6f01b82158c67ec89a026cfb5700788..0000000000000000000000000000000000000000
--- a/contrib/lbfgs/solvers.h
+++ /dev/null
@@ -1,1353 +0,0 @@
-/*************************************************************************
-Copyright (c) Sergey Bochkanov (ALGLIB project).
-
->>> SOURCE LICENSE >>>
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation (www.fsf.org); either version 2 of the
-License, or (at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-GNU General Public License for more details.
-
-A copy of the GNU General Public License is available at
-http://www.fsf.org/licensing/licenses
->>> END OF LICENSE >>>
-*************************************************************************/
-#ifndef _solvers_pkg_h
-#define _solvers_pkg_h
-#include "ap.h"
-#include "alglibinternal.h"
-#include "linalg.h"
-#include "alglibmisc.h"
-
-/////////////////////////////////////////////////////////////////////////
-//
-// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
-//
-/////////////////////////////////////////////////////////////////////////
-namespace alglib_impl
-{
-typedef struct
-{
- double r1;
- double rinf;
-} densesolverreport;
-typedef struct
-{
- double r2;
- ae_matrix cx;
- ae_int_t n;
- ae_int_t k;
-} densesolverlsreport;
-typedef struct
-{
- ae_int_t n;
- ae_int_t m;
- double epsf;
- ae_int_t maxits;
- ae_bool xrep;
- double stpmax;
- ae_vector x;
- double f;
- ae_vector fi;
- ae_matrix j;
- ae_bool needf;
- ae_bool needfij;
- ae_bool xupdated;
- rcommstate rstate;
- ae_int_t repiterationscount;
- ae_int_t repnfunc;
- ae_int_t repnjac;
- ae_int_t repterminationtype;
- ae_vector xbase;
- double fbase;
- double fprev;
- ae_vector candstep;
- ae_vector rightpart;
- ae_vector cgbuf;
-} nleqstate;
-typedef struct
-{
- ae_int_t iterationscount;
- ae_int_t nfunc;
- ae_int_t njac;
- ae_int_t terminationtype;
-} nleqreport;
-
-}
-
-/////////////////////////////////////////////////////////////////////////
-//
-// THIS SECTION CONTAINS C++ INTERFACE
-//
-/////////////////////////////////////////////////////////////////////////
-namespace alglib
-{
-
-/*************************************************************************
-
-*************************************************************************/
-class _densesolverreport_owner
-{
-public:
- _densesolverreport_owner();
- _densesolverreport_owner(const _densesolverreport_owner &rhs);
- _densesolverreport_owner& operator=(const _densesolverreport_owner &rhs);
- virtual ~_densesolverreport_owner();
- alglib_impl::densesolverreport* c_ptr();
- alglib_impl::densesolverreport* c_ptr() const;
-protected:
- alglib_impl::densesolverreport *p_struct;
-};
-class densesolverreport : public _densesolverreport_owner
-{
-public:
- densesolverreport();
- densesolverreport(const densesolverreport &rhs);
- densesolverreport& operator=(const densesolverreport &rhs);
- virtual ~densesolverreport();
- double &r1;
- double &rinf;
-
-};
-
-
-/*************************************************************************
-
-*************************************************************************/
-class _densesolverlsreport_owner
-{
-public:
- _densesolverlsreport_owner();
- _densesolverlsreport_owner(const _densesolverlsreport_owner &rhs);
- _densesolverlsreport_owner& operator=(const _densesolverlsreport_owner &rhs);
- virtual ~_densesolverlsreport_owner();
- alglib_impl::densesolverlsreport* c_ptr();
- alglib_impl::densesolverlsreport* c_ptr() const;
-protected:
- alglib_impl::densesolverlsreport *p_struct;
-};
-class densesolverlsreport : public _densesolverlsreport_owner
-{
-public:
- densesolverlsreport();
- densesolverlsreport(const densesolverlsreport &rhs);
- densesolverlsreport& operator=(const densesolverlsreport &rhs);
- virtual ~densesolverlsreport();
- double &r2;
- real_2d_array cx;
- ae_int_t &n;
- ae_int_t &k;
-
-};
-
-/*************************************************************************
-
-*************************************************************************/
-class _nleqstate_owner
-{
-public:
- _nleqstate_owner();
- _nleqstate_owner(const _nleqstate_owner &rhs);
- _nleqstate_owner& operator=(const _nleqstate_owner &rhs);
- virtual ~_nleqstate_owner();
- alglib_impl::nleqstate* c_ptr();
- alglib_impl::nleqstate* c_ptr() const;
-protected:
- alglib_impl::nleqstate *p_struct;
-};
-class nleqstate : public _nleqstate_owner
-{
-public:
- nleqstate();
- nleqstate(const nleqstate &rhs);
- nleqstate& operator=(const nleqstate &rhs);
- virtual ~nleqstate();
- ae_bool &needf;
- ae_bool &needfij;
- ae_bool &xupdated;
- double &f;
- real_1d_array fi;
- real_2d_array j;
- real_1d_array x;
-
-};
-
-
-/*************************************************************************
-
-*************************************************************************/
-class _nleqreport_owner
-{
-public:
- _nleqreport_owner();
- _nleqreport_owner(const _nleqreport_owner &rhs);
- _nleqreport_owner& operator=(const _nleqreport_owner &rhs);
- virtual ~_nleqreport_owner();
- alglib_impl::nleqreport* c_ptr();
- alglib_impl::nleqreport* c_ptr() const;
-protected:
- alglib_impl::nleqreport *p_struct;
-};
-class nleqreport : public _nleqreport_owner
-{
-public:
- nleqreport();
- nleqreport(const nleqreport &rhs);
- nleqreport& operator=(const nleqreport &rhs);
- virtual ~nleqreport();
- ae_int_t &iterationscount;
- ae_int_t &nfunc;
- ae_int_t &njac;
- ae_int_t &terminationtype;
-
-};
-
-/*************************************************************************
-Dense solver.
-
-This subroutine solves a system A*x=b, where A is NxN non-denegerate
-real matrix, x and b are vectors.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* iterative refinement
-* O(N^3) complexity
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- N - size of A
- B - array[0..N-1], right part
-
-OUTPUT PARAMETERS
- Info - return code:
- * -3 A is singular, or VERY close to singular.
- X is filled by zeros in such cases.
- * -1 N<=0 was passed
- * 1 task is solved (but matrix A may be ill-conditioned,
- check R1/RInf parameters for condition numbers).
- Rep - solver report, see below for more info
- X - array[0..N-1], it contains:
- * solution of A*x=b if A is non-singular (well-conditioned
- or ill-conditioned, but not very close to singular)
- * zeros, if A is singular or VERY close to singular
- (in this case Info=-3).
-
-SOLVER REPORT
-
-Subroutine sets following fields of the Rep structure:
-* R1 reciprocal of condition number: 1/cond(A), 1-norm.
-* RInf reciprocal of condition number: 1/cond(A), inf-norm.
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void rmatrixsolve(const real_2d_array &a, const ae_int_t n, const real_1d_array &b, ae_int_t &info, densesolverreport &rep, real_1d_array &x);
-
-
-/*************************************************************************
-Dense solver.
-
-Similar to RMatrixSolve() but solves task with multiple right parts (where
-b and x are NxM matrices).
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* optional iterative refinement
-* O(N^3+M*N^2) complexity
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- N - size of A
- B - array[0..N-1,0..M-1], right part
- M - right part size
- RFS - iterative refinement switch:
- * True - refinement is used.
- Less performance, more precision.
- * False - refinement is not used.
- More performance, less precision.
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void rmatrixsolvem(const real_2d_array &a, const ae_int_t n, const real_2d_array &b, const ae_int_t m, const bool rfs, ae_int_t &info, densesolverreport &rep, real_2d_array &x);
-
-
-/*************************************************************************
-Dense solver.
-
-This subroutine solves a system A*X=B, where A is NxN non-denegerate
-real matrix given by its LU decomposition, X and B are NxM real matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* O(N^2) complexity
-* condition number estimation
-
-No iterative refinement is provided because exact form of original matrix
-is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
- P - array[0..N-1], pivots array, RMatrixLU result
- N - size of A
- B - array[0..N-1], right part
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void rmatrixlusolve(const real_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const real_1d_array &b, ae_int_t &info, densesolverreport &rep, real_1d_array &x);
-
-
-/*************************************************************************
-Dense solver.
-
-Similar to RMatrixLUSolve() but solves task with multiple right parts
-(where b and x are NxM matrices).
-
-Algorithm features:
-* automatic detection of degenerate cases
-* O(M*N^2) complexity
-* condition number estimation
-
-No iterative refinement is provided because exact form of original matrix
-is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
- P - array[0..N-1], pivots array, RMatrixLU result
- N - size of A
- B - array[0..N-1,0..M-1], right part
- M - right part size
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void rmatrixlusolvem(const real_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const real_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, real_2d_array &x);
-
-
-/*************************************************************************
-Dense solver.
-
-This subroutine solves a system A*x=b, where BOTH ORIGINAL A AND ITS
-LU DECOMPOSITION ARE KNOWN. You can use it if for some reasons you have
-both A and its LU decomposition.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* iterative refinement
-* O(N^2) complexity
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
- P - array[0..N-1], pivots array, RMatrixLU result
- N - size of A
- B - array[0..N-1], right part
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolveM
- Rep - same as in RMatrixSolveM
- X - same as in RMatrixSolveM
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void rmatrixmixedsolve(const real_2d_array &a, const real_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const real_1d_array &b, ae_int_t &info, densesolverreport &rep, real_1d_array &x);
-
-
-/*************************************************************************
-Dense solver.
-
-Similar to RMatrixMixedSolve() but solves task with multiple right parts
-(where b and x are NxM matrices).
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* iterative refinement
-* O(M*N^2) complexity
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
- P - array[0..N-1], pivots array, RMatrixLU result
- N - size of A
- B - array[0..N-1,0..M-1], right part
- M - right part size
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolveM
- Rep - same as in RMatrixSolveM
- X - same as in RMatrixSolveM
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void rmatrixmixedsolvem(const real_2d_array &a, const real_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const real_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, real_2d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixSolveM(), but for complex matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* iterative refinement
-* O(N^3+M*N^2) complexity
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- N - size of A
- B - array[0..N-1,0..M-1], right part
- M - right part size
- RFS - iterative refinement switch:
- * True - refinement is used.
- Less performance, more precision.
- * False - refinement is not used.
- More performance, less precision.
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void cmatrixsolvem(const complex_2d_array &a, const ae_int_t n, const complex_2d_array &b, const ae_int_t m, const bool rfs, ae_int_t &info, densesolverreport &rep, complex_2d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixSolve(), but for complex matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* iterative refinement
-* O(N^3) complexity
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- N - size of A
- B - array[0..N-1], right part
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void cmatrixsolve(const complex_2d_array &a, const ae_int_t n, const complex_1d_array &b, ae_int_t &info, densesolverreport &rep, complex_1d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixLUSolveM(), but for complex matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* O(M*N^2) complexity
-* condition number estimation
-
-No iterative refinement is provided because exact form of original matrix
-is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
- P - array[0..N-1], pivots array, RMatrixLU result
- N - size of A
- B - array[0..N-1,0..M-1], right part
- M - right part size
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void cmatrixlusolvem(const complex_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const complex_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, complex_2d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixLUSolve(), but for complex matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* O(N^2) complexity
-* condition number estimation
-
-No iterative refinement is provided because exact form of original matrix
-is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
- P - array[0..N-1], pivots array, CMatrixLU result
- N - size of A
- B - array[0..N-1], right part
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void cmatrixlusolve(const complex_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const complex_1d_array &b, ae_int_t &info, densesolverreport &rep, complex_1d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixMixedSolveM(), but for complex matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* iterative refinement
-* O(M*N^2) complexity
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
- P - array[0..N-1], pivots array, CMatrixLU result
- N - size of A
- B - array[0..N-1,0..M-1], right part
- M - right part size
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolveM
- Rep - same as in RMatrixSolveM
- X - same as in RMatrixSolveM
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void cmatrixmixedsolvem(const complex_2d_array &a, const complex_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const complex_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, complex_2d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixMixedSolve(), but for complex matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* iterative refinement
-* O(N^2) complexity
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
- P - array[0..N-1], pivots array, CMatrixLU result
- N - size of A
- B - array[0..N-1], right part
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolveM
- Rep - same as in RMatrixSolveM
- X - same as in RMatrixSolveM
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void cmatrixmixedsolve(const complex_2d_array &a, const complex_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const complex_1d_array &b, ae_int_t &info, densesolverreport &rep, complex_1d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixSolveM(), but for symmetric positive definite
-matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* O(N^3+M*N^2) complexity
-* matrix is represented by its upper or lower triangle
-
-No iterative refinement is provided because such partial representation of
-matrix does not allow efficient calculation of extra-precise matrix-vector
-products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- N - size of A
- IsUpper - what half of A is provided
- B - array[0..N-1,0..M-1], right part
- M - right part size
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve.
- Returns -3 for non-SPD matrices.
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void spdmatrixsolvem(const real_2d_array &a, const ae_int_t n, const bool isupper, const real_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, real_2d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixSolve(), but for SPD matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* O(N^3) complexity
-* matrix is represented by its upper or lower triangle
-
-No iterative refinement is provided because such partial representation of
-matrix does not allow efficient calculation of extra-precise matrix-vector
-products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- N - size of A
- IsUpper - what half of A is provided
- B - array[0..N-1], right part
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Returns -3 for non-SPD matrices.
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void spdmatrixsolve(const real_2d_array &a, const ae_int_t n, const bool isupper, const real_1d_array &b, ae_int_t &info, densesolverreport &rep, real_1d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixLUSolveM(), but for SPD matrices represented
-by their Cholesky decomposition.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* O(M*N^2) complexity
-* condition number estimation
-* matrix is represented by its upper or lower triangle
-
-No iterative refinement is provided because such partial representation of
-matrix does not allow efficient calculation of extra-precise matrix-vector
-products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- CHA - array[0..N-1,0..N-1], Cholesky decomposition,
- SPDMatrixCholesky result
- N - size of CHA
- IsUpper - what half of CHA is provided
- B - array[0..N-1,0..M-1], right part
- M - right part size
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void spdmatrixcholeskysolvem(const real_2d_array &cha, const ae_int_t n, const bool isupper, const real_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, real_2d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixLUSolve(), but for SPD matrices represented
-by their Cholesky decomposition.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* O(N^2) complexity
-* condition number estimation
-* matrix is represented by its upper or lower triangle
-
-No iterative refinement is provided because such partial representation of
-matrix does not allow efficient calculation of extra-precise matrix-vector
-products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- CHA - array[0..N-1,0..N-1], Cholesky decomposition,
- SPDMatrixCholesky result
- N - size of A
- IsUpper - what half of CHA is provided
- B - array[0..N-1], right part
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void spdmatrixcholeskysolve(const real_2d_array &cha, const ae_int_t n, const bool isupper, const real_1d_array &b, ae_int_t &info, densesolverreport &rep, real_1d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixSolveM(), but for Hermitian positive definite
-matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* O(N^3+M*N^2) complexity
-* matrix is represented by its upper or lower triangle
-
-No iterative refinement is provided because such partial representation of
-matrix does not allow efficient calculation of extra-precise matrix-vector
-products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- N - size of A
- IsUpper - what half of A is provided
- B - array[0..N-1,0..M-1], right part
- M - right part size
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve.
- Returns -3 for non-HPD matrices.
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void hpdmatrixsolvem(const complex_2d_array &a, const ae_int_t n, const bool isupper, const complex_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, complex_2d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixSolve(), but for Hermitian positive definite
-matrices.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* condition number estimation
-* O(N^3) complexity
-* matrix is represented by its upper or lower triangle
-
-No iterative refinement is provided because such partial representation of
-matrix does not allow efficient calculation of extra-precise matrix-vector
-products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- A - array[0..N-1,0..N-1], system matrix
- N - size of A
- IsUpper - what half of A is provided
- B - array[0..N-1], right part
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Returns -3 for non-HPD matrices.
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void hpdmatrixsolve(const complex_2d_array &a, const ae_int_t n, const bool isupper, const complex_1d_array &b, ae_int_t &info, densesolverreport &rep, complex_1d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixLUSolveM(), but for HPD matrices represented
-by their Cholesky decomposition.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* O(M*N^2) complexity
-* condition number estimation
-* matrix is represented by its upper or lower triangle
-
-No iterative refinement is provided because such partial representation of
-matrix does not allow efficient calculation of extra-precise matrix-vector
-products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- CHA - array[0..N-1,0..N-1], Cholesky decomposition,
- HPDMatrixCholesky result
- N - size of CHA
- IsUpper - what half of CHA is provided
- B - array[0..N-1,0..M-1], right part
- M - right part size
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void hpdmatrixcholeskysolvem(const complex_2d_array &cha, const ae_int_t n, const bool isupper, const complex_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, complex_2d_array &x);
-
-
-/*************************************************************************
-Dense solver. Same as RMatrixLUSolve(), but for HPD matrices represented
-by their Cholesky decomposition.
-
-Algorithm features:
-* automatic detection of degenerate cases
-* O(N^2) complexity
-* condition number estimation
-* matrix is represented by its upper or lower triangle
-
-No iterative refinement is provided because such partial representation of
-matrix does not allow efficient calculation of extra-precise matrix-vector
-products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
-need iterative refinement.
-
-INPUT PARAMETERS
- CHA - array[0..N-1,0..N-1], Cholesky decomposition,
- SPDMatrixCholesky result
- N - size of A
- IsUpper - what half of CHA is provided
- B - array[0..N-1], right part
-
-OUTPUT PARAMETERS
- Info - same as in RMatrixSolve
- Rep - same as in RMatrixSolve
- X - same as in RMatrixSolve
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-void hpdmatrixcholeskysolve(const complex_2d_array &cha, const ae_int_t n, const bool isupper, const complex_1d_array &b, ae_int_t &info, densesolverreport &rep, complex_1d_array &x);
-
-
-/*************************************************************************
-Dense solver.
-
-This subroutine finds solution of the linear system A*X=B with non-square,
-possibly degenerate A. System is solved in the least squares sense, and
-general least squares solution X = X0 + CX*y which minimizes |A*X-B| is
-returned. If A is non-degenerate, solution in the usual sense is returned
-
-Algorithm features:
-* automatic detection of degenerate cases
-* iterative refinement
-* O(N^3) complexity
-
-INPUT PARAMETERS
- A - array[0..NRows-1,0..NCols-1], system matrix
- NRows - vertical size of A
- NCols - horizontal size of A
- B - array[0..NCols-1], right part
- Threshold- a number in [0,1]. Singular values beyond Threshold are
- considered zero. Set it to 0.0, if you don't understand
- what it means, so the solver will choose good value on its
- own.
-
-OUTPUT PARAMETERS
- Info - return code:
- * -4 SVD subroutine failed
- * -1 if NRows<=0 or NCols<=0 or Threshold<0 was passed
- * 1 if task is solved
- Rep - solver report, see below for more info
- X - array[0..N-1,0..M-1], it contains:
- * solution of A*X=B if A is non-singular (well-conditioned
- or ill-conditioned, but not very close to singular)
- * zeros, if A is singular or VERY close to singular
- (in this case Info=-3).
-
-SOLVER REPORT
-
-Subroutine sets following fields of the Rep structure:
-* R2 reciprocal of condition number: 1/cond(A), 2-norm.
-* N = NCols
-* K dim(Null(A))
-* CX array[0..N-1,0..K-1], kernel of A.
- Columns of CX store such vectors that A*CX[i]=0.
-
- -- ALGLIB --
- Copyright 24.08.2009 by Bochkanov Sergey
-*************************************************************************/
-void rmatrixsolvels(const real_2d_array &a, const ae_int_t nrows, const ae_int_t ncols, const real_1d_array &b, const double threshold, ae_int_t &info, densesolverlsreport &rep, real_1d_array &x);
-
-/*************************************************************************
- LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER
-
-DESCRIPTION:
-This algorithm solves system of nonlinear equations
- F[0](x[0], ..., x[n-1]) = 0
- F[1](x[0], ..., x[n-1]) = 0
- ...
- F[M-1](x[0], ..., x[n-1]) = 0
-with M/N do not necessarily coincide. Algorithm converges quadratically
-under following conditions:
- * the solution set XS is nonempty
- * for some xs in XS there exist such neighbourhood N(xs) that:
- * vector function F(x) and its Jacobian J(x) are continuously
- differentiable on N
- * ||F(x)|| provides local error bound on N, i.e. there exists such
- c1, that ||F(x)||>c1*distance(x,XS)
-Note that these conditions are much more weaker than usual non-singularity
-conditions. For example, algorithm will converge for any affine function
-F (whether its Jacobian singular or not).
-
-
-REQUIREMENTS:
-Algorithm will request following information during its operation:
-* function vector F[] and Jacobian matrix at given point X
-* value of merit function f(x)=F[0]^2(x)+...+F[M-1]^2(x) at given point X
-
-
-USAGE:
-1. User initializes algorithm state with NLEQCreateLM() call
-2. User tunes solver parameters with NLEQSetCond(), NLEQSetStpMax() and
- other functions
-3. User calls NLEQSolve() function which takes algorithm state and
- pointers (delegates, etc.) to callback functions which calculate merit
- function value and Jacobian.
-4. User calls NLEQResults() to get solution
-5. Optionally, user may call NLEQRestartFrom() to solve another problem
- with same parameters (N/M) but another starting point and/or another
- function vector. NLEQRestartFrom() allows to reuse already initialized
- structure.
-
-
-INPUT PARAMETERS:
- N - space dimension, N>1:
- * if provided, only leading N elements of X are used
- * if not provided, determined automatically from size of X
- M - system size
- X - starting point
-
-
-OUTPUT PARAMETERS:
- State - structure which stores algorithm state
-
-
-NOTES:
-1. you may tune stopping conditions with NLEQSetCond() function
-2. if target function contains exp() or other fast growing functions, and
- optimization algorithm makes too large steps which leads to overflow,
- use NLEQSetStpMax() function to bound algorithm's steps.
-3. this algorithm is a slightly modified implementation of the method
- described in 'Levenberg-Marquardt method for constrained nonlinear
- equations with strong local convergence properties' by Christian Kanzow
- Nobuo Yamashita and Masao Fukushima and further developed in 'On the
- convergence of a New Levenberg-Marquardt Method' by Jin-yan Fan and
- Ya-Xiang Yuan.
-
-
- -- ALGLIB --
- Copyright 20.08.2009 by Bochkanov Sergey
-*************************************************************************/
-void nleqcreatelm(const ae_int_t n, const ae_int_t m, const real_1d_array &x, nleqstate &state);
-void nleqcreatelm(const ae_int_t m, const real_1d_array &x, nleqstate &state);
-
-
-/*************************************************************************
-This function sets stopping conditions for the nonlinear solver
-
-INPUT PARAMETERS:
- State - structure which stores algorithm state
- EpsF - >=0
- The subroutine finishes its work if on k+1-th iteration
- the condition ||F||<=EpsF is satisfied
- MaxIts - maximum number of iterations. If MaxIts=0, the number of
- iterations is unlimited.
-
-Passing EpsF=0 and MaxIts=0 simultaneously will lead to automatic
-stopping criterion selection (small EpsF).
-
-NOTES:
-
- -- ALGLIB --
- Copyright 20.08.2010 by Bochkanov Sergey
-*************************************************************************/
-void nleqsetcond(const nleqstate &state, const double epsf, const ae_int_t maxits);
-
-
-/*************************************************************************
-This function turns on/off reporting.
-
-INPUT PARAMETERS:
- State - structure which stores algorithm state
- NeedXRep- whether iteration reports are needed or not
-
-If NeedXRep is True, algorithm will call rep() callback function if it is
-provided to NLEQSolve().
-
- -- ALGLIB --
- Copyright 20.08.2010 by Bochkanov Sergey
-*************************************************************************/
-void nleqsetxrep(const nleqstate &state, const bool needxrep);
-
-
-/*************************************************************************
-This function sets maximum step length
-
-INPUT PARAMETERS:
- State - structure which stores algorithm state
- StpMax - maximum step length, >=0. Set StpMax to 0.0, if you don't
- want to limit step length.
-
-Use this subroutine when target function contains exp() or other fast
-growing functions, and algorithm makes too large steps which lead to
-overflow. This function allows us to reject steps that are too large (and
-therefore expose us to the possible overflow) without actually calculating
-function value at the x+stp*d.
-
- -- ALGLIB --
- Copyright 20.08.2010 by Bochkanov Sergey
-*************************************************************************/
-void nleqsetstpmax(const nleqstate &state, const double stpmax);
-
-
-/*************************************************************************
-This function provides reverse communication interface
-Reverse communication interface is not documented or recommended to use.
-See below for functions which provide better documented API
-*************************************************************************/
-bool nleqiteration(const nleqstate &state);
-
-
-/*************************************************************************
-This family of functions is used to launcn iterations of nonlinear solver
-
-These functions accept following parameters:
- state - algorithm state
- func - callback which calculates function (or merit function)
- value func at given point x
- jac - callback which calculates function vector fi[]
- and Jacobian jac at given point x
- rep - optional callback which is called after each iteration
- can be NULL
- ptr - optional pointer which is passed to func/grad/hess/jac/rep
- can be NULL
-
-
- -- ALGLIB --
- Copyright 20.03.2009 by Bochkanov Sergey
-
-*************************************************************************/
-void nleqsolve(nleqstate &state,
- void (*func)(const real_1d_array &x, double &func, void *ptr),
- void (*jac)(const real_1d_array &x, real_1d_array &fi, real_2d_array &jac, void *ptr),
- void (*rep)(const real_1d_array &x, double func, void *ptr) = NULL,
- void *ptr = NULL);
-
-
-/*************************************************************************
-NLEQ solver results
-
-INPUT PARAMETERS:
- State - algorithm state.
-
-OUTPUT PARAMETERS:
- X - array[0..N-1], solution
- Rep - optimization report:
- * Rep.TerminationType completetion code:
- * -4 ERROR: algorithm has converged to the
- stationary point Xf which is local minimum of
- f=F[0]^2+...+F[m-1]^2, but is not solution of
- nonlinear system.
- * 1 sqrt(f)<=EpsF.
- * 5 MaxIts steps was taken
- * 7 stopping conditions are too stringent,
- further improvement is impossible
- * Rep.IterationsCount contains iterations count
- * NFEV countains number of function calculations
- * ActiveConstraints contains number of active constraints
-
- -- ALGLIB --
- Copyright 20.08.2009 by Bochkanov Sergey
-*************************************************************************/
-void nleqresults(const nleqstate &state, real_1d_array &x, nleqreport &rep);
-
-
-/*************************************************************************
-NLEQ solver results
-
-Buffered implementation of NLEQResults(), which uses pre-allocated buffer
-to store X[]. If buffer size is too small, it resizes buffer. It is
-intended to be used in the inner cycles of performance critical algorithms
-where array reallocation penalty is too large to be ignored.
-
- -- ALGLIB --
- Copyright 20.08.2009 by Bochkanov Sergey
-*************************************************************************/
-void nleqresultsbuf(const nleqstate &state, real_1d_array &x, nleqreport &rep);
-
-
-/*************************************************************************
-This subroutine restarts CG algorithm from new point. All optimization
-parameters are left unchanged.
-
-This function allows to solve multiple optimization problems (which
-must have same number of dimensions) without object reallocation penalty.
-
-INPUT PARAMETERS:
- State - structure used for reverse communication previously
- allocated with MinCGCreate call.
- X - new starting point.
- BndL - new lower bounds
- BndU - new upper bounds
-
- -- ALGLIB --
- Copyright 30.07.2010 by Bochkanov Sergey
-*************************************************************************/
-void nleqrestartfrom(const nleqstate &state, const real_1d_array &x);
-}
-
-/////////////////////////////////////////////////////////////////////////
-//
-// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
-//
-/////////////////////////////////////////////////////////////////////////
-namespace alglib_impl
-{
-void rmatrixsolve(/* Real */ ae_matrix* a,
- ae_int_t n,
- /* Real */ ae_vector* b,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_vector* x,
- ae_state *_state);
-void rmatrixsolvem(/* Real */ ae_matrix* a,
- ae_int_t n,
- /* Real */ ae_matrix* b,
- ae_int_t m,
- ae_bool rfs,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_matrix* x,
- ae_state *_state);
-void rmatrixlusolve(/* Real */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- ae_int_t n,
- /* Real */ ae_vector* b,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_vector* x,
- ae_state *_state);
-void rmatrixlusolvem(/* Real */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- ae_int_t n,
- /* Real */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_matrix* x,
- ae_state *_state);
-void rmatrixmixedsolve(/* Real */ ae_matrix* a,
- /* Real */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- ae_int_t n,
- /* Real */ ae_vector* b,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_vector* x,
- ae_state *_state);
-void rmatrixmixedsolvem(/* Real */ ae_matrix* a,
- /* Real */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- ae_int_t n,
- /* Real */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_matrix* x,
- ae_state *_state);
-void cmatrixsolvem(/* Complex */ ae_matrix* a,
- ae_int_t n,
- /* Complex */ ae_matrix* b,
- ae_int_t m,
- ae_bool rfs,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_matrix* x,
- ae_state *_state);
-void cmatrixsolve(/* Complex */ ae_matrix* a,
- ae_int_t n,
- /* Complex */ ae_vector* b,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_vector* x,
- ae_state *_state);
-void cmatrixlusolvem(/* Complex */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- ae_int_t n,
- /* Complex */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_matrix* x,
- ae_state *_state);
-void cmatrixlusolve(/* Complex */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- ae_int_t n,
- /* Complex */ ae_vector* b,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_vector* x,
- ae_state *_state);
-void cmatrixmixedsolvem(/* Complex */ ae_matrix* a,
- /* Complex */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- ae_int_t n,
- /* Complex */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_matrix* x,
- ae_state *_state);
-void cmatrixmixedsolve(/* Complex */ ae_matrix* a,
- /* Complex */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- ae_int_t n,
- /* Complex */ ae_vector* b,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_vector* x,
- ae_state *_state);
-void spdmatrixsolvem(/* Real */ ae_matrix* a,
- ae_int_t n,
- ae_bool isupper,
- /* Real */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_matrix* x,
- ae_state *_state);
-void spdmatrixsolve(/* Real */ ae_matrix* a,
- ae_int_t n,
- ae_bool isupper,
- /* Real */ ae_vector* b,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_vector* x,
- ae_state *_state);
-void spdmatrixcholeskysolvem(/* Real */ ae_matrix* cha,
- ae_int_t n,
- ae_bool isupper,
- /* Real */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_matrix* x,
- ae_state *_state);
-void spdmatrixcholeskysolve(/* Real */ ae_matrix* cha,
- ae_int_t n,
- ae_bool isupper,
- /* Real */ ae_vector* b,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_vector* x,
- ae_state *_state);
-void hpdmatrixsolvem(/* Complex */ ae_matrix* a,
- ae_int_t n,
- ae_bool isupper,
- /* Complex */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_matrix* x,
- ae_state *_state);
-void hpdmatrixsolve(/* Complex */ ae_matrix* a,
- ae_int_t n,
- ae_bool isupper,
- /* Complex */ ae_vector* b,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_vector* x,
- ae_state *_state);
-void hpdmatrixcholeskysolvem(/* Complex */ ae_matrix* cha,
- ae_int_t n,
- ae_bool isupper,
- /* Complex */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_matrix* x,
- ae_state *_state);
-void hpdmatrixcholeskysolve(/* Complex */ ae_matrix* cha,
- ae_int_t n,
- ae_bool isupper,
- /* Complex */ ae_vector* b,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_vector* x,
- ae_state *_state);
-void rmatrixsolvels(/* Real */ ae_matrix* a,
- ae_int_t nrows,
- ae_int_t ncols,
- /* Real */ ae_vector* b,
- double threshold,
- ae_int_t* info,
- densesolverlsreport* rep,
- /* Real */ ae_vector* x,
- ae_state *_state);
-ae_bool _densesolverreport_init(densesolverreport* p, ae_state *_state, ae_bool make_automatic);
-ae_bool _densesolverreport_init_copy(densesolverreport* dst, densesolverreport* src, ae_state *_state, ae_bool make_automatic);
-void _densesolverreport_clear(densesolverreport* p);
-ae_bool _densesolverlsreport_init(densesolverlsreport* p, ae_state *_state, ae_bool make_automatic);
-ae_bool _densesolverlsreport_init_copy(densesolverlsreport* dst, densesolverlsreport* src, ae_state *_state, ae_bool make_automatic);
-void _densesolverlsreport_clear(densesolverlsreport* p);
-void nleqcreatelm(ae_int_t n,
- ae_int_t m,
- /* Real */ ae_vector* x,
- nleqstate* state,
- ae_state *_state);
-void nleqsetcond(nleqstate* state,
- double epsf,
- ae_int_t maxits,
- ae_state *_state);
-void nleqsetxrep(nleqstate* state, ae_bool needxrep, ae_state *_state);
-void nleqsetstpmax(nleqstate* state, double stpmax, ae_state *_state);
-ae_bool nleqiteration(nleqstate* state, ae_state *_state);
-void nleqresults(nleqstate* state,
- /* Real */ ae_vector* x,
- nleqreport* rep,
- ae_state *_state);
-void nleqresultsbuf(nleqstate* state,
- /* Real */ ae_vector* x,
- nleqreport* rep,
- ae_state *_state);
-void nleqrestartfrom(nleqstate* state,
- /* Real */ ae_vector* x,
- ae_state *_state);
-ae_bool _nleqstate_init(nleqstate* p, ae_state *_state, ae_bool make_automatic);
-ae_bool _nleqstate_init_copy(nleqstate* dst, nleqstate* src, ae_state *_state, ae_bool make_automatic);
-void _nleqstate_clear(nleqstate* p);
-ae_bool _nleqreport_init(nleqreport* p, ae_state *_state, ae_bool make_automatic);
-ae_bool _nleqreport_init_copy(nleqreport* dst, nleqreport* src, ae_state *_state, ae_bool make_automatic);
-void _nleqreport_clear(nleqreport* p);
-
-}
-#endif
-