1 files changed, 95 insertions, 0 deletions
diff --git a/math/minpack/qform.f b/math/minpack/qform.f
new file mode 100644
index 00000000..087b2478
--- /dev/null
+++ b/math/minpack/qform.f
@@ -0,0 +1,95 @@
+      subroutine qform(m,n,q,ldq,wa)
+      integer m,n,ldq
+      double precision q(ldq,m),wa(m)
+c     **********
+c
+c     subroutine qform
+c
+c     this subroutine proceeds from the computed qr factorization of
+c     an m by n matrix a to accumulate the m by m orthogonal matrix
+c     q from its factored form.
+c
+c     the subroutine statement is
+c
+c       subroutine qform(m,n,q,ldq,wa)
+c
+c     where
+c
+c       m is a positive integer input variable set to the number
+c         of rows of a and the order of q.
+c
+c       n is a positive integer input variable set to the number
+c         of columns of a.
+c
+c       q is an m by m array. on input the full lower trapezoid in
+c         the first min(m,n) columns of q contains the factored form.
+c         on output q has been accumulated into a square matrix.
+c
+c       ldq is a positive integer input variable not less than m
+c         which specifies the leading dimension of the array q.
+c
+c       wa is a work array of length m.
+c
+c     subprograms called
+c
+c       fortran-supplied ... min0
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer i,j,jm1,k,l,minmn,np1
+      double precision one,sum,temp,zero
+      data one,zero /1.0d0,0.0d0/
+c
+c     zero out upper triangle of q in the first min(m,n) columns.
+c
+      minmn = min0(m,n)
+      if (minmn .lt. 2) go to 30
+      do 20 j = 2, minmn
+         jm1 = j - 1
+         do 10 i = 1, jm1
+            q(i,j) = zero
+   10       continue
+   20    continue
+   30 continue
+c
+c     initialize remaining columns to those of the identity matrix.
+c
+      np1 = n + 1
+      if (m .lt. np1) go to 60
+      do 50 j = np1, m
+         do 40 i = 1, m
+            q(i,j) = zero
+   40       continue
+         q(j,j) = one
+   50    continue
+   60 continue
+c
+c     accumulate q from its factored form.
+c
+      do 120 l = 1, minmn
+         k = minmn - l + 1
+         do 70 i = k, m
+            wa(i) = q(i,k)
+            q(i,k) = zero
+   70       continue
+         q(k,k) = one
+         if (wa(k) .eq. zero) go to 110
+         do 100 j = k, m
+            sum = zero
+            do 80 i = k, m
+               sum = sum + q(i,j)*wa(i)
+   80          continue
+            temp = sum/wa(k)
+            do 90 i = k, m
+               q(i,j) = q(i,j) - temp*wa(i)
+   90          continue
+  100       continue
+  110    continue
+  120    continue
+      return
+c
+c     last card of subroutine qform.
+c
+      end