Return to sqrsl.f CVS log

Up to [Research Unix] / researchv10no / cmd / view2d / num

Annotation of researchv10no/cmd/view2d/num/sqrsl.f, revision 1.1

1.1     ! root        1:       subroutine sqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info)
        !             2:       integer ldx,n,k,job,info
        !             3:       real x(ldx,1),qraux(1),y(1),qy(1),qty(1),b(1),rsd(1),xb(1)
        !             4: c
        !             5: c     sqrsl applies the output of sqrdc to compute coordinate
        !             6: c     transformations, projections, and least squares solutions.
        !             7: c     for k .le. min(n,p), let xk be the matrix
        !             8: c
        !             9: c            xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k)))
        !            10: c
        !            11: c     formed from columnns jpvt(1), ... ,jpvt(k) of the original
        !            12: c     n x p matrix x that was input to sqrdc (if no pivoting was
        !            13: c     done, xk consists of the first k columns of x in their
        !            14: c     original order).  sqrdc produces a factored orthogonal matrix q
        !            15: c     and an upper triangular matrix r such that
        !            16: c
        !            17: c              xk = q * (r)
        !            18: c                       (0)
        !            19: c
        !            20: c     this information is contained in coded form in the arrays
        !            21: c     x and qraux.
        !            22: c
        !            23: c     on entry
        !            24: c
        !            25: c        x      real(ldx,p).
        !            26: c               x contains the output of sqrdc.
        !            27: c
        !            28: c        ldx    integer.
        !            29: c               ldx is the leading dimension of the array x.
        !            30: c
        !            31: c        n      integer.
        !            32: c               n is the number of rows of the matrix xk.  it must
        !            33: c               have the same value as n in sqrdc.
        !            34: c
        !            35: c        k      integer.
        !            36: c               k is the number of columns of the matrix xk.  k
        !            37: c               must nnot be greater than min(n,p), where p is the
        !            38: c               same as in the calling sequence to sqrdc.
        !            39: c
        !            40: c        qraux  real(p).
        !            41: c               qraux contains the auxiliary output from sqrdc.
        !            42: c
        !            43: c        y      real(n)
        !            44: c               y contains an n-vector that is to be manipulated
        !            45: c               by sqrsl.
        !            46: c
        !            47: c        job    integer.
        !            48: c               job specifies what is to be computed.  job has
        !            49: c               the decimal expansion abcde, with the following
        !            50: c               meaning.
        !            51: c
        !            52: c                    if a.ne.0, compute qy.
        !            53: c                    if b,c,d, or e .ne. 0, compute qty.
        !            54: c                    if c.ne.0, compute b.
        !            55: c                    if d.ne.0, compute rsd.
        !            56: c                    if e.ne.0, compute xb.
        !            57: c
        !            58: c               note that a request to compute b, rsd, or xb
        !            59: c               automatically triggers the computation of qty, for
        !            60: c               which an array must be provided in the calling
        !            61: c               sequence.
        !            62: c
        !            63: c     on return
        !            64: c
        !            65: c        qy     real(n).
        !            66: c               qy conntains q*y, if its computation has been
        !            67: c               requested.
        !            68: c
        !            69: c        qty    real(n).
        !            70: c               qty contains trans(q)*y, if its computation has
        !            71: c               been requested.  here trans(q) is the
        !            72: c               transpose of the matrix q.
        !            73: c
        !            74: c        b      real(k)
        !            75: c               b contains the solution of the least squares problem
        !            76: c
        !            77: c                    minimize norm2(y - xk*b),
        !            78: c
        !            79: c               if its computation has been requested.  (note that
        !            80: c               if pivoting was requested in sqrdc, the j-th
        !            81: c               component of b will be associated with column jpvt(j)
        !            82: c               of the original matrix x that was input into sqrdc.)
        !            83: c
        !            84: c        rsd    real(n).
        !            85: c               rsd contains the least squares residual y - xk*b,
        !            86: c               if its computation has been requested.  rsd is
        !            87: c               also the orthogonal projection of y onto the
        !            88: c               orthogonal complement of the column space of xk.
        !            89: c
        !            90: c        xb     real(n).
        !            91: c               xb contains the least squares approximation xk*b,
        !            92: c               if its computation has been requested.  xb is also
        !            93: c               the orthogonal projection of y onto the column space
        !            94: c               of x.
        !            95: c
        !            96: c        info   integer.
        !            97: c               info is zero unless the computation of b has
        !            98: c               been requested and r is exactly singular.  in
        !            99: c               this case, info is the index of the first zero
        !           100: c               diagonal element of r and b is left unaltered.
        !           101: c
        !           102: c     the parameters qy, qty, b, rsd, and xb are not referenced
        !           103: c     if their computation is not requested and in this case
        !           104: c     can be replaced by dummy variables in the calling program.
        !           105: c     to save storage, the user may in some cases use the same
        !           106: c     array for different parameters in the calling sequence.  a
        !           107: c     frequently occuring example is when one wishes to compute
        !           108: c     any of b, rsd, or xb and does not need y or qty.  in this
        !           109: c     case one may identify y, qty, and one of b, rsd, or xb, while
        !           110: c     providing separate arrays for anything else that is to be
        !           111: c     computed.  thus the calling sequence
        !           112: c
        !           113: c          call sqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info)
        !           114: c
        !           115: c     will result in the computation of b and rsd, with rsd
        !           116: c     overwriting y.  more generally, each item in the following
        !           117: c     list contains groups of permissible identifications for
        !           118: c     a single callinng sequence.
        !           119: c
        !           120: c          1. (y,qty,b) (rsd) (xb) (qy)
        !           121: c
        !           122: c          2. (y,qty,rsd) (b) (xb) (qy)
        !           123: c
        !           124: c          3. (y,qty,xb) (b) (rsd) (qy)
        !           125: c
        !           126: c          4. (y,qy) (qty,b) (rsd) (xb)
        !           127: c
        !           128: c          5. (y,qy) (qty,rsd) (b) (xb)
        !           129: c
        !           130: c          6. (y,qy) (qty,xb) (b) (rsd)
        !           131: c
        !           132: c     in any group the value returned in the array allocated to
        !           133: c     the group corresponds to the last member of the group.
        !           134: c
        !           135: c     linpack. this version dated 08/14/78 .
        !           136: c     g.w. stewart, university of maryland, argonne national lab.
        !           137: c
        !           138: c     sqrsl uses the following functions and subprograms.
        !           139: c
        !           140: c     blas saxpy,scopy,sdot
        !           141: c     fortran abs,min0,mod
        !           142: c
        !           143: c     internal variables
        !           144: c
        !           145:       integer i,j,jj,ju,kp1
        !           146:       real sdot,t,temp
        !           147:       logical cb,cqy,cqty,cr,cxb
        !           148: c
        !           149: c
        !           150: c     set info flag.
        !           151: c
        !           152:       info = 0
        !           153: c
        !           154: c     determine what is to be computed.
        !           155: c
        !           156:       cqy = job/10000 .ne. 0
        !           157:       cqty = mod(job,10000) .ne. 0
        !           158:       cb = mod(job,1000)/100 .ne. 0
        !           159:       cr = mod(job,100)/10 .ne. 0
        !           160:       cxb = mod(job,10) .ne. 0
        !           161:       ju = min0(k,n-1)
        !           162: c
        !           163: c     special action when n=1.
        !           164: c
        !           165:       if (ju .ne. 0) go to 40
        !           166:          if (cqy) qy(1) = y(1)
        !           167:          if (cqty) qty(1) = y(1)
        !           168:          if (cxb) xb(1) = y(1)
        !           169:          if (.not.cb) go to 30
        !           170:             if (x(1,1) .ne. 0.0e0) go to 10
        !           171:                info = 1
        !           172:             go to 20
        !           173:    10       continue
        !           174:                b(1) = y(1)/x(1,1)
        !           175:    20       continue
        !           176:    30    continue
        !           177:          if (cr) rsd(1) = 0.0e0
        !           178:       go to 250
        !           179:    40 continue
        !           180: c
        !           181: c        set up to compute qy or qty.
        !           182: c
        !           183:          if (cqy) call scopy(n,y,1,qy,1)
        !           184:          if (cqty) call scopy(n,y,1,qty,1)
        !           185:          if (.not.cqy) go to 70
        !           186: c
        !           187: c           compute qy.
        !           188: c
        !           189:             do 60 jj = 1, ju
        !           190:                j = ju - jj + 1
        !           191:                if (qraux(j) .eq. 0.0e0) go to 50
        !           192:                   temp = x(j,j)
        !           193:                   x(j,j) = qraux(j)
        !           194:                   t = -sdot(n-j+1,x(j,j),1,qy(j),1)/x(j,j)
        !           195:                   call saxpy(n-j+1,t,x(j,j),1,qy(j),1)
        !           196:                   x(j,j) = temp
        !           197:    50          continue
        !           198:    60       continue
        !           199:    70    continue
        !           200:          if (.not.cqty) go to 100
        !           201: c
        !           202: c           compute trans(q)*y.
        !           203: c
        !           204:             do 90 j = 1, ju
        !           205:                if (qraux(j) .eq. 0.0e0) go to 80
        !           206:                   temp = x(j,j)
        !           207:                   x(j,j) = qraux(j)
        !           208:                   t = -sdot(n-j+1,x(j,j),1,qty(j),1)/x(j,j)
        !           209:                   call saxpy(n-j+1,t,x(j,j),1,qty(j),1)
        !           210:                   x(j,j) = temp
        !           211:    80          continue
        !           212:    90       continue
        !           213:   100    continue
        !           214: c
        !           215: c        set up to compute b, rsd, or xb.
        !           216: c
        !           217:          if (cb) call scopy(k,qty,1,b,1)
        !           218:          kp1 = k + 1
        !           219:          if (cxb) call scopy(k,qty,1,xb,1)
        !           220:          if (cr .and. k .lt. n) call scopy(n-k,qty(kp1),1,rsd(kp1),1)
        !           221:          if (.not.cxb .or. kp1 .gt. n) go to 120
        !           222:             do 110 i = kp1, n
        !           223:                xb(i) = 0.0e0
        !           224:   110       continue
        !           225:   120    continue
        !           226:          if (.not.cr) go to 140
        !           227:             do 130 i = 1, k
        !           228:                rsd(i) = 0.0e0
        !           229:   130       continue
        !           230:   140    continue
        !           231:          if (.not.cb) go to 190
        !           232: c
        !           233: c           compute b.
        !           234: c
        !           235:             do 170 jj = 1, k
        !           236:                j = k - jj + 1
        !           237:                if (x(j,j) .ne. 0.0e0) go to 150
        !           238:                   info = j
        !           239: c           ......exit
        !           240:                   go to 180
        !           241:   150          continue
        !           242:                b(j) = b(j)/x(j,j)
        !           243:                if (j .eq. 1) go to 160
        !           244:                   t = -b(j)
        !           245:                   call saxpy(j-1,t,x(1,j),1,b,1)
        !           246:   160          continue
        !           247:   170       continue
        !           248:   180       continue
        !           249:   190    continue
        !           250:          if (.not.cr .and. .not.cxb) go to 240
        !           251: c
        !           252: c           compute rsd or xb as required.
        !           253: c
        !           254:             do 230 jj = 1, ju
        !           255:                j = ju - jj + 1
        !           256:                if (qraux(j) .eq. 0.0e0) go to 220
        !           257:                   temp = x(j,j)
        !           258:                   x(j,j) = qraux(j)
        !           259:                   if (.not.cr) go to 200
        !           260:                      t = -sdot(n-j+1,x(j,j),1,rsd(j),1)/x(j,j)
        !           261:                      call saxpy(n-j+1,t,x(j,j),1,rsd(j),1)
        !           262:   200             continue
        !           263:                   if (.not.cxb) go to 210
        !           264:                      t = -sdot(n-j+1,x(j,j),1,xb(j),1)/x(j,j)
        !           265:                      call saxpy(n-j+1,t,x(j,j),1,xb(j),1)
        !           266:   210             continue
        !           267:                   x(j,j) = temp
        !           268:   220          continue
        !           269:   230       continue
        !           270:   240    continue
        !           271:   250 continue
        !           272:       return
        !           273:       end

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.