computes column norms outside of ei_qrfac()

This commit is contained in:
Thomas Capricelli 2009-11-26 02:53:58 +01:00
parent 9cbfdbad22
commit 746c787a76
3 changed files with 27 additions and 36 deletions

View File

@ -226,15 +226,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
return UserAksed;
++njev;
/* compute the qr factorization of the jacobian. */
ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, wa1.data(), wa2.data());
wa2 = fjac.colwise().blueNorm();
/* on the first iteration and if mode is 1, scale according */
/* to the norms of the columns of the initial jacobian. */
if (iter == 1) {
/* on the first iteration and if mode is 1, scale according */
/* to the norms of the columns of the initial jacobian. */
if (mode != 2)
for (j = 0; j < n; ++j) {
diag[j] = wa2[j];
@ -251,6 +247,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
delta = parameters.factor;
}
/* compute the qr factorization of the jacobian. */
ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, wa1.data());
/* form (q transpose)*fvec and store in qtf. */
qtf = fvec;
@ -269,18 +268,16 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
sing = false;
for (j = 0; j < n; ++j) {
l = j;
if (j)
for (i = 0; i < j; ++i) {
R[l] = fjac(i,j);
l = l + n - i -1;
}
for (i = 0; i < j; ++i) {
R[l] = fjac(i,j);
l = l + n - i -1;
}
R[l] = wa1[j];
if (wa1[j] == 0.)
sing = true;
}
/* accumulate the orthogonal factor in fjac. */
ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
/* rescale if necessary. */
@ -543,13 +540,10 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
return UserAksed;
nfev += std::min(parameters.nb_of_subdiagonals+parameters.nb_of_superdiagonals+ 1, n);
/* compute the qr factorization of the jacobian. */
ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, wa1.data(), wa2.data());
wa2 = fjac.colwise().blueNorm();
/* on the first iteration and if mode is 1, scale according */
/* to the norms of the columns of the initial jacobian. */
if (iter == 1) {
if (mode != 2)
for (j = 0; j < n; ++j) {
@ -560,7 +554,6 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
/* on the first iteration, calculate the norm of the scaled x */
/* and initialize the step bound delta. */
wa3 = diag.cwise() * x;
xnorm = wa3.stableNorm();
delta = parameters.factor * xnorm;
@ -568,6 +561,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
delta = parameters.factor;
}
/* compute the qr factorization of the jacobian. */
ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, wa1.data());
/* form (q transpose)*fvec and store in qtf. */
qtf = fvec;
@ -586,18 +582,16 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
sing = false;
for (j = 0; j < n; ++j) {
l = j;
if (j)
for (i = 0; i < j; ++i) {
R[l] = fjac(i,j);
l = l + n - i -1;
}
for (i = 0; i < j; ++i) {
R[l] = fjac(i,j);
l = l + n - i -1;
}
R[l] = wa1[j];
if (wa1[j] == 0.)
sing = true;
}
/* accumulate the orthogonal factor in fjac. */
ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
/* rescale if necessary. */

View File

@ -248,8 +248,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
/* compute the qr factorization of the jacobian. */
ei_qrfac<Scalar>(m, n, fjac.data(), fjac.rows(), true, ipvt.data(), wa1.data(), wa2.data());
ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
wa2 = fjac.colwise().blueNorm();
ei_qrfac<Scalar>(m, n, fjac.data(), fjac.rows(), true, ipvt.data(), wa1.data());
ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convention (1->n), convert it to c (0->n-1)
/* on the first iteration and if mode is 1, scale according */
/* to the norms of the columns of the initial jacobian. */
@ -537,8 +538,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
}
if (sing) {
ipvt.cwise()+=1;
ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), true, ipvt.data(), wa1.data(), wa2.data());
ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
wa2 = fjac.colwise().blueNorm();
ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), true, ipvt.data(), wa1.data());
ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convention (1->n), convert it to c (0->n-1)
for (j = 0; j < n; ++j) {
if (fjac(j,j) != 0.) {
sum = 0.;

View File

@ -1,8 +1,7 @@
template <typename Scalar>
void ei_qrfac(int m, int n, Scalar *a, int
lda, int pivot, int *ipvt, Scalar *rdiag,
Scalar *acnorm)
lda, int pivot, int *ipvt, Scalar *rdiag)
{
/* System generated locals */
int a_dim1, a_offset;
@ -18,7 +17,6 @@ void ei_qrfac(int m, int n, Scalar *a, int
Matrix< Scalar, Dynamic, 1 > wa(n+1);
/* Parameter adjustments */
--acnorm;
--rdiag;
a_dim1 = lda;
a_offset = 1 + a_dim1 * 1;
@ -31,13 +29,10 @@ void ei_qrfac(int m, int n, Scalar *a, int
/* compute the initial column norms and initialize several arrays. */
for (j = 1; j <= n; ++j) {
acnorm[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm();
rdiag[j] = acnorm[j];
rdiag[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm();
wa[j] = rdiag[j];
if (pivot) {
if (pivot)
ipvt[j] = j;
}
/* L10: */
}
/* reduce a to r with householder transformations. */