mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-08-14 20:56:00 +08:00
raw import of covar() : this is the last one, and we now do not depend on
the cminpack library anymore.
This commit is contained in:
parent
16061a46db
commit
aa3a7c3303
@ -36,10 +36,6 @@ namespace Eigen {
|
|||||||
*/
|
*/
|
||||||
//@{
|
//@{
|
||||||
|
|
||||||
#if 1
|
|
||||||
#include <cminpack.h>
|
|
||||||
#else
|
|
||||||
|
|
||||||
/* Declarations for minpack */
|
/* Declarations for minpack */
|
||||||
typedef int (*minpack_func_nn)(void *p, int n, const double *x, double *fvec, int iflag );
|
typedef int (*minpack_func_nn)(void *p, int n, const double *x, double *fvec, int iflag );
|
||||||
typedef int (*minpack_funcder_nn)(void *p, int n, const double *x, double *fvec, double *fjac, int ldfjac, int iflag );
|
typedef int (*minpack_funcder_nn)(void *p, int n, const double *x, double *fvec, double *fjac, int ldfjac, int iflag );
|
||||||
@ -47,9 +43,6 @@ typedef int (*minpack_func_mn)(void *p, int m, int n, const double *x, double *f
|
|||||||
typedef int (*minpack_funcder_mn)(void *p, int m, int n, const double *x, double *fvec, double *fjac, int ldfjac, int iflag );
|
typedef int (*minpack_funcder_mn)(void *p, int m, int n, const double *x, double *fvec, double *fjac, int ldfjac, int iflag );
|
||||||
typedef int (*minpack_funcderstr_mn)(void *p, int m, int n, const double *x, double *fvec, double *fjrow, int iflag );
|
typedef int (*minpack_funcderstr_mn)(void *p, int m, int n, const double *x, double *fvec, double *fjrow, int iflag );
|
||||||
|
|
||||||
void covar(int n, double *r__, int ldr, const int *ipvt, double tol, double *wa);
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#include "src/NonLinear/qrsolv.h"
|
#include "src/NonLinear/qrsolv.h"
|
||||||
#include "src/NonLinear/r1updt.h"
|
#include "src/NonLinear/r1updt.h"
|
||||||
#include "src/NonLinear/r1mpyq.h"
|
#include "src/NonLinear/r1mpyq.h"
|
||||||
|
129
unsupported/Eigen/src/NonLinear/covar.h
Normal file
129
unsupported/Eigen/src/NonLinear/covar.h
Normal file
@ -0,0 +1,129 @@
|
|||||||
|
|
||||||
|
template <typename Scalar>
|
||||||
|
void ei_covar(int n, Scalar *r__, int ldr,
|
||||||
|
const int *ipvt, Scalar tol, Scalar *wa)
|
||||||
|
{
|
||||||
|
/* System generated locals */
|
||||||
|
int r_dim1, r_offset, i__1, i__2, i__3;
|
||||||
|
|
||||||
|
/* Local variables */
|
||||||
|
int i__, j, k, l, ii, jj, km1;
|
||||||
|
int sing;
|
||||||
|
Scalar temp, tolr;
|
||||||
|
|
||||||
|
/* Parameter adjustments */
|
||||||
|
--wa;
|
||||||
|
--ipvt;
|
||||||
|
tolr = tol * ei_abs(r__[0]);
|
||||||
|
r_dim1 = ldr;
|
||||||
|
r_offset = 1 + r_dim1;
|
||||||
|
r__ -= r_offset;
|
||||||
|
|
||||||
|
/* Function Body */
|
||||||
|
|
||||||
|
/* form the inverse of r in the full upper triangle of r. */
|
||||||
|
|
||||||
|
l = 0;
|
||||||
|
i__1 = n;
|
||||||
|
for (k = 1; k <= i__1; ++k) {
|
||||||
|
if (ei_abs(r__[k + k * r_dim1]) <= tolr) {
|
||||||
|
goto L50;
|
||||||
|
}
|
||||||
|
r__[k + k * r_dim1] = 1. / r__[k + k * r_dim1];
|
||||||
|
km1 = k - 1;
|
||||||
|
if (km1 < 1) {
|
||||||
|
goto L30;
|
||||||
|
}
|
||||||
|
i__2 = km1;
|
||||||
|
for (j = 1; j <= i__2; ++j) {
|
||||||
|
temp = r__[k + k * r_dim1] * r__[j + k * r_dim1];
|
||||||
|
r__[j + k * r_dim1] = 0.;
|
||||||
|
i__3 = j;
|
||||||
|
for (i__ = 1; i__ <= i__3; ++i__) {
|
||||||
|
r__[i__ + k * r_dim1] -= temp * r__[i__ + j * r_dim1];
|
||||||
|
/* L10: */
|
||||||
|
}
|
||||||
|
/* L20: */
|
||||||
|
}
|
||||||
|
L30:
|
||||||
|
l = k;
|
||||||
|
/* L40: */
|
||||||
|
}
|
||||||
|
L50:
|
||||||
|
|
||||||
|
/* form the full upper triangle of the inverse of (r transpose)*r */
|
||||||
|
/* in the full upper triangle of r. */
|
||||||
|
|
||||||
|
if (l < 1) {
|
||||||
|
goto L110;
|
||||||
|
}
|
||||||
|
i__1 = l;
|
||||||
|
for (k = 1; k <= i__1; ++k) {
|
||||||
|
km1 = k - 1;
|
||||||
|
if (km1 < 1) {
|
||||||
|
goto L80;
|
||||||
|
}
|
||||||
|
i__2 = km1;
|
||||||
|
for (j = 1; j <= i__2; ++j) {
|
||||||
|
temp = r__[j + k * r_dim1];
|
||||||
|
i__3 = j;
|
||||||
|
for (i__ = 1; i__ <= i__3; ++i__) {
|
||||||
|
r__[i__ + j * r_dim1] += temp * r__[i__ + k * r_dim1];
|
||||||
|
/* L60: */
|
||||||
|
}
|
||||||
|
/* L70: */
|
||||||
|
}
|
||||||
|
L80:
|
||||||
|
temp = r__[k + k * r_dim1];
|
||||||
|
i__2 = k;
|
||||||
|
for (i__ = 1; i__ <= i__2; ++i__) {
|
||||||
|
r__[i__ + k * r_dim1] = temp * r__[i__ + k * r_dim1];
|
||||||
|
/* L90: */
|
||||||
|
}
|
||||||
|
/* L100: */
|
||||||
|
}
|
||||||
|
L110:
|
||||||
|
|
||||||
|
/* form the full lower triangle of the covariance matrix */
|
||||||
|
/* in the strict lower triangle of r and in wa. */
|
||||||
|
|
||||||
|
i__1 = n;
|
||||||
|
for (j = 1; j <= i__1; ++j) {
|
||||||
|
jj = ipvt[j];
|
||||||
|
sing = j > l;
|
||||||
|
i__2 = j;
|
||||||
|
for (i__ = 1; i__ <= i__2; ++i__) {
|
||||||
|
if (sing) {
|
||||||
|
r__[i__ + j * r_dim1] = 0.;
|
||||||
|
}
|
||||||
|
ii = ipvt[i__];
|
||||||
|
if (ii > jj) {
|
||||||
|
r__[ii + jj * r_dim1] = r__[i__ + j * r_dim1];
|
||||||
|
}
|
||||||
|
if (ii < jj) {
|
||||||
|
r__[jj + ii * r_dim1] = r__[i__ + j * r_dim1];
|
||||||
|
}
|
||||||
|
/* L120: */
|
||||||
|
}
|
||||||
|
wa[jj] = r__[j + j * r_dim1];
|
||||||
|
/* L130: */
|
||||||
|
}
|
||||||
|
|
||||||
|
/* symmetrize the covariance matrix in r. */
|
||||||
|
|
||||||
|
i__1 = n;
|
||||||
|
for (j = 1; j <= i__1; ++j) {
|
||||||
|
i__2 = j;
|
||||||
|
for (i__ = 1; i__ <= i__2; ++i__) {
|
||||||
|
r__[i__ + j * r_dim1] = r__[j + i__ * r_dim1];
|
||||||
|
/* L140: */
|
||||||
|
}
|
||||||
|
r__[j + j * r_dim1] = wa[j];
|
||||||
|
/* L150: */
|
||||||
|
}
|
||||||
|
/*return 0;*/
|
||||||
|
|
||||||
|
/* last card of subroutine covar. */
|
||||||
|
|
||||||
|
} /* covar_ */
|
||||||
|
|
@ -15,10 +15,6 @@ else(ADOLC_FOUND)
|
|||||||
ei_add_property(EIGEN_MISSING_BACKENDS "Adolc")
|
ei_add_property(EIGEN_MISSING_BACKENDS "Adolc")
|
||||||
endif(ADOLC_FOUND)
|
endif(ADOLC_FOUND)
|
||||||
|
|
||||||
# temporary :
|
|
||||||
include_directories(/home/orzel/tmp/cminpack-1.0.2/)
|
|
||||||
LINK_LIBRARIES(/home/orzel/tmp/cminpack-1.0.2/libminpack.a)
|
|
||||||
|
|
||||||
ei_add_test(NonLinear)
|
ei_add_test(NonLinear)
|
||||||
ei_add_test(autodiff)
|
ei_add_test(autodiff)
|
||||||
ei_add_test(BVH)
|
ei_add_test(BVH)
|
||||||
|
@ -247,7 +247,7 @@ void testLmder()
|
|||||||
covfac = fnorm*fnorm/(m-n);
|
covfac = fnorm*fnorm/(m-n);
|
||||||
VectorXd wa(n);
|
VectorXd wa(n);
|
||||||
ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides)
|
ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides)
|
||||||
covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
|
ei_covar<double>(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
|
||||||
|
|
||||||
MatrixXd cov_ref(n,n);
|
MatrixXd cov_ref(n,n);
|
||||||
cov_ref <<
|
cov_ref <<
|
||||||
@ -763,7 +763,7 @@ void testLmdif()
|
|||||||
covfac = fnorm*fnorm/(m-n);
|
covfac = fnorm*fnorm/(m-n);
|
||||||
VectorXd wa(n);
|
VectorXd wa(n);
|
||||||
ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides)
|
ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides)
|
||||||
covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
|
ei_covar<double>(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
|
||||||
|
|
||||||
MatrixXd cov_ref(n,n);
|
MatrixXd cov_ref(n,n);
|
||||||
cov_ref <<
|
cov_ref <<
|
||||||
|
Loading…
x
Reference in New Issue
Block a user