CFPRIME DGEFA2                                                          DGEFA2      1
C***********************************************************************RIGHTE      1
C*                                                                     *RIGHTE      2
C*        Extended Transient-Midterm Stability Program (ETMSP)         *RIGHTE      3
C*        ====================================================         *RIGHTE      4
C*                                                                     *RIGHTE      5
C*       Copyright (C) (1993) Electric Power Research Institute        *RIGHTE      6
C*                                                                     *RIGHTE      7
C* EPRI reserves all rights in the software and associated             *RIGHTE      8
C* documentation as delivered. The software and its associated         *RIGHTE      9
C* documention or any portion thereof may not be reproduced in any form*RIGHTE     10
C* whatsoever without the consent of EPRI, except as explicity provided*RIGHTE     11
C* for by license. Licensee shall submit a detailed description to EPRI*RIGHTE     12
C* of any changes or modifications to the software within two (2)      *RIGHTE     13
C* months after such changes or modifications are made. At EPRI's sole *RIGHTE     14
C* discretion, the software, as changed or modified, shall be given a  *RIGHTE     15
C* new designation sufficiently different from its then current        *RIGHTE     16
C* designation as to prevent mistake,confusion, or deception as between*RIGHTE     17
C* the current software as changed or modified.                        *RIGHTE     18
C*                                                                     *RIGHTE     19
C* A license under EPRI's rights in the software can be obtained       *RIGHTE     20
C* directly from EPRI.                                                 *RIGHTE     21
C*                                                                     *RIGHTE     22
C* This program and supporting materials may be independently available*RIGHTE     23
C* from Electric Power Software centre (EPSC) for an appropriate       *RIGHTE     24
C* distribution fee. The EPSC is a service of                          *RIGHTE     25
C*                                                                     *RIGHTE     26
C*                    Power Computing Company                          *RIGHTE     27
C*                    1930 Hi Line                                     *RIGHTE     28
C*                    Dallas, Tx 75207                                 *RIGHTE     29
C*                                                                     *RIGHTE     30
C* NEITHER EPRI, ANY MEMBER OF EPRI, NOR ANY PERSON OR  ORGANIZATION   *RIGHTE     31
C* ACTING ON BEHAVE OF THEM:                                           *RIGHTE     32
C*                                                                     *RIGHTE     33
C* 1. MAKE ANY WARRANTY OR REPRESENTATION WHATSOEVER, EXPRESS OR       *RIGHTE     34
C*    IMPLIED, INCLUDING ANY WARRANTY OF MERCHANTABILITY OR FITNESS OR *RIGHTE     35
C*    ANY PURPOSE WITH RESPECT TO THE SOFTWARE OR ITS ASSOCIATED       *RIGHTE     36
C*    DOCUMENTATION; OR                                                *RIGHTE     37
C*                                                                     *RIGHTE     38
C* 2. ASSUMES ANY LIABILITY WHATSOEVER WITH RESPECT TO ANY USE OF THE  *RIGHTE     39
C*    SOFTWARE OR ITS ASSOCIATED DOCUMENTATION OR ANY PORTION THEREOF  *RIGHTE     40
C*    OR WITH RESPECT TO ANY DAMAGES WHICH MAY RESULT FROM SUCH USE.   *RIGHTE     41
C*                                                                     *RIGHTE     42
C* THIS NOTICE MAY NOT BE REMOVED FROM THE SOFTWARE BY ANY USER        *RIGHTE     43
C* THEREOF.                                                            *RIGHTE     44
C*                                                                     *RIGHTE     45
C***********************************************************************RIGHTE     46
      SUBROUTINE DGEFA2 (A, LN, N, IPVT, INFO)                          DGEFA2      3
