CFPRIME LUDCM                                                           LUDCM       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 LUDCM(ISTATE)                                          LUDCM       3
CD                                                                      LUDCM       4
CD    TO CALCULATE THE L-U DECOMPSITION OF THE JACOBIAN MATRIX A FROM   LUDCM       5
CD    THE TRAPEZOIDAL METHOD.                                           LUDCM       6
CD                                                                      LUDCM       7
CD    ALGORITHM : SEE ALGORITHM 4.2-1, PAGE 57 OF MATRIX COMPUTATION    LUDCM       8
CD                BY G. H. GOLUB AND C. F. VAN LOAN                     LUDCM       9
CD                                                                      LUDCM      10
CD    NOTE      : THERE IS NO PIVOTING IN ELEMINATION IN ORDER TO KEEP  LUDCM      11
CD                LIFE SIMPLE (ASSUMING THAT THE DIAGONAL ELEMENTS ARE  LUDCM      12
CD                NON-ZERO IN ELEMINATION)                              LUDCM      13
CD                                                                      LUDCM      14
      INCLUDE 'PARX.INS'                                                LUDCM      15
      INCLUDE 'CNABAK.INS'                                              LUDCM      16
C                                                                       LUDCM      17
      DO 10 K=1,ISTATE-1                                                LUDCM      18
         DO 20 I=K+1,ISTATE                                             LUDCM      19
            IF (IA(I,K).EQ.0) GOTO 20                                   LUDCM      20
            A(I,K)=A(I,K)/A(K,K)                                        LUDCM      21
            DO 30 J=K+1,ISTATE                                          LUDCM      22
               A(I,J)=A(I,J)-A(I,K)*A(K,J)                              LUDCM      23
  30        CONTINUE                                                    LUDCM      24
  20     CONTINUE                                                       LUDCM      25
  10  CONTINUE                                                          LUDCM      26
C                                                                       LUDCM      27
      RETURN                                                            LUDCM      28
      END                                                               LUDCM      29
C                                                                       LUDCM      30
      SUBROUTINE FBSUB(NAIDEX,ISTATE,IRHS)                              LUDCM      31
CD                                                                      LUDCM      32
CD    TO PERFORM FORWARD AND BACKWARD SUBSTITUTIONS TO SOLVE THE LINEAR LUDCM      33
CD    EQUATION FROM THE TRAPEZOIDAL METHOD.                             LUDCM      34
CD                                                                      LUDCM      35
CD    ALGORITHM : STRAIGHTFORWARD                                       LUDCM      36
CD                                                                      LUDCM      37
CD    NOTE:                                                             LUDCM      38
CD                                                                      LUDCM      39
CD    1. THE L-U FACTORIZED MATRIX A IS STORED IN ARRAY AJCOB WITH      LUDCM      40
CD       OFFSET VALUE IAIDEX. ON RETURN, IAIDEX IS ADVANCED TO THE      LUDCM      41
CD       OFFSET VALUE OF THE NEXT DEVICE.                               LUDCM      42
CD                                                                      LUDCM      43
CD    2. THE RIGHT HAND SIDE VECTOR IS IN ARRAY XRES WITH OFFSET VALUE  LUDCM      44
CD       IRHS. ON RETURN, XRES CONTAINS THE SOLUTION VECTOR.            LUDCM      45
CD                                                                      LUDCM      46
      INCLUDE 'PARX.INS'                                                LUDCM      47
      INCLUDE 'CNSVS.INS'                                               LUDCM      48
      INCLUDE 'CNJCOB.INS'                                              LUDCM      49
      INCLUDE 'CNXDOT.INS'                                              LUDCM      50
      INCLUDE 'CNABAK.INS'                                              LUDCM      51
C                                                                       LUDCM      52
      REAL WV1(NSID)                                                    LUDCM      53
C                                                                       LUDCM      54
C     FORWARD SUBSTITUTION                                              LUDCM      55
C                                                                       LUDCM      56
      DO 10 I=2,ISTATE                                                  LUDCM      57
         IWS1=IRHS+I                                                    LUDCM      58
         DO 20 J=1,I-1                                                  LUDCM      59
            XRES(IWS1)=XRES(IWS1)-AJCOB(NAIDEX+J)*XRES(IRHS+J)          LUDCM      60
  20     CONTINUE                                                       LUDCM      61
         NAIDEX=NAIDEX+ISTATE                                           LUDCM      62
  10  CONTINUE                                                          LUDCM      63
C                                                                       LUDCM      64
      IAIDEX=NAIDEX                                                     LUDCM      65
