From aa3a7c33030857b781f0bd6392932290a16f3dba Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Sat, 22 Aug 2009 06:44:41 +0200 Subject: [PATCH] raw import of covar() : this is the last one, and we now do not depend on the cminpack library anymore. --- unsupported/Eigen/NonLinear | 7 -- unsupported/Eigen/src/NonLinear/covar.h | 129 ++++++++++++++++++++++++ unsupported/test/CMakeLists.txt | 4 - unsupported/test/NonLinear.cpp | 4 +- 4 files changed, 131 insertions(+), 13 deletions(-) create mode 100644 unsupported/Eigen/src/NonLinear/covar.h diff --git a/unsupported/Eigen/NonLinear b/unsupported/Eigen/NonLinear index 79c7e4c4c..e904844aa 100644 --- a/unsupported/Eigen/NonLinear +++ b/unsupported/Eigen/NonLinear @@ -36,10 +36,6 @@ namespace Eigen { */ //@{ -#if 1 -#include -#else - /* Declarations for minpack */ 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 ); @@ -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_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/r1updt.h" #include "src/NonLinear/r1mpyq.h" diff --git a/unsupported/Eigen/src/NonLinear/covar.h b/unsupported/Eigen/src/NonLinear/covar.h new file mode 100644 index 000000000..69e649d63 --- /dev/null +++ b/unsupported/Eigen/src/NonLinear/covar.h @@ -0,0 +1,129 @@ + +template +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_ */ + diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index 89b91ca29..007a9c4c2 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -15,10 +15,6 @@ else(ADOLC_FOUND) ei_add_property(EIGEN_MISSING_BACKENDS "Adolc") 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(autodiff) ei_add_test(BVH) diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp index 6092ad59b..299ee9019 100644 --- a/unsupported/test/NonLinear.cpp +++ b/unsupported/test/NonLinear.cpp @@ -247,7 +247,7 @@ void testLmder() covfac = fnorm*fnorm/(m-n); VectorXd wa(n); ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides) - covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data()); + ei_covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data()); MatrixXd cov_ref(n,n); cov_ref << @@ -763,7 +763,7 @@ void testLmdif() covfac = fnorm*fnorm/(m-n); VectorXd wa(n); ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides) - covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data()); + ei_covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data()); MatrixXd cov_ref(n,n); cov_ref <<