Skip to content
Snippets Groups Projects
Commit ecc979a0 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

Revert "3D lloyd et rtree" 

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