diff --git a/contrib/lbfgs/solvers.h b/contrib/lbfgs/solvers.h
new file mode 100755
index 0000000000000000000000000000000000000000..2f7e1b5ae6f01b82158c67ec89a026cfb5700788
--- /dev/null
+++ b/contrib/lbfgs/solvers.h
@@ -0,0 +1,1353 @@
+/*************************************************************************
+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
+