Added code for solving nonlinear systems.

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@1665 f25d12d1-65f4-0310-ae8a-bbce733d8d8e

c_runtime/dogleg.f

@@ -0,0 +1,177 @@
+      subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
+      integer n,lr
+      double precision delta
+      double precision r(lr),diag(n),qtb(n),x(n),wa1(n),wa2(n)
+c     **********
+c
+c     subroutine dogleg
+c
+c     given an m by n matrix a, an n by n nonsingular diagonal
+c     matrix d, an m-vector b, and a positive number delta, the
+c     problem is to determine the convex combination x of the
+c     gauss-newton and scaled gradient directions that minimizes
+c     (a*x - b) in the least squares sense, subject to the
+c     restriction that the euclidean norm of d*x be at most delta.
+c
+c     this subroutine completes the solution of the problem
+c     if it is provided with the necessary information from the
+c     qr factorization of a. that is, if a = q*r, where q has
+c     orthogonal columns and r is an upper triangular matrix,
+c     then dogleg expects the full upper triangle of r and
+c     the first n components of (q transpose)*b.
+c
+c     the subroutine statement is
+c
+c       subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
+c
+c     where
+c
+c       n is a positive integer input variable set to the order of r.
+c
+c       r is an input array of length lr which must contain the upper
+c         triangular matrix r stored by rows.
+c
+c       lr is a positive integer input variable not less than
+c         (n*(n+1))/2.
+c
+c       diag is an input array of length n which must contain the
+c         diagonal elements of the matrix d.
+c
+c       qtb is an input array of length n which must contain the first
+c         n elements of the vector (q transpose)*b.
+c
+c       delta is a positive input variable which specifies an upper
+c         bound on the euclidean norm of d*x.
+c
+c       x is an output array of length n which contains the desired
+c         convex combination of the gauss-newton direction and the
+c         scaled gradient direction.
+c
+c       wa1 and wa2 are work arrays of length n.
+c
+c     subprograms called
+c
+c       minpack-supplied ... dpmpar,enorm
+c
+c       fortran-supplied ... dabs,dmax1,dmin1,dsqrt
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer i,j,jj,jp1,k,l
+      double precision alpha,bnorm,epsmch,gnorm,one,qnorm,sgnorm,sum,
+     *                 temp,zero
+      double precision dpmpar,enorm
+      data one,zero /1.0d0,0.0d0/
+c
+c     epsmch is the machine precision.
+c
+      epsmch = dpmpar(1)
+c
+c     first, calculate the gauss-newton direction.
+c
+      jj = (n*(n + 1))/2 + 1
+      do 50 k = 1, n
+         j = n - k + 1
+         jp1 = j + 1
+         jj = jj - k
+         l = jj + 1
+         sum = zero
+         if (n .lt. jp1) go to 20
+         do 10 i = jp1, n
+            sum = sum + r(l)*x(i)
+            l = l + 1
+   10       continue
+   20    continue
+         temp = r(jj)
+         if (temp .ne. zero) go to 40
+         l = j
+         do 30 i = 1, j
+            temp = dmax1(temp,dabs(r(l)))
+            l = l + n - i
+   30       continue
+         temp = epsmch*temp
+         if (temp .eq. zero) temp = epsmch
+   40    continue
+         x(j) = (qtb(j) - sum)/temp
+   50    continue
+c
+c     test whether the gauss-newton direction is acceptable.
+c
+      do 60 j = 1, n
+         wa1(j) = zero
+         wa2(j) = diag(j)*x(j)
+   60    continue
+      qnorm = enorm(n,wa2)
+      if (qnorm .le. delta) go to 140
+c
+c     the gauss-newton direction is not acceptable.
+c     next, calculate the scaled gradient direction.
+c
+      l = 1
+      do 80 j = 1, n
+         temp = qtb(j)
+         do 70 i = j, n
+            wa1(i) = wa1(i) + r(l)*temp
+            l = l + 1
+   70       continue
+         wa1(j) = wa1(j)/diag(j)
+   80    continue
+c
+c     calculate the norm of the scaled gradient and test for
+c     the special case in which the scaled gradient is zero.
+c
+      gnorm = enorm(n,wa1)
+      sgnorm = zero
+      alpha = delta/qnorm
+      if (gnorm .eq. zero) go to 120
+c
+c     calculate the point along the scaled gradient
+c     at which the quadratic is minimized.
+c
+      do 90 j = 1, n
+         wa1(j) = (wa1(j)/gnorm)/diag(j)
+   90    continue
+      l = 1
+      do 110 j = 1, n
+         sum = zero
+         do 100 i = j, n
+            sum = sum + r(l)*wa1(i)
+            l = l + 1
+  100       continue
+         wa2(j) = sum
+  110    continue
+      temp = enorm(n,wa2)
+      sgnorm = (gnorm/temp)/temp
+c
+c     test whether the scaled gradient direction is acceptable.
+c
+      alpha = zero
+      if (sgnorm .ge. delta) go to 120
+c
+c     the scaled gradient direction is not acceptable.
+c     finally, calculate the point along the dogleg
+c     at which the quadratic is minimized.
+c
+      bnorm = enorm(n,qtb)
+      temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
+      temp = temp - (delta/qnorm)*(sgnorm/delta)**2
+     *       + dsqrt((temp-(delta/qnorm))**2
+     *               +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2))
+      alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp
+  120 continue
+c
+c     form appropriate convex combination of the gauss-newton
+c     direction and the scaled gradient direction.
+c
+      temp = (one - alpha)*dmin1(sgnorm,delta)
+      do 130 j = 1, n
+         x(j) = temp*wa1(j) + alpha*x(j)
+  130    continue
+  140 continue
+      return
+c
+c     last card of subroutine dogleg.
+c
+      end