C                                                                       LUDCM      66
C     BACKWARD SUBSTITUTION                                             LUDCM      67
C                                                                       LUDCM      68
      XRES(IRHS+ISTATE)=XRES(IRHS+ISTATE)/AJCOB(NAIDEX)                 LUDCM      69
      NAIDEX=NAIDEX-ISTATE                                              LUDCM      70
      DO 30 I=ISTATE-1,1,-1                                             LUDCM      71
         IWS1=IRHS+I                                                    LUDCM      72
         IWS2=NAIDEX-ISTATE                                             LUDCM      73
         DO 40 J=ISTATE,I+1,-1                                          LUDCM      74
            XRES(IWS1)=XRES(IWS1)-AJCOB(IWS2+J)*XRES(IRHS+J)            LUDCM      75
  40     CONTINUE                                                       LUDCM      76
         XRES(IWS1)=XRES(IWS1)/AJCOB(IWS2+I)                            LUDCM      77
         NAIDEX=IWS2                                                    LUDCM      78
  30  CONTINUE                                                          LUDCM      79
C                                                                       LUDCM      80
      RETURN                                                            LUDCM      81
C                                                                       LUDCM      82
      ENTRY FBSUBA(ISTATE)                                              LUDCM      83
CD                                                                      LUDCM      84
CD    SIMILAR TO FBSUB EXCEPT THAT INVERSE(A)*B IS SOLVED AND RESULTS   LUDCM      85
CD    ARE STILL STORED IN B ARRAY.                                      LUDCM      86
CD                                                                      LUDCM      87
      DO 50 K=1,ICOLBV*2                                                LUDCM      88
C                                                                       LUDCM      89
C     FORWARD SUBSTITUTION                                              LUDCM      90
C                                                                       LUDCM      91
         DO 60 I=2,ISTATE                                               LUDCM      92
            DO 70 J=1,I-1                                               LUDCM      93
               B(I,K)=B(I,K)-A(I,J)*B(J,K)                              LUDCM      94
  70        CONTINUE                                                    LUDCM      95
  60     CONTINUE                                                       LUDCM      96
C                                                                       LUDCM      97
C     BACKWARD SUBSTITUTION                                             LUDCM      98
C                                                                       LUDCM      99
         B(ISTATE,K)=B(ISTATE,K)/A(ISTATE,ISTATE)                       LUDCM     100
         DO 80 I=ISTATE-1,1,-1                                          LUDCM     101
            DO 90 J=ISTATE,I+1,-1                                       LUDCM     102
               B(I,K)=B(I,K)-A(I,J)*B(J,K)                              LUDCM     103
  90        CONTINUE                                                    LUDCM     104
            B(I,K)=B(I,K)/A(I,I)                                        LUDCM     105
  80     CONTINUE                                                       LUDCM     106
C                                                                       LUDCM     107
  50  CONTINUE                                                          LUDCM     108
C                                                                       LUDCM     109
      RETURN                                                            LUDCM     110
C                                                                       LUDCM     111
      ENTRY FBSUBB(IOFSET,ISTATE)                                       LUDCM     112
CD                                                                      LUDCM     113
CD    SIMILAR TO FBSUBA EXCEPT THAT MATRIX A IS IN A DIFFERENT          LUDCM     114
CD    ORDERING.                                                         LUDCM     115
CD                                                                      LUDCM     116
      DO 100 K=1,ICOLBV*2                                               LUDCM     117
         DO 105 I=1,ISTATE                                              LUDCM     118
            WV1(I)=B(JCBPEM(IOFSET+I),K)                                LUDCM     119
 105     CONTINUE                                                       LUDCM     120
C                                                                       LUDCM     121
C     FORWARD SUBSTITUTION                                              LUDCM     122
C                                                                       LUDCM     123
         DO 110 I=2,ISTATE                                              LUDCM     124
            DO 120 J=1,I-1                                              LUDCM     125
               IF (IA(I,J).NE.0) WV1(I)=WV1(I)-A(I,J)*WV1(J)            LUDCM     126
 120        CONTINUE                                                    LUDCM     127
 110     CONTINUE                                                       LUDCM     128
C                                                                       LUDCM     129
C     BACKWARD SUBSTITUTION                                             LUDCM     130
C                                                                       LUDCM     131
         WV1(ISTATE)=WV1(ISTATE)/A(ISTATE,ISTATE)                       LUDCM     132
         DO 130 I=ISTATE-1,1,-1                                         LUDCM     133
            DO 140 J=ISTATE,I+1,-1                                      LUDCM     134
               IF (IA(I,J).NE.0) WV1(I)=WV1(I)-A(I,J)*WV1(J)            LUDCM     135
 140        CONTINUE                                                    LUDCM     136
            WV1(I)=WV1(I)/A(I,I)                                        LUDCM     137
 130     CONTINUE                                                       LUDCM     138
