From 113995355b1063a8d2971f50df0df66c3721097d Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Tue, 26 Jan 2010 08:42:48 +0100 Subject: get rid of ei_qform + lot of other cleaning, now that we do not depend on minpack qr factorization. --- .../NonLinearOptimization/HybridNonLinearSolver.h | 51 ++++-------------- .../Eigen/src/NonLinearOptimization/qform.h | 60 ---------------------- 2 files changed, 10 insertions(+), 101 deletions(-) delete mode 100644 unsupported/Eigen/src/NonLinearOptimization/qform.h (limited to 'unsupported/Eigen/src') diff --git a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h index ee5b63c44..4b985db74 100644 --- a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h +++ b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h @@ -217,7 +217,7 @@ HybridNonLinearSolver::solveOneStep( const int mode ) { - int i, j; + int j; jeval = true; /* calculate the jacobian matrix. */ @@ -246,30 +246,15 @@ HybridNonLinearSolver::solveOneStep( /* compute the qr factorization of the jacobian. */ wa2 = fjac.colwise().blueNorm(); HouseholderQR qrfac(fjac); // no pivoting: - fjac = qrfac.matrixQR(); - wa1 = fjac.diagonal(); - fjac.diagonal() = qrfac.hCoeffs(); - // TODO : avoid this: - for(int ii=0; ii< fjac.cols(); ii++) fjac.col(ii).segment(ii+1, fjac.rows()-ii-1) *= fjac(ii,ii); // rescale vectors - - /* form (q transpose)*fvec and store in qtf. */ - qtf = fvec; - for (j = 0; j < n; ++j) - if (fjac(j,j) != 0.) { - sum = 0.; - for (i = j; i < n; ++i) - sum += fjac(i,j) * qtf[i]; - temp = -sum / fjac(j,j); - for (i = j; i < n; ++i) - qtf[i] += fjac(i,j) * temp; - } /* copy the triangular factor of the qr factorization into r. */ R = qrfac.matrixQR(); - sing = wa1.cwiseAbs().minCoeff()==0.; /* accumulate the orthogonal factor in fjac. */ - ei_qform(n, n, fjac.data(), fjac.rows(), wa1.data()); + fjac = qrfac.householderQ(); + + /* form (q transpose)*fvec and store in qtf. */ + qtf = fjac.transpose() * fvec; /* rescale if necessary. */ if (mode != 2) @@ -480,13 +465,12 @@ HybridNonLinearSolver::solveNumericalDiffOneStep( const int mode ) { - int i, j; + int j; jeval = true; if (parameters.nb_of_subdiagonals<0) parameters.nb_of_subdiagonals= n-1; if (parameters.nb_of_superdiagonals<0) parameters.nb_of_superdiagonals= n-1; /* calculate the jacobian matrix. */ - if (ei_fdjac1(functor, x, fvec, fjac, parameters.nb_of_subdiagonals, parameters.nb_of_superdiagonals, parameters.epsfcn) <0) return UserAksed; nfev += std::min(parameters.nb_of_subdiagonals+parameters.nb_of_superdiagonals+ 1, n); @@ -512,30 +496,15 @@ HybridNonLinearSolver::solveNumericalDiffOneStep( /* compute the qr factorization of the jacobian. */ wa2 = fjac.colwise().blueNorm(); HouseholderQR qrfac(fjac); // no pivoting: - fjac = qrfac.matrixQR(); - wa1 = fjac.diagonal(); - fjac.diagonal() = qrfac.hCoeffs(); - // TODO : avoid this: - for(int ii=0; ii< fjac.cols(); ii++) fjac.col(ii).segment(ii+1, fjac.rows()-ii-1) *= fjac(ii,ii); // rescale vectors - - /* form (q transpose)*fvec and store in qtf. */ - qtf = fvec; - for (j = 0; j < n; ++j) - if (fjac(j,j) != 0.) { - sum = 0.; - for (i = j; i < n; ++i) - sum += fjac(i,j) * qtf[i]; - temp = -sum / fjac(j,j); - for (i = j; i < n; ++i) - qtf[i] += fjac(i,j) * temp; - } /* copy the triangular factor of the qr factorization into r. */ R = qrfac.matrixQR(); - sing = wa1.cwiseAbs().minCoeff()==0.; /* accumulate the orthogonal factor in fjac. */ - ei_qform(n, n, fjac.data(), fjac.rows(), wa1.data()); + fjac = qrfac.householderQ(); + + /* form (q transpose)*fvec and store in qtf. */ + qtf = fjac.transpose() * fvec; /* rescale if necessary. */ if (mode != 2) diff --git a/unsupported/Eigen/src/NonLinearOptimization/qform.h b/unsupported/Eigen/src/NonLinearOptimization/qform.h deleted file mode 100644 index 1c3e3e267..000000000 --- a/unsupported/Eigen/src/NonLinearOptimization/qform.h +++ /dev/null @@ -1,60 +0,0 @@ - -template -void ei_qform(int m, int n, Scalar *q, int - ldq, Scalar *wa) -{ - /* System generated locals */ - int q_dim1, q_offset; - - /* Local variables */ - int i, j, k, l; - Scalar sum, temp; - int minmn; - - /* Parameter adjustments */ - --wa; - q_dim1 = ldq; - q_offset = 1 + q_dim1 * 1; - q -= q_offset; - - /* Function Body */ - - /* zero out upper triangle of q in the first min(m,n) columns. */ - - minmn = std::min(m,n); - for (j = 2; j <= minmn; ++j) { - for (i = 1; i <= j-1; ++i) - q[i + j * q_dim1] = 0.; - } - - /* initialize remaining columns to those of the identity matrix. */ - - for (j = n+1; j <= m; ++j) { - for (i = 1; i <= m; ++i) - q[i + j * q_dim1] = 0.; - q[j + j * q_dim1] = 1.; - } - - /* accumulate q from its factored form. */ - - for (l = 1; l <= minmn; ++l) { - k = minmn - l + 1; - for (i = k; i <= m; ++i) { - wa[i] = q[i + k * q_dim1]; - q[i + k * q_dim1] = 0.; - } - q[k + k * q_dim1] = 1.; - if (wa[k] == 0.) - continue; - for (j = k; j <= m; ++j) { - sum = 0.; - for (i = k; i <= m; ++i) - sum += q[i + j * q_dim1] * wa[i]; - temp = sum / wa[k]; - for (i = k; i <= m; ++i) - q[i + j * q_dim1] -= temp * wa[i]; - } - } - -} - -- cgit v1.2.3