c_runtime/dpmpar.f

@@ -0,0 +1,177 @@
+      double precision function dpmpar(i)
+      integer i
+c     **********
+c
+c     Function dpmpar
+c
+c     This function provides double precision machine parameters
+c     when the appropriate set of data statements is activated (by
+c     removing the c from column 1) and all other data statements are
+c     rendered inactive. Most of the parameter values were obtained
+c     from the corresponding Bell Laboratories Port Library function.
+c
+c     The function statement is
+c
+c       double precision function dpmpar(i)
+c
+c     where
+c
+c       i is an integer input variable set to 1, 2, or 3 which
+c         selects the desired machine parameter. If the machine has
+c         t base b digits and its smallest and largest exponents are
+c         emin and emax, respectively, then these parameters are
+c
+c         dpmpar(1) = b**(1 - t), the machine precision,
+c
+c         dpmpar(2) = b**(emin - 1), the smallest magnitude,
+c
+c         dpmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude.
+c
+c     Argonne National Laboratory. MINPACK Project. November 1996.
+c     Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More'
+c
+c     **********
+      integer mcheps(4)
+      integer minmag(4)
+      integer maxmag(4)
+      double precision dmach(3)
+      equivalence (dmach(1),mcheps(1))
+      equivalence (dmach(2),minmag(1))
+      equivalence (dmach(3),maxmag(1))
+c
+c     Machine constants for the IBM 360/370 series,
+c     the Amdahl 470/V6, the ICL 2900, the Itel AS/6,
+c     the Xerox Sigma 5/7/9 and the Sel systems 85/86.
+c
+c     data mcheps(1),mcheps(2) / z34100000, z00000000 /
+c     data minmag(1),minmag(2) / z00100000, z00000000 /
+c     data maxmag(1),maxmag(2) / z7fffffff, zffffffff /
+c
+c     Machine constants for the Honeywell 600/6000 series.
+c
+c     data mcheps(1),mcheps(2) / o606400000000, o000000000000 /
+c     data minmag(1),minmag(2) / o402400000000, o000000000000 /
+c     data maxmag(1),maxmag(2) / o376777777777, o777777777777 /
+c
+c     Machine constants for the CDC 6000/7000 series.
+c
+c     data mcheps(1) / 15614000000000000000b /
+c     data mcheps(2) / 15010000000000000000b /
+c
+c     data minmag(1) / 00604000000000000000b /
+c     data minmag(2) / 00000000000000000000b /
+c
+c     data maxmag(1) / 37767777777777777777b /
+c     data maxmag(2) / 37167777777777777777b /
+c
+c     Machine constants for the PDP-10 (KA processor).
+c
+c     data mcheps(1),mcheps(2) / "114400000000, "000000000000 /
+c     data minmag(1),minmag(2) / "033400000000, "000000000000 /
+c     data maxmag(1),maxmag(2) / "377777777777, "344777777777 /
+c
+c     Machine constants for the PDP-10 (KI processor).
+c
+c     data mcheps(1),mcheps(2) / "104400000000, "000000000000 /
+c     data minmag(1),minmag(2) / "000400000000, "000000000000 /
+c     data maxmag(1),maxmag(2) / "377777777777, "377777777777 /
+c
+c     Machine constants for the PDP-11. 
+c
+c     data mcheps(1),mcheps(2) /   9472,      0 /
+c     data mcheps(3),mcheps(4) /      0,      0 /
+c
+c     data minmag(1),minmag(2) /    128,      0 /
+c     data minmag(3),minmag(4) /      0,      0 /
+c
+c     data maxmag(1),maxmag(2) /  32767,     -1 /
+c     data maxmag(3),maxmag(4) /     -1,     -1 /
+c
+c     Machine constants for the Burroughs 6700/7700 systems.
+c
+c     data mcheps(1) / o1451000000000000 /
+c     data mcheps(2) / o0000000000000000 /
+c
+c     data minmag(1) / o1771000000000000 /
+c     data minmag(2) / o7770000000000000 /
+c
+c     data maxmag(1) / o0777777777777777 /
+c     data maxmag(2) / o7777777777777777 /
+c
+c     Machine constants for the Burroughs 5700 system.
+c
+c     data mcheps(1) / o1451000000000000 /
+c     data mcheps(2) / o0000000000000000 /
+c
+c     data minmag(1) / o1771000000000000 /
+c     data minmag(2) / o0000000000000000 /
+c
+c     data maxmag(1) / o0777777777777777 /
+c     data maxmag(2) / o0007777777777777 /
+c
+c     Machine constants for the Burroughs 1700 system.
+c
+c     data mcheps(1) / zcc6800000 /
+c     data mcheps(2) / z000000000 /
+c
+c     data minmag(1) / zc00800000 /
+c     data minmag(2) / z000000000 /
+c
+c     data maxmag(1) / zdffffffff /
+c     data maxmag(2) / zfffffffff /
+c
+c     Machine constants for the Univac 1100 series.
+c
+c     data mcheps(1),mcheps(2) / o170640000000, o000000000000 /
+c     data minmag(1),minmag(2) / o000040000000, o000000000000 /
+c     data maxmag(1),maxmag(2) / o377777777777, o777777777777 /
+c
+c     Machine constants for the Data General Eclipse S/200.
+c
+c     Note - it may be appropriate to include the following card -
+c     static dmach(3)
+c
+c     data minmag/20k,3*0/,maxmag/77777k,3*177777k/
+c     data mcheps/32020k,3*0/
+c
+c     Machine constants for the Harris 220.
+c
+c     data mcheps(1),mcheps(2) / '20000000, '00000334 /
+c     data minmag(1),minmag(2) / '20000000, '00000201 /
+c     data maxmag(1),maxmag(2) / '37777777, '37777577 /
+c
+c     Machine constants for the Cray-1.
+c
+c     data mcheps(1) / 0376424000000000000000b /
+c     data mcheps(2) / 0000000000000000000000b /
+c
+c     data minmag(1) / 0200034000000000000000b /
+c     data minmag(2) / 0000000000000000000000b /
+c
+c     data maxmag(1) / 0577777777777777777777b /
+c     data maxmag(2) / 0000007777777777777776b /
+c
+c     Machine constants for the Prime 400.
+c
+c     data mcheps(1),mcheps(2) / :10000000000, :00000000123 /
+c     data minmag(1),minmag(2) / :10000000000, :00000100000 /
+c     data maxmag(1),maxmag(2) / :17777777777, :37777677776 /
+c
+c     Machine constants for the VAX-11.
+c
+c     data mcheps(1),mcheps(2) /   9472,  0 /
+c     data minmag(1),minmag(2) /    128,  0 /
+c     data maxmag(1),maxmag(2) / -32769, -1 /
+c
+c     Machine constants for IEEE machines.
+c
+      data dmach(1) /2.22044604926d-16/
+      data dmach(2) /2.22507385852d-308/
+      data dmach(3) /1.79769313485d+308/
+c
+      dpmpar = dmach(i)
+      return
+c
+c     Last card of function dpmpar.
+c
+      end