CD -------------------------------------------------------------------  DGEFA2      4
CD                                                                      DGEFA2      5
CD              *****     DGEFA2  ROUTINE   *****                       DGEFA2      6
CD                                                                      DGEFA2      7
CD -------------------------------------------------------------------  DGEFA2      8
CD                                                                      DGEFA2      9
CD    DGEFA2 factors a double precision matrix by Gaussian Elimination. DGEFA2     10
CD                                                                      DGEFA2     11
CD    DGEFA2 is usually called by DGECO, but it can be called           DGEFA2     12
CD    directly with a saving in time if  RCOND  is not needed.          DGEFA2     13
CD    (TIME FOR DGECO) = (1 + 9 / N) * (TIME FOR DGEFA2).               DGEFA2     14
CD                                                                      DGEFA2     15
CD    On Entry:                                                         DGEFA2     16
CD                                                                      DGEFA2     17
CD       A       DOUBLE PRECISION(LN, N)                                DGEFA2     18
CD               The matrix to be factored.                             DGEFA2     19
CD                                                                      DGEFA2     20
CD       LN     INTEGER                                                 DGEFA2     21
CD               The leading dimension of the array  A.                 DGEFA2     22
CD                                                                      DGEFA2     23
CD       N       INTEGER                                                DGEFA2     24
CD               The order of the matrix  A.                            DGEFA2     25
CD                                                                      DGEFA2     26
CD    On Return:                                                        DGEFA2     27
CD                                                                      DGEFA2     28
CD       A       An upper triangular matrix and the multipliers         DGEFA2     29
CD               which were used to obtain it.                          DGEFA2     30
CD               The factorization can be written  A = L * U  where     DGEFA2     31
CD               L  is a product of permutation and unit lower          DGEFA2     32
CD               triangular matrices and  U  is upper triangular.       DGEFA2     33
CD                                                                      DGEFA2     34
CD       IPVT    INTEGER(N)                                             DGEFA2     35
CD               An integer vector of pivot indices.                    DGEFA2     36
CD                                                                      DGEFA2     37
CD       INFO    INTEGER                                                DGEFA2     38
CD               = 0  Normal value.                                     DGEFA2     39
CD               = K  If  U(K,K) .EQ. 0.0 .  This is not an error       DGEFA2     40
CD                    condition for this subroutine, but it does        DGEFA2     41
CD                    indicate that DGESL2 or DGEDI will divide by zero DGEFA2     42
CD                    if called.  Use  RCOND  in DGECO for a reliable   DGEFA2     43
CD                    indication of singularity.                        DGEFA2     44
CD                                                                      DGEFA2     45
CD    LINPACK. THIS VERSION DATED 08/14/78 .                            DGEFA2     46
CD    CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.      DGEFA2     47
CD                                                                      DGEFA2     48
CD -------------------------------------------------------------------  DGEFA2     49
C                                                                       DGEFA2     50
C                                                                       DGEFA2     51
      DOUBLE PRECISION A(LN,LN), T                                      DGEFA2     52
C                                                                       DGEFA2     53
      DIMENSION          IPVT(*)                                        DGEFA2     54
C                                                                       DGEFA2     55
C                                                                       DGEFA2     56
C---------------------------------------------------------------------  DGEFA2     57
C                                                                       DGEFA2     58
C                                                                       DGEFA2     59
C     Gaussian Elimination with partial pivoting.                       DGEFA2     60
C                                                                       DGEFA2     61
      INFO = 0                                                          DGEFA2     62
      NM1 = N - 1                                                       DGEFA2     63
      IF (NM1 .LT. 1) GO TO 70                                          DGEFA2     64
      DO 60 K = 1, NM1                                                  DGEFA2     65
         KP1 = K + 1                                                    DGEFA2     66
C                                                                       DGEFA2     67
C        Find L = Pivot index.                                          DGEFA2     68
C                                                                       DGEFA2     69
         L = IDAMAX(N-K+1,A(K,K),1) + K - 1                             DGEFA2     70
         IPVT(K) = L                                                    DGEFA2     71
C                                                                       DGEFA2     72
C        Zero pivot implies this column already triangularized.         DGEFA2     73
C                                                                       DGEFA2     74
         IF (ABS(A(L,K)) .LE. 0.0D0) GO TO 40                           DGEFA2     75
C                                                                       DGEFA2     76
C           Interchange if necessary.                                   DGEFA2     77
C                                                                       DGEFA2     78
            IF (L .EQ. K) GO TO 10                                      DGEFA2     79
               T = A(L,K)                                               DGEFA2     80
               A(L,K) = A(K,K)                                          DGEFA2     81
               A(K,K) = T                                               DGEFA2     82
10          CONTINUE                                                    DGEFA2     83
C                                                                       DGEFA2     84
C           Compute multipliers.                                        DGEFA2     85
C                                                                       DGEFA2     86
            T = -1.0D0 / A(K,K)                                         DGEFA2     87
            CALL DSCAL (N-K, T, A(K+1,K), 1)                            DGEFA2     88
C                                                                       DGEFA2     89
C           Row Elimination with column indexing.                       DGEFA2     90
C                                                                       DGEFA2     91
            DO 30 J = KP1, N                                            DGEFA2     92
               T = A(L,J)                                               DGEFA2     93
               IF (L .EQ. K) GO TO 20                                   DGEFA2     94
                  A(L,J) = A(K,J)                                       DGEFA2     95
                  A(K,J) = T                                            DGEFA2     96
20             CONTINUE                                                 DGEFA2     97
               CALL DAXPY (N-K, T, A(K+1,K), 1, A(K+1,J), 1)            DGEFA2     98
30          CONTINUE                                                    DGEFA2     99
         GO TO 50                                                       DGEFA2    100
40       CONTINUE                                                       DGEFA2    101
            INFO = K                                                    DGEFA2    102
50       CONTINUE                                                       DGEFA2    103
60    CONTINUE                                                          DGEFA2    104
70    CONTINUE                                                          DGEFA2    105
      IPVT(N) = N                                                       DGEFA2    106
      IF (ABS(A(N,N)) .LE. 0.0D0) INFO = N                              DGEFA2    107
      RETURN                                                            DGEFA2    108
      END                                                               DGEFA2    109
