Return to sgesl.f CVS log

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

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

1.1     ! root        1:       subroutine sgesl(a,lda,n,ipvt,b,job)
        !             2:       integer lda,n,ipvt(1),job
        !             3:       real a(lda,1),b(1)
        !             4: c
        !             5: c     sgesl solves the real system
        !             6: c     a * x = b  or  trans(a) * x = b
        !             7: c     using the factors computed by sgeco or sgefa.
        !             8: c
        !             9: c     on entry
        !            10: c
        !            11: c        a       real(lda, n)
        !            12: c                the output from sgeco or sgefa.
        !            13: c
        !            14: c        lda     integer
        !            15: c                the leading dimension of the array  a .
        !            16: c
        !            17: c        n       integer
        !            18: c                the order of the matrix  a .
        !            19: c
        !            20: c        ipvt    integer(n)
        !            21: c                the pivot vector from sgeco or sgefa.
        !            22: c
        !            23: c        b       real(n)
        !            24: c                the right hand side vector.
        !            25: c
        !            26: c        job     integer
        !            27: c                = 0         to solve  a*x = b ,
        !            28: c                = nonzero   to solve  trans(a)*x = b  where
        !            29: c                            trans(a)  is the transpose.
        !            30: c
        !            31: c     on return
        !            32: c
        !            33: c        b       the solution vector  x .
        !            34: c
        !            35: c     error condition
        !            36: c
        !            37: c        a division by zero will occur if the input factor contains a
        !            38: c        zero on the diagonal.  technically this indicates singularity
        !            39: c        but it is often caused by improper arguments or improper
        !            40: c        setting of lda .  it will not occur if the subroutines are
        !            41: c        called correctly and if sgeco has set rcond .gt. 0.0
        !            42: c        or sgefa has set info .eq. 0 .
        !            43: c
        !            44: c     to compute  inverse(a) * c  where  c  is a matrix
        !            45: c     with  p  columns
        !            46: c           call sgeco(a,lda,n,ipvt,rcond,z)
        !            47: c           if (rcond is too small) go to ...
        !            48: c           do 10 j = 1, p
        !            49: c              call sgesl(a,lda,n,ipvt,c(1,j),0)
        !            50: c        10 continue
        !            51: c
        !            52: c     linpack. this version dated 08/14/78 .
        !            53: c     cleve moler, university of new mexico, argonne national lab.
        !            54: c
        !            55: c     subroutines and functions
        !            56: c
        !            57: c     blas saxpy,sdot
        !            58: c
        !            59: c     internal variables
        !            60: c
        !            61:       real sdot,t
        !            62:       integer k,kb,l,nm1
        !            63: c
        !            64:       nm1 = n - 1
        !            65:       if (job .ne. 0) go to 50
        !            66: c
        !            67: c        job = 0 , solve  a * x = b
        !            68: c        first solve  l*y = b
        !            69: c
        !            70:          if (nm1 .lt. 1) go to 30
        !            71:          do 20 k = 1, nm1
        !            72:             l = ipvt(k)
        !            73:             t = b(l)
        !            74:             if (l .eq. k) go to 10
        !            75:                b(l) = b(k)
        !            76:                b(k) = t
        !            77:    10       continue
        !            78:             call saxpy(n-k,t,a(k+1,k),1,b(k+1),1)
        !            79:    20    continue
        !            80:    30    continue
        !            81: c
        !            82: c        now solve  u*x = y
        !            83: c
        !            84:          do 40 kb = 1, n
        !            85:             k = n + 1 - kb
        !            86:             b(k) = b(k)/a(k,k)
        !            87:             t = -b(k)
        !            88:             call saxpy(k-1,t,a(1,k),1,b(1),1)
        !            89:    40    continue
        !            90:       go to 100
        !            91:    50 continue
        !            92: c
        !            93: c        job = nonzero, solve  trans(a) * x = b
        !            94: c        first solve  trans(u)*y = b
        !            95: c
        !            96:          do 60 k = 1, n
        !            97:             t = sdot(k-1,a(1,k),1,b(1),1)
        !            98:             b(k) = (b(k) - t)/a(k,k)
        !            99:    60    continue
        !           100: c
        !           101: c        now solve trans(l)*x = y
        !           102: c
        !           103:          if (nm1 .lt. 1) go to 90
        !           104:          do 80 kb = 1, nm1
        !           105:             k = n - kb
        !           106:             b(k) = b(k) + sdot(n-k,a(k+1,k),1,b(k+1),1)
        !           107:             l = ipvt(k)
        !           108:             if (l .eq. k) go to 70
        !           109:                t = b(l)
        !           110:                b(l) = b(k)
        !           111:                b(k) = t
        !           112:    70       continue
        !           113:    80    continue
        !           114:    90    continue
        !           115:   100 continue
        !           116:       return
        !           117:       end

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