Return to sqrdc.f CVS log

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

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

1.1     ! root        1:       subroutine sqrdc(x,ldx,n,p,qraux,jpvt,work,job)
        !             2:       integer ldx,n,p,job
        !             3:       integer jpvt(1)
        !             4:       real x(ldx,1),qraux(1),work(1)
        !             5: c
        !             6: c     sqrdc uses householder transformations to compute the qr
        !             7: c     factorization of an n by p matrix x.  column pivoting
        !             8: c     based on the 2-norms of the reduced columns may be
        !             9: c     performed at the users option.
        !            10: c
        !            11: c     on entry
        !            12: c
        !            13: c        x       real(ldx,p), where ldx .ge. n.
        !            14: c                x contains the matrix whose decomposition is to be
        !            15: c                computed.
        !            16: c
        !            17: c        ldx     integer.
        !            18: c                ldx is the leading dimension of the array x.
        !            19: c
        !            20: c        n       integer.
        !            21: c                n is the number of rows of the matrix x.
        !            22: c
        !            23: c        p       integer.
        !            24: c                p is the number of columns of the matrix x.
        !            25: c
        !            26: c        jpvt    integer(p).
        !            27: c                jpvt contains integers that control the selection
        !            28: c                of the pivot columns.  the k-th column x(k) of x
        !            29: c                is placed in one of three classes according to the
        !            30: c                value of jpvt(k).
        !            31: c
        !            32: c                   if jpvt(k) .gt. 0, then x(k) is an initial
        !            33: c                                      column.
        !            34: c
        !            35: c                   if jpvt(k) .eq. 0, then x(k) is a free column.
        !            36: c
        !            37: c                   if jpvt(k) .lt. 0, then x(k) is a final column.
        !            38: c
        !            39: c                before the decomposition is computed, initial columns
        !            40: c                are moved to the beginning of the array x and final
        !            41: c                columns to the end.  both initial and final columns
        !            42: c                are frozen in place during the computation and only
        !            43: c                free columns are moved.  at the k-th stage of the
        !            44: c                reduction, if x(k) is occupied by a free column
        !            45: c                it is interchanged with the free column of largest
        !            46: c                reduced norm.  jpvt is not referenced if
        !            47: c                job .eq. 0.
        !            48: c
        !            49: c        work    real(p).
        !            50: c                work is a work array.  work is not referenced if
        !            51: c                job .eq. 0.
        !            52: c
        !            53: c        job     integer.
        !            54: c                job is an integer that initiates column pivoting.
        !            55: c                if job .eq. 0, no pivoting is done.
        !            56: c                if job .ne. 0, pivoting is done.
        !            57: c
        !            58: c     on return
        !            59: c
        !            60: c        x       x contains in its upper triangle the upper
        !            61: c                triangular matrix r of the qr factorization.
        !            62: c                below its diagonal x contains information from
        !            63: c                which the orthogonal part of the decomposition
        !            64: c                can be recovered.  note that if pivoting has
        !            65: c                been requested, the decomposition is not that
        !            66: c                of the original matrix x but that of x
        !            67: c                with its columns permuted as described by jpvt.
        !            68: c
        !            69: c        qraux   real(p).
        !            70: c                qraux contains further information required to recover
        !            71: c                the orthogonal part of the decomposition.
        !            72: c
        !            73: c        jpvt    jpvt(k) contains the index of the column of the
        !            74: c                original matrix that has been interchanged into
        !            75: c                the k-th column, if pivoting was requested.
        !            76: c
        !            77: c     linpack. this version dated 08/14/78 .
        !            78: c     g.w. stewart, university of maryland, argonne national lab.
        !            79: c
        !            80: c     sqrdc uses the following functions and subprograms.
        !            81: c
        !            82: c     blas saxpy,sdot,sscal,sswap,snrm2
        !            83: c     fortran abs,amax1,min0,sqrt
        !            84: c
        !            85: c     internal variables
        !            86: c
        !            87:       integer j,jp,l,lp1,lup,maxj,pl,pu
        !            88:       real maxnrm,snrm2,tt
        !            89:       real sdot,nrmxl,t
        !            90:       logical negj,swapj
        !            91: c
        !            92: c
        !            93:       pl = 1
        !            94:       pu = 0
        !            95:       if (job .eq. 0) go to 60
        !            96: c
        !            97: c        pivoting has been requested.  rearrange the columns
        !            98: c        according to jpvt.
        !            99: c
        !           100:          do 20 j = 1, p
        !           101:             swapj = jpvt(j) .gt. 0
        !           102:             negj = jpvt(j) .lt. 0
        !           103:             jpvt(j) = j
        !           104:             if (negj) jpvt(j) = -j
        !           105:             if (.not.swapj) go to 10
        !           106:                if (j .ne. pl) call sswap(n,x(1,pl),1,x(1,j),1)
        !           107:                jpvt(j) = jpvt(pl)
        !           108:                jpvt(pl) = j
        !           109:                pl = pl + 1
        !           110:    10       continue
        !           111:    20    continue
        !           112:          pu = p
        !           113:          do 50 jj = 1, p
        !           114:             j = p - jj + 1
        !           115:             if (jpvt(j) .ge. 0) go to 40
        !           116:                jpvt(j) = -jpvt(j)
        !           117:                if (j .eq. pu) go to 30
        !           118:                   call sswap(n,x(1,pu),1,x(1,j),1)
        !           119:                   jp = jpvt(pu)
        !           120:                   jpvt(pu) = jpvt(j)
        !           121:                   jpvt(j) = jp
        !           122:    30          continue
        !           123:                pu = pu - 1
        !           124:    40       continue
        !           125:    50    continue
        !           126:    60 continue
        !           127: c
        !           128: c     compute the norms of the free columns.
        !           129: c
        !           130:       if (pu .lt. pl) go to 80
        !           131:       do 70 j = pl, pu
        !           132:          qraux(j) = snrm2(n,x(1,j),1)
        !           133:          work(j) = qraux(j)
        !           134:    70 continue
        !           135:    80 continue
        !           136: c
        !           137: c     perform the householder reduction of x.
        !           138: c
        !           139:       lup = min0(n,p)
        !           140:       do 200 l = 1, lup
        !           141:          if (l .lt. pl .or. l .ge. pu) go to 120
        !           142: c
        !           143: c           locate the column of largest norm and bring it
        !           144: c           into the pivot position.
        !           145: c
        !           146:             maxnrm = 0.0e0
        !           147:             maxj = l
        !           148:             do 100 j = l, pu
        !           149:                if (qraux(j) .le. maxnrm) go to 90
        !           150:                   maxnrm = qraux(j)
        !           151:                   maxj = j
        !           152:    90          continue
        !           153:   100       continue
        !           154:             if (maxj .eq. l) go to 110
        !           155:                call sswap(n,x(1,l),1,x(1,maxj),1)
        !           156:                qraux(maxj) = qraux(l)
        !           157:                work(maxj) = work(l)
        !           158:                jp = jpvt(maxj)
        !           159:                jpvt(maxj) = jpvt(l)
        !           160:                jpvt(l) = jp
        !           161:   110       continue
        !           162:   120    continue
        !           163:          qraux(l) = 0.0e0
        !           164:          if (l .eq. n) go to 190
        !           165: c
        !           166: c           compute the householder transformation for column l.
        !           167: c
        !           168:             nrmxl = snrm2(n-l+1,x(l,l),1)
        !           169:             if (nrmxl .eq. 0.0e0) go to 180
        !           170:                if (x(l,l) .ne. 0.0e0) nrmxl = sign(nrmxl,x(l,l))
        !           171:                call sscal(n-l+1,1.0e0/nrmxl,x(l,l),1)
        !           172:                x(l,l) = 1.0e0 + x(l,l)
        !           173: c
        !           174: c              apply the transformation to the remaining columns,
        !           175: c              updating the norms.
        !           176: c
        !           177:                lp1 = l + 1
        !           178:                if (p .lt. lp1) go to 170
        !           179:                do 160 j = lp1, p
        !           180:                   t = -sdot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
        !           181:                   call saxpy(n-l+1,t,x(l,l),1,x(l,j),1)
        !           182:                   if (j .lt. pl .or. j .gt. pu) go to 150
        !           183:                   if (qraux(j) .eq. 0.0e0) go to 150
        !           184:                      tt = 1.0e0 - (abs(x(l,j))/qraux(j))**2
        !           185:                      tt = amax1(tt,0.0e0)
        !           186:                      t = tt
        !           187:                      tt = 1.0e0 + 0.05e0*tt*(qraux(j)/work(j))**2
        !           188:                      if (tt .eq. 1.0e0) go to 130
        !           189:                         qraux(j) = qraux(j)*sqrt(t)
        !           190:                      go to 140
        !           191:   130                continue
        !           192:                         qraux(j) = snrm2(n-l,x(l+1,j),1)
        !           193:                         work(j) = qraux(j)
        !           194:   140                continue
        !           195:   150             continue
        !           196:   160          continue
        !           197:   170          continue
        !           198: c
        !           199: c              save the transformation.
        !           200: c
        !           201:                qraux(l) = x(l,l)
        !           202:                x(l,l) = -nrmxl
        !           203:   180       continue
        !           204:   190    continue
        !           205:   200 continue
        !           206:       return
        !           207:       end

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