Repatch (from linux) of OSX IRAF

1 files changed, 440 insertions, 0 deletions

diff --git a/math/minpack/hybrj.f b/math/minpack/hybrj.f
new file mode 100644
index 00000000..3070dad3
--- /dev/null
+++ b/math/minpack/hybrj.f
@@ -0,0 +1,440 @@
+      subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode,
+     *                 factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2,
+     *                 wa3,wa4)
+      integer n,ldfjac,maxfev,mode,nprint,info,nfev,njev,lr
+      double precision xtol,factor
+      double precision x(n),fvec(n),fjac(ldfjac,n),diag(n),r(lr),
+     *                 qtf(n),wa1(n),wa2(n),wa3(n),wa4(n)
+c     **********
+c
+c     subroutine hybrj
+c
+c     the purpose of hybrj is to find a zero of a system of
+c     n nonlinear functions in n variables by a modification
+c     of the powell hybrid method. the user must provide a
+c     subroutine which calculates the functions and the jacobian.
+c
+c     the subroutine statement is
+c
+c       subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,
+c                        mode,factor,nprint,info,nfev,njev,r,lr,qtf,
+c                        wa1,wa2,wa3,wa4)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions and the jacobian. fcn must
+c         be declared in an external statement in the user
+c         calling program, and should be written as follows.
+c
+c         subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
+c         integer n,ldfjac,iflag
+c         double precision x(n),fvec(n),fjac(ldfjac,n)
+c         ----------
+c         if iflag = 1 calculate the functions at x and
+c         return this vector in fvec. do not alter fjac.
+c         if iflag = 2 calculate the jacobian at x and
+c         return this matrix in fjac. do not alter fvec.
+c         ---------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of hybrj.
+c         in this case set iflag to a negative integer.
+c
+c       n is a positive integer input variable set to the number
+c         of functions and variables.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length n which contains
+c         the functions evaluated at the output x.
+c
+c       fjac is an output n by n array which contains the
+c         orthogonal matrix q produced by the qr factorization
+c         of the final approximate jacobian.
+c
+c       ldfjac is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array fjac.
+c
+c       xtol is a nonnegative input variable. termination
+c         occurs when the relative error between two consecutive
+c         iterates is at most xtol.
+c
+c       maxfev is a positive integer input variable. termination
+c         occurs when the number of calls to fcn with iflag = 1
+c         has reached maxfev.
+c
+c       diag is an array of length n. if mode = 1 (see
+c         below), diag is internally set. if mode = 2, diag
+c         must contain positive entries that serve as
+c         multiplicative scale factors for the variables.
+c
+c       mode is an integer input variable. if mode = 1, the
+c         variables will be scaled internally. if mode = 2,
+c         the scaling is specified by the input diag. other
+c         values of mode are equivalent to mode = 1.
+c
+c       factor is a positive input variable used in determining the
+c         initial step bound. this bound is set to the product of
+c         factor and the euclidean norm of diag*x if nonzero, or else
+c         to factor itself. in most cases factor should lie in the
+c         interval (.1,100.). 100. is a generally recommended value.
+c
+c       nprint is an integer input variable that enables controlled
+c         printing of iterates if it is positive. in this case,
+c         fcn is called with iflag = 0 at the beginning of the first
+c         iteration and every nprint iterations thereafter and
+c         immediately prior to return, with x and fvec available
+c         for printing. fvec and fjac should not be altered.
+c         if nprint is not positive, no special calls of fcn
+c         with iflag = 0 are made.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0   improper input parameters.
+c
+c         info = 1   relative error between two consecutive iterates
+c                    is at most xtol.
+c
+c         info = 2   number of calls to fcn with iflag = 1 has
+c                    reached maxfev.
+c
+c         info = 3   xtol is too small. no further improvement in
+c                    the approximate solution x is possible.
+c
+c         info = 4   iteration is not making good progress, as
+c                    measured by the improvement from the last
+c                    five jacobian evaluations.
+c
+c         info = 5   iteration is not making good progress, as
+c                    measured by the improvement from the last
+c                    ten iterations.
+c
+c       nfev is an integer output variable set to the number of
+c         calls to fcn with iflag = 1.
+c
+c       njev is an integer output variable set to the number of
+c         calls to fcn with iflag = 2.
+c
+c       r is an output array of length lr which contains the
+c         upper triangular matrix produced by the qr factorization
+c         of the final approximate jacobian, stored rowwise.
+c
+c       lr is a positive integer input variable not less than
+c         (n*(n+1))/2.
+c
+c       qtf is an output array of length n which contains
+c         the vector (q transpose)*fvec.
+c
+c       wa1, wa2, wa3, and wa4 are work arrays of length n.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... dogleg,dpmpar,enorm,
+c                            qform,qrfac,r1mpyq,r1updt
+c
+c       fortran-supplied ... dabs,dmax1,dmin1,mod
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer i,iflag,iter,j,jm1,l,ncfail,ncsuc,nslow1,nslow2
+      integer iwa(1)
+      logical jeval,sing
+      double precision actred,delta,epsmch,fnorm,fnorm1,one,pnorm,
+     *                 prered,p1,p5,p001,p0001,ratio,sum,temp,xnorm,
+     *                 zero
+      double precision dpmpar,enorm
+      data one,p1,p5,p001,p0001,zero
+     *     /1.0d0,1.0d-1,5.0d-1,1.0d-3,1.0d-4,0.0d0/
+c
+c     epsmch is the machine precision.
+c
+      epsmch = dpmpar(1)
+c
+      info = 0
+      iflag = 0
+      nfev = 0
+      njev = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. ldfjac .lt. n .or. xtol .lt. zero
+     *    .or. maxfev .le. 0 .or. factor .le. zero
+     *    .or. lr .lt. (n*(n + 1))/2) go to 300
+      if (mode .ne. 2) go to 20
+      do 10 j = 1, n
+         if (diag(j) .le. zero) go to 300
+   10    continue
+   20 continue
+c
+c     evaluate the function at the starting point
+c     and calculate its norm.
+c
+      iflag = 1
+      call fcn(n,x,fvec,fjac,ldfjac,iflag)
+      nfev = 1
+      if (iflag .lt. 0) go to 300
+      fnorm = enorm(n,fvec)
+c
+c     initialize iteration counter and monitors.
+c
+      iter = 1
+      ncsuc = 0
+      ncfail = 0
+      nslow1 = 0
+      nslow2 = 0
+c
+c     beginning of the outer loop.
+c
+   30 continue
+         jeval = .true.
+c
+c        calculate the jacobian matrix.
+c
+         iflag = 2
+         call fcn(n,x,fvec,fjac,ldfjac,iflag)
+         njev = njev + 1
+         if (iflag .lt. 0) go to 300
+c
+c        compute the qr factorization of the jacobian.
+c
+         call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
+c
+c        on the first iteration and if mode is 1, scale according
+c        to the norms of the columns of the initial jacobian.
+c
+         if (iter .ne. 1) go to 70
+         if (mode .eq. 2) go to 50
+         do 40 j = 1, n
+            diag(j) = wa2(j)
+            if (wa2(j) .eq. zero) diag(j) = one
+   40       continue
+   50    continue
+c
+c        on the first iteration, calculate the norm of the scaled x
+c        and initialize the step bound delta.
+c
+         do 60 j = 1, n
+            wa3(j) = diag(j)*x(j)
+   60       continue
+         xnorm = enorm(n,wa3)
+         delta = factor*xnorm
+         if (delta .eq. zero) delta = factor
+   70    continue
+c
+c        form (q transpose)*fvec and store in qtf.
+c
+         do 80 i = 1, n
+            qtf(i) = fvec(i)
+   80       continue
+         do 120 j = 1, n
+            if (fjac(j,j) .eq. zero) go to 110
+            sum = zero
+            do 90 i = j, n
+               sum = sum + fjac(i,j)*qtf(i)
+   90          continue
+            temp = -sum/fjac(j,j)
+            do 100 i = j, n
+               qtf(i) = qtf(i) + fjac(i,j)*temp
+  100          continue
+  110       continue
+  120       continue
+c
+c        copy the triangular factor of the qr factorization into r.
+c
+         sing = .false.
+         do 150 j = 1, n
+            l = j
+            jm1 = j - 1
+            if (jm1 .lt. 1) go to 140
+            do 130 i = 1, jm1
+               r(l) = fjac(i,j)
+               l = l + n - i
+  130          continue
+  140       continue
+            r(l) = wa1(j)
+            if (wa1(j) .eq. zero) sing = .true.
+  150       continue
+c
+c        accumulate the orthogonal factor in fjac.
+c
+         call qform(n,n,fjac,ldfjac,wa1)
+c
+c        rescale if necessary.
+c
+         if (mode .eq. 2) go to 170
+         do 160 j = 1, n
+            diag(j) = dmax1(diag(j),wa2(j))
+  160       continue
+  170    continue
+c
+c        beginning of the inner loop.
+c
+  180    continue
+c
+c           if requested, call fcn to enable printing of iterates.
+c
+            if (nprint .le. 0) go to 190
+            iflag = 0
+            if (mod(iter-1,nprint) .eq. 0)
+     *         call fcn(n,x,fvec,fjac,ldfjac,iflag)
+            if (iflag .lt. 0) go to 300
+  190       continue
+c
+c           determine the direction p.
+c
+            call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
+c
+c           store the direction p and x + p. calculate the norm of p.
+c
+            do 200 j = 1, n
+               wa1(j) = -wa1(j)
+               wa2(j) = x(j) + wa1(j)
+               wa3(j) = diag(j)*wa1(j)
+  200          continue
+            pnorm = enorm(n,wa3)
+c
+c           on the first iteration, adjust the initial step bound.
+c
+            if (iter .eq. 1) delta = dmin1(delta,pnorm)
+c
+c           evaluate the function at x + p and calculate its norm.
+c
+            iflag = 1
+            call fcn(n,wa2,wa4,fjac,ldfjac,iflag)
+            nfev = nfev + 1
+            if (iflag .lt. 0) go to 300
+            fnorm1 = enorm(n,wa4)
+c
+c           compute the scaled actual reduction.
+c
+            actred = -one
+            if (fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
+c
+c           compute the scaled predicted reduction.
+c
+            l = 1
+            do 220 i = 1, n
+               sum = zero
+               do 210 j = i, n
+                  sum = sum + r(l)*wa1(j)
+                  l = l + 1
+  210             continue
+               wa3(i) = qtf(i) + sum
+  220          continue
+            temp = enorm(n,wa3)
+            prered = zero
+            if (temp .lt. fnorm) prered = one - (temp/fnorm)**2
+c
+c           compute the ratio of the actual to the predicted
+c           reduction.
+c
+            ratio = zero
+            if (prered .gt. zero) ratio = actred/prered
+c
+c           update the step bound.
+c
+            if (ratio .ge. p1) go to 230
+               ncsuc = 0
+               ncfail = ncfail + 1
+               delta = p5*delta
+               go to 240
+  230       continue
+               ncfail = 0
+               ncsuc = ncsuc + 1
+               if (ratio .ge. p5 .or. ncsuc .gt. 1)
+     *            delta = dmax1(delta,pnorm/p5)
+               if (dabs(ratio-one) .le. p1) delta = pnorm/p5
+  240       continue
+c
+c           test for successful iteration.
+c
+            if (ratio .lt. p0001) go to 260
+c
+c           successful iteration. update x, fvec, and their norms.
+c
+            do 250 j = 1, n
+               x(j) = wa2(j)
+               wa2(j) = diag(j)*x(j)
+               fvec(j) = wa4(j)
+  250          continue
+            xnorm = enorm(n,wa2)
+            fnorm = fnorm1
+            iter = iter + 1
+  260       continue
+c
+c           determine the progress of the iteration.
+c
+            nslow1 = nslow1 + 1
+            if (actred .ge. p001) nslow1 = 0
+            if (jeval) nslow2 = nslow2 + 1
+            if (actred .ge. p1) nslow2 = 0
+c
+c           test for convergence.
+c
+            if (delta .le. xtol*xnorm .or. fnorm .eq. zero) info = 1
+            if (info .ne. 0) go to 300
+c
+c           tests for termination and stringent tolerances.
+c
+            if (nfev .ge. maxfev) info = 2
+            if (p1*dmax1(p1*delta,pnorm) .le. epsmch*xnorm) info = 3
+            if (nslow2 .eq. 5) info = 4
+            if (nslow1 .eq. 10) info = 5
+            if (info .ne. 0) go to 300
+c
+c           criterion for recalculating jacobian.
+c
+            if (ncfail .eq. 2) go to 290
+c
+c           calculate the rank one modification to the jacobian
+c           and update qtf if necessary.
+c
+            do 280 j = 1, n
+               sum = zero
+               do 270 i = 1, n
+                  sum = sum + fjac(i,j)*wa4(i)
+  270             continue
+               wa2(j) = (sum - wa3(j))/pnorm
+               wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
+               if (ratio .ge. p0001) qtf(j) = sum
+  280          continue
+c
+c           compute the qr factorization of the updated jacobian.
+c
+            call r1updt(n,n,r,lr,wa1,wa2,wa3,sing)
+            call r1mpyq(n,n,fjac,ldfjac,wa2,wa3)
+            call r1mpyq(1,n,qtf,1,wa2,wa3)
+c
+c           end of the inner loop.
+c
+            jeval = .false.
+            go to 180
+  290    continue
+c
+c        end of the outer loop.
+c
+         go to 30
+  300 continue
+c
+c     termination, either normal or user imposed.
+c
+      if (iflag .lt. 0) info = iflag
+      iflag = 0
+      if (nprint .gt. 0) call fcn(n,x,fvec,fjac,ldfjac,iflag)
+      return
+c
+c     last card of subroutine hybrj.
+c
+      end