| dataanalysis.h | | dataanalysis.h | |
|---|
| | | | |
| skipping to change at line 25 | | skipping to change at line 25 | |
|---|
| A copy of the GNU General Public License is available at | | A copy of the GNU General Public License is available at | |
| http://www.fsf.org/licensing/licenses | | http://www.fsf.org/licensing/licenses | |
| >>> END OF LICENSE >>> | | >>> END OF LICENSE >>> | |
| *************************************************************************/ | | *************************************************************************/ | |
| #ifndef _dataanalysis_pkg_h | | #ifndef _dataanalysis_pkg_h | |
| #define _dataanalysis_pkg_h | | #define _dataanalysis_pkg_h | |
| #include "ap.h" | | #include "ap.h" | |
| #include "alglibinternal.h" | | #include "alglibinternal.h" | |
| #include "linalg.h" | | #include "linalg.h" | |
| #include "statistics.h" | | #include "statistics.h" | |
|
| #include "alglibmisc.h" | | | |
| #include "specialfunctions.h" | | #include "specialfunctions.h" | |
|
| | | #include "alglibmisc.h" | |
| #include "solvers.h" | | #include "solvers.h" | |
| #include "optimization.h" | | #include "optimization.h" | |
| | | | |
| ///////////////////////////////////////////////////////////////////////// | | ///////////////////////////////////////////////////////////////////////// | |
| // | | // | |
| // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES) | | // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES) | |
| // | | // | |
| ///////////////////////////////////////////////////////////////////////// | | ///////////////////////////////////////////////////////////////////////// | |
| namespace alglib_impl | | namespace alglib_impl | |
| { | | { | |
| | | | |
| skipping to change at line 207 | | skipping to change at line 207 | |
|---|
| typedef struct | | typedef struct | |
| { | | { | |
| double relclserror; | | double relclserror; | |
| double avgce; | | double avgce; | |
| double rmserror; | | double rmserror; | |
| double avgerror; | | double avgerror; | |
| double avgrelerror; | | double avgrelerror; | |
| } mlpcvreport; | | } mlpcvreport; | |
| typedef struct | | typedef struct | |
| { | | { | |
|
| ae_vector structinfo; | | | |
| ae_int_t ensemblesize; | | ae_int_t ensemblesize; | |
|
| ae_int_t nin; | | | |
| ae_int_t nout; | | | |
| ae_int_t wcount; | | | |
| ae_bool issoftmax; | | | |
| ae_bool postprocessing; | | | |
| ae_vector weights; | | ae_vector weights; | |
| ae_vector columnmeans; | | ae_vector columnmeans; | |
| ae_vector columnsigmas; | | ae_vector columnsigmas; | |
|
| ae_int_t serializedlen; | | multilayerperceptron network; | |
| ae_vector serializedmlp; | | | |
| ae_vector tmpweights; | | | |
| ae_vector tmpmeans; | | | |
| ae_vector tmpsigmas; | | | |
| ae_vector neurons; | | | |
| ae_vector dfdnet; | | | |
| ae_vector y; | | ae_vector y; | |
| } mlpensemble; | | } mlpensemble; | |
| | | | |
| } | | } | |
| | | | |
| ///////////////////////////////////////////////////////////////////////// | | ///////////////////////////////////////////////////////////////////////// | |
| // | | // | |
| // THIS SECTION CONTAINS C++ INTERFACE | | // THIS SECTION CONTAINS C++ INTERFACE | |
| // | | // | |
| ///////////////////////////////////////////////////////////////////////// | | ///////////////////////////////////////////////////////////////////////// | |
| | | | |
| skipping to change at line 882 | | skipping to change at line 870 | |
|---|
| Its meaning for regression task is obvious. As for | | Its meaning for regression task is obvious. As for | |
| classification task, it means average relative error when estimating | | classification task, it means average relative error when estimating | |
| posterior probability of belonging to the correct class. | | posterior probability of belonging to the correct class. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 16.02.2009 by Bochkanov Sergey | | Copyright 16.02.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| double dfavgrelerror(const decisionforest &df, const real_2d_array &xy, con
st ae_int_t npoints); | | double dfavgrelerror(const decisionforest &df, const real_2d_array &xy, con
st ae_int_t npoints); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| k-means++ clusterization | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| XY - dataset, array [0..NPoints-1,0..NVars-1]. | | | |
| NPoints - dataset size, NPoints>=K | | | |
| NVars - number of variables, NVars>=1 | | | |
| K - desired number of clusters, K>=1 | | | |
| Restarts - number of restarts, Restarts>=1 | | | |
| | | | |
| OUTPUT PARAMETERS: | | | |
| Info - return code: | | | |
| * -3, if task is degenerate (number of distinct points | | | |
| is | | | |
| less than K) | | | |
| * -1, if incorrect NPoints/NFeatures/K/Restarts was pas | | | |
| sed | | | |
| * 1, if subroutine finished successfully | | | |
| C - array[0..NVars-1,0..K-1].matrix whose columns store | | | |
| cluster's centers | | | |
| XYC - array[NPoints], which contains cluster indexes | | | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 21.03.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void kmeansgenerate(const real_2d_array &xy, const ae_int_t npoints, const | | | |
| ae_int_t nvars, const ae_int_t k, const ae_int_t restarts, ae_int_t &info, | | | |
| real_2d_array &c, integer_1d_array &xyc); | | | |
| | | | |
| /************************************************************************* | | | |
| Multiclass Fisher LDA | | | |
| | | | |
| Subroutine finds coefficients of linear combination which optimally separat | | | |
| es | | | |
| training set on classes. | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| XY - training set, array[0..NPoints-1,0..NVars]. | | | |
| First NVars columns store values of independent | | | |
| variables, next column stores number of class (from 0 | | | |
| to NClasses-1) which dataset element belongs to. Fracti | | | |
| onal | | | |
| values are rounded to nearest integer. | | | |
| NPoints - training set size, NPoints>=0 | | | |
| NVars - number of independent variables, NVars>=1 | | | |
| NClasses - number of classes, NClasses>=2 | | | |
| | | | |
| OUTPUT PARAMETERS: | | | |
| Info - return code: | | | |
| * -4, if internal EVD subroutine hasn't converged | | | |
| * -2, if there is a point with class number | | | |
| outside of [0..NClasses-1]. | | | |
| * -1, if incorrect parameters was passed (NPoints<0, | | | |
| NVars<1, NClasses<2) | | | |
| * 1, if task has been solved | | | |
| * 2, if there was a multicollinearity in training set, | | | |
| but task has been solved. | | | |
| W - linear combination coefficients, array[0..NVars-1] | | | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 31.05.2008 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void fisherlda(const real_2d_array &xy, const ae_int_t npoints, const ae_in | | | |
| t_t nvars, const ae_int_t nclasses, ae_int_t &info, real_1d_array &w); | | | |
| | | | |
| /************************************************************************* | | | |
| N-dimensional multiclass Fisher LDA | | | |
| | | | |
| Subroutine finds coefficients of linear combinations which optimally separa | | | |
| tes | | | |
| training set on classes. It returns N-dimensional basis whose vector are so | | | |
| rted | | | |
| by quality of training set separation (in descending order). | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| XY - training set, array[0..NPoints-1,0..NVars]. | | | |
| First NVars columns store values of independent | | | |
| variables, next column stores number of class (from 0 | | | |
| to NClasses-1) which dataset element belongs to. Fracti | | | |
| onal | | | |
| values are rounded to nearest integer. | | | |
| NPoints - training set size, NPoints>=0 | | | |
| NVars - number of independent variables, NVars>=1 | | | |
| NClasses - number of classes, NClasses>=2 | | | |
| | | | |
| OUTPUT PARAMETERS: | | | |
| Info - return code: | | | |
| * -4, if internal EVD subroutine hasn't converged | | | |
| * -2, if there is a point with class number | | | |
| outside of [0..NClasses-1]. | | | |
| * -1, if incorrect parameters was passed (NPoints<0, | | | |
| NVars<1, NClasses<2) | | | |
| * 1, if task has been solved | | | |
| * 2, if there was a multicollinearity in training set, | | | |
| but task has been solved. | | | |
| W - basis, array[0..NVars-1,0..NVars-1] | | | |
| columns of matrix stores basis vectors, sorted by | | | |
| quality of training set separation (in descending order | | | |
| ) | | | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 31.05.2008 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void fisherldan(const real_2d_array &xy, const ae_int_t npoints, const ae_i | | | |
| nt_t nvars, const ae_int_t nclasses, ae_int_t &info, real_2d_array &w); | | | |
| | | | |
| /************************************************************************* | | | |
| Linear regression | | Linear regression | |
| | | | |
| Subroutine builds model: | | Subroutine builds model: | |
| | | | |
| Y = A(0)*X[0] + ... + A(N-1)*X[N-1] + A(N) | | Y = A(0)*X[0] + ... + A(N-1)*X[N-1] + A(N) | |
| | | | |
| and model found in ALGLIB format, covariation matrix, training set errors | | and model found in ALGLIB format, covariation matrix, training set errors | |
| (rms, average, average relative) and leave-one-out cross-validation | | (rms, average, average relative) and leave-one-out cross-validation | |
| estimate of the generalization error. CV estimate calculated using fast | | estimate of the generalization error. CV estimate calculated using fast | |
| algorithm with O(NPoints*NVars) complexity. | | algorithm with O(NPoints*NVars) complexity. | |
| | | | |
| skipping to change at line 1164 | | skipping to change at line 1058 | |
|---|
| | | | |
| RESULT: | | RESULT: | |
| average relative error. | | average relative error. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 30.08.2008 by Bochkanov Sergey | | Copyright 30.08.2008 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| double lravgrelerror(const linearmodel &lm, const real_2d_array &xy, const
ae_int_t npoints); | | double lravgrelerror(const linearmodel &lm, const real_2d_array &xy, const
ae_int_t npoints); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | Filters: simple moving averages (unsymmetric). | |
| | | | |
| | | This filter replaces array by results of SMA(K) filter. SMA(K) is defined | |
| | | as filter which averages at most K previous points (previous - not points | |
| | | AROUND central point) - or less, in case of the first K-1 points. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - array[N], array to process. It can be larger than N, | |
| | | in this case only first N points are processed. | |
| | | N - points count, N>=0 | |
| | | K - K>=1 (K can be larger than N , such cases will be | |
| | | correctly handled). Window width. K=1 corresponds to | |
| | | identity transformation (nothing changes). | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | X - array, whose first N elements were processed with SMA(K | |
| | | ) | |
| | | | |
| | | NOTE 1: this function uses efficient in-place algorithm which does not | |
| | | allocate temporary arrays. | |
| | | | |
| | | NOTE 2: this algorithm makes only one pass through array and uses running | |
| | | sum to speed-up calculation of the averages. Additional measures | |
| | | are taken to ensure that running sum on a long sequence of zero | |
| | | elements will be correctly reset to zero even in the presence of | |
| | | round-off error. | |
| | | | |
| | | NOTE 3: this is unsymmetric version of the algorithm, which does NOT | |
| | | averages points after the current one. Only X[i], X[i-1], ... are | |
| | | used when calculating new value of X[i]. We should also note that | |
| | | this algorithm uses BOTH previous points and current one, i.e. | |
| | | new value of X[i] depends on BOTH previous point and X[i] itself. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 25.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void filtersma(real_1d_array &x, const ae_int_t n, const ae_int_t k); | |
| | | void filtersma(real_1d_array &x, const ae_int_t k); | |
| | | | |
| | | /************************************************************************* | |
| | | Filters: exponential moving averages. | |
| | | | |
| | | This filter replaces array by results of EMA(alpha) filter. EMA(alpha) is | |
| | | defined as filter which replaces X[] by S[]: | |
| | | S[0] = X[0] | |
| | | S[t] = alpha*X[t] + (1-alpha)*S[t-1] | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - array[N], array to process. It can be larger than N, | |
| | | in this case only first N points are processed. | |
| | | N - points count, N>=0 | |
| | | alpha - 0<alpha<=1, smoothing parameter. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | X - array, whose first N elements were processed | |
| | | with EMA(alpha) | |
| | | | |
| | | NOTE 1: this function uses efficient in-place algorithm which does not | |
| | | allocate temporary arrays. | |
| | | | |
| | | NOTE 2: this algorithm uses BOTH previous points and current one, i.e. | |
| | | new value of X[i] depends on BOTH previous point and X[i] itself. | |
| | | | |
| | | NOTE 3: technical analytis users quite often work with EMA coefficient | |
| | | expressed in DAYS instead of fractions. If you want to calculate | |
| | | EMA(N), where N is a number of days, you can use alpha=2/(N+1). | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 25.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void filterema(real_1d_array &x, const ae_int_t n, const double alpha); | |
| | | void filterema(real_1d_array &x, const double alpha); | |
| | | | |
| | | /************************************************************************* | |
| | | Filters: linear regression moving averages. | |
| | | | |
| | | This filter replaces array by results of LRMA(K) filter. | |
| | | | |
| | | LRMA(K) is defined as filter which, for each data point, builds linear | |
| | | regression model using K prevous points (point itself is included in | |
| | | these K points) and calculates value of this linear model at the point in | |
| | | question. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - array[N], array to process. It can be larger than N, | |
| | | in this case only first N points are processed. | |
| | | N - points count, N>=0 | |
| | | K - K>=1 (K can be larger than N , such cases will be | |
| | | correctly handled). Window width. K=1 corresponds to | |
| | | identity transformation (nothing changes). | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | X - array, whose first N elements were processed with SMA(K | |
| | | ) | |
| | | | |
| | | NOTE 1: this function uses efficient in-place algorithm which does not | |
| | | allocate temporary arrays. | |
| | | | |
| | | NOTE 2: this algorithm makes only one pass through array and uses running | |
| | | sum to speed-up calculation of the averages. Additional measures | |
| | | are taken to ensure that running sum on a long sequence of zero | |
| | | elements will be correctly reset to zero even in the presence of | |
| | | round-off error. | |
| | | | |
| | | NOTE 3: this is unsymmetric version of the algorithm, which does NOT | |
| | | averages points after the current one. Only X[i], X[i-1], ... are | |
| | | used when calculating new value of X[i]. We should also note that | |
| | | this algorithm uses BOTH previous points and current one, i.e. | |
| | | new value of X[i] depends on BOTH previous point and X[i] itself. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 25.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void filterlrma(real_1d_array &x, const ae_int_t n, const ae_int_t k); | |
| | | void filterlrma(real_1d_array &x, const ae_int_t k); | |
| | | | |
| | | /************************************************************************* | |
| | | k-means++ clusterization | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | XY - dataset, array [0..NPoints-1,0..NVars-1]. | |
| | | NPoints - dataset size, NPoints>=K | |
| | | NVars - number of variables, NVars>=1 | |
| | | K - desired number of clusters, K>=1 | |
| | | Restarts - number of restarts, Restarts>=1 | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info - return code: | |
| | | * -3, if task is degenerate (number of distinct points | |
| | | is | |
| | | less than K) | |
| | | * -1, if incorrect NPoints/NFeatures/K/Restarts was pas | |
| | | sed | |
| | | * 1, if subroutine finished successfully | |
| | | C - array[0..NVars-1,0..K-1].matrix whose columns store | |
| | | cluster's centers | |
| | | XYC - array[NPoints], which contains cluster indexes | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 21.03.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void kmeansgenerate(const real_2d_array &xy, const ae_int_t npoints, const | |
| | | ae_int_t nvars, const ae_int_t k, const ae_int_t restarts, ae_int_t &info, | |
| | | real_2d_array &c, integer_1d_array &xyc); | |
| | | | |
| | | /************************************************************************* | |
| | | Multiclass Fisher LDA | |
| | | | |
| | | Subroutine finds coefficients of linear combination which optimally separat | |
| | | es | |
| | | training set on classes. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | XY - training set, array[0..NPoints-1,0..NVars]. | |
| | | First NVars columns store values of independent | |
| | | variables, next column stores number of class (from 0 | |
| | | to NClasses-1) which dataset element belongs to. Fracti | |
| | | onal | |
| | | values are rounded to nearest integer. | |
| | | NPoints - training set size, NPoints>=0 | |
| | | NVars - number of independent variables, NVars>=1 | |
| | | NClasses - number of classes, NClasses>=2 | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info - return code: | |
| | | * -4, if internal EVD subroutine hasn't converged | |
| | | * -2, if there is a point with class number | |
| | | outside of [0..NClasses-1]. | |
| | | * -1, if incorrect parameters was passed (NPoints<0, | |
| | | NVars<1, NClasses<2) | |
| | | * 1, if task has been solved | |
| | | * 2, if there was a multicollinearity in training set, | |
| | | but task has been solved. | |
| | | W - linear combination coefficients, array[0..NVars-1] | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 31.05.2008 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void fisherlda(const real_2d_array &xy, const ae_int_t npoints, const ae_in | |
| | | t_t nvars, const ae_int_t nclasses, ae_int_t &info, real_1d_array &w); | |
| | | | |
| | | /************************************************************************* | |
| | | N-dimensional multiclass Fisher LDA | |
| | | | |
| | | Subroutine finds coefficients of linear combinations which optimally separa | |
| | | tes | |
| | | training set on classes. It returns N-dimensional basis whose vector are so | |
| | | rted | |
| | | by quality of training set separation (in descending order). | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | XY - training set, array[0..NPoints-1,0..NVars]. | |
| | | First NVars columns store values of independent | |
| | | variables, next column stores number of class (from 0 | |
| | | to NClasses-1) which dataset element belongs to. Fracti | |
| | | onal | |
| | | values are rounded to nearest integer. | |
| | | NPoints - training set size, NPoints>=0 | |
| | | NVars - number of independent variables, NVars>=1 | |
| | | NClasses - number of classes, NClasses>=2 | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info - return code: | |
| | | * -4, if internal EVD subroutine hasn't converged | |
| | | * -2, if there is a point with class number | |
| | | outside of [0..NClasses-1]. | |
| | | * -1, if incorrect parameters was passed (NPoints<0, | |
| | | NVars<1, NClasses<2) | |
| | | * 1, if task has been solved | |
| | | * 2, if there was a multicollinearity in training set, | |
| | | but task has been solved. | |
| | | W - basis, array[0..NVars-1,0..NVars-1] | |
| | | columns of matrix stores basis vectors, sorted by | |
| | | quality of training set separation (in descending order | |
| | | ) | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 31.05.2008 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void fisherldan(const real_2d_array &xy, const ae_int_t npoints, const ae_i | |
| | | nt_t nvars, const ae_int_t nclasses, ae_int_t &info, real_2d_array &w); | |
| | | | |
| | | /************************************************************************* | |
| This function serializes data structure to string. | | This function serializes data structure to string. | |
| | | | |
| Important properties of s_out: | | Important properties of s_out: | |
| * it contains alphanumeric characters, dots, underscores, minus signs | | * it contains alphanumeric characters, dots, underscores, minus signs | |
| * these symbols are grouped into words, which are separated by spaces | | * these symbols are grouped into words, which are separated by spaces | |
| and Windows-style (CR+LF) newlines | | and Windows-style (CR+LF) newlines | |
| * although serializer uses spaces and CR+LF as separators, you can | | * although serializer uses spaces and CR+LF as separators, you can | |
| replace any separator character by arbitrary combination of spaces, | | replace any separator character by arbitrary combination of spaces, | |
| tabs, Windows or Unix newlines. It allows flexible reformatting of | | tabs, Windows or Unix newlines. It allows flexible reformatting of | |
| the string in case you want to include it into text or XML file. | | the string in case you want to include it into text or XML file. | |
| | | | |
| skipping to change at line 1332 | | skipping to change at line 1435 | |
|---|
| /************************************************************************* | | /************************************************************************* | |
| Returns information about initialized network: number of inputs, outputs, | | Returns information about initialized network: number of inputs, outputs, | |
| weights. | | weights. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 04.11.2007 by Bochkanov Sergey | | Copyright 04.11.2007 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void mlpproperties(const multilayerperceptron &network, ae_int_t &nin, ae_i
nt_t &nout, ae_int_t &wcount); | | void mlpproperties(const multilayerperceptron &network, ae_int_t &nin, ae_i
nt_t &nout, ae_int_t &wcount); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | Returns number of inputs. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 19.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | ae_int_t mlpgetinputscount(const multilayerperceptron &network); | |
| | | | |
| | | /************************************************************************* | |
| | | Returns number of outputs. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 19.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | ae_int_t mlpgetoutputscount(const multilayerperceptron &network); | |
| | | | |
| | | /************************************************************************* | |
| | | Returns number of weights. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 19.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | ae_int_t mlpgetweightscount(const multilayerperceptron &network); | |
| | | | |
| | | /************************************************************************* | |
| Tells whether network is SOFTMAX-normalized (i.e. classifier) or not. | | Tells whether network is SOFTMAX-normalized (i.e. classifier) or not. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 04.11.2007 by Bochkanov Sergey | | Copyright 04.11.2007 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| bool mlpissoftmax(const multilayerperceptron &network); | | bool mlpissoftmax(const multilayerperceptron &network); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This function returns total number of layers (including input, hidden and | | This function returns total number of layers (including input, hidden and | |
| output layers). | | output layers). | |
| | | | |
| skipping to change at line 2700 | | skipping to change at line 2827 | |
|---|
| *************************************************************************/ | | *************************************************************************/ | |
| void mlptrainlbfgs(const multilayerperceptron &network, const real_2d_array
&xy, const ae_int_t npoints, const double decay, const ae_int_t restarts,
const double wstep, const ae_int_t maxits, ae_int_t &info, mlpreport &rep); | | void mlptrainlbfgs(const multilayerperceptron &network, const real_2d_array
&xy, const ae_int_t npoints, const double decay, const ae_int_t restarts,
const double wstep, const ae_int_t maxits, ae_int_t &info, mlpreport &rep); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| Neural network training using early stopping (base algorithm - L-BFGS with | | Neural network training using early stopping (base algorithm - L-BFGS with | |
| regularization). | | regularization). | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| Network - neural network with initialized geometry | | Network - neural network with initialized geometry | |
| TrnXY - training set | | TrnXY - training set | |
|
| TrnSize - training set size | | TrnSize - training set size, TrnSize>0 | |
| ValXY - validation set | | ValXY - validation set | |
|
| ValSize - validation set size | | ValSize - validation set size, ValSize>0 | |
| Decay - weight decay constant, >=0.001 | | Decay - weight decay constant, >=0.001 | |
| Decay term 'Decay*||Weights||^2' is added to error | | Decay term 'Decay*||Weights||^2' is added to error | |
| function. | | function. | |
| If you don't know what Decay to choose, use 0.001. | | If you don't know what Decay to choose, use 0.001. | |
|
| Restarts - number of restarts from random position, >0. | | Restarts - number of restarts, either: | |
| If you don't know what Restarts to choose, use 2. | | * strictly positive number - algorithm make specified | |
| | | number of restarts from random position. | |
| | | * -1, in which case algorithm makes exactly one run | |
| | | from the initial state of the network (no randomizati | |
| | | on). | |
| | | If you don't know what Restarts to choose, choose one | |
| | | one the following: | |
| | | * -1 (deterministic start) | |
| | | * +1 (one random restart) | |
| | | * +5 (moderate amount of random restarts) | |
| | | | |
| OUTPUT PARAMETERS: | | OUTPUT PARAMETERS: | |
| Network - trained neural network. | | Network - trained neural network. | |
| Info - return code: | | Info - return code: | |
| * -2, if there is a point with class number | | * -2, if there is a point with class number | |
| outside of [0..NOut-1]. | | outside of [0..NOut-1]. | |
| * -1, if wrong parameters specified | | * -1, if wrong parameters specified | |
| (NPoints<0, Restarts<1, ...). | | (NPoints<0, Restarts<1, ...). | |
| * 2, task has been solved, stopping criterion met - | | * 2, task has been solved, stopping criterion met - | |
| sufficiently small step size. Not expected (we | | sufficiently small step size. Not expected (we | |
| | | | |
| skipping to change at line 2798 | | skipping to change at line 2933 | |
|---|
| Info - return code, same as in MLPTrainLBFGS | | Info - return code, same as in MLPTrainLBFGS | |
| Rep - report, same as in MLPTrainLM/MLPTrainLBFGS | | Rep - report, same as in MLPTrainLM/MLPTrainLBFGS | |
| CVRep - generalization error estimates | | CVRep - generalization error estimates | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 09.12.2007 by Bochkanov Sergey | | Copyright 09.12.2007 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void mlpkfoldcvlm(const multilayerperceptron &network, const real_2d_array
&xy, const ae_int_t npoints, const double decay, const ae_int_t restarts, c
onst ae_int_t foldscount, ae_int_t &info, mlpreport &rep, mlpcvreport &cvre
p); | | void mlpkfoldcvlm(const multilayerperceptron &network, const real_2d_array
&xy, const ae_int_t npoints, const double decay, const ae_int_t restarts, c
onst ae_int_t foldscount, ae_int_t &info, mlpreport &rep, mlpcvreport &cvre
p); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | This function serializes data structure to string. | |
| | | | |
| | | Important properties of s_out: | |
| | | * it contains alphanumeric characters, dots, underscores, minus signs | |
| | | * these symbols are grouped into words, which are separated by spaces | |
| | | and Windows-style (CR+LF) newlines | |
| | | * although serializer uses spaces and CR+LF as separators, you can | |
| | | replace any separator character by arbitrary combination of spaces, | |
| | | tabs, Windows or Unix newlines. It allows flexible reformatting of | |
| | | the string in case you want to include it into text or XML file. | |
| | | But you should not insert separators into the middle of the "words" | |
| | | nor you should change case of letters. | |
| | | * s_out can be freely moved between 32-bit and 64-bit systems, little | |
| | | and big endian machines, and so on. You can serialize structure on | |
| | | 32-bit machine and unserialize it on 64-bit one (or vice versa), or | |
| | | serialize it on SPARC and unserialize on x86. You can also | |
| | | serialize it in C++ version of ALGLIB and unserialize in C# one, | |
| | | and vice versa. | |
| | | *************************************************************************/ | |
| | | void mlpeserialize(mlpensemble &obj, std::string &s_out); | |
| | | | |
| | | /************************************************************************* | |
| | | This function unserializes data structure from string. | |
| | | *************************************************************************/ | |
| | | void mlpeunserialize(std::string &s_in, mlpensemble &obj); | |
| | | | |
| | | /************************************************************************* | |
| Like MLPCreate0, but for ensembles. | | Like MLPCreate0, but for ensembles. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 18.02.2009 by Bochkanov Sergey | | Copyright 18.02.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void mlpecreate0(const ae_int_t nin, const ae_int_t nout, const ae_int_t en
semblesize, mlpensemble &ensemble); | | void mlpecreate0(const ae_int_t nin, const ae_int_t nout, const ae_int_t en
semblesize, mlpensemble &ensemble); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| Like MLPCreate1, but for ensembles. | | Like MLPCreate1, but for ensembles. | |
| | | | |
| | | | |
| skipping to change at line 3334 | | skipping to change at line 3496 | |
|---|
| ae_state *_state); | | ae_state *_state); | |
| ae_bool _decisionforest_init(decisionforest* p, ae_state *_state, ae_bool m
ake_automatic); | | ae_bool _decisionforest_init(decisionforest* p, ae_state *_state, ae_bool m
ake_automatic); | |
| ae_bool _decisionforest_init_copy(decisionforest* dst, decisionforest* src,
ae_state *_state, ae_bool make_automatic); | | ae_bool _decisionforest_init_copy(decisionforest* dst, decisionforest* src,
ae_state *_state, ae_bool make_automatic); | |
| void _decisionforest_clear(decisionforest* p); | | void _decisionforest_clear(decisionforest* p); | |
| ae_bool _dfreport_init(dfreport* p, ae_state *_state, ae_bool make_automati
c); | | ae_bool _dfreport_init(dfreport* p, ae_state *_state, ae_bool make_automati
c); | |
| ae_bool _dfreport_init_copy(dfreport* dst, dfreport* src, ae_state *_state,
ae_bool make_automatic); | | ae_bool _dfreport_init_copy(dfreport* dst, dfreport* src, ae_state *_state,
ae_bool make_automatic); | |
| void _dfreport_clear(dfreport* p); | | void _dfreport_clear(dfreport* p); | |
| ae_bool _dfinternalbuffers_init(dfinternalbuffers* p, ae_state *_state, ae_
bool make_automatic); | | ae_bool _dfinternalbuffers_init(dfinternalbuffers* p, ae_state *_state, ae_
bool make_automatic); | |
| ae_bool _dfinternalbuffers_init_copy(dfinternalbuffers* dst, dfinternalbuff
ers* src, ae_state *_state, ae_bool make_automatic); | | ae_bool _dfinternalbuffers_init_copy(dfinternalbuffers* dst, dfinternalbuff
ers* src, ae_state *_state, ae_bool make_automatic); | |
| void _dfinternalbuffers_clear(dfinternalbuffers* p); | | void _dfinternalbuffers_clear(dfinternalbuffers* p); | |
|
| void kmeansgenerate(/* Real */ ae_matrix* xy, | | | |
| ae_int_t npoints, | | | |
| ae_int_t nvars, | | | |
| ae_int_t k, | | | |
| ae_int_t restarts, | | | |
| ae_int_t* info, | | | |
| /* Real */ ae_matrix* c, | | | |
| /* Integer */ ae_vector* xyc, | | | |
| ae_state *_state); | | | |
| void fisherlda(/* Real */ ae_matrix* xy, | | | |
| ae_int_t npoints, | | | |
| ae_int_t nvars, | | | |
| ae_int_t nclasses, | | | |
| ae_int_t* info, | | | |
| /* Real */ ae_vector* w, | | | |
| ae_state *_state); | | | |
| void fisherldan(/* Real */ ae_matrix* xy, | | | |
| ae_int_t npoints, | | | |
| ae_int_t nvars, | | | |
| ae_int_t nclasses, | | | |
| ae_int_t* info, | | | |
| /* Real */ ae_matrix* w, | | | |
| ae_state *_state); | | | |
| void lrbuild(/* Real */ ae_matrix* xy, | | void lrbuild(/* Real */ ae_matrix* xy, | |
| ae_int_t npoints, | | ae_int_t npoints, | |
| ae_int_t nvars, | | ae_int_t nvars, | |
| ae_int_t* info, | | ae_int_t* info, | |
| linearmodel* lm, | | linearmodel* lm, | |
| lrreport* ar, | | lrreport* ar, | |
| ae_state *_state); | | ae_state *_state); | |
| void lrbuilds(/* Real */ ae_matrix* xy, | | void lrbuilds(/* Real */ ae_matrix* xy, | |
| /* Real */ ae_vector* s, | | /* Real */ ae_vector* s, | |
| ae_int_t npoints, | | ae_int_t npoints, | |
| | | | |
| skipping to change at line 3435 | | skipping to change at line 3574 | |
|---|
| ae_int_t* info, | | ae_int_t* info, | |
| double* a, | | double* a, | |
| double* b, | | double* b, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_bool _linearmodel_init(linearmodel* p, ae_state *_state, ae_bool make_au
tomatic); | | ae_bool _linearmodel_init(linearmodel* p, ae_state *_state, ae_bool make_au
tomatic); | |
| ae_bool _linearmodel_init_copy(linearmodel* dst, linearmodel* src, ae_state
*_state, ae_bool make_automatic); | | ae_bool _linearmodel_init_copy(linearmodel* dst, linearmodel* src, ae_state
*_state, ae_bool make_automatic); | |
| void _linearmodel_clear(linearmodel* p); | | void _linearmodel_clear(linearmodel* p); | |
| ae_bool _lrreport_init(lrreport* p, ae_state *_state, ae_bool make_automati
c); | | ae_bool _lrreport_init(lrreport* p, ae_state *_state, ae_bool make_automati
c); | |
| ae_bool _lrreport_init_copy(lrreport* dst, lrreport* src, ae_state *_state,
ae_bool make_automatic); | | ae_bool _lrreport_init_copy(lrreport* dst, lrreport* src, ae_state *_state,
ae_bool make_automatic); | |
| void _lrreport_clear(lrreport* p); | | void _lrreport_clear(lrreport* p); | |
|
| | | void filtersma(/* Real */ ae_vector* x, | |
| | | ae_int_t n, | |
| | | ae_int_t k, | |
| | | ae_state *_state); | |
| | | void filterema(/* Real */ ae_vector* x, | |
| | | ae_int_t n, | |
| | | double alpha, | |
| | | ae_state *_state); | |
| | | void filterlrma(/* Real */ ae_vector* x, | |
| | | ae_int_t n, | |
| | | ae_int_t k, | |
| | | ae_state *_state); | |
| | | void kmeansgenerate(/* Real */ ae_matrix* xy, | |
| | | ae_int_t npoints, | |
| | | ae_int_t nvars, | |
| | | ae_int_t k, | |
| | | ae_int_t restarts, | |
| | | ae_int_t* info, | |
| | | /* Real */ ae_matrix* c, | |
| | | /* Integer */ ae_vector* xyc, | |
| | | ae_state *_state); | |
| | | void fisherlda(/* Real */ ae_matrix* xy, | |
| | | ae_int_t npoints, | |
| | | ae_int_t nvars, | |
| | | ae_int_t nclasses, | |
| | | ae_int_t* info, | |
| | | /* Real */ ae_vector* w, | |
| | | ae_state *_state); | |
| | | void fisherldan(/* Real */ ae_matrix* xy, | |
| | | ae_int_t npoints, | |
| | | ae_int_t nvars, | |
| | | ae_int_t nclasses, | |
| | | ae_int_t* info, | |
| | | /* Real */ ae_matrix* w, | |
| | | ae_state *_state); | |
| void mlpcreate0(ae_int_t nin, | | void mlpcreate0(ae_int_t nin, | |
| ae_int_t nout, | | ae_int_t nout, | |
| multilayerperceptron* network, | | multilayerperceptron* network, | |
| ae_state *_state); | | ae_state *_state); | |
| void mlpcreate1(ae_int_t nin, | | void mlpcreate1(ae_int_t nin, | |
| ae_int_t nhid, | | ae_int_t nhid, | |
| ae_int_t nout, | | ae_int_t nout, | |
| multilayerperceptron* network, | | multilayerperceptron* network, | |
| ae_state *_state); | | ae_state *_state); | |
| void mlpcreate2(ae_int_t nin, | | void mlpcreate2(ae_int_t nin, | |
| | | | |
| skipping to change at line 3528 | | skipping to change at line 3702 | |
|---|
| void mlprandomizefull(multilayerperceptron* network, ae_state *_state); | | void mlprandomizefull(multilayerperceptron* network, ae_state *_state); | |
| void mlpinitpreprocessor(multilayerperceptron* network, | | void mlpinitpreprocessor(multilayerperceptron* network, | |
| /* Real */ ae_matrix* xy, | | /* Real */ ae_matrix* xy, | |
| ae_int_t ssize, | | ae_int_t ssize, | |
| ae_state *_state); | | ae_state *_state); | |
| void mlpproperties(multilayerperceptron* network, | | void mlpproperties(multilayerperceptron* network, | |
| ae_int_t* nin, | | ae_int_t* nin, | |
| ae_int_t* nout, | | ae_int_t* nout, | |
| ae_int_t* wcount, | | ae_int_t* wcount, | |
| ae_state *_state); | | ae_state *_state); | |
|
| | | ae_int_t mlpgetinputscount(multilayerperceptron* network, | |
| | | ae_state *_state); | |
| | | ae_int_t mlpgetoutputscount(multilayerperceptron* network, | |
| | | ae_state *_state); | |
| | | ae_int_t mlpgetweightscount(multilayerperceptron* network, | |
| | | ae_state *_state); | |
| ae_bool mlpissoftmax(multilayerperceptron* network, ae_state *_state); | | ae_bool mlpissoftmax(multilayerperceptron* network, ae_state *_state); | |
| ae_int_t mlpgetlayerscount(multilayerperceptron* network, | | ae_int_t mlpgetlayerscount(multilayerperceptron* network, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_int_t mlpgetlayersize(multilayerperceptron* network, | | ae_int_t mlpgetlayersize(multilayerperceptron* network, | |
| ae_int_t k, | | ae_int_t k, | |
| ae_state *_state); | | ae_state *_state); | |
| void mlpgetinputscaling(multilayerperceptron* network, | | void mlpgetinputscaling(multilayerperceptron* network, | |
| ae_int_t i, | | ae_int_t i, | |
| double* mean, | | double* mean, | |
| double* sigma, | | double* sigma, | |
| | | | |
| skipping to change at line 3950 | | skipping to change at line 4130 | |
|---|
| ae_int_t ensemblesize, | | ae_int_t ensemblesize, | |
| mlpensemble* ensemble, | | mlpensemble* ensemble, | |
| ae_state *_state); | | ae_state *_state); | |
| void mlpecreatefromnetwork(multilayerperceptron* network, | | void mlpecreatefromnetwork(multilayerperceptron* network, | |
| ae_int_t ensemblesize, | | ae_int_t ensemblesize, | |
| mlpensemble* ensemble, | | mlpensemble* ensemble, | |
| ae_state *_state); | | ae_state *_state); | |
| void mlpecopy(mlpensemble* ensemble1, | | void mlpecopy(mlpensemble* ensemble1, | |
| mlpensemble* ensemble2, | | mlpensemble* ensemble2, | |
| ae_state *_state); | | ae_state *_state); | |
|
| void mlpeserialize(mlpensemble* ensemble, | | | |
| /* Real */ ae_vector* ra, | | | |
| ae_int_t* rlen, | | | |
| ae_state *_state); | | | |
| void mlpeunserialize(/* Real */ ae_vector* ra, | | | |
| mlpensemble* ensemble, | | | |
| ae_state *_state); | | | |
| void mlperandomize(mlpensemble* ensemble, ae_state *_state); | | void mlperandomize(mlpensemble* ensemble, ae_state *_state); | |
| void mlpeproperties(mlpensemble* ensemble, | | void mlpeproperties(mlpensemble* ensemble, | |
| ae_int_t* nin, | | ae_int_t* nin, | |
| ae_int_t* nout, | | ae_int_t* nout, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_bool mlpeissoftmax(mlpensemble* ensemble, ae_state *_state); | | ae_bool mlpeissoftmax(mlpensemble* ensemble, ae_state *_state); | |
| void mlpeprocess(mlpensemble* ensemble, | | void mlpeprocess(mlpensemble* ensemble, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| ae_state *_state); | | ae_state *_state); | |
| | | | |
| skipping to change at line 4019 | | skipping to change at line 4192 | |
|---|
| mlpcvreport* ooberrors, | | mlpcvreport* ooberrors, | |
| ae_state *_state); | | ae_state *_state); | |
| void mlpetraines(mlpensemble* ensemble, | | void mlpetraines(mlpensemble* ensemble, | |
| /* Real */ ae_matrix* xy, | | /* Real */ ae_matrix* xy, | |
| ae_int_t npoints, | | ae_int_t npoints, | |
| double decay, | | double decay, | |
| ae_int_t restarts, | | ae_int_t restarts, | |
| ae_int_t* info, | | ae_int_t* info, | |
| mlpreport* rep, | | mlpreport* rep, | |
| ae_state *_state); | | ae_state *_state); | |
|
| | | void mlpealloc(ae_serializer* s, mlpensemble* ensemble, ae_state *_state); | |
| | | void mlpeserialize(ae_serializer* s, | |
| | | mlpensemble* ensemble, | |
| | | ae_state *_state); | |
| | | void mlpeunserialize(ae_serializer* s, | |
| | | mlpensemble* ensemble, | |
| | | ae_state *_state); | |
| ae_bool _mlpensemble_init(mlpensemble* p, ae_state *_state, ae_bool make_au
tomatic); | | ae_bool _mlpensemble_init(mlpensemble* p, ae_state *_state, ae_bool make_au
tomatic); | |
| ae_bool _mlpensemble_init_copy(mlpensemble* dst, mlpensemble* src, ae_state
*_state, ae_bool make_automatic); | | ae_bool _mlpensemble_init_copy(mlpensemble* dst, mlpensemble* src, ae_state
*_state, ae_bool make_automatic); | |
| void _mlpensemble_clear(mlpensemble* p); | | void _mlpensemble_clear(mlpensemble* p); | |
| void pcabuildbasis(/* Real */ ae_matrix* x, | | void pcabuildbasis(/* Real */ ae_matrix* x, | |
| ae_int_t npoints, | | ae_int_t npoints, | |
| ae_int_t nvars, | | ae_int_t nvars, | |
| ae_int_t* info, | | ae_int_t* info, | |
| /* Real */ ae_vector* s2, | | /* Real */ ae_vector* s2, | |
| /* Real */ ae_matrix* v, | | /* Real */ ae_matrix* v, | |
| ae_state *_state); | | ae_state *_state); | |
| | | | |
| End of changes. 17 change blocks. |
|---|
| 154 lines changed or deleted | | 337 lines changed or added | |
|---|
|
| interpolation.h | | interpolation.h | |
|---|
| | | | |
| skipping to change at line 68 | | skipping to change at line 68 | |
|---|
| double sy; | | double sy; | |
| ae_vector x; | | ae_vector x; | |
| ae_vector y; | | ae_vector y; | |
| ae_vector w; | | ae_vector w; | |
| } barycentricinterpolant; | | } barycentricinterpolant; | |
| typedef struct | | typedef struct | |
| { | | { | |
| ae_bool periodic; | | ae_bool periodic; | |
| ae_int_t n; | | ae_int_t n; | |
| ae_int_t k; | | ae_int_t k; | |
|
| | | ae_int_t continuity; | |
| ae_vector x; | | ae_vector x; | |
| ae_vector c; | | ae_vector c; | |
| } spline1dinterpolant; | | } spline1dinterpolant; | |
| typedef struct | | typedef struct | |
| { | | { | |
| double taskrcond; | | double taskrcond; | |
| double rmserror; | | double rmserror; | |
| double avgerror; | | double avgerror; | |
| double avgrelerror; | | double avgrelerror; | |
| double maxerror; | | double maxerror; | |
| | | | |
| skipping to change at line 168 | | skipping to change at line 169 | |
|---|
| { | | { | |
| ae_int_t n; | | ae_int_t n; | |
| ae_bool periodic; | | ae_bool periodic; | |
| ae_vector p; | | ae_vector p; | |
| spline1dinterpolant x; | | spline1dinterpolant x; | |
| spline1dinterpolant y; | | spline1dinterpolant y; | |
| spline1dinterpolant z; | | spline1dinterpolant z; | |
| } pspline3interpolant; | | } pspline3interpolant; | |
| typedef struct | | typedef struct | |
| { | | { | |
|
| | | ae_int_t ny; | |
| | | ae_int_t nx; | |
| | | ae_int_t nc; | |
| | | ae_int_t nl; | |
| | | kdtree tree; | |
| | | ae_matrix xc; | |
| | | ae_matrix wr; | |
| | | double rmax; | |
| | | ae_matrix v; | |
| | | ae_int_t gridtype; | |
| | | ae_bool fixrad; | |
| | | double lambdav; | |
| | | double radvalue; | |
| | | double radzvalue; | |
| | | ae_int_t nlayers; | |
| | | ae_int_t aterm; | |
| | | ae_int_t algorithmtype; | |
| | | double epsort; | |
| | | double epserr; | |
| | | ae_int_t maxits; | |
| | | double h; | |
| | | ae_int_t n; | |
| | | ae_matrix x; | |
| | | ae_matrix y; | |
| | | ae_vector calcbufxcx; | |
| | | ae_matrix calcbufx; | |
| | | ae_vector calcbuftags; | |
| | | } rbfmodel; | |
| | | typedef struct | |
| | | { | |
| | | ae_int_t arows; | |
| | | ae_int_t acols; | |
| | | ae_int_t annz; | |
| | | ae_int_t iterationscount; | |
| | | ae_int_t nmv; | |
| | | ae_int_t terminationtype; | |
| | | } rbfreport; | |
| | | typedef struct | |
| | | { | |
| ae_int_t k; | | ae_int_t k; | |
| ae_vector c; | | ae_vector c; | |
| } spline2dinterpolant; | | } spline2dinterpolant; | |
| | | | |
| } | | } | |
| | | | |
| ///////////////////////////////////////////////////////////////////////// | | ///////////////////////////////////////////////////////////////////////// | |
| // | | // | |
| // THIS SECTION CONTAINS C++ INTERFACE | | // THIS SECTION CONTAINS C++ INTERFACE | |
| // | | // | |
| | | | |
| skipping to change at line 233 | | skipping to change at line 273 | |
|---|
| { | | { | |
| public: | | public: | |
| barycentricinterpolant(); | | barycentricinterpolant(); | |
| barycentricinterpolant(const barycentricinterpolant &rhs); | | barycentricinterpolant(const barycentricinterpolant &rhs); | |
| barycentricinterpolant& operator=(const barycentricinterpolant &rhs); | | barycentricinterpolant& operator=(const barycentricinterpolant &rhs); | |
| virtual ~barycentricinterpolant(); | | virtual ~barycentricinterpolant(); | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| 1-dimensional spline inteprolant | | 1-dimensional spline interpolant | |
| *************************************************************************/ | | *************************************************************************/ | |
| class _spline1dinterpolant_owner | | class _spline1dinterpolant_owner | |
| { | | { | |
| public: | | public: | |
| _spline1dinterpolant_owner(); | | _spline1dinterpolant_owner(); | |
| _spline1dinterpolant_owner(const _spline1dinterpolant_owner &rhs); | | _spline1dinterpolant_owner(const _spline1dinterpolant_owner &rhs); | |
| _spline1dinterpolant_owner& operator=(const _spline1dinterpolant_owner
&rhs); | | _spline1dinterpolant_owner& operator=(const _spline1dinterpolant_owner
&rhs); | |
| virtual ~_spline1dinterpolant_owner(); | | virtual ~_spline1dinterpolant_owner(); | |
| alglib_impl::spline1dinterpolant* c_ptr(); | | alglib_impl::spline1dinterpolant* c_ptr(); | |
| alglib_impl::spline1dinterpolant* c_ptr() const; | | alglib_impl::spline1dinterpolant* c_ptr() const; | |
| | | | |
| skipping to change at line 501 | | skipping to change at line 541 | |
|---|
| { | | { | |
| public: | | public: | |
| pspline3interpolant(); | | pspline3interpolant(); | |
| pspline3interpolant(const pspline3interpolant &rhs); | | pspline3interpolant(const pspline3interpolant &rhs); | |
| pspline3interpolant& operator=(const pspline3interpolant &rhs); | | pspline3interpolant& operator=(const pspline3interpolant &rhs); | |
| virtual ~pspline3interpolant(); | | virtual ~pspline3interpolant(); | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | RBF model. | |
| | | | |
| | | Never try to directly work with fields of this object - always use ALGLIB | |
| | | functions to use this object. | |
| | | *************************************************************************/ | |
| | | class _rbfmodel_owner | |
| | | { | |
| | | public: | |
| | | _rbfmodel_owner(); | |
| | | _rbfmodel_owner(const _rbfmodel_owner &rhs); | |
| | | _rbfmodel_owner& operator=(const _rbfmodel_owner &rhs); | |
| | | virtual ~_rbfmodel_owner(); | |
| | | alglib_impl::rbfmodel* c_ptr(); | |
| | | alglib_impl::rbfmodel* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::rbfmodel *p_struct; | |
| | | }; | |
| | | class rbfmodel : public _rbfmodel_owner | |
| | | { | |
| | | public: | |
| | | rbfmodel(); | |
| | | rbfmodel(const rbfmodel &rhs); | |
| | | rbfmodel& operator=(const rbfmodel &rhs); | |
| | | virtual ~rbfmodel(); | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| | | RBF solution report: | |
| | | * TerminationType - termination type, positive values - success, | |
| | | non-positive - failure. | |
| | | *************************************************************************/ | |
| | | class _rbfreport_owner | |
| | | { | |
| | | public: | |
| | | _rbfreport_owner(); | |
| | | _rbfreport_owner(const _rbfreport_owner &rhs); | |
| | | _rbfreport_owner& operator=(const _rbfreport_owner &rhs); | |
| | | virtual ~_rbfreport_owner(); | |
| | | alglib_impl::rbfreport* c_ptr(); | |
| | | alglib_impl::rbfreport* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::rbfreport *p_struct; | |
| | | }; | |
| | | class rbfreport : public _rbfreport_owner | |
| | | { | |
| | | public: | |
| | | rbfreport(); | |
| | | rbfreport(const rbfreport &rhs); | |
| | | rbfreport& operator=(const rbfreport &rhs); | |
| | | virtual ~rbfreport(); | |
| | | ae_int_t &arows; | |
| | | ae_int_t &acols; | |
| | | ae_int_t &annz; | |
| | | ae_int_t &iterationscount; | |
| | | ae_int_t &nmv; | |
| | | ae_int_t &terminationtype; | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| 2-dimensional spline inteprolant | | 2-dimensional spline inteprolant | |
| *************************************************************************/ | | *************************************************************************/ | |
| class _spline2dinterpolant_owner | | class _spline2dinterpolant_owner | |
| { | | { | |
| public: | | public: | |
| _spline2dinterpolant_owner(); | | _spline2dinterpolant_owner(); | |
| _spline2dinterpolant_owner(const _spline2dinterpolant_owner &rhs); | | _spline2dinterpolant_owner(const _spline2dinterpolant_owner &rhs); | |
| _spline2dinterpolant_owner& operator=(const _spline2dinterpolant_owner
&rhs); | | _spline2dinterpolant_owner& operator=(const _spline2dinterpolant_owner
&rhs); | |
| virtual ~_spline2dinterpolant_owner(); | | virtual ~_spline2dinterpolant_owner(); | |
| alglib_impl::spline2dinterpolant* c_ptr(); | | alglib_impl::spline2dinterpolant* c_ptr(); | |
| | | | |
| skipping to change at line 3358 | | skipping to change at line 3459 | |
|---|
| | | | |
| RESULT: | | RESULT: | |
| length of arc starting at T=A and ending at T=B. | | length of arc starting at T=A and ending at T=B. | |
| | | | |
| -- ALGLIB PROJECT -- | | -- ALGLIB PROJECT -- | |
| Copyright 30.05.2010 by Bochkanov Sergey | | Copyright 30.05.2010 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| double pspline3arclength(const pspline3interpolant &p, const double a, cons
t double b); | | double pspline3arclength(const pspline3interpolant &p, const double a, cons
t double b); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | This function serializes data structure to string. | |
| | | | |
| | | Important properties of s_out: | |
| | | * it contains alphanumeric characters, dots, underscores, minus signs | |
| | | * these symbols are grouped into words, which are separated by spaces | |
| | | and Windows-style (CR+LF) newlines | |
| | | * although serializer uses spaces and CR+LF as separators, you can | |
| | | replace any separator character by arbitrary combination of spaces, | |
| | | tabs, Windows or Unix newlines. It allows flexible reformatting of | |
| | | the string in case you want to include it into text or XML file. | |
| | | But you should not insert separators into the middle of the "words" | |
| | | nor you should change case of letters. | |
| | | * s_out can be freely moved between 32-bit and 64-bit systems, little | |
| | | and big endian machines, and so on. You can serialize structure on | |
| | | 32-bit machine and unserialize it on 64-bit one (or vice versa), or | |
| | | serialize it on SPARC and unserialize on x86. You can also | |
| | | serialize it in C++ version of ALGLIB and unserialize in C# one, | |
| | | and vice versa. | |
| | | *************************************************************************/ | |
| | | void rbfserialize(rbfmodel &obj, std::string &s_out); | |
| | | | |
| | | /************************************************************************* | |
| | | This function unserializes data structure from string. | |
| | | *************************************************************************/ | |
| | | void rbfunserialize(std::string &s_in, rbfmodel &obj); | |
| | | | |
| | | /************************************************************************* | |
| | | This function creates RBF model for a scalar (NY=1) or vector (NY>1) | |
| | | function in a NX-dimensional space (NX=2 or NX=3). | |
| | | | |
| | | Newly created model is empty. It can be used for interpolation right after | |
| | | creation, but it just returns zeros. You have to add points to the model, | |
| | | tune interpolation settings, and then call model construction function | |
| | | RBFBuildModel() which will update model according to your specification. | |
| | | | |
| | | USAGE: | |
| | | 1. User creates model with RBFCreate() | |
| | | 2. User adds dataset with RBFSetPoints() (points do NOT have to be on a | |
| | | regular grid) | |
| | | 3. (OPTIONAL) User chooses polynomial term by calling: | |
| | | * RBFLinTerm() to set linear term | |
| | | * RBFConstTerm() to set constant term | |
| | | * RBFZeroTerm() to set zero term | |
| | | By default, linear term is used. | |
| | | 4. User chooses specific RBF algorithm to use: either QNN (RBFSetAlgoQNN) | |
| | | or ML (RBFSetAlgoMultiLayer). | |
| | | 5. User calls RBFBuildModel() function which rebuilds model according to | |
| | | the specification | |
| | | 6. User may call RBFCalc() to calculate model value at the specified point, | |
| | | RBFGridCalc() to calculate model values at the points of the regular | |
| | | grid. User may extract model coefficients with RBFUnpack() call. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | NX - dimension of the space, NX=2 or NX=3 | |
| | | NY - function dimension, NY>=1 | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | S - RBF model (initially equals to zero) | |
| | | | |
| | | NOTE 1: memory requirements. RBF models require amount of memory which is | |
| | | proportional to the number of data points. Memory is allocated | |
| | | during model construction, but most of this memory is freed after | |
| | | model coefficients are calculated. | |
| | | | |
| | | Some approximate estimates for N centers with default settings are | |
| | | given below: | |
| | | * about 250*N*(sizeof(double)+2*sizeof(int)) bytes of memory is | |
| | | needed during model construction stage. | |
| | | * about 15*N*sizeof(double) bytes is needed after model is built. | |
| | | For example, for N=100000 we may need 0.6 GB of memory to build | |
| | | model, but just about 0.012 GB to store it. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfcreate(const ae_int_t nx, const ae_int_t ny, rbfmodel &s); | |
| | | | |
| | | /************************************************************************* | |
| | | This function adds dataset. | |
| | | | |
| | | This function overrides results of the previous calls, i.e. multiple calls | |
| | | of this function will result in only the last set being added. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model, initialized by RBFCreate() call. | |
| | | XY - points, array[N,NX+NY]. One row corresponds to one point | |
| | | in the dataset. First NX elements are coordinates, next | |
| | | NY elements are function values. Array may be larger than | |
| | | specific, in this case only leading [N,NX+NY] elements | |
| | | will be used. | |
| | | N - number of points in the dataset | |
| | | | |
| | | After you've added dataset and (optionally) tuned algorithm settings you | |
| | | should call RBFBuildModel() in order to build a model for you. | |
| | | | |
| | | NOTE: this function has some serialization-related subtleties. We | |
| | | recommend you to study serialization examples from ALGLIB Reference | |
| | | Manual if you want to perform serialization of your models. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfsetpoints(const rbfmodel &s, const real_2d_array &xy, const ae_int_ | |
| | | t n); | |
| | | void rbfsetpoints(const rbfmodel &s, const real_2d_array &xy); | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets RBF interpolation algorithm. ALGLIB supports several | |
| | | RBF algorithms with different properties. | |
| | | | |
| | | This algorithm is called RBF-QNN and it is good for point sets with | |
| | | following properties: | |
| | | a) all points are distinct | |
| | | b) all points are well separated. | |
| | | c) points distribution is approximately uniform. There is no "contour | |
| | | lines", clusters of points, or other small-scale structures. | |
| | | | |
| | | Algorithm description: | |
| | | 1) interpolation centers are allocated to data points | |
| | | 2) interpolation radii are calculated as distances to the nearest centers | |
| | | times Q coefficient (where Q is a value from [0.75,1.50]). | |
| | | 3) after performing (2) radii are transformed in order to avoid situation | |
| | | when single outlier has very large radius and influences many points | |
| | | across all dataset. Transformation has following form: | |
| | | new_r[i] = min(r[i],Z*median(r[])) | |
| | | where r[i] is I-th radius, median() is a median radius across entire | |
| | | dataset, Z is user-specified value which controls amount of deviation | |
| | | from median radius. | |
| | | | |
| | | When (a) is violated, we will be unable to build RBF model. When (b) or | |
| | | (c) are violated, model will be built, but interpolation quality will be | |
| | | low. See http://www.alglib.net/interpolation/ for more information on this | |
| | | subject. | |
| | | | |
| | | This algorithm is used by default. | |
| | | | |
| | | Additional Q parameter controls smoothness properties of the RBF basis: | |
| | | * Q<0.75 will give perfectly conditioned basis, but terrible smoothness | |
| | | properties (RBF interpolant will have sharp peaks around function values) | |
| | | * Q around 1.0 gives good balance between smoothness and condition number | |
| | | * Q>1.5 will lead to badly conditioned systems and slow convergence of the | |
| | | underlying linear solver (although smoothness will be very good) | |
| | | * Q>2.0 will effectively make optimizer useless because it won't converge | |
| | | within reasonable amount of iterations. It is possible to set such large | |
| | | Q, but it is advised not to do so. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model, initialized by RBFCreate() call | |
| | | Q - Q parameter, Q>0, recommended value - 1.0 | |
| | | Z - Z parameter, Z>0, recommended value - 5.0 | |
| | | | |
| | | NOTE: this function has some serialization-related subtleties. We | |
| | | recommend you to study serialization examples from ALGLIB Reference | |
| | | Manual if you want to perform serialization of your models. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfsetalgoqnn(const rbfmodel &s, const double q, const double z); | |
| | | void rbfsetalgoqnn(const rbfmodel &s); | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets RBF interpolation algorithm. ALGLIB supports several | |
| | | RBF algorithms with different properties. | |
| | | | |
| | | This algorithm is called RBF-ML. It builds multilayer RBF model, i.e. | |
| | | model with subsequently decreasing radii, which allows us to combine | |
| | | smoothness (due to large radii of the first layers) with exactness (due | |
| | | to small radii of the last layers) and fast convergence. | |
| | | | |
| | | Internally RBF-ML uses many different means of acceleration, from sparse | |
| | | matrices to KD-trees, which results in algorithm whose working time is | |
| | | roughly proportional to N*log(N)*Density*RBase^2*NLayers, where N is a | |
| | | number of points, Density is an average density if points per unit of the | |
| | | interpolation space, RBase is an initial radius, NLayers is a number of | |
| | | layers. | |
| | | | |
| | | RBF-ML is good for following kinds of interpolation problems: | |
| | | 1. "exact" problems (perfect fit) with well separated points | |
| | | 2. least squares problems with arbitrary distribution of points (algorithm | |
| | | gives perfect fit where it is possible, and resorts to least squares | |
| | | fit in the hard areas). | |
| | | 3. noisy problems where we want to apply some controlled amount of | |
| | | smoothing. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model, initialized by RBFCreate() call | |
| | | RBase - RBase parameter, RBase>0 | |
| | | NLayers - NLayers parameter, NLayers>0, recommended value to start | |
| | | with - about 5. | |
| | | LambdaV - regularization value, can be useful when solving problem | |
| | | in the least squares sense. Optimal lambda is problem- | |
| | | dependent and require trial and error. In our experience, | |
| | | good lambda can be as large as 0.1, and you can use 0.001 | |
| | | as initial guess. | |
| | | Default value - 0.01, which is used when LambdaV is not | |
| | | given. You can specify zero value, but it is not | |
| | | recommended to do so. | |
| | | | |
| | | TUNING ALGORITHM | |
| | | | |
| | | In order to use this algorithm you have to choose three parameters: | |
| | | * initial radius RBase | |
| | | * number of layers in the model NLayers | |
| | | * regularization coefficient LambdaV | |
| | | | |
| | | Initial radius is easy to choose - you can pick any number several times | |
| | | larger than the average distance between points. Algorithm won't break | |
| | | down if you choose radius which is too large (model construction time will | |
| | | increase, but model will be built correctly). | |
| | | | |
| | | Choose such number of layers that RLast=RBase/2^(NLayers-1) (radius used | |
| | | by the last layer) will be smaller than the typical distance between | |
| | | points. In case model error is too large, you can increase number of | |
| | | layers. Having more layers will make model construction and evaluation | |
| | | proportionally slower, but it will allow you to have model which precisely | |
| | | fits your data. From the other side, if you want to suppress noise, you | |
| | | can DECREASE number of layers to make your model less flexible. | |
| | | | |
| | | Regularization coefficient LambdaV controls smoothness of the individual | |
| | | models built for each layer. We recommend you to use default value in case | |
| | | you don't want to tune this parameter, because having non-zero LambdaV | |
| | | accelerates and stabilizes internal iterative algorithm. In case you want | |
| | | to suppress noise you can use LambdaV as additional parameter (larger | |
| | | value = more smoothness) to tune. | |
| | | | |
| | | TYPICAL ERRORS | |
| | | | |
| | | 1. Using initial radius which is too large. Memory requirements of the | |
| | | RBF-ML are roughly proportional to N*Density*RBase^2 (where Density is | |
| | | an average density of points per unit of the interpolation space). In | |
| | | the extreme case of the very large RBase we will need O(N^2) units of | |
| | | memory - and many layers in order to decrease radius to some reasonably | |
| | | small value. | |
| | | | |
| | | 2. Using too small number of layers - RBF models with large radius are not | |
| | | flexible enough to reproduce small variations in the target function. | |
| | | You need many layers with different radii, from large to small, in | |
| | | order to have good model. | |
| | | | |
| | | 3. Using initial radius which is too small. You will get model with | |
| | | "holes" in the areas which are too far away from interpolation centers. | |
| | | However, algorithm will work correctly (and quickly) in this case. | |
| | | | |
| | | 4. Using too many layers - you will get too large and too slow model. This | |
| | | model will perfectly reproduce your function, but maybe you will be | |
| | | able to achieve similar results with less layers (and less memory). | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 02.03.2012 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfsetalgomultilayer(const rbfmodel &s, const double rbase, const ae_i | |
| | | nt_t nlayers, const double lambdav); | |
| | | void rbfsetalgomultilayer(const rbfmodel &s, const double rbase, const ae_i | |
| | | nt_t nlayers); | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets linear term (model is a sum of radial basis functions | |
| | | plus linear polynomial). This function won't have effect until next call | |
| | | to RBFBuildModel(). | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model, initialized by RBFCreate() call | |
| | | | |
| | | NOTE: this function has some serialization-related subtleties. We | |
| | | recommend you to study serialization examples from ALGLIB Reference | |
| | | Manual if you want to perform serialization of your models. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfsetlinterm(const rbfmodel &s); | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets constant term (model is a sum of radial basis functions | |
| | | plus constant). This function won't have effect until next call to | |
| | | RBFBuildModel(). | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model, initialized by RBFCreate() call | |
| | | | |
| | | NOTE: this function has some serialization-related subtleties. We | |
| | | recommend you to study serialization examples from ALGLIB Reference | |
| | | Manual if you want to perform serialization of your models. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfsetconstterm(const rbfmodel &s); | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets zero term (model is a sum of radial basis functions | |
| | | without polynomial term). This function won't have effect until next call | |
| | | to RBFBuildModel(). | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model, initialized by RBFCreate() call | |
| | | | |
| | | NOTE: this function has some serialization-related subtleties. We | |
| | | recommend you to study serialization examples from ALGLIB Reference | |
| | | Manual if you want to perform serialization of your models. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfsetzeroterm(const rbfmodel &s); | |
| | | | |
| | | /************************************************************************* | |
| | | This function builds RBF model and returns report (contains some | |
| | | information which can be used for evaluation of the algorithm properties). | |
| | | | |
| | | Call to this function modifies RBF model by calculating its centers/radii/ | |
| | | weights and saving them into RBFModel structure. Initially RBFModel | |
| | | contain zero coefficients, but after call to this function we will have | |
| | | coefficients which were calculated in order to fit our dataset. | |
| | | | |
| | | After you called this function you can call RBFCalc(), RBFGridCalc() and | |
| | | other model calculation functions. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model, initialized by RBFCreate() call | |
| | | Rep - report: | |
| | | * Rep.TerminationType: | |
| | | * -5 - non-distinct basis function centers were detected, | |
| | | interpolation aborted | |
| | | * -4 - nonconvergence of the internal SVD solver | |
| | | * 1 - successful termination | |
| | | Fields are used for debugging purposes: | |
| | | * Rep.IterationsCount - iterations count of the LSQR solver | |
| | | * Rep.NMV - number of matrix-vector products | |
| | | * Rep.ARows - rows count for the system matrix | |
| | | * Rep.ACols - columns count for the system matrix | |
| | | * Rep.ANNZ - number of significantly non-zero elements | |
| | | (elements above some algorithm-determined threshold) | |
| | | | |
| | | NOTE: failure to build model will leave current state of the structure | |
| | | unchanged. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfbuildmodel(const rbfmodel &s, rbfreport &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates values of the RBF model in the given point. | |
| | | | |
| | | This function should be used when we have NY=1 (scalar function) and NX=2 | |
| | | (2-dimensional space). If you have 3-dimensional space, use RBFCalc3(). If | |
| | | you have general situation (NX-dimensional space, NY-dimensional function) | |
| | | you should use general, less efficient implementation RBFCalc(). | |
| | | | |
| | | If you want to calculate function values many times, consider using | |
| | | RBFGridCalc2(), which is far more efficient than many subsequent calls to | |
| | | RBFCalc2(). | |
| | | | |
| | | This function returns 0.0 when: | |
| | | * model is not initialized | |
| | | * NX<>2 | |
| | | *NY<>1 | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model | |
| | | X0 - first coordinate, finite number | |
| | | X1 - second coordinate, finite number | |
| | | | |
| | | RESULT: | |
| | | value of the model or 0.0 (as defined above) | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | double rbfcalc2(const rbfmodel &s, const double x0, const double x1); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates values of the RBF model in the given point. | |
| | | | |
| | | This function should be used when we have NY=1 (scalar function) and NX=3 | |
| | | (3-dimensional space). If you have 2-dimensional space, use RBFCalc2(). If | |
| | | you have general situation (NX-dimensional space, NY-dimensional function) | |
| | | you should use general, less efficient implementation RBFCalc(). | |
| | | | |
| | | This function returns 0.0 when: | |
| | | * model is not initialized | |
| | | * NX<>3 | |
| | | *NY<>1 | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model | |
| | | X0 - first coordinate, finite number | |
| | | X1 - second coordinate, finite number | |
| | | X2 - third coordinate, finite number | |
| | | | |
| | | RESULT: | |
| | | value of the model or 0.0 (as defined above) | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | double rbfcalc3(const rbfmodel &s, const double x0, const double x1, const | |
| | | double x2); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates values of the RBF model at the given point. | |
| | | | |
| | | This is general function which can be used for arbitrary NX (dimension of | |
| | | the space of arguments) and NY (dimension of the function itself). However | |
| | | when you have NY=1 you may find more convenient to use RBFCalc2() or | |
| | | RBFCalc3(). | |
| | | | |
| | | This function returns 0.0 when model is not initialized. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model | |
| | | X - coordinates, array[NX]. | |
| | | X may have more than NX elements, in this case only | |
| | | leading NX will be used. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Y - function value, array[NY]. Y is out-parameter and | |
| | | reallocated after call to this function. In case you want | |
| | | to reuse previously allocated Y, you may use RBFCalcBuf(), | |
| | | which reallocates Y only when it is too small. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfcalc(const rbfmodel &s, const real_1d_array &x, real_1d_array &y); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates values of the RBF model at the given point. | |
| | | | |
| | | Same as RBFCalc(), but does not reallocate Y when in is large enough to | |
| | | store function values. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model | |
| | | X - coordinates, array[NX]. | |
| | | X may have more than NX elements, in this case only | |
| | | leading NX will be used. | |
| | | Y - possibly preallocated array | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Y - function value, array[NY]. Y is not reallocated when it | |
| | | is larger than NY. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfcalcbuf(const rbfmodel &s, const real_1d_array &x, real_1d_array &y | |
| | | ); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates values of the RBF model at the regular grid. | |
| | | | |
| | | Grid have N0*N1 points, with Point[I,J] = (X0[I], X1[J]) | |
| | | | |
| | | This function returns 0.0 when: | |
| | | * model is not initialized | |
| | | * NX<>2 | |
| | | *NY<>1 | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model | |
| | | X0 - array of grid nodes, first coordinates, array[N0] | |
| | | N0 - grid size (number of nodes) in the first dimension | |
| | | X1 - array of grid nodes, second coordinates, array[N1] | |
| | | N1 - grid size (number of nodes) in the second dimension | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Y - function values, array[N0,N1]. Y is out-variable and | |
| | | is reallocated by this function. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfgridcalc2(const rbfmodel &s, const real_1d_array &x0, const ae_int_ | |
| | | t n0, const real_1d_array &x1, const ae_int_t n1, real_2d_array &y); | |
| | | | |
| | | /************************************************************************* | |
| | | This function "unpacks" RBF model by extracting its coefficients. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | S - RBF model | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | NX - dimensionality of argument | |
| | | NY - dimensionality of the target function | |
| | | XWR - model information, array[NC,NX+NY+1]. | |
| | | One row of the array corresponds to one basis function: | |
| | | * first NX columns - coordinates of the center | |
| | | * next NY columns - weights, one per dimension of the | |
| | | function being modelled | |
| | | * last column - radius, same for all dimensions of | |
| | | the function being modelled | |
| | | NC - number of the centers | |
| | | V - polynomial term , array[NY,NX+1]. One row per one | |
| | | dimension of the function being modelled. First NX | |
| | | elements are linear coefficients, V[NX] is equal to the | |
| | | constant part. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void rbfunpack(const rbfmodel &s, ae_int_t &nx, ae_int_t &ny, real_2d_array | |
| | | &xwr, ae_int_t &nc, real_2d_array &v); | |
| | | | |
| | | /************************************************************************* | |
| This subroutine builds bilinear spline coefficients table. | | This subroutine builds bilinear spline coefficients table. | |
| | | | |
| Input parameters: | | Input parameters: | |
| X - spline abscissas, array[0..N-1] | | X - spline abscissas, array[0..N-1] | |
| Y - spline ordinates, array[0..M-1] | | Y - spline ordinates, array[0..M-1] | |
| F - function values, array[0..M-1,0..N-1] | | F - function values, array[0..M-1,0..N-1] | |
| M,N - grid size, M>=2, N>=2 | | M,N - grid size, M>=2, N>=2 | |
| | | | |
| Output parameters: | | Output parameters: | |
| C - spline interpolant | | C - spline interpolant | |
| | | | |
| skipping to change at line 3803 | | skipping to change at line 4404 | |
|---|
| ae_int_t n, | | ae_int_t n, | |
| /* Real */ ae_vector* x2, | | /* Real */ ae_vector* x2, | |
| ae_int_t n2, | | ae_int_t n2, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| ae_bool needy, | | ae_bool needy, | |
| /* Real */ ae_vector* d1, | | /* Real */ ae_vector* d1, | |
| ae_bool needd1, | | ae_bool needd1, | |
| /* Real */ ae_vector* d2, | | /* Real */ ae_vector* d2, | |
| ae_bool needd2, | | ae_bool needd2, | |
| ae_state *_state); | | ae_state *_state); | |
|
| | | void spline1drootsandextrema(spline1dinterpolant* c, | |
| | | /* Real */ ae_vector* r, | |
| | | ae_int_t* nr, | |
| | | ae_bool* dr, | |
| | | /* Real */ ae_vector* e, | |
| | | /* Integer */ ae_vector* et, | |
| | | ae_int_t* ne, | |
| | | ae_bool* de, | |
| | | ae_state *_state); | |
| void heapsortdpoints(/* Real */ ae_vector* x, | | void heapsortdpoints(/* Real */ ae_vector* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| /* Real */ ae_vector* d, | | /* Real */ ae_vector* d, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_state *_state); | | ae_state *_state); | |
|
| | | void solvepolinom2(double p0, | |
| | | double m0, | |
| | | double p1, | |
| | | double m1, | |
| | | double* x0, | |
| | | double* x1, | |
| | | ae_int_t* nr, | |
| | | ae_state *_state); | |
| | | void solvecubicpolinom(double pa, | |
| | | double ma, | |
| | | double pb, | |
| | | double mb, | |
| | | double a, | |
| | | double b, | |
| | | double* x0, | |
| | | double* x1, | |
| | | double* x2, | |
| | | double* ex0, | |
| | | double* ex1, | |
| | | ae_int_t* nr, | |
| | | ae_int_t* ne, | |
| | | /* Real */ ae_vector* tempdata, | |
| | | ae_state *_state); | |
| | | ae_int_t bisectmethod(double pa, | |
| | | double ma, | |
| | | double pb, | |
| | | double mb, | |
| | | double a, | |
| | | double b, | |
| | | double* x, | |
| | | ae_state *_state); | |
| ae_bool _spline1dinterpolant_init(spline1dinterpolant* p, ae_state *_state,
ae_bool make_automatic); | | ae_bool _spline1dinterpolant_init(spline1dinterpolant* p, ae_state *_state,
ae_bool make_automatic); | |
| ae_bool _spline1dinterpolant_init_copy(spline1dinterpolant* dst, spline1din
terpolant* src, ae_state *_state, ae_bool make_automatic); | | ae_bool _spline1dinterpolant_init_copy(spline1dinterpolant* dst, spline1din
terpolant* src, ae_state *_state, ae_bool make_automatic); | |
| void _spline1dinterpolant_clear(spline1dinterpolant* p); | | void _spline1dinterpolant_clear(spline1dinterpolant* p); | |
| void polynomialfit(/* Real */ ae_vector* x, | | void polynomialfit(/* Real */ ae_vector* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_int_t m, | | ae_int_t m, | |
| ae_int_t* info, | | ae_int_t* info, | |
| barycentricinterpolant* p, | | barycentricinterpolant* p, | |
| polynomialfitreport* rep, | | polynomialfitreport* rep, | |
| | | | |
| skipping to change at line 4162 | | skipping to change at line 4803 | |
|---|
| double pspline3arclength(pspline3interpolant* p, | | double pspline3arclength(pspline3interpolant* p, | |
| double a, | | double a, | |
| double b, | | double b, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_bool _pspline2interpolant_init(pspline2interpolant* p, ae_state *_state,
ae_bool make_automatic); | | ae_bool _pspline2interpolant_init(pspline2interpolant* p, ae_state *_state,
ae_bool make_automatic); | |
| ae_bool _pspline2interpolant_init_copy(pspline2interpolant* dst, pspline2in
terpolant* src, ae_state *_state, ae_bool make_automatic); | | ae_bool _pspline2interpolant_init_copy(pspline2interpolant* dst, pspline2in
terpolant* src, ae_state *_state, ae_bool make_automatic); | |
| void _pspline2interpolant_clear(pspline2interpolant* p); | | void _pspline2interpolant_clear(pspline2interpolant* p); | |
| ae_bool _pspline3interpolant_init(pspline3interpolant* p, ae_state *_state,
ae_bool make_automatic); | | ae_bool _pspline3interpolant_init(pspline3interpolant* p, ae_state *_state,
ae_bool make_automatic); | |
| ae_bool _pspline3interpolant_init_copy(pspline3interpolant* dst, pspline3in
terpolant* src, ae_state *_state, ae_bool make_automatic); | | ae_bool _pspline3interpolant_init_copy(pspline3interpolant* dst, pspline3in
terpolant* src, ae_state *_state, ae_bool make_automatic); | |
| void _pspline3interpolant_clear(pspline3interpolant* p); | | void _pspline3interpolant_clear(pspline3interpolant* p); | |
|
| | | void rbfcreate(ae_int_t nx, ae_int_t ny, rbfmodel* s, ae_state *_state); | |
| | | void rbfsetpoints(rbfmodel* s, | |
| | | /* Real */ ae_matrix* xy, | |
| | | ae_int_t n, | |
| | | ae_state *_state); | |
| | | void rbfsetalgoqnn(rbfmodel* s, double q, double z, ae_state *_state); | |
| | | void rbfsetalgomultilayer(rbfmodel* s, | |
| | | double rbase, | |
| | | ae_int_t nlayers, | |
| | | double lambdav, | |
| | | ae_state *_state); | |
| | | void rbfsetlinterm(rbfmodel* s, ae_state *_state); | |
| | | void rbfsetconstterm(rbfmodel* s, ae_state *_state); | |
| | | void rbfsetzeroterm(rbfmodel* s, ae_state *_state); | |
| | | void rbfsetcond(rbfmodel* s, | |
| | | double epsort, | |
| | | double epserr, | |
| | | ae_int_t maxits, | |
| | | ae_state *_state); | |
| | | void rbfbuildmodel(rbfmodel* s, rbfreport* rep, ae_state *_state); | |
| | | double rbfcalc2(rbfmodel* s, double x0, double x1, ae_state *_state); | |
| | | double rbfcalc3(rbfmodel* s, | |
| | | double x0, | |
| | | double x1, | |
| | | double x2, | |
| | | ae_state *_state); | |
| | | void rbfcalc(rbfmodel* s, | |
| | | /* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_state *_state); | |
| | | void rbfcalcbuf(rbfmodel* s, | |
| | | /* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_state *_state); | |
| | | void rbfgridcalc2(rbfmodel* s, | |
| | | /* Real */ ae_vector* x0, | |
| | | ae_int_t n0, | |
| | | /* Real */ ae_vector* x1, | |
| | | ae_int_t n1, | |
| | | /* Real */ ae_matrix* y, | |
| | | ae_state *_state); | |
| | | void rbfunpack(rbfmodel* s, | |
| | | ae_int_t* nx, | |
| | | ae_int_t* ny, | |
| | | /* Real */ ae_matrix* xwr, | |
| | | ae_int_t* nc, | |
| | | /* Real */ ae_matrix* v, | |
| | | ae_state *_state); | |
| | | void rbfalloc(ae_serializer* s, rbfmodel* model, ae_state *_state); | |
| | | void rbfserialize(ae_serializer* s, rbfmodel* model, ae_state *_state); | |
| | | void rbfunserialize(ae_serializer* s, rbfmodel* model, ae_state *_state); | |
| | | ae_bool _rbfmodel_init(rbfmodel* p, ae_state *_state, ae_bool make_automati | |
| | | c); | |
| | | ae_bool _rbfmodel_init_copy(rbfmodel* dst, rbfmodel* src, ae_state *_state, | |
| | | ae_bool make_automatic); | |
| | | void _rbfmodel_clear(rbfmodel* p); | |
| | | ae_bool _rbfreport_init(rbfreport* p, ae_state *_state, ae_bool make_automa | |
| | | tic); | |
| | | ae_bool _rbfreport_init_copy(rbfreport* dst, rbfreport* src, ae_state *_sta | |
| | | te, ae_bool make_automatic); | |
| | | void _rbfreport_clear(rbfreport* p); | |
| void spline2dbuildbilinear(/* Real */ ae_vector* x, | | void spline2dbuildbilinear(/* Real */ ae_vector* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| /* Real */ ae_matrix* f, | | /* Real */ ae_matrix* f, | |
| ae_int_t m, | | ae_int_t m, | |
| ae_int_t n, | | ae_int_t n, | |
| spline2dinterpolant* c, | | spline2dinterpolant* c, | |
| ae_state *_state); | | ae_state *_state); | |
| void spline2dbuildbicubic(/* Real */ ae_vector* x, | | void spline2dbuildbicubic(/* Real */ ae_vector* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| /* Real */ ae_matrix* f, | | /* Real */ ae_matrix* f, | |
| | | | |
| End of changes. 8 change blocks. |
|---|
| 1 lines changed or deleted | | 710 lines changed or added | |
|---|
|
| linalg.h | | linalg.h | |
|---|
| | | | |
| skipping to change at line 55 | | skipping to change at line 55 | |
|---|
| ae_vector rk; | | ae_vector rk; | |
| ae_vector rk1; | | ae_vector rk1; | |
| ae_vector xk; | | ae_vector xk; | |
| ae_vector xk1; | | ae_vector xk1; | |
| ae_vector pk; | | ae_vector pk; | |
| ae_vector pk1; | | ae_vector pk1; | |
| ae_vector b; | | ae_vector b; | |
| rcommstate rstate; | | rcommstate rstate; | |
| ae_vector tmp2; | | ae_vector tmp2; | |
| } fblslincgstate; | | } fblslincgstate; | |
|
| | | typedef struct | |
| | | { | |
| | | ae_vector vals; | |
| | | ae_vector idx; | |
| | | ae_vector ridx; | |
| | | ae_vector didx; | |
| | | ae_vector uidx; | |
| | | ae_int_t matrixtype; | |
| | | ae_int_t m; | |
| | | ae_int_t n; | |
| | | ae_int_t nfree; | |
| | | ae_int_t ninitialized; | |
| | | } sparsematrix; | |
| | | typedef struct | |
| | | { | |
| | | ae_int_t n; | |
| | | ae_int_t m; | |
| | | ae_int_t nstart; | |
| | | ae_int_t nits; | |
| | | ae_int_t seedval; | |
| | | ae_vector x0; | |
| | | ae_vector x1; | |
| | | ae_vector t; | |
| | | ae_vector xbest; | |
| | | hqrndstate r; | |
| | | ae_vector x; | |
| | | ae_vector mv; | |
| | | ae_vector mtv; | |
| | | ae_bool needmv; | |
| | | ae_bool needmtv; | |
| | | double repnorm; | |
| | | rcommstate rstate; | |
| | | } normestimatorstate; | |
| | | | |
| } | | } | |
| | | | |
| ///////////////////////////////////////////////////////////////////////// | | ///////////////////////////////////////////////////////////////////////// | |
| // | | // | |
| // THIS SECTION CONTAINS C++ INTERFACE | | // THIS SECTION CONTAINS C++ INTERFACE | |
| // | | // | |
| ///////////////////////////////////////////////////////////////////////// | | ///////////////////////////////////////////////////////////////////////// | |
| namespace alglib | | namespace alglib | |
| { | | { | |
| | | | |
| skipping to change at line 96 | | skipping to change at line 129 | |
|---|
| matinvreport(); | | matinvreport(); | |
| matinvreport(const matinvreport &rhs); | | matinvreport(const matinvreport &rhs); | |
| matinvreport& operator=(const matinvreport &rhs); | | matinvreport& operator=(const matinvreport &rhs); | |
| virtual ~matinvreport(); | | virtual ~matinvreport(); | |
| double &r1; | | double &r1; | |
| double &rinf; | | double &rinf; | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | Sparse matrix | |
| | | | |
| | | You should use ALGLIB functions to work with sparse matrix. | |
| | | Never try to access its fields directly! | |
| | | *************************************************************************/ | |
| | | class _sparsematrix_owner | |
| | | { | |
| | | public: | |
| | | _sparsematrix_owner(); | |
| | | _sparsematrix_owner(const _sparsematrix_owner &rhs); | |
| | | _sparsematrix_owner& operator=(const _sparsematrix_owner &rhs); | |
| | | virtual ~_sparsematrix_owner(); | |
| | | alglib_impl::sparsematrix* c_ptr(); | |
| | | alglib_impl::sparsematrix* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::sparsematrix *p_struct; | |
| | | }; | |
| | | class sparsematrix : public _sparsematrix_owner | |
| | | { | |
| | | public: | |
| | | sparsematrix(); | |
| | | sparsematrix(const sparsematrix &rhs); | |
| | | sparsematrix& operator=(const sparsematrix &rhs); | |
| | | virtual ~sparsematrix(); | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| | | This object stores state of the iterative norm estimation algorithm. | |
| | | | |
| | | You should use ALGLIB functions to work with this object. | |
| | | *************************************************************************/ | |
| | | class _normestimatorstate_owner | |
| | | { | |
| | | public: | |
| | | _normestimatorstate_owner(); | |
| | | _normestimatorstate_owner(const _normestimatorstate_owner &rhs); | |
| | | _normestimatorstate_owner& operator=(const _normestimatorstate_owner &r | |
| | | hs); | |
| | | virtual ~_normestimatorstate_owner(); | |
| | | alglib_impl::normestimatorstate* c_ptr(); | |
| | | alglib_impl::normestimatorstate* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::normestimatorstate *p_struct; | |
| | | }; | |
| | | class normestimatorstate : public _normestimatorstate_owner | |
| | | { | |
| | | public: | |
| | | normestimatorstate(); | |
| | | normestimatorstate(const normestimatorstate &rhs); | |
| | | normestimatorstate& operator=(const normestimatorstate &rhs); | |
| | | virtual ~normestimatorstate(); | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| Cache-oblivous complex "copy-and-transpose" | | Cache-oblivous complex "copy-and-transpose" | |
| | | | |
| Input parameters: | | Input parameters: | |
| M - number of rows | | M - number of rows | |
| N - number of columns | | N - number of columns | |
| A - source matrix, MxN submatrix is copied and transposed | | A - source matrix, MxN submatrix is copied and transposed | |
| IA - submatrix offset (row index) | | IA - submatrix offset (row index) | |
| JA - submatrix offset (column index) | | JA - submatrix offset (column index) | |
|
| A - destination matrix | | B - destination matrix, must be large enough to store result | |
| IB - submatrix offset (row index) | | IB - submatrix offset (row index) | |
| JB - submatrix offset (column index) | | JB - submatrix offset (column index) | |
| *************************************************************************/ | | *************************************************************************/ | |
| void cmatrixtranspose(const ae_int_t m, const ae_int_t n, const complex_2d_
array &a, const ae_int_t ia, const ae_int_t ja, complex_2d_array &b, const
ae_int_t ib, const ae_int_t jb); | | void cmatrixtranspose(const ae_int_t m, const ae_int_t n, const complex_2d_
array &a, const ae_int_t ia, const ae_int_t ja, complex_2d_array &b, const
ae_int_t ib, const ae_int_t jb); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| Cache-oblivous real "copy-and-transpose" | | Cache-oblivous real "copy-and-transpose" | |
| | | | |
| Input parameters: | | Input parameters: | |
| M - number of rows | | M - number of rows | |
| N - number of columns | | N - number of columns | |
| A - source matrix, MxN submatrix is copied and transposed | | A - source matrix, MxN submatrix is copied and transposed | |
| IA - submatrix offset (row index) | | IA - submatrix offset (row index) | |
| JA - submatrix offset (column index) | | JA - submatrix offset (column index) | |
|
| A - destination matrix | | B - destination matrix, must be large enough to store result | |
| IB - submatrix offset (row index) | | IB - submatrix offset (row index) | |
| JB - submatrix offset (column index) | | JB - submatrix offset (column index) | |
| *************************************************************************/ | | *************************************************************************/ | |
| void rmatrixtranspose(const ae_int_t m, const ae_int_t n, const real_2d_arr
ay &a, const ae_int_t ia, const ae_int_t ja, real_2d_array &b, const ae_int
_t ib, const ae_int_t jb); | | void rmatrixtranspose(const ae_int_t m, const ae_int_t n, const real_2d_arr
ay &a, const ae_int_t ia, const ae_int_t ja, real_2d_array &b, const ae_int
_t ib, const ae_int_t jb); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| Copy | | Copy | |
| | | | |
| Input parameters: | | Input parameters: | |
| M - number of rows | | M - number of rows | |
| N - number of columns | | N - number of columns | |
| A - source matrix, MxN submatrix is copied and transposed | | A - source matrix, MxN submatrix is copied and transposed | |
| IA - submatrix offset (row index) | | IA - submatrix offset (row index) | |
| JA - submatrix offset (column index) | | JA - submatrix offset (column index) | |
|
| B - destination matrix | | B - destination matrix, must be large enough to store result | |
| IB - submatrix offset (row index) | | IB - submatrix offset (row index) | |
| JB - submatrix offset (column index) | | JB - submatrix offset (column index) | |
| *************************************************************************/ | | *************************************************************************/ | |
| void cmatrixcopy(const ae_int_t m, const ae_int_t n, const complex_2d_array
&a, const ae_int_t ia, const ae_int_t ja, complex_2d_array &b, const ae_in
t_t ib, const ae_int_t jb); | | void cmatrixcopy(const ae_int_t m, const ae_int_t n, const complex_2d_array
&a, const ae_int_t ia, const ae_int_t ja, complex_2d_array &b, const ae_in
t_t ib, const ae_int_t jb); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| Copy | | Copy | |
| | | | |
| Input parameters: | | Input parameters: | |
| M - number of rows | | M - number of rows | |
| N - number of columns | | N - number of columns | |
| A - source matrix, MxN submatrix is copied and transposed | | A - source matrix, MxN submatrix is copied and transposed | |
| IA - submatrix offset (row index) | | IA - submatrix offset (row index) | |
| JA - submatrix offset (column index) | | JA - submatrix offset (column index) | |
|
| B - destination matrix | | B - destination matrix, must be large enough to store result | |
| IB - submatrix offset (row index) | | IB - submatrix offset (row index) | |
| JB - submatrix offset (column index) | | JB - submatrix offset (column index) | |
| *************************************************************************/ | | *************************************************************************/ | |
| void rmatrixcopy(const ae_int_t m, const ae_int_t n, const real_2d_array &a
, const ae_int_t ia, const ae_int_t ja, real_2d_array &b, const ae_int_t ib
, const ae_int_t jb); | | void rmatrixcopy(const ae_int_t m, const ae_int_t n, const real_2d_array &a
, const ae_int_t ia, const ae_int_t ja, real_2d_array &b, const ae_int_t ib
, const ae_int_t jb); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| Rank-1 correction: A := A + u*v' | | Rank-1 correction: A := A + u*v' | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| M - number of rows | | M - number of rows | |
| | | | |
| skipping to change at line 205 | | skipping to change at line 293 | |
|---|
| A - target matrix | | A - target matrix | |
| IA - submatrix offset (row index) | | IA - submatrix offset (row index) | |
| JA - submatrix offset (column index) | | JA - submatrix offset (column index) | |
| OpA - operation type: | | OpA - operation type: | |
| * OpA=0 => op(A) = A | | * OpA=0 => op(A) = A | |
| * OpA=1 => op(A) = A^T | | * OpA=1 => op(A) = A^T | |
| * OpA=2 => op(A) = A^H | | * OpA=2 => op(A) = A^H | |
| X - input vector | | X - input vector | |
| IX - subvector offset | | IX - subvector offset | |
| IY - subvector offset | | IY - subvector offset | |
|
| | | Y - preallocated matrix, must be large enough to store result | |
| | | | |
| OUTPUT PARAMETERS: | | OUTPUT PARAMETERS: | |
| Y - vector which stores result | | Y - vector which stores result | |
| | | | |
| if M=0, then subroutine does nothing. | | if M=0, then subroutine does nothing. | |
| if N=0, Y is filled by zeros. | | if N=0, Y is filled by zeros. | |
| | | | |
| -- ALGLIB routine -- | | -- ALGLIB routine -- | |
| | | | |
| 28.01.2010 | | 28.01.2010 | |
| | | | |
| skipping to change at line 234 | | skipping to change at line 323 | |
|---|
| N - number of columns of op(A) | | N - number of columns of op(A) | |
| A - target matrix | | A - target matrix | |
| IA - submatrix offset (row index) | | IA - submatrix offset (row index) | |
| JA - submatrix offset (column index) | | JA - submatrix offset (column index) | |
| OpA - operation type: | | OpA - operation type: | |
| * OpA=0 => op(A) = A | | * OpA=0 => op(A) = A | |
| * OpA=1 => op(A) = A^T | | * OpA=1 => op(A) = A^T | |
| X - input vector | | X - input vector | |
| IX - subvector offset | | IX - subvector offset | |
| IY - subvector offset | | IY - subvector offset | |
|
| | | Y - preallocated matrix, must be large enough to store result | |
| | | | |
| OUTPUT PARAMETERS: | | OUTPUT PARAMETERS: | |
| Y - vector which stores result | | Y - vector which stores result | |
| | | | |
| if M=0, then subroutine does nothing. | | if M=0, then subroutine does nothing. | |
| if N=0, Y is filled by zeros. | | if N=0, Y is filled by zeros. | |
| | | | |
| -- ALGLIB routine -- | | -- ALGLIB routine -- | |
| | | | |
| 28.01.2010 | | 28.01.2010 | |
| | | | |
| skipping to change at line 269 | | skipping to change at line 359 | |
|---|
| M - matrix size, N>=0 | | M - matrix size, N>=0 | |
| A - matrix, actial matrix is stored in A[I1:I1+N-1,J1:J1+N-1] | | A - matrix, actial matrix is stored in A[I1:I1+N-1,J1:J1+N-1] | |
| I1 - submatrix offset | | I1 - submatrix offset | |
| J1 - submatrix offset | | J1 - submatrix offset | |
| IsUpper - whether matrix is upper triangular | | IsUpper - whether matrix is upper triangular | |
| IsUnit - whether matrix is unitriangular | | IsUnit - whether matrix is unitriangular | |
| OpType - transformation type: | | OpType - transformation type: | |
| * 0 - no transformation | | * 0 - no transformation | |
| * 1 - transposition | | * 1 - transposition | |
| * 2 - conjugate transposition | | * 2 - conjugate transposition | |
|
| C - matrix, actial matrix is stored in C[I2:I2+M-1,J2:J2+N-1] | | X - matrix, actial matrix is stored in X[I2:I2+M-1,J2:J2+N-1] | |
| I2 - submatrix offset | | I2 - submatrix offset | |
| J2 - submatrix offset | | J2 - submatrix offset | |
| | | | |
| -- ALGLIB routine -- | | -- ALGLIB routine -- | |
| 15.12.2009 | | 15.12.2009 | |
| Bochkanov Sergey | | Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void cmatrixrighttrsm(const ae_int_t m, const ae_int_t n, const complex_2d_
array &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const b
ool isunit, const ae_int_t optype, complex_2d_array &x, const ae_int_t i2,
const ae_int_t j2); | | void cmatrixrighttrsm(const ae_int_t m, const ae_int_t n, const complex_2d_
array &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const b
ool isunit, const ae_int_t optype, complex_2d_array &x, const ae_int_t i2,
const ae_int_t j2); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| | | | |
| skipping to change at line 300 | | skipping to change at line 390 | |
|---|
| M - matrix size, N>=0 | | M - matrix size, N>=0 | |
| A - matrix, actial matrix is stored in A[I1:I1+M-1,J1:J1+M-1] | | A - matrix, actial matrix is stored in A[I1:I1+M-1,J1:J1+M-1] | |
| I1 - submatrix offset | | I1 - submatrix offset | |
| J1 - submatrix offset | | J1 - submatrix offset | |
| IsUpper - whether matrix is upper triangular | | IsUpper - whether matrix is upper triangular | |
| IsUnit - whether matrix is unitriangular | | IsUnit - whether matrix is unitriangular | |
| OpType - transformation type: | | OpType - transformation type: | |
| * 0 - no transformation | | * 0 - no transformation | |
| * 1 - transposition | | * 1 - transposition | |
| * 2 - conjugate transposition | | * 2 - conjugate transposition | |
|
| C - matrix, actial matrix is stored in C[I2:I2+M-1,J2:J2+N-1] | | X - matrix, actial matrix is stored in X[I2:I2+M-1,J2:J2+N-1] | |
| I2 - submatrix offset | | I2 - submatrix offset | |
| J2 - submatrix offset | | J2 - submatrix offset | |
| | | | |
| -- ALGLIB routine -- | | -- ALGLIB routine -- | |
| 15.12.2009 | | 15.12.2009 | |
| Bochkanov Sergey | | Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void cmatrixlefttrsm(const ae_int_t m, const ae_int_t n, const complex_2d_a
rray &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bo
ol isunit, const ae_int_t optype, complex_2d_array &x, const ae_int_t i2, c
onst ae_int_t j2); | | void cmatrixlefttrsm(const ae_int_t m, const ae_int_t n, const complex_2d_a
rray &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bo
ol isunit, const ae_int_t optype, complex_2d_array &x, const ae_int_t i2, c
onst ae_int_t j2); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Same as CMatrixRightTRSM, but for real matrices | | This subroutine calculates X*op(A^-1) where: | |
| | | * X is MxN general matrix | |
| | | * A is NxN upper/lower triangular/unitriangular matrix | |
| | | * "op" may be identity transformation, transposition | |
| | | | |
| | | Multiplication result replaces X. | |
| | | Cache-oblivious algorithm is used. | |
| | | | |
|
| OpType may be only 0 or 1. | | INPUT PARAMETERS | |
| | | N - matrix size, N>=0 | |
| | | M - matrix size, N>=0 | |
| | | A - matrix, actial matrix is stored in A[I1:I1+N-1,J1:J1+N-1] | |
| | | I1 - submatrix offset | |
| | | J1 - submatrix offset | |
| | | IsUpper - whether matrix is upper triangular | |
| | | IsUnit - whether matrix is unitriangular | |
| | | OpType - transformation type: | |
| | | * 0 - no transformation | |
| | | * 1 - transposition | |
| | | X - matrix, actial matrix is stored in X[I2:I2+M-1,J2:J2+N-1] | |
| | | I2 - submatrix offset | |
| | | J2 - submatrix offset | |
| | | | |
| -- ALGLIB routine -- | | -- ALGLIB routine -- | |
| 15.12.2009 | | 15.12.2009 | |
| Bochkanov Sergey | | Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void rmatrixrighttrsm(const ae_int_t m, const ae_int_t n, const real_2d_arr
ay &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bool
isunit, const ae_int_t optype, real_2d_array &x, const ae_int_t i2, const
ae_int_t j2); | | void rmatrixrighttrsm(const ae_int_t m, const ae_int_t n, const real_2d_arr
ay &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bool
isunit, const ae_int_t optype, real_2d_array &x, const ae_int_t i2, const
ae_int_t j2); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Same as CMatrixLeftTRSM, but for real matrices | | This subroutine calculates op(A^-1)*X where: | |
| | | * X is MxN general matrix | |
| | | * A is MxM upper/lower triangular/unitriangular matrix | |
| | | * "op" may be identity transformation, transposition | |
| | | | |
| | | Multiplication result replaces X. | |
| | | Cache-oblivious algorithm is used. | |
| | | | |
|
| OpType may be only 0 or 1. | | INPUT PARAMETERS | |
| | | N - matrix size, N>=0 | |
| | | M - matrix size, N>=0 | |
| | | A - matrix, actial matrix is stored in A[I1:I1+M-1,J1:J1+M-1] | |
| | | I1 - submatrix offset | |
| | | J1 - submatrix offset | |
| | | IsUpper - whether matrix is upper triangular | |
| | | IsUnit - whether matrix is unitriangular | |
| | | OpType - transformation type: | |
| | | * 0 - no transformation | |
| | | * 1 - transposition | |
| | | X - matrix, actial matrix is stored in X[I2:I2+M-1,J2:J2+N-1] | |
| | | I2 - submatrix offset | |
| | | J2 - submatrix offset | |
| | | | |
| -- ALGLIB routine -- | | -- ALGLIB routine -- | |
| 15.12.2009 | | 15.12.2009 | |
| Bochkanov Sergey | | Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void rmatrixlefttrsm(const ae_int_t m, const ae_int_t n, const real_2d_arra
y &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bool
isunit, const ae_int_t optype, real_2d_array &x, const ae_int_t i2, const a
e_int_t j2); | | void rmatrixlefttrsm(const ae_int_t m, const ae_int_t n, const real_2d_arra
y &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bool
isunit, const ae_int_t optype, real_2d_array &x, const ae_int_t i2, const a
e_int_t j2); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This subroutine calculates C=alpha*A*A^H+beta*C or C=alpha*A^H*A+beta*C | | This subroutine calculates C=alpha*A*A^H+beta*C or C=alpha*A^H*A+beta*C | |
| where: | | where: | |
| | | | |
| skipping to change at line 368 | | skipping to change at line 496 | |
|---|
| JC - submatrix offset | | JC - submatrix offset | |
| IsUpper - whether C is upper triangular or lower triangular | | IsUpper - whether C is upper triangular or lower triangular | |
| | | | |
| -- ALGLIB routine -- | | -- ALGLIB routine -- | |
| 16.12.2009 | | 16.12.2009 | |
| Bochkanov Sergey | | Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void cmatrixsyrk(const ae_int_t n, const ae_int_t k, const double alpha, co
nst complex_2d_array &a, const ae_int_t ia, const ae_int_t ja, const ae_int
_t optypea, const double beta, complex_2d_array &c, const ae_int_t ic, cons
t ae_int_t jc, const bool isupper); | | void cmatrixsyrk(const ae_int_t n, const ae_int_t k, const double alpha, co
nst complex_2d_array &a, const ae_int_t ia, const ae_int_t ja, const ae_int
_t optypea, const double beta, complex_2d_array &c, const ae_int_t ic, cons
t ae_int_t jc, const bool isupper); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Same as CMatrixSYRK, but for real matrices | | This subroutine calculates C=alpha*A*A^T+beta*C or C=alpha*A^T*A+beta*C | |
| | | where: | |
| | | * C is NxN symmetric matrix given by its upper/lower triangle | |
| | | * A is NxK matrix when A*A^T is calculated, KxN matrix otherwise | |
| | | | |
| | | Additional info: | |
| | | * cache-oblivious algorithm is used. | |
| | | * multiplication result replaces C. If Beta=0, C elements are not used in | |
| | | calculations (not multiplied by zero - just not referenced) | |
| | | * if Alpha=0, A is not used (not multiplied by zero - just not referenced) | |
| | | * if both Beta and Alpha are zero, C is filled by zeros. | |
| | | | |
|
| OpType may be only 0 or 1. | | INPUT PARAMETERS | |
| | | N - matrix size, N>=0 | |
| | | K - matrix size, K>=0 | |
| | | Alpha - coefficient | |
| | | A - matrix | |
| | | IA - submatrix offset | |
| | | JA - submatrix offset | |
| | | OpTypeA - multiplication type: | |
| | | * 0 - A*A^T is calculated | |
| | | * 2 - A^T*A is calculated | |
| | | Beta - coefficient | |
| | | C - matrix | |
| | | IC - submatrix offset | |
| | | JC - submatrix offset | |
| | | IsUpper - whether C is upper triangular or lower triangular | |
| | | | |
| -- ALGLIB routine -- | | -- ALGLIB routine -- | |
| 16.12.2009 | | 16.12.2009 | |
| Bochkanov Sergey | | Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void rmatrixsyrk(const ae_int_t n, const ae_int_t k, const double alpha, co
nst real_2d_array &a, const ae_int_t ia, const ae_int_t ja, const ae_int_t
optypea, const double beta, real_2d_array &c, const ae_int_t ic, const ae_i
nt_t jc, const bool isupper); | | void rmatrixsyrk(const ae_int_t n, const ae_int_t k, const double alpha, co
nst real_2d_array &a, const ae_int_t ia, const ae_int_t ja, const ae_int_t
optypea, const double beta, real_2d_array &c, const ae_int_t ic, const ae_i
nt_t jc, const bool isupper); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This subroutine calculates C = alpha*op1(A)*op2(B) +beta*C where: | | This subroutine calculates C = alpha*op1(A)*op2(B) +beta*C where: | |
| * C is MxN general matrix | | * C is MxN general matrix | |
| | | | |
| skipping to change at line 393 | | skipping to change at line 545 | |
|---|
| * "op" may be identity transformation, transposition, conjugate transpositi
on | | * "op" may be identity transformation, transposition, conjugate transpositi
on | |
| | | | |
| Additional info: | | Additional info: | |
| * cache-oblivious algorithm is used. | | * cache-oblivious algorithm is used. | |
| * multiplication result replaces C. If Beta=0, C elements are not used in | | * multiplication result replaces C. If Beta=0, C elements are not used in | |
| calculations (not multiplied by zero - just not referenced) | | calculations (not multiplied by zero - just not referenced) | |
| * if Alpha=0, A is not used (not multiplied by zero - just not referenced) | | * if Alpha=0, A is not used (not multiplied by zero - just not referenced) | |
| * if both Beta and Alpha are zero, C is filled by zeros. | | * if both Beta and Alpha are zero, C is filled by zeros. | |
| | | | |
| INPUT PARAMETERS | | INPUT PARAMETERS | |
|
| | | M - matrix size, M>0 | |
| N - matrix size, N>0 | | N - matrix size, N>0 | |
|
| M - matrix size, N>0 | | | |
| K - matrix size, K>0 | | K - matrix size, K>0 | |
| Alpha - coefficient | | Alpha - coefficient | |
| A - matrix | | A - matrix | |
| IA - submatrix offset | | IA - submatrix offset | |
| JA - submatrix offset | | JA - submatrix offset | |
| OpTypeA - transformation type: | | OpTypeA - transformation type: | |
| * 0 - no transformation | | * 0 - no transformation | |
| * 1 - transposition | | * 1 - transposition | |
| * 2 - conjugate transposition | | * 2 - conjugate transposition | |
| B - matrix | | B - matrix | |
| | | | |
| skipping to change at line 423 | | skipping to change at line 575 | |
|---|
| IC - submatrix offset | | IC - submatrix offset | |
| JC - submatrix offset | | JC - submatrix offset | |
| | | | |
| -- ALGLIB routine -- | | -- ALGLIB routine -- | |
| 16.12.2009 | | 16.12.2009 | |
| Bochkanov Sergey | | Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void cmatrixgemm(const ae_int_t m, const ae_int_t n, const ae_int_t k, cons
t alglib::complex alpha, const complex_2d_array &a, const ae_int_t ia, cons
t ae_int_t ja, const ae_int_t optypea, const complex_2d_array &b, const ae_
int_t ib, const ae_int_t jb, const ae_int_t optypeb, const alglib::complex
beta, complex_2d_array &c, const ae_int_t ic, const ae_int_t jc); | | void cmatrixgemm(const ae_int_t m, const ae_int_t n, const ae_int_t k, cons
t alglib::complex alpha, const complex_2d_array &a, const ae_int_t ia, cons
t ae_int_t ja, const ae_int_t optypea, const complex_2d_array &b, const ae_
int_t ib, const ae_int_t jb, const ae_int_t optypeb, const alglib::complex
beta, complex_2d_array &c, const ae_int_t ic, const ae_int_t jc); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Same as CMatrixGEMM, but for real numbers. | | This subroutine calculates C = alpha*op1(A)*op2(B) +beta*C where: | |
| OpType may be only 0 or 1. | | * C is MxN general matrix | |
| | | * op1(A) is MxK matrix | |
| | | * op2(B) is KxN matrix | |
| | | * "op" may be identity transformation, transposition | |
| | | | |
| | | Additional info: | |
| | | * cache-oblivious algorithm is used. | |
| | | * multiplication result replaces C. If Beta=0, C elements are not used in | |
| | | calculations (not multiplied by zero - just not referenced) | |
| | | * if Alpha=0, A is not used (not multiplied by zero - just not referenced) | |
| | | * if both Beta and Alpha are zero, C is filled by zeros. | |
| | | | |
| | | INPUT PARAMETERS | |
| | | M - matrix size, M>0 | |
| | | N - matrix size, N>0 | |
| | | K - matrix size, K>0 | |
| | | Alpha - coefficient | |
| | | A - matrix | |
| | | IA - submatrix offset | |
| | | JA - submatrix offset | |
| | | OpTypeA - transformation type: | |
| | | * 0 - no transformation | |
| | | * 1 - transposition | |
| | | B - matrix | |
| | | IB - submatrix offset | |
| | | JB - submatrix offset | |
| | | OpTypeB - transformation type: | |
| | | * 0 - no transformation | |
| | | * 1 - transposition | |
| | | Beta - coefficient | |
| | | C - matrix | |
| | | IC - submatrix offset | |
| | | JC - submatrix offset | |
| | | | |
| -- ALGLIB routine -- | | -- ALGLIB routine -- | |
| 16.12.2009 | | 16.12.2009 | |
| Bochkanov Sergey | | Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void rmatrixgemm(const ae_int_t m, const ae_int_t n, const ae_int_t k, cons
t double alpha, const real_2d_array &a, const ae_int_t ia, const ae_int_t j
a, const ae_int_t optypea, const real_2d_array &b, const ae_int_t ib, const
ae_int_t jb, const ae_int_t optypeb, const double beta, real_2d_array &c,
const ae_int_t ic, const ae_int_t jc); | | void rmatrixgemm(const ae_int_t m, const ae_int_t n, const ae_int_t k, cons
t double alpha, const real_2d_array &a, const ae_int_t ia, const ae_int_t j
a, const ae_int_t optypea, const real_2d_array &b, const ae_int_t ib, const
ae_int_t jb, const ae_int_t optypeb, const double beta, real_2d_array &c,
const ae_int_t ic, const ae_int_t jc); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| QR decomposition of a rectangular matrix of size MxN | | QR decomposition of a rectangular matrix of size MxN | |
| | | | |
| | | | |
| skipping to change at line 1162 | | skipping to change at line 1346 | |
|---|
| Output parameters: | | Output parameters: | |
| Q - transformation matrix. | | Q - transformation matrix. | |
| array with elements [0..N-1, 0..N-1]. | | array with elements [0..N-1, 0..N-1]. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 2005-2010 by Bochkanov Sergey | | Copyright 2005-2010 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void hmatrixtdunpackq(const complex_2d_array &a, const ae_int_t n, const bo
ol isupper, const complex_1d_array &tau, complex_2d_array &q); | | void hmatrixtdunpackq(const complex_2d_array &a, const ae_int_t n, const bo
ol isupper, const complex_1d_array &tau, complex_2d_array &q); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | Singular value decomposition of a bidiagonal matrix (extended algorithm) | |
| | | | |
| | | The algorithm performs the singular value decomposition of a bidiagonal | |
| | | matrix B (upper or lower) representing it as B = Q*S*P^T, where Q and P - | |
| | | orthogonal matrices, S - diagonal matrix with non-negative elements on the | |
| | | main diagonal, in descending order. | |
| | | | |
| | | The algorithm finds singular values. In addition, the algorithm can | |
| | | calculate matrices Q and P (more precisely, not the matrices, but their | |
| | | product with given matrices U and VT - U*Q and (P^T)*VT)). Of course, | |
| | | matrices U and VT can be of any type, including identity. Furthermore, the | |
| | | algorithm can calculate Q'*C (this product is calculated more effectively | |
| | | than U*Q, because this calculation operates with rows instead of matrix | |
| | | columns). | |
| | | | |
| | | The feature of the algorithm is its ability to find all singular values | |
| | | including those which are arbitrarily close to 0 with relative accuracy | |
| | | close to machine precision. If the parameter IsFractionalAccuracyRequired | |
| | | is set to True, all singular values will have high relative accuracy close | |
| | | to machine precision. If the parameter is set to False, only the biggest | |
| | | singular value will have relative accuracy close to machine precision. | |
| | | The absolute error of other singular values is equal to the absolute error | |
| | | of the biggest singular value. | |
| | | | |
| | | Input parameters: | |
| | | D - main diagonal of matrix B. | |
| | | Array whose index ranges within [0..N-1]. | |
| | | E - superdiagonal (or subdiagonal) of matrix B. | |
| | | Array whose index ranges within [0..N-2]. | |
| | | N - size of matrix B. | |
| | | IsUpper - True, if the matrix is upper bidiagonal. | |
| | | IsFractionalAccuracyRequired - | |
| | | THIS PARAMETER IS IGNORED SINCE ALGLIB 3.5.0 | |
| | | SINGULAR VALUES ARE ALWAYS SEARCHED WITH HIGH ACCURACY. | |
| | | U - matrix to be multiplied by Q. | |
| | | Array whose indexes range within [0..NRU-1, 0..N-1]. | |
| | | The matrix can be bigger, in that case only the submatrix | |
| | | [0..NRU-1, 0..N-1] will be multiplied by Q. | |
| | | NRU - number of rows in matrix U. | |
| | | C - matrix to be multiplied by Q'. | |
| | | Array whose indexes range within [0..N-1, 0..NCC-1]. | |
| | | The matrix can be bigger, in that case only the submatrix | |
| | | [0..N-1, 0..NCC-1] will be multiplied by Q'. | |
| | | NCC - number of columns in matrix C. | |
| | | VT - matrix to be multiplied by P^T. | |
| | | Array whose indexes range within [0..N-1, 0..NCVT-1]. | |
| | | The matrix can be bigger, in that case only the submatrix | |
| | | [0..N-1, 0..NCVT-1] will be multiplied by P^T. | |
| | | NCVT - number of columns in matrix VT. | |
| | | | |
| | | Output parameters: | |
| | | D - singular values of matrix B in descending order. | |
| | | U - if NRU>0, contains matrix U*Q. | |
| | | VT - if NCVT>0, contains matrix (P^T)*VT. | |
| | | C - if NCC>0, contains matrix Q'*C. | |
| | | | |
| | | Result: | |
| | | True, if the algorithm has converged. | |
| | | False, if the algorithm hasn't converged (rare case). | |
| | | | |
| | | Additional information: | |
| | | The type of convergence is controlled by the internal parameter TOL. | |
| | | If the parameter is greater than 0, the singular values will have | |
| | | relative accuracy TOL. If TOL<0, the singular values will have | |
| | | absolute accuracy ABS(TOL)*norm(B). | |
| | | By default, |TOL| falls within the range of 10*Epsilon and 100*Epsilon, | |
| | | where Epsilon is the machine precision. It is not recommended to use | |
| | | TOL less than 10*Epsilon since this will considerably slow down the | |
| | | algorithm and may not lead to error decreasing. | |
| | | History: | |
| | | * 31 March, 2007. | |
| | | changed MAXITR from 6 to 12. | |
| | | | |
| | | -- LAPACK routine (version 3.0) -- | |
| | | Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., | |
| | | Courant Institute, Argonne National Lab, and Rice University | |
| | | October 31, 1999. | |
| | | *************************************************************************/ | |
| | | bool rmatrixbdsvd(real_1d_array &d, const real_1d_array &e, const ae_int_t | |
| | | n, const bool isupper, const bool isfractionalaccuracyrequired, real_2d_arr | |
| | | ay &u, const ae_int_t nru, real_2d_array &c, const ae_int_t ncc, real_2d_ar | |
| | | ray &vt, const ae_int_t ncvt); | |
| | | | |
| | | /************************************************************************* | |
| | | Singular value decomposition of a rectangular matrix. | |
| | | | |
| | | The algorithm calculates the singular value decomposition of a matrix of | |
| | | size MxN: A = U * S * V^T | |
| | | | |
| | | The algorithm finds the singular values and, optionally, matrices U and V^T | |
| | | . | |
| | | The algorithm can find both first min(M,N) columns of matrix U and rows of | |
| | | matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM | |
| | | and NxN respectively). | |
| | | | |
| | | Take into account that the subroutine does not return matrix V but V^T. | |
| | | | |
| | | Input parameters: | |
| | | A - matrix to be decomposed. | |
| | | Array whose indexes range within [0..M-1, 0..N-1]. | |
| | | M - number of rows in matrix A. | |
| | | N - number of columns in matrix A. | |
| | | UNeeded - 0, 1 or 2. See the description of the parameter U. | |
| | | VTNeeded - 0, 1 or 2. See the description of the parameter VT. | |
| | | AdditionalMemory - | |
| | | If the parameter: | |
| | | * equals 0, the algorithm doesn | |
| | | memory (lower requirements, lower performance). | |
| | | * equals 1, the algorithm uses additional | |
| | | memory of size min(M,N)*min(M,N) of real numbers. | |
| | | It often speeds up the algorithm. | |
| | | * equals 2, the algorithm uses additional | |
| | | memory of size M*min(M,N) of real numbers. | |
| | | It allows to get a maximum performance. | |
| | | The recommended value of the parameter is 2. | |
| | | | |
| | | Output parameters: | |
| | | W - contains singular values in descending order. | |
| | | U - if UNeeded=0, U isn't changed, the left singular vector | |
| | | s | |
| | | are not calculated. | |
| | | if Uneeded=1, U contains left singular vectors (first | |
| | | min(M,N) columns of matrix U). Array whose indexes rang | |
| | | e | |
| | | within [0..M-1, 0..Min(M,N)-1]. | |
| | | if UNeeded=2, U contains matrix U wholly. Array whose | |
| | | indexes range within [0..M-1, 0..M-1]. | |
| | | VT - if VTNeeded=0, VT isn | |
| | | are not calculated. | |
| | | if VTNeeded=1, VT contains right singular vectors (firs | |
| | | t | |
| | | min(M,N) rows of matrix V^T). Array whose indexes range | |
| | | within [0..min(M,N)-1, 0..N-1]. | |
| | | if VTNeeded=2, VT contains matrix V^T wholly. Array who | |
| | | se | |
| | | indexes range within [0..N-1, 0..N-1]. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 2005 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | bool rmatrixsvd(const real_2d_array &a, const ae_int_t m, const ae_int_t n, | |
| | | const ae_int_t uneeded, const ae_int_t vtneeded, const ae_int_t additional | |
| | | memory, real_1d_array &w, real_2d_array &u, real_2d_array &vt); | |
| | | | |
| | | /************************************************************************* | |
| Finding the eigenvalues and eigenvectors of a symmetric matrix | | Finding the eigenvalues and eigenvectors of a symmetric matrix | |
| | | | |
| The algorithm finds eigen pairs of a symmetric matrix by reducing it to | | The algorithm finds eigen pairs of a symmetric matrix by reducing it to | |
| tridiagonal form and using the QL/QR algorithm. | | tridiagonal form and using the QL/QR algorithm. | |
| | | | |
| Input parameters: | | Input parameters: | |
| A - symmetric matrix which is given by its upper or lower | | A - symmetric matrix which is given by its upper or lower | |
| triangular part. | | triangular part. | |
| Array whose indexes range within [0..N-1, 0..N-1]. | | Array whose indexes range within [0..N-1, 0..N-1]. | |
| N - size of matrix A. | | N - size of matrix A. | |
| | | | |
| skipping to change at line 2679 | | skipping to change at line 2998 | |
|---|
| Rep - same as for RMatrixLUInverse | | Rep - same as for RMatrixLUInverse | |
| A - same as for RMatrixLUInverse. | | A - same as for RMatrixLUInverse. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 05.02.2010 by Bochkanov Sergey | | Copyright 05.02.2010 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void cmatrixtrinverse(complex_2d_array &a, const ae_int_t n, const bool isu
pper, const bool isunit, ae_int_t &info, matinvreport &rep); | | void cmatrixtrinverse(complex_2d_array &a, const ae_int_t n, const bool isu
pper, const bool isunit, ae_int_t &info, matinvreport &rep); | |
| void cmatrixtrinverse(complex_2d_array &a, const bool isupper, ae_int_t &in
fo, matinvreport &rep); | | void cmatrixtrinverse(complex_2d_array &a, const bool isupper, ae_int_t &in
fo, matinvreport &rep); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Singular value decomposition of a bidiagonal matrix (extended algorithm) | | This function creates sparse matrix in a Hash-Table format. | |
| | | | |
|
| The algorithm performs the singular value decomposition of a bidiagonal | | This function creates Hast-Table matrix, which can be converted to CRS | |
| matrix B (upper or lower) representing it as B = Q*S*P^T, where Q and P - | | format after its initialization is over. Typical usage scenario for a | |
| orthogonal matrices, S - diagonal matrix with non-negative elements on the | | sparse matrix is: | |
| main diagonal, in descending order. | | 1. creation in a Hash-Table format | |
| | | 2. insertion of the matrix elements | |
| | | 3. conversion to the CRS representation | |
| | | 4. matrix is passed to some linear algebra algorithm | |
| | | | |
|
| The algorithm finds singular values. In addition, the algorithm can | | Some information about different matrix formats can be found below, in | |
| calculate matrices Q and P (more precisely, not the matrices, but their | | the "NOTES" section. | |
| product with given matrices U and VT - U*Q and (P^T)*VT)). Of course, | | | |
| matrices U and VT can be of any type, including identity. Furthermore, the | | | |
| algorithm can calculate Q'*C (this product is calculated more effectively | | | |
| than U*Q, because this calculation operates with rows instead of matrix | | | |
| columns). | | | |
| | | | |
|
| The feature of the algorithm is its ability to find all singular values | | INPUT PARAMETERS | |
| including those which are arbitrarily close to 0 with relative accuracy | | M - number of rows in a matrix, M>=1 | |
| close to machine precision. If the parameter IsFractionalAccuracyRequired | | N - number of columns in a matrix, N>=1 | |
| is set to True, all singular values will have high relative accuracy close | | K - K>=0, expected number of non-zero elements in a matrix. | |
| to machine precision. If the parameter is set to False, only the biggest | | K can be inexact approximation, can be less than actual | |
| singular value will have relative accuracy close to machine precision. | | number of elements (table will grow when needed) or | |
| The absolute error of other singular values is equal to the absolute error | | even zero). | |
| of the biggest singular value. | | It is important to understand that although hash-table | |
| | | may grow automatically, it is better to provide good | |
| | | estimate of data size. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | S - sparse M*N matrix in Hash-Table representation. | |
| | | All elements of the matrix are zero. | |
| | | | |
| | | NOTE 1. | |
| | | | |
| | | Sparse matrices can be stored using either Hash-Table representation or | |
| | | Compressed Row Storage representation. Hast-table is better suited for | |
| | | querying and dynamic operations (thus, it is used for matrix | |
| | | initialization), but it is inefficient when you want to make some linear | |
| | | algebra operations. | |
| | | | |
| | | From the other side, CRS is better suited for linear algebra operations, | |
| | | but initialization is less convenient - you have to tell row sizes at the | |
| | | initialization, and you can fill matrix only row by row, from left to | |
| | | right. CRS is also very inefficient when you want to find matrix element | |
| | | by its index. | |
| | | | |
| | | Thus, Hash-Table representation does not support linear algebra | |
| | | operations, while CRS format does not support modification of the table. | |
| | | Tables below outline information about these two formats: | |
| | | | |
| | | OPERATIONS WITH MATRIX HASH CRS | |
| | | create + + | |
| | | read element + + | |
| | | modify element + | |
| | | add value to element + | |
| | | A*x (dense vector) + | |
| | | A'*x (dense vector) + | |
| | | A*X (dense matrix) + | |
| | | A'*X (dense matrix) + | |
| | | | |
| | | NOTE 2. | |
| | | | |
| | | Hash-tables use memory inefficiently, and they have to keep some amount | |
| | | of the "spare memory" in order to have good performance. Hash table for | |
| | | matrix with K non-zero elements will need C*K*(8+2*sizeof(int)) bytes, | |
| | | where C is a small constant, about 1.5-2 in magnitude. | |
| | | | |
| | | CRS storage, from the other side, is more memory-efficient, and needs | |
| | | just K*(8+sizeof(int))+M*sizeof(int) bytes, where M is a number of rows | |
| | | in a matrix. | |
| | | | |
| | | When you convert from the Hash-Table to CRS representation, all unneeded | |
| | | memory will be freed. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparsecreate(const ae_int_t m, const ae_int_t n, const ae_int_t k, spa | |
| | | rsematrix &s); | |
| | | void sparsecreate(const ae_int_t m, const ae_int_t n, sparsematrix &s); | |
| | | | |
| | | /************************************************************************* | |
| | | This function creates sparse matrix in a CRS format (expert function for | |
| | | situations when you are running out of memory). | |
| | | | |
| | | This function creates CRS matrix. Typical usage scenario for a CRS matrix | |
| | | is: | |
| | | 1. creation (you have to tell number of non-zero elements at each row at | |
| | | this moment) | |
| | | 2. insertion of the matrix elements (row by row, from left to right) | |
| | | 3. matrix is passed to some linear algebra algorithm | |
| | | | |
| | | This function is a memory-efficient alternative to SparseCreate(), but it | |
| | | is more complex because it requires you to know in advance how large your | |
| | | matrix is. Some information about different matrix formats can be found | |
| | | below, in the "NOTES" section. | |
| | | | |
|
| Input parameters: | | INPUT PARAMETERS | |
| D - main diagonal of matrix B. | | M - number of rows in a matrix, M>=1 | |
| Array whose index ranges within [0..N-1]. | | N - number of columns in a matrix, N>=1 | |
| E - superdiagonal (or subdiagonal) of matrix B. | | NER - number of elements at each row, array[M], NER[i]>=0 | |
| Array whose index ranges within [0..N-2]. | | | |
| N - size of matrix B. | | | |
| IsUpper - True, if the matrix is upper bidiagonal. | | | |
| IsFractionalAccuracyRequired - | | | |
| accuracy to search singular values with. | | | |
| U - matrix to be multiplied by Q. | | | |
| Array whose indexes range within [0..NRU-1, 0..N-1]. | | | |
| The matrix can be bigger, in that case only the submatrix | | | |
| [0..NRU-1, 0..N-1] will be multiplied by Q. | | | |
| NRU - number of rows in matrix U. | | | |
| C - matrix to be multiplied by Q'. | | | |
| Array whose indexes range within [0..N-1, 0..NCC-1]. | | | |
| The matrix can be bigger, in that case only the submatrix | | | |
| [0..N-1, 0..NCC-1] will be multiplied by Q'. | | | |
| NCC - number of columns in matrix C. | | | |
| VT - matrix to be multiplied by P^T. | | | |
| Array whose indexes range within [0..N-1, 0..NCVT-1]. | | | |
| The matrix can be bigger, in that case only the submatrix | | | |
| [0..N-1, 0..NCVT-1] will be multiplied by P^T. | | | |
| NCVT - number of columns in matrix VT. | | | |
| | | | |
|
| Output parameters: | | OUTPUT PARAMETERS | |
| D - singular values of matrix B in descending order. | | S - sparse M*N matrix in CRS representation. | |
| U - if NRU>0, contains matrix U*Q. | | You have to fill ALL non-zero elements by calling | |
| VT - if NCVT>0, contains matrix (P^T)*VT. | | SparseSet() BEFORE you try to use this matrix. | |
| C - if NCC>0, contains matrix Q'*C. | | | |
| | | | |
|
| Result: | | NOTE 1. | |
| True, if the algorithm has converged. | | | |
| False, if the algorithm hasn't converged (rare case). | | | |
| | | | |
|
| Additional information: | | Sparse matrices can be stored using either Hash-Table representation or | |
| The type of convergence is controlled by the internal parameter TOL. | | Compressed Row Storage representation. Hast-table is better suited for | |
| If the parameter is greater than 0, the singular values will have | | querying and dynamic operations (thus, it is used for matrix | |
| relative accuracy TOL. If TOL<0, the singular values will have | | initialization), but it is inefficient when you want to make some linear | |
| absolute accuracy ABS(TOL)*norm(B). | | algebra operations. | |
| By default, |TOL| falls within the range of 10*Epsilon and 100*Epsilon, | | | |
| where Epsilon is the machine precision. It is not recommended to use | | | |
| TOL less than 10*Epsilon since this will considerably slow down the | | | |
| algorithm and may not lead to error decreasing. | | | |
| History: | | | |
| * 31 March, 2007. | | | |
| changed MAXITR from 6 to 12. | | | |
| | | | |
|
| -- LAPACK routine (version 3.0) -- | | From the other side, CRS is better suited for linear algebra operations, | |
| Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., | | but initialization is less convenient - you have to tell row sizes at the | |
| Courant Institute, Argonne National Lab, and Rice University | | initialization, and you can fill matrix only row by row, from left to | |
| October 31, 1999. | | right. CRS is also very inefficient when you want to find matrix element | |
| | | by its index. | |
| | | | |
| | | Thus, Hash-Table representation does not support linear algebra | |
| | | operations, while CRS format does not support modification of the table. | |
| | | Tables below outline information about these two formats: | |
| | | | |
| | | OPERATIONS WITH MATRIX HASH CRS | |
| | | create + + | |
| | | read element + + | |
| | | modify element + | |
| | | add value to element + | |
| | | A*x (dense vector) + | |
| | | A'*x (dense vector) + | |
| | | A*X (dense matrix) + | |
| | | A'*X (dense matrix) + | |
| | | | |
| | | NOTE 2. | |
| | | | |
| | | Hash-tables use memory inefficiently, and they have to keep some amount | |
| | | of the "spare memory" in order to have good performance. Hash table for | |
| | | matrix with K non-zero elements will need C*K*(8+2*sizeof(int)) bytes, | |
| | | where C is a small constant, about 1.5-2 in magnitude. | |
| | | | |
| | | CRS storage, from the other side, is more memory-efficient, and needs | |
| | | just K*(8+sizeof(int))+M*sizeof(int) bytes, where M is a number of rows | |
| | | in a matrix. | |
| | | | |
| | | When you convert from the Hash-Table to CRS representation, all unneeded | |
| | | memory will be freed. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| bool rmatrixbdsvd(real_1d_array &d, const real_1d_array &e, const ae_int_t
n, const bool isupper, const bool isfractionalaccuracyrequired, real_2d_arr
ay &u, const ae_int_t nru, real_2d_array &c, const ae_int_t ncc, real_2d_ar
ray &vt, const ae_int_t ncvt); | | void sparsecreatecrs(const ae_int_t m, const ae_int_t n, const integer_1d_a
rray &ner, sparsematrix &s); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Singular value decomposition of a rectangular matrix. | | This function copies S0 to S1. | |
| | | | |
|
| The algorithm calculates the singular value decomposition of a matrix of | | NOTE: this function does not verify its arguments, it just copies all | |
| size MxN: A = U * S * V^T | | fields of the structure. | |
| | | | |
|
| The algorithm finds the singular values and, optionally, matrices U and V^T | | -- ALGLIB PROJECT -- | |
| . | | Copyright 14.10.2011 by Bochkanov Sergey | |
| The algorithm can find both first min(M,N) columns of matrix U and rows of | | *************************************************************************/ | |
| matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM | | void sparsecopy(const sparsematrix &s0, sparsematrix &s1); | |
| and NxN respectively). | | | |
| | | | |
|
| Take into account that the subroutine does not return matrix V but V^T. | | /************************************************************************* | |
| | | This function adds value to S[i,j] - element of the sparse matrix. Matrix | |
| | | must be in a Hash-Table mode. | |
| | | | |
|
| Input parameters: | | In case S[i,j] already exists in the table, V i added to its value. In | |
| A - matrix to be decomposed. | | case S[i,j] is non-existent, it is inserted in the table. Table | |
| Array whose indexes range within [0..M-1, 0..N-1]. | | automatically grows when necessary. | |
| M - number of rows in matrix A. | | | |
| N - number of columns in matrix A. | | | |
| UNeeded - 0, 1 or 2. See the description of the parameter U. | | | |
| VTNeeded - 0, 1 or 2. See the description of the parameter VT. | | | |
| AdditionalMemory - | | | |
| If the parameter: | | | |
| * equals 0, the algorithm doesn | | | |
| memory (lower requirements, lower performance). | | | |
| * equals 1, the algorithm uses additional | | | |
| memory of size min(M,N)*min(M,N) of real numbers. | | | |
| It often speeds up the algorithm. | | | |
| * equals 2, the algorithm uses additional | | | |
| memory of size M*min(M,N) of real numbers. | | | |
| It allows to get a maximum performance. | | | |
| The recommended value of the parameter is 2. | | | |
| | | | |
|
| Output parameters: | | INPUT PARAMETERS | |
| W - contains singular values in descending order. | | S - sparse M*N matrix in Hash-Table representation. | |
| U - if UNeeded=0, U isn't changed, the left singular vector | | Exception will be thrown for CRS matrix. | |
| s | | I - row index of the element to modify, 0<=I<M | |
| are not calculated. | | J - column index of the element to modify, 0<=J<N | |
| if Uneeded=1, U contains left singular vectors (first | | V - value to add, must be finite number | |
| min(M,N) columns of matrix U). Array whose indexes rang | | | |
| e | | OUTPUT PARAMETERS | |
| within [0..M-1, 0..Min(M,N)-1]. | | S - modified matrix | |
| if UNeeded=2, U contains matrix U wholly. Array whose | | | |
| indexes range within [0..M-1, 0..M-1]. | | NOTE 1: when S[i,j] is exactly zero after modification, it is deleted | |
| VT - if VTNeeded=0, VT isn | | from the table. | |
| are not calculated. | | | |
| if VTNeeded=1, VT contains right singular vectors (firs | | -- ALGLIB PROJECT -- | |
| t | | Copyright 14.10.2011 by Bochkanov Sergey | |
| min(M,N) rows of matrix V^T). Array whose indexes range | | *************************************************************************/ | |
| within [0..min(M,N)-1, 0..N-1]. | | void sparseadd(const sparsematrix &s, const ae_int_t i, const ae_int_t j, c | |
| if VTNeeded=2, VT contains matrix V^T wholly. Array who | | onst double v); | |
| se | | | |
| indexes range within [0..N-1, 0..N-1]. | | /************************************************************************* | |
| | | This function modifies S[i,j] - element of the sparse matrix. Matrix must | |
| | | be in a Hash-Table mode. | |
| | | | |
| | | In case new value of S[i,j] is zero, this element is deleted from the | |
| | | table. | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse M*N matrix in Hash-Table representation. | |
| | | Exception will be thrown for CRS matrix. | |
| | | I - row index of the element to modify, 0<=I<M | |
| | | J - column index of the element to modify, 0<=J<N | |
| | | V - value to set, must be finite number, can be zero | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | S - modified matrix | |
| | | | |
| | | NOTE: this function has no effect when called with zero V for non- | |
| | | existent element. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparseset(const sparsematrix &s, const ae_int_t i, const ae_int_t j, c | |
| | | onst double v); | |
| | | | |
| | | /************************************************************************* | |
| | | This function returns S[i,j] - element of the sparse matrix. Matrix can | |
| | | be in any mode (Hash-Table or CRS), but this function is less efficient | |
| | | for CRS matrices. Hash-Table matrices can find element in O(1) time, | |
| | | while CRS matrices need O(RS) time, where RS is an number of non-zero | |
| | | elements in a row. | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse M*N matrix in Hash-Table representation. | |
| | | Exception will be thrown for CRS matrix. | |
| | | I - row index of the element to modify, 0<=I<M | |
| | | J - column index of the element to modify, 0<=J<N | |
| | | | |
| | | RESULT | |
| | | value of S[I,J] or zero (in case no element with such index is found) | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | double sparseget(const sparsematrix &s, const ae_int_t i, const ae_int_t j) | |
| | | ; | |
| | | | |
| | | /************************************************************************* | |
| | | This function converts matrix to CRS format. | |
| | | | |
| | | Some algorithms (linear algebra ones, for example) require matrices in | |
| | | CRS format. | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse M*N matrix. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | S - matrix in CRS format | |
| | | | |
| | | NOTE: this function has no effect when called with matrix which is | |
| | | already in CRS mode. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparseconverttocrs(const sparsematrix &s); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates matrix-vector product S*x. Matrix S must be | |
| | | stored in CRS format (exception will be thrown otherwise). | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse M*N matrix in CRS format (you MUST convert it | |
| | | to CRS before calling this function). | |
| | | X - array[N], input vector. For performance reasons we | |
| | | make only quick checks - we check that array size is | |
| | | at least N, but we do not check for NAN's or INF's. | |
| | | Y - output buffer, possibly preallocated. In case buffer | |
| | | size is too small to store result, this buffer is | |
| | | automatically resized. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | Y - array[M], S*x | |
| | | | |
| | | NOTE: this function throws exception when called for non-CRS matrix. You | |
| | | must convert your matrix with SparseConvertToCRS() before using this | |
| | | function. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparsemv(const sparsematrix &s, const real_1d_array &x, real_1d_array | |
| | | &y); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates matrix-vector product S^T*x. Matrix S must be | |
| | | stored in CRS format (exception will be thrown otherwise). | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse M*N matrix in CRS format (you MUST convert it | |
| | | to CRS before calling this function). | |
| | | X - array[M], input vector. For performance reasons we | |
| | | make only quick checks - we check that array size is | |
| | | at least M, but we do not check for NAN's or INF's. | |
| | | Y - output buffer, possibly preallocated. In case buffer | |
| | | size is too small to store result, this buffer is | |
| | | automatically resized. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | Y - array[N], S^T*x | |
| | | | |
| | | NOTE: this function throws exception when called for non-CRS matrix. You | |
| | | must convert your matrix with SparseConvertToCRS() before using this | |
| | | function. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparsemtv(const sparsematrix &s, const real_1d_array &x, real_1d_array | |
| | | &y); | |
| | | | |
| | | /************************************************************************* | |
| | | This function simultaneously calculates two matrix-vector products: | |
| | | S*x and S^T*x. | |
| | | S must be square (non-rectangular) matrix stored in CRS format (exception | |
| | | will be thrown otherwise). | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse N*N matrix in CRS format (you MUST convert it | |
| | | to CRS before calling this function). | |
| | | X - array[N], input vector. For performance reasons we | |
| | | make only quick checks - we check that array size is | |
| | | at least N, but we do not check for NAN's or INF's. | |
| | | Y0 - output buffer, possibly preallocated. In case buffer | |
| | | size is too small to store result, this buffer is | |
| | | automatically resized. | |
| | | Y1 - output buffer, possibly preallocated. In case buffer | |
| | | size is too small to store result, this buffer is | |
| | | automatically resized. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | Y0 - array[N], S*x | |
| | | Y1 - array[N], S^T*x | |
| | | | |
| | | NOTE: this function throws exception when called for non-CRS matrix. You | |
| | | must convert your matrix with SparseConvertToCRS() before using this | |
| | | function. It also throws exception when S is non-square. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparsemv2(const sparsematrix &s, const real_1d_array &x, real_1d_array | |
| | | &y0, real_1d_array &y1); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates matrix-vector product S*x, when S is symmetric | |
| | | matrix. Matrix S must be stored in CRS format (exception will be | |
| | | thrown otherwise). | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse M*M matrix in CRS format (you MUST convert it | |
| | | to CRS before calling this function). | |
| | | X - array[N], input vector. For performance reasons we | |
| | | make only quick checks - we check that array size is | |
| | | at least N, but we do not check for NAN's or INF's. | |
| | | Y - output buffer, possibly preallocated. In case buffer | |
| | | size is too small to store result, this buffer is | |
| | | automatically resized. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | Y - array[M], S*x | |
| | | | |
| | | NOTE: this function throws exception when called for non-CRS matrix. You | |
| | | must convert your matrix with SparseConvertToCRS() before using this | |
| | | function. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparsesmv(const sparsematrix &s, const bool isupper, const real_1d_arr | |
| | | ay &x, real_1d_array &y); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates matrix-matrix product S*A. Matrix S must be | |
| | | stored in CRS format (exception will be thrown otherwise). | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse M*N matrix in CRS format (you MUST convert it | |
| | | to CRS before calling this function). | |
| | | A - array[N][K], input dense matrix. For performance reaso | |
| | | ns | |
| | | we make only quick checks - we check that array size | |
| | | is at least N, but we do not check for NAN's or INF's. | |
| | | K - number of columns of matrix (A). | |
| | | B - output buffer, possibly preallocated. In case buffer | |
| | | size is too small to store result, this buffer is | |
| | | automatically resized. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | B - array[M][K], S*A | |
| | | | |
| | | NOTE: this function throws exception when called for non-CRS matrix. You | |
| | | must convert your matrix with SparseConvertToCRS() before using this | |
| | | function. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparsemm(const sparsematrix &s, const real_2d_array &a, const ae_int_t | |
| | | k, real_2d_array &b); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates matrix-matrix product S^T*A. Matrix S must be | |
| | | stored in CRS format (exception will be thrown otherwise). | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse M*N matrix in CRS format (you MUST convert it | |
| | | to CRS before calling this function). | |
| | | A - array[M][K], input dense matrix. For performance reason | |
| | | s | |
| | | we make only quick checks - we check that array size i | |
| | | s | |
| | | at least M, but we do not check for NAN's or INF's. | |
| | | K - number of columns of matrix (A). | |
| | | B - output buffer, possibly preallocated. In case buffer | |
| | | size is too small to store result, this buffer is | |
| | | automatically resized. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | B - array[N][K], S^T*A | |
| | | | |
| | | NOTE: this function throws exception when called for non-CRS matrix. You | |
| | | must convert your matrix with SparseConvertToCRS() before using this | |
| | | function. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparsemtm(const sparsematrix &s, const real_2d_array &a, const ae_int_ | |
| | | t k, real_2d_array &b); | |
| | | | |
| | | /************************************************************************* | |
| | | This function simultaneously calculates two matrix-matrix products: | |
| | | S*A and S^T*A. | |
| | | S must be square (non-rectangular) matrix stored in CRS format (exception | |
| | | will be thrown otherwise). | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse N*N matrix in CRS format (you MUST convert it | |
| | | to CRS before calling this function). | |
| | | A - array[N][K], input dense matrix. For performance reason | |
| | | s | |
| | | we make only quick checks - we check that array size i | |
| | | s | |
| | | at least N, but we do not check for NAN's or INF's. | |
| | | K - number of columns of matrix (A). | |
| | | B0 - output buffer, possibly preallocated. In case buffer | |
| | | size is too small to store result, this buffer is | |
| | | automatically resized. | |
| | | B1 - output buffer, possibly preallocated. In case buffer | |
| | | size is too small to store result, this buffer is | |
| | | automatically resized. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | B0 - array[N][K], S*A | |
| | | B1 - array[N][K], S^T*A | |
| | | | |
| | | NOTE: this function throws exception when called for non-CRS matrix. You | |
| | | must convert your matrix with SparseConvertToCRS() before using this | |
| | | function. It also throws exception when S is non-square. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparsemm2(const sparsematrix &s, const real_2d_array &a, const ae_int_ | |
| | | t k, real_2d_array &b0, real_2d_array &b1); | |
| | | | |
| | | /************************************************************************* | |
| | | This function calculates matrix-matrix product S*A, when S is symmetric | |
| | | matrix. Matrix S must be stored in CRS format (exception will be | |
| | | thrown otherwise). | |
| | | | |
| | | INPUT PARAMETERS | |
| | | S - sparse M*M matrix in CRS format (you MUST convert it | |
| | | to CRS before calling this function). | |
| | | A - array[N][K], input dense matrix. For performance reason | |
| | | s | |
| | | we make only quick checks - we check that array size is | |
| | | at least N, but we do not check for NAN's or INF's. | |
| | | K - number of columns of matrix (A). | |
| | | B - output buffer, possibly preallocated. In case buffer | |
| | | size is too small to store result, this buffer is | |
| | | automatically resized. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | B - array[M][K], S*A | |
| | | | |
| | | NOTE: this function throws exception when called for non-CRS matrix. You | |
| | | must convert your matrix with SparseConvertToCRS() before using this | |
| | | function. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparsesmm(const sparsematrix &s, const bool isupper, const real_2d_arr | |
| | | ay &a, const ae_int_t k, real_2d_array &b); | |
| | | | |
| | | /************************************************************************* | |
| | | This procedure resizes Hash-Table matrix. It can be called when you have | |
| | | deleted too many elements from the matrix, and you want to free unneeded | |
| | | memory. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 14.10.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void sparseresizematrix(const sparsematrix &s); | |
| | | | |
| | | /************************************************************************* | |
| | | This procedure initializes matrix norm estimator. | |
| | | | |
| | | USAGE: | |
| | | 1. User initializes algorithm state with NormEstimatorCreate() call | |
| | | 2. User calls NormEstimatorEstimateSparse() (or NormEstimatorIteration()) | |
| | | 3. User calls NormEstimatorResults() to get solution. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | M - number of rows in the matrix being estimated, M>0 | |
| | | N - number of columns in the matrix being estimated, N>0 | |
| | | NStart - number of random starting vectors | |
| | | recommended value - at least 5. | |
| | | NIts - number of iterations to do with best starting vector | |
| | | recommended value - at least 5. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | NOTE: this algorithm is effectively deterministic, i.e. it always returns | |
| | | same result when repeatedly called for the same matrix. In fact, algorithm | |
| | | uses randomized starting vectors, but internal random numbers generator | |
| | | always generates same sequence of the random values (it is a feature, not | |
| | | bug). | |
| | | | |
| | | Algorithm can be made non-deterministic with NormEstimatorSetSeed(0) call. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
|
| Copyright 2005 by Bochkanov Sergey | | Copyright 06.12.2011 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| bool rmatrixsvd(const real_2d_array &a, const ae_int_t m, const ae_int_t n, | | void normestimatorcreate(const ae_int_t m, const ae_int_t n, const ae_int_t | |
| const ae_int_t uneeded, const ae_int_t vtneeded, const ae_int_t additional | | nstart, const ae_int_t nits, normestimatorstate &state); | |
| memory, real_1d_array &w, real_2d_array &u, real_2d_array &vt); | | | |
| | | /************************************************************************* | |
| | | This function changes seed value used by algorithm. In some cases we need | |
| | | deterministic processing, i.e. subsequent calls must return equal results, | |
| | | in other cases we need non-deterministic algorithm which returns different | |
| | | results for the same matrix on every pass. | |
| | | | |
| | | Setting zero seed will lead to non-deterministic algorithm, while non-zero | |
| | | value will make our algorithm deterministic. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | State - norm estimator state, must be initialized with a call | |
| | | to NormEstimatorCreate() | |
| | | SeedVal - seed value, >=0. Zero value = non-deterministic algo. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 06.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void normestimatorsetseed(const normestimatorstate &state, const ae_int_t s | |
| | | eedval); | |
| | | | |
| | | /************************************************************************* | |
| | | This function estimates norm of the sparse M*N matrix A. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | State - norm estimator state, must be initialized with a call | |
| | | to NormEstimatorCreate() | |
| | | A - sparse M*N matrix, must be converted to CRS format | |
| | | prior to calling this function. | |
| | | | |
| | | After this function is over you can call NormEstimatorResults() to get | |
| | | estimate of the norm(A). | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 06.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void normestimatorestimatesparse(const normestimatorstate &state, const spa | |
| | | rsematrix &a); | |
| | | | |
| | | /************************************************************************* | |
| | | Matrix norm estimation results | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | State - algorithm state | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Nrm - estimate of the matrix norm, Nrm>=0 | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 06.12.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void normestimatorresults(const normestimatorstate &state, double &nrm); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| Determinant calculation of the matrix given by its LU decomposition. | | Determinant calculation of the matrix given by its LU decomposition. | |
| | | | |
| Input parameters: | | Input parameters: | |
| A - LU decomposition of the matrix (output of | | A - LU decomposition of the matrix (output of | |
| RMatrixLU subroutine). | | RMatrixLU subroutine). | |
| Pivots - table of permutations which were made during | | Pivots - table of permutations which were made during | |
| the LU decomposition. | | the LU decomposition. | |
| Output of RMatrixLU subroutine. | | Output of RMatrixLU subroutine. | |
| | | | |
| skipping to change at line 3465 | | skipping to change at line 4221 | |
|---|
| ae_int_t n, | | ae_int_t n, | |
| /* Complex */ ae_vector* tau, | | /* Complex */ ae_vector* tau, | |
| ae_int_t qrows, | | ae_int_t qrows, | |
| /* Complex */ ae_matrix* q, | | /* Complex */ ae_matrix* q, | |
| ae_state *_state); | | ae_state *_state); | |
| void cmatrixlqunpackl(/* Complex */ ae_matrix* a, | | void cmatrixlqunpackl(/* Complex */ ae_matrix* a, | |
| ae_int_t m, | | ae_int_t m, | |
| ae_int_t n, | | ae_int_t n, | |
| /* Complex */ ae_matrix* l, | | /* Complex */ ae_matrix* l, | |
| ae_state *_state); | | ae_state *_state); | |
|
| | | void rmatrixqrbasecase(/* Real */ ae_matrix* a, | |
| | | ae_int_t m, | |
| | | ae_int_t n, | |
| | | /* Real */ ae_vector* work, | |
| | | /* Real */ ae_vector* t, | |
| | | /* Real */ ae_vector* tau, | |
| | | ae_state *_state); | |
| | | void rmatrixlqbasecase(/* Real */ ae_matrix* a, | |
| | | ae_int_t m, | |
| | | ae_int_t n, | |
| | | /* Real */ ae_vector* work, | |
| | | /* Real */ ae_vector* t, | |
| | | /* Real */ ae_vector* tau, | |
| | | ae_state *_state); | |
| void rmatrixbd(/* Real */ ae_matrix* a, | | void rmatrixbd(/* Real */ ae_matrix* a, | |
| ae_int_t m, | | ae_int_t m, | |
| ae_int_t n, | | ae_int_t n, | |
| /* Real */ ae_vector* tauq, | | /* Real */ ae_vector* tauq, | |
| /* Real */ ae_vector* taup, | | /* Real */ ae_vector* taup, | |
| ae_state *_state); | | ae_state *_state); | |
| void rmatrixbdunpackq(/* Real */ ae_matrix* qp, | | void rmatrixbdunpackq(/* Real */ ae_matrix* qp, | |
| ae_int_t m, | | ae_int_t m, | |
| ae_int_t n, | | ae_int_t n, | |
| /* Real */ ae_vector* tauq, | | /* Real */ ae_vector* tauq, | |
| | | | |
| skipping to change at line 3551 | | skipping to change at line 4321 | |
|---|
| /* Complex */ ae_vector* tau, | | /* Complex */ ae_vector* tau, | |
| /* Real */ ae_vector* d, | | /* Real */ ae_vector* d, | |
| /* Real */ ae_vector* e, | | /* Real */ ae_vector* e, | |
| ae_state *_state); | | ae_state *_state); | |
| void hmatrixtdunpackq(/* Complex */ ae_matrix* a, | | void hmatrixtdunpackq(/* Complex */ ae_matrix* a, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_bool isupper, | | ae_bool isupper, | |
| /* Complex */ ae_vector* tau, | | /* Complex */ ae_vector* tau, | |
| /* Complex */ ae_matrix* q, | | /* Complex */ ae_matrix* q, | |
| ae_state *_state); | | ae_state *_state); | |
|
| | | ae_bool rmatrixbdsvd(/* Real */ ae_vector* d, | |
| | | /* Real */ ae_vector* e, | |
| | | ae_int_t n, | |
| | | ae_bool isupper, | |
| | | ae_bool isfractionalaccuracyrequired, | |
| | | /* Real */ ae_matrix* u, | |
| | | ae_int_t nru, | |
| | | /* Real */ ae_matrix* c, | |
| | | ae_int_t ncc, | |
| | | /* Real */ ae_matrix* vt, | |
| | | ae_int_t ncvt, | |
| | | ae_state *_state); | |
| | | ae_bool bidiagonalsvddecomposition(/* Real */ ae_vector* d, | |
| | | /* Real */ ae_vector* e, | |
| | | ae_int_t n, | |
| | | ae_bool isupper, | |
| | | ae_bool isfractionalaccuracyrequired, | |
| | | /* Real */ ae_matrix* u, | |
| | | ae_int_t nru, | |
| | | /* Real */ ae_matrix* c, | |
| | | ae_int_t ncc, | |
| | | /* Real */ ae_matrix* vt, | |
| | | ae_int_t ncvt, | |
| | | ae_state *_state); | |
| | | ae_bool rmatrixsvd(/* Real */ ae_matrix* a, | |
| | | ae_int_t m, | |
| | | ae_int_t n, | |
| | | ae_int_t uneeded, | |
| | | ae_int_t vtneeded, | |
| | | ae_int_t additionalmemory, | |
| | | /* Real */ ae_vector* w, | |
| | | /* Real */ ae_matrix* u, | |
| | | /* Real */ ae_matrix* vt, | |
| | | ae_state *_state); | |
| ae_bool smatrixevd(/* Real */ ae_matrix* a, | | ae_bool smatrixevd(/* Real */ ae_matrix* a, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_int_t zneeded, | | ae_int_t zneeded, | |
| ae_bool isupper, | | ae_bool isupper, | |
| /* Real */ ae_vector* d, | | /* Real */ ae_vector* d, | |
| /* Real */ ae_matrix* z, | | /* Real */ ae_matrix* z, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_bool smatrixevdr(/* Real */ ae_matrix* a, | | ae_bool smatrixevdr(/* Real */ ae_matrix* a, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_int_t zneeded, | | ae_int_t zneeded, | |
| | | | |
| skipping to change at line 3854 | | skipping to change at line 4658 | |
|---|
| void cmatrixtrinverse(/* Complex */ ae_matrix* a, | | void cmatrixtrinverse(/* Complex */ ae_matrix* a, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_bool isupper, | | ae_bool isupper, | |
| ae_bool isunit, | | ae_bool isunit, | |
| ae_int_t* info, | | ae_int_t* info, | |
| matinvreport* rep, | | matinvreport* rep, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_bool _matinvreport_init(matinvreport* p, ae_state *_state, ae_bool make_
automatic); | | ae_bool _matinvreport_init(matinvreport* p, ae_state *_state, ae_bool make_
automatic); | |
| ae_bool _matinvreport_init_copy(matinvreport* dst, matinvreport* src, ae_st
ate *_state, ae_bool make_automatic); | | ae_bool _matinvreport_init_copy(matinvreport* dst, matinvreport* src, ae_st
ate *_state, ae_bool make_automatic); | |
| void _matinvreport_clear(matinvreport* p); | | void _matinvreport_clear(matinvreport* p); | |
|
| ae_bool rmatrixbdsvd(/* Real */ ae_vector* d, | | | |
| /* Real */ ae_vector* e, | | | |
| ae_int_t n, | | | |
| ae_bool isupper, | | | |
| ae_bool isfractionalaccuracyrequired, | | | |
| /* Real */ ae_matrix* u, | | | |
| ae_int_t nru, | | | |
| /* Real */ ae_matrix* c, | | | |
| ae_int_t ncc, | | | |
| /* Real */ ae_matrix* vt, | | | |
| ae_int_t ncvt, | | | |
| ae_state *_state); | | | |
| ae_bool bidiagonalsvddecomposition(/* Real */ ae_vector* d, | | | |
| /* Real */ ae_vector* e, | | | |
| ae_int_t n, | | | |
| ae_bool isupper, | | | |
| ae_bool isfractionalaccuracyrequired, | | | |
| /* Real */ ae_matrix* u, | | | |
| ae_int_t nru, | | | |
| /* Real */ ae_matrix* c, | | | |
| ae_int_t ncc, | | | |
| /* Real */ ae_matrix* vt, | | | |
| ae_int_t ncvt, | | | |
| ae_state *_state); | | | |
| ae_bool rmatrixsvd(/* Real */ ae_matrix* a, | | | |
| ae_int_t m, | | | |
| ae_int_t n, | | | |
| ae_int_t uneeded, | | | |
| ae_int_t vtneeded, | | | |
| ae_int_t additionalmemory, | | | |
| /* Real */ ae_vector* w, | | | |
| /* Real */ ae_matrix* u, | | | |
| /* Real */ ae_matrix* vt, | | | |
| ae_state *_state); | | | |
| void fblscholeskysolve(/* Real */ ae_matrix* cha, | | void fblscholeskysolve(/* Real */ ae_matrix* cha, | |
| double sqrtscalea, | | double sqrtscalea, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_bool isupper, | | ae_bool isupper, | |
| /* Real */ ae_vector* xb, | | /* Real */ ae_vector* xb, | |
| /* Real */ ae_vector* tmp, | | /* Real */ ae_vector* tmp, | |
| ae_state *_state); | | ae_state *_state); | |
| void fblssolvecgx(/* Real */ ae_matrix* a, | | void fblssolvecgx(/* Real */ ae_matrix* a, | |
| ae_int_t m, | | ae_int_t m, | |
| ae_int_t n, | | ae_int_t n, | |
| | | | |
| skipping to change at line 3909 | | skipping to change at line 4679 | |
|---|
| /* Real */ ae_vector* b, | | /* Real */ ae_vector* b, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| /* Real */ ae_vector* buf, | | /* Real */ ae_vector* buf, | |
| ae_state *_state); | | ae_state *_state); | |
| void fblscgcreate(/* Real */ ae_vector* x, | | void fblscgcreate(/* Real */ ae_vector* x, | |
| /* Real */ ae_vector* b, | | /* Real */ ae_vector* b, | |
| ae_int_t n, | | ae_int_t n, | |
| fblslincgstate* state, | | fblslincgstate* state, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_bool fblscgiteration(fblslincgstate* state, ae_state *_state); | | ae_bool fblscgiteration(fblslincgstate* state, ae_state *_state); | |
|
| | | void fblssolvels(/* Real */ ae_matrix* a, | |
| | | /* Real */ ae_vector* b, | |
| | | ae_int_t m, | |
| | | ae_int_t n, | |
| | | /* Real */ ae_vector* tmp0, | |
| | | /* Real */ ae_vector* tmp1, | |
| | | /* Real */ ae_vector* tmp2, | |
| | | ae_state *_state); | |
| ae_bool _fblslincgstate_init(fblslincgstate* p, ae_state *_state, ae_bool m
ake_automatic); | | ae_bool _fblslincgstate_init(fblslincgstate* p, ae_state *_state, ae_bool m
ake_automatic); | |
| ae_bool _fblslincgstate_init_copy(fblslincgstate* dst, fblslincgstate* src,
ae_state *_state, ae_bool make_automatic); | | ae_bool _fblslincgstate_init_copy(fblslincgstate* dst, fblslincgstate* src,
ae_state *_state, ae_bool make_automatic); | |
| void _fblslincgstate_clear(fblslincgstate* p); | | void _fblslincgstate_clear(fblslincgstate* p); | |
|
| | | void sparsecreate(ae_int_t m, | |
| | | ae_int_t n, | |
| | | ae_int_t k, | |
| | | sparsematrix* s, | |
| | | ae_state *_state); | |
| | | void sparsecreatecrs(ae_int_t m, | |
| | | ae_int_t n, | |
| | | /* Integer */ ae_vector* ner, | |
| | | sparsematrix* s, | |
| | | ae_state *_state); | |
| | | void sparsecopy(sparsematrix* s0, sparsematrix* s1, ae_state *_state); | |
| | | void sparseadd(sparsematrix* s, | |
| | | ae_int_t i, | |
| | | ae_int_t j, | |
| | | double v, | |
| | | ae_state *_state); | |
| | | void sparseset(sparsematrix* s, | |
| | | ae_int_t i, | |
| | | ae_int_t j, | |
| | | double v, | |
| | | ae_state *_state); | |
| | | double sparseget(sparsematrix* s, | |
| | | ae_int_t i, | |
| | | ae_int_t j, | |
| | | ae_state *_state); | |
| | | void sparseconverttocrs(sparsematrix* s, ae_state *_state); | |
| | | void sparsemv(sparsematrix* s, | |
| | | /* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_state *_state); | |
| | | void sparsemtv(sparsematrix* s, | |
| | | /* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_state *_state); | |
| | | void sparsemv2(sparsematrix* s, | |
| | | /* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y0, | |
| | | /* Real */ ae_vector* y1, | |
| | | ae_state *_state); | |
| | | void sparsesmv(sparsematrix* s, | |
| | | ae_bool isupper, | |
| | | /* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_state *_state); | |
| | | void sparsemm(sparsematrix* s, | |
| | | /* Real */ ae_matrix* a, | |
| | | ae_int_t k, | |
| | | /* Real */ ae_matrix* b, | |
| | | ae_state *_state); | |
| | | void sparsemtm(sparsematrix* s, | |
| | | /* Real */ ae_matrix* a, | |
| | | ae_int_t k, | |
| | | /* Real */ ae_matrix* b, | |
| | | ae_state *_state); | |
| | | void sparsemm2(sparsematrix* s, | |
| | | /* Real */ ae_matrix* a, | |
| | | ae_int_t k, | |
| | | /* Real */ ae_matrix* b0, | |
| | | /* Real */ ae_matrix* b1, | |
| | | ae_state *_state); | |
| | | void sparsesmm(sparsematrix* s, | |
| | | ae_bool isupper, | |
| | | /* Real */ ae_matrix* a, | |
| | | ae_int_t k, | |
| | | /* Real */ ae_matrix* b, | |
| | | ae_state *_state); | |
| | | void sparseresizematrix(sparsematrix* s, ae_state *_state); | |
| | | double sparsegetaveragelengthofchain(sparsematrix* s, ae_state *_state); | |
| | | ae_bool _sparsematrix_init(sparsematrix* p, ae_state *_state, ae_bool make_ | |
| | | automatic); | |
| | | ae_bool _sparsematrix_init_copy(sparsematrix* dst, sparsematrix* src, ae_st | |
| | | ate *_state, ae_bool make_automatic); | |
| | | void _sparsematrix_clear(sparsematrix* p); | |
| | | void normestimatorcreate(ae_int_t m, | |
| | | ae_int_t n, | |
| | | ae_int_t nstart, | |
| | | ae_int_t nits, | |
| | | normestimatorstate* state, | |
| | | ae_state *_state); | |
| | | void normestimatorsetseed(normestimatorstate* state, | |
| | | ae_int_t seedval, | |
| | | ae_state *_state); | |
| | | ae_bool normestimatoriteration(normestimatorstate* state, | |
| | | ae_state *_state); | |
| | | void normestimatorestimatesparse(normestimatorstate* state, | |
| | | sparsematrix* a, | |
| | | ae_state *_state); | |
| | | void normestimatorresults(normestimatorstate* state, | |
| | | double* nrm, | |
| | | ae_state *_state); | |
| | | void normestimatorrestart(normestimatorstate* state, ae_state *_state); | |
| | | ae_bool _normestimatorstate_init(normestimatorstate* p, ae_state *_state, a | |
| | | e_bool make_automatic); | |
| | | ae_bool _normestimatorstate_init_copy(normestimatorstate* dst, normestimato | |
| | | rstate* src, ae_state *_state, ae_bool make_automatic); | |
| | | void _normestimatorstate_clear(normestimatorstate* p); | |
| double rmatrixludet(/* Real */ ae_matrix* a, | | double rmatrixludet(/* Real */ ae_matrix* a, | |
| /* Integer */ ae_vector* pivots, | | /* Integer */ ae_vector* pivots, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_state *_state); | | ae_state *_state); | |
| double rmatrixdet(/* Real */ ae_matrix* a, | | double rmatrixdet(/* Real */ ae_matrix* a, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_complex cmatrixludet(/* Complex */ ae_matrix* a, | | ae_complex cmatrixludet(/* Complex */ ae_matrix* a, | |
| /* Integer */ ae_vector* pivots, | | /* Integer */ ae_vector* pivots, | |
| ae_int_t n, | | ae_int_t n, | |
| | | | |
| End of changes. 43 change blocks. |
|---|
| 169 lines changed or deleted | | 1068 lines changed or added | |
|---|
|
| solvers.h | | solvers.h | |
|---|
| | | | |
| skipping to change at line 47 | | skipping to change at line 47 | |
|---|
| } densesolverreport; | | } densesolverreport; | |
| typedef struct | | typedef struct | |
| { | | { | |
| double r2; | | double r2; | |
| ae_matrix cx; | | ae_matrix cx; | |
| ae_int_t n; | | ae_int_t n; | |
| ae_int_t k; | | ae_int_t k; | |
| } densesolverlsreport; | | } densesolverlsreport; | |
| typedef struct | | typedef struct | |
| { | | { | |
|
| | | normestimatorstate nes; | |
| | | ae_vector rx; | |
| | | ae_vector b; | |
| | | ae_int_t n; | |
| | | ae_int_t m; | |
| | | ae_vector ui; | |
| | | ae_vector uip1; | |
| | | ae_vector vi; | |
| | | ae_vector vip1; | |
| | | ae_vector omegai; | |
| | | ae_vector omegaip1; | |
| | | double alphai; | |
| | | double alphaip1; | |
| | | double betai; | |
| | | double betaip1; | |
| | | double phibari; | |
| | | double phibarip1; | |
| | | double phii; | |
| | | double rhobari; | |
| | | double rhobarip1; | |
| | | double rhoi; | |
| | | double ci; | |
| | | double si; | |
| | | double theta; | |
| | | double lambdai; | |
| | | ae_vector d; | |
| | | double anorm; | |
| | | double bnorm2; | |
| | | double dnorm; | |
| | | double r2; | |
| | | ae_vector x; | |
| | | ae_vector mv; | |
| | | ae_vector mtv; | |
| | | double epsa; | |
| | | double epsb; | |
| | | double epsc; | |
| | | ae_int_t maxits; | |
| | | ae_bool xrep; | |
| | | ae_bool xupdated; | |
| | | ae_bool needmv; | |
| | | ae_bool needmtv; | |
| | | ae_bool needmv2; | |
| | | ae_bool needvmv; | |
| | | ae_bool needprec; | |
| | | ae_int_t repiterationscount; | |
| | | ae_int_t repnmv; | |
| | | ae_int_t repterminationtype; | |
| | | ae_bool running; | |
| | | rcommstate rstate; | |
| | | } linlsqrstate; | |
| | | typedef struct | |
| | | { | |
| | | ae_int_t iterationscount; | |
| | | ae_int_t nmv; | |
| | | ae_int_t terminationtype; | |
| | | } linlsqrreport; | |
| | | typedef struct | |
| | | { | |
| | | ae_vector rx; | |
| | | ae_vector b; | |
| | | ae_int_t n; | |
| | | ae_vector cx; | |
| | | ae_vector cr; | |
| | | ae_vector cz; | |
| | | ae_vector p; | |
| | | ae_vector r; | |
| | | ae_vector z; | |
| | | double alpha; | |
| | | double beta; | |
| | | double r2; | |
| | | double meritfunction; | |
| | | ae_vector x; | |
| | | ae_vector mv; | |
| | | ae_vector pv; | |
| | | double vmv; | |
| | | ae_vector startx; | |
| | | double epsf; | |
| | | ae_int_t maxits; | |
| | | ae_int_t itsbeforerestart; | |
| | | ae_int_t itsbeforerupdate; | |
| | | ae_bool xrep; | |
| | | ae_bool xupdated; | |
| | | ae_bool needmv; | |
| | | ae_bool needmtv; | |
| | | ae_bool needmv2; | |
| | | ae_bool needvmv; | |
| | | ae_bool needprec; | |
| | | ae_int_t repiterationscount; | |
| | | ae_int_t repnmv; | |
| | | ae_int_t repterminationtype; | |
| | | ae_bool running; | |
| | | rcommstate rstate; | |
| | | } lincgstate; | |
| | | typedef struct | |
| | | { | |
| | | ae_int_t iterationscount; | |
| | | ae_int_t nmv; | |
| | | ae_int_t terminationtype; | |
| | | double r2; | |
| | | } lincgreport; | |
| | | typedef struct | |
| | | { | |
| ae_int_t n; | | ae_int_t n; | |
| ae_int_t m; | | ae_int_t m; | |
| double epsf; | | double epsf; | |
| ae_int_t maxits; | | ae_int_t maxits; | |
| ae_bool xrep; | | ae_bool xrep; | |
| double stpmax; | | double stpmax; | |
| ae_vector x; | | ae_vector x; | |
| double f; | | double f; | |
| ae_vector fi; | | ae_vector fi; | |
| ae_matrix j; | | ae_matrix j; | |
| | | | |
| skipping to change at line 147 | | skipping to change at line 249 | |
|---|
| densesolverlsreport& operator=(const densesolverlsreport &rhs); | | densesolverlsreport& operator=(const densesolverlsreport &rhs); | |
| virtual ~densesolverlsreport(); | | virtual ~densesolverlsreport(); | |
| double &r2; | | double &r2; | |
| real_2d_array cx; | | real_2d_array cx; | |
| ae_int_t &n; | | ae_int_t &n; | |
| ae_int_t &k; | | ae_int_t &k; | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | This object stores state of the LinLSQR method. | |
| | | | |
| | | You should use ALGLIB functions to work with this object. | |
| | | *************************************************************************/ | |
| | | class _linlsqrstate_owner | |
| | | { | |
| | | public: | |
| | | _linlsqrstate_owner(); | |
| | | _linlsqrstate_owner(const _linlsqrstate_owner &rhs); | |
| | | _linlsqrstate_owner& operator=(const _linlsqrstate_owner &rhs); | |
| | | virtual ~_linlsqrstate_owner(); | |
| | | alglib_impl::linlsqrstate* c_ptr(); | |
| | | alglib_impl::linlsqrstate* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::linlsqrstate *p_struct; | |
| | | }; | |
| | | class linlsqrstate : public _linlsqrstate_owner | |
| | | { | |
| | | public: | |
| | | linlsqrstate(); | |
| | | linlsqrstate(const linlsqrstate &rhs); | |
| | | linlsqrstate& operator=(const linlsqrstate &rhs); | |
| | | virtual ~linlsqrstate(); | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| | | | |
| | | *************************************************************************/ | |
| | | class _linlsqrreport_owner | |
| | | { | |
| | | public: | |
| | | _linlsqrreport_owner(); | |
| | | _linlsqrreport_owner(const _linlsqrreport_owner &rhs); | |
| | | _linlsqrreport_owner& operator=(const _linlsqrreport_owner &rhs); | |
| | | virtual ~_linlsqrreport_owner(); | |
| | | alglib_impl::linlsqrreport* c_ptr(); | |
| | | alglib_impl::linlsqrreport* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::linlsqrreport *p_struct; | |
| | | }; | |
| | | class linlsqrreport : public _linlsqrreport_owner | |
| | | { | |
| | | public: | |
| | | linlsqrreport(); | |
| | | linlsqrreport(const linlsqrreport &rhs); | |
| | | linlsqrreport& operator=(const linlsqrreport &rhs); | |
| | | virtual ~linlsqrreport(); | |
| | | ae_int_t &iterationscount; | |
| | | ae_int_t &nmv; | |
| | | ae_int_t &terminationtype; | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| | | This object stores state of the linear CG method. | |
| | | | |
| | | You should use ALGLIB functions to work with this object. | |
| | | Never try to access its fields directly! | |
| | | *************************************************************************/ | |
| | | class _lincgstate_owner | |
| | | { | |
| | | public: | |
| | | _lincgstate_owner(); | |
| | | _lincgstate_owner(const _lincgstate_owner &rhs); | |
| | | _lincgstate_owner& operator=(const _lincgstate_owner &rhs); | |
| | | virtual ~_lincgstate_owner(); | |
| | | alglib_impl::lincgstate* c_ptr(); | |
| | | alglib_impl::lincgstate* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::lincgstate *p_struct; | |
| | | }; | |
| | | class lincgstate : public _lincgstate_owner | |
| | | { | |
| | | public: | |
| | | lincgstate(); | |
| | | lincgstate(const lincgstate &rhs); | |
| | | lincgstate& operator=(const lincgstate &rhs); | |
| | | virtual ~lincgstate(); | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| | | | |
| | | *************************************************************************/ | |
| | | class _lincgreport_owner | |
| | | { | |
| | | public: | |
| | | _lincgreport_owner(); | |
| | | _lincgreport_owner(const _lincgreport_owner &rhs); | |
| | | _lincgreport_owner& operator=(const _lincgreport_owner &rhs); | |
| | | virtual ~_lincgreport_owner(); | |
| | | alglib_impl::lincgreport* c_ptr(); | |
| | | alglib_impl::lincgreport* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::lincgreport *p_struct; | |
| | | }; | |
| | | class lincgreport : public _lincgreport_owner | |
| | | { | |
| | | public: | |
| | | lincgreport(); | |
| | | lincgreport(const lincgreport &rhs); | |
| | | lincgreport& operator=(const lincgreport &rhs); | |
| | | virtual ~lincgreport(); | |
| | | ae_int_t &iterationscount; | |
| | | ae_int_t &nmv; | |
| | | ae_int_t &terminationtype; | |
| | | double &r2; | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| | | | |
| *************************************************************************/ | | *************************************************************************/ | |
| class _nleqstate_owner | | class _nleqstate_owner | |
| { | | { | |
| public: | | public: | |
| _nleqstate_owner(); | | _nleqstate_owner(); | |
| _nleqstate_owner(const _nleqstate_owner &rhs); | | _nleqstate_owner(const _nleqstate_owner &rhs); | |
| _nleqstate_owner& operator=(const _nleqstate_owner &rhs); | | _nleqstate_owner& operator=(const _nleqstate_owner &rhs); | |
| virtual ~_nleqstate_owner(); | | virtual ~_nleqstate_owner(); | |
| alglib_impl::nleqstate* c_ptr(); | | alglib_impl::nleqstate* c_ptr(); | |
| | | | |
| skipping to change at line 878 | | skipping to change at line 1092 | |
|---|
| * K dim(Null(A)) | | * K dim(Null(A)) | |
| * CX array[0..N-1,0..K-1], kernel of A. | | * CX array[0..N-1,0..K-1], kernel of A. | |
| Columns of CX store such vectors that A*CX[i]=0. | | Columns of CX store such vectors that A*CX[i]=0. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 24.08.2009 by Bochkanov Sergey | | 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); | | 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); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | This function initializes linear LSQR Solver. This solver is used to solve | |
| | | non-symmetric (and, possibly, non-square) problems. Least squares solution | |
| | | is returned for non-compatible systems. | |
| | | | |
| | | USAGE: | |
| | | 1. User initializes algorithm state with LinLSQRCreate() call | |
| | | 2. User tunes solver parameters with LinLSQRSetCond() and other functions | |
| | | 3. User calls LinLSQRSolveSparse() function which takes algorithm state | |
| | | and SparseMatrix object. | |
| | | 4. User calls LinLSQRResults() to get solution | |
| | | 5. Optionally, user may call LinLSQRSolveSparse() again to solve another | |
| | | problem with different matrix and/or right part without reinitializing | |
| | | LinLSQRState structure. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | M - number of rows in A | |
| | | N - number of variables, N>0 | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 30.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void linlsqrcreate(const ae_int_t m, const ae_int_t n, linlsqrstate &state) | |
| | | ; | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets optional Tikhonov regularization coefficient. | |
| | | It is zero by default. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | LambdaI - regularization factor, LambdaI>=0 | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 30.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void linlsqrsetlambdai(const linlsqrstate &state, const double lambdai); | |
| | | | |
| | | /************************************************************************* | |
| | | Procedure for solution of A*x=b with sparse A. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | State - algorithm state | |
| | | A - sparse M*N matrix in the CRS format (you MUST contvert it | |
| | | to CRS format by calling SparseConvertToCRS() function | |
| | | BEFORE you pass it to this function). | |
| | | B - right part, array[M] | |
| | | | |
| | | RESULT: | |
| | | This function returns no result. | |
| | | You can get solution by calling LinCGResults() | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 30.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void linlsqrsolvesparse(const linlsqrstate &state, const sparsematrix &a, c | |
| | | onst real_1d_array &b); | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets stopping criteria. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | EpsA - algorithm will be stopped if ||A^T*Rk||/(||A||*||Rk||)<=Eps | |
| | | A. | |
| | | EpsB - algorithm will be stopped if ||Rk||<=EpsB*||B|| | |
| | | MaxIts - algorithm will be stopped if number of iterations | |
| | | more than MaxIts. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | NOTE: if EpsA,EpsB,EpsC and MaxIts are zero then these variables will | |
| | | be setted as default values. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 30.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void linlsqrsetcond(const linlsqrstate &state, const double epsa, const dou | |
| | | ble epsb, const ae_int_t maxits); | |
| | | | |
| | | /************************************************************************* | |
| | | LSQR solver: results. | |
| | | | |
| | | This function must be called after LinLSQRSolve | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | State - algorithm state | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | X - array[N], solution | |
| | | Rep - optimization report: | |
| | | * Rep.TerminationType completetion code: | |
| | | * 1 ||Rk||<=EpsB*||B|| | |
| | | * 4 ||A^T*Rk||/(||A||*||Rk||)<=EpsA | |
| | | * 5 MaxIts steps was taken | |
| | | * 7 rounding errors prevent further progress, | |
| | | X contains best point found so far. | |
| | | (sometimes returned on singular systems) | |
| | | * Rep.IterationsCount contains iterations count | |
| | | * NMV countains number of matrix-vector calculations | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 30.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void linlsqrresults(const linlsqrstate &state, real_1d_array &x, linlsqrrep | |
| | | ort &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | 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 MinCGOptimize(). | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 30.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void linlsqrsetxrep(const linlsqrstate &state, const bool needxrep); | |
| | | | |
| | | /************************************************************************* | |
| | | This function initializes linear CG Solver. This solver is used to solve | |
| | | symmetric positive definite problems. If you want to solve nonsymmetric | |
| | | (or non-positive definite) problem you may use LinLSQR solver provided by | |
| | | ALGLIB. | |
| | | | |
| | | USAGE: | |
| | | 1. User initializes algorithm state with LinCGCreate() call | |
| | | 2. User tunes solver parameters with LinCGSetCond() and other functions | |
| | | 3. Optionally, user sets starting point with LinCGSetStartingPoint() | |
| | | 4. User calls LinCGSolveSparse() function which takes algorithm state and | |
| | | SparseMatrix object. | |
| | | 5. User calls LinCGResults() to get solution | |
| | | 6. Optionally, user may call LinCGSolveSparse() again to solve another | |
| | | problem with different matrix and/or right part without reinitializing | |
| | | LinCGState structure. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | N - problem dimension, N>0 | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 14.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lincgcreate(const ae_int_t n, lincgstate &state); | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets starting point. | |
| | | By default, zero starting point is used. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - starting point, array[N] | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 14.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lincgsetstartingpoint(const lincgstate &state, const real_1d_array &x) | |
| | | ; | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets stopping criteria. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | EpsF - algorithm will be stopped if norm of residual is less than | |
| | | EpsF*||b||. | |
| | | MaxIts - algorithm will be stopped if number of iterations is more | |
| | | than MaxIts. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | NOTES: | |
| | | If both EpsF and MaxIts are zero then small EpsF will be set to small | |
| | | value. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 14.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lincgsetcond(const lincgstate &state, const double epsf, const ae_int_ | |
| | | t maxits); | |
| | | | |
| | | /************************************************************************* | |
| | | Procedure for solution of A*x=b with sparse A. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | State - algorithm state | |
| | | A - sparse matrix in the CRS format (you MUST contvert it to | |
| | | CRS format by calling SparseConvertToCRS() function). | |
| | | IsUpper - whether upper or lower triangle of A is used: | |
| | | * IsUpper=True => only upper triangle is used and lower | |
| | | triangle is not referenced at all | |
| | | * IsUpper=False => only lower triangle is used and upper | |
| | | triangle is not referenced at all | |
| | | B - right part, array[N] | |
| | | | |
| | | RESULT: | |
| | | This function returns no result. | |
| | | You can get solution by calling LinCGResults() | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 14.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lincgsolvesparse(const lincgstate &state, const sparsematrix &a, const | |
| | | bool isupper, const real_1d_array &b); | |
| | | | |
| | | /************************************************************************* | |
| | | CG-solver: results. | |
| | | | |
| | | This function must be called after LinCGSolve | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | State - algorithm state | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | X - array[N], solution | |
| | | Rep - optimization report: | |
| | | * Rep.TerminationType completetion code: | |
| | | * -5 input matrix is either not positive definite, | |
| | | too large or too small | |
| | | * -4 overflow/underflow during solution | |
| | | (ill conditioned problem) | |
| | | * 1 ||residual||<=EpsF*||b|| | |
| | | * 5 MaxIts steps was taken | |
| | | * 7 rounding errors prevent further progress, | |
| | | best point found is returned | |
| | | * Rep.IterationsCount contains iterations count | |
| | | * NMV countains number of matrix-vector calculations | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 14.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lincgresults(const lincgstate &state, real_1d_array &x, lincgreport &r | |
| | | ep); | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets restart frequency. By default, algorithm is restarted | |
| | | after N subsequent iterations. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 14.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lincgsetrestartfreq(const lincgstate &state, const ae_int_t srf); | |
| | | | |
| | | /************************************************************************* | |
| | | This function sets frequency of residual recalculations. | |
| | | | |
| | | Algorithm updates residual r_k using iterative formula, but recalculates | |
| | | it from scratch after each 10 iterations. It is done to avoid accumulation | |
| | | of numerical errors and to stop algorithm when r_k starts to grow. | |
| | | | |
| | | Such low update frequence (1/10) gives very little overhead, but makes | |
| | | algorithm a bit more robust against numerical errors. However, you may | |
| | | change it | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | Freq - desired update frequency, Freq>=0. | |
| | | Zero value means that no updates will be done. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 14.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lincgsetrupdatefreq(const lincgstate &state, const ae_int_t freq); | |
| | | | |
| | | /************************************************************************* | |
| | | 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 MinCGOptimize(). | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 14.11.2011 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lincgsetxrep(const lincgstate &state, const bool needxrep); | |
| | | | |
| | | /************************************************************************* | |
| LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER | | LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER | |
| | | | |
| DESCRIPTION: | | DESCRIPTION: | |
| This algorithm solves system of nonlinear equations | | This algorithm solves system of nonlinear equations | |
| F[0](x[0], ..., x[n-1]) = 0 | | F[0](x[0], ..., x[n-1]) = 0 | |
| F[1](x[0], ..., x[n-1]) = 0 | | F[1](x[0], ..., x[n-1]) = 0 | |
| ... | | ... | |
| F[M-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 | | with M/N do not necessarily coincide. Algorithm converges quadratically | |
| under following conditions: | | under following conditions: | |
| | | | |
| skipping to change at line 1284 | | skipping to change at line 1779 | |
|---|
| ae_int_t* info, | | ae_int_t* info, | |
| densesolverlsreport* rep, | | densesolverlsreport* rep, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_bool _densesolverreport_init(densesolverreport* p, ae_state *_state, ae_
bool make_automatic); | | ae_bool _densesolverreport_init(densesolverreport* p, ae_state *_state, ae_
bool make_automatic); | |
| ae_bool _densesolverreport_init_copy(densesolverreport* dst, densesolverrep
ort* src, ae_state *_state, ae_bool make_automatic); | | ae_bool _densesolverreport_init_copy(densesolverreport* dst, densesolverrep
ort* src, ae_state *_state, ae_bool make_automatic); | |
| void _densesolverreport_clear(densesolverreport* p); | | void _densesolverreport_clear(densesolverreport* p); | |
| ae_bool _densesolverlsreport_init(densesolverlsreport* p, ae_state *_state,
ae_bool make_automatic); | | ae_bool _densesolverlsreport_init(densesolverlsreport* p, ae_state *_state,
ae_bool make_automatic); | |
| ae_bool _densesolverlsreport_init_copy(densesolverlsreport* dst, densesolve
rlsreport* src, ae_state *_state, ae_bool make_automatic); | | ae_bool _densesolverlsreport_init_copy(densesolverlsreport* dst, densesolve
rlsreport* src, ae_state *_state, ae_bool make_automatic); | |
| void _densesolverlsreport_clear(densesolverlsreport* p); | | void _densesolverlsreport_clear(densesolverlsreport* p); | |
|
| | | void linlsqrcreate(ae_int_t m, | |
| | | ae_int_t n, | |
| | | linlsqrstate* state, | |
| | | ae_state *_state); | |
| | | void linlsqrsetb(linlsqrstate* state, | |
| | | /* Real */ ae_vector* b, | |
| | | ae_state *_state); | |
| | | void linlsqrsetlambdai(linlsqrstate* state, | |
| | | double lambdai, | |
| | | ae_state *_state); | |
| | | ae_bool linlsqriteration(linlsqrstate* state, ae_state *_state); | |
| | | void linlsqrsolvesparse(linlsqrstate* state, | |
| | | sparsematrix* a, | |
| | | /* Real */ ae_vector* b, | |
| | | ae_state *_state); | |
| | | void linlsqrsetcond(linlsqrstate* state, | |
| | | double epsa, | |
| | | double epsb, | |
| | | ae_int_t maxits, | |
| | | ae_state *_state); | |
| | | void linlsqrresults(linlsqrstate* state, | |
| | | /* Real */ ae_vector* x, | |
| | | linlsqrreport* rep, | |
| | | ae_state *_state); | |
| | | void linlsqrsetxrep(linlsqrstate* state, | |
| | | ae_bool needxrep, | |
| | | ae_state *_state); | |
| | | void linlsqrrestart(linlsqrstate* state, ae_state *_state); | |
| | | ae_bool _linlsqrstate_init(linlsqrstate* p, ae_state *_state, ae_bool make_ | |
| | | automatic); | |
| | | ae_bool _linlsqrstate_init_copy(linlsqrstate* dst, linlsqrstate* src, ae_st | |
| | | ate *_state, ae_bool make_automatic); | |
| | | void _linlsqrstate_clear(linlsqrstate* p); | |
| | | ae_bool _linlsqrreport_init(linlsqrreport* p, ae_state *_state, ae_bool mak | |
| | | e_automatic); | |
| | | ae_bool _linlsqrreport_init_copy(linlsqrreport* dst, linlsqrreport* src, ae | |
| | | _state *_state, ae_bool make_automatic); | |
| | | void _linlsqrreport_clear(linlsqrreport* p); | |
| | | void lincgcreate(ae_int_t n, lincgstate* state, ae_state *_state); | |
| | | void lincgsetstartingpoint(lincgstate* state, | |
| | | /* Real */ ae_vector* x, | |
| | | ae_state *_state); | |
| | | void lincgsetb(lincgstate* state, | |
| | | /* Real */ ae_vector* b, | |
| | | ae_state *_state); | |
| | | void lincgsetcond(lincgstate* state, | |
| | | double epsf, | |
| | | ae_int_t maxits, | |
| | | ae_state *_state); | |
| | | ae_bool lincgiteration(lincgstate* state, ae_state *_state); | |
| | | void lincgsolvesparse(lincgstate* state, | |
| | | sparsematrix* a, | |
| | | ae_bool isupper, | |
| | | /* Real */ ae_vector* b, | |
| | | ae_state *_state); | |
| | | void lincgresults(lincgstate* state, | |
| | | /* Real */ ae_vector* x, | |
| | | lincgreport* rep, | |
| | | ae_state *_state); | |
| | | void lincgsetrestartfreq(lincgstate* state, | |
| | | ae_int_t srf, | |
| | | ae_state *_state); | |
| | | void lincgsetrupdatefreq(lincgstate* state, | |
| | | ae_int_t freq, | |
| | | ae_state *_state); | |
| | | void lincgsetxrep(lincgstate* state, ae_bool needxrep, ae_state *_state); | |
| | | void lincgrestart(lincgstate* state, ae_state *_state); | |
| | | ae_bool _lincgstate_init(lincgstate* p, ae_state *_state, ae_bool make_auto | |
| | | matic); | |
| | | ae_bool _lincgstate_init_copy(lincgstate* dst, lincgstate* src, ae_state *_ | |
| | | state, ae_bool make_automatic); | |
| | | void _lincgstate_clear(lincgstate* p); | |
| | | ae_bool _lincgreport_init(lincgreport* p, ae_state *_state, ae_bool make_au | |
| | | tomatic); | |
| | | ae_bool _lincgreport_init_copy(lincgreport* dst, lincgreport* src, ae_state | |
| | | *_state, ae_bool make_automatic); | |
| | | void _lincgreport_clear(lincgreport* p); | |
| void nleqcreatelm(ae_int_t n, | | void nleqcreatelm(ae_int_t n, | |
| ae_int_t m, | | ae_int_t m, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| nleqstate* state, | | nleqstate* state, | |
| ae_state *_state); | | ae_state *_state); | |
| void nleqsetcond(nleqstate* state, | | void nleqsetcond(nleqstate* state, | |
| double epsf, | | double epsf, | |
| ae_int_t maxits, | | ae_int_t maxits, | |
| ae_state *_state); | | ae_state *_state); | |
| void nleqsetxrep(nleqstate* state, ae_bool needxrep, ae_state *_state); | | void nleqsetxrep(nleqstate* state, ae_bool needxrep, ae_state *_state); | |
| | | | |
| End of changes. 4 change blocks. |
|---|
| 0 lines changed or deleted | | 581 lines changed or added | |
|---|
|