share/minpack/fortran/qrsolv.f

   1       subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
   2       integer n,ldr
   3       integer ipvt(n)
   4       double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa(n)
   5 c     **********
   6 c
   7 c     subroutine qrsolv
   8 c
   9 c     given an m by n matrix a, an n by n diagonal matrix d,
  10 c     and an m-vector b, the problem is to determine an x which
  11 c     solves the system
  12 c
  13 c           a*x = b ,     d*x = 0 ,
  14 c
  15 c     in the least squares sense.
  16 c
  17 c     this subroutine completes the solution of the problem
  18 c     if it is provided with the necessary information from the
  19 c     qr factorization, with column pivoting, of a. that is, if
  20 c     a*p = q*r, where p is a permutation matrix, q has orthogonal
  21 c     columns, and r is an upper triangular matrix with diagonal
  22 c     elements of nonincreasing magnitude, then qrsolv expects
  23 c     the full upper triangle of r, the permutation matrix p,
  24 c     and the first n components of (q transpose)*b. the system
  25 c     a*x = b, d*x = 0, is then equivalent to
  26 c
  27 c                  t       t
  28 c           r*z = q *b ,  p *d*p*z = 0 ,
  29 c
  30 c     where x = p*z. if this system does not have full rank,
  31 c     then a least squares solution is obtained. on output qrsolv
  32 c     also provides an upper triangular matrix s such that
  33 c
  34 c            t   t               t
  35 c           p *(a *a + d*d)*p = s *s .
  36 c
  37 c     s is computed within qrsolv and may be of separate interest.
  38 c
  39 c     the subroutine statement is
  40 c
  41 c       subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
  42 c
  43 c     where
  44 c
  45 c       n is a positive integer input variable set to the order of r.
  46 c
  47 c       r is an n by n array. on input the full upper triangle
  48 c         must contain the full upper triangle of the matrix r.
  49 c         on output the full upper triangle is unaltered, and the
  50 c         strict lower triangle contains the strict upper triangle
  51 c         (transposed) of the upper triangular matrix s.
  52 c
  53 c       ldr is a positive integer input variable not less than n
  54 c         which specifies the leading dimension of the array r.
  55 c
  56 c       ipvt is an integer input array of length n which defines the
  57 c         permutation matrix p such that a*p = q*r. column j of p
  58 c         is column ipvt(j) of the identity matrix.
  59 c
  60 c       diag is an input array of length n which must contain the
  61 c         diagonal elements of the matrix d.
  62 c
  63 c       qtb is an input array of length n which must contain the first
  64 c         n elements of the vector (q transpose)*b.
  65 c
  66 c       x is an output array of length n which contains the least
  67 c         squares solution of the system a*x = b, d*x = 0.
  68 c
  69 c       sdiag is an output array of length n which contains the
  70 c         diagonal elements of the upper triangular matrix s.
  71 c
  72 c       wa is a work array of length n.
  73 c
  74 c     subprograms called
  75 c
  76 c       fortran-supplied ... dabs,dsqrt
  77 c
  78 c     argonne national laboratory. minpack project. march 1980.
  79 c     burton s. garbow, kenneth e. hillstrom, jorge j. more
  80 c
  81 c     **********
  82       integer i,j,jp1,k,kp1,l,nsing
  83       double precision cos,cotan,p5,p25,qtbpj,sin,sum,tan,temp,zero
  84       data p5,p25,zero /5.0d-1,2.5d-1,0.0d0/
  85 c
  86 c     copy r and (q transpose)*b to preserve input and initialize s.
  87 c     in particular, save the diagonal elements of r in x.
  88 c
  89       do 20 j = 1, n
  90          do 10 i = j, n
  91             r(i,j) = r(j,i)
  92    10       continue
  93          x(j) = r(j,j)
  94          wa(j) = qtb(j)
  95    20    continue
  96 c
  97 c     eliminate the diagonal matrix d using a givens rotation.
  98 c
  99       do 100 j = 1, n
 100 c
 101 c        prepare the row of d to be eliminated, locating the
 102 c        diagonal element using p from the qr factorization.
 103 c
 104          l = ipvt(j)
 105          if (diag(l) .eq. zero) go to 90
 106          do 30 k = j, n
 107             sdiag(k) = zero
 108    30       continue
 109          sdiag(j) = diag(l)
 110 c
 111 c        the transformations to eliminate the row of d
 112 c        modify only a single element of (q transpose)*b
 113 c        beyond the first n, which is initially zero.
 114 c
 115          qtbpj = zero
 116          do 80 k = j, n
 117 c
 118 c           determine a givens rotation which eliminates the
 119 c           appropriate element in the current row of d.
 120 c
 121             if (sdiag(k) .eq. zero) go to 70
 122             if (dabs(r(k,k)) .ge. dabs(sdiag(k))) go to 40
 123                cotan = r(k,k)/sdiag(k)
 124                sin = p5/dsqrt(p25+p25*cotan**2)
 125                cos = sin*cotan
 126                go to 50
 127    40       continue
 128                tan = sdiag(k)/r(k,k)
 129                cos = p5/dsqrt(p25+p25*tan**2)
 130                sin = cos*tan
 131    50       continue
 132 c
 133 c           compute the modified diagonal element of r and
 134 c           the modified element of ((q transpose)*b,0).
 135 c
 136             r(k,k) = cos*r(k,k) + sin*sdiag(k)
 137             temp = cos*wa(k) + sin*qtbpj
 138             qtbpj = -sin*wa(k) + cos*qtbpj
 139             wa(k) = temp
 140 c
 141 c           accumulate the transformation in the row of s.
 142 c
 143             kp1 = k + 1
 144             if (n .lt. kp1) go to 70
 145             do 60 i = kp1, n
 146                temp = cos*r(i,k) + sin*sdiag(i)
 147                sdiag(i) = -sin*r(i,k) + cos*sdiag(i)
 148                r(i,k) = temp
 149    60          continue
 150    70       continue
 151    80       continue
 152    90    continue
 153 c
 154 c        store the diagonal element of s and restore
 155 c        the corresponding diagonal element of r.
 156 c
 157          sdiag(j) = r(j,j)
 158          r(j,j) = x(j)
 159   100    continue
 160 c
 161 c     solve the triangular system for z. if the system is
 162 c     singular, then obtain a least squares solution.
 163 c
 164       nsing = n
 165       do 110 j = 1, n
 166          if (sdiag(j) .eq. zero .and. nsing .eq. n) nsing = j - 1
 167          if (nsing .lt. n) wa(j) = zero
 168   110    continue
 169       if (nsing .lt. 1) go to 150
 170       do 140 k = 1, nsing
 171          j = nsing - k + 1
 172          sum = zero
 173          jp1 = j + 1
 174          if (nsing .lt. jp1) go to 130
 175          do 120 i = jp1, nsing
 176             sum = sum + r(i,j)*wa(i)
 177   120       continue
 178   130    continue
 179          wa(j) = (wa(j) - sum)/sdiag(j)
 180   140    continue
 181   150 continue
 182 c
 183 c     permute the components of z back to components of x.
 184 c
 185       do 160 j = 1, n
 186          l = ipvt(j)
 187          x(l) = wa(j)
 188   160    continue
 189       return
 190 c
 191 c     last card of subroutine qrsolv.
 192 c
 193       end