share/minpack/fortran/qrfac.f

   1       subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
   2       integer m,n,lda,lipvt
   3       integer ipvt(lipvt)
   4       logical pivot
   5       double precision a(lda,n),rdiag(n),acnorm(n),wa(n)
   6 c     **********
   7 c
   8 c     subroutine qrfac
   9 c
  10 c     this subroutine uses householder transformations with column
  11 c     pivoting (optional) to compute a qr factorization of the
  12 c     m by n matrix a. that is, qrfac determines an orthogonal
  13 c     matrix q, a permutation matrix p, and an upper trapezoidal
  14 c     matrix r with diagonal elements of nonincreasing magnitude,
  15 c     such that a*p = q*r. the householder transformation for
  16 c     column k, k = 1,2,...,min(m,n), is of the form
  17 c
  18 c                           t
  19 c           i - (1/u(k))*u*u
  20 c
  21 c     where u has zeros in the first k-1 positions. the form of
  22 c     this transformation and the method of pivoting first
  23 c     appeared in the corresponding linpack subroutine.
  24 c
  25 c     the subroutine statement is
  26 c
  27 c       subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
  28 c
  29 c     where
  30 c
  31 c       m is a positive integer input variable set to the number
  32 c         of rows of a.
  33 c
  34 c       n is a positive integer input variable set to the number
  35 c         of columns of a.
  36 c
  37 c       a is an m by n array. on input a contains the matrix for
  38 c         which the qr factorization is to be computed. on output
  39 c         the strict upper trapezoidal part of a contains the strict
  40 c         upper trapezoidal part of r, and the lower trapezoidal
  41 c         part of a contains a factored form of q (the non-trivial
  42 c         elements of the u vectors described above).
  43 c
  44 c       lda is a positive integer input variable not less than m
  45 c         which specifies the leading dimension of the array a.
  46 c
  47 c       pivot is a logical input variable. if pivot is set true,
  48 c         then column pivoting is enforced. if pivot is set false,
  49 c         then no column pivoting is done.
  50 c
  51 c       ipvt is an integer output array of length lipvt. ipvt
  52 c         defines the permutation matrix p such that a*p = q*r.
  53 c         column j of p is column ipvt(j) of the identity matrix.
  54 c         if pivot is false, ipvt is not referenced.
  55 c
  56 c       lipvt is a positive integer input variable. if pivot is false,
  57 c         then lipvt may be as small as 1. if pivot is true, then
  58 c         lipvt must be at least n.
  59 c
  60 c       rdiag is an output array of length n which contains the
  61 c         diagonal elements of r.
  62 c
  63 c       acnorm is an output array of length n which contains the
  64 c         norms of the corresponding columns of the input matrix a.
  65 c         if this information is not needed, then acnorm can coincide
  66 c         with rdiag.
  67 c
  68 c       wa is a work array of length n. if pivot is false, then wa
  69 c         can coincide with rdiag.
  70 c
  71 c     subprograms called
  72 c
  73 c       minpack-supplied ... dpmpar,enorm
  74 c
  75 c       fortran-supplied ... dmax1,dsqrt,min0
  76 c
  77 c     argonne national laboratory. minpack project. march 1980.
  78 c     burton s. garbow, kenneth e. hillstrom, jorge j. more
  79 c
  80 c     **********
  81       integer i,j,jp1,k,kmax,minmn
  82       double precision ajnorm,epsmch,one,p05,sum,temp,zero
  83       double precision dpmpar,enorm
  84       data one,p05,zero /1.0d0,5.0d-2,0.0d0/
  85 c
  86 c     epsmch is the machine precision.
  87 c
  88       epsmch = dpmpar(1)
  89 c
  90 c     compute the initial column norms and initialize several arrays.
  91 c
  92       do 10 j = 1, n
  93          acnorm(j) = enorm(m,a(1,j))
  94          rdiag(j) = acnorm(j)
  95          wa(j) = rdiag(j)
  96          if (pivot) ipvt(j) = j
  97    10    continue
  98 c
  99 c     reduce a to r with householder transformations.
 100 c
 101       minmn = min0(m,n)
 102       do 110 j = 1, minmn
 103          if (.not.pivot) go to 40
 104 c
 105 c        bring the column of largest norm into the pivot position.
 106 c
 107          kmax = j
 108          do 20 k = j, n
 109             if (rdiag(k) .gt. rdiag(kmax)) kmax = k
 110    20       continue
 111          if (kmax .eq. j) go to 40
 112          do 30 i = 1, m
 113             temp = a(i,j)
 114             a(i,j) = a(i,kmax)
 115             a(i,kmax) = temp
 116    30       continue
 117          rdiag(kmax) = rdiag(j)
 118          wa(kmax) = wa(j)
 119          k = ipvt(j)
 120          ipvt(j) = ipvt(kmax)
 121          ipvt(kmax) = k
 122    40    continue
 123 c
 124 c        compute the householder transformation to reduce the
 125 c        j-th column of a to a multiple of the j-th unit vector.
 126 c
 127          ajnorm = enorm(m-j+1,a(j,j))
 128          if (ajnorm .eq. zero) go to 100
 129          if (a(j,j) .lt. zero) ajnorm = -ajnorm
 130          do 50 i = j, m
 131             a(i,j) = a(i,j)/ajnorm
 132    50       continue
 133          a(j,j) = a(j,j) + one
 134 c
 135 c        apply the transformation to the remaining columns
 136 c        and update the norms.
 137 c
 138          jp1 = j + 1
 139          if (n .lt. jp1) go to 100
 140          do 90 k = jp1, n
 141             sum = zero
 142             do 60 i = j, m
 143                sum = sum + a(i,j)*a(i,k)
 144    60          continue
 145             temp = sum/a(j,j)
 146             do 70 i = j, m
 147                a(i,k) = a(i,k) - temp*a(i,j)
 148    70          continue
 149             if (.not.pivot .or. rdiag(k) .eq. zero) go to 80
 150             temp = a(j,k)/rdiag(k)
 151             rdiag(k) = rdiag(k)*dsqrt(dmax1(zero,one-temp**2))
 152             if (p05*(rdiag(k)/wa(k))**2 .gt. epsmch) go to 80
 153             rdiag(k) = enorm(m-j,a(jp1,k))
 154             wa(k) = rdiag(k)
 155    80       continue
 156    90       continue
 157   100    continue
 158          rdiag(j) = -ajnorm
 159   110    continue
 160       return
 161 c
 162 c     last card of subroutine qrfac.
 163 c
 164       end