share/colnew/fortran/dgefa.f

   1       subroutine dgefa(a,lda,n,ipvt,info)
   2       integer lda,n,ipvt(1),info
   3       double precision a(lda,1)
   4 c
   5 c     dgefa factors a double precision matrix by gaussian elimination.
   6 c
   7 c     dgefa is usually called by dgeco, but it can be called
   8 c     directly with a saving in time if  rcond  is not needed.
   9 c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
  10 c
  11 c     on entry
  12 c
  13 c        a       double precision(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 dgesl or dgedi will divide by zero
  38 c                     if called.  use  rcond  in dgeco 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 daxpy,dscal,idamax
  47 c
  48 c     internal variables
  49 c
  50       double precision t
  51       integer idamax,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 = idamax(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.0d0) 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.0d0/a(k,k)
  82             call dscal(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 daxpy(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.0d0) info = n
 102       return
 103       end