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

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.