Return to sgefa.f CVS log

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

Annotation of researchv10no/cmd/view2d/num/sgefa.f, revision 1.1.1.1

1.1       root        1:       subroutine sgefa(a,lda,n,ipvt,info)
                      2:       integer lda,n,ipvt(1),info
                      3:       real a(lda,1)
                      4: c
                      5: c     sgefa factors a real matrix by gaussian elimination.
                      6: c
                      7: c     sgefa is usually called by sgeco, but it can be called
                      8: c     directly with a saving in time if  rcond  is not needed.
                      9: c     (time for sgeco) = (1 + 9/n)*(time for sgefa) .
                     10: c
                     11: c     on entry
                     12: c
                     13: c        a       real(lda, n)
                     14: c                the matrix to be factored.
                     15: c
                     16: c        lda     integer
                     17: c                the leading dimension of the array  a .
                     18: c
                     19: c        n       integer
                     20: c                the order of the matrix  a .
                     21: c
                     22: c     on return
                     23: c
                     24: c        a       an upper triangular matrix and the multipliers
                     25: c                which were used to obtain it.
                     26: c                the factorization can be written  a = l*u  where
                     27: c                l  is a product of permutation and unit lower
                     28: c                triangular matrices and  u  is upper triangular.
                     29: c
                     30: c        ipvt    integer(n)
                     31: c                an integer vector of pivot indices.
                     32: c
                     33: c        info    integer
                     34: c                = 0  normal value.
                     35: c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
                     36: c                     condition for this subroutine, but it does
                     37: c                     indicate that sgesl or sgedi will divide by zero
                     38: c                     if called.  use  rcond  in sgeco for a reliable
                     39: c                     indication of singularity.
                     40: c
                     41: c     linpack. this version dated 08/14/78 .
                     42: c     cleve moler, university of new mexico, argonne national lab.
                     43: c
                     44: c     subroutines and functions
                     45: c
                     46: c     blas saxpy,sscal,isamax
                     47: c
                     48: c     internal variables
                     49: c
                     50:       real t
                     51:       integer isamax,j,k,kp1,l,nm1
                     52: c
                     53: c
                     54: c     gaussian elimination with partial pivoting
                     55: c
                     56:       info = 0
                     57:       nm1 = n - 1
                     58:       if (nm1 .lt. 1) go to 70
                     59:       do 60 k = 1, nm1
                     60:          kp1 = k + 1
                     61: c
                     62: c        find l = pivot index
                     63: c
                     64:          l = isamax(n-k+1,a(k,k),1) + k - 1
                     65:          ipvt(k) = l
                     66: c
                     67: c        zero pivot implies this column already triangularized
                     68: c
                     69:          if (a(l,k) .eq. 0.0e0) go to 40
                     70: c
                     71: c           interchange if necessary
                     72: c
                     73:             if (l .eq. k) go to 10
                     74:                t = a(l,k)
                     75:                a(l,k) = a(k,k)
                     76:                a(k,k) = t
                     77:    10       continue
                     78: c
                     79: c           compute multipliers
                     80: c
                     81:             t = -1.0e0/a(k,k)
                     82:             call sscal(n-k,t,a(k+1,k),1)
                     83: c
                     84: c           row elimination with column indexing
                     85: c
                     86:             do 30 j = kp1, n
                     87:                t = a(l,j)
                     88:                if (l .eq. k) go to 20
                     89:                   a(l,j) = a(k,j)
                     90:                   a(k,j) = t
                     91:    20          continue
                     92:                call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
                     93:    30       continue
                     94:          go to 50
                     95:    40    continue
                     96:             info = k
                     97:    50    continue
                     98:    60 continue
                     99:    70 continue
                    100:       ipvt(n) = n
                    101:       if (a(n,n) .eq. 0.0e0) info = n
                    102:       return
                    103:       end

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