mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

get rid of ei_qform + lot of other cleaning, now that we do not depend on

Browse Source 
minpack qr factorization.
This commit is contained in:
Thomas Capricelli 2010-01-26 08:42:48 +01:00
parent ba2a9cce03
commit 113995355b
2 changed files with 10 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
View File

@ -217,7 +217,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
const int mode
)
{
int i, j;
int j;
jeval = true;
/* calculate the jacobian matrix. */
@ -246,30 +246,15 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
/* compute the qr factorization of the jacobian. */
wa2 = fjac.colwise().blueNorm();
HouseholderQR<JacobianType> 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<Scalar>(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<FunctorType,Scalar>::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<FunctorType,Scalar>::solveNumericalDiffOneStep(
/* compute the qr factorization of the jacobian. */
wa2 = fjac.colwise().blueNorm();
HouseholderQR<JacobianType> 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<Scalar>(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)

60
unsupported/Eigen/src/NonLinearOptimization/qform.h
View File

@ -1,60 +0,0 @@
template <typename Scalar>
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];
}
}
}