TABLE OF CONTENTS
- ABINIT/m_levenberg_marquardt
- m_levenberg_marquardt/enorm
- m_levenberg_marquardt/fdjac2
- m_levenberg_marquardt/lm_fit_print_info
- m_levenberg_marquardt/lmdif
- m_levenberg_marquardt/lmdif1
- m_levenberg_marquardt/lmpar
- m_levenberg_marquardt/qrfac
- m_levenberg_marquardt/qrsolv
ABINIT/m_levenberg_marquardt [ Modules ]
NAME
m_levenberg_marquardt
FUNCTION
Module containing routines for performing Levenberg-Marquardt nonlinear least-squares fit. These routines are the conversions to F90 of the Minpack F77 routines by Alan Miller TODO: For some reason, the routine lmder, which uses analytic evaluation of the Jacobian, did not work in testing (bug?). lmdif proved stable, but potentially uses more funtion evaluations.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MS) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
23 #if defined HAVE_CONFIG_H 24 #include "config.h" 25 #endif 26 27 #include "abi_common.h" 28 29 module m_levenberg_marquardt 30 31 use defs_basis 32 use m_abicore 33 ! MINPACK routines which are used by both LMDIF & LMDER 34 ! 25 October 2001: 35 ! Changed INTENT of iflag in several places to IN OUT. 36 ! Changed INTENT of fvec to IN OUT in user routine FCN. 37 ! Removed arguments diag and qtv from LMDIF & LMDER. 38 ! Replaced several DO loops with array operations. 39 ! amiller @ bigpond.net.au 40 41 implicit none 42 43 private 44 45 public :: lmdif1 46 public :: lmdif 47 public :: lm_fit_print_info 48 49 CONTAINS !==============================================================================
m_levenberg_marquardt/enorm [ Functions ]
[ Top ] [ m_levenberg_marquardt ] [ Functions ]
NAME
FUNCTION
INPUTS
OUTPUT
NOTES
SOURCE
1317 function enorm(n,x) result(fn_val) 1318 1319 ! Code converted using TO_F90 by Alan Miller 1320 ! Date: 1999-12-09 Time: 12:45:34 1321 1322 !Arguments ------------------------------------ 1323 !scalars 1324 real (dp) :: fn_val 1325 integer, intent(in) :: n 1326 !arrays 1327 real (dp), intent(in) :: x(n) 1328 1329 !Local variables------------------------------- 1330 integer :: i 1331 real (dp) :: agiant, floatn, s1, s2, s3, xabs, x1max, x3max 1332 real (dp), parameter :: one = 1.0_dp, zero = 0.0_dp, rdwarf = 3.834e-20_dp, rgiant = 1.304e+19_dp 1333 ! ********************************************************************* 1334 1335 ! function enorm 1336 1337 ! given an n-vector x, this function calculates the euclidean norm of x. 1338 1339 ! the euclidean norm is computed by accumulating the sum of squares in 1340 ! three different sums. The sums of squares for the small and large 1341 ! components are scaled so that no overflows occur. Non-destructive 1342 ! underflows are permitted. Underflows and overflows do not occur in the 1343 ! computation of the unscaled sum of squares for the intermediate 1344 ! components. The definitions of small, intermediate and large components 1345 ! depend on two constants, rdwarf and rgiant. The main restrictions on 1346 ! these constants are that rdwarf**2 not underflow and rgiant**2 not 1347 ! overflow. The constants given here are suitable for every known computer. 1348 1349 ! the function statement is 1350 1351 ! REAL (dp) function enorm(n,x) 1352 1353 ! where 1354 1355 ! n is a positive integer input variable. 1356 1357 ! x is an input array of length n. 1358 1359 ! subprograms called 1360 1361 ! fortran-supplied ... ABS,SQRT 1362 1363 ! argonne national laboratory. minpack project. march 1980. 1364 ! burton s. garbow, kenneth e. hillstrom, jorge j. more 1365 1366 ! ********** 1367 1368 1369 s1 = zero 1370 s2 = zero 1371 s3 = zero 1372 x1max = zero 1373 x3max = zero 1374 floatn = n 1375 agiant = rgiant/floatn 1376 DO i = 1, n 1377 xabs = ABS(x(i)) 1378 IF (xabs > rdwarf .AND. xabs < agiant) GO TO 70 1379 IF (xabs <= rdwarf) GO TO 30 1380 1381 ! sum for large components. 1382 1383 IF (xabs <= x1max) GO TO 10 1384 s1 = one + s1*(x1max/xabs)**2 1385 x1max = xabs 1386 GO TO 20 1387 1388 10 s1 = s1 + (xabs/x1max)**2 1389 1390 20 GO TO 60 1391 1392 ! sum for small components. 1393 1394 30 IF (xabs <= x3max) GO TO 40 1395 s3 = one + s3*(x3max/xabs)**2 1396 x3max = xabs 1397 GO TO 60 1398 1399 40 IF (xabs /= zero) s3 = s3 + (xabs/x3max)**2 1400 1401 60 CYCLE 1402 1403 ! sum for intermediate components. 1404 1405 70 s2 = s2 + xabs**2 1406 END DO 1407 1408 ! calculation of norm. 1409 1410 IF (s1 == zero) GO TO 100 1411 fn_val = x1max*SQRT(s1 + (s2/x1max)/x1max) 1412 GO TO 120 1413 1414 100 IF (s2 == zero) GO TO 110 1415 IF (s2 >= x3max) fn_val = SQRT(s2*(one + (x3max/s2)*(x3max*s3))) 1416 IF (s2 < x3max) fn_val = SQRT(x3max*((s2/x3max) + (x3max*s3))) 1417 GO TO 120 1418 1419 110 fn_val = x3max*SQRT(s3) 1420 1421 120 RETURN 1422 1423 ! last card of function enorm. 1424 1425 END FUNCTION enorm
m_levenberg_marquardt/fdjac2 [ Functions ]
[ Top ] [ m_levenberg_marquardt ] [ Functions ]
NAME
FUNCTION
INPUTS
OUTPUT
NOTES
SOURCE
1443 SUBROUTINE fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn) 1444 1445 ! Code converted using TO_F90 by Alan Miller 1446 ! Date: 1999-12-09 Time: 12:45:44 1447 1448 ! N.B. Arguments LDFJAC & WA have been removed. 1449 integer, intent(in) :: m 1450 integer, intent(in) :: n 1451 real (dp), intent(in out) :: x(n) 1452 real (dp), intent(in) :: fvec(m) 1453 real (dp), intent(out) :: fjac(:,:) ! fjac(ldfjac,n) 1454 integer, intent(in out) :: iflag 1455 real (dp), intent(in) :: epsfcn 1456 1457 interface 1458 subroutine fcn(m, n, x, fvec, iflag, y, nfqre, nfqim, bsign, csign, multi_x_exp) 1459 implicit none 1460 integer, parameter :: dp = KIND(1.0d0) 1461 integer, intent(in) :: m, n 1462 real (dp), intent(in) :: x(n) 1463 real (dp), intent(in out) :: fvec(m) 1464 integer, intent(in out) :: iflag 1465 integer, optional, intent(in) :: nfqre,nfqim,bsign,csign,multi_x_exp 1466 real (dp), optional, intent(in) :: y(m) 1467 end subroutine fcn 1468 end interface 1469 1470 ! ********** 1471 1472 ! subroutine fdjac2 1473 1474 ! this subroutine computes a forward-difference approximation 1475 ! to the m by n jacobian matrix associated with a specified 1476 ! problem of m functions in n variables. 1477 1478 ! the subroutine statement is 1479 1480 ! subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa) 1481 1482 ! where 1483 1484 ! fcn is the name of the user-supplied subroutine which calculates the 1485 ! functions. fcn must be declared in an external statement in the user 1486 ! calling program, and should be written as follows. 1487 1488 ! subroutine fcn(m,n,x,fvec,iflag) 1489 ! integer m,n,iflag 1490 ! REAL (dp) x(n),fvec(m) 1491 ! ---------- 1492 ! calculate the functions at x and 1493 ! return this vector in fvec. 1494 ! ---------- 1495 ! return 1496 ! end 1497 1498 ! the value of iflag should not be changed by fcn unless 1499 ! the user wants to terminate execution of fdjac2. 1500 ! in this case set iflag to a negative integer. 1501 1502 ! m is a positive integer input variable set to the number of functions. 1503 1504 ! n is a positive integer input variable set to the number of variables. 1505 ! n must not exceed m. 1506 1507 ! x is an input array of length n. 1508 1509 ! fvec is an input array of length m which must contain the 1510 ! functions evaluated at x. 1511 1512 ! fjac is an output m by n array which contains the 1513 ! approximation to the jacobian matrix evaluated at x. 1514 1515 ! ldfjac is a positive integer input variable not less than m 1516 ! which specifies the leading dimension of the array fjac. 1517 1518 ! iflag is an integer variable which can be used to terminate 1519 ! the execution of fdjac2. see description of fcn. 1520 1521 ! epsfcn is an input variable used in determining a suitable step length 1522 ! for the forward-difference approximation. This approximation assumes 1523 ! that the relative errors in the functions are of the order of epsfcn. 1524 ! If epsfcn is less than the machine precision, it is assumed that the 1525 ! relative errors in the functions are of the order of the machine 1526 ! precision. 1527 1528 ! wa is a work array of length m. 1529 1530 ! subprograms called 1531 1532 ! user-supplied ...... fcn 1533 1534 ! minpack-supplied ... dpmpar 1535 1536 ! fortran-supplied ... ABS,MAX,SQRT 1537 1538 ! argonne national laboratory. minpack project. march 1980. 1539 ! burton s. garbow, kenneth e. hillstrom, jorge j. more 1540 1541 ! ********** 1542 integer :: j 1543 real (dp) :: eps, epsmch, h, temp, wa(m) 1544 real (dp), parameter :: zero = 0.0_dp 1545 1546 ! epsmch is the machine precision. 1547 1548 epsmch = EPSILON(zero) 1549 1550 eps = SQRT(MAX(epsfcn, epsmch)) 1551 DO j = 1, n 1552 temp = x(j) 1553 h = eps*ABS(temp) 1554 IF (h == zero) h = eps 1555 x(j) = temp + h 1556 CALL fcn(m, n, x, wa, iflag) 1557 IF (iflag < 0) EXIT 1558 x(j) = temp 1559 fjac(1:m,j) = (wa(1:m) - fvec(1:m))/h 1560 END DO 1561 1562 RETURN 1563 1564 ! last card of subroutine fdjac2. 1565 1566 end subroutine fdjac2
m_levenberg_marquardt/lm_fit_print_info [ Functions ]
[ Top ] [ m_levenberg_marquardt ] [ Functions ]
NAME
FUNCTION
INPUTS
OUTPUT
NOTES
SOURCE
1581 subroutine lm_fit_print_info(info,msg) 1582 1583 integer, intent(in) :: info 1584 character(len=500) :: msg 1585 1586 select case (info) 1587 case (:-1) 1588 write(msg,'(a,i0)') 'STOP LM fit: fcn in LM fit returned info = ', -info 1589 case (0) 1590 write(msg,'(a)') 'ERROR LM fit: improper values for input parameters in LM fit' 1591 case (1:3) 1592 write(msg,'(a)') 'LM fit converged' 1593 case (4) 1594 write(msg,'(a)') 'ERROR LM fit: residuals orthogonal to the jacobian, there may be an error in fcn.' 1595 case (5) 1596 write(msg,'(a)') 'ERROR LM fit: too many calls to fcn, either slow convergence, or an error in opt.' 1597 case (6:7) 1598 write(msg,'(a)') 'ERROR LM fit: tol was set too small' 1599 case default 1600 write(msg,'(a,i0)') 'ERROR LM fit: unknown info = ',info 1601 end select 1602 1603 return 1604 1605 end subroutine lm_fit_print_info 1606 1607 end module m_levenberg_marquardt
m_levenberg_marquardt/lmdif [ Functions ]
[ Top ] [ m_levenberg_marquardt ] [ Functions ]
NAME
FUNCTION
INPUTS
OUTPUT
NOTES
SOURCE
231 SUBROUTINE lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, & 232 mode, factor, nprint, info, nfev, fjac, ipvt) 233 234 ! N.B. Arguments LDFJAC, DIAG, QTF, WA1, WA2, WA3 & WA4 have been removed. 235 integer, intent(in) :: m 236 integer, intent(in) :: n 237 real (dp), intent(in out) :: x(:) 238 real (dp), intent(out) :: fvec(:) 239 real (dp), intent(in) :: ftol 240 real (dp), intent(in) :: xtol 241 real (dp), intent(in out) :: gtol 242 integer, intent(in out) :: maxfev 243 real (dp), intent(in out) :: epsfcn 244 integer, intent(in) :: mode 245 real (dp), intent(in) :: factor 246 integer, intent(in) :: nprint 247 integer, intent(out) :: info 248 integer, intent(out) :: nfev 249 real (dp), intent(out) :: fjac(:,:) ! fjac(ldfjac,n) 250 integer, intent(out) :: ipvt(:) 251 252 ! external fcn 253 254 interface 255 subroutine fcn(m, n, x, fvec, iflag, y, nfqre, nfqim, bsign, csign, multi_x_exp) 256 implicit none 257 integer, parameter :: dp = KIND(1.0d0) 258 integer, intent(in) :: m, n 259 real (dp), intent(in) :: x(n) 260 real (dp), intent(in out) :: fvec(m) 261 integer, intent(in out) :: iflag 262 integer, optional, intent(in) :: nfqre,nfqim,bsign,csign,multi_x_exp 263 real (dp), optional, intent(in) :: y(m) 264 end subroutine fcn 265 end interface 266 267 ! ********** 268 269 ! subroutine lmdif 270 271 ! The purpose of lmdif is to minimize the sum of the squares of m nonlinear 272 ! functions in n variables by a modification of the Levenberg-Marquardt 273 ! algorithm. The user must provide a subroutine which calculates the 274 ! functions. The jacobian is then calculated by a forward-difference 275 ! approximation. 276 277 ! the subroutine statement is 278 279 ! subroutine lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, 280 ! diag, mode, factor, nprint, info, nfev, fjac, 281 ! ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4) 282 283 ! N.B. 7 of these arguments have been removed in this version. 284 285 ! where 286 287 ! fcn is the name of the user-supplied subroutine which calculates the 288 ! functions. fcn must be declared in an external statement in the user 289 ! calling program, and should be written as follows. 290 291 ! subroutine fcn(m, n, x, fvec, iflag) 292 ! integer m, n, iflag 293 ! REAL (dp) x(:), fvec(m) 294 ! ---------- 295 ! calculate the functions at x and return this vector in fvec. 296 ! ---------- 297 ! return 298 ! end 299 300 ! the value of iflag should not be changed by fcn unless 301 ! the user wants to terminate execution of lmdif. 302 ! in this case set iflag to a negative integer. 303 304 ! m is a positive integer input variable set to the number of functions. 305 306 ! n is a positive integer input variable set to the number of variables. 307 ! n must not exceed m. 308 309 ! x is an array of length n. On input x must contain an initial estimate 310 ! of the solution vector. On output x contains the final estimate of the 311 ! solution vector. 312 313 ! fvec is an output array of length m which contains 314 ! the functions evaluated at the output x. 315 316 ! ftol is a nonnegative input variable. Termination occurs when both the 317 ! actual and predicted relative reductions in the sum of squares are at 318 ! most ftol. Therefore, ftol measures the relative error desired 319 ! in the sum of squares. 320 321 ! xtol is a nonnegative input variable. Termination occurs when the 322 ! relative error between two consecutive iterates is at most xtol. 323 ! Therefore, xtol measures the relative error desired in the approximate 324 ! solution. 325 326 ! gtol is a nonnegative input variable. Termination occurs when the cosine 327 ! of the angle between fvec and any column of the jacobian is at most 328 ! gtol in absolute value. Therefore, gtol measures the orthogonality 329 ! desired between the function vector and the columns of the jacobian. 330 331 ! maxfev is a positive integer input variable. Termination occurs when the 332 ! number of calls to fcn is at least maxfev by the end of an iteration. 333 334 ! epsfcn is an input variable used in determining a suitable step length 335 ! for the forward-difference approximation. This approximation assumes 336 ! that the relative errors in the functions are of the order of epsfcn. 337 ! If epsfcn is less than the machine precision, it is assumed that the 338 ! relative errors in the functions are of the order of the machine 339 ! precision. 340 341 ! diag is an array of length n. If mode = 1 (see below), diag is 342 ! internally set. If mode = 2, diag must contain positive entries that 343 ! serve as multiplicative scale factors for the variables. 344 345 ! mode is an integer input variable. If mode = 1, the variables will be 346 ! scaled internally. If mode = 2, the scaling is specified by the input 347 ! diag. other values of mode are equivalent to mode = 1. 348 349 ! factor is a positive input variable used in determining the initial step 350 ! bound. This bound is set to the product of factor and the euclidean 351 ! norm of diag*x if nonzero, or else to factor itself. In most cases 352 ! factor should lie in the interval (.1,100.). 100. is a generally 353 ! recommended value. 354 355 ! nprint is an integer input variable that enables controlled printing of 356 ! iterates if it is positive. In this case, fcn is called with iflag = 0 357 ! at the beginning of the first iteration and every nprint iterations 358 ! thereafter and immediately prior to return, with x and fvec available 359 ! for printing. If nprint is not positive, no special calls 360 ! of fcn with iflag = 0 are made. 361 362 ! info is an integer output variable. If the user has terminated 363 ! execution, info is set to the (negative) value of iflag. 364 ! See description of fcn. Otherwise, info is set as follows. 365 366 ! info = 0 improper input parameters. 367 368 ! info = 1 both actual and predicted relative reductions 369 ! in the sum of squares are at most ftol. 370 371 ! info = 2 relative error between two consecutive iterates <= xtol. 372 373 ! info = 3 conditions for info = 1 and info = 2 both hold. 374 375 ! info = 4 the cosine of the angle between fvec and any column of 376 ! the Jacobian is at most gtol in absolute value. 377 378 ! info = 5 number of calls to fcn has reached or exceeded maxfev. 379 380 ! info = 6 ftol is too small. no further reduction in 381 ! the sum of squares is possible. 382 383 ! info = 7 xtol is too small. no further improvement in 384 ! the approximate solution x is possible. 385 386 ! info = 8 gtol is too small. fvec is orthogonal to the 387 ! columns of the jacobian to machine precision. 388 389 ! nfev is an integer output variable set to the number of calls to fcn. 390 391 ! fjac is an output m by n array. the upper n by n submatrix 392 ! of fjac contains an upper triangular matrix r with 393 ! diagonal elements of nonincreasing magnitude such that 394 395 ! t t t 396 ! p *(jac *jac)*p = r *r, 397 398 ! where p is a permutation matrix and jac is the final calculated 399 ! Jacobian. Column j of p is column ipvt(j) (see below) of the 400 ! identity matrix. the lower trapezoidal part of fjac contains 401 ! information generated during the computation of r. 402 403 ! ldfjac is a positive integer input variable not less than m 404 ! which specifies the leading dimension of the array fjac. 405 406 ! ipvt is an integer output array of length n. ipvt defines a permutation 407 ! matrix p such that jac*p = q*r, where jac is the final calculated 408 ! jacobian, q is orthogonal (not stored), and r is upper triangular 409 ! with diagonal elements of nonincreasing magnitude. 410 ! Column j of p is column ipvt(j) of the identity matrix. 411 412 ! qtf is an output array of length n which contains 413 ! the first n elements of the vector (q transpose)*fvec. 414 415 ! wa1, wa2, and wa3 are work arrays of length n. 416 417 ! wa4 is a work array of length m. 418 419 ! subprograms called 420 421 ! user-supplied ...... fcn 422 423 ! minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac 424 425 ! fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod 426 427 ! argonne national laboratory. minpack project. march 1980. 428 ! burton s. garbow, kenneth e. hillstrom, jorge j. more 429 430 ! ********** 431 integer :: i, iflag, iter, j, l 432 real (dp) :: actred, delta, dirder, epsmch, fnorm, fnorm1, gnorm 433 real(dp) :: par, pnorm, prered, ratio, sum, temp, temp1, temp2, xnorm 434 real (dp) :: diag(n), qtf(n), wa1(n), wa2(n), wa3(n), wa4(m) 435 real (dp), parameter :: one = 1.0_dp, p1 = 0.1_dp, p5 = 0.5_dp 436 real(dp),parameter :: p25 = 0.25_dp, p75 = 0.75_dp, p0001 = 0.0001_dp, zero = 0.0_dp 437 438 ! epsmch is the machine precision. 439 440 epsmch = EPSILON(zero) 441 442 info = 0 443 iflag = 0 444 nfev = 0 445 446 ! check the input parameters for errors. 447 448 IF (n <= 0 .OR. m < n .OR. ftol < zero .OR. xtol < zero .OR. gtol < zero & 449 & .OR. maxfev <= 0 .OR. factor <= zero) GO TO 300 450 IF (mode /= 2) GO TO 20 451 DO j = 1, n 452 IF (diag(j) <= zero) GO TO 300 453 END DO 454 455 ! evaluate the function at the starting point and calculate its norm. 456 457 20 iflag = 1 458 CALL fcn(m, n, x, fvec, iflag) 459 nfev = 1 460 IF (iflag < 0) GO TO 300 461 fnorm = enorm(m, fvec) 462 463 ! initialize levenberg-marquardt parameter and iteration counter. 464 465 par = zero 466 iter = 1 467 468 ! beginning of the outer loop. 469 470 ! calculate the jacobian matrix. 471 472 30 iflag = 2 473 CALL fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn) 474 nfev = nfev + n 475 IF (iflag < 0) GO TO 300 476 477 ! If requested, call fcn to enable printing of iterates. 478 479 IF (nprint <= 0) GO TO 40 480 iflag = 0 481 IF (MOD(iter-1,nprint) == 0) THEN 482 CALL fcn(m, n, x, fvec, iflag) 483 END IF 484 IF (iflag < 0) GO TO 300 485 486 ! Compute the qr factorization of the jacobian. 487 488 40 CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2) 489 490 ! On the first iteration and if mode is 1, scale according 491 ! to the norms of the columns of the initial jacobian. 492 493 IF (iter /= 1) GO TO 80 494 IF (mode == 2) GO TO 60 495 DO j = 1, n 496 diag(j) = wa2(j) 497 IF (wa2(j) == zero) diag(j) = one 498 END DO 499 500 ! On the first iteration, calculate the norm of the scaled x 501 ! and initialize the step bound delta. 502 503 60 wa3(1:n) = diag(1:n)*x(1:n) 504 xnorm = enorm(n, wa3) 505 delta = factor*xnorm 506 IF (delta == zero) delta = factor 507 508 ! Form (q transpose)*fvec and store the first n components in qtf. 509 510 80 wa4(1:m) = fvec(1:m) 511 DO j = 1, n 512 IF (fjac(j,j) == zero) GO TO 120 513 sum = DOT_PRODUCT( fjac(j:m,j), wa4(j:m) ) 514 temp = -sum/fjac(j,j) 515 DO i = j, m 516 wa4(i) = wa4(i) + fjac(i,j)*temp 517 END DO 518 120 fjac(j,j) = wa1(j) 519 qtf(j) = wa4(j) 520 END DO 521 522 ! compute the norm of the scaled gradient. 523 524 gnorm = zero 525 IF (fnorm == zero) GO TO 170 526 DO j = 1, n 527 l = ipvt(j) 528 IF (wa2(l) == zero) CYCLE 529 sum = zero 530 DO i = 1, j 531 sum = sum + fjac(i,j)*(qtf(i)/fnorm) 532 END DO 533 gnorm = MAX(gnorm, ABS(sum/wa2(l))) 534 END DO 535 536 ! test for convergence of the gradient norm. 537 538 170 IF (gnorm <= gtol) info = 4 539 IF (info /= 0) GO TO 300 540 541 ! rescale if necessary. 542 543 IF (mode == 2) GO TO 200 544 DO j = 1, n 545 diag(j) = MAX(diag(j), wa2(j)) 546 END DO 547 548 ! beginning of the inner loop. 549 550 ! determine the Levenberg-Marquardt parameter. 551 552 200 CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2) 553 554 ! store the direction p and x + p. calculate the norm of p. 555 556 DO j = 1, n 557 wa1(j) = -wa1(j) 558 wa2(j) = x(j) + wa1(j) 559 wa3(j) = diag(j)*wa1(j) 560 END DO 561 pnorm = enorm(n, wa3) 562 563 ! on the first iteration, adjust the initial step bound. 564 565 IF (iter == 1) delta = MIN(delta, pnorm) 566 567 ! evaluate the function at x + p and calculate its norm. 568 569 iflag = 1 570 CALL fcn(m, n, wa2, wa4, iflag) 571 nfev = nfev + 1 572 IF (iflag < 0) GO TO 300 573 fnorm1 = enorm(m, wa4) 574 575 ! compute the scaled actual reduction. 576 577 actred = -one 578 IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2 579 580 ! Compute the scaled predicted reduction and 581 ! the scaled directional derivative. 582 583 DO j = 1, n 584 wa3(j) = zero 585 l = ipvt(j) 586 temp = wa1(l) 587 DO i = 1, j 588 wa3(i) = wa3(i) + fjac(i,j)*temp 589 END DO 590 END DO 591 temp1 = enorm(n,wa3)/fnorm 592 temp2 = (SQRT(par)*pnorm)/fnorm 593 prered = temp1**2 + temp2**2/p5 594 dirder = -(temp1**2 + temp2**2) 595 596 ! compute the ratio of the actual to the predicted reduction. 597 598 ratio = zero 599 IF (prered /= zero) ratio = actred/prered 600 601 ! update the step bound. 602 603 IF (ratio <= p25) THEN 604 IF (actred >= zero) temp = p5 605 IF (actred < zero) temp = p5*dirder/(dirder + p5*actred) 606 IF (p1*fnorm1 >= fnorm .OR. temp < p1) temp = p1 607 delta = temp*MIN(delta,pnorm/p1) 608 par = par/temp 609 ELSE 610 IF (par /= zero .AND. ratio < p75) GO TO 260 611 delta = pnorm/p5 612 par = p5*par 613 END IF 614 615 ! test for successful iteration. 616 617 260 IF (ratio < p0001) GO TO 290 618 619 ! successful iteration. update x, fvec, and their norms. 620 621 DO j = 1, n 622 x(j) = wa2(j) 623 wa2(j) = diag(j)*x(j) 624 END DO 625 fvec(1:m) = wa4(1:m) 626 xnorm = enorm(n, wa2) 627 fnorm = fnorm1 628 iter = iter + 1 629 630 ! tests for convergence. 631 632 290 IF (ABS(actred) <= ftol .AND. prered <= ftol .AND. p5*ratio <= one) info = 1 633 IF (delta <= xtol*xnorm) info = 2 634 IF (ABS(actred) <= ftol .AND. prered <= ftol & 635 & .AND. p5*ratio <= one .AND. info == 2) info = 3 636 IF (info /= 0) GO TO 300 637 638 ! tests for termination and stringent tolerances. 639 640 IF (nfev >= maxfev) info = 5 641 IF (ABS(actred) <= epsmch .AND. prered <= epsmch & 642 & .AND. p5*ratio <= one) info = 6 643 IF (delta <= epsmch*xnorm) info = 7 644 IF (gnorm <= epsmch) info = 8 645 IF (info /= 0) GO TO 300 646 647 ! end of the inner loop. repeat if iteration unsuccessful. 648 649 IF (ratio < p0001) GO TO 200 650 651 ! end of the outer loop. 652 653 GO TO 30 654 655 ! termination, either normal or user imposed. 656 657 300 IF (iflag < 0) info = iflag 658 iflag = 0 659 IF (nprint > 0) THEN 660 CALL fcn(m, n, x, fvec, iflag) 661 END IF 662 RETURN 663 664 ! last card of subroutine lmdif. 665 666 end subroutine lmdif
m_levenberg_marquardt/lmdif1 [ Functions ]
[ Top ] [ m_levenberg_marquardt ] [ Functions ]
NAME
FUNCTION
INPUTS
OUTPUT
NOTES
SOURCE
65 subroutine lmdif1(fcn, m, n, x, fvec, tol, info, iwa) 66 67 ! code converted using to_f90 by alan miller 68 ! date: 1999-12-11 time: 00:51:44 69 70 ! n.b. arguments wa & lwa have been removed. 71 integer, intent(in) :: m 72 integer, intent(in) :: n 73 real (dp), intent(in out) :: x(:) 74 real (dp), intent(out) :: fvec(:) 75 real (dp), intent(in) :: tol 76 integer, intent(out) :: info 77 integer, intent(out) :: iwa(:) 78 79 ! external fcn 80 81 interface 82 subroutine fcn(m, n, x, fvec, iflag, y, nfqre, nfqim, bsign, csign, multi_x_exp) 83 implicit none 84 integer, parameter :: dp = KIND(1.0d0) 85 integer, intent(in) :: m, n 86 real (dp), intent(in) :: x(n) 87 real (dp), intent(in out) :: fvec(m) 88 integer, intent(in out) :: iflag 89 integer, optional, intent(in) :: nfqre,nfqim,bsign,csign,multi_x_exp 90 real (dp), optional, intent(in) :: y(m) 91 end subroutine fcn 92 end interface 93 94 ! ********** 95 96 ! subroutine lmdif1 97 98 ! The purpose of lmdif1 is to minimize the sum of the squares of m nonlinear 99 ! functions in n variables by a modification of the Levenberg-Marquardt 100 ! algorithm. This is done by using the more general least-squares 101 ! solver lmdif. The user must provide a subroutine which calculates the 102 ! functions. The jacobian is then calculated by a forward-difference 103 ! approximation. 104 105 ! the subroutine statement is 106 107 ! subroutine lmdif1(fcn, m, n, x, fvec, tol, info, iwa) 108 109 ! where 110 111 ! fcn is the name of the user-supplied subroutine which calculates 112 ! the functions. fcn must be declared in an external statement in the 113 ! user calling program, and should be written as follows. 114 115 ! subroutine fcn(m, n, x, fvec, iflag) 116 ! integer m, n, iflag 117 ! REAL (dp) x(n), fvec(m) 118 ! ---------- 119 ! calculate the functions at x and return this vector in fvec. 120 ! ---------- 121 ! return 122 ! end 123 124 ! the value of iflag should not be changed by fcn unless 125 ! the user wants to terminate execution of lmdif1. 126 ! In this case set iflag to a negative integer. 127 128 ! m is a positive integer input variable set to the number of functions. 129 130 ! n is a positive integer input variable set to the number of variables. 131 ! n must not exceed m. 132 133 ! x is an array of length n. On input x must contain an initial estimate 134 ! of the solution vector. On output x contains the final estimate of 135 ! the solution vector. 136 137 ! fvec is an output array of length m which contains 138 ! the functions evaluated at the output x. 139 140 ! tol is a nonnegative input variable. Termination occurs when the 141 ! algorithm estimates either that the relative error in the sum of 142 ! squares is at most tol or that the relative error between x and the 143 ! solution is at most tol. 144 145 ! info is an integer output variable. If the user has terminated execution, 146 ! info is set to the (negative) value of iflag. See description of fcn. 147 ! Otherwise, info is set as follows. 148 149 ! info = 0 improper input parameters. 150 151 ! info = 1 algorithm estimates that the relative error 152 ! in the sum of squares is at most tol. 153 154 ! info = 2 algorithm estimates that the relative error 155 ! between x and the solution is at most tol. 156 157 ! info = 3 conditions for info = 1 and info = 2 both hold. 158 159 ! info = 4 fvec is orthogonal to the columns of the 160 ! jacobian to machine precision. 161 162 ! info = 5 number of calls to fcn has reached or exceeded 200*(n+1). 163 164 ! info = 6 tol is too small. no further reduction in 165 ! the sum of squares is possible. 166 167 ! info = 7 tol is too small. No further improvement in 168 ! the approximate solution x is possible. 169 170 ! iwa is an integer work array of length n. 171 172 ! wa is a work array of length lwa. 173 174 ! lwa is a positive integer input variable not less than m*n+5*n+m. 175 176 ! subprograms called 177 178 ! user-supplied ...... fcn 179 180 ! minpack-supplied ... lmdif 181 182 ! argonne national laboratory. minpack project. march 1980. 183 ! burton s. garbow, kenneth e. hillstrom, jorge j. more 184 185 ! ********** 186 integer :: maxfev, mode, nfev, nprint 187 real (dp) :: epsfcn, ftol, gtol, xtol, fjac(m,n) 188 real (dp), parameter :: factor = 100._dp , zero = 0.0_dp 189 190 info = 0 191 192 ! check the input parameters for errors. 193 194 IF (n <= 0 .OR. m < n .OR. tol < zero) GO TO 10 195 196 ! call lmdif. 197 198 maxfev = 200*(n + 1) 199 ftol = tol 200 xtol = tol 201 gtol = zero 202 epsfcn = zero 203 mode = 1 204 nprint = 0 205 CALL lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, & 206 & mode, factor, nprint, info, nfev, fjac, iwa) 207 IF (info == 8) info = 4 208 209 10 RETURN 210 211 ! last card of subroutine lmdif1. 212 213 end subroutine lmdif1
m_levenberg_marquardt/lmpar [ Functions ]
[ Top ] [ m_levenberg_marquardt ] [ Functions ]
NAME
FUNCTION
INPUTS
OUTPUT
NOTES
SOURCE
684 subroutine lmpar(n, r, ipvt, diag, qtb, delta, par, x, sdiag) 685 686 ! Code converted using TO_F90 by Alan Miller 687 ! Date: 1999-12-09 Time: 12:46:12 688 689 ! N.B. Arguments LDR, WA1 & WA2 have been removed. 690 integer, intent(in) :: n 691 real (dp), intent(in out) :: r(:,:) 692 integer, intent(in) :: ipvt(:) 693 real (dp), intent(in) :: diag(:) 694 real (dp), intent(in) :: qtb(:) 695 real (dp), intent(in) :: delta 696 real (dp), intent(out) :: par 697 real (dp), intent(out) :: x(:) 698 real (dp), intent(out) :: sdiag(:) 699 700 ! ********** 701 702 ! subroutine lmpar 703 704 ! Given an m by n matrix a, an n by n nonsingular diagonal matrix d, 705 ! an m-vector b, and a positive number delta, the problem is to determine a 706 ! value for the parameter par such that if x solves the system 707 708 ! a*x = b , sqrt(par)*d*x = 0 , 709 710 ! in the least squares sense, and dxnorm is the euclidean 711 ! norm of d*x, then either par is zero and 712 713 ! (dxnorm-delta) <= 0.1*delta , 714 715 ! or par is positive and 716 717 ! abs(dxnorm-delta) <= 0.1*delta . 718 719 ! This subroutine completes the solution of the problem if it is provided 720 ! with the necessary information from the r factorization, with column 721 ! qpivoting, of a. That is, if a*p = q*r, where p is a permutation matrix, 722 ! q has orthogonal columns, and r is an upper triangular matrix with diagonal 723 ! elements of nonincreasing magnitude, then lmpar expects the full upper 724 ! triangle of r, the permutation matrix p, and the first n components of 725 ! (q transpose)*b. 726 ! On output lmpar also provides an upper triangular matrix s such that 727 728 ! t t t 729 ! p *(a *a + par*d*d)*p = s *s . 730 731 ! s is employed within lmpar and may be of separate interest. 732 733 ! Only a few iterations are generally needed for convergence of the algorithm. 734 ! If, however, the limit of 10 iterations is reached, then the output par 735 ! will contain the best value obtained so far. 736 737 ! the subroutine statement is 738 739 ! subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, wa1,wa2) 740 741 ! where 742 743 ! n is a positive integer input variable set to the order of r. 744 745 ! r is an n by n array. on input the full upper triangle 746 ! must contain the full upper triangle of the matrix r. 747 ! On output the full upper triangle is unaltered, and the 748 ! strict lower triangle contains the strict upper triangle 749 ! (transposed) of the upper triangular matrix s. 750 751 ! ldr is a positive integer input variable not less than n 752 ! which specifies the leading dimension of the array r. 753 754 ! ipvt is an integer input array of length n which defines the 755 ! permutation matrix p such that a*p = q*r. column j of p 756 ! is column ipvt(j) of the identity matrix. 757 758 ! diag is an input array of length n which must contain the 759 ! diagonal elements of the matrix d. 760 761 ! qtb is an input array of length n which must contain the first 762 ! n elements of the vector (q transpose)*b. 763 764 ! delta is a positive input variable which specifies an upper 765 ! bound on the euclidean norm of d*x. 766 767 ! par is a nonnegative variable. on input par contains an 768 ! initial estimate of the levenberg-marquardt parameter. 769 ! on output par contains the final estimate. 770 771 ! x is an output array of length n which contains the least 772 ! squares solution of the system a*x = b, sqrt(par)*d*x = 0, 773 ! for the output par. 774 775 ! sdiag is an output array of length n which contains the 776 ! diagonal elements of the upper triangular matrix s. 777 778 ! wa1 and wa2 are work arrays of length n. 779 780 ! subprograms called 781 782 ! minpack-supplied ... dpmpar,enorm,qrsolv 783 784 ! fortran-supplied ... ABS,MAX,MIN,SQRT 785 786 ! argonne national laboratory. minpack project. march 1980. 787 ! burton s. garbow, kenneth e. hillstrom, jorge j. more 788 789 ! ********** 790 integer :: iter, j, jm1, jp1, k, l, nsing 791 real (dp) :: dxnorm, dwarf, fp, gnorm, parc, parl, paru, sum, temp 792 real (dp) :: wa1(n), wa2(n) 793 real (dp), parameter :: p1 = 0.1_dp, p001 = 0.001_dp, zero = 0.0_dp 794 795 ! dwarf is the smallest positive magnitude. 796 797 dwarf = TINY(zero) 798 799 ! compute and store in x the gauss-newton direction. if the 800 ! jacobian is rank-deficient, obtain a least squares solution. 801 802 nsing = n 803 DO j = 1, n 804 wa1(j) = qtb(j) 805 IF (r(j,j) == zero .AND. nsing == n) nsing = j - 1 806 IF (nsing < n) wa1(j) = zero 807 END DO 808 809 DO k = 1, nsing 810 j = nsing - k + 1 811 wa1(j) = wa1(j)/r(j,j) 812 temp = wa1(j) 813 jm1 = j - 1 814 wa1(1:jm1) = wa1(1:jm1) - r(1:jm1,j)*temp 815 END DO 816 817 DO j = 1, n 818 l = ipvt(j) 819 x(l) = wa1(j) 820 END DO 821 822 ! initialize the iteration counter. 823 ! evaluate the function at the origin, and test 824 ! for acceptance of the gauss-newton direction. 825 826 iter = 0 827 wa2(1:n) = diag(1:n)*x(1:n) 828 dxnorm = enorm(n, wa2) 829 fp = dxnorm - delta 830 IF (fp <= p1*delta) GO TO 220 831 832 ! if the jacobian is not rank deficient, the newton 833 ! step provides a lower bound, parl, for the zero of 834 ! the function. Otherwise set this bound to zero. 835 836 parl = zero 837 IF (nsing < n) GO TO 120 838 DO j = 1, n 839 l = ipvt(j) 840 wa1(j) = diag(l)*(wa2(l)/dxnorm) 841 END DO 842 DO j = 1, n 843 sum = DOT_PRODUCT( r(1:j-1,j), wa1(1:j-1) ) 844 wa1(j) = (wa1(j) - sum)/r(j,j) 845 END DO 846 temp = enorm(n,wa1) 847 parl = ((fp/delta)/temp)/temp 848 849 ! calculate an upper bound, paru, for the zero of the function. 850 851 120 DO j = 1, n 852 sum = DOT_PRODUCT( r(1:j,j), qtb(1:j) ) 853 l = ipvt(j) 854 wa1(j) = sum/diag(l) 855 END DO 856 gnorm = enorm(n,wa1) 857 paru = gnorm/delta 858 IF (paru == zero) paru = dwarf/MIN(delta,p1) 859 860 ! if the input par lies outside of the interval (parl,paru), 861 ! set par to the closer endpoint. 862 863 par = MAX(par,parl) 864 par = MIN(par,paru) 865 IF (par == zero) par = gnorm/dxnorm 866 867 ! beginning of an iteration. 868 869 150 iter = iter + 1 870 871 ! evaluate the function at the current value of par. 872 873 IF (par == zero) par = MAX(dwarf, p001*paru) 874 temp = SQRT(par) 875 wa1(1:n) = temp*diag(1:n) 876 CALL qrsolv(n, r, ipvt, wa1, qtb, x, sdiag) 877 wa2(1:n) = diag(1:n)*x(1:n) 878 dxnorm = enorm(n, wa2) 879 temp = fp 880 fp = dxnorm - delta 881 882 ! if the function is small enough, accept the current value 883 ! of par. also test for the exceptional cases where parl 884 ! is zero or the number of iterations has reached 10. 885 886 IF (ABS(fp) <= p1*delta .OR. parl == zero .AND. fp <= temp & 887 & .AND. temp < zero .OR. iter == 10) GO TO 220 888 889 ! compute the newton correction. 890 891 DO j = 1, n 892 l = ipvt(j) 893 wa1(j) = diag(l)*(wa2(l)/dxnorm) 894 END DO 895 DO j = 1, n 896 wa1(j) = wa1(j)/sdiag(j) 897 temp = wa1(j) 898 jp1 = j + 1 899 wa1(jp1:n) = wa1(jp1:n) - r(jp1:n,j)*temp 900 END DO 901 temp = enorm(n,wa1) 902 parc = ((fp/delta)/temp)/temp 903 904 ! depending on the sign of the function, update parl or paru. 905 906 IF (fp > zero) parl = MAX(parl,par) 907 IF (fp < zero) paru = MIN(paru,par) 908 909 ! compute an improved estimate for par. 910 911 par = MAX(parl, par+parc) 912 913 ! end of an iteration. 914 915 GO TO 150 916 917 ! termination. 918 919 220 IF (iter == 0) par = zero 920 RETURN 921 922 ! last card of subroutine lmpar. 923 924 end subroutine lmpar
m_levenberg_marquardt/qrfac [ Functions ]
[ Top ] [ m_levenberg_marquardt ] [ Functions ]
NAME
FUNCTION
INPUTS
OUTPUT
NOTES
SOURCE
942 subroutine qrfac(m, n, a, pivot, ipvt, rdiag, acnorm) 943 944 ! Code converted using TO_F90 by Alan Miller 945 ! Date: 1999-12-09 Time: 12:46:17 946 947 ! N.B. Arguments LDA, LIPVT & WA have been removed. 948 integer, intent(in) :: m 949 integer, intent(in) :: n 950 real (dp), intent(in out) :: a(:,:) 951 logical, intent(in) :: pivot 952 integer, intent(out) :: ipvt(:) 953 real (dp), intent(out) :: rdiag(:) 954 real (dp), intent(out) :: acnorm(:) 955 956 ! ********** 957 958 ! subroutine qrfac 959 960 ! This subroutine uses Householder transformations with column pivoting 961 ! (optional) to compute a qr factorization of the m by n matrix a. 962 ! That is, qrfac determines an orthogonal matrix q, a permutation matrix p, 963 ! and an upper trapezoidal matrix r with diagonal elements of nonincreasing 964 ! magnitude, such that a*p = q*r. The householder transformation for 965 ! column k, k = 1,2,...,min(m,n), is of the form 966 967 ! t 968 ! i - (1/u(k))*u*u 969 970 ! where u has zeros in the first k-1 positions. The form of this 971 ! transformation and the method of pivoting first appeared in the 972 ! corresponding linpack subroutine. 973 974 ! the subroutine statement is 975 976 ! subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa) 977 978 ! N.B. 3 of these arguments have been omitted in this version. 979 980 ! where 981 982 ! m is a positive integer input variable set to the number of rows of a. 983 984 ! n is a positive integer input variable set to the number of columns of a. 985 986 ! a is an m by n array. On input a contains the matrix for 987 ! which the qr factorization is to be computed. On output 988 ! the strict upper trapezoidal part of a contains the strict 989 ! upper trapezoidal part of r, and the lower trapezoidal 990 ! part of a contains a factored form of q (the non-trivial 991 ! elements of the u vectors described above). 992 993 ! lda is a positive integer input variable not less than m 994 ! which specifies the leading dimension of the array a. 995 996 ! pivot is a logical input variable. If pivot is set true, 997 ! then column pivoting is enforced. If pivot is set false, 998 ! then no column pivoting is done. 999 1000 ! ipvt is an integer output array of length lipvt. ipvt 1001 ! defines the permutation matrix p such that a*p = q*r. 1002 ! Column j of p is column ipvt(j) of the identity matrix. 1003 ! If pivot is false, ipvt is not referenced. 1004 1005 ! lipvt is a positive integer input variable. If pivot is false, 1006 ! then lipvt may be as small as 1. If pivot is true, then 1007 ! lipvt must be at least n. 1008 1009 ! rdiag is an output array of length n which contains the 1010 ! diagonal elements of r. 1011 1012 ! acnorm is an output array of length n which contains the norms of the 1013 ! corresponding columns of the input matrix a. 1014 ! If this information is not needed, then acnorm can coincide with rdiag. 1015 1016 ! wa is a work array of length n. If pivot is false, then wa 1017 ! can coincide with rdiag. 1018 1019 ! subprograms called 1020 1021 ! minpack-supplied ... dpmpar,enorm 1022 1023 ! fortran-supplied ... MAX,SQRT,MIN 1024 1025 ! argonne national laboratory. minpack project. march 1980. 1026 ! burton s. garbow, kenneth e. hillstrom, jorge j. more 1027 1028 ! ********** 1029 integer :: i, j, jp1, k, kmax, minmn 1030 real (dp) :: ajnorm, epsmch, sum, temp, wa(n) 1031 real (dp), parameter :: one = 1.0_dp, p05 = 0.05_dp, zero = 0.0_dp 1032 1033 ! epsmch is the machine precision. 1034 1035 epsmch = EPSILON(zero) 1036 1037 ! compute the initial column norms and initialize several arrays. 1038 1039 DO j = 1, n 1040 acnorm(j) = enorm(m,a(1:,j)) 1041 rdiag(j) = acnorm(j) 1042 wa(j) = rdiag(j) 1043 IF (pivot) ipvt(j) = j 1044 END DO 1045 1046 ! Reduce a to r with Householder transformations. 1047 1048 minmn = MIN(m,n) 1049 DO j = 1, minmn 1050 IF (.NOT.pivot) GO TO 40 1051 1052 ! Bring the column of largest norm into the pivot position. 1053 1054 kmax = j 1055 DO k = j, n 1056 IF (rdiag(k) > rdiag(kmax)) kmax = k 1057 END DO 1058 IF (kmax == j) GO TO 40 1059 DO i = 1, m 1060 temp = a(i,j) 1061 a(i,j) = a(i,kmax) 1062 a(i,kmax) = temp 1063 END DO 1064 rdiag(kmax) = rdiag(j) 1065 wa(kmax) = wa(j) 1066 k = ipvt(j) 1067 ipvt(j) = ipvt(kmax) 1068 ipvt(kmax) = k 1069 1070 ! Compute the Householder transformation to reduce the 1071 ! j-th column of a to a multiple of the j-th unit vector. 1072 1073 40 ajnorm = enorm(m-j+1, a(j:,j)) 1074 IF (ajnorm == zero) CYCLE 1075 IF (a(j,j) < zero) ajnorm = -ajnorm 1076 a(j:m,j) = a(j:m,j)/ajnorm 1077 a(j,j) = a(j,j) + one 1078 1079 ! Apply the transformation to the remaining columns and update the norms. 1080 1081 jp1 = j + 1 1082 DO k = jp1, n 1083 sum = DOT_PRODUCT( a(j:m,j), a(j:m,k) ) 1084 temp = sum/a(j,j) 1085 a(j:m,k) = a(j:m,k) - temp*a(j:m,j) 1086 IF (.NOT.pivot .OR. rdiag(k) == zero) CYCLE 1087 temp = a(j,k)/rdiag(k) 1088 rdiag(k) = rdiag(k)*SQRT(MAX(zero, one-temp**2)) 1089 IF (p05*(rdiag(k)/wa(k))**2 > epsmch) CYCLE 1090 rdiag(k) = enorm(m-j, a(jp1:,k)) 1091 wa(k) = rdiag(k) 1092 END DO 1093 rdiag(j) = -ajnorm 1094 END DO 1095 RETURN 1096 1097 ! last card of subroutine qrfac. 1098 1099 end subroutine qrfac
m_levenberg_marquardt/qrsolv [ Functions ]
[ Top ] [ m_levenberg_marquardt ] [ Functions ]
NAME
FUNCTION
INPUTS
OUTPUT
NOTES
SOURCE
1117 SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag) 1118 1119 ! N.B. Arguments LDR & WA have been removed. 1120 integer, intent(in) :: n 1121 real (dp), intent(in out) :: r(:,:) 1122 integer, intent(in) :: ipvt(:) 1123 real (dp), intent(in) :: diag(:) 1124 real (dp), intent(in) :: qtb(:) 1125 real (dp), intent(out) :: x(:) 1126 real (dp), intent(out) :: sdiag(:) 1127 1128 ! ********** 1129 1130 ! subroutine qrsolv 1131 1132 ! Given an m by n matrix a, an n by n diagonal matrix d, and an m-vector b, 1133 ! the problem is to determine an x which solves the system 1134 1135 ! a*x = b , d*x = 0 , 1136 1137 ! in the least squares sense. 1138 1139 ! This subroutine completes the solution of the problem if it is provided 1140 ! with the necessary information from the qr factorization, with column 1141 ! pivoting, of a. That is, if a*p = q*r, where p is a permutation matrix, 1142 ! q has orthogonal columns, and r is an upper triangular matrix with diagonal 1143 ! elements of nonincreasing magnitude, then qrsolv expects the full upper 1144 ! triangle of r, the permutation matrix p, and the first n components of 1145 ! (q transpose)*b. The system a*x = b, d*x = 0, is then equivalent to 1146 1147 ! t t 1148 ! r*z = q *b , p *d*p*z = 0 , 1149 1150 ! where x = p*z. if this system does not have full rank, 1151 ! then a least squares solution is obtained. On output qrsolv 1152 ! also provides an upper triangular matrix s such that 1153 1154 ! t t t 1155 ! p *(a *a + d*d)*p = s *s . 1156 1157 ! s is computed within qrsolv and may be of separate interest. 1158 1159 ! the subroutine statement is 1160 1161 ! subroutine qrsolv(n, r, ldr, ipvt, diag, qtb, x, sdiag, wa) 1162 1163 ! N.B. Arguments LDR and WA have been removed in this version. 1164 1165 ! where 1166 1167 ! n is a positive integer input variable set to the order of r. 1168 1169 ! r is an n by n array. On input the full upper triangle must contain 1170 ! the full upper triangle of the matrix r. 1171 ! On output the full upper triangle is unaltered, and the strict lower 1172 ! triangle contains the strict upper triangle (transposed) of the 1173 ! upper triangular matrix s. 1174 1175 ! ldr is a positive integer input variable not less than n 1176 ! which specifies the leading dimension of the array r. 1177 1178 ! ipvt is an integer input array of length n which defines the 1179 ! permutation matrix p such that a*p = q*r. Column j of p 1180 ! is column ipvt(j) of the identity matrix. 1181 1182 ! diag is an input array of length n which must contain the 1183 ! diagonal elements of the matrix d. 1184 1185 ! qtb is an input array of length n which must contain the first 1186 ! n elements of the vector (q transpose)*b. 1187 1188 ! x is an output array of length n which contains the least 1189 ! squares solution of the system a*x = b, d*x = 0. 1190 1191 ! sdiag is an output array of length n which contains the 1192 ! diagonal elements of the upper triangular matrix s. 1193 1194 ! wa is a work array of length n. 1195 1196 ! subprograms called 1197 1198 ! fortran-supplied ... ABS,SQRT 1199 1200 ! argonne national laboratory. minpack project. march 1980. 1201 ! burton s. garbow, kenneth e. hillstrom, jorge j. more 1202 1203 ! ********** 1204 integer :: i, j, k, kp1, l, nsing 1205 real (dp) :: cos, cotan, qtbpj, sin, sum, tan, temp, wa(n) 1206 real (dp), parameter :: p5 = 0.5_dp, p25 = 0.25_dp, zero = 0.0_dp 1207 1208 ! Copy r and (q transpose)*b to preserve input and initialize s. 1209 ! In particular, save the diagonal elements of r in x. 1210 1211 DO j = 1, n 1212 r(j:n,j) = r(j,j:n) 1213 x(j) = r(j,j) 1214 wa(j) = qtb(j) 1215 END DO 1216 1217 ! Eliminate the diagonal matrix d using a givens rotation. 1218 1219 DO j = 1, n 1220 1221 ! Prepare the row of d to be eliminated, locating the 1222 ! diagonal element using p from the qr factorization. 1223 1224 l = ipvt(j) 1225 IF (diag(l) == zero) CYCLE 1226 sdiag(j:n) = zero 1227 sdiag(j) = diag(l) 1228 1229 ! The transformations to eliminate the row of d modify only a single 1230 ! element of (q transpose)*b beyond the first n, which is initially zero. 1231 1232 qtbpj = zero 1233 DO k = j, n 1234 1235 ! Determine a givens rotation which eliminates the 1236 ! appropriate element in the current row of d. 1237 1238 IF (sdiag(k) == zero) CYCLE 1239 IF (ABS(r(k,k)) < ABS(sdiag(k))) THEN 1240 cotan = r(k,k)/sdiag(k) 1241 SIN = p5/SQRT(p25 + p25*cotan**2) 1242 COS = SIN*cotan 1243 ELSE 1244 TAN = sdiag(k)/r(k,k) 1245 COS = p5/SQRT(p25 + p25*TAN**2) 1246 SIN = COS*TAN 1247 END IF 1248 1249 ! Compute the modified diagonal element of r and 1250 ! the modified element of ((q transpose)*b,0). 1251 1252 r(k,k) = COS*r(k,k) + SIN*sdiag(k) 1253 temp = COS*wa(k) + SIN*qtbpj 1254 qtbpj = -SIN*wa(k) + COS*qtbpj 1255 wa(k) = temp 1256 1257 ! Accumulate the tranformation in the row of s. 1258 1259 kp1 = k + 1 1260 DO i = kp1, n 1261 temp = COS*r(i,k) + SIN*sdiag(i) 1262 sdiag(i) = -SIN*r(i,k) + COS*sdiag(i) 1263 r(i,k) = temp 1264 END DO 1265 END DO 1266 1267 ! Store the diagonal element of s and restore 1268 ! the corresponding diagonal element of r. 1269 1270 sdiag(j) = r(j,j) 1271 r(j,j) = x(j) 1272 END DO 1273 1274 ! Solve the triangular system for z. If the system is singular, 1275 ! then obtain a least squares solution. 1276 1277 nsing = n 1278 DO j = 1, n 1279 IF (sdiag(j) == zero .AND. nsing == n) nsing = j - 1 1280 IF (nsing < n) wa(j) = zero 1281 END DO 1282 1283 DO k = 1, nsing 1284 j = nsing - k + 1 1285 sum = DOT_PRODUCT( r(j+1:nsing,j), wa(j+1:nsing) ) 1286 wa(j) = (wa(j) - sum)/sdiag(j) 1287 END DO 1288 1289 ! Permute the components of z back to components of x. 1290 1291 DO j = 1, n 1292 l = ipvt(j) 1293 x(l) = wa(j) 1294 END DO 1295 RETURN 1296 1297 ! last card of subroutine qrsolv. 1298 1299 END SUBROUTINE qrsolv