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     ! 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.