C                                                                       LUDCM     139
         DO 150 I=1,ISTATE                                              LUDCM     140
            B(JCBPEM(IOFSET+I),K)=WV1(I)                                LUDCM     141
 150     CONTINUE                                                       LUDCM     142
 100  CONTINUE                                                          LUDCM     143
C                                                                       LUDCM     144
      RETURN                                                            LUDCM     145
      END                                                               LUDCM     146
C                                                                       LUDCM     147
      SUBROUTINE FBSUBS(ISTATE,IRHS)                                    LUDCM     148
CD                                                                      LUDCM     149
CD    SIMILAR TO SUBROUTINE FBSUB EXCEPT THAT THE SPARSE A STRUCTURE    LUDCM     150
CD    IS USED TO SOLVE THE EQUATIONS. THIS SUBROUTINE IS FOR SOLVING    LUDCM     151
CD    DETAILED GENERATOR EQUATIONS ONLY.                                LUDCM     152
CD                                                                      LUDCM     153
CD    COMPUTATION STEPS:                                                LUDCM     154
CD                                                                      LUDCM     155
      INCLUDE 'PARX.INS'                                                LUDCM     156
      INCLUDE 'CNJCOB.INS'                                              LUDCM     157
      INCLUDE 'CNXDOT.INS'                                              LUDCM     158
      INCLUDE 'CNABAK.INS'                                              LUDCM     159
C                                                                       LUDCM     160
      REAL WV1(NSID)                                                    LUDCM     161
C                                                                       LUDCM     162
      DO 10 I=1,ISTATE                                                  LUDCM     163
         IWS1=IRHS+JCBPEM(IRHS+I)                                       LUDCM     164
         WV1(I)=XRES(IWS1)                                              LUDCM     165
  10  CONTINUE                                                          LUDCM     166
CD                                                                      LUDCM     167
CD    STEP 1 PERFORM FORWARD SUBSTITUTION.                              LUDCM     168
CD                                                                      LUDCM     169
      NAIDEX=IAIDEX+2                                                   LUDCM     170
      NALUFT=IAJCOB(1,IAIDEX+1)                                         LUDCM     171
      DO 100 I=2,ISTATE                                                 LUDCM     172
         NROW=IAJCOB(1,NAIDEX+I)                                        LUDCM     173
         IF (NROW.EQ.0) GOTO 100                                        LUDCM     174
         DO 110 J=1,NROW                                                LUDCM     175
            NALUFT=NALUFT+1                                             LUDCM     176
            NCOL=IALUFT(NALUFT)                                         LUDCM     177
            WV1(I)=WV1(I)-ALUFT(NALUFT)*WV1(NCOL)                       LUDCM     178
  110    CONTINUE                                                       LUDCM     179
  100 CONTINUE                                                          LUDCM     180
CD                                                                      LUDCM     181
CD    STEP 2 PERFORM BACKWARD SUBSTITUTION.                             LUDCM     182
CD                                                                      LUDCM     183
      WV1(ISTATE)=WV1(ISTATE)*AJCOB(NAIDEX+ISTATE)                      LUDCM     184
      DO 120 I=ISTATE-1,1,-1                                            LUDCM     185
         NROW=IAJCOB(2,NAIDEX+I)                                        LUDCM     186
         IF (NROW.EQ.0) GOTO 140                                        LUDCM     187
         DO 130 J=1,NROW                                                LUDCM     188
            NALUFT=NALUFT+1                                             LUDCM     189
            NCOL=IALUFT(NALUFT)                                         LUDCM     190
            WV1(I)=WV1(I)-ALUFT(NALUFT)*WV1(NCOL)                       LUDCM     191
  130    CONTINUE                                                       LUDCM     192
  140    WV1(I)=WV1(I)*AJCOB(NAIDEX+I)                                  LUDCM     193
  120 CONTINUE                                                          LUDCM     194
C                                                                       LUDCM     195
      IAIDEX=NAIDEX+ISTATE                                              LUDCM     196
C                                                                       LUDCM     197
      DO 160 I=1,ISTATE                                                 LUDCM     198
         IWS1=IRHS+JCBPEM(IRHS+I)                                       LUDCM     199
         XRES(IWS1)=WV1(I)                                              LUDCM     200
 160  CONTINUE                                                          LUDCM     201
C                                                                       LUDCM     202
      RETURN                                                            LUDCM     203
      END                                                               LUDCM     204
