Repatch (from linux) of OSX IRAF

1 files changed, 113 insertions, 0 deletions

diff --git a/math/minpack/rwupdt.f b/math/minpack/rwupdt.f
new file mode 100644
index 00000000..05282b55
--- /dev/null
+++ b/math/minpack/rwupdt.f
@@ -0,0 +1,113 @@
+      subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
+      integer n,ldr
+      double precision alpha
+      double precision r(ldr,n),w(n),b(n),cos(n),sin(n)
+c     **********
+c
+c     subroutine rwupdt
+c
+c     given an n by n upper triangular matrix r, this subroutine
+c     computes the qr decomposition of the matrix formed when a row
+c     is added to r. if the row is specified by the vector w, then
+c     rwupdt determines an orthogonal matrix q such that when the
+c     n+1 by n matrix composed of r augmented by w is premultiplied
+c     by (q transpose), the resulting matrix is upper trapezoidal.
+c     the matrix (q transpose) is the product of n transformations
+c
+c           g(n)*g(n-1)* ... *g(1)
+c
+c     where g(i) is a givens rotation in the (i,n+1) plane which
+c     eliminates elements in the (n+1)-st plane. rwupdt also
+c     computes the product (q transpose)*c where c is the
+c     (n+1)-vector (b,alpha). q itself is not accumulated, rather
+c     the information to recover the g rotations is supplied.
+c
+c     the subroutine statement is
+c
+c       subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
+c
+c     where
+c
+c       n is a positive integer input variable set to the order of r.
+c
+c       r is an n by n array. on input the upper triangular part of
+c         r must contain the matrix to be updated. on output r
+c         contains the updated triangular matrix.
+c
+c       ldr is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array r.
+c
+c       w is an input array of length n which must contain the row
+c         vector to be added to r.
+c
+c       b is an array of length n. on input b must contain the
+c         first n elements of the vector c. on output b contains
+c         the first n elements of the vector (q transpose)*c.
+c
+c       alpha is a variable. on input alpha must contain the
+c         (n+1)-st element of the vector c. on output alpha contains
+c         the (n+1)-st element of the vector (q transpose)*c.
+c
+c       cos is an output array of length n which contains the
+c         cosines of the transforming givens rotations.
+c
+c       sin is an output array of length n which contains the
+c         sines of the transforming givens rotations.
+c
+c     subprograms called
+c
+c       fortran-supplied ... dabs,dsqrt
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, dudley v. goetschel, kenneth e. hillstrom,
+c     jorge j. more
+c
+c     **********
+      integer i,j,jm1
+      double precision cotan,one,p5,p25,rowj,tan,temp,zero
+      data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
+c
+      do 60 j = 1, n
+         rowj = w(j)
+         jm1 = j - 1
+c
+c        apply the previous transformations to
+c        r(i,j), i=1,2,...,j-1, and to w(j).
+c
+         if (jm1 .lt. 1) go to 20
+         do 10 i = 1, jm1
+            temp = cos(i)*r(i,j) + sin(i)*rowj
+            rowj = -sin(i)*r(i,j) + cos(i)*rowj
+            r(i,j) = temp
+   10       continue
+   20    continue
+c
+c        determine a givens rotation which eliminates w(j).
+c
+         cos(j) = one
+         sin(j) = zero
+         if (rowj .eq. zero) go to 50
+         if (dabs(r(j,j)) .ge. dabs(rowj)) go to 30
+            cotan = r(j,j)/rowj
+            sin(j) = p5/dsqrt(p25+p25*cotan**2)
+            cos(j) = sin(j)*cotan
+            go to 40
+   30    continue
+            tan = rowj/r(j,j)
+            cos(j) = p5/dsqrt(p25+p25*tan**2)
+            sin(j) = cos(j)*tan
+   40    continue
+c
+c        apply the current transformation to r(j,j), b(j), and alpha.
+c
+         r(j,j) = cos(j)*r(j,j) + sin(j)*rowj
+         temp = cos(j)*b(j) + sin(j)*alpha
+         alpha = -sin(j)*b(j) + cos(j)*alpha
+         b(j) = temp
+   50    continue
+   60    continue
+      return
+c
+c     last card of subroutine rwupdt.
+c
+      end