share/colnew/fortran/factrb.f

   1       SUBROUTINE FACTRB ( W, IPIVOT, D, NROW, NCOL, LAST, INFO)
   2 C
   3 C********************************************************************
   4 C
   5 C     adapted from p.132 of  element.numer.analysis  by conte-de boor
   6 C
   7 C     constructs a partial plu factorization, corresponding to steps
   8 C      1,..., last   in gauss elimination, for the matrix  w  of
   9 C      order ( nrow ,  ncol ), using pivoting of scaled rows.
  10 C
  11 C     parameters
  12 C       w       contains the (nrow,ncol) matrix to be partially factored
  13 C               on input, and the partial factorization on output.
  14 C       ipivot  an integer array of length last containing a record of
  15 C               the pivoting strategy used; explicit interchanges
  16 C               are used for pivoting.
  17 C       d       a work array of length nrow used to store row sizes
  18 C               temporarily.
  19 C       nrow    number of rows of w.
  20 C       ncol    number of columns of w.
  21 C       last    number of elimination steps to be carried out.
  22 C       info    on output, zero if the matrix is found to be non-
  23 C               singular, in case a zero pivot was encountered in row
  24 C               n,  info = n on output.
  25 C
  26 C**********************************************************************
  27 C
  28       INTEGER IPIVOT(NROW),NCOL,LAST,INFO, I,J,K,L,KP1
  29       DOUBLE PRECISION W(NROW,NCOL),D(NROW), COLMAX,T,S
  30       DOUBLE PRECISION DABS,DMAX1
  31 C
  32 C...  initialize  d
  33 C
  34       DO 10 I = 1, NROW
  35         D(I) = 0.D0
  36    10 CONTINUE
  37       DO 20 J = 1, NCOL
  38         DO 20 I = 1, NROW
  39           D(I) = DMAX1( D(I) , DABS(W(I,J)))
  40    20 CONTINUE
  41 C
  42 C...  gauss elimination with pivoting of scaled rows, loop over
  43 C...  k=1,.,last
  44 C
  45       K = 1
  46 C
  47 C...  as pivot row for k-th step, pick among the rows not yet used,
  48 C...  i.e., from rows  k ,..., nrow , the one whose k-th entry
  49 C...  (compared to the row size) is largest. then, if this row
  50 C...  does not turn out to be row k, interchange row k with this
  51 C...  particular row and redefine ipivot(k).
  52 C
  53    30      CONTINUE
  54            IF ( D(K) .EQ. 0.D0 )                    GO TO 90
  55            IF (K .EQ. NROW)                         GO TO 80
  56            L = K
  57            KP1 = K+1
  58            COLMAX = DABS(W(K,K)) / D(K)
  59 C
  60 C...       find the (relatively) largest pivot
  61 C
  62            DO 40 I = KP1, NROW
  63              IF ( DABS(W(I,K)) .LE. COLMAX * D(I) ) GO TO 40
  64              COLMAX = DABS(W(I,K)) / D(I)
  65              L = I
  66    40      CONTINUE
  67            IPIVOT(K) = L
  68            T = W(L,K)
  69            S = D(L)
  70            IF ( L .EQ. K )                          GO TO 50
  71              W(L,K) = W(K,K)
  72              W(K,K) = T
  73              D(L) = D(K)
  74              D(K) = S
  75    50      CONTINUE
  76 C
  77 C...       if pivot element is too small in absolute value, declare
  78 C...       matrix to be noninvertible and quit.
  79 C
  80            IF ( DABS(T)+D(K) .LE. D(K) )            GO TO 90
  81 C
  82 C...       otherwise, subtract the appropriate multiple of the pivot
  83 C...       row from remaining rows, i.e., the rows (k+1),..., (nrow)
  84 C...       to make k-th entry zero. save the multiplier in its place.
  85 C...       for high performance do this operations column oriented.
  86 C
  87            T = -1.0D0/T
  88            DO 60 I = KP1, NROW
  89    60        W(I,K) = W(I,K) * T
  90            DO 70 J=KP1,NCOL
  91              T = W(L,J)
  92              IF ( L .EQ. K )                        GO TO 62
  93                W(L,J) = W(K,J)
  94                W(K,J) = T
  95    62        IF ( T .EQ. 0.D0 )                     GO TO 70
  96              DO 64 I = KP1, NROW
  97    64           W(I,J) = W(I,J) + W(I,K) * T
  98    70      CONTINUE
  99            K = KP1
 100 C
 101 C...       check for having reached the next block.
 102 C
 103            IF ( K .LE. LAST )                       GO TO 30
 104       RETURN
 105 C
 106 C...  if  last  .eq. nrow , check now that pivot element in last row
 107 C...  is nonzero.
 108 C
 109    80 IF( DABS(W(NROW,NROW))+D(NROW) .GT. D(NROW) ) RETURN
 110 C
 111 C...  singularity flag set
 112 C
 113    90 INFO = K
 114       RETURN
 115       END