diff --git a/contrib/lbfgs/solvers.cpp b/contrib/lbfgs/solvers.cpp
deleted file mode 100755
index 595486b309f676d0fd0c3d51cfe3d2b8e1564ca1..0000000000000000000000000000000000000000
--- a/contrib/lbfgs/solvers.cpp
+++ /dev/null
@@ -1,5665 +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 >>>
-*************************************************************************/
-#include "stdafx.h"
-#include "solvers.h"
-
-// disable some irrelevant warnings
-#if (AE_COMPILER==AE_MSVC)
-#pragma warning(disable:4100)
-#pragma warning(disable:4127)
-#pragma warning(disable:4702)
-#pragma warning(disable:4996)
-#endif
-using namespace std;
-
-/////////////////////////////////////////////////////////////////////////
-//
-// THIS SECTION CONTAINS IMPLEMENTATION OF C++ INTERFACE
-//
-/////////////////////////////////////////////////////////////////////////
-namespace alglib
-{
-
-
-/*************************************************************************
-
-*************************************************************************/
-_densesolverreport_owner::_densesolverreport_owner()
-{
- p_struct = (alglib_impl::densesolverreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::densesolverreport), NULL);
- if( p_struct==NULL )
- throw ap_error("ALGLIB: malloc error");
- if( !alglib_impl::_densesolverreport_init(p_struct, NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
-}
-
-_densesolverreport_owner::_densesolverreport_owner(const _densesolverreport_owner &rhs)
-{
- p_struct = (alglib_impl::densesolverreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::densesolverreport), NULL);
- if( p_struct==NULL )
- throw ap_error("ALGLIB: malloc error");
- if( !alglib_impl::_densesolverreport_init_copy(p_struct, const_cast<alglib_impl::densesolverreport*>(rhs.p_struct), NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
-}
-
-_densesolverreport_owner& _densesolverreport_owner::operator=(const _densesolverreport_owner &rhs)
-{
- if( this==&rhs )
- return *this;
- alglib_impl::_densesolverreport_clear(p_struct);
- if( !alglib_impl::_densesolverreport_init_copy(p_struct, const_cast<alglib_impl::densesolverreport*>(rhs.p_struct), NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
- return *this;
-}
-
-_densesolverreport_owner::~_densesolverreport_owner()
-{
- alglib_impl::_densesolverreport_clear(p_struct);
- ae_free(p_struct);
-}
-
-alglib_impl::densesolverreport* _densesolverreport_owner::c_ptr()
-{
- return p_struct;
-}
-
-alglib_impl::densesolverreport* _densesolverreport_owner::c_ptr() const
-{
- return const_cast<alglib_impl::densesolverreport*>(p_struct);
-}
-densesolverreport::densesolverreport() : _densesolverreport_owner() ,r1(p_struct->r1),rinf(p_struct->rinf)
-{
-}
-
-densesolverreport::densesolverreport(const densesolverreport &rhs):_densesolverreport_owner(rhs) ,r1(p_struct->r1),rinf(p_struct->rinf)
-{
-}
-
-densesolverreport& densesolverreport::operator=(const densesolverreport &rhs)
-{
- if( this==&rhs )
- return *this;
- _densesolverreport_owner::operator=(rhs);
- return *this;
-}
-
-densesolverreport::~densesolverreport()
-{
-}
-
-
-/*************************************************************************
-
-*************************************************************************/
-_densesolverlsreport_owner::_densesolverlsreport_owner()
-{
- p_struct = (alglib_impl::densesolverlsreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::densesolverlsreport), NULL);
- if( p_struct==NULL )
- throw ap_error("ALGLIB: malloc error");
- if( !alglib_impl::_densesolverlsreport_init(p_struct, NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
-}
-
-_densesolverlsreport_owner::_densesolverlsreport_owner(const _densesolverlsreport_owner &rhs)
-{
- p_struct = (alglib_impl::densesolverlsreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::densesolverlsreport), NULL);
- if( p_struct==NULL )
- throw ap_error("ALGLIB: malloc error");
- if( !alglib_impl::_densesolverlsreport_init_copy(p_struct, const_cast<alglib_impl::densesolverlsreport*>(rhs.p_struct), NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
-}
-
-_densesolverlsreport_owner& _densesolverlsreport_owner::operator=(const _densesolverlsreport_owner &rhs)
-{
- if( this==&rhs )
- return *this;
- alglib_impl::_densesolverlsreport_clear(p_struct);
- if( !alglib_impl::_densesolverlsreport_init_copy(p_struct, const_cast<alglib_impl::densesolverlsreport*>(rhs.p_struct), NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
- return *this;
-}
-
-_densesolverlsreport_owner::~_densesolverlsreport_owner()
-{
- alglib_impl::_densesolverlsreport_clear(p_struct);
- ae_free(p_struct);
-}
-
-alglib_impl::densesolverlsreport* _densesolverlsreport_owner::c_ptr()
-{
- return p_struct;
-}
-
-alglib_impl::densesolverlsreport* _densesolverlsreport_owner::c_ptr() const
-{
- return const_cast<alglib_impl::densesolverlsreport*>(p_struct);
-}
-densesolverlsreport::densesolverlsreport() : _densesolverlsreport_owner() ,r2(p_struct->r2),cx(&p_struct->cx),n(p_struct->n),k(p_struct->k)
-{
-}
-
-densesolverlsreport::densesolverlsreport(const densesolverlsreport &rhs):_densesolverlsreport_owner(rhs) ,r2(p_struct->r2),cx(&p_struct->cx),n(p_struct->n),k(p_struct->k)
-{
-}
-
-densesolverlsreport& densesolverlsreport::operator=(const densesolverlsreport &rhs)
-{
- if( this==&rhs )
- return *this;
- _densesolverlsreport_owner::operator=(rhs);
- return *this;
-}
-
-densesolverlsreport::~densesolverlsreport()
-{
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::rmatrixsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::rmatrixsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, rfs, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::rmatrixlusolve(const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::rmatrixlusolvem(const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::rmatrixmixedsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::rmatrixmixedsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::cmatrixsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, rfs, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::cmatrixsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::cmatrixlusolvem(const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::cmatrixlusolve(const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::cmatrixmixedsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::cmatrixmixedsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::spdmatrixsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, isupper, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::spdmatrixsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, isupper, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::spdmatrixcholeskysolvem(const_cast<alglib_impl::ae_matrix*>(cha.c_ptr()), n, isupper, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::spdmatrixcholeskysolve(const_cast<alglib_impl::ae_matrix*>(cha.c_ptr()), n, isupper, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::hpdmatrixsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, isupper, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::hpdmatrixsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, isupper, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::hpdmatrixcholeskysolvem(const_cast<alglib_impl::ae_matrix*>(cha.c_ptr()), n, isupper, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::hpdmatrixcholeskysolve(const_cast<alglib_impl::ae_matrix*>(cha.c_ptr()), n, isupper, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::rmatrixsolvels(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), nrows, ncols, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), threshold, &info, const_cast<alglib_impl::densesolverlsreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-
-*************************************************************************/
-_nleqstate_owner::_nleqstate_owner()
-{
- p_struct = (alglib_impl::nleqstate*)alglib_impl::ae_malloc(sizeof(alglib_impl::nleqstate), NULL);
- if( p_struct==NULL )
- throw ap_error("ALGLIB: malloc error");
- if( !alglib_impl::_nleqstate_init(p_struct, NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
-}
-
-_nleqstate_owner::_nleqstate_owner(const _nleqstate_owner &rhs)
-{
- p_struct = (alglib_impl::nleqstate*)alglib_impl::ae_malloc(sizeof(alglib_impl::nleqstate), NULL);
- if( p_struct==NULL )
- throw ap_error("ALGLIB: malloc error");
- if( !alglib_impl::_nleqstate_init_copy(p_struct, const_cast<alglib_impl::nleqstate*>(rhs.p_struct), NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
-}
-
-_nleqstate_owner& _nleqstate_owner::operator=(const _nleqstate_owner &rhs)
-{
- if( this==&rhs )
- return *this;
- alglib_impl::_nleqstate_clear(p_struct);
- if( !alglib_impl::_nleqstate_init_copy(p_struct, const_cast<alglib_impl::nleqstate*>(rhs.p_struct), NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
- return *this;
-}
-
-_nleqstate_owner::~_nleqstate_owner()
-{
- alglib_impl::_nleqstate_clear(p_struct);
- ae_free(p_struct);
-}
-
-alglib_impl::nleqstate* _nleqstate_owner::c_ptr()
-{
- return p_struct;
-}
-
-alglib_impl::nleqstate* _nleqstate_owner::c_ptr() const
-{
- return const_cast<alglib_impl::nleqstate*>(p_struct);
-}
-nleqstate::nleqstate() : _nleqstate_owner() ,needf(p_struct->needf),needfij(p_struct->needfij),xupdated(p_struct->xupdated),f(p_struct->f),fi(&p_struct->fi),j(&p_struct->j),x(&p_struct->x)
-{
-}
-
-nleqstate::nleqstate(const nleqstate &rhs):_nleqstate_owner(rhs) ,needf(p_struct->needf),needfij(p_struct->needfij),xupdated(p_struct->xupdated),f(p_struct->f),fi(&p_struct->fi),j(&p_struct->j),x(&p_struct->x)
-{
-}
-
-nleqstate& nleqstate::operator=(const nleqstate &rhs)
-{
- if( this==&rhs )
- return *this;
- _nleqstate_owner::operator=(rhs);
- return *this;
-}
-
-nleqstate::~nleqstate()
-{
-}
-
-
-/*************************************************************************
-
-*************************************************************************/
-_nleqreport_owner::_nleqreport_owner()
-{
- p_struct = (alglib_impl::nleqreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::nleqreport), NULL);
- if( p_struct==NULL )
- throw ap_error("ALGLIB: malloc error");
- if( !alglib_impl::_nleqreport_init(p_struct, NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
-}
-
-_nleqreport_owner::_nleqreport_owner(const _nleqreport_owner &rhs)
-{
- p_struct = (alglib_impl::nleqreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::nleqreport), NULL);
- if( p_struct==NULL )
- throw ap_error("ALGLIB: malloc error");
- if( !alglib_impl::_nleqreport_init_copy(p_struct, const_cast<alglib_impl::nleqreport*>(rhs.p_struct), NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
-}
-
-_nleqreport_owner& _nleqreport_owner::operator=(const _nleqreport_owner &rhs)
-{
- if( this==&rhs )
- return *this;
- alglib_impl::_nleqreport_clear(p_struct);
- if( !alglib_impl::_nleqreport_init_copy(p_struct, const_cast<alglib_impl::nleqreport*>(rhs.p_struct), NULL, ae_false) )
- throw ap_error("ALGLIB: malloc error");
- return *this;
-}
-
-_nleqreport_owner::~_nleqreport_owner()
-{
- alglib_impl::_nleqreport_clear(p_struct);
- ae_free(p_struct);
-}
-
-alglib_impl::nleqreport* _nleqreport_owner::c_ptr()
-{
- return p_struct;
-}
-
-alglib_impl::nleqreport* _nleqreport_owner::c_ptr() const
-{
- return const_cast<alglib_impl::nleqreport*>(p_struct);
-}
-nleqreport::nleqreport() : _nleqreport_owner() ,iterationscount(p_struct->iterationscount),nfunc(p_struct->nfunc),njac(p_struct->njac),terminationtype(p_struct->terminationtype)
-{
-}
-
-nleqreport::nleqreport(const nleqreport &rhs):_nleqreport_owner(rhs) ,iterationscount(p_struct->iterationscount),nfunc(p_struct->nfunc),njac(p_struct->njac),terminationtype(p_struct->terminationtype)
-{
-}
-
-nleqreport& nleqreport::operator=(const nleqreport &rhs)
-{
- if( this==&rhs )
- return *this;
- _nleqreport_owner::operator=(rhs);
- return *this;
-}
-
-nleqreport::~nleqreport()
-{
-}
-
-/*************************************************************************
- 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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::nleqcreatelm(n, m, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::nleqstate*>(state.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
- 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 m, const real_1d_array &x, nleqstate &state)
-{
- alglib_impl::ae_state _alglib_env_state;
- ae_int_t n;
-
- n = x.length();
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::nleqcreatelm(n, m, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::nleqstate*>(state.c_ptr()), &_alglib_env_state);
-
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::nleqsetcond(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), epsf, maxits, &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::nleqsetxrep(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), needxrep, &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::nleqsetstpmax(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), stpmax, &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- ae_bool result = alglib_impl::nleqiteration(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return *(reinterpret_cast<bool*>(&result));
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-
-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),
- void *ptr)
-{
- alglib_impl::ae_state _alglib_env_state;
- if( func==NULL )
- throw ap_error("ALGLIB: error in 'nleqsolve()' (func is NULL)");
- if( jac==NULL )
- throw ap_error("ALGLIB: error in 'nleqsolve()' (jac is NULL)");
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- while( alglib_impl::nleqiteration(state.c_ptr(), &_alglib_env_state) )
- {
- if( state.needf )
- {
- func(state.x, state.f, ptr);
- continue;
- }
- if( state.needfij )
- {
- jac(state.x, state.fi, state.j, ptr);
- continue;
- }
- if( state.xupdated )
- {
- if( rep!=NULL )
- rep(state.x, state.f, ptr);
- continue;
- }
- throw ap_error("ALGLIB: error in 'nleqsolve' (some derivatives were not provided?)");
- }
- alglib_impl::ae_state_clear(&_alglib_env_state);
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::nleqresults(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::nleqreport*>(rep.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::nleqresultsbuf(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::nleqreport*>(rep.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-
-/*************************************************************************
-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)
-{
- alglib_impl::ae_state _alglib_env_state;
- alglib_impl::ae_state_init(&_alglib_env_state);
- try
- {
- alglib_impl::nleqrestartfrom(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
- alglib_impl::ae_state_clear(&_alglib_env_state);
- return;
- }
- catch(alglib_impl::ae_error_type)
- {
- throw ap_error(_alglib_env_state.error_msg);
- }
- catch(...)
- {
- throw;
- }
-}
-}
-
-/////////////////////////////////////////////////////////////////////////
-//
-// THIS SECTION CONTAINS IMPLEMENTATION OF COMPUTATIONAL CORE
-//
-/////////////////////////////////////////////////////////////////////////
-namespace alglib_impl
-{
-static void densesolver_rmatrixlusolveinternal(/* Real */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- double scalea,
- ae_int_t n,
- /* Real */ ae_matrix* a,
- ae_bool havea,
- /* Real */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_matrix* x,
- ae_state *_state);
-static void densesolver_spdmatrixcholeskysolveinternal(/* Real */ ae_matrix* cha,
- double sqrtscalea,
- ae_int_t n,
- ae_bool isupper,
- /* Real */ ae_matrix* a,
- ae_bool havea,
- /* Real */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_matrix* x,
- ae_state *_state);
-static void densesolver_cmatrixlusolveinternal(/* Complex */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- double scalea,
- ae_int_t n,
- /* Complex */ ae_matrix* a,
- ae_bool havea,
- /* Complex */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_matrix* x,
- ae_state *_state);
-static void densesolver_hpdmatrixcholeskysolveinternal(/* Complex */ ae_matrix* cha,
- double sqrtscalea,
- ae_int_t n,
- ae_bool isupper,
- /* Complex */ ae_matrix* a,
- ae_bool havea,
- /* Complex */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_matrix* x,
- ae_state *_state);
-static ae_int_t densesolver_densesolverrfsmax(ae_int_t n,
- double r1,
- double rinf,
- ae_state *_state);
-static ae_int_t densesolver_densesolverrfsmaxv2(ae_int_t n,
- double r2,
- ae_state *_state);
-static void densesolver_rbasiclusolve(/* Real */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- double scalea,
- ae_int_t n,
- /* Real */ ae_vector* xb,
- /* Real */ ae_vector* tmp,
- ae_state *_state);
-static void densesolver_spdbasiccholeskysolve(/* Real */ ae_matrix* cha,
- double sqrtscalea,
- ae_int_t n,
- ae_bool isupper,
- /* Real */ ae_vector* xb,
- /* Real */ ae_vector* tmp,
- ae_state *_state);
-static void densesolver_cbasiclusolve(/* Complex */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- double scalea,
- ae_int_t n,
- /* Complex */ ae_vector* xb,
- /* Complex */ ae_vector* tmp,
- ae_state *_state);
-static void densesolver_hpdbasiccholeskysolve(/* Complex */ ae_matrix* cha,
- double sqrtscalea,
- ae_int_t n,
- ae_bool isupper,
- /* Complex */ ae_vector* xb,
- /* Complex */ ae_vector* tmp,
- ae_state *_state);
-
-
-static ae_int_t nleq_armijomaxfev = 20;
-static void nleq_clearrequestfields(nleqstate* state, ae_state *_state);
-static ae_bool nleq_increaselambda(double* lambdav,
- double* nu,
- double lambdaup,
- ae_state *_state);
-static void nleq_decreaselambda(double* lambdav,
- double* nu,
- double lambdadown,
- ae_state *_state);
-
-
-
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix bm;
- ae_matrix xm;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_vector_clear(x);
- ae_matrix_init(&bm, 0, 0, DT_REAL, _state, ae_true);
- ae_matrix_init(&xm, 0, 0, DT_REAL, _state, ae_true);
-
- if( n<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&bm, n, 1, _state);
- ae_v_move(&bm.ptr.pp_double[0][0], bm.stride, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
- rmatrixsolvem(a, n, &bm, 1, ae_true, info, rep, &xm, _state);
- ae_vector_set_length(x, n, _state);
- ae_v_move(&x->ptr.p_double[0], 1, &xm.ptr.pp_double[0][0], xm.stride, ae_v_len(0,n-1));
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix da;
- ae_matrix emptya;
- ae_vector p;
- double scalea;
- ae_int_t i;
- ae_int_t j;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_matrix_init(&da, 0, 0, DT_REAL, _state, ae_true);
- ae_matrix_init(&emptya, 0, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&p, 0, DT_INT, _state, ae_true);
-
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&da, n, n, _state);
-
- /*
- * 1. scale matrix, max(|A[i,j]|)
- * 2. factorize scaled matrix
- * 3. solve
- */
- scalea = 0;
- for(i=0; i<=n-1; i++)
- {
- for(j=0; j<=n-1; j++)
- {
- scalea = ae_maxreal(scalea, ae_fabs(a->ptr.pp_double[i][j], _state), _state);
- }
- }
- if( ae_fp_eq(scalea,0) )
- {
- scalea = 1;
- }
- scalea = 1/scalea;
- for(i=0; i<=n-1; i++)
- {
- ae_v_move(&da.ptr.pp_double[i][0], 1, &a->ptr.pp_double[i][0], 1, ae_v_len(0,n-1));
- }
- rmatrixlu(&da, n, n, &p, _state);
- if( rfs )
- {
- densesolver_rmatrixlusolveinternal(&da, &p, scalea, n, a, ae_true, b, m, info, rep, x, _state);
- }
- else
- {
- densesolver_rmatrixlusolveinternal(&da, &p, scalea, n, &emptya, ae_false, b, m, info, rep, x, _state);
- }
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix bm;
- ae_matrix xm;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_vector_clear(x);
- ae_matrix_init(&bm, 0, 0, DT_REAL, _state, ae_true);
- ae_matrix_init(&xm, 0, 0, DT_REAL, _state, ae_true);
-
- if( n<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&bm, n, 1, _state);
- ae_v_move(&bm.ptr.pp_double[0][0], bm.stride, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
- rmatrixlusolvem(lua, p, n, &bm, 1, info, rep, &xm, _state);
- ae_vector_set_length(x, n, _state);
- ae_v_move(&x->ptr.p_double[0], 1, &xm.ptr.pp_double[0][0], xm.stride, ae_v_len(0,n-1));
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix emptya;
- ae_int_t i;
- ae_int_t j;
- double scalea;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_matrix_init(&emptya, 0, 0, DT_REAL, _state, ae_true);
-
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
-
- /*
- * 1. scale matrix, max(|U[i,j]|)
- * we assume that LU is in its normal form, i.e. |L[i,j]|<=1
- * 2. solve
- */
- scalea = 0;
- for(i=0; i<=n-1; i++)
- {
- for(j=i; j<=n-1; j++)
- {
- scalea = ae_maxreal(scalea, ae_fabs(lua->ptr.pp_double[i][j], _state), _state);
- }
- }
- if( ae_fp_eq(scalea,0) )
- {
- scalea = 1;
- }
- scalea = 1/scalea;
- densesolver_rmatrixlusolveinternal(lua, p, scalea, n, &emptya, ae_false, b, m, info, rep, x, _state);
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix bm;
- ae_matrix xm;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_vector_clear(x);
- ae_matrix_init(&bm, 0, 0, DT_REAL, _state, ae_true);
- ae_matrix_init(&xm, 0, 0, DT_REAL, _state, ae_true);
-
- if( n<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&bm, n, 1, _state);
- ae_v_move(&bm.ptr.pp_double[0][0], bm.stride, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
- rmatrixmixedsolvem(a, lua, p, n, &bm, 1, info, rep, &xm, _state);
- ae_vector_set_length(x, n, _state);
- ae_v_move(&x->ptr.p_double[0], 1, &xm.ptr.pp_double[0][0], xm.stride, ae_v_len(0,n-1));
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- double scalea;
- ae_int_t i;
- ae_int_t j;
-
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
-
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- return;
- }
-
- /*
- * 1. scale matrix, max(|A[i,j]|)
- * 2. factorize scaled matrix
- * 3. solve
- */
- scalea = 0;
- for(i=0; i<=n-1; i++)
- {
- for(j=0; j<=n-1; j++)
- {
- scalea = ae_maxreal(scalea, ae_fabs(a->ptr.pp_double[i][j], _state), _state);
- }
- }
- if( ae_fp_eq(scalea,0) )
- {
- scalea = 1;
- }
- scalea = 1/scalea;
- densesolver_rmatrixlusolveinternal(lua, p, scalea, n, a, ae_true, b, m, info, rep, x, _state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix da;
- ae_matrix emptya;
- ae_vector p;
- double scalea;
- ae_int_t i;
- ae_int_t j;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_matrix_init(&da, 0, 0, DT_COMPLEX, _state, ae_true);
- ae_matrix_init(&emptya, 0, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&p, 0, DT_INT, _state, ae_true);
-
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&da, n, n, _state);
-
- /*
- * 1. scale matrix, max(|A[i,j]|)
- * 2. factorize scaled matrix
- * 3. solve
- */
- scalea = 0;
- for(i=0; i<=n-1; i++)
- {
- for(j=0; j<=n-1; j++)
- {
- scalea = ae_maxreal(scalea, ae_c_abs(a->ptr.pp_complex[i][j], _state), _state);
- }
- }
- if( ae_fp_eq(scalea,0) )
- {
- scalea = 1;
- }
- scalea = 1/scalea;
- for(i=0; i<=n-1; i++)
- {
- ae_v_cmove(&da.ptr.pp_complex[i][0], 1, &a->ptr.pp_complex[i][0], 1, "N", ae_v_len(0,n-1));
- }
- cmatrixlu(&da, n, n, &p, _state);
- if( rfs )
- {
- densesolver_cmatrixlusolveinternal(&da, &p, scalea, n, a, ae_true, b, m, info, rep, x, _state);
- }
- else
- {
- densesolver_cmatrixlusolveinternal(&da, &p, scalea, n, &emptya, ae_false, b, m, info, rep, x, _state);
- }
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix bm;
- ae_matrix xm;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_vector_clear(x);
- ae_matrix_init(&bm, 0, 0, DT_COMPLEX, _state, ae_true);
- ae_matrix_init(&xm, 0, 0, DT_COMPLEX, _state, ae_true);
-
- if( n<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&bm, n, 1, _state);
- ae_v_cmove(&bm.ptr.pp_complex[0][0], bm.stride, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
- cmatrixsolvem(a, n, &bm, 1, ae_true, info, rep, &xm, _state);
- ae_vector_set_length(x, n, _state);
- ae_v_cmove(&x->ptr.p_complex[0], 1, &xm.ptr.pp_complex[0][0], xm.stride, "N", ae_v_len(0,n-1));
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix emptya;
- ae_int_t i;
- ae_int_t j;
- double scalea;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_matrix_init(&emptya, 0, 0, DT_COMPLEX, _state, ae_true);
-
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
-
- /*
- * 1. scale matrix, max(|U[i,j]|)
- * we assume that LU is in its normal form, i.e. |L[i,j]|<=1
- * 2. solve
- */
- scalea = 0;
- for(i=0; i<=n-1; i++)
- {
- for(j=i; j<=n-1; j++)
- {
- scalea = ae_maxreal(scalea, ae_c_abs(lua->ptr.pp_complex[i][j], _state), _state);
- }
- }
- if( ae_fp_eq(scalea,0) )
- {
- scalea = 1;
- }
- scalea = 1/scalea;
- densesolver_cmatrixlusolveinternal(lua, p, scalea, n, &emptya, ae_false, b, m, info, rep, x, _state);
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix bm;
- ae_matrix xm;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_vector_clear(x);
- ae_matrix_init(&bm, 0, 0, DT_COMPLEX, _state, ae_true);
- ae_matrix_init(&xm, 0, 0, DT_COMPLEX, _state, ae_true);
-
- if( n<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&bm, n, 1, _state);
- ae_v_cmove(&bm.ptr.pp_complex[0][0], bm.stride, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
- cmatrixlusolvem(lua, p, n, &bm, 1, info, rep, &xm, _state);
- ae_vector_set_length(x, n, _state);
- ae_v_cmove(&x->ptr.p_complex[0], 1, &xm.ptr.pp_complex[0][0], xm.stride, "N", ae_v_len(0,n-1));
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- double scalea;
- ae_int_t i;
- ae_int_t j;
-
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
-
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- return;
- }
-
- /*
- * 1. scale matrix, max(|A[i,j]|)
- * 2. factorize scaled matrix
- * 3. solve
- */
- scalea = 0;
- for(i=0; i<=n-1; i++)
- {
- for(j=0; j<=n-1; j++)
- {
- scalea = ae_maxreal(scalea, ae_c_abs(a->ptr.pp_complex[i][j], _state), _state);
- }
- }
- if( ae_fp_eq(scalea,0) )
- {
- scalea = 1;
- }
- scalea = 1/scalea;
- densesolver_cmatrixlusolveinternal(lua, p, scalea, n, a, ae_true, b, m, info, rep, x, _state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix bm;
- ae_matrix xm;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_vector_clear(x);
- ae_matrix_init(&bm, 0, 0, DT_COMPLEX, _state, ae_true);
- ae_matrix_init(&xm, 0, 0, DT_COMPLEX, _state, ae_true);
-
- if( n<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&bm, n, 1, _state);
- ae_v_cmove(&bm.ptr.pp_complex[0][0], bm.stride, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
- cmatrixmixedsolvem(a, lua, p, n, &bm, 1, info, rep, &xm, _state);
- ae_vector_set_length(x, n, _state);
- ae_v_cmove(&x->ptr.p_complex[0], 1, &xm.ptr.pp_complex[0][0], xm.stride, "N", ae_v_len(0,n-1));
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix da;
- double sqrtscalea;
- ae_int_t i;
- ae_int_t j;
- ae_int_t j1;
- ae_int_t j2;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_matrix_init(&da, 0, 0, DT_REAL, _state, ae_true);
-
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&da, n, n, _state);
-
- /*
- * 1. scale matrix, max(|A[i,j]|)
- * 2. factorize scaled matrix
- * 3. solve
- */
- sqrtscalea = 0;
- for(i=0; i<=n-1; i++)
- {
- if( isupper )
- {
- j1 = i;
- j2 = n-1;
- }
- else
- {
- j1 = 0;
- j2 = i;
- }
- for(j=j1; j<=j2; j++)
- {
- sqrtscalea = ae_maxreal(sqrtscalea, ae_fabs(a->ptr.pp_double[i][j], _state), _state);
- }
- }
- if( ae_fp_eq(sqrtscalea,0) )
- {
- sqrtscalea = 1;
- }
- sqrtscalea = 1/sqrtscalea;
- sqrtscalea = ae_sqrt(sqrtscalea, _state);
- for(i=0; i<=n-1; i++)
- {
- if( isupper )
- {
- j1 = i;
- j2 = n-1;
- }
- else
- {
- j1 = 0;
- j2 = i;
- }
- ae_v_move(&da.ptr.pp_double[i][j1], 1, &a->ptr.pp_double[i][j1], 1, ae_v_len(j1,j2));
- }
- if( !spdmatrixcholesky(&da, n, isupper, _state) )
- {
- ae_matrix_set_length(x, n, m, _state);
- for(i=0; i<=n-1; i++)
- {
- for(j=0; j<=m-1; j++)
- {
- x->ptr.pp_double[i][j] = 0;
- }
- }
- rep->r1 = 0;
- rep->rinf = 0;
- *info = -3;
- ae_frame_leave(_state);
- return;
- }
- *info = 1;
- densesolver_spdmatrixcholeskysolveinternal(&da, sqrtscalea, n, isupper, a, ae_true, b, m, info, rep, x, _state);
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix bm;
- ae_matrix xm;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_vector_clear(x);
- ae_matrix_init(&bm, 0, 0, DT_REAL, _state, ae_true);
- ae_matrix_init(&xm, 0, 0, DT_REAL, _state, ae_true);
-
- if( n<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&bm, n, 1, _state);
- ae_v_move(&bm.ptr.pp_double[0][0], bm.stride, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
- spdmatrixsolvem(a, n, isupper, &bm, 1, info, rep, &xm, _state);
- ae_vector_set_length(x, n, _state);
- ae_v_move(&x->ptr.p_double[0], 1, &xm.ptr.pp_double[0][0], xm.stride, ae_v_len(0,n-1));
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix emptya;
- double sqrtscalea;
- ae_int_t i;
- ae_int_t j;
- ae_int_t j1;
- ae_int_t j2;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_matrix_init(&emptya, 0, 0, DT_REAL, _state, ae_true);
-
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
-
- /*
- * 1. scale matrix, max(|U[i,j]|)
- * 2. factorize scaled matrix
- * 3. solve
- */
- sqrtscalea = 0;
- for(i=0; i<=n-1; i++)
- {
- if( isupper )
- {
- j1 = i;
- j2 = n-1;
- }
- else
- {
- j1 = 0;
- j2 = i;
- }
- for(j=j1; j<=j2; j++)
- {
- sqrtscalea = ae_maxreal(sqrtscalea, ae_fabs(cha->ptr.pp_double[i][j], _state), _state);
- }
- }
- if( ae_fp_eq(sqrtscalea,0) )
- {
- sqrtscalea = 1;
- }
- sqrtscalea = 1/sqrtscalea;
- densesolver_spdmatrixcholeskysolveinternal(cha, sqrtscalea, n, isupper, &emptya, ae_false, b, m, info, rep, x, _state);
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix bm;
- ae_matrix xm;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_vector_clear(x);
- ae_matrix_init(&bm, 0, 0, DT_REAL, _state, ae_true);
- ae_matrix_init(&xm, 0, 0, DT_REAL, _state, ae_true);
-
- if( n<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&bm, n, 1, _state);
- ae_v_move(&bm.ptr.pp_double[0][0], bm.stride, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
- spdmatrixcholeskysolvem(cha, n, isupper, &bm, 1, info, rep, &xm, _state);
- ae_vector_set_length(x, n, _state);
- ae_v_move(&x->ptr.p_double[0], 1, &xm.ptr.pp_double[0][0], xm.stride, ae_v_len(0,n-1));
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix da;
- double sqrtscalea;
- ae_int_t i;
- ae_int_t j;
- ae_int_t j1;
- ae_int_t j2;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_matrix_init(&da, 0, 0, DT_COMPLEX, _state, ae_true);
-
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&da, n, n, _state);
-
- /*
- * 1. scale matrix, max(|A[i,j]|)
- * 2. factorize scaled matrix
- * 3. solve
- */
- sqrtscalea = 0;
- for(i=0; i<=n-1; i++)
- {
- if( isupper )
- {
- j1 = i;
- j2 = n-1;
- }
- else
- {
- j1 = 0;
- j2 = i;
- }
- for(j=j1; j<=j2; j++)
- {
- sqrtscalea = ae_maxreal(sqrtscalea, ae_c_abs(a->ptr.pp_complex[i][j], _state), _state);
- }
- }
- if( ae_fp_eq(sqrtscalea,0) )
- {
- sqrtscalea = 1;
- }
- sqrtscalea = 1/sqrtscalea;
- sqrtscalea = ae_sqrt(sqrtscalea, _state);
- for(i=0; i<=n-1; i++)
- {
- if( isupper )
- {
- j1 = i;
- j2 = n-1;
- }
- else
- {
- j1 = 0;
- j2 = i;
- }
- ae_v_cmove(&da.ptr.pp_complex[i][j1], 1, &a->ptr.pp_complex[i][j1], 1, "N", ae_v_len(j1,j2));
- }
- if( !hpdmatrixcholesky(&da, n, isupper, _state) )
- {
- ae_matrix_set_length(x, n, m, _state);
- for(i=0; i<=n-1; i++)
- {
- for(j=0; j<=m-1; j++)
- {
- x->ptr.pp_complex[i][j] = ae_complex_from_d(0);
- }
- }
- rep->r1 = 0;
- rep->rinf = 0;
- *info = -3;
- ae_frame_leave(_state);
- return;
- }
- *info = 1;
- densesolver_hpdmatrixcholeskysolveinternal(&da, sqrtscalea, n, isupper, a, ae_true, b, m, info, rep, x, _state);
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix bm;
- ae_matrix xm;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_vector_clear(x);
- ae_matrix_init(&bm, 0, 0, DT_COMPLEX, _state, ae_true);
- ae_matrix_init(&xm, 0, 0, DT_COMPLEX, _state, ae_true);
-
- if( n<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&bm, n, 1, _state);
- ae_v_cmove(&bm.ptr.pp_complex[0][0], bm.stride, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
- hpdmatrixsolvem(a, n, isupper, &bm, 1, info, rep, &xm, _state);
- ae_vector_set_length(x, n, _state);
- ae_v_cmove(&x->ptr.p_complex[0], 1, &xm.ptr.pp_complex[0][0], xm.stride, "N", ae_v_len(0,n-1));
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix emptya;
- double sqrtscalea;
- ae_int_t i;
- ae_int_t j;
- ae_int_t j1;
- ae_int_t j2;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_matrix_init(&emptya, 0, 0, DT_COMPLEX, _state, ae_true);
-
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
-
- /*
- * 1. scale matrix, max(|U[i,j]|)
- * 2. factorize scaled matrix
- * 3. solve
- */
- sqrtscalea = 0;
- for(i=0; i<=n-1; i++)
- {
- if( isupper )
- {
- j1 = i;
- j2 = n-1;
- }
- else
- {
- j1 = 0;
- j2 = i;
- }
- for(j=j1; j<=j2; j++)
- {
- sqrtscalea = ae_maxreal(sqrtscalea, ae_c_abs(cha->ptr.pp_complex[i][j], _state), _state);
- }
- }
- if( ae_fp_eq(sqrtscalea,0) )
- {
- sqrtscalea = 1;
- }
- sqrtscalea = 1/sqrtscalea;
- densesolver_hpdmatrixcholeskysolveinternal(cha, sqrtscalea, n, isupper, &emptya, ae_false, b, m, info, rep, x, _state);
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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)
-{
- ae_frame _frame_block;
- ae_matrix bm;
- ae_matrix xm;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_vector_clear(x);
- ae_matrix_init(&bm, 0, 0, DT_COMPLEX, _state, ae_true);
- ae_matrix_init(&xm, 0, 0, DT_COMPLEX, _state, ae_true);
-
- if( n<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(&bm, n, 1, _state);
- ae_v_cmove(&bm.ptr.pp_complex[0][0], bm.stride, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
- hpdmatrixcholeskysolvem(cha, n, isupper, &bm, 1, info, rep, &xm, _state);
- ae_vector_set_length(x, n, _state);
- ae_v_cmove(&x->ptr.p_complex[0], 1, &xm.ptr.pp_complex[0][0], xm.stride, "N", ae_v_len(0,n-1));
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-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(/* 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_frame _frame_block;
- ae_vector sv;
- ae_matrix u;
- ae_matrix vt;
- ae_vector rp;
- ae_vector utb;
- ae_vector sutb;
- ae_vector tmp;
- ae_vector ta;
- ae_vector tx;
- ae_vector buf;
- ae_vector w;
- ae_int_t i;
- ae_int_t j;
- ae_int_t nsv;
- ae_int_t kernelidx;
- double v;
- double verr;
- ae_bool svdfailed;
- ae_bool zeroa;
- ae_int_t rfs;
- ae_int_t nrfs;
- ae_bool terminatenexttime;
- ae_bool smallerr;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverlsreport_clear(rep);
- ae_vector_clear(x);
- ae_vector_init(&sv, 0, DT_REAL, _state, ae_true);
- ae_matrix_init(&u, 0, 0, DT_REAL, _state, ae_true);
- ae_matrix_init(&vt, 0, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&rp, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&utb, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&sutb, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&tmp, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&ta, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&tx, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&w, 0, DT_REAL, _state, ae_true);
-
- if( (nrows<=0||ncols<=0)||ae_fp_less(threshold,0) )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- if( ae_fp_eq(threshold,0) )
- {
- threshold = 1000*ae_machineepsilon;
- }
-
- /*
- * Factorize A first
- */
- svdfailed = !rmatrixsvd(a, nrows, ncols, 1, 2, 2, &sv, &u, &vt, _state);
- zeroa = ae_fp_eq(sv.ptr.p_double[0],0);
- if( svdfailed||zeroa )
- {
- if( svdfailed )
- {
- *info = -4;
- }
- else
- {
- *info = 1;
- }
- ae_vector_set_length(x, ncols, _state);
- for(i=0; i<=ncols-1; i++)
- {
- x->ptr.p_double[i] = 0;
- }
- rep->n = ncols;
- rep->k = ncols;
- ae_matrix_set_length(&rep->cx, ncols, ncols, _state);
- for(i=0; i<=ncols-1; i++)
- {
- for(j=0; j<=ncols-1; j++)
- {
- if( i==j )
- {
- rep->cx.ptr.pp_double[i][j] = 1;
- }
- else
- {
- rep->cx.ptr.pp_double[i][j] = 0;
- }
- }
- }
- rep->r2 = 0;
- ae_frame_leave(_state);
- return;
- }
- nsv = ae_minint(ncols, nrows, _state);
- if( nsv==ncols )
- {
- rep->r2 = sv.ptr.p_double[nsv-1]/sv.ptr.p_double[0];
- }
- else
- {
- rep->r2 = 0;
- }
- rep->n = ncols;
- *info = 1;
-
- /*
- * Iterative refinement of xc combined with solution:
- * 1. xc = 0
- * 2. calculate r = bc-A*xc using extra-precise dot product
- * 3. solve A*y = r
- * 4. update x:=x+r
- * 5. goto 2
- *
- * This cycle is executed until one of two things happens:
- * 1. maximum number of iterations reached
- * 2. last iteration decreased error to the lower limit
- */
- ae_vector_set_length(&utb, nsv, _state);
- ae_vector_set_length(&sutb, nsv, _state);
- ae_vector_set_length(x, ncols, _state);
- ae_vector_set_length(&tmp, ncols, _state);
- ae_vector_set_length(&ta, ncols+1, _state);
- ae_vector_set_length(&tx, ncols+1, _state);
- ae_vector_set_length(&buf, ncols+1, _state);
- for(i=0; i<=ncols-1; i++)
- {
- x->ptr.p_double[i] = 0;
- }
- kernelidx = nsv;
- for(i=0; i<=nsv-1; i++)
- {
- if( ae_fp_less_eq(sv.ptr.p_double[i],threshold*sv.ptr.p_double[0]) )
- {
- kernelidx = i;
- break;
- }
- }
- rep->k = ncols-kernelidx;
- nrfs = densesolver_densesolverrfsmaxv2(ncols, rep->r2, _state);
- terminatenexttime = ae_false;
- ae_vector_set_length(&rp, nrows, _state);
- for(rfs=0; rfs<=nrfs; rfs++)
- {
- if( terminatenexttime )
- {
- break;
- }
-
- /*
- * calculate right part
- */
- if( rfs==0 )
- {
- ae_v_move(&rp.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,nrows-1));
- }
- else
- {
- smallerr = ae_true;
- for(i=0; i<=nrows-1; i++)
- {
- ae_v_move(&ta.ptr.p_double[0], 1, &a->ptr.pp_double[i][0], 1, ae_v_len(0,ncols-1));
- ta.ptr.p_double[ncols] = -1;
- ae_v_move(&tx.ptr.p_double[0], 1, &x->ptr.p_double[0], 1, ae_v_len(0,ncols-1));
- tx.ptr.p_double[ncols] = b->ptr.p_double[i];
- xdot(&ta, &tx, ncols+1, &buf, &v, &verr, _state);
- rp.ptr.p_double[i] = -v;
- smallerr = smallerr&&ae_fp_less(ae_fabs(v, _state),4*verr);
- }
- if( smallerr )
- {
- terminatenexttime = ae_true;
- }
- }
-
- /*
- * solve A*dx = rp
- */
- for(i=0; i<=ncols-1; i++)
- {
- tmp.ptr.p_double[i] = 0;
- }
- for(i=0; i<=nsv-1; i++)
- {
- utb.ptr.p_double[i] = 0;
- }
- for(i=0; i<=nrows-1; i++)
- {
- v = rp.ptr.p_double[i];
- ae_v_addd(&utb.ptr.p_double[0], 1, &u.ptr.pp_double[i][0], 1, ae_v_len(0,nsv-1), v);
- }
- for(i=0; i<=nsv-1; i++)
- {
- if( i<kernelidx )
- {
- sutb.ptr.p_double[i] = utb.ptr.p_double[i]/sv.ptr.p_double[i];
- }
- else
- {
- sutb.ptr.p_double[i] = 0;
- }
- }
- for(i=0; i<=nsv-1; i++)
- {
- v = sutb.ptr.p_double[i];
- ae_v_addd(&tmp.ptr.p_double[0], 1, &vt.ptr.pp_double[i][0], 1, ae_v_len(0,ncols-1), v);
- }
-
- /*
- * update x: x:=x+dx
- */
- ae_v_add(&x->ptr.p_double[0], 1, &tmp.ptr.p_double[0], 1, ae_v_len(0,ncols-1));
- }
-
- /*
- * fill CX
- */
- if( rep->k>0 )
- {
- ae_matrix_set_length(&rep->cx, ncols, rep->k, _state);
- for(i=0; i<=rep->k-1; i++)
- {
- ae_v_move(&rep->cx.ptr.pp_double[0][i], rep->cx.stride, &vt.ptr.pp_double[kernelidx+i][0], 1, ae_v_len(0,ncols-1));
- }
- }
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-Internal LU solver
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-static void densesolver_rmatrixlusolveinternal(/* Real */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- double scalea,
- ae_int_t n,
- /* Real */ ae_matrix* a,
- ae_bool havea,
- /* Real */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_matrix* x,
- ae_state *_state)
-{
- ae_frame _frame_block;
- ae_int_t i;
- ae_int_t j;
- ae_int_t k;
- ae_int_t rfs;
- ae_int_t nrfs;
- ae_vector xc;
- ae_vector y;
- ae_vector bc;
- ae_vector xa;
- ae_vector xb;
- ae_vector tx;
- double v;
- double verr;
- double mxb;
- double scaleright;
- ae_bool smallerr;
- ae_bool terminatenexttime;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_vector_init(&xc, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&y, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&bc, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&xa, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&xb, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&tx, 0, DT_REAL, _state, ae_true);
-
- ae_assert(ae_fp_greater(scalea,0), "Assertion failed", _state);
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- for(i=0; i<=n-1; i++)
- {
- if( p->ptr.p_int[i]>n-1||p->ptr.p_int[i]<i )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- }
- ae_matrix_set_length(x, n, m, _state);
- ae_vector_set_length(&y, n, _state);
- ae_vector_set_length(&xc, n, _state);
- ae_vector_set_length(&bc, n, _state);
- ae_vector_set_length(&tx, n+1, _state);
- ae_vector_set_length(&xa, n+1, _state);
- ae_vector_set_length(&xb, n+1, _state);
-
- /*
- * estimate condition number, test for near singularity
- */
- rep->r1 = rmatrixlurcond1(lua, n, _state);
- rep->rinf = rmatrixlurcondinf(lua, n, _state);
- if( ae_fp_less(rep->r1,rcondthreshold(_state))||ae_fp_less(rep->rinf,rcondthreshold(_state)) )
- {
- for(i=0; i<=n-1; i++)
- {
- for(j=0; j<=m-1; j++)
- {
- x->ptr.pp_double[i][j] = 0;
- }
- }
- rep->r1 = 0;
- rep->rinf = 0;
- *info = -3;
- ae_frame_leave(_state);
- return;
- }
- *info = 1;
-
- /*
- * solve
- */
- for(k=0; k<=m-1; k++)
- {
-
- /*
- * copy B to contiguous storage
- */
- ae_v_move(&bc.ptr.p_double[0], 1, &b->ptr.pp_double[0][k], b->stride, ae_v_len(0,n-1));
-
- /*
- * Scale right part:
- * * MX stores max(|Bi|)
- * * ScaleRight stores actual scaling applied to B when solving systems
- * it is chosen to make |scaleRight*b| close to 1.
- */
- mxb = 0;
- for(i=0; i<=n-1; i++)
- {
- mxb = ae_maxreal(mxb, ae_fabs(bc.ptr.p_double[i], _state), _state);
- }
- if( ae_fp_eq(mxb,0) )
- {
- mxb = 1;
- }
- scaleright = 1/mxb;
-
- /*
- * First, non-iterative part of solution process.
- * We use separate code for this task because
- * XDot is quite slow and we want to save time.
- */
- ae_v_moved(&xc.ptr.p_double[0], 1, &bc.ptr.p_double[0], 1, ae_v_len(0,n-1), scaleright);
- densesolver_rbasiclusolve(lua, p, scalea, n, &xc, &tx, _state);
-
- /*
- * Iterative refinement of xc:
- * * calculate r = bc-A*xc using extra-precise dot product
- * * solve A*y = r
- * * update x:=x+r
- *
- * This cycle is executed until one of two things happens:
- * 1. maximum number of iterations reached
- * 2. last iteration decreased error to the lower limit
- */
- if( havea )
- {
- nrfs = densesolver_densesolverrfsmax(n, rep->r1, rep->rinf, _state);
- terminatenexttime = ae_false;
- for(rfs=0; rfs<=nrfs-1; rfs++)
- {
- if( terminatenexttime )
- {
- break;
- }
-
- /*
- * generate right part
- */
- smallerr = ae_true;
- ae_v_move(&xb.ptr.p_double[0], 1, &xc.ptr.p_double[0], 1, ae_v_len(0,n-1));
- for(i=0; i<=n-1; i++)
- {
- ae_v_moved(&xa.ptr.p_double[0], 1, &a->ptr.pp_double[i][0], 1, ae_v_len(0,n-1), scalea);
- xa.ptr.p_double[n] = -1;
- xb.ptr.p_double[n] = scaleright*bc.ptr.p_double[i];
- xdot(&xa, &xb, n+1, &tx, &v, &verr, _state);
- y.ptr.p_double[i] = -v;
- smallerr = smallerr&&ae_fp_less(ae_fabs(v, _state),4*verr);
- }
- if( smallerr )
- {
- terminatenexttime = ae_true;
- }
-
- /*
- * solve and update
- */
- densesolver_rbasiclusolve(lua, p, scalea, n, &y, &tx, _state);
- ae_v_add(&xc.ptr.p_double[0], 1, &y.ptr.p_double[0], 1, ae_v_len(0,n-1));
- }
- }
-
- /*
- * Store xc.
- * Post-scale result.
- */
- v = scalea*mxb;
- ae_v_moved(&x->ptr.pp_double[0][k], x->stride, &xc.ptr.p_double[0], 1, ae_v_len(0,n-1), v);
- }
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-Internal Cholesky solver
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-static void densesolver_spdmatrixcholeskysolveinternal(/* Real */ ae_matrix* cha,
- double sqrtscalea,
- ae_int_t n,
- ae_bool isupper,
- /* Real */ ae_matrix* a,
- ae_bool havea,
- /* Real */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Real */ ae_matrix* x,
- ae_state *_state)
-{
- ae_frame _frame_block;
- ae_int_t i;
- ae_int_t j;
- ae_int_t k;
- ae_vector xc;
- ae_vector y;
- ae_vector bc;
- ae_vector xa;
- ae_vector xb;
- ae_vector tx;
- double v;
- double mxb;
- double scaleright;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_vector_init(&xc, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&y, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&bc, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&xa, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&xb, 0, DT_REAL, _state, ae_true);
- ae_vector_init(&tx, 0, DT_REAL, _state, ae_true);
-
- ae_assert(ae_fp_greater(sqrtscalea,0), "Assertion failed", _state);
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(x, n, m, _state);
- ae_vector_set_length(&y, n, _state);
- ae_vector_set_length(&xc, n, _state);
- ae_vector_set_length(&bc, n, _state);
- ae_vector_set_length(&tx, n+1, _state);
- ae_vector_set_length(&xa, n+1, _state);
- ae_vector_set_length(&xb, n+1, _state);
-
- /*
- * estimate condition number, test for near singularity
- */
- rep->r1 = spdmatrixcholeskyrcond(cha, n, isupper, _state);
- rep->rinf = rep->r1;
- if( ae_fp_less(rep->r1,rcondthreshold(_state)) )
- {
- for(i=0; i<=n-1; i++)
- {
- for(j=0; j<=m-1; j++)
- {
- x->ptr.pp_double[i][j] = 0;
- }
- }
- rep->r1 = 0;
- rep->rinf = 0;
- *info = -3;
- ae_frame_leave(_state);
- return;
- }
- *info = 1;
-
- /*
- * solve
- */
- for(k=0; k<=m-1; k++)
- {
-
- /*
- * copy B to contiguous storage
- */
- ae_v_move(&bc.ptr.p_double[0], 1, &b->ptr.pp_double[0][k], b->stride, ae_v_len(0,n-1));
-
- /*
- * Scale right part:
- * * MX stores max(|Bi|)
- * * ScaleRight stores actual scaling applied to B when solving systems
- * it is chosen to make |scaleRight*b| close to 1.
- */
- mxb = 0;
- for(i=0; i<=n-1; i++)
- {
- mxb = ae_maxreal(mxb, ae_fabs(bc.ptr.p_double[i], _state), _state);
- }
- if( ae_fp_eq(mxb,0) )
- {
- mxb = 1;
- }
- scaleright = 1/mxb;
-
- /*
- * First, non-iterative part of solution process.
- * We use separate code for this task because
- * XDot is quite slow and we want to save time.
- */
- ae_v_moved(&xc.ptr.p_double[0], 1, &bc.ptr.p_double[0], 1, ae_v_len(0,n-1), scaleright);
- densesolver_spdbasiccholeskysolve(cha, sqrtscalea, n, isupper, &xc, &tx, _state);
-
- /*
- * Store xc.
- * Post-scale result.
- */
- v = ae_sqr(sqrtscalea, _state)*mxb;
- ae_v_moved(&x->ptr.pp_double[0][k], x->stride, &xc.ptr.p_double[0], 1, ae_v_len(0,n-1), v);
- }
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-Internal LU solver
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-static void densesolver_cmatrixlusolveinternal(/* Complex */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- double scalea,
- ae_int_t n,
- /* Complex */ ae_matrix* a,
- ae_bool havea,
- /* Complex */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_matrix* x,
- ae_state *_state)
-{
- ae_frame _frame_block;
- ae_int_t i;
- ae_int_t j;
- ae_int_t k;
- ae_int_t rfs;
- ae_int_t nrfs;
- ae_vector xc;
- ae_vector y;
- ae_vector bc;
- ae_vector xa;
- ae_vector xb;
- ae_vector tx;
- ae_vector tmpbuf;
- ae_complex v;
- double verr;
- double mxb;
- double scaleright;
- ae_bool smallerr;
- ae_bool terminatenexttime;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_vector_init(&xc, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&y, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&bc, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&xa, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&xb, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&tx, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&tmpbuf, 0, DT_REAL, _state, ae_true);
-
- ae_assert(ae_fp_greater(scalea,0), "Assertion failed", _state);
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- for(i=0; i<=n-1; i++)
- {
- if( p->ptr.p_int[i]>n-1||p->ptr.p_int[i]<i )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- }
- ae_matrix_set_length(x, n, m, _state);
- ae_vector_set_length(&y, n, _state);
- ae_vector_set_length(&xc, n, _state);
- ae_vector_set_length(&bc, n, _state);
- ae_vector_set_length(&tx, n, _state);
- ae_vector_set_length(&xa, n+1, _state);
- ae_vector_set_length(&xb, n+1, _state);
- ae_vector_set_length(&tmpbuf, 2*n+2, _state);
-
- /*
- * estimate condition number, test for near singularity
- */
- rep->r1 = cmatrixlurcond1(lua, n, _state);
- rep->rinf = cmatrixlurcondinf(lua, n, _state);
- if( ae_fp_less(rep->r1,rcondthreshold(_state))||ae_fp_less(rep->rinf,rcondthreshold(_state)) )
- {
- for(i=0; i<=n-1; i++)
- {
- for(j=0; j<=m-1; j++)
- {
- x->ptr.pp_complex[i][j] = ae_complex_from_d(0);
- }
- }
- rep->r1 = 0;
- rep->rinf = 0;
- *info = -3;
- ae_frame_leave(_state);
- return;
- }
- *info = 1;
-
- /*
- * solve
- */
- for(k=0; k<=m-1; k++)
- {
-
- /*
- * copy B to contiguous storage
- */
- ae_v_cmove(&bc.ptr.p_complex[0], 1, &b->ptr.pp_complex[0][k], b->stride, "N", ae_v_len(0,n-1));
-
- /*
- * Scale right part:
- * * MX stores max(|Bi|)
- * * ScaleRight stores actual scaling applied to B when solving systems
- * it is chosen to make |scaleRight*b| close to 1.
- */
- mxb = 0;
- for(i=0; i<=n-1; i++)
- {
- mxb = ae_maxreal(mxb, ae_c_abs(bc.ptr.p_complex[i], _state), _state);
- }
- if( ae_fp_eq(mxb,0) )
- {
- mxb = 1;
- }
- scaleright = 1/mxb;
-
- /*
- * First, non-iterative part of solution process.
- * We use separate code for this task because
- * XDot is quite slow and we want to save time.
- */
- ae_v_cmoved(&xc.ptr.p_complex[0], 1, &bc.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1), scaleright);
- densesolver_cbasiclusolve(lua, p, scalea, n, &xc, &tx, _state);
-
- /*
- * Iterative refinement of xc:
- * * calculate r = bc-A*xc using extra-precise dot product
- * * solve A*y = r
- * * update x:=x+r
- *
- * This cycle is executed until one of two things happens:
- * 1. maximum number of iterations reached
- * 2. last iteration decreased error to the lower limit
- */
- if( havea )
- {
- nrfs = densesolver_densesolverrfsmax(n, rep->r1, rep->rinf, _state);
- terminatenexttime = ae_false;
- for(rfs=0; rfs<=nrfs-1; rfs++)
- {
- if( terminatenexttime )
- {
- break;
- }
-
- /*
- * generate right part
- */
- smallerr = ae_true;
- ae_v_cmove(&xb.ptr.p_complex[0], 1, &xc.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
- for(i=0; i<=n-1; i++)
- {
- ae_v_cmoved(&xa.ptr.p_complex[0], 1, &a->ptr.pp_complex[i][0], 1, "N", ae_v_len(0,n-1), scalea);
- xa.ptr.p_complex[n] = ae_complex_from_d(-1);
- xb.ptr.p_complex[n] = ae_c_mul_d(bc.ptr.p_complex[i],scaleright);
- xcdot(&xa, &xb, n+1, &tmpbuf, &v, &verr, _state);
- y.ptr.p_complex[i] = ae_c_neg(v);
- smallerr = smallerr&&ae_fp_less(ae_c_abs(v, _state),4*verr);
- }
- if( smallerr )
- {
- terminatenexttime = ae_true;
- }
-
- /*
- * solve and update
- */
- densesolver_cbasiclusolve(lua, p, scalea, n, &y, &tx, _state);
- ae_v_cadd(&xc.ptr.p_complex[0], 1, &y.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
- }
- }
-
- /*
- * Store xc.
- * Post-scale result.
- */
- v = ae_complex_from_d(scalea*mxb);
- ae_v_cmovec(&x->ptr.pp_complex[0][k], x->stride, &xc.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1), v);
- }
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-Internal Cholesky solver
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-static void densesolver_hpdmatrixcholeskysolveinternal(/* Complex */ ae_matrix* cha,
- double sqrtscalea,
- ae_int_t n,
- ae_bool isupper,
- /* Complex */ ae_matrix* a,
- ae_bool havea,
- /* Complex */ ae_matrix* b,
- ae_int_t m,
- ae_int_t* info,
- densesolverreport* rep,
- /* Complex */ ae_matrix* x,
- ae_state *_state)
-{
- ae_frame _frame_block;
- ae_int_t i;
- ae_int_t j;
- ae_int_t k;
- ae_vector xc;
- ae_vector y;
- ae_vector bc;
- ae_vector xa;
- ae_vector xb;
- ae_vector tx;
- double v;
- double mxb;
- double scaleright;
-
- ae_frame_make(_state, &_frame_block);
- *info = 0;
- _densesolverreport_clear(rep);
- ae_matrix_clear(x);
- ae_vector_init(&xc, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&y, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&bc, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&xa, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&xb, 0, DT_COMPLEX, _state, ae_true);
- ae_vector_init(&tx, 0, DT_COMPLEX, _state, ae_true);
-
- ae_assert(ae_fp_greater(sqrtscalea,0), "Assertion failed", _state);
-
- /*
- * prepare: check inputs, allocate space...
- */
- if( n<=0||m<=0 )
- {
- *info = -1;
- ae_frame_leave(_state);
- return;
- }
- ae_matrix_set_length(x, n, m, _state);
- ae_vector_set_length(&y, n, _state);
- ae_vector_set_length(&xc, n, _state);
- ae_vector_set_length(&bc, n, _state);
- ae_vector_set_length(&tx, n+1, _state);
- ae_vector_set_length(&xa, n+1, _state);
- ae_vector_set_length(&xb, n+1, _state);
-
- /*
- * estimate condition number, test for near singularity
- */
- rep->r1 = hpdmatrixcholeskyrcond(cha, n, isupper, _state);
- rep->rinf = rep->r1;
- if( ae_fp_less(rep->r1,rcondthreshold(_state)) )
- {
- for(i=0; i<=n-1; i++)
- {
- for(j=0; j<=m-1; j++)
- {
- x->ptr.pp_complex[i][j] = ae_complex_from_d(0);
- }
- }
- rep->r1 = 0;
- rep->rinf = 0;
- *info = -3;
- ae_frame_leave(_state);
- return;
- }
- *info = 1;
-
- /*
- * solve
- */
- for(k=0; k<=m-1; k++)
- {
-
- /*
- * copy B to contiguous storage
- */
- ae_v_cmove(&bc.ptr.p_complex[0], 1, &b->ptr.pp_complex[0][k], b->stride, "N", ae_v_len(0,n-1));
-
- /*
- * Scale right part:
- * * MX stores max(|Bi|)
- * * ScaleRight stores actual scaling applied to B when solving systems
- * it is chosen to make |scaleRight*b| close to 1.
- */
- mxb = 0;
- for(i=0; i<=n-1; i++)
- {
- mxb = ae_maxreal(mxb, ae_c_abs(bc.ptr.p_complex[i], _state), _state);
- }
- if( ae_fp_eq(mxb,0) )
- {
- mxb = 1;
- }
- scaleright = 1/mxb;
-
- /*
- * First, non-iterative part of solution process.
- * We use separate code for this task because
- * XDot is quite slow and we want to save time.
- */
- ae_v_cmoved(&xc.ptr.p_complex[0], 1, &bc.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1), scaleright);
- densesolver_hpdbasiccholeskysolve(cha, sqrtscalea, n, isupper, &xc, &tx, _state);
-
- /*
- * Store xc.
- * Post-scale result.
- */
- v = ae_sqr(sqrtscalea, _state)*mxb;
- ae_v_cmoved(&x->ptr.pp_complex[0][k], x->stride, &xc.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1), v);
- }
- ae_frame_leave(_state);
-}
-
-
-/*************************************************************************
-Internal subroutine.
-Returns maximum count of RFS iterations as function of:
-1. machine epsilon
-2. task size.
-3. condition number
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-static ae_int_t densesolver_densesolverrfsmax(ae_int_t n,
- double r1,
- double rinf,
- ae_state *_state)
-{
- ae_int_t result;
-
-
- result = 5;
- return result;
-}
-
-
-/*************************************************************************
-Internal subroutine.
-Returns maximum count of RFS iterations as function of:
-1. machine epsilon
-2. task size.
-3. norm-2 condition number
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-static ae_int_t densesolver_densesolverrfsmaxv2(ae_int_t n,
- double r2,
- ae_state *_state)
-{
- ae_int_t result;
-
-
- result = densesolver_densesolverrfsmax(n, 0, 0, _state);
- return result;
-}
-
-
-/*************************************************************************
-Basic LU solver for ScaleA*PLU*x = y.
-
-This subroutine assumes that:
-* L is well-scaled, and it is U which needs scaling by ScaleA.
-* A=PLU is well-conditioned, so no zero divisions or overflow may occur
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-static void densesolver_rbasiclusolve(/* Real */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- double scalea,
- ae_int_t n,
- /* Real */ ae_vector* xb,
- /* Real */ ae_vector* tmp,
- ae_state *_state)
-{
- ae_int_t i;
- double v;
-
-
- for(i=0; i<=n-1; i++)
- {
- if( p->ptr.p_int[i]!=i )
- {
- v = xb->ptr.p_double[i];
- xb->ptr.p_double[i] = xb->ptr.p_double[p->ptr.p_int[i]];
- xb->ptr.p_double[p->ptr.p_int[i]] = v;
- }
- }
- for(i=1; i<=n-1; i++)
- {
- v = ae_v_dotproduct(&lua->ptr.pp_double[i][0], 1, &xb->ptr.p_double[0], 1, ae_v_len(0,i-1));
- xb->ptr.p_double[i] = xb->ptr.p_double[i]-v;
- }
- xb->ptr.p_double[n-1] = xb->ptr.p_double[n-1]/(scalea*lua->ptr.pp_double[n-1][n-1]);
- for(i=n-2; i>=0; i--)
- {
- ae_v_moved(&tmp->ptr.p_double[i+1], 1, &lua->ptr.pp_double[i][i+1], 1, ae_v_len(i+1,n-1), scalea);
- v = ae_v_dotproduct(&tmp->ptr.p_double[i+1], 1, &xb->ptr.p_double[i+1], 1, ae_v_len(i+1,n-1));
- xb->ptr.p_double[i] = (xb->ptr.p_double[i]-v)/(scalea*lua->ptr.pp_double[i][i]);
- }
-}
-
-
-/*************************************************************************
-Basic Cholesky solver for ScaleA*Cholesky(A)'*x = y.
-
-This subroutine assumes that:
-* A*ScaleA is well scaled
-* A is well-conditioned, so no zero divisions or overflow may occur
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-static void densesolver_spdbasiccholeskysolve(/* Real */ ae_matrix* cha,
- double sqrtscalea,
- ae_int_t n,
- ae_bool isupper,
- /* Real */ ae_vector* xb,
- /* Real */ ae_vector* tmp,
- ae_state *_state)
-{
- ae_int_t i;
- double v;
-
-
-
- /*
- * A = L*L' or A=U'*U
- */
- if( isupper )
- {
-
- /*
- * Solve U'*y=b first.
- */
- for(i=0; i<=n-1; i++)
- {
- xb->ptr.p_double[i] = xb->ptr.p_double[i]/(sqrtscalea*cha->ptr.pp_double[i][i]);
- if( i<n-1 )
- {
- v = xb->ptr.p_double[i];
- ae_v_moved(&tmp->ptr.p_double[i+1], 1, &cha->ptr.pp_double[i][i+1], 1, ae_v_len(i+1,n-1), sqrtscalea);
- ae_v_subd(&xb->ptr.p_double[i+1], 1, &tmp->ptr.p_double[i+1], 1, ae_v_len(i+1,n-1), v);
- }
- }
-
- /*
- * Solve U*x=y then.
- */
- for(i=n-1; i>=0; i--)
- {
- if( i<n-1 )
- {
- ae_v_moved(&tmp->ptr.p_double[i+1], 1, &cha->ptr.pp_double[i][i+1], 1, ae_v_len(i+1,n-1), sqrtscalea);
- v = ae_v_dotproduct(&tmp->ptr.p_double[i+1], 1, &xb->ptr.p_double[i+1], 1, ae_v_len(i+1,n-1));
- xb->ptr.p_double[i] = xb->ptr.p_double[i]-v;
- }
- xb->ptr.p_double[i] = xb->ptr.p_double[i]/(sqrtscalea*cha->ptr.pp_double[i][i]);
- }
- }
- else
- {
-
- /*
- * Solve L*y=b first
- */
- for(i=0; i<=n-1; i++)
- {
- if( i>0 )
- {
- ae_v_moved(&tmp->ptr.p_double[0], 1, &cha->ptr.pp_double[i][0], 1, ae_v_len(0,i-1), sqrtscalea);
- v = ae_v_dotproduct(&tmp->ptr.p_double[0], 1, &xb->ptr.p_double[0], 1, ae_v_len(0,i-1));
- xb->ptr.p_double[i] = xb->ptr.p_double[i]-v;
- }
- xb->ptr.p_double[i] = xb->ptr.p_double[i]/(sqrtscalea*cha->ptr.pp_double[i][i]);
- }
-
- /*
- * Solve L'*x=y then.
- */
- for(i=n-1; i>=0; i--)
- {
- xb->ptr.p_double[i] = xb->ptr.p_double[i]/(sqrtscalea*cha->ptr.pp_double[i][i]);
- if( i>0 )
- {
- v = xb->ptr.p_double[i];
- ae_v_moved(&tmp->ptr.p_double[0], 1, &cha->ptr.pp_double[i][0], 1, ae_v_len(0,i-1), sqrtscalea);
- ae_v_subd(&xb->ptr.p_double[0], 1, &tmp->ptr.p_double[0], 1, ae_v_len(0,i-1), v);
- }
- }
- }
-}
-
-
-/*************************************************************************
-Basic LU solver for ScaleA*PLU*x = y.
-
-This subroutine assumes that:
-* L is well-scaled, and it is U which needs scaling by ScaleA.
-* A=PLU is well-conditioned, so no zero divisions or overflow may occur
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-static void densesolver_cbasiclusolve(/* Complex */ ae_matrix* lua,
- /* Integer */ ae_vector* p,
- double scalea,
- ae_int_t n,
- /* Complex */ ae_vector* xb,
- /* Complex */ ae_vector* tmp,
- ae_state *_state)
-{
- ae_int_t i;
- ae_complex v;
-
-
- for(i=0; i<=n-1; i++)
- {
- if( p->ptr.p_int[i]!=i )
- {
- v = xb->ptr.p_complex[i];
- xb->ptr.p_complex[i] = xb->ptr.p_complex[p->ptr.p_int[i]];
- xb->ptr.p_complex[p->ptr.p_int[i]] = v;
- }
- }
- for(i=1; i<=n-1; i++)
- {
- v = ae_v_cdotproduct(&lua->ptr.pp_complex[i][0], 1, "N", &xb->ptr.p_complex[0], 1, "N", ae_v_len(0,i-1));
- xb->ptr.p_complex[i] = ae_c_sub(xb->ptr.p_complex[i],v);
- }
- xb->ptr.p_complex[n-1] = ae_c_div(xb->ptr.p_complex[n-1],ae_c_mul_d(lua->ptr.pp_complex[n-1][n-1],scalea));
- for(i=n-2; i>=0; i--)
- {
- ae_v_cmoved(&tmp->ptr.p_complex[i+1], 1, &lua->ptr.pp_complex[i][i+1], 1, "N", ae_v_len(i+1,n-1), scalea);
- v = ae_v_cdotproduct(&tmp->ptr.p_complex[i+1], 1, "N", &xb->ptr.p_complex[i+1], 1, "N", ae_v_len(i+1,n-1));
- xb->ptr.p_complex[i] = ae_c_div(ae_c_sub(xb->ptr.p_complex[i],v),ae_c_mul_d(lua->ptr.pp_complex[i][i],scalea));
- }
-}
-
-
-/*************************************************************************
-Basic Cholesky solver for ScaleA*Cholesky(A)'*x = y.
-
-This subroutine assumes that:
-* A*ScaleA is well scaled
-* A is well-conditioned, so no zero divisions or overflow may occur
-
- -- ALGLIB --
- Copyright 27.01.2010 by Bochkanov Sergey
-*************************************************************************/
-static void densesolver_hpdbasiccholeskysolve(/* Complex */ ae_matrix* cha,
- double sqrtscalea,
- ae_int_t n,
- ae_bool isupper,
- /* Complex */ ae_vector* xb,
- /* Complex */ ae_vector* tmp,
- ae_state *_state)
-{
- ae_int_t i;
- ae_complex v;
-
-
-
- /*
- * A = L*L' or A=U'*U
- */
- if( isupper )
- {
-
- /*
- * Solve U'*y=b first.
- */
- for(i=0; i<=n-1; i++)
- {
- xb->ptr.p_complex[i] = ae_c_div(xb->ptr.p_complex[i],ae_c_mul_d(ae_c_conj(cha->ptr.pp_complex[i][i], _state),sqrtscalea));
- if( i<n-1 )
- {
- v = xb->ptr.p_complex[i];
- ae_v_cmoved(&tmp->ptr.p_complex[i+1], 1, &cha->ptr.pp_complex[i][i+1], 1, "Conj", ae_v_len(i+1,n-1), sqrtscalea);
- ae_v_csubc(&xb->ptr.p_complex[i+1], 1, &tmp->ptr.p_complex[i+1], 1, "N", ae_v_len(i+1,n-1), v);
- }
- }
-
- /*
- * Solve U*x=y then.
- */
- for(i=n-1; i>=0; i--)
- {
- if( i<n-1 )
- {
- ae_v_cmoved(&tmp->ptr.p_complex[i+1], 1, &cha->ptr.pp_complex[i][i+1], 1, "N", ae_v_len(i+1,n-1), sqrtscalea);
- v = ae_v_cdotproduct(&tmp->ptr.p_complex[i+1], 1, "N", &xb->ptr.p_complex[i+1], 1, "N", ae_v_len(i+1,n-1));
- xb->ptr.p_complex[i] = ae_c_sub(xb->ptr.p_complex[i],v);
- }
- xb->ptr.p_complex[i] = ae_c_div(xb->ptr.p_complex[i],ae_c_mul_d(cha->ptr.pp_complex[i][i],sqrtscalea));
- }
- }
- else
- {
-
- /*
- * Solve L*y=b first
- */
- for(i=0; i<=n-1; i++)
- {
- if( i>0 )
- {
- ae_v_cmoved(&tmp->ptr.p_complex[0], 1, &cha->ptr.pp_complex[i][0], 1, "N", ae_v_len(0,i-1), sqrtscalea);
- v = ae_v_cdotproduct(&tmp->ptr.p_complex[0], 1, "N", &xb->ptr.p_complex[0], 1, "N", ae_v_len(0,i-1));
- xb->ptr.p_complex[i] = ae_c_sub(xb->ptr.p_complex[i],v);
- }
- xb->ptr.p_complex[i] = ae_c_div(xb->ptr.p_complex[i],ae_c_mul_d(cha->ptr.pp_complex[i][i],sqrtscalea));
- }
-
- /*
- * Solve L'*x=y then.
- */
- for(i=n-1; i>=0; i--)
- {
- xb->ptr.p_complex[i] = ae_c_div(xb->ptr.p_complex[i],ae_c_mul_d(ae_c_conj(cha->ptr.pp_complex[i][i], _state),sqrtscalea));
- if( i>0 )
- {
- v = xb->ptr.p_complex[i];
- ae_v_cmoved(&tmp->ptr.p_complex[0], 1, &cha->ptr.pp_complex[i][0], 1, "Conj", ae_v_len(0,i-1), sqrtscalea);
- ae_v_csubc(&xb->ptr.p_complex[0], 1, &tmp->ptr.p_complex[0], 1, "N", ae_v_len(0,i-1), v);
- }
- }
- }
-}
-
-
-ae_bool _densesolverreport_init(densesolverreport* p, ae_state *_state, ae_bool make_automatic)
-{
- return ae_true;
-}
-
-
-ae_bool _densesolverreport_init_copy(densesolverreport* dst, densesolverreport* src, ae_state *_state, ae_bool make_automatic)
-{
- dst->r1 = src->r1;
- dst->rinf = src->rinf;
- return ae_true;
-}
-
-
-void _densesolverreport_clear(densesolverreport* p)
-{
-}
-
-
-ae_bool _densesolverlsreport_init(densesolverlsreport* p, ae_state *_state, ae_bool make_automatic)
-{
- if( !ae_matrix_init(&p->cx, 0, 0, DT_REAL, _state, make_automatic) )
- return ae_false;
- return ae_true;
-}
-
-
-ae_bool _densesolverlsreport_init_copy(densesolverlsreport* dst, densesolverlsreport* src, ae_state *_state, ae_bool make_automatic)
-{
- dst->r2 = src->r2;
- if( !ae_matrix_init_copy(&dst->cx, &src->cx, _state, make_automatic) )
- return ae_false;
- dst->n = src->n;
- dst->k = src->k;
- return ae_true;
-}
-
-
-void _densesolverlsreport_clear(densesolverlsreport* p)
-{
- ae_matrix_clear(&p->cx);
-}
-
-
-
-
-/*************************************************************************
- 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(ae_int_t n,
- ae_int_t m,
- /* Real */ ae_vector* x,
- nleqstate* state,
- ae_state *_state)
-{
-
- _nleqstate_clear(state);
-
- ae_assert(n>=1, "NLEQCreateLM: N<1!", _state);
- ae_assert(m>=1, "NLEQCreateLM: M<1!", _state);
- ae_assert(x->cnt>=n, "NLEQCreateLM: Length(X)<N!", _state);
- ae_assert(isfinitevector(x, n, _state), "NLEQCreateLM: X contains infinite or NaN values!", _state);
-
- /*
- * Initialize
- */
- state->n = n;
- state->m = m;
- nleqsetcond(state, 0, 0, _state);
- nleqsetxrep(state, ae_false, _state);
- nleqsetstpmax(state, 0, _state);
- ae_vector_set_length(&state->x, n, _state);
- ae_vector_set_length(&state->xbase, n, _state);
- ae_matrix_set_length(&state->j, m, n, _state);
- ae_vector_set_length(&state->fi, m, _state);
- ae_vector_set_length(&state->rightpart, n, _state);
- ae_vector_set_length(&state->candstep, n, _state);
- nleqrestartfrom(state, x, _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(nleqstate* state,
- double epsf,
- ae_int_t maxits,
- ae_state *_state)
-{
-
-
- ae_assert(ae_isfinite(epsf, _state), "NLEQSetCond: EpsF is not finite number!", _state);
- ae_assert(ae_fp_greater_eq(epsf,0), "NLEQSetCond: negative EpsF!", _state);
- ae_assert(maxits>=0, "NLEQSetCond: negative MaxIts!", _state);
- if( ae_fp_eq(epsf,0)&&maxits==0 )
- {
- epsf = 1.0E-6;
- }
- state->epsf = epsf;
- state->maxits = 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(nleqstate* state, ae_bool needxrep, ae_state *_state)
-{
-
-
- state->xrep = 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(nleqstate* state, double stpmax, ae_state *_state)
-{
-
-
- ae_assert(ae_isfinite(stpmax, _state), "NLEQSetStpMax: StpMax is not finite!", _state);
- ae_assert(ae_fp_greater_eq(stpmax,0), "NLEQSetStpMax: StpMax<0!", _state);
- state->stpmax = stpmax;
-}
-
-
-/*************************************************************************
-
- -- ALGLIB --
- Copyright 20.03.2009 by Bochkanov Sergey
-*************************************************************************/
-ae_bool nleqiteration(nleqstate* state, ae_state *_state)
-{
- ae_int_t n;
- ae_int_t m;
- ae_int_t i;
- double lambdaup;
- double lambdadown;
- double lambdav;
- double rho;
- double mu;
- double stepnorm;
- ae_bool b;
- ae_bool result;
-
-
-
- /*
- * Reverse communication preparations
- * I know it looks ugly, but it works the same way
- * anywhere from C++ to Python.
- *
- * This code initializes locals by:
- * * random values determined during code
- * generation - on first subroutine call
- * * values from previous call - on subsequent calls
- */
- if( state->rstate.stage>=0 )
- {
- n = state->rstate.ia.ptr.p_int[0];
- m = state->rstate.ia.ptr.p_int[1];
- i = state->rstate.ia.ptr.p_int[2];
- b = state->rstate.ba.ptr.p_bool[0];
- lambdaup = state->rstate.ra.ptr.p_double[0];
- lambdadown = state->rstate.ra.ptr.p_double[1];
- lambdav = state->rstate.ra.ptr.p_double[2];
- rho = state->rstate.ra.ptr.p_double[3];
- mu = state->rstate.ra.ptr.p_double[4];
- stepnorm = state->rstate.ra.ptr.p_double[5];
- }
- else
- {
- n = -983;
- m = -989;
- i = -834;
- b = ae_false;
- lambdaup = -287;
- lambdadown = 364;
- lambdav = 214;
- rho = -338;
- mu = -686;
- stepnorm = 912;
- }
- if( state->rstate.stage==0 )
- {
- goto lbl_0;
- }
- if( state->rstate.stage==1 )
- {
- goto lbl_1;
- }
- if( state->rstate.stage==2 )
- {
- goto lbl_2;
- }
- if( state->rstate.stage==3 )
- {
- goto lbl_3;
- }
- if( state->rstate.stage==4 )
- {
- goto lbl_4;
- }
-
- /*
- * Routine body
- */
-
- /*
- * Prepare
- */
- n = state->n;
- m = state->m;
- state->repterminationtype = 0;
- state->repiterationscount = 0;
- state->repnfunc = 0;
- state->repnjac = 0;
-
- /*
- * Calculate F/G, initialize algorithm
- */
- nleq_clearrequestfields(state, _state);
- state->needf = ae_true;
- state->rstate.stage = 0;
- goto lbl_rcomm;
-lbl_0:
- state->needf = ae_false;
- state->repnfunc = state->repnfunc+1;
- ae_v_move(&state->xbase.ptr.p_double[0], 1, &state->x.ptr.p_double[0], 1, ae_v_len(0,n-1));
- state->fbase = state->f;
- state->fprev = ae_maxrealnumber;
- if( !state->xrep )
- {
- goto lbl_5;
- }
-
- /*
- * progress report
- */
- nleq_clearrequestfields(state, _state);
- state->xupdated = ae_true;
- state->rstate.stage = 1;
- goto lbl_rcomm;
-lbl_1:
- state->xupdated = ae_false;
-lbl_5:
- if( ae_fp_less_eq(state->f,ae_sqr(state->epsf, _state)) )
- {
- state->repterminationtype = 1;
- result = ae_false;
- return result;
- }
-
- /*
- * Main cycle
- */
- lambdaup = 10;
- lambdadown = 0.3;
- lambdav = 0.001;
- rho = 1;
-lbl_7:
- if( ae_false )
- {
- goto lbl_8;
- }
-
- /*
- * Get Jacobian;
- * before we get to this point we already have State.XBase filled
- * with current point and State.FBase filled with function value
- * at XBase
- */
- nleq_clearrequestfields(state, _state);
- state->needfij = ae_true;
- ae_v_move(&state->x.ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,n-1));
- state->rstate.stage = 2;
- goto lbl_rcomm;
-lbl_2:
- state->needfij = ae_false;
- state->repnfunc = state->repnfunc+1;
- state->repnjac = state->repnjac+1;
- rmatrixmv(n, m, &state->j, 0, 0, 1, &state->fi, 0, &state->rightpart, 0, _state);
- ae_v_muld(&state->rightpart.ptr.p_double[0], 1, ae_v_len(0,n-1), -1);
-
- /*
- * Inner cycle: find good lambda
- */
-lbl_9:
- if( ae_false )
- {
- goto lbl_10;
- }
-
- /*
- * Solve (J^T*J + (Lambda+Mu)*I)*y = J^T*F
- * to get step d=-y where:
- * * Mu=||F|| - is damping parameter for nonlinear system
- * * Lambda - is additional Levenberg-Marquardt parameter
- * for better convergence when far away from minimum
- */
- for(i=0; i<=n-1; i++)
- {
- state->candstep.ptr.p_double[i] = 0;
- }
- fblssolvecgx(&state->j, m, n, lambdav, &state->rightpart, &state->candstep, &state->cgbuf, _state);
-
- /*
- * Normalize step (it must be no more than StpMax)
- */
- stepnorm = 0;
- for(i=0; i<=n-1; i++)
- {
- if( ae_fp_neq(state->candstep.ptr.p_double[i],0) )
- {
- stepnorm = 1;
- break;
- }
- }
- linminnormalized(&state->candstep, &stepnorm, n, _state);
- if( ae_fp_neq(state->stpmax,0) )
- {
- stepnorm = ae_minreal(stepnorm, state->stpmax, _state);
- }
-
- /*
- * Test new step - is it good enough?
- * * if not, Lambda is increased and we try again.
- * * if step is good, we decrease Lambda and move on.
- *
- * We can break this cycle on two occasions:
- * * step is so small that x+step==x (in floating point arithmetics)
- * * lambda is so large
- */
- ae_v_move(&state->x.ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,n-1));
- ae_v_addd(&state->x.ptr.p_double[0], 1, &state->candstep.ptr.p_double[0], 1, ae_v_len(0,n-1), stepnorm);
- b = ae_true;
- for(i=0; i<=n-1; i++)
- {
- if( ae_fp_neq(state->x.ptr.p_double[i],state->xbase.ptr.p_double[i]) )
- {
- b = ae_false;
- break;
- }
- }
- if( b )
- {
-
- /*
- * Step is too small, force zero step and break
- */
- stepnorm = 0;
- ae_v_move(&state->x.ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,n-1));
- state->f = state->fbase;
- goto lbl_10;
- }
- nleq_clearrequestfields(state, _state);
- state->needf = ae_true;
- state->rstate.stage = 3;
- goto lbl_rcomm;
-lbl_3:
- state->needf = ae_false;
- state->repnfunc = state->repnfunc+1;
- if( ae_fp_less(state->f,state->fbase) )
- {
-
- /*
- * function value decreased, move on
- */
- nleq_decreaselambda(&lambdav, &rho, lambdadown, _state);
- goto lbl_10;
- }
- if( !nleq_increaselambda(&lambdav, &rho, lambdaup, _state) )
- {
-
- /*
- * Lambda is too large (near overflow), force zero step and break
- */
- stepnorm = 0;
- ae_v_move(&state->x.ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,n-1));
- state->f = state->fbase;
- goto lbl_10;
- }
- goto lbl_9;
-lbl_10:
-
- /*
- * Accept step:
- * * new position
- * * new function value
- */
- state->fbase = state->f;
- ae_v_addd(&state->xbase.ptr.p_double[0], 1, &state->candstep.ptr.p_double[0], 1, ae_v_len(0,n-1), stepnorm);
- state->repiterationscount = state->repiterationscount+1;
-
- /*
- * Report new iteration
- */
- if( !state->xrep )
- {
- goto lbl_11;
- }
- nleq_clearrequestfields(state, _state);
- state->xupdated = ae_true;
- state->f = state->fbase;
- ae_v_move(&state->x.ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,n-1));
- state->rstate.stage = 4;
- goto lbl_rcomm;
-lbl_4:
- state->xupdated = ae_false;
-lbl_11:
-
- /*
- * Test stopping conditions on F, step (zero/non-zero) and MaxIts;
- * If one of the conditions is met, RepTerminationType is changed.
- */
- if( ae_fp_less_eq(ae_sqrt(state->f, _state),state->epsf) )
- {
- state->repterminationtype = 1;
- }
- if( ae_fp_eq(stepnorm,0)&&state->repterminationtype==0 )
- {
- state->repterminationtype = -4;
- }
- if( state->repiterationscount>=state->maxits&&state->maxits>0 )
- {
- state->repterminationtype = 5;
- }
- if( state->repterminationtype!=0 )
- {
- goto lbl_8;
- }
-
- /*
- * Now, iteration is finally over
- */
- goto lbl_7;
-lbl_8:
- result = ae_false;
- return result;
-
- /*
- * Saving state
- */
-lbl_rcomm:
- result = ae_true;
- state->rstate.ia.ptr.p_int[0] = n;
- state->rstate.ia.ptr.p_int[1] = m;
- state->rstate.ia.ptr.p_int[2] = i;
- state->rstate.ba.ptr.p_bool[0] = b;
- state->rstate.ra.ptr.p_double[0] = lambdaup;
- state->rstate.ra.ptr.p_double[1] = lambdadown;
- state->rstate.ra.ptr.p_double[2] = lambdav;
- state->rstate.ra.ptr.p_double[3] = rho;
- state->rstate.ra.ptr.p_double[4] = mu;
- state->rstate.ra.ptr.p_double[5] = stepnorm;
- return result;
-}
-
-
-/*************************************************************************
-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(nleqstate* state,
- /* Real */ ae_vector* x,
- nleqreport* rep,
- ae_state *_state)
-{
-
- ae_vector_clear(x);
- _nleqreport_clear(rep);
-
- nleqresultsbuf(state, x, rep, _state);
-}
-
-
-/*************************************************************************
-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(nleqstate* state,
- /* Real */ ae_vector* x,
- nleqreport* rep,
- ae_state *_state)
-{
-
-
- if( x->cnt<state->n )
- {
- ae_vector_set_length(x, state->n, _state);
- }
- ae_v_move(&x->ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,state->n-1));
- rep->iterationscount = state->repiterationscount;
- rep->nfunc = state->repnfunc;
- rep->njac = state->repnjac;
- rep->terminationtype = state->repterminationtype;
-}
-
-
-/*************************************************************************
-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(nleqstate* state,
- /* Real */ ae_vector* x,
- ae_state *_state)
-{
-
-
- ae_assert(x->cnt>=state->n, "NLEQRestartFrom: Length(X)<N!", _state);
- ae_assert(isfinitevector(x, state->n, _state), "NLEQRestartFrom: X contains infinite or NaN values!", _state);
- ae_v_move(&state->x.ptr.p_double[0], 1, &x->ptr.p_double[0], 1, ae_v_len(0,state->n-1));
- ae_vector_set_length(&state->rstate.ia, 2+1, _state);
- ae_vector_set_length(&state->rstate.ba, 0+1, _state);
- ae_vector_set_length(&state->rstate.ra, 5+1, _state);
- state->rstate.stage = -1;
- nleq_clearrequestfields(state, _state);
-}
-
-
-/*************************************************************************
-Clears request fileds (to be sure that we don't forgot to clear something)
-*************************************************************************/
-static void nleq_clearrequestfields(nleqstate* state, ae_state *_state)
-{
-
-
- state->needf = ae_false;
- state->needfij = ae_false;
- state->xupdated = ae_false;
-}
-
-
-/*************************************************************************
-Increases lambda, returns False when there is a danger of overflow
-*************************************************************************/
-static ae_bool nleq_increaselambda(double* lambdav,
- double* nu,
- double lambdaup,
- ae_state *_state)
-{
- double lnlambda;
- double lnnu;
- double lnlambdaup;
- double lnmax;
- ae_bool result;
-
-
- result = ae_false;
- lnlambda = ae_log(*lambdav, _state);
- lnlambdaup = ae_log(lambdaup, _state);
- lnnu = ae_log(*nu, _state);
- lnmax = 0.5*ae_log(ae_maxrealnumber, _state);
- if( ae_fp_greater(lnlambda+lnlambdaup+lnnu,lnmax) )
- {
- return result;
- }
- if( ae_fp_greater(lnnu+ae_log(2, _state),lnmax) )
- {
- return result;
- }
- *lambdav = *lambdav*lambdaup*(*nu);
- *nu = *nu*2;
- result = ae_true;
- return result;
-}
-
-
-/*************************************************************************
-Decreases lambda, but leaves it unchanged when there is danger of underflow.
-*************************************************************************/
-static void nleq_decreaselambda(double* lambdav,
- double* nu,
- double lambdadown,
- ae_state *_state)
-{
-
-
- *nu = 1;
- if( ae_fp_less(ae_log(*lambdav, _state)+ae_log(lambdadown, _state),ae_log(ae_minrealnumber, _state)) )
- {
- *lambdav = ae_minrealnumber;
- }
- else
- {
- *lambdav = *lambdav*lambdadown;
- }
-}
-
-
-ae_bool _nleqstate_init(nleqstate* p, ae_state *_state, ae_bool make_automatic)
-{
- if( !ae_vector_init(&p->x, 0, DT_REAL, _state, make_automatic) )
- return ae_false;
- if( !ae_vector_init(&p->fi, 0, DT_REAL, _state, make_automatic) )
- return ae_false;
- if( !ae_matrix_init(&p->j, 0, 0, DT_REAL, _state, make_automatic) )
- return ae_false;
- if( !_rcommstate_init(&p->rstate, _state, make_automatic) )
- return ae_false;
- if( !ae_vector_init(&p->xbase, 0, DT_REAL, _state, make_automatic) )
- return ae_false;
- if( !ae_vector_init(&p->candstep, 0, DT_REAL, _state, make_automatic) )
- return ae_false;
- if( !ae_vector_init(&p->rightpart, 0, DT_REAL, _state, make_automatic) )
- return ae_false;
- if( !ae_vector_init(&p->cgbuf, 0, DT_REAL, _state, make_automatic) )
- return ae_false;
- return ae_true;
-}
-
-
-ae_bool _nleqstate_init_copy(nleqstate* dst, nleqstate* src, ae_state *_state, ae_bool make_automatic)
-{
- dst->n = src->n;
- dst->m = src->m;
- dst->epsf = src->epsf;
- dst->maxits = src->maxits;
- dst->xrep = src->xrep;
- dst->stpmax = src->stpmax;
- if( !ae_vector_init_copy(&dst->x, &src->x, _state, make_automatic) )
- return ae_false;
- dst->f = src->f;
- if( !ae_vector_init_copy(&dst->fi, &src->fi, _state, make_automatic) )
- return ae_false;
- if( !ae_matrix_init_copy(&dst->j, &src->j, _state, make_automatic) )
- return ae_false;
- dst->needf = src->needf;
- dst->needfij = src->needfij;
- dst->xupdated = src->xupdated;
- if( !_rcommstate_init_copy(&dst->rstate, &src->rstate, _state, make_automatic) )
- return ae_false;
- dst->repiterationscount = src->repiterationscount;
- dst->repnfunc = src->repnfunc;
- dst->repnjac = src->repnjac;
- dst->repterminationtype = src->repterminationtype;
- if( !ae_vector_init_copy(&dst->xbase, &src->xbase, _state, make_automatic) )
- return ae_false;
- dst->fbase = src->fbase;
- dst->fprev = src->fprev;
- if( !ae_vector_init_copy(&dst->candstep, &src->candstep, _state, make_automatic) )
- return ae_false;
- if( !ae_vector_init_copy(&dst->rightpart, &src->rightpart, _state, make_automatic) )
- return ae_false;
- if( !ae_vector_init_copy(&dst->cgbuf, &src->cgbuf, _state, make_automatic) )
- return ae_false;
- return ae_true;
-}
-
-
-void _nleqstate_clear(nleqstate* p)
-{
- ae_vector_clear(&p->x);
- ae_vector_clear(&p->fi);
- ae_matrix_clear(&p->j);
- _rcommstate_clear(&p->rstate);
- ae_vector_clear(&p->xbase);
- ae_vector_clear(&p->candstep);
- ae_vector_clear(&p->rightpart);
- ae_vector_clear(&p->cgbuf);
-}
-
-
-ae_bool _nleqreport_init(nleqreport* p, ae_state *_state, ae_bool make_automatic)
-{
- return ae_true;
-}
-
-
-ae_bool _nleqreport_init_copy(nleqreport* dst, nleqreport* src, ae_state *_state, ae_bool make_automatic)
-{
- dst->iterationscount = src->iterationscount;
- dst->nfunc = src->nfunc;
- dst->njac = src->njac;
- dst->terminationtype = src->terminationtype;
- return ae_true;
-}
-
-
-void _nleqreport_clear(nleqreport* p)
-{
-}
-
-
-
-}
-