CFPRIME DGESL2                                                          DGESL2      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 DGESL2 (A, LDA, N, IPVT, B, JOB)                       DGESL2      3
CD -------------------------------------------------------------------  DGESL2      4
CD                                                                      DGESL2      5
CD              *****     DGESL2  ROUTINE   *****                       DGESL2      6
CD                                                                      DGESL2      7
CD -------------------------------------------------------------------  DGESL2      8
CD                                                                      DGESL2      9
CD    DGESL2 solves the double precision system                         DGESL2     10
CD    A * X = B  or  TRANS(A) * X = B                                   DGESL2     11
CD    using the factors computed by DGECO or DGEFA2.                    DGESL2     12
CD                                                                      DGESL2     13
CD    On Entry:                                                         DGESL2     14
CD                                                                      DGESL2     15
CD       A       DOUBLE PRECISION(LDA, N)                               DGESL2     16
CD               The output from DGECO or DGEFA2.                       DGESL2     17
CD                                                                      DGESL2     18
CD       LDA     INTEGER                                                DGESL2     19
CD               The leading dimension of the array  A.                 DGESL2     20
CD                                                                      DGESL2     21
CD       N       INTEGER                                                DGESL2     22
CD               The order of the matrix  A.                            DGESL2     23
CD                                                                      DGESL2     24
CD       IPVT    INTEGER(N)                                             DGESL2     25
CD               The pivot vector from DGECO or DGEFA2.                 DGESL2     26
CD                                                                      DGESL2     27
CD       B       DOUBLE PRECISION(N)                                    DGESL2     28
CD               The right hand side vector.                            DGESL2     29
CD                                                                      DGESL2     30
CD       JOB     INTEGER                                                DGESL2     31
CD               = 0         to solve  A * X = B ,                      DGESL2     32
CD               = NONZERO   to solve  TRANS(A) * X = B  where          DGESL2     33
CD                           TRANS(A)  is the transpose.                DGESL2     34
CD                                                                      DGESL2     35
CD    On Return:                                                        DGESL2     36
CD                                                                      DGESL2     37
CD       B       The solution vector  X.                                DGESL2     38
CD                                                                      DGESL2     39
CD    Error Condition:                                                  DGESL2     40
CD                                                                      DGESL2     41
CD       A division by zero will occur if the input factor contains a   DGESL2     42
CD       zero on the diagonal.  Technically this indicates singularity  DGESL2     43
CD       but it is often caused by improper arguments or improper       DGESL2     44
CD       setting of LDA.  It will not occur if the subroutines are      DGESL2     45
CD       called correctly and if DGECO has set RCOND .GT. 0.0           DGESL2     46
CD       or DGEFA2 has set INFO .EQ. 0 .                                DGESL2     47
CD                                                                      DGESL2     48
CD    To compute  INVERSE(A) * CD where  CD is a matrix                 DGESL2     49
CD    with  P  columns                                                  DGESL2     50
CD                                                                      DGESL2     51
CD          call DGECO (A, LDA, N, IPVT, RCOND, Z)                      DGESL2     52
CD          if (RCOND is too small) go to ...                           DGESL2     53
CD          do 10 J = 1, P                                              DGESL2     54
CD             call DGESL2 (A, LDA, N, IPVT, C(1,J), 0)                 DGESL2     55
CD       10 continue                                                    DGESL2     56
CD                                                                      DGESL2     57
CD    LINPACK. This version dated 08/14/78.                             DGESL2     58
CD    CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.      DGESL2     59
CD                                                                      DGESL2     60
CD -------------------------------------------------------------------  DGESL2     61
C                                                                       DGESL2     62
C                                                                       DGESL2     63
      DOUBLE PRECISION A(LDA,LDA), B(*), DDOT, T                        DGESL2     64
C                                                                       DGESL2     65
      DIMENSION          IPVT(*)                                        DGESL2     66
C                                                                       DGESL2     67
C                                                                       DGESL2     68
C---------------------------------------------------------------------  DGESL2     69
C                                                                       DGESL2     70
      NM1 = N - 1                                                       DGESL2     71
      IF (JOB .NE. 0) GO TO 50                                          DGESL2     72
C                                                                       DGESL2     73
C        JOB = 0 , solve  A * X = B                                     DGESL2     74
C        First solve  L * Y = B                                         DGESL2     75
C                                                                       DGESL2     76
         IF (NM1 .LT. 1) GO TO 30                                       DGESL2     77
         DO 20 K = 1, NM1                                               DGESL2     78
            L = IPVT(K)                                                 DGESL2     79
            T = B(L)                                                    DGESL2     80
            IF (L .EQ. K) GO TO 10                                      DGESL2     81
               B(L) = B(K)                                              DGESL2     82
               B(K) = T                                                 DGESL2     83
10          CONTINUE                                                    DGESL2     84
            CALL DAXPY (N-K, T, A(K+1,K), 1, B(K+1), 1)                 DGESL2     85
20       CONTINUE                                                       DGESL2     86
30       CONTINUE                                                       DGESL2     87
C                                                                       DGESL2     88
C        Now solve  U * X = Y.                                          DGESL2     89
C                                                                       DGESL2     90
         DO 40 KB = 1, N                                                DGESL2     91
            K = N + 1 - KB                                              DGESL2     92
            B(K) = B(K) / A(K,K)                                        DGESL2     93
            T = -B(K)                                                   DGESL2     94
            CALL DAXPY (K-1, T, A(1,K), 1, B(1), 1)                     DGESL2     95
40       CONTINUE                                                       DGESL2     96
      GO TO 100                                                         DGESL2     97
50    CONTINUE                                                          DGESL2     98
C                                                                       DGESL2     99
C        JOB = nonzero, solve  TRANS(A) * X = B                         DGESL2    100
C        First solve  TRANS(U) * Y = B                                  DGESL2    101
C                                                                       DGESL2    102
         DO 60 K = 1, N                                                 DGESL2    103
            T = DDOT (K-1, A(1,K), 1, B(1), 1)                          DGESL2    104
            B(K) = (B(K) - T) / A(K,K)                                  DGESL2    105
60       CONTINUE                                                       DGESL2    106
C                                                                       DGESL2    107
C        Now solve TRANS(L) * X = Y                                     DGESL2    108
C                                                                       DGESL2    109
         IF (NM1 .LT. 1) GO TO 90                                       DGESL2    110
         DO 80 KB = 1, NM1                                              DGESL2    111
            K = N - KB                                                  DGESL2    112
            B(K) = B(K) + DDOT (N-K, A(K+1,K), 1, B(K+1), 1)            DGESL2    113
            L = IPVT(K)                                                 DGESL2    114
            IF (L .EQ. K) GO TO 70                                      DGESL2    115
               T = B(L)                                                 DGESL2    116
               B(L) = B(K)                                              DGESL2    117
               B(K) = T                                                 DGESL2    118
70          CONTINUE                                                    DGESL2    119
80       CONTINUE                                                       DGESL2    120
90       CONTINUE                                                       DGESL2    121
100   CONTINUE                                                          DGESL2    122
      RETURN                                                            DGESL2    123
      END                                                               DGESL2    124
