CFPRIME PRONYX                                                          PRONYX      1
C***********************************************************************RIGHTE      1
C*                                                                     *RIGHTE      2
C*               Output Analysis Program for ETMSP (OAP)               *RIGHTE      4
C*               =======================================               *RIGHTE      5
C*                                                                     *RIGHTE     11
C*        copyright @(1993) Electric Power Research Institute          *RIGHTE     12
C*                                                                     *RIGHTE     13
C* EPRI reserves all rights in the software and associated             *RIGHTE     14
C* documentation as delivered. The software and its associated         *RIGHTE     15
C* documention or any portion thereof may not be reproduced in any form*RIGHTE     16
C* whatsoever without the consent of EPRI, except as explicity provided*RIGHTE     17
C* for by license. Licensee shall submit a detailed description to EPRI*RIGHTE     18
C* of any changes or modifications to the software within two (2)      *RIGHTE     19
C* months after such changes or modifications are made. At EPRI's sole *RIGHTE     20
C* discretion, the software, as changed or modified, shall be given a  *RIGHTE     21
C* new designation sufficiently different from its then current        *RIGHTE     22
C* designation as to prevent mistake,confusion, or deception as between*RIGHTE     23
C* the current software as changed or modified.                        *RIGHTE     24
C*                                                                     *RIGHTE     25
C* A license under EPRI's rights in the software can be obtained       *RIGHTE     26
C* directly from EPRI.                                                 *RIGHTE     27
C*                                                                     *RIGHTE     28
C* This program and supporting materials may be independently available*RIGHTE     29
C* from Electric Power Software centre (EPSC) for an appropriate       *RIGHTE     30
C* distribution fee. The EPSC is a service of                          *RIGHTE     31
C*                                                                     *RIGHTE     32
C*                    Power Computing Company                          *RIGHTE     33
C*                    1930 Hi Line                                     *RIGHTE     34
C*                    Dallas, Tx 75207                                 *RIGHTE     35
C*                                                                     *RIGHTE     36
C* NEITHER EPRI, ANY MEMBER OF EPRI, NOR ANY PERSON OR  ORGANIZATION   *RIGHTE     37
C* ACTING ON BEHAVE OF THEM:                                           *RIGHTE     38
C*                                                                     *RIGHTE     39
C* 1. MAKE ANY WARRANTY OR REPRESENTATION WHATSOEVER, EXPRESS OR       *RIGHTE     40
C*    IMPLIED, INCLUDING ANY WARRANTY OF MERCHANTABILITY OR FITNESS OR *RIGHTE     41
C*    ANY PURPOSE WITH RESPECT TO THE SOFTWARE OR ITS ASSOCIATED       *RIGHTE     42
C*    DOCUMENTATION; OR                                                *RIGHTE     43
C*                                                                     *RIGHTE     44
C* 2. ASSUMES ANY LIABILITY WHATSOEVER WITH RESPECT TO ANY USE OF THE  *RIGHTE     45
C*    SOFTWARE OR ITS ASSOCIATED DOCUMENTATION OR ANY PORTION THEREOF  *RIGHTE     46
C*    OR WITH RESPECT TO ANY DAMAGES WHICH MAY RESULT FROM SUCH USE.   *RIGHTE     47
C*                                                                     *RIGHTE     48
C* THIS NOTICE MAY NOT BE REMOVED FROM THE SOFTWARE BY ANY USER        *RIGHTE     49
C* THEREOF.                                                            *RIGHTE     50
C*                                                                     *RIGHTE     51
C***********************************************************************RIGHTE     52
      SUBROUTINE PRONYX(NFSAVE,TSTEP,NOGO,ANS,LU15,DBP,                 PRONYX      3
     1 FREQ,RESP,ANGLE,XSIG,SHAT)                                       PRONYX      4
CD                                                                      PRONYX      5
CD       THIS ROUTINE IS THE MAIN DRIVER FOR THE SPECTRAL ANALYSIS.     PRONYX      6
CD       IT FUNCTIONS LISTED AS FOLLOWS :                               PRONYX      7
CD                                                                      PRONYX      8
CD           1) SET MATHEMATICAL CONSTANTS                              PRONYX      9
CD           2) SET LOGICAL UNIT POINTERS                               PRONYX     10
CD           3) INITIALIZE OTHER CONTROLS                               PRONYX     11
CD           4) INTIALIZE COUNTERS AND DEFAULT CONTROLS                 PRONYX     12
CD           5) ENTER PRONY LOGIC                                       PRONYX     13
CD           6) CALCULATE THE FREQUENCY RESPONSE                        PRONYX     14
CD                                                                      PRONYX     15
CD       CALLING PARAMETERS :                                           PRONYX     16
CD                                                                      PRONYX     17
CD           NFSAVE - NUMBER OF DATA POINT                              PRONYX     18
CD           TSTEP  - TIME INTERVAL                                     PRONYX     19
CD           NOGO   - I/O ERROR FLAG                                    PRONYX     20
CD           ANS    - MODE TRIM VARIABLE                                PRONYX     21
CD           LU15   - LOGICAL UNIT NUMBER FOR PRONY ANALYSIS OUTPUT     PRONYX     22
CD           DBP    - LOGICAL FLAG FOR DEBUG                            PRONYX     23
CD           FREQ   - FREQUENCY ARRAY                                   PRONYX     24
CD           RESP   - MAGNITUDE ARRAY (FUNCTION OF FREQUENCY)           PRONYX     25
CD           ANGLE  - ANGLE ARRAY (FUNCTION OF FREQUENCY)               PRONYX     26
CD           XSIG   - ORIGINAL Y-AXIS ARRAY                             PRONYX     27
CD           SHAT   - CALCULATED Y-AXIS ARRAY AFTER PRONY ANALYSIS      PRONYX     28
CD                                                                      PRONYX     29
      INCLUDE 'SOAPC1.INS'                                              PRONYX     30
      INTEGER NSIZE ,NTAIL                                              PRONYX     31
      PARAMETER (NSIZE=128)                                             PRONYX     32
      PARAMETER (NTAIL=1024)                                            PRONYX     33
      PARAMETER (NSTEP=100)                                             PRONYX     34
      DOUBLE PRECISION DAMP(1:NSIZE),FRQHZ(1:NSIZE),RELMAG(1:NSIZE),    PRONYX     35
     1                 PHASE(1:NSIZE),AMPOPT(1:NSIZE),PHIOPT(1:NSIZE)   PRONYX     36
      DOUBLE PRECISION FTRIMH,FTRIML                                    PRONYX     37
      INTEGER LU5,LU6,LU8,LU15,MRET                                     PRONYX     38
      INTEGER IZERO,IONE                                                PRONYX     39
      INTEGER NFSAVE,NSTART,LPOORG,LPMORG,FBSAVE,PIRORG,QSAVE           PRONYX     40
      INTEGER NNFIT,NNHAT,QCON                                          PRONYX     41
      DOUBLE PRECISION XSIG(0:NTAIL-1),SHAT(0:NTAIL-1)                  PRONYX     42
      DOUBLE PRECISION TRMLPM,SNR                                       PRONYX     43
      REAL TSTEP,ANS                                                    PRONYX     44
      DIMENSION FREQ(NSTEP+2),RESP(NSTEP+2),ANGLE(NSTEP+2)              PRONYX     45
      COMPLEX K1,K2,RES(NSTEP)                                          PRONYX     46
                                                                        PRONYX     47
      LOGICAL DBP                                                       PRONYX     48
      INTEGER NOGO                                                      PRONYX     49
C      CHARACTER*80 MSG,STRING                                          PRONYX     50
      NOGO=0                                                            PRONYX     51
C                                                                       PRONYX     52
C *** SET MATHEMATICAL CONSTANTS                                        PRONYX     53
      IZERO=0                                                           PRONYX     54
      IONE=1                                                            PRONYX     55
      TWOPI=8.0*ATAN(1.0)                                               PRONYX     56
      DEG=360.0/TWOPI                                                   PRONYX     57
      DEGI=TWOPI/360.0                                                  PRONYX     58
C                                                                       PRONYX     59
C *** SET LOGICAL UNIT POINTERS                                         PRONYX     60
      LU5=LIN                                                           PRONYX     61
      LU6=LOUT                                                          PRONYX     62
      LU8=8                                                             PRONYX     63
C                                                                       PRONYX     64
C *** INITIALIZE OTHER CONTROLS                                         PRONYX     65
      TRMLPM=DBLE(ANS/100.0)                                            PRONYX     66
C                                                                       PRONYX     67
C *** INTIALIZE COUNTERS AND DEFAULT CONTROLS                           PRONYX     68
      NNFIT=0                                                           PRONYX     69
      NSTART=0                                                          PRONYX     70
      LPOORG=0                                                          PRONYX     71
      LPMORG=3                                                          PRONYX     72
      FBSAVE=1                                                          PRONYX     73
      PIRORG=0                                                          PRONYX     74
      QSAVE=0                                                           PRONYX     75
C                                                                       PRONYX     76
C *** ENTER PRONY LOGIC                                                 PRONYX     77
      NNHAT=NNFIT                                                       PRONYX     78
      FTRIMH=0.0                                                        PRONYX     79
      FTRIML=0.0                                                        PRONYX     80
      CALL PRSPAK (XSIG,TSTEP,NFSAVE,NSTART                             PRONYX     81
     *,QCON,SNR,SHAT,DAMP,FRQHZ,RELMAG,PHASE,AMPOPT,PHIOPT              PRONYX     82
     *,LPOORG,LPMORG,FBSAVE,PIRORG,QSAVE,NNHAT                          PRONYX     83
     *,TRMLPM,FTRIMH,FTRIML                                             PRONYX     84
     *,LU5,LU6,LU15,MRET,DBP,NOGO)                                      PRONYX     85
                                                                        PRONYX     86
      IF(NOGO.LT.0) GOTO 9999                                           PRONYX     87
                                                                        PRONYX     88
C *** CALCULATE THE FREQUENCY RESPONSE                                  PRONYX     89
      FRQNYQ=1.0/(2.0*TSTEP)                                            PRONYX     90
      FRQ=FRQNYQ/REAL(NSTEP)                                            PRONYX     91
      DO 10 N=1,NSTEP                                                   PRONYX     92
        FREQ(N)=REAL(N)*FRQ                                             PRONYX     93
        RES(N)=CMPLX(0.,0.)                                             PRONYX     94
 10   CONTINUE                                                          PRONYX     95
C                                                                       PRONYX     96
      DO 30 K=1,QCON                                                    PRONYX     97
        AMPK=0.5*AMPOPT(K)                                              PRONYX     98
        PHIOPK=PHIOPT(K)*DEGI                                           PRONYX     99
        K1=CMPLX(AMPK*COS(PHIOPK),AMPK*SIN(PHIOPK))                     PRONYX    100
        K2=CMPLX(AMPK*COS(PHIOPK),-AMPK*SIN(PHIOPK))                    PRONYX    101
        DO 20 N=1,NSTEP                                                 PRONYX    102
          RES(N)=RES(N)+K1/CMPLX(-DAMP(K),(FREQ(N)-FRQHZ(K))*TWOPI) +   PRONYX    103
     +                  K2/CMPLX(-DAMP(K),(FREQ(N)+FRQHZ(K))*TWOPI)     PRONYX    104
 20     CONTINUE                                                        PRONYX    105
 30   CONTINUE                                                          PRONYX    106
      DO 40 N=1,NSTEP                                                   PRONYX    107
        RESP(N)=CABS(RES(N))                                            PRONYX    108
        RESNI=AIMAG(RES(N))                                             PRONYX    109
        RESNR=REAL(RES(N))                                              PRONYX    110
        ANGLE(N)=ATAN2(RESNI,RESNR) * DEG                               PRONYX    111
C        WRITE (LUT,'(3E21.10)') FREQ(N),RESP(N),ANGLE(N)               PRONYX    112
 40   CONTINUE                                                          PRONYX    113
C                                                                       PRONYX    114
 9999 CONTINUE                                                          PRONYX    115
C                                                                       PRONYX    116
      RETURN                                                            PRONYX    117
C                                                                       PRONYX    118
      END                                                               PRONYX    119
C                                                                       PRONYX    120
C                                                                       PRONYX    121
        SUBROUTINE COMBIN (N,M,IC)                                      PRONYX    122
CD                                                                      PRONYX    123
CD         This routine performs complementing two bits to get          PRONYX    124
CD         the next combination                                         PRONYX    125
CD                                                                      PRONYX    126
CD         CALLING PARAMETERS :                                         PRONYX    127
CD                                                                      PRONYX    128
CD                         N  -  ACTUAL DIMENSION OF THE VECTOR         PRONYX    129
CD                         M  -  MAX DIMENSION OF THE VECTOR            PRONYX    130
CD                         IC -  INTEGER ARRAY                          PRONYX    131
CD                                                                      PRONYX    132
        INTEGER IC(1:*) ,N ,N1 ,NJ ,M ,I ,I1 ,J ,J1 ,JP ,K1 ,K2         PRONYX    133
                                                                        PRONYX    134
C Check ending pattern of vector                                        PRONYX    135
        IF (M.GE.N .OR. M.EQ.0) GOTO 140                                PRONYX    136
        N1 = N - 1                                                      PRONYX    137
        DO 10 J=1,N1                                                    PRONYX    138
          NJ = N - J                                                    PRONYX    139
          IF (IC(N).EQ.IC(NJ)) GOTO 10                                  PRONYX    140
          J1 = J                                                        PRONYX    141
          GOTO 20                                                       PRONYX    142
10      CONTINUE                                                        PRONYX    143
20      IF (MOD(M,2).EQ.1) GOTO 90                                      PRONYX    144
C For M even                                                            PRONYX    145
        IF (IC(N).EQ.1) GOTO 30                                         PRONYX    146
        K1 = N - J1                                                     PRONYX    147
        K2 = K1 + 1                                                     PRONYX    148
        GOTO 130                                                        PRONYX    149
30      IF (MOD(J1,2).EQ.1) GOTO 40                                     PRONYX    150
        GOTO 120                                                        PRONYX    151
C Scan from right to left                                               PRONYX    152
40      JP = (N-J1) - 1                                                 PRONYX    153
        DO 50 I=1,JP                                                    PRONYX    154
          I1 = JP + 2 - I                                               PRONYX    155
          IF (IC(I1).EQ.0) GOTO 50                                      PRONYX    156
          IF (IC(I1-1).EQ.1) GOTO 70                                    PRONYX    157
          GOTO 80                                                       PRONYX    158
50      CONTINUE                                                        PRONYX    159
60      K1 = 1                                                          PRONYX    160
        K2 = (N+1) - M                                                  PRONYX    161
        GOTO 130                                                        PRONYX    162
70      K1 = I1 - 1                                                     PRONYX    163
        K2 = N - J1                                                     PRONYX    164
        GOTO 130                                                        PRONYX    165
80      K1 = I1 - 1                                                     PRONYX    166
        K2 = (N+1) - J1                                                 PRONYX    167
        GOTO 130                                                        PRONYX    168
C For M odd                                                             PRONYX    169
90      IF (IC(N).EQ.1) GOTO 110                                        PRONYX    170
        K2 = (N-J1) - 1                                                 PRONYX    171
        IF (K2.EQ.0) GOTO 60                                            PRONYX    172
        IF (IC(K2+1).EQ.1 .AND. IC(K2).EQ.1) GOTO 100                   PRONYX    173
        K1 = K2 + 1                                                     PRONYX    174
        GOTO 130                                                        PRONYX    175
100     K1 = N                                                          PRONYX    176
        GOTO 130                                                        PRONYX    177
110     IF (MOD(J1,2).EQ.1) GOTO 120                                    PRONYX    178
        GOTO 40                                                         PRONYX    179
120     K1 = N - J1                                                     PRONYX    180
        K2 = MIN0((K1+2),N)                                             PRONYX    181
C Complementing two bits to get the next combination                    PRONYX    182
130     IC(K1) = 1 - IC(K1)                                             PRONYX    183
        IC(K2) = 1 - IC(K2)                                             PRONYX    184
140     RETURN                                                          PRONYX    185
        END                                                             PRONYX    186
C                                                                       PRONYX    187
        SUBROUTINE DSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)       PRONYX    188
CD                                                                      PRONYX    189
CD         This routine reduces a double precision NXP matrix X by      PRONYX    190
CD         orthronal transformations U and V to diag form. The diag     PRONYX    191
CD         elements S(I) are the singular values of X. The col of U     PRONYX    192
CD         are the corresponding to the left singular vectors, and      PRONYX    193
Cd         the col. of V the right singular vectors                     PRONYX    194
CD                                                                      PRONYX    195
CD  PARAMETERS :                                                        PRONYX    196
CD                                                                      PRONYX    197
CD   ON ENTRY :                                                         PRONYX    198
CD       X  - DOUBLE PREC(LDX,P), where LDX.GE.N. The matrix whose      PRONYX    199
CD            singular value decomp. is to be computed. X is destroyed  PRONYX    200
CD            by DSVDC.                                                 PRONYX    201
CD       LDX- INTEGER, The leading dimension of the array X.            PRONYX    202
CD       N  - INTEGER, The number of columns of the matrix X.           PRONYX    203
CD       P  - INTEGER, The number of rows of the matrix X.              PRONYX    204
CD       LDU- INTEGER, The leading dimen of the array U. (see below).   PRONYX    205
CD       LDV- INTEGER, The leading dimen of the array V. (see below).   PRONYX    206
CD      WORK- DOUBLE PRECISION(N), WORK is a scratch array.             PRONYX    207
CD       JOB- INTEGER, JOB controls the comp of the singular vectors.   PRONYX    208
CD            It has the decimal expansion AB with the  meaning below   PRONYX    209
CD            A.EQ.0  Do not compute the left singular vectors.         PRONYX    210
CD            A.EQ.1  Return the N left singular vectors in U.          PRONYX    211
CD            A.GE.2  Return the first MIN(N,P) singular vectors in U.  PRONYX    212
CD            B.EQ.0  Do not compute the right singular vectors.        PRONYX    213
CD            B.EQ.1  Return the right singular vectors in V.           PRONYX    214
CD ON RETURN :                                                          PRONYX    215
CD       S  - DOUBLE PREC(MM), where MM=MIN(N+1,P).The first MIN(N,P).  PRONYX    216
CD            S contain the singular values of X arranged in,           PRONYX    217
CD            descending order of magnitude.                            PRONYX    218
CD       E  - DOUBLE PREC(P). E ordinarily contains zeros.However see.  PRONYX    219
CD            the discussion of INFO for exceptions.                    PRONYX    220
CD       U  - DOUBLE PREC(LDU,K), where LDU.GE.N. If JOBA.EQ.1 then     PRONYX    221
CD            K.EQ.N,If JOBA.GE.2 then K.EQ.MIN(N,P).U IS the matrix    PRONYX    222
CD            of right singular vectors. U is not ref if JOBA.EQ.0. If  PRONYX    223
CD            N.LE.P or if JOBA=.2 ,then U may be ident. with X in the  PRONYX    224
CD            subroutine call.                                          PRONYX    225
CD       V  - DOUBLE PREC(LDV,P), where LDV.GE.P.V contains the matrix  PRONYX    226
CD            of right singular vectors.V is not ref if JOB.EQ.0. Iff   PRONYX    227
CD            P.LE.N,then V may be ident with X in the subroutine call  PRONYX    228
CD     INFO - INTEGER, The singular values (and their corresp singular  PRONYX    229
CD            vectors) S(INFO+1),S(INFO+2),...,S(M) are correct (here   PRONYX    230
CD            M=MIN(N,P)).Thus if INFO=0, all the singular values andl  PRONYX    231
CD            their vectors are correct. In any event, the matrix       PRONYX    232
CD            B=TRANS(U)*X*V is the bidiag matrix with the elements of  PRONYX    233
CD            S on its diag and the elements of E on its super-diag.    PRONYX    234
CD            Thus the singular values of X and B are the same.         PRONYX    235
CD                                                                      PRONYX    236
CD LINPACK. THIS VERSION DATED 03/19/79 .                               PRONYX    237
CD G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.          PRONYX    238
CD                                                                      PRONYX    239
CD SUBROUTINES AND FUNCTIONS : EXTERNAL DROT                            PRONYX    240
CD                             BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2,DROTG  PRONYX    241
CD                             FORTRAN DABS,DMAX1,MAX0,MIN0,MOD,DSQRT   PRONYX    242
C Argument variables                                                    PRONYX    243
        INTEGER LDX,N,P,LDU,LDV,JOB,INFO                                PRONYX    244
        DOUBLE PRECISION X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1)   PRONYX    245
C Internal variables                                                    PRONYX    246
        INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT, PRONYX    247
     *        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1                        PRONYX    248
        DOUBLE PRECISION DDOT,T,R                                       PRONYX    249
        DOUBLE PRECISION B,C,CS,EL,EMM1,F,G,DNRM2,SCALE,SHIFT,SL,SM,SN, PRONYX    250
     *                 SMM1,T1,TEST,ZTEST                               PRONYX    251
        LOGICAL WANTU,WANTV                                             PRONYX    252
C                                                                       PRONYX    253
C set the maximum number of iterations.                                 PRONYX    254
C                                                                       PRONYX    255
        MAXIT = 30                                                      PRONYX    256
C                                                                       PRONYX    257
C determine what is to be computed.                                     PRONYX    258
C                                                                       PRONYX    259
        WANTU = .FALSE.                                                 PRONYX    260
        WANTV = .FALSE.                                                 PRONYX    261
        JOBU = MOD(JOB,100)/10                                          PRONYX    262
        NCU = N                                                         PRONYX    263
        IF (JOBU .GT. 1) NCU = MIN0(N,P)                                PRONYX    264
        IF (JOBU .NE. 0) WANTU = .TRUE.                                 PRONYX    265
        IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.                          PRONYX    266
C                                                                       PRONYX    267
C reduce X to bidiagonal form, storing the diagonal elements in S and thPRONYX    268
C super-diagonal elements in E.                                         PRONYX    269
C                                                                       PRONYX    270
        INFO = 0                                                        PRONYX    271
        NCT = MIN0(N-1,P)                                               PRONYX    272
        NRT = MAX0(0,MIN0(P-2,N))                                       PRONYX    273
        LU = MAX0(NCT,NRT)                                              PRONYX    274
        IF (LU .LT. 1) GO TO 170                                        PRONYX    275
        DO 160 L = 1, LU                                                PRONYX    276
           LP1 = L + 1                                                  PRONYX    277
           IF (L .GT. NCT) GO TO 20                                     PRONYX    278
C                                                                       PRONYX    279
C compute the transformation for the L-TH column and place the L-TH diagPRONYX    280
C S(L).                                                                 PRONYX    281
C                                                                       PRONYX    282
              S(L) = DNRM2(N-L+1,X(L,L),1)                              PRONYX    283
              IF (S(L) .EQ. 0.0D0) GO TO 10                             PRONYX    284
                 IF (X(L,L) .NE. 0.0D0) S(L) = DSIGN(S(L),X(L,L))       PRONYX    285
                 CALL DSCAL(N-L+1,1.0D0/S(L),X(L,L),1)                  PRONYX    286
                 X(L,L) = 1.0D0 + X(L,L)                                PRONYX    287
10            CONTINUE                                                  PRONYX    288
              S(L) = -S(L)                                              PRONYX    289
20         CONTINUE                                                     PRONYX    290
           IF (P .LT. LP1) GO TO 50                                     PRONYX    291
           DO 40 J = LP1, P                                             PRONYX    292
              IF (L .GT. NCT) GO TO 30                                  PRONYX    293
              IF (S(L) .EQ. 0.0D0) GO TO 30                             PRONYX    294
C                                                                       PRONYX    295
C apply the transformation.                                             PRONYX    296
C                                                                       PRONYX    297
                 T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)              PRONYX    298
                 CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1)                  PRONYX    299
30            CONTINUE                                                  PRONYX    300
C                                                                       PRONYX    301
C place the L-TH row of X into E for the subsequent calculation of the rPRONYX    302
C transformation.                                                       PRONYX    303
C                                                                       PRONYX    304
              E(J) = X(L,J)                                             PRONYX    305
40         CONTINUE                                                     PRONYX    306
50         CONTINUE                                                     PRONYX    307
           IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70                     PRONYX    308
C                                                                       PRONYX    309
C place the transformation in U for subsequent back multiplication.     PRONYX    310
C                                                                       PRONYX    311
              DO 60 I = L, N                                            PRONYX    312
                 U(I,L) = X(I,L)                                        PRONYX    313
60            CONTINUE                                                  PRONYX    314
70         CONTINUE                                                     PRONYX    315
           IF (L .GT. NRT) GO TO 150                                    PRONYX    316
C                                                                       PRONYX    317
C compute the L-TH row transformation and place the L-TH super-diagonal PRONYX    318
C                                                                       PRONYX    319
              E(L) = DNRM2(P-L,E(LP1),1)                                PRONYX    320
              IF (E(L) .EQ. 0.0D0) GO TO 80                             PRONYX    321
                 IF (E(LP1) .NE. 0.0D0) E(L) = DSIGN(E(L),E(LP1))       PRONYX    322
                 CALL DSCAL(P-L,1.0D0/E(L),E(LP1),1)                    PRONYX    323
                 E(LP1) = 1.0D0 + E(LP1)                                PRONYX    324
80            CONTINUE                                                  PRONYX    325
              E(L) = -E(L)                                              PRONYX    326
              IF (LP1 .GT. N .OR. E(L) .EQ. 0.0D0) GO TO 120            PRONYX    327
C                                                                       PRONYX    328
C apply the transformation.                                             PRONYX    329
C                                                                       PRONYX    330
                 DO 90 I = LP1, N                                       PRONYX    331
                    WORK(I) = 0.0D0                                     PRONYX    332
90               CONTINUE                                               PRONYX    333
                 DO 100 J = LP1, P                                      PRONYX    334
                    CALL DAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1)         PRONYX    335
100              CONTINUE                                               PRONYX    336
                 DO 110 J = LP1, P                                      PRONYX    337
                    CALL DAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1) PRONYX    338
110              CONTINUE                                               PRONYX    339
120           CONTINUE                                                  PRONYX    340
              IF (.NOT.WANTV) GO TO 140                                 PRONYX    341
C                                                                       PRONYX    342
C place the transformation in V for subsequent back multiplication.     PRONYX    343
C                                                                       PRONYX    344
                 DO 130 I = LP1, P                                      PRONYX    345
                    V(I,L) = E(I)                                       PRONYX    346
130              CONTINUE                                               PRONYX    347
140           CONTINUE                                                  PRONYX    348
150        CONTINUE                                                     PRONYX    349
160     CONTINUE                                                        PRONYX    350
170     CONTINUE                                                        PRONYX    351
C                                                                       PRONYX    352
C set up the final bidiagonal matrix of order M.                        PRONYX    353
C                                                                       PRONYX    354
        M = MIN0(P,N+1)                                                 PRONYX    355
        NCTP1 = NCT + 1                                                 PRONYX    356
        NRTP1 = NRT + 1                                                 PRONYX    357
        IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)                       PRONYX    358
        IF (N .LT. M) S(M) = 0.0D0                                      PRONYX    359
        IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)                         PRONYX    360
        E(M) = 0.0D0                                                    PRONYX    361
C                                                                       PRONYX    362
C if required, generate U.                                              PRONYX    363
C                                                                       PRONYX    364
        IF (.NOT.WANTU) GO TO 300                                       PRONYX    365
           IF (NCU .LT. NCTP1) GO TO 200                                PRONYX    366
           DO 190 J = NCTP1, NCU                                        PRONYX    367
              DO 180 I = 1, N                                           PRONYX    368
                 U(I,J) = 0.0D0                                         PRONYX    369
180           CONTINUE                                                  PRONYX    370
              U(J,J) = 1.0D0                                            PRONYX    371
190        CONTINUE                                                     PRONYX    372
200        CONTINUE                                                     PRONYX    373
           IF (NCT .LT. 1) GO TO 290                                    PRONYX    374
           DO 280 LL = 1, NCT                                           PRONYX    375
              L = NCT - LL + 1                                          PRONYX    376
              IF (S(L) .EQ. 0.0D0) GO TO 250                            PRONYX    377
                 LP1 = L + 1                                            PRONYX    378
                 IF (NCU .LT. LP1) GO TO 220                            PRONYX    379
                 DO 210 J = LP1, NCU                                    PRONYX    380
                    T = -DDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)           PRONYX    381
                    CALL DAXPY(N-L+1,T,U(L,L),1,U(L,J),1)               PRONYX    382
210              CONTINUE                                               PRONYX    383
220              CONTINUE                                               PRONYX    384
                 CALL DSCAL(N-L+1,-1.0D0,U(L,L),1)                      PRONYX    385
                 U(L,L) = 1.0D0 + U(L,L)                                PRONYX    386
                 LM1 = L - 1                                            PRONYX    387
                 IF (LM1 .LT. 1) GO TO 240                              PRONYX    388
                 DO 230 I = 1, LM1                                      PRONYX    389
                    U(I,L) = 0.0D0                                      PRONYX    390
230              CONTINUE                                               PRONYX    391
240              CONTINUE                                               PRONYX    392
              GO TO 270                                                 PRONYX    393
250           CONTINUE                                                  PRONYX    394
                 DO 260 I = 1, N                                        PRONYX    395
                    U(I,L) = 0.0D0                                      PRONYX    396
260              CONTINUE                                               PRONYX    397
                 U(L,L) = 1.0D0                                         PRONYX    398
270           CONTINUE                                                  PRONYX    399
280        CONTINUE                                                     PRONYX    400
290        CONTINUE                                                     PRONYX    401
300     CONTINUE                                                        PRONYX    402
C                                                                       PRONYX    403
C if it is required, generate V.                                        PRONYX    404
C                                                                       PRONYX    405
        IF (.NOT.WANTV) GO TO 350                                       PRONYX    406
           DO 340 LL = 1, P                                             PRONYX    407
              L = P - LL + 1                                            PRONYX    408
              LP1 = L + 1                                               PRONYX    409
              IF (L .GT. NRT) GO TO 320                                 PRONYX    410
              IF (E(L) .EQ. 0.0D0) GO TO 320                            PRONYX    411
                 DO 310 J = LP1, P                                      PRONYX    412
                    T = -DDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)       PRONYX    413
                    CALL DAXPY(P-L,T,V(LP1,L),1,V(LP1,J),1)             PRONYX    414
310              CONTINUE                                               PRONYX    415
320           CONTINUE                                                  PRONYX    416
              DO 330 I = 1, P                                           PRONYX    417
                 V(I,L) = 0.0D0                                         PRONYX    418
330           CONTINUE                                                  PRONYX    419
              V(L,L) = 1.0D0                                            PRONYX    420
340        CONTINUE                                                     PRONYX    421
350     CONTINUE                                                        PRONYX    422
C                                                                       PRONYX    423
C main iteration loop for the singular values.                          PRONYX    424
C                                                                       PRONYX    425
        MM = M                                                          PRONYX    426
        ITER = 0                                                        PRONYX    427
360     CONTINUE                                                        PRONYX    428
C                                                                       PRONYX    429
C quit if all the singular values have been found.                      PRONYX    430
C                                                                       PRONYX    431
C ...EXIT                                                               PRONYX    432
           IF (M .EQ. 0) GO TO 620                                      PRONYX    433
C                                                                       PRONYX    434
C if too many iterations have been performed, set flag and return.      PRONYX    435
C                                                                       PRONYX    436
           IF (ITER .LT. MAXIT) GO TO 370                               PRONYX    437
              INFO = M                                                  PRONYX    438
C ...EXIT                                                               PRONYX    439
              GO TO 620                                                 PRONYX    440
370        CONTINUE                                                     PRONYX    441
C                                                                       PRONYX    442
C This section of the program inspects for negligible elements in the S PRONYX    443
C arrays. On completion the variables KASE and L are set as follows.    PRONYX    444
C       KASE = 1  If S(M) and E(L-1) are negligible and L.LT.M          PRONYX    445
C       KASE = 2  If S(L) is negligible and L.LT.M                      PRONYX    446
C       KASE = 3  If E(L-1) is SMALL, L.LT.M, and S(L),., S(M) are not. PRONYX    447
C                 negligible (QR step).                                 PRONYX    448
C       KASE = 4  If E(M-1) is negligible (convergence).                PRONYX    449
C                                                                       PRONYX    450
           DO 390 LL = 1, M                                             PRONYX    451
              L = M - LL                                                PRONYX    452
C ...EXIT                                                               PRONYX    453
              IF (L .EQ. 0) GO TO 400                                   PRONYX    454
              TEST = DABS(S(L)) + DABS(S(L+1))                          PRONYX    455
              ZTEST = TEST + DABS(E(L))                                 PRONYX    456
              IF (ZTEST .NE. TEST) GO TO 380                            PRONYX    457
                 E(L) = 0.0D0                                           PRONYX    458
C ...EXIT                                                               PRONYX    459
                 GO TO 400                                              PRONYX    460
380           CONTINUE                                                  PRONYX    461
390        CONTINUE                                                     PRONYX    462
400        CONTINUE                                                     PRONYX    463
           IF (L .NE. M - 1) GO TO 410                                  PRONYX    464
              KASE = 4                                                  PRONYX    465
           GO TO 480                                                    PRONYX    466
410        CONTINUE                                                     PRONYX    467
              LP1 = L + 1                                               PRONYX    468
              MP1 = M + 1                                               PRONYX    469
              DO 430 LLS = LP1, MP1                                     PRONYX    470
                 LS = M - LLS + LP1                                     PRONYX    471
C ...EXIT                                                               PRONYX    472
                 IF (LS .EQ. L) GO TO 440                               PRONYX    473
                 TEST = 0.0D0                                           PRONYX    474
                 IF (LS .NE. M) TEST = TEST + DABS(E(LS))               PRONYX    475
                 IF (LS .NE. L + 1) TEST = TEST + DABS(E(LS-1))         PRONYX    476
                 ZTEST = TEST + DABS(S(LS))                             PRONYX    477
                 IF (ZTEST .NE. TEST) GO TO 420                         PRONYX    478
                    S(LS) = 0.0D0                                       PRONYX    479
C ..EXIT                                                                PRONYX    480
                    GO TO 440                                           PRONYX    481
420              CONTINUE                                               PRONYX    482
430           CONTINUE                                                  PRONYX    483
440           CONTINUE                                                  PRONYX    484
              IF (LS .NE. L) GO TO 450                                  PRONYX    485
                 KASE = 3                                               PRONYX    486
              GO TO 470                                                 PRONYX    487
450           CONTINUE                                                  PRONYX    488
              IF (LS .NE. M) GO TO 460                                  PRONYX    489
                 KASE = 1                                               PRONYX    490
              GO TO 470                                                 PRONYX    491
460           CONTINUE                                                  PRONYX    492
                 KASE = 2                                               PRONYX    493
                 L = LS                                                 PRONYX    494
470           CONTINUE                                                  PRONYX    495
480        CONTINUE                                                     PRONYX    496
           L = L + 1                                                    PRONYX    497
C                                                                       PRONYX    498
C perform the task indicated by KASE.                                   PRONYX    499
C                                                                       PRONYX    500
           GO TO (490,520,540,570), KASE                                PRONYX    501
C                                                                       PRONYX    502
C deflate negligible S(M).                                              PRONYX    503
C                                                                       PRONYX    504
490        CONTINUE                                                     PRONYX    505
              MM1 = M - 1                                               PRONYX    506
              F = E(M-1)                                                PRONYX    507
              E(M-1) = 0.0D0                                            PRONYX    508
              DO 510 KK = L, MM1                                        PRONYX    509
                 K = MM1 - KK + L                                       PRONYX    510
                 T1 = S(K)                                              PRONYX    511
                 CALL DROTG(T1,F,CS,SN)                                 PRONYX    512
                 S(K) = T1                                              PRONYX    513
                 IF (K .EQ. L) GO TO 500                                PRONYX    514
                    F = -SN*E(K-1)                                      PRONYX    515
                    E(K-1) = CS*E(K-1)                                  PRONYX    516
500              CONTINUE                                               PRONYX    517
                 IF (WANTV) CALL DROT(P,V(1,K),1,V(1,M),1,CS,SN)        PRONYX    518
510           CONTINUE                                                  PRONYX    519
           GO TO 610                                                    PRONYX    520
C                                                                       PRONYX    521
C split at negligible S(L).                                             PRONYX    522
C                                                                       PRONYX    523
520        CONTINUE                                                     PRONYX    524
              F = E(L-1)                                                PRONYX    525
              E(L-1) = 0.0D0                                            PRONYX    526
              DO 530 K = L, M                                           PRONYX    527
                 T1 = S(K)                                              PRONYX    528
                 CALL DROTG(T1,F,CS,SN)                                 PRONYX    529
                 S(K) = T1                                              PRONYX    530
                 F = -SN*E(K)                                           PRONYX    531
                 E(K) = CS*E(K)                                         PRONYX    532
                 IF (WANTU) CALL DROT(N,U(1,K),1,U(1,L-1),1,CS,SN)      PRONYX    533
530           CONTINUE                                                  PRONYX    534
           GO TO 610                                                    PRONYX    535
C                                                                       PRONYX    536
C perform one QR step.                                                  PRONYX    537
C                                                                       PRONYX    538
540        CONTINUE                                                     PRONYX    539
C                                                                       PRONYX    540
C calculate the shift.                                                  PRONYX    541
C                                                                       PRONYX    542
              SCALE = DMAX1(DABS(S(M)),DABS(S(M-1)),DABS(E(M-1)),       PRONYX    543
     *                    DABS(S(L)),DABS(E(L)))                        PRONYX    544
              SM = S(M)/SCALE                                           PRONYX    545
              SMM1 = S(M-1)/SCALE                                       PRONYX    546
              EMM1 = E(M-1)/SCALE                                       PRONYX    547
              SL = S(L)/SCALE                                           PRONYX    548
              EL = E(L)/SCALE                                           PRONYX    549
              B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0             PRONYX    550
              C = (SM*EMM1)**2                                          PRONYX    551
              SHIFT = 0.0D0                                             PRONYX    552
              IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 550            PRONYX    553
                 SHIFT = DSQRT(B**2+C)                                  PRONYX    554
                 IF (B .LT. 0.0D0) SHIFT = -SHIFT                       PRONYX    555
                 SHIFT = C/(B + SHIFT)                                  PRONYX    556
550           CONTINUE                                                  PRONYX    557
              F = (SL + SM)*(SL - SM) - SHIFT                           PRONYX    558
              G = SL*EL                                                 PRONYX    559
C                                                                       PRONYX    560
C chase zeros.                                                          PRONYX    561
C                                                                       PRONYX    562
              MM1 = M - 1                                               PRONYX    563
              DO 560 K = L, MM1                                         PRONYX    564
                 CALL DROTG(F,G,CS,SN)                                  PRONYX    565
                 IF (K .NE. L) E(K-1) = F                               PRONYX    566
                 F = CS*S(K) + SN*E(K)                                  PRONYX    567
                 E(K) = CS*E(K) - SN*S(K)                               PRONYX    568
                 G = SN*S(K+1)                                          PRONYX    569
                 S(K+1) = CS*S(K+1)                                     PRONYX    570
                 IF (WANTV) CALL DROT(P,V(1,K),1,V(1,K+1),1,CS,SN)      PRONYX    571
                 CALL DROTG(F,G,CS,SN)                                  PRONYX    572
                 S(K) = F                                               PRONYX    573
                 F = CS*E(K) + SN*S(K+1)                                PRONYX    574
                 S(K+1) = -SN*E(K) + CS*S(K+1)                          PRONYX    575
                 G = SN*E(K+1)                                          PRONYX    576
                 E(K+1) = CS*E(K+1)                                     PRONYX    577
                 IF (WANTU .AND. K .LT. N)                              PRONYX    578
     *            CALL DROT(N,U(1,K),1,U(1,K+1),1,CS,SN)                PRONYX    579
560           CONTINUE                                                  PRONYX    580
              E(M-1) = F                                                PRONYX    581
              ITER = ITER + 1                                           PRONYX    582
           GO TO 610                                                    PRONYX    583
C                                                                       PRONYX    584
C convergence.                                                          PRONYX    585
C                                                                       PRONYX    586
570        CONTINUE                                                     PRONYX    587
C                                                                       PRONYX    588
C make the singular value positive.                                     PRONYX    589
C                                                                       PRONYX    590
              IF (S(L) .GE. 0.0D0) GO TO 580                            PRONYX    591
                 S(L) = -S(L)                                           PRONYX    592
                 IF (WANTV) CALL DSCAL(P,-1.0D0,V(1,L),1)               PRONYX    593
580           CONTINUE                                                  PRONYX    594
C                                                                       PRONYX    595
C order the singular value.                                             PRONYX    596
C                                                                       PRONYX    597
590           IF (L .EQ. MM) GO TO 600                                  PRONYX    598
C ...EXIT                                                               PRONYX    599
                 IF (S(L) .GE. S(L+1)) GO TO 600                        PRONYX    600
                 T = S(L)                                               PRONYX    601
                 S(L) = S(L+1)                                          PRONYX    602
                 S(L+1) = T                                             PRONYX    603
                 IF (WANTV .AND. L .LT. P)                              PRONYX    604
     *            CALL DSWAP(P,V(1,L),1,V(1,L+1),1)                     PRONYX    605
                 IF (WANTU .AND. L .LT. N)                              PRONYX    606
     *            CALL DSWAP(N,U(1,L),1,U(1,L+1),1)                     PRONYX    607
                 L = L + 1                                              PRONYX    608
              GO TO 590                                                 PRONYX    609
600           CONTINUE                                                  PRONYX    610
              ITER = 0                                                  PRONYX    611
              M = M - 1                                                 PRONYX    612
610        CONTINUE                                                     PRONYX    613
        GO TO 360                                                       PRONYX    614
620     CONTINUE                                                        PRONYX    615
        RETURN                                                          PRONYX    616
        END                                                             PRONYX    617
C                                                                       PRONYX    618
        SUBROUTINE DQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB)                 PRONYX    619
CD           This routine uses householder transformations to compute   PRONYX    620
CD           the QR factorization of N by P matx X. Col pivoting based  PRONYX    621
CD           on the 2-NORMS of the reduced columns may be performed at  PRONYX    622
CD           the users option.                                          PRONYX    623
CD                                                                      PRONYX    624
CD     PARAMETERS :                                                     PRONYX    625
CD                                                                      PRONYX    626
CD      ON ENTRY :                                                      PRONYX    627
CD                                                                      PRONYX    628
CD       X   - DOUBLE PREC(LDX,P), where LDX .GE. N. The matrix whosed  PRONYX    629
CD             decomposition is to be computed.                         PRONYX    630
CD       LDX - INTEGER, The leading dimension of the array X.           PRONYX    631
CD       N   - INTEGER, The number of rows of the matrix X.             PRONYX    632
CD       P   - INTEGER, The number of columns of the matrix X.          PRONYX    633
CD       JPVT-  INT(P),JPVT contains integers controling  selection of  PRONYX    634
CD             the pivot cols.The K-TH col X(K) of X is placed in one   PRONYX    635
CD             of three classes according to the value of JPVT(K).      PRONYX    636
CD             If JPVT(K) .GT. 0, then X(K) is an initial column.       PRONYX    637
CD             If JPVT(K) .EQ. 0, then X(K) is a free column.           PRONYX    638
CD             If JPVT(K) .LT. 0, then X(K) is a final column.          PRONYX    639
CD                                                                      PRONYX    640
CD             Before the decomp is computed, initial cols are moved to PRONYX    641
CD             the beginning of the array X and final cols to the end.  PRONYX    642
CD             Both init and final cols are frozen in place during the  PRONYX    643
CD             comp and only free cols are moved. At the K-TH staged    PRONYX    644
CD             of reduction,If X(k) is filled by a free column it is    PRONYX    645
CD             interch with the free col of largest reduced norm. JPVT  PRONYX    646
CD             is not referenced if JOB .EQ. 0.                         PRONYX    647
CD       WORK- DOUBLE PREC(P), Work array. WORK is not referenced if J  PRONYX    648
CD             JOB .EQ. 0.                                              PRONYX    649
CD       JOB - INTEGER, Initiates column pivoting.                      PRONYX    650
CD             If JOB .EQ. 0, no pivoting is done.                      PRONYX    651
CD             If JOB .NE. 0, pivoting is done.                         PRONYX    652
CD ON RETURN :                                                          PRONYX    653
CD                                                                      PRONYX    654
CD       X   - X contains in its upper triangle  matrix RT              PRONYX    655
CD             of the QR factorization. Below its diagonal X contains   PRONYX    656
CD             information from which the orthogonal part of the decomp PRONYX    657
CD             can be recovered.Note if pivoting has been requested,the PRONYX    658
CD             decomp is not that of the original matrix  but that of X PRONYX    659
CD             with its columns permuted as described by JPVT.          PRONYX    660
CD     QRAUX - DOUBLE PRECISION(P). QRAUX contains further information  PRONYX    661
CD             required to recover the orthogonal part of the decomp.   PRONYX    662
CD      JPVT - JPVT(K) contains the index of the col of the orig. matrX PRONYX    663
CD             that has been interchged into the K-TH col, if pivoting  PRONYX    664
CD             was requested.                                           PRONYX    665
CD                                                                      PRONYX    666
CD LINPACK. THIS VERSION DATED 08/14/78 .                               PRONYX    667
CD G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.          PRONYX    668
CD                                                                      PRONYX    669
CD SUBROUTINES AND FUNCTIONS :   BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2      PRONYX    670
CD                               FORTRAN DABS,DMAX1,MIN0,DSQRT          PRONYX    671
C Argument variables                                                    PRONYX    672
        INTEGER LDX,N,P,JOB                                             PRONYX    673
        INTEGER JPVT(1)                                                 PRONYX    674
        DOUBLE PRECISION X(LDX,1),QRAUX(1),WORK(1)                      PRONYX    675
C Internal variables                                                    PRONYX    676
        INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU                               PRONYX    677
        DOUBLE PRECISION MAXNRM,DNRM2,TT                                PRONYX    678
        DOUBLE PRECISION DDOT,NRMXL,T                                   PRONYX    679
        LOGICAL NEGJ,SWAPJ                                              PRONYX    680
C                                                                       PRONYX    681
        PL = 1                                                          PRONYX    682
        PU = 0                                                          PRONYX    683
        IF (JOB .EQ. 0) GO TO 60                                        PRONYX    684
C                                                                       PRONYX    685
C Pivoting has been requested.  rearrange the columns according to JPVT.PRONYX    686
C                                                                       PRONYX    687
           DO 20 J = 1, P                                               PRONYX    688
              SWAPJ = JPVT(J) .GT. 0                                    PRONYX    689
              NEGJ = JPVT(J) .LT. 0                                     PRONYX    690
              JPVT(J) = J                                               PRONYX    691
              IF (NEGJ) JPVT(J) = -J                                    PRONYX    692
              IF (.NOT.SWAPJ) GO TO 10                                  PRONYX    693
                 IF (J .NE. PL) CALL DSWAP(N,X(1,PL),1,X(1,J),1)        PRONYX    694
                 JPVT(J) = JPVT(PL)                                     PRONYX    695
                 JPVT(PL) = J                                           PRONYX    696
                 PL = PL + 1                                            PRONYX    697
10            CONTINUE                                                  PRONYX    698
20         CONTINUE                                                     PRONYX    699
           PU = P                                                       PRONYX    700
           DO 50 JJ = 1, P                                              PRONYX    701
              J = P - JJ + 1                                            PRONYX    702
              IF (JPVT(J) .GE. 0) GO TO 40                              PRONYX    703
                 JPVT(J) = -JPVT(J)                                     PRONYX    704
                 IF (J .EQ. PU) GO TO 30                                PRONYX    705
                    CALL DSWAP(N,X(1,PU),1,X(1,J),1)                    PRONYX    706
                    JP = JPVT(PU)                                       PRONYX    707
                    JPVT(PU) = JPVT(J)                                  PRONYX    708
                    JPVT(J) = JP                                        PRONYX    709
30               CONTINUE                                               PRONYX    710
                 PU = PU - 1                                            PRONYX    711
40            CONTINUE                                                  PRONYX    712
50         CONTINUE                                                     PRONYX    713
60      CONTINUE                                                        PRONYX    714
C                                                                       PRONYX    715
C compute the norms of the free columns.                                PRONYX    716
C                                                                       PRONYX    717
        IF (PU .LT. PL) GO TO 80                                        PRONYX    718
        DO 70 J = PL, PU                                                PRONYX    719
           QRAUX(J) = DNRM2(N,X(1,J),1)                                 PRONYX    720
           WORK(J) = QRAUX(J)                                           PRONYX    721
70      CONTINUE                                                        PRONYX    722
80      CONTINUE                                                        PRONYX    723
C                                                                       PRONYX    724
C perform the HOUSEHOLDER reduction of X.                               PRONYX    725
C                                                                       PRONYX    726
        LUP = MIN0(N,P)                                                 PRONYX    727
        DO 200 L = 1, LUP                                               PRONYX    728
           IF (L .LT. PL .OR. L .GE. PU) GO TO 120                      PRONYX    729
C                                                                       PRONYX    730
C locate the column of largest norm and bring it into the pivot positionPRONYX    731
C                                                                       PRONYX    732
              MAXNRM = 0.0D0                                            PRONYX    733
              MAXJ = L                                                  PRONYX    734
              DO 100 J = L, PU                                          PRONYX    735
                 IF (QRAUX(J) .LE. MAXNRM) GO TO 90                     PRONYX    736
                    MAXNRM = QRAUX(J)                                   PRONYX    737
                    MAXJ = J                                            PRONYX    738
90               CONTINUE                                               PRONYX    739
100           CONTINUE                                                  PRONYX    740
              IF (MAXJ .EQ. L) GO TO 110                                PRONYX    741
                 CALL DSWAP(N,X(1,L),1,X(1,MAXJ),1)                     PRONYX    742
                 QRAUX(MAXJ) = QRAUX(L)                                 PRONYX    743
                 WORK(MAXJ) = WORK(L)                                   PRONYX    744
                 JP = JPVT(MAXJ)                                        PRONYX    745
                 JPVT(MAXJ) = JPVT(L)                                   PRONYX    746
                 JPVT(L) = JP                                           PRONYX    747
110          CONTINUE                                                   PRONYX    748
120        CONTINUE                                                     PRONYX    749
           QRAUX(L) = 0.0D0                                             PRONYX    750
           IF (L .EQ. N) GO TO 190                                      PRONYX    751
C                                                                       PRONYX    752
C compute the HOUSEHOLDER transformation for column L.                  PRONYX    753
C                                                                       PRONYX    754
              NRMXL = DNRM2(N-L+1,X(L,L),1)                             PRONYX    755
              IF (NRMXL .EQ. 0.0D0) GO TO 180                           PRONYX    756
                 IF (X(L,L) .NE. 0.0D0) NRMXL = DSIGN(NRMXL,X(L,L))     PRONYX    757
                 CALL DSCAL(N-L+1,1.0D0/NRMXL,X(L,L),1)                 PRONYX    758
                 X(L,L) = 1.0D0 + X(L,L)                                PRONYX    759
C                                                                       PRONYX    760
C apply the transformation to the remaining columns, updating the norms.PRONYX    761
C                                                                       PRONYX    762
                 LP1 = L + 1                                            PRONYX    763
                 IF (P .LT. LP1) GO TO 170                              PRONYX    764
                 DO 160 J = LP1, P                                      PRONYX    765
                    T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)           PRONYX    766
                    CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1)               PRONYX    767
                    IF (J .LT. PL .OR. J .GT. PU) GO TO 150             PRONYX    768
                    IF (QRAUX(J) .EQ. 0.0D0) GO TO 150                  PRONYX    769
                       TT = 1.0D0 - (DABS(X(L,J))/QRAUX(J))**2          PRONYX    770
                       TT = DMAX1(TT,0.0D0)                             PRONYX    771
                       T = TT                                           PRONYX    772
                       TT = 1.0D0 + 0.05D0*TT*(QRAUX(J)/WORK(J))**2     PRONYX    773
                       IF (TT .EQ. 1.0D0) GO TO 130                     PRONYX    774
                          QRAUX(J) = QRAUX(J)*DSQRT(T)                  PRONYX    775
                       GO TO 140                                        PRONYX    776
130                    CONTINUE                                         PRONYX    777
                          QRAUX(J) = DNRM2(N-L,X(L+1,J),1)              PRONYX    778
                          WORK(J) = QRAUX(J)                            PRONYX    779
140                    CONTINUE                                         PRONYX    780
150                 CONTINUE                                            PRONYX    781
160              CONTINUE                                               PRONYX    782
170              CONTINUE                                               PRONYX    783
C                                                                       PRONYX    784
C save the transformation.                                              PRONYX    785
C                                                                       PRONYX    786
                 QRAUX(L) = X(L,L)                                      PRONYX    787
                 X(L,L) = -NRMXL                                        PRONYX    788
180           CONTINUE                                                  PRONYX    789
190        CONTINUE                                                     PRONYX    790
200     CONTINUE                                                        PRONYX    791
        RETURN                                                          PRONYX    792
        END                                                             PRONYX    793
        SUBROUTINE DQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO)    PRONYX    794
CD                                                                      PRONYX    795
CD           This routine applies the output of DQRDC to compute        PRONYX    796
CD           coordinate transformation projections, and least sqs       PRONYX    797
CD           solutions. For K .LE. MIN(N,P),let XK be the matrix        PRONYX    798
CD               XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K)))          PRONYX    799
CD           formed from cols JPVT(1),.,JPVT(K) of the original N X P   PRONYX    800
CD           matrix X that was input to DQRDC(if no pivoting was done,  PRONYX    801
CD           XK consists of first K columns of X in their original      PRONYX    802
CD           order). DQRDC produces a factored orthogonal matrix Q and  PRONYX    803
CD           an upper triang matrix R such that                         PRONYX    804
CD               XK = Q * (R)                                           PRONYX    805
CD                        (0)                                           PRONYX    806
CD     This info is contained in coded form in the arrays X and QRAUX.  PRONYX    807
CD                                                                      PRONYX    808
CD  PARAMETERS :                                                        PRONYX    809
CD                                                                      PRONYX    810
CD   ON ENTRY :                                                         PRONYX    811
CD       X     DOUBLE PRECISION(LDX,P), X contains the output of DQRDC. PRONYX    812
CD       LDX   INTEGER, The leading dimension of the array X.           PRONYX    813
CD       N     INT,The no. of rows of the matrix XK. It must have the   PRONYX    814
CD             same value as N in DQRDC.                                PRONYX    815
CD       K     INT, K is the no. of cols of the matrix XK. K must not   PRONYX    816
CD             be > MIN(N,P), where P is the same as in the calling     PRONYX    817
CD             sequence to DQRDC.                                       PRONYX    818
CD       QRAUX DOUBLE PRECISION(P), The auxiliary output from DQRDC.    PRONYX    819
CD       Y     DOUBLE PRECISION(N), contains an N-VECTOR that is to be  PRONYX    820
CD             manipulated by DQRSL.                                    PRONYX    821
CD       JOB   INT, specifies what is to be computed.JOB has the deci   PRONYX    822
CD             expansion ABCDE, with the following meaning.             PRONYX    823
CD                       If A.NE.0, compute QY.                         PRONYX    824
CD                       If B,C,D, or E .NE. 0, compute QTY.            PRONYX    825
CD                       If C.NE.0, compute B.                          PRONYX    826
CD                       If D.NE.0, compute RSD.                        PRONYX    827
CD                       If E.NE.0, compute XB.                         PRONYX    828
CD             Note that a request to comp B, RSD, or XB automat        PRONYX    829
CD             triggers the comp of QTY, for which an array must be     PRONYX    830
CD             provided in the calling sequence.                        PRONYX    831
CD   ON RETURN :                                                        PRONYX    832
CD       QY    DOUBLE PREC(N), contains Q*Y, if its comp has been       PRONYX    833
CD             requested.                                               PRONYX    834
CD       QTY   DOUBLE PREC(N),contains TRANS(Q)*Y, if its comp has      PRONYX    835
CD             been requested. trans(q) IS THE transpose of the matx q. PRONYX    836
CD       B     DOUBLE PREC(K),contains the soln of the least squaress   PRONYX    837
CD             problem                                                  PRONYX    838
CD             MINIMIZE NORM2(Y - XK*B),                                PRONYX    839
CD             if its comp has been requested. (note that if pivoting   PRONYX    840
CD             was requested in DQRDC, the J-TH component of B will be  PRONYX    841
CD             associated with col JPVT(J) of the orig matrix X that    PRONYX    842
CD             was input into DQRDC.)                                   PRONYX    843
CD       RSD   DOUBLE PRE(N), The least sqs residual Y - XK*B, if its   PRONYX    844
CD             comp has been requested. RSD is also the orthogonal      PRONYX    845
CD             projection of Y onto the orthogonal complement of the co PRONYX    846
CD             space of XK.                                             PRONYX    847
CD       XB    DOUBLE PREC(N), The least sqs approximation XK*B, ifits  PRONYX    848
CD             its comp has been requested. XB is also the orthogonal   PRONYX    849
CD             projection of Y onto the column space of X.              PRONYX    850
CD       INFO  INT, INFO is zero unless the computation of B has been   PRONYX    851
CD             requested and R is  singular. In this case, INFO is the  PRONYX    852
CD             index of the first zero diag element of R and B is left  PRONYX    853
CD             unaltered.                                               PRONYX    854
CD                                                                      PRONYX    855
CD The parameters QY, QTY, B, RSD, and XB are not referenced if their   PRONYX    856
CD computation is not requested and in this case can be replaced by     PRONYX    857
CD dummy variables in the calling program. To save storage, the user    PRONYX    858
CD may in some use the same array for different parameters in the       PRONYX    859
CD calling sequence. A frequently occuring example is when one wishes   PRONYX    860
CD to compute any of B, RS and does not need Y or QTY. In this case one PRONYX    861
CD may identify Y, QTY, and  B, RSD, or XB, while providing separate    PRONYX    862
CD arrays for anything else that computed. Thus the calling sequence    PRONYX    863
CD       CALL DQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO)           PRONYX    864
CD will result in the computation of B and RSD, with RSD overwriting Y. PRONYX    865
CD generally, each item in the following list contains groups of permissPRONYX    866
CD identifications for a single callinng sequence.                      PRONYX    867
CD                                                                      PRONYX    868
CD       1. (Y,QTY,B) (RSD) (XB) (QY)                                   PRONYX    869
CD       2. (Y,QTY,RSD) (B) (XB) (QY)                                   PRONYX    870
CD       3. (Y,QTY,XB) (B) (RSD) (QY)                                   PRONYX    871
CD       4. (Y,QY) (QTY,B) (RSD) (XB)                                   PRONYX    872
CD       5. (Y,QY) (QTY,RSD) (B) (XB)                                   PRONYX    873
CD       6. (Y,QY) (QTY,XB) (B) (RSD)                                   PRONYX    874
CD                                                                      PRONYX    875
CD In any group the value returned in the array allocated to the group  PRONYX    876
CD corresponds to the last member of the group.                         PRONYX    877
CD                                                                      PRONYX    878
CD LINPACK. THIS VERSION DATED 08/14/78 .                               PRONYX    879
CD G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.          PRONYX    880
CD                                                                      PRONYX    881
CD SUBROUTINES AND FUNCTIONS :   BLAS DAXPY,DCOPY,DDOT                  PRONYX    882
CD                               FORTRAN DABS,MIN0,MOD                  PRONYX    883
CD                                                                      PRONYX    884
C Argument variables                                                    PRONYX    885
        INTEGER LDX,N,K,JOB,INFO                                        PRONYX    886
        DOUBLE PRECISION X(LDX,1),QRAUX(1),Y(1),QY(1),QTY(1),B(1),      PRONYX    887
     *                RSD(1),XB(1)                                      PRONYX    888
C Internal variables                                                    PRONYX    889
        INTEGER I,J,JJ,JU,KP1                                           PRONYX    890
        DOUBLE PRECISION DDOT,T,TEMP                                    PRONYX    891
        LOGICAL CB,CQY,CQTY,CR,CXB                                      PRONYX    892
C                                                                       PRONYX    893
C set INFO flag.                                                        PRONYX    894
C                                                                       PRONYX    895
        INFO = 0                                                        PRONYX    896
C                                                                       PRONYX    897
C determine what is to be computed.                                     PRONYX    898
C                                                                       PRONYX    899
        CQY = JOB/10000 .NE. 0                                          PRONYX    900
        CQTY = MOD(JOB,10000) .NE. 0                                    PRONYX    901
        CB = MOD(JOB,1000)/100 .NE. 0                                   PRONYX    902
        CR = MOD(JOB,100)/10 .NE. 0                                     PRONYX    903
        CXB = MOD(JOB,10) .NE. 0                                        PRONYX    904
        JU = MIN0(K,N-1)                                                PRONYX    905
C                                                                       PRONYX    906
C special action when N=1.                                              PRONYX    907
C                                                                       PRONYX    908
        IF (JU .NE. 0) GO TO 40                                         PRONYX    909
           IF (CQY) QY(1) = Y(1)                                        PRONYX    910
           IF (CQTY) QTY(1) = Y(1)                                      PRONYX    911
           IF (CXB) XB(1) = Y(1)                                        PRONYX    912
           IF (.NOT.CB) GO TO 30                                        PRONYX    913
              IF (X(1,1) .NE. 0.0D0) GO TO 10                           PRONYX    914
                 INFO = 1                                               PRONYX    915
              GO TO 20                                                  PRONYX    916
10            CONTINUE                                                  PRONYX    917
                 B(1) = Y(1)/X(1,1)                                     PRONYX    918
20            CONTINUE                                                  PRONYX    919
30         CONTINUE                                                     PRONYX    920
           IF (CR) RSD(1) = 0.0D0                                       PRONYX    921
        GO TO 250                                                       PRONYX    922
40      CONTINUE                                                        PRONYX    923
C                                                                       PRONYX    924
C set up to compute QY or QTY.                                          PRONYX    925
C                                                                       PRONYX    926
           IF (CQY) CALL DCOPY(N,Y,1,QY,1)                              PRONYX    927
           IF (CQTY) CALL DCOPY(N,Y,1,QTY,1)                            PRONYX    928
           IF (.NOT.CQY) GO TO 70                                       PRONYX    929
C                                                                       PRONYX    930
C compute QY.                                                           PRONYX    931
C                                                                       PRONYX    932
              DO 60 JJ = 1, JU                                          PRONYX    933
                 J = JU - JJ + 1                                        PRONYX    934
                 IF (QRAUX(J) .EQ. 0.0D0) GO TO 50                      PRONYX    935
                    TEMP = X(J,J)                                       PRONYX    936
                    X(J,J) = QRAUX(J)                                   PRONYX    937
                    T = -DDOT(N-J+1,X(J,J),1,QY(J),1)/X(J,J)            PRONYX    938
                    CALL DAXPY(N-J+1,T,X(J,J),1,QY(J),1)                PRONYX    939
                    X(J,J) = TEMP                                       PRONYX    940
50               CONTINUE                                               PRONYX    941
60            CONTINUE                                                  PRONYX    942
70         CONTINUE                                                     PRONYX    943
           IF (.NOT.CQTY) GO TO 100                                     PRONYX    944
C                                                                       PRONYX    945
C compute TRANS(Q)*Y.                                                   PRONYX    946
C                                                                       PRONYX    947
              DO 90 J = 1, JU                                           PRONYX    948
                 IF (QRAUX(J) .EQ. 0.0D0) GO TO 80                      PRONYX    949
                    TEMP = X(J,J)                                       PRONYX    950
                    X(J,J) = QRAUX(J)                                   PRONYX    951
                    T = -DDOT(N-J+1,X(J,J),1,QTY(J),1)/X(J,J)           PRONYX    952
                    CALL DAXPY(N-J+1,T,X(J,J),1,QTY(J),1)               PRONYX    953
                    X(J,J) = TEMP                                       PRONYX    954
80               CONTINUE                                               PRONYX    955
90            CONTINUE                                                  PRONYX    956
100        CONTINUE                                                     PRONYX    957
C                                                                       PRONYX    958
C set up to compute B, RSD, or XB.                                      PRONYX    959
C                                                                       PRONYX    960
           IF (CB) CALL DCOPY(K,QTY,1,B,1)                              PRONYX    961
           KP1 = K + 1                                                  PRONYX    962
           IF (CXB) CALL DCOPY(K,QTY,1,XB,1)                            PRONYX    963
           IF (CR .AND. K .LT. N) CALL DCOPY(N-K,QTY(KP1),1,RSD(KP1),1) PRONYX    964
           IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 120                      PRONYX    965
              DO 110 I = KP1, N                                         PRONYX    966
                 XB(I) = 0.0D0                                          PRONYX    967
110           CONTINUE                                                  PRONYX    968
120        CONTINUE                                                     PRONYX    969
           IF (.NOT.CR) GO TO 140                                       PRONYX    970
              DO 130 I = 1, K                                           PRONYX    971
                 RSD(I) = 0.0D0                                         PRONYX    972
130           CONTINUE                                                  PRONYX    973
140        CONTINUE                                                     PRONYX    974
           IF (.NOT.CB) GO TO 190                                       PRONYX    975
C                                                                       PRONYX    976
C compute B.                                                            PRONYX    977
C                                                                       PRONYX    978
              DO 170 JJ = 1, K                                          PRONYX    979
                 J = K - JJ + 1                                         PRONYX    980
                 IF (X(J,J) .NE. 0.0D0) GO TO 150                       PRONYX    981
                    INFO = J                                            PRONYX    982
C ...EXIT                                                               PRONYX    983
                    GO TO 180                                           PRONYX    984
150              CONTINUE                                               PRONYX    985
                 B(J) = B(J)/X(J,J)                                     PRONYX    986
                 IF (J .EQ. 1) GO TO 160                                PRONYX    987
                    T = -B(J)                                           PRONYX    988
                    CALL DAXPY(J-1,T,X(1,J),1,B,1)                      PRONYX    989
160              CONTINUE                                               PRONYX    990
170           CONTINUE                                                  PRONYX    991
180           CONTINUE                                                  PRONYX    992
190        CONTINUE                                                     PRONYX    993
           IF (.NOT.CR .AND. .NOT.CXB) GO TO 240                        PRONYX    994
C                                                                       PRONYX    995
C compute RSD or XB as required.                                        PRONYX    996
C                                                                       PRONYX    997
              DO 230 JJ = 1, JU                                         PRONYX    998
                 J = JU - JJ + 1                                        PRONYX    999
                 IF (QRAUX(J) .EQ. 0.0D0) GO TO 220                     PRONYX   1000
                    TEMP = X(J,J)                                       PRONYX   1001
                    X(J,J) = QRAUX(J)                                   PRONYX   1002
                    IF (.NOT.CR) GO TO 200                              PRONYX   1003
                       T = -DDOT(N-J+1,X(J,J),1,RSD(J),1)/X(J,J)        PRONYX   1004
                       CALL DAXPY(N-J+1,T,X(J,J),1,RSD(J),1)            PRONYX   1005
200                 CONTINUE                                            PRONYX   1006
                    IF (.NOT.CXB) GO TO 210                             PRONYX   1007
                       T = -DDOT(N-J+1,X(J,J),1,XB(J),1)/X(J,J)         PRONYX   1008
                       CALL DAXPY(N-J+1,T,X(J,J),1,XB(J),1)             PRONYX   1009
210                 CONTINUE                                            PRONYX   1010
                    X(J,J) = TEMP                                       PRONYX   1011
220              CONTINUE                                               PRONYX   1012
230           CONTINUE                                                  PRONYX   1013
240        CONTINUE                                                     PRONYX   1014
250     CONTINUE                                                        PRONYX   1015
        RETURN                                                          PRONYX   1016
        END                                                             PRONYX   1017
C***********************************************************************PRONYX   1018
C                                                                      *PRONYX   1019
C       Functions and subroutines used by double precision LINPACK     *PRONYX   1020
C                                                                      *PRONYX   1021
C       DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX)                     *PRONYX   1022
C       SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)                         *PRONYX   1023
C       SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)                            *PRONYX   1024
C       DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)              *PRONYX   1025
C       DOUBLE PRECISION FUNCTION DNRM2 (N,DX,INCX)                    *PRONYX   1026
C       SUBROUTINE DROT (N,DX,INCX,DY,INCY,C,S)                        *PRONYX   1027
C       SUBROUTINE DROTG(DA,DB,C,S)                                    *PRONYX   1028
C       SUBROUTINE DSCAL(N,DA,DX,INCX)                                 *PRONYX   1029
C       SUBROUTINE DSWAP (N,DX,INCX,DY,INCY)                           *PRONYX   1030
C       INTEGER FUNCTION IDAMAX(N,DX,INCX)                             *PRONYX   1031
C                                                                      *PRONYX   1032
C***********************************************************************PRONYX   1033
                                                                        PRONYX   1034
C                                                                       PRONYX   1035
        DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX)                      PRONYX   1036
CD                                                                      PRONYX   1037
CD      This routine takes the sum of the absolute values of an array.  PRONYX   1038
CD                                                                      PRONYX   1039
CD      CALLING PARAMETERS :                                            PRONYX   1040
CD                                                                      PRONYX   1041
CD                     N    - NO. OF PONIT                              PRONYX   1042
CD                     DX   - INCOMING ARRAY                            PRONYX   1043
CD                     INCX - INCREMENT                                 PRONYX   1044
C JACK DONGARRA, LINPACK, 3/11/78.                                      PRONYX   1045
        DOUBLE PRECISION DX(1),DTEMP                                    PRONYX   1046
        INTEGER I,INCX,M,MP1,N,NINCX                                    PRONYX   1047
        DASUM = 0.0D0                                                   PRONYX   1048
        DTEMP = 0.0D0                                                   PRONYX   1049
        IF(N.LE.0)RETURN                                                PRONYX   1050
        IF(INCX.EQ.1)GO TO 20                                           PRONYX   1051
C                                                                       PRONYX   1052
C code for increment not equal to 1                                     PRONYX   1053
C                                                                       PRONYX   1054
        NINCX = N*INCX                                                  PRONYX   1055
        DO 10 I = 1,NINCX,INCX                                          PRONYX   1056
            DTEMP = DTEMP + DABS(DX(I))                                 PRONYX   1057
10      CONTINUE                                                        PRONYX   1058
        DASUM = DTEMP                                                   PRONYX   1059
        RETURN                                                          PRONYX   1060
C                                                                       PRONYX   1061
C code for increment equal to 1                                         PRONYX   1062
C                                                                       PRONYX   1063
C clean-up loop                                                         PRONYX   1064
C                                                                       PRONYX   1065
20      M = MOD(N,6)                                                    PRONYX   1066
        IF( M .EQ. 0 ) GO TO 40                                         PRONYX   1067
        DO 30 I = 1,M                                                   PRONYX   1068
          DTEMP = DTEMP + DABS(DX(I))                                   PRONYX   1069
30      CONTINUE                                                        PRONYX   1070
        IF( N .LT. 6 ) GO TO 60                                         PRONYX   1071
40      MP1 = M + 1                                                     PRONYX   1072
        DO 50 I = MP1,N,6                                               PRONYX   1073
          DTEMP = DTEMP+DABS(DX(I))+DABS(DX(I + 1))+DABS(DX(I + 2))     PRONYX   1074
     *  + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5))           PRONYX   1075
50      CONTINUE                                                        PRONYX   1076
60      DASUM = DTEMP                                                   PRONYX   1077
        RETURN                                                          PRONYX   1078
        END                                                             PRONYX   1079
                                                                        PRONYX   1080
C                                                                       PRONYX   1081
        SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)                          PRONYX   1082
CD                                                                      PRONYX   1083
CD            This routine performs a function of constant times        PRONYX   1084
CD            a vector plus a vector. It uses unrolled loops for        PRONYX   1085
CD            increments equal to one. I.E.                             PRONYX   1086
CD                                                                      PRONYX   1087
CD                DY(.) = DY(.) + DA*DX(.)                              PRONYX   1088
CD                                                                      PRONYX   1089
CD            CALLING PARAMETERS :                                      PRONYX   1090
CD                                                                      PRONYX   1091
CD                    N   - NO. OF POINTS IN THE ARRAY                  PRONYX   1092
CD                    DA  - INCOMING CONSTANT                           PRONYX   1093
CD                    DX  - INCOMING X-ARRAY                            PRONYX   1094
CD                  INCX  - INCREMENT OF X                              PRONYX   1095
CD                    DY  - INCOMING Y-ARRAY                            PRONYX   1096
CD                  INCY  - INCREMENT OF Y                              PRONYX   1097
CD                                                                      PRONYX   1098
C JACK DONGARRA, LINPACK, 3/11/78.                                      PRONYX   1099
        DOUBLE PRECISION DX(1),DY(1),DA                                 PRONYX   1100
        INTEGER I,INCX,INCY,M,MP1,N                                     PRONYX   1101
        IF(N.LE.0)RETURN                                                PRONYX   1102
        IF (DA .EQ. 0.0D0) RETURN                                       PRONYX   1103
        IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20                             PRONYX   1104
C                                                                       PRONYX   1105
C code for unequal increments or equal increments not equal to 1        PRONYX   1106
C                                                                       PRONYX   1107
        IX = 1                                                          PRONYX   1108
        IY = 1                                                          PRONYX   1109
        IF(INCX.LT.0)IX = (-N+1)*INCX + 1                               PRONYX   1110
        IF(INCY.LT.0)IY = (-N+1)*INCY + 1                               PRONYX   1111
        DO 10 I = 1,N                                                   PRONYX   1112
          DY(IY) = DY(IY) + DA*DX(IX)                                   PRONYX   1113
          IX = IX + INCX                                                PRONYX   1114
          IY = IY + INCY                                                PRONYX   1115
10      CONTINUE                                                        PRONYX   1116
        RETURN                                                          PRONYX   1117
C                                                                       PRONYX   1118
C code for both increments equal to 1                                   PRONYX   1119
C                                                                       PRONYX   1120
C clean-up loop                                                         PRONYX   1121
C                                                                       PRONYX   1122
20      M = MOD(N,4)                                                    PRONYX   1123
        IF( M .EQ. 0 ) GO TO 40                                         PRONYX   1124
        DO 30 I = 1,M                                                   PRONYX   1125
          DY(I) = DY(I) + DA*DX(I)                                      PRONYX   1126
30      CONTINUE                                                        PRONYX   1127
        IF( N .LT. 4 ) RETURN                                           PRONYX   1128
40      MP1 = M + 1                                                     PRONYX   1129
        DO 50 I = MP1,N,4                                               PRONYX   1130
          DY(I) = DY(I) + DA*DX(I)                                      PRONYX   1131
          DY(I + 1) = DY(I + 1) + DA*DX(I + 1)                          PRONYX   1132
          DY(I + 2) = DY(I + 2) + DA*DX(I + 2)                          PRONYX   1133
          DY(I + 3) = DY(I + 3) + DA*DX(I + 3)                          PRONYX   1134
50      CONTINUE                                                        PRONYX   1135
        RETURN                                                          PRONYX   1136
        END                                                             PRONYX   1137
C                                                                       PRONYX   1138
        SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)                             PRONYX   1139
CD                                                                      PRONYX   1140
CD            This routine copies a vector, X, to a vector, Y.          PRONYX   1141
CD            Uses unrolled loops for increments equal to one.          PRONYX   1142
CD                                                                      PRONYX   1143
CD            CALLING PARAMETERS :                                      PRONYX   1144
CD                                                                      PRONYX   1145
CD                    N   - NO. OF POINTS IN THE ARRAY                  PRONYX   1146
CD                    DX  - INCOMING X-ARRAY                            PRONYX   1147
CD                  INCX  - INCREMENT OF X                              PRONYX   1148
CD                    DY  - INCOMING Y-ARRAY                            PRONYX   1149
CD                  INCY  - INCREMENT OF Y                              PRONYX   1150
CD                                                                      PRONYX   1151
C JACK DONGARRA, LINPACK, 3/11/78.                                      PRONYX   1152
        DOUBLE PRECISION DX(1),DY(1)                                    PRONYX   1153
        INTEGER I,INCX,INCY,IX,IY,M,MP1,N                               PRONYX   1154
        IF(N.LE.0)RETURN                                                PRONYX   1155
        IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20                             PRONYX   1156
C                                                                       PRONYX   1157
C code for unequal increments or equal increments not equal to 1        PRONYX   1158
C                                                                       PRONYX   1159
        IX = 1                                                          PRONYX   1160
        IY = 1                                                          PRONYX   1161
        IF(INCX.LT.0)IX = (-N+1)*INCX + 1                               PRONYX   1162
        IF(INCY.LT.0)IY = (-N+1)*INCY + 1                               PRONYX   1163
        DO 10 I = 1,N                                                   PRONYX   1164
          DY(IY) = DX(IX)                                               PRONYX   1165
          IX = IX + INCX                                                PRONYX   1166
          IY = IY + INCY                                                PRONYX   1167
10      CONTINUE                                                        PRONYX   1168
        RETURN                                                          PRONYX   1169
C                                                                       PRONYX   1170
C code for both increments equal to 1                                   PRONYX   1171
C                                                                       PRONYX   1172
C clean-up loop                                                         PRONYX   1173
C                                                                       PRONYX   1174
20      M = MOD(N,7)                                                    PRONYX   1175
        IF( M .EQ. 0 ) GO TO 40                                         PRONYX   1176
        DO 30 I = 1,M                                                   PRONYX   1177
          DY(I) = DX(I)                                                 PRONYX   1178
30      CONTINUE                                                        PRONYX   1179
        IF( N .LT. 7 ) RETURN                                           PRONYX   1180
40      MP1 = M + 1                                                     PRONYX   1181
        DO 50 I = MP1,N,7                                               PRONYX   1182
          DY(I) = DX(I)                                                 PRONYX   1183
          DY(I + 1) = DX(I + 1)                                         PRONYX   1184
          DY(I + 2) = DX(I + 2)                                         PRONYX   1185
          DY(I + 3) = DX(I + 3)                                         PRONYX   1186
          DY(I + 4) = DX(I + 4)                                         PRONYX   1187
          DY(I + 5) = DX(I + 5)                                         PRONYX   1188
          DY(I + 6) = DX(I + 6)                                         PRONYX   1189
50      CONTINUE                                                        PRONYX   1190
        RETURN                                                          PRONYX   1191
        END                                                             PRONYX   1192
C                                                                       PRONYX   1193
        DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)               PRONYX   1194
CD                                                                      PRONYX   1195
CD            The routine forms the dot product of two vectors.         PRONYX   1196
CD            Uses unrolled loops for increments equal to one.          PRONYX   1197
CD                                                                      PRONYX   1198
CD            CALLING PARAMETERS :                                      PRONYX   1199
CD                                                                      PRONYX   1200
CD                    N   - NO. OF POINTS IN THE ARRAY                  PRONYX   1201
CD                    DX  - INCOMING X-ARRAY                            PRONYX   1202
CD                  INCX  - INCREMENT OF X                              PRONYX   1203
CD                    DY  - INCOMING Y-ARRAY                            PRONYX   1204
CD                  INCY  - INCREMENT OF Y                              PRONYX   1205
CD                                                                      PRONYX   1206
C JACK DONGARRA, LINPACK, 3/11/78.                                      PRONYX   1207
        DOUBLE PRECISION DX(1),DY(1),DTEMP                              PRONYX   1208
        INTEGER I,INCX,INCY,IX,IY,M,MP1,N                               PRONYX   1209
        DDOT = 0.0D0                                                    PRONYX   1210
        DTEMP = 0.0D0                                                   PRONYX   1211
        IF(N.LE.0)RETURN                                                PRONYX   1212
        IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20                             PRONYX   1213
C                                                                       PRONYX   1214
C code for unequal increments or equal increments not equal to 1        PRONYX   1215
C                                                                       PRONYX   1216
        IX = 1                                                          PRONYX   1217
        IY = 1                                                          PRONYX   1218
        IF(INCX.LT.0)IX = (-N+1)*INCX + 1                               PRONYX   1219
        IF(INCY.LT.0)IY = (-N+1)*INCY + 1                               PRONYX   1220
        DO 10 I = 1,N                                                   PRONYX   1221
          DTEMP = DTEMP + DX(IX)*DY(IY)                                 PRONYX   1222
          IX = IX + INCX                                                PRONYX   1223
          IY = IY + INCY                                                PRONYX   1224
10      CONTINUE                                                        PRONYX   1225
        DDOT = DTEMP                                                    PRONYX   1226
        RETURN                                                          PRONYX   1227
C                                                                       PRONYX   1228
C code for both increments equal to 1                                   PRONYX   1229
C                                                                       PRONYX   1230
C clean-up loop                                                         PRONYX   1231
C                                                                       PRONYX   1232
20      M = MOD(N,5)                                                    PRONYX   1233
        IF( M .EQ. 0 ) GO TO 40                                         PRONYX   1234
        DO 30 I = 1,M                                                   PRONYX   1235
          DTEMP = DTEMP + DX(I)*DY(I)                                   PRONYX   1236
30      CONTINUE                                                        PRONYX   1237
        IF( N .LT. 5 ) GO TO 60                                         PRONYX   1238
40      MP1 = M + 1                                                     PRONYX   1239
        DO 50 I = MP1,N,5                                               PRONYX   1240
          DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +           PRONYX   1241
     *   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)PRONYX   1242
50      CONTINUE                                                        PRONYX   1243
60      DDOT = DTEMP                                                    PRONYX   1244
        RETURN                                                          PRONYX   1245
        END                                                             PRONYX   1246
        DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX)                  PRONYX   1247
CD                                                                      PRONYX   1248
CD            This rotunie calculates  Euclidean norm of the N-VECTOR   PRONYX   1249
CD            stored in DX() with storage include INCX                  PRONYX   1250
CD                                                                      PRONYX   1251
CD            CALLING PARAMETERS :                                      PRONYX   1252
CD                                                                      PRONYX   1253
CD                    N   - NO. OF POINTS IN THE ARRAY                  PRONYX   1254
CD                    DX  - INCOMING X-ARRAY                            PRONYX   1255
CD                  INCX  - INCREMENT OF X                              PRONYX   1256
CD                                                                      PRONYX   1257
CD                  If N .LE. 0 return with RESULT = 0.                 PRONYX   1258
CD                  If N .GE. 1 then INCX must be .GE. 1                PRONYX   1259
CD                                                                      PRONYX   1260
C                                                                       PRONYX   1261
C C.L.LAWSON, 1978 JAN 08                                               PRONYX   1262
C                                                                       PRONYX   1263
C Four phase method : using two built-in constants that are hopefully   PRONYX   1264
C applicable to all machines.                                           PRONYX   1265
C  CUTLO = MAXIMUM OF  DSQRT(U/EPS)  OVER ALL KNOWN MACHINES.           PRONYX   1266
C  CUTHI = MINIMUM OF  DSQRT(V)      OVER ALL KNOWN MACHINES.           PRONYX   1267
C where                                                                 PRONYX   1268
C  EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.                        PRONYX   1269
C  U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)                      PRONYX   1270
C  V   = LARGEST  NO.            (OVERFLOW  LIMIT)                      PRONYX   1271
C                                                                       PRONYX   1272
C Brief outline of algorithm..                                          PRONYX   1273
C                                                                       PRONYX   1274
C PHASE 1    scans zero components.                                     PRONYX   1275
C move TO PHASE 2 when a component is nonzero and .LE. CUTLO            PRONYX   1276
C move TO PHASE 3 when a component is .GT. CUTLO                        PRONYX   1277
C move TO PHASE 4 when a component is .GE. CUTHI/M                      PRONYX   1278
C where M = N for X() real and M = 2*N for complex                      PRONYX   1279
C                                                                       PRONYX   1280
C Values for CUTLO and CUTHI..                                          PRONYX   1281
C From the environmental parameters listed in the IMSL converter        PRONYX   1282
C document the limiting values are as follows..                         PRONYX   1283
C CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS AR     PRONYX   1284
C               UNIVAC AND DEC AT 2**(-103)                             PRONYX   1285
C               THUS CUTLO = 2**(-51) = 4.44089E-16                     PRONYX   1286
C CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.              PRONYX   1287
C               THUS CUTHI = 2**(63.5) = 1.30438E19                     PRONYX   1288
C CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.                 PRONYX   1289
C               THUS CUTLO = 2**(-33.5) = 8.23181D-11                   PRONYX   1290
C CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19                        PRONYX   1291
C DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /                            PRONYX   1292
C DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /                            PRONYX   1293
        INTEGER NEXT                                                    PRONYX   1294
        DOUBLE PRECISION DX(1),CUTLO,CUTHI,HITEST,SUM,XMAX,ZERO,ONE     PRONYX   1295
        DATA ZERO, ONE /0.0D0, 1.0D0/                                   PRONYX   1296
        DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /                      PRONYX   1297
C                                                                       PRONYX   1298
        IF(N .GT. 0) GO TO 10                                           PRONYX   1299
           DNRM2  = ZERO                                                PRONYX   1300
           GO TO 300                                                    PRONYX   1301
C                                                                       PRONYX   1302
10      ASSIGN 30 TO NEXT                                               PRONYX   1303
        SUM = ZERO                                                      PRONYX   1304
        NN = N * INCX                                                   PRONYX   1305
C begin main loop                                                       PRONYX   1306
        I = 1                                                           PRONYX   1307
20      GO TO NEXT,(30, 50, 70, 110)                                    PRONYX   1308
30      IF( DABS(DX(I)) .GT. CUTLO) GO TO 85                            PRONYX   1309
        ASSIGN 50 TO NEXT                                               PRONYX   1310
        XMAX = ZERO                                                     PRONYX   1311
C                                                                       PRONYX   1312
C PHASE 1. sum is zero                                                  PRONYX   1313
C                                                                       PRONYX   1314
50      IF( DX(I) .EQ. ZERO) GO TO 200                                  PRONYX   1315
        IF( DABS(DX(I)) .GT. CUTLO) GO TO 85                            PRONYX   1316
C                                                                       PRONYX   1317
C prepare for PHASE 2.                                                  PRONYX   1318
        ASSIGN 70 TO NEXT                                               PRONYX   1319
        GO TO 105                                                       PRONYX   1320
C                                                                       PRONYX   1321
C prepare for PHASE 4.                                                  PRONYX   1322
C                                                                       PRONYX   1323
100     I = J                                                           PRONYX   1324
        ASSIGN 110 TO NEXT                                              PRONYX   1325
        SUM = (SUM / DX(I)) / DX(I)                                     PRONYX   1326
105     XMAX = DABS(DX(I))                                              PRONYX   1327
        GO TO 115                                                       PRONYX   1328
C                                                                       PRONYX   1329
C PHASE 2.  sum is small. scale to avoid destructive underflow.         PRONYX   1330
C                                                                       PRONYX   1331
70      IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75                           PRONYX   1332
C                                                                       PRONYX   1333
C common code for PHASES 2 and 4.                                       PRONYX   1334
C In phase 4 SUM is large. Scale to avoid overflow                      PRONYX   1335
C                                                                       PRONYX   1336
110     IF( DABS(DX(I)) .LE. XMAX ) GO TO 115                           PRONYX   1337
           SUM = ONE + SUM * (XMAX / DX(I))**2                          PRONYX   1338
           XMAX = DABS(DX(I))                                           PRONYX   1339
           GO TO 200                                                    PRONYX   1340
115     SUM = SUM + (DX(I)/XMAX)**2                                     PRONYX   1341
        GO TO 200                                                       PRONYX   1342
C                                                                       PRONYX   1343
C Prepare for PHASE 3.                                                  PRONYX   1344
C                                                                       PRONYX   1345
75      SUM = (SUM * XMAX) * XMAX                                       PRONYX   1346
C                                                                       PRONYX   1347
C For REAL or D.P. set HITEST = CUTHI/N                                 PRONYX   1348
C For COMPLEX      set HITEST = CUTHI/(2*N)                             PRONYX   1349
C                                                                       PRONYX   1350
85      HITEST = CUTHI/REAL( N )                                        PRONYX   1351
C                                                                       PRONYX   1352
C PHASE 3.  Sum is mid-range.  no scaling.                              PRONYX   1353
C                                                                       PRONYX   1354
        DO 95 J =I,NN,INCX                                              PRONYX   1355
        IF(DABS(DX(J)) .GE. HITEST) GO TO 100                           PRONYX   1356
95      SUM = SUM + DX(J)**2                                            PRONYX   1357
        DNRM2 = DSQRT( SUM )                                            PRONYX   1358
        GO TO 300                                                       PRONYX   1359
200     CONTINUE                                                        PRONYX   1360
        I = I + INCX                                                    PRONYX   1361
        IF ( I .LE. NN ) GO TO 20                                       PRONYX   1362
C                                                                       PRONYX   1363
C End of main loop.                                                     PRONYX   1364
C                                                                       PRONYX   1365
C Compute square root and adjust for scaling.                           PRONYX   1366
C                                                                       PRONYX   1367
        DNRM2 = XMAX * DSQRT(SUM)                                       PRONYX   1368
300     CONTINUE                                                        PRONYX   1369
        RETURN                                                          PRONYX   1370
        END                                                             PRONYX   1371
        SUBROUTINE DROT (N,DX,INCX,DY,INCY,C,S)                         PRONYX   1372
CD                                                                      PRONYX   1373
CD       This routine Applies a plane rotation.                         PRONYX   1374
CD       Uses unrolled loops for increments equal to one.               PRONYX   1375
CD                                                                      PRONYX   1376
CD            CALLING PARAMETERS :                                      PRONYX   1377
CD                                                                      PRONYX   1378
CD                    N   - NO. OF POINTS IN THE ARRAY                  PRONYX   1379
CD                    DX  - INCOMING X-ARRAY                            PRONYX   1380
CD                  INCX  - INCREMENT OF X                              PRONYX   1381
CD                    DY  - INCOMING Y-ARRAY                            PRONYX   1382
CD                  INCY  - INCREMENT OF Y                              PRONYX   1383
CD                     C  - EQUIVALENT TO COSINE(ANGLE OF ROTATION)     PRONYX   1384
CD                     S  - EQUIVALENT TO SINE(ANGLE OF ROTATION)       PRONYX   1385
CD                                                                      PRONYX   1386
C      JACK DONGARRA, LINPACK, 3/11/78.                                 PRONYX   1387
C                                                                       PRONYX   1388
        DOUBLE PRECISION DX(1),DY(1),DTEMP,C,S                          PRONYX   1389
        INTEGER I,INCX,INCY,IX,IY,N                                     PRONYX   1390
        IF(N.LE.0)RETURN                                                PRONYX   1391
        IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20                             PRONYX   1392
C                                                                       PRONYX   1393
C Code for unequal increments or equal increments not equal to 1        PRONYX   1394
C                                                                       PRONYX   1395
        IX = 1                                                          PRONYX   1396
        IY = 1                                                          PRONYX   1397
        IF(INCX.LT.0)IX = (-N+1)*INCX + 1                               PRONYX   1398
        IF(INCY.LT.0)IY = (-N+1)*INCY + 1                               PRONYX   1399
        DO 10 I = 1,N                                                   PRONYX   1400
          DTEMP = C*DX(IX) + S*DY(IY)                                   PRONYX   1401
          DY(IY) = C*DY(IY) - S*DX(IX)                                  PRONYX   1402
          DX(IX) = DTEMP                                                PRONYX   1403
          IX = IX + INCX                                                PRONYX   1404
          IY = IY + INCY                                                PRONYX   1405
10      CONTINUE                                                        PRONYX   1406
        RETURN                                                          PRONYX   1407
C                                                                       PRONYX   1408
C Code for both increments equal to 1                                   PRONYX   1409
C                                                                       PRONYX   1410
20      DO 30 I = 1,N                                                   PRONYX   1411
          DTEMP = C*DX(I) + S*DY(I)                                     PRONYX   1412
          DY(I) = C*DY(I) - S*DX(I)                                     PRONYX   1413
          DX(I) = DTEMP                                                 PRONYX   1414
30      CONTINUE                                                        PRONYX   1415
        RETURN                                                          PRONYX   1416
        END                                                             PRONYX   1417
        SUBROUTINE DROTG(DA,DB,C,S)                                     PRONYX   1418
CD                                                                      PRONYX   1419
CD       This routine construct givens plane rotation.                  PRONYX   1420
CD                                                                      PRONYX   1421
CD       CALLING PARAMETERS :                                           PRONYX   1422
CD                                                                      PRONYX   1423
CD                    DA,DB - INCOMING VARIABLES WHICH IS USED TO       PRONYX   1424
CD                            CALCULATE PLANE ROTATION                  PRONYX   1425
CD                        C - ROTATION FACTOR IN THE X-AXIS             PRONYX   1426
CD                        S - ROTATION FACTOR IN THE Y-AXIS             PRONYX   1427
CD                                                                      PRONYX   1428
C           JACK DONGARRA, LINPACK, 3/11/                               PRONYX   1429
        DOUBLE PRECISION DA,DB,C,S,ROE,SCALE,R,Z                        PRONYX   1430
        ROE = DB                                                        PRONYX   1431
        IF( DABS(DA) .GT. DABS(DB) ) ROE = DA                           PRONYX   1432
        SCALE = DABS(DA) + DABS(DB)                                     PRONYX   1433
        IF( SCALE .NE. 0.0D0 ) GO TO 10                                 PRONYX   1434
           C = 1.0D0                                                    PRONYX   1435
           S = 0.0D0                                                    PRONYX   1436
           R = 0.0D0                                                    PRONYX   1437
           GO TO 20                                                     PRONYX   1438
10      R = SCALE*DSQRT((DA/SCALE)**2 + (DB/SCALE)**2)                  PRONYX   1439
        R = DSIGN(1.0D0,ROE)*R                                          PRONYX   1440
        C = DA/R                                                        PRONYX   1441
        S = DB/R                                                        PRONYX   1442
20      Z = 1.0D0                                                       PRONYX   1443
        IF( DABS(DA) .GT. DABS(DB) ) Z = S                              PRONYX   1444
        IF( DABS(DB) .GE. DABS(DA) .AND. C .NE. 0.0D0 ) Z = 1.0D0/C     PRONYX   1445
        DA = R                                                          PRONYX   1446
        DB = Z                                                          PRONYX   1447
        RETURN                                                          PRONYX   1448
        END                                                             PRONYX   1449
C                                                                       PRONYX   1450
        SUBROUTINE DSCAL(N,DA,DX,INCX)                                  PRONYX   1451
CD                                                                      PRONYX   1452
CD             This routine scales a vector by a constant.              PRONYX   1453
CD             Uses unrolled loops for increment equal to one.          PRONYX   1454
CD                                                                      PRONYX   1455
CD            CALLING PARAMETERS :                                      PRONYX   1456
CD                                                                      PRONYX   1457
CD                    N   - NO. OF POINTS IN THE ARRAY                  PRONYX   1458
CD                    DA  - INCOMING CONSTANT                           PRONYX   1459
CD                    DX  - INCOMING X-ARRAY                            PRONYX   1460
CD                  INCX  - INCREMENT OF X                              PRONYX   1461
CD                                                                      PRONYX   1462
C JACK DONGARRA, LINPACK, 3/11/78.                                      PRONYX   1463
        DOUBLE PRECISION DA,DX(1)                                       PRONYX   1464
        INTEGER I,INCX,M,MP1,N,NINCX                                    PRONYX   1465
        IF(N.LE.0)RETURN                                                PRONYX   1466
        IF(INCX.EQ.1)GO TO 20                                           PRONYX   1467
C                                                                       PRONYX   1468
C Code for increment not equal to 1                                     PRONYX   1469
C                                                                       PRONYX   1470
        NINCX = N*INCX                                                  PRONYX   1471
        DO 10 I = 1,NINCX,INCX                                          PRONYX   1472
          DX(I) = DA*DX(I)                                              PRONYX   1473
10      CONTINUE                                                        PRONYX   1474
        RETURN                                                          PRONYX   1475
C                                                                       PRONYX   1476
C Code for increment equal to 1                                         PRONYX   1477
C                                                                       PRONYX   1478
C Clean-up loop                                                         PRONYX   1479
C                                                                       PRONYX   1480
20      M = MOD(N,5)                                                    PRONYX   1481
        IF( M .EQ. 0 ) GO TO 40                                         PRONYX   1482
        DO 30 I = 1,M                                                   PRONYX   1483
          DX(I) = DA*DX(I)                                              PRONYX   1484
30      CONTINUE                                                        PRONYX   1485
        IF( N .LT. 5 ) RETURN                                           PRONYX   1486
40      MP1 = M + 1                                                     PRONYX   1487
        DO 50 I = MP1,N,5                                               PRONYX   1488
          DX(I) = DA*DX(I)                                              PRONYX   1489
          DX(I + 1) = DA*DX(I + 1)                                      PRONYX   1490
          DX(I + 2) = DA*DX(I + 2)                                      PRONYX   1491
          DX(I + 3) = DA*DX(I + 3)                                      PRONYX   1492
          DX(I + 4) = DA*DX(I + 4)                                      PRONYX   1493
50      CONTINUE                                                        PRONYX   1494
        RETURN                                                          PRONYX   1495
        END                                                             PRONYX   1496
C                                                                       PRONYX   1497
        SUBROUTINE DSWAP (N,DX,INCX,DY,INCY)                            PRONYX   1498
CD                                                                      PRONYX   1499
CD             This routine Interchanges two vectors.                   PRONYX   1500
CD             Uses unrolled loops for increments equal one.            PRONYX   1501
CD                                                                      PRONYX   1502
CD            CALLING PARAMETERS :                                      PRONYX   1503
CD                                                                      PRONYX   1504
CD                    N   - NO. OF POINTS IN THE ARRAY                  PRONYX   1505
CD                    DX  - INCOMING X-ARRAY                            PRONYX   1506
CD                  INCX  - INCREMENT OF X                              PRONYX   1507
CD                    DY  - INCOMING Y-ARRAY                            PRONYX   1508
CD                  INCY  - INCREMENT OF Y                              PRONYX   1509
CD                                                                      PRONYX   1510
C              JACK DONGARRA, LINPACK, 3/11/78.                         PRONYX   1511
        DOUBLE PRECISION DX(1),DY(1),DTEMP                              PRONYX   1512
        INTEGER I,INCX,INCY,IX,IY,M,MP1,N                               PRONYX   1513
        IF(N.LE.0)RETURN                                                PRONYX   1514
        IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20                             PRONYX   1515
C                                                                       PRONYX   1516
C Code for unequal increments or equal increments not equal to 1        PRONYX   1517
C                                                                       PRONYX   1518
        IX = 1                                                          PRONYX   1519
        IY = 1                                                          PRONYX   1520
        IF(INCX.LT.0)IX = (-N+1)*INCX + 1                               PRONYX   1521
        IF(INCY.LT.0)IY = (-N+1)*INCY + 1                               PRONYX   1522
        DO 10 I = 1,N                                                   PRONYX   1523
          DTEMP = DX(IX)                                                PRONYX   1524
          DX(IX) = DY(IY)                                               PRONYX   1525
          DY(IY) = DTEMP                                                PRONYX   1526
          IX = IX + INCX                                                PRONYX   1527
          IY = IY + INCY                                                PRONYX   1528
10      CONTINUE                                                        PRONYX   1529
        RETURN                                                          PRONYX   1530
C                                                                       PRONYX   1531
C Code for both increments equal to 1                                   PRONYX   1532
C                                                                       PRONYX   1533
C Clean-up loop                                                         PRONYX   1534
C                                                                       PRONYX   1535
20      M = MOD(N,3)                                                    PRONYX   1536
        IF( M .EQ. 0 ) GO TO 40                                         PRONYX   1537
        DO 30 I = 1,M                                                   PRONYX   1538
          DTEMP = DX(I)                                                 PRONYX   1539
          DX(I) = DY(I)                                                 PRONYX   1540
          DY(I) = DTEMP                                                 PRONYX   1541
30      CONTINUE                                                        PRONYX   1542
        IF( N .LT. 3 ) RETURN                                           PRONYX   1543
40      MP1 = M + 1                                                     PRONYX   1544
        DO 50 I = MP1,N,3                                               PRONYX   1545
          DTEMP = DX(I)                                                 PRONYX   1546
          DX(I) = DY(I)                                                 PRONYX   1547
          DY(I) = DTEMP                                                 PRONYX   1548
          DTEMP = DX(I + 1)                                             PRONYX   1549
          DX(I + 1) = DY(I + 1)                                         PRONYX   1550
          DY(I + 1) = DTEMP                                             PRONYX   1551
          DTEMP = DX(I + 2)                                             PRONYX   1552
          DX(I + 2) = DY(I + 2)                                         PRONYX   1553
          DY(I + 2) = DTEMP                                             PRONYX   1554
50      CONTINUE                                                        PRONYX   1555
        RETURN                                                          PRONYX   1556
        END                                                             PRONYX   1557
C                                                                       PRONYX   1558
        INTEGER FUNCTION IDAMAX(N,DX,INCX)                              PRONYX   1559
CD                                                                      PRONYX   1560
CD       This routine finds the index of element having maximum         PRONYX   1561
CD       absolute value.                                                PRONYX   1562
CD                                                                      PRONYX   1563
CD            CALLING PARAMETERS :                                      PRONYX   1564
CD                                                                      PRONYX   1565
CD                    N   - NO. OF POINTS IN THE ARRAY                  PRONYX   1566
CD                    DX  - INCOMING X-ARRAY                            PRONYX   1567
CD                  INCX  - INCREMENT OF X                              PRONYX   1568
CD                                                                      PRONYX   1569
C         JACK DONGARRA, LINPACK, 3/11/78.                              PRONYX   1570
        DOUBLE PRECISION DX(1),DMAX                                     PRONYX   1571
        INTEGER I,INCX,IX,N                                             PRONYX   1572
        IDAMAX = 0                                                      PRONYX   1573
        IF( N .LT. 1 ) RETURN                                           PRONYX   1574
        IDAMAX = 1                                                      PRONYX   1575
        IF(N.EQ.1)RETURN                                                PRONYX   1576
        IF(INCX.EQ.1)GO TO 20                                           PRONYX   1577
C                                                                       PRONYX   1578
C Code for increment not equal to 1                                     PRONYX   1579
C                                                                       PRONYX   1580
        IX = 1                                                          PRONYX   1581
        DMAX = DABS(DX(1))                                              PRONYX   1582
        IX = IX + INCX                                                  PRONYX   1583
        DO 10 I = 2,N                                                   PRONYX   1584
           IF(DABS(DX(IX)).LE.DMAX) GO TO 5                             PRONYX   1585
           IDAMAX = I                                                   PRONYX   1586
           DMAX = DABS(DX(IX))                                          PRONYX   1587
5          IX = IX + INCX                                               PRONYX   1588
10      CONTINUE                                                        PRONYX   1589
        RETURN                                                          PRONYX   1590
C                                                                       PRONYX   1591
C Code for increment equal to 1                                         PRONYX   1592
C                                                                       PRONYX   1593
20      DMAX = DABS(DX(1))                                              PRONYX   1594
        DO 30 I = 2,N                                                   PRONYX   1595
           IF(DABS(DX(I)).LE.DMAX) GO TO 30                             PRONYX   1596
           IDAMAX = I                                                   PRONYX   1597
           DMAX = DABS(DX(I))                                           PRONYX   1598
30      CONTINUE                                                        PRONYX   1599
        RETURN                                                          PRONYX   1600
        END                                                             PRONYX   1601
                                                                        PRONYX   1602
C                                                                       PRONYX   1603
        SUBROUTINE DSQSVD(A,U,V,D,C,PIRCON,N1,LPOCON,LUOUT,BATCH,NOGO)  PRONYX   1604
CD           This routine computes the Least Squares solution to        PRONYX   1605
CD           Y A = C  where a rank PIRCON version of the width P        PRONYX   1606
CD           matrix Y is used, computed using a singular value decompn  PRONYX   1607
CD           Y=U Sigma Trans(V).                                        PRONYX   1608
CD           The diagonal matrix Sigma is stored in D. Only the first   PRONYX   1609
CD           PIRCON singular values are used in computing the least     PRONYX   1610
CD           sqs soln.                                                  PRONYX   1611
CD                                                                      PRONYX   1612
C LAST MODIFIED 07/05/89 [JFH]                                          PRONYX   1613
C                                                                       PRONYX   1614
CD  Calling Parameters :                                                PRONYX   1615
CD                                                                      PRONYX   1616
CD    Inputs :   U,D,V   : Matrices of the SVD of Y                     PRONYX   1617
CD               C       : Right hand side (already with sign changed)  PRONYX   1618
CD               PIRCON  : Forced rank limit of Y                       PRONYX   1619
CD               N1      : Number of rows of Y                          PRONYX   1620
CD               P       : Number of columns of Y [LPOCON]              PRONYX   1621
CD               BATCH   : FLAG FOR BATCH EXECUTION                     PRONYX   1622
CD    Outputs :                                                         PRONYX   1623
CD               A       : Solution vector [0,...,P], A(0)=1.           PRONYX   1624
        INTEGER NSIZE ,NTAIL                                            PRONYX   1625
        PARAMETER (NSIZE=128)                                           PRONYX   1626
        PARAMETER (NTAIL=1024)                                          PRONYX   1627
        INTEGER LPOCON,PIRCON,N1,I,J,LUOUT                              PRONYX   1628
        DOUBLE PRECISION U(1:NTAIL,1:NSIZE),V(1:NSIZE,1:NSIZE)          PRONYX   1629
        DOUBLE PRECISION A(0:*),C(1:*),D(1:*),E(1:NSIZE),ALPHA          PRONYX   1630
        DOUBLE PRECISION EPS,DTEST                                      PRONYX   1631
        LOGICAL BATCH,TRACK                                             PRONYX   1632
        INTEGER NOGO                                                    PRONYX   1633
C                                                                       PRONYX   1634
      TRACK=.TRUE.                                                      PRONYX   1635
C                                                                       PRONYX   1636
C *** DIAGNOSTIC                                                        PRONYX   1637
      IF (TRACK) THEN                                                   PRONYX   1638
C#      IF (BATCH) WRITE (*,9000) PIRCON,N1,LPOCON                      PRONYX   1639
C#        WRITE (LUOUT,9000) PIRCON,N1,LPOCON                           PRONYX   1640
      END IF                                                            PRONYX   1641
C                                                                       PRONYX   1642
C *** EXAMINE SINGULAR VALUES                                           PRONYX   1643
      EPS=1.D-10                                                        PRONYX   1644
      DO 10 I=2,LPOCON                                                  PRONYX   1645
        IF (D(I).GT.D(I-1)) THEN                                        PRONYX   1646
          IF (BATCH) WRITE (*,9020) I                                   PRONYX   1647
          WRITE (LUOUT,9020) I                                          PRONYX   1648
        END IF                                                          PRONYX   1649
        DTEST=D(I)/D(1)                                                 PRONYX   1650
        IF (DTEST.LT.EPS) THEN                                          PRONYX   1651
          IF (BATCH) WRITE (*,9040) I,DTEST                             PRONYX   1652
          WRITE (LUOUT,9040) I,DTEST                                    PRONYX   1653
          IF (I.GE.PIRCON) THEN                                         PRONYX   1654
            PIRCON=I-1                                                  PRONYX   1655
C#          IF (BATCH) WRITE (*,9042) PIRCON                            PRONYX   1656
C#            WRITE (LUOUT,9042) PIRCON                                 PRONYX   1657
          END IF                                                        PRONYX   1658
          GO TO 200                                                     PRONYX   1659
        END IF                                                          PRONYX   1660
 10   CONTINUE                                                          PRONYX   1661
C                                                                       PRONYX   1662
C *** COMPUTE E = LOW RANK INVERSE (SIGMA) * TRANS(U) * C               PRONYX   1663
  200 CONTINUE                                                          PRONYX   1664
C FLF                                                                   PRONYX   1665
      IF (DTEST.LT.EPS) THEN                                            PRONYX   1666
          NOGO=-1                                                       PRONYX   1667
          RETURN                                                        PRONYX   1668
      ENDIF                                                             PRONYX   1669
C FLF                                                                   PRONYX   1670
      DO 30 I=1,PIRCON                                                  PRONYX   1671
        ALPHA = 0.D0                                                    PRONYX   1672
        DO 20 J=1,N1                                                    PRONYX   1673
          ALPHA = ALPHA + U(J,I) * C(J)                                 PRONYX   1674
 20     CONTINUE                                                        PRONYX   1675
        E(I) = ALPHA / D(I)                                             PRONYX   1676
 30   CONTINUE                                                          PRONYX   1677
      DO 40 I=PIRCON+1,LPOCON                                           PRONYX   1678
        E(I) = 0.D0                                                     PRONYX   1679
 40   CONTINUE                                                          PRONYX   1680
C                                                                       PRONYX   1681
C *** COMPUTE A = V * E                                                 PRONYX   1682
      A(0) = 1.D0                                                       PRONYX   1683
      DO 60 I=1,LPOCON                                                  PRONYX   1684
        ALPHA = 0.D0                                                    PRONYX   1685
        DO 50 J=1,LPOCON                                                PRONYX   1686
          ALPHA = ALPHA + V(I,J) * E(J)                                 PRONYX   1687
 50     CONTINUE                                                        PRONYX   1688
        A(I) = ALPHA                                                    PRONYX   1689
 60   CONTINUE                                                          PRONYX   1690
C                                                                       PRONYX   1691
C 9000 FORMAT (2X,'#IN DLEASTSVD: PIRCON,N1,LPOCON =',3I4)              PRONYX   1692
 9020 FORMAT (2X,'#IN DLEASTSVD: SIGMA',I3,' OUT OF ORDER')             PRONYX   1693
 9040 FORMAT (2X,'#IN DLEASTSVD: SIGMA',I3,' TOO SMALL: DTEST = ',D15.7)PRONYX   1694
C 9042 FORMAT (2X,'#IN DLEASTSVD: REDUCING PIRCON TO',I4)               PRONYX   1695
C                                                                       PRONYX   1696
      RETURN                                                            PRONYX   1697
      END                                                               PRONYX   1698
      SUBROUTINE DSQVDM (X,ROOTR,ROOTI,MAG,PHA,AMP,PHI,N,Q,IND,IC       PRONYX   1699
     1,ERR,NOGO)                                                        PRONYX   1700
CD            This routine forms the real version of the complex        PRONYX   1701
CD            Van Der Monde system M containing the modes.              PRONYX   1702
CD            Solve the Least Squares problem M D = A, for D            PRONYX   1703
CD            where A is a real vector, containing the data. The soln   PRONYX   1704
CD            as well as the residual vector are computed using the     PRONYX   1705
CD            LINPACK routines and SQRSL.                               PRONYX   1706
CD            Use the soln D to get the amplitude and initial phase     PRONYX   1707
CD            of the cosines.                                           PRONYX   1708
C                                                                       PRONYX   1709
C Date : Feb 1 1988                                                     PRONYX   1710
CD                                                                      PRONYX   1711
CD Calling Parameters :                                                 PRONYX   1712
CD                                                                      PRONYX   1713
CD Inputs :      X       : Data vector [0,...,N-1]                      PRONYX   1714
CD               MAG     : Magnitude of the modes [0,...,Q]             PRONYX   1715
CD               PHA     : Phase of the modes [0,...,Q] (in degrees)    PRONYX   1716
CD               N       : Number of data values                        PRONYX   1717
CD               Q       : Number of modes                              PRONYX   1718
CD               IND     : Total number of available modes              PRONYX   1719
CD               IC      : Binary array,indicate posn of desired modes  PRONYX   1720
CD Outputs :                                                            PRONYX   1721
CD               AMP     : Amplitude of the modes [0,...,Q].            PRONYX   1722
CD               PHI     : Initial phase of the modes [0,...,Q].        PRONYX   1723
CD               ERR     : Sum of Squared residuals (prediction error). PRONYX   1724
CD                                                                      PRONYX   1725
        INTEGER NSIZE ,NTAIL                                            PRONYX   1726
        PARAMETER (NSIZE=128)                                           PRONYX   1727
        PARAMETER (NTAIL=1024)                                          PRONYX   1728
        INTEGER N ,Q ,I ,J ,INFO ,DEGREE ,IND ,IC(1:*) ,JDUM(1)         PRONYX   1729
        DOUBLE PRECISION X(0:*) ,ROOTR(1:*) ,ROOTI(1:*) ,MAG(1:*) ,     PRONYX   1730
     +                   PHA(1:*),DUM(1)                                PRONYX   1731
        DOUBLE PRECISION EPS ,AMP(1:*) ,PHI(1:*) ,PI ,ERR               PRONYX   1732
        DOUBLE PRECISION MAT(1:NTAIL,1:NSIZE) ,QRAUX(1:NSIZE)           PRONYX   1733
        DOUBLE PRECISION A(1:NTAIL) ,B(1:NTAIL) ,D(1:NSIZE) ,AA         PRONYX   1734
        INTEGER NOGO                                                    PRONYX   1735
        COMMON / DSQCOM /A,B,D,QRAUX,MAT                                PRONYX   1736
C                                                                       PRONYX   1737
        PI = 4.0D0 * DATAN (1.0D0)                                      PRONYX   1738
        EPS = 1.0D-12                                                   PRONYX   1739
C                                                                       PRONYX   1740
C Set the "Van der Monde" least squares system                          PRONYX   1741
        CALL DSTVDM (MAT,MAG,PHA,N,Q,IND,IC,DEGREE,NOGO)                PRONYX   1742
C FLF                                                                   PRONYX   1743
        IF(NOGO.LT.0) RETURN                                            PRONYX   1744
C FLF                                                                   PRONYX   1745
C                                                                       PRONYX   1746
C Compute its QR factorization                                          PRONYX   1747
        CALL DQRDC (MAT,NTAIL,N,DEGREE,QRAUX,JDUM,DUM,0)                PRONYX   1748
C                                                                       PRONYX   1749
C Set up the right hand side                                            PRONYX   1750
        DO 10 I=0,N-1                                                   PRONYX   1751
          A(I+1) = X(I)                                                 PRONYX   1752
 10     CONTINUE                                                        PRONYX   1753
C                                                                       PRONYX   1754
C Compute the least square solution vector D and the residual vector B  PRONYX   1755
        CALL DQRSL (MAT,NTAIL,N,DEGREE,QRAUX,A,DUM,B,D,B,DUM,110,INFO)  PRONYX   1756
C                                                                       PRONYX   1757
C Compute the residual sum of squares                                   PRONYX   1758
        ERR = 0.D0                                                      PRONYX   1759
        DO 20 I=1,N                                                     PRONYX   1760
          ERR = ERR + B(I) * B(I)                                       PRONYX   1761
 20     CONTINUE                                                        PRONYX   1762
C                                                                       PRONYX   1763
C From those complex amplitude, compute the power and the initial phase PRONYX   1764
        J = 1                                                           PRONYX   1765
        DO 30 I=1,IND                                                   PRONYX   1766
          IF (IC(I).EQ.1) THEN                                          PRONYX   1767
C Treat the roots on the real axis separately                           PRONYX   1768
            AA = MAG(I) * SIN(PI * PHA(I) / 180.D0)                     PRONYX   1769
            IF (DABS(AA).LT.EPS) THEN                                   PRONYX   1770
              PHI(I) = 0.D0                                             PRONYX   1771
              AMP(I) = D(J)                                             PRONYX   1772
              J = J + 1                                                 PRONYX   1773
            ELSE                                                        PRONYX   1774
              AMP(I) = DSQRT( D(J)*D(J) + D(J+1)*D(J+1))                PRONYX   1775
              PHI(I) = DATAN2( - D(J+1) , D(J) )                        PRONYX   1776
              PHI(I) = PHI(I) * 180.D0 / PI                             PRONYX   1777
              J = J + 2                                                 PRONYX   1778
            END IF                                                      PRONYX   1779
          END IF                                                        PRONYX   1780
 30     CONTINUE                                                        PRONYX   1781
C                                                                       PRONYX   1782
        RETURN                                                          PRONYX   1783
        END                                                             PRONYX   1784
      SUBROUTINE DPOLRT (ADDP,NDEG,ZERORQ,ZEROIQ,IFAIL)                 PRONYX   1785
CD                                                                      PRONYX   1786
CD     THIS IS A DUMMY ROUTINE ,THE ROOTING IS                          PRONYX   1787
CD     ACTUALLY VIA 'QPOLRT'                                            PRONYX   1788
CD                                                                      PRONYX   1789
*IF -VAX                                                                PRONYX   1790
CSE                                                                     PRONYX   1791
CPS      PURPOSE:                                                       PRONYX   1792
C                                                                       PRONYX   1793
C           FINDS THE ZEROS OF A REAL POLYNOMIAL OF THE FORM            PRONYX   1794
C                                                                       PRONYX   1795
C             X**N + A(N) * X**(N-1) + .... + A(2) * X + A(1) = 0       PRONYX   1796
C                                                                       PRONYX   1797
C           THAT IS, MONIC POLYNOMIALS ONLY ARE CONSIDERED.             PRONYX   1798
C                                                                       PRONYX   1799
C           (CAUTION: ROOT CRACKING FOR LARGE DEGREE POLYNOMIALS IS     PRONYX   1800
C            AN ILL-CONDITIONED COMPUTING PROBLEM!!!!)                  PRONYX   1801
C                                                                       PRONYX   1802
CPE                                                                     PRONYX   1803
CAS           ***  ARGUMENT LIST  ***                                   PRONYX   1804
C                                                                       PRONYX   1805
C  AD   -  COEFFICIENTS IN ORDER OF INCREASING POWER, QUAD PRECISION    PRONYX   1806
C                                                                       PRONYX   1807
C  NDEG -  INTEGER DEGREE OF POLYNOMIAL                                 PRONYX   1808
C                                                                       PRONYX   1809
C  ZEROR-  REAL PARTS OF ZEROS, DOUBLE PRECISION                        PRONYX   1810
C          (SEE REMARKS BELOW)                                          PRONYX   1811
C  ZEROI-  IMAGINARY PARTS OF ZEROS, DOUBLE PRECISION                   PRONYX   1812
C          (SEE REMARKS BELOW)                                          PRONYX   1813
C                                                                       PRONYX   1814
C  IFAIL-  OUTPUT FLAG, SET TO 0 IF NO TROUBLE, SET TO 1 IF FEWER THAN  PRONYX   1815
C          NDEG ROOTS ARE FOUND, IN WHICH CASE, NDEG IS RESET TO THE    PRONYX   1816
C          NUMBER OF ROOTS FOUND.  IFAIL IS SET TO -1 IF ALL ROOTS ARE  PRONYX   1817
C          FOUND, BUT ACCURACY OF RECONSTRUCTED COEFFICIENTS DOES NOT   PRONYX   1818
C          SPECIFIED TOLERANCE OF 'ERRTOL'.  SEE BELOW.                 PRONYX   1819
C                                                                       PRONYX   1820
C  THIS PROGRAM IS AN ADAPTATION OF SUBROUTINE RPOLY, DUE TO JENKINS    PRONYX   1821
C  AND TRAUB, AND GIVEN IN  'TRANSACTIONS ON MATHEMATICAL SOFTWARE?'    PRONYX   1822
C  VOLUME 1, !2 (JUNE, 1975).  THE PROGRAM HAS BEEN REWORKED TO PROVIDE PRONYX   1823
C  GREATER USER CONVENIENCE, WITH ESSENTIALLY NO CHANGE IN THE LOGICAL  PRONYX   1824
C  STRUCTURE.  THE EXCEPTION IS THAT THE PROGRAM DOES NOT TREAT         PRONYX   1825
C  UNNORMALIZED POLYNOMIALS ANY LONGER; HOWEVER, THE RESTRICTION TO     PRONYX   1826
C  MONIC POLYNOMIALS CAN BE REMOVED WITH A SIMPLE CHANGE TO THE SECTION PRONYX   1827
C  BELOW LABELLED 'MAKE A COPY OF THE COEFFICIENTS?'                    PRONYX   1828
C  THE MAIN CHANGES HAVE BEEN TO ELIMINATE ALL SUPPORTING SUBROUTINES   PRONYX   1829
C  IN FAVOR OF IN-LINE CODE (FOR THOSE SUBROUTINES WHICH WERE CALLED    PRONYX   1830
C  ONLY ONCE; FXSHFR, QUADIT, REALIT) OR SEPARATE SECTIONS OF CODE      PRONYX   1831
C  WHICH ARE BRANCHED OUT TO FROM THE MAIN SECTION.  THE ORIGINAL       PRONYX   1832
C  QUADRATIC SYNTHETIC DIVISION SUBROUTINE, QUADSD, WAS ALSO REPLACED   PRONYX   1833
C  WITH IN-LINE CODE BECAUSE IT IS SHORT AND WAS CALLED ONLY FIVE TIMES.PRONYX   1834
C  SINGLE PRECISION CALCULATIONS ARE USED FOR SCALING, BOUNDS, AND      PRONYX   1835
C  ERROR CALCULATIONS.  ALL ITERATIONS ARE DONE IN DOUBLE PRECISION.    PRONYX   1836
C                                                                       PRONYX   1837
C  IF ALL ROOTS ARE FOUND, THE PROGRAM CHECKS THE ACCURACY OF THE       PRONYX   1838
C  ROOTING PROCEDURE BY RECONSTRUCTING THE COEFFICIENTS OF POLYNOMIAL   PRONYX   1839
C  FROM THE ROOTS.  THE ERROR PARAMETER USED IS THE EUCLIDIAN NORM      PRONYX   1840
C  OF ERROR VECTOR (ORIGINAL COEFFICIENT MINUS RECONSTRUCTED) DIVIDED   PRONYX   1841
C  BY THE EUCLIDIAN NORM OF THE ORIGINAL COEFFICIENT VECTOR.  THIS      PRONYX   1842
C  SHOULD BE SMALLER THAN 'ERRTOL', OR IFAIL IS SET TO -1.              PRONYX   1843
C                                                                       PRONYX   1844
C  EXTRACTED FROM NASA 'SAMSAN' CODES AND REVISED TO QUAD PRECISION BY  PRONYX   1845
C  J.F. HAUER OF THE BONNEVILLE POWER ADMINISTRATION, PORTLAND, OREGON. PRONYX   1846
C  THIS IS A REVERSION TO DOUBLE PRECISION.                             PRONYX   1847
C                                                                       PRONYX   1848
C *** LAST JFH REVISION 10/26/88                                        PRONYX   1849
C                                                                       PRONYX   1850
C                                                                       PRONYX   1851
CAE                                                                     PRONYX   1852
      IMPLICIT DOUBLE PRECISION(A-H,O-Z),INTEGER(I-N)                   PRONYX   1853
      DOUBLE PRECISION ADDP,ZERORQ,ZEROIQ                               PRONYX   1854
      DIMENSION ADDP(1),ZERORQ(1),ZEROIQ(1)                             PRONYX   1855
      DOUBLE PRECISION K,KV,LZI,LZR                                     PRONYX   1856
      LOGICAL  SPASS,STRY,TRIED,VPASS,VTRY,ZEROK                        PRONYX   1857
      LOGICAL TRACK                                                     PRONYX   1858
      DIMENSION  AD(121),PT(121),ZEROR(121),ZEROI(121)                  PRONYX   1859
      DIMENSION  K(121),P(121),QK(121),QP(121),SVK(121),TEMP(121)       PRONYX   1860
      EQUIVALENCE  (PT(1),K(1))                                         PRONYX   1861
      DATA  ERRTOL/1.0D-07/                                             PRONYX   1862
      DATA NBITSD/55/,BASE/2.D0/,RINFIN/.17D+38/,SMALNO/.29D-37/        PRONYX   1863
C     DATA NBITSD/112/,BASE/2.Q0/,RINFIN/1.Q4930/,SMALNO/1.Q-4930/      PRONYX   1864
C                                                                       PRONYX   1865
C  THE DATA STATEMENT ABOVE AND THE NEXT FOUR STATEMENTS SET MACHINE-   PRONYX   1866
C  DEPENDENT CONSTANTS USED IN VARIOUS PARTS OF THE PROGRAM.  THE       PRONYX   1867
C  CONSTANTS ARE:                                                       PRONYX   1868
C                                                                       PRONYX   1869
C  NBITSD -  # OF SIGNIFICANT BITS IN THE EXTENDED-PRECISON MANTSSA OF  PRONYX   1870
C            MACHINE BEING USED.                                        PRONYX   1871
C                                                                       PRONYX   1872
C  BASE   -  BASE OF MACHINE'S FLOATING POINT NUMBER SYSTEM.            PRONYX   1873
C                                                                       PRONYX   1874
C  RINFIN -  A LARGE FLOATING POINT NUMBER NEAR THE UPPER END OF THE    PRONYX   1875
C            MACHINE'S PRECISION RANGE.                                 PRONYX   1876
C                                                                       PRONYX   1877
C  SMALNO -  A SMALL FLOATING POINT NUMBER NEAR THE LOWER END OF THE    PRONYX   1878
C            MACHINE'S PRECISION RANGE.                                 PRONYX   1879
C                                                                       PRONYX   1880
C  ETA    -  THE MAXIMUM RELATIVE REPRESENTATION ERROR; I.E., SMALLEST  PRONYX   1881
C            POSITIVE FLOATING POINT NUMBER WHICH CAN BE ADDED TO 1.0   PRONYX   1882
C            AND PRODUCE A RESULT LARGER THAN 1.0                       PRONYX   1883
C                                                                       PRONYX   1884
C  ARE    -  UNIT ERROR IN THE ADDITION OPERATION, ASSUMED EQUAL TO ETA.PRONYX   1885
C                                                                       PRONYX   1886
C  XMRE   -  UNIT ERROR IN THE MULTIPLY OPERATION, ASSUMED EQUAL TO ETA.PRONYX   1887
C                                                                       PRONYX   1888
      ETA = .5**(NBITSD-1)                                              PRONYX   1889
      ARE = ETA                                                         PRONYX   1890
      XMRE = ETA                                                        PRONYX   1891
      RLO = SMALNO/ETA                                                  PRONYX   1892
C                                                                       PRONYX   1893
C *** LOAD DOUBLE-PRECISION DATA VECTOR [NOT ACTUALLY NECESSARY]        PRONYX   1894
      DO 1002 I=1,NDEG                                                  PRONYX   1895
 1002 AD(I)=ADDP(I)                                                     PRONYX   1896
C                                                                       PRONYX   1897
C *** CLEAR WORKING ARRAYS                                              PRONYX   1898
      DO 1008 M=1,121                                                   PRONYX   1899
      K(M)=0.D0                                                         PRONYX   1900
      P(M)=0.D0                                                         PRONYX   1901
      QK(M)=0.D0                                                        PRONYX   1902
      QP(M)=0.D0                                                        PRONYX   1903
      SVK(M)=0.D0                                                       PRONYX   1904
      TEMP(M)=0.D0                                                      PRONYX   1905
 1008 CONTINUE                                                          PRONYX   1906
C                                                                       PRONYX   1907
C *** INITIALIZE CONSTANTS FOR SHIFT ROTATION                           PRONYX   1908
      XX = 1.D0/DSQRT(2.D0)                                             PRONYX   1909
      YY = -XX                                                          PRONYX   1910
      TWOPI=4.D0*DASIN(1.D0)                                            PRONYX   1911
      COSR=-DCOS(86.D0/TWOPI)                                           PRONYX   1912
      SINR=DSIN(86.D0/TWOPI)                                            PRONYX   1913
C                                                                       PRONYX   1914
      IFAIL = 0                                                         PRONYX   1915
C                                                                       PRONYX   1916
C *** MAKE A COPY OF THE COEFFICIENTS.  NOTE WE'RE REVERSING THE ORDER  PRONYX   1917
C     AND FORCING THE POLYNOMIAL MONIC                                  PRONYX   1918
      N = NDEG                                                          PRONYX   1919
      NN = N + 1                                                        PRONYX   1920
      P(1) = 1.D0                                                       PRONYX   1921
      DO 1010 I = 2, NN                                                 PRONYX   1922
      P(I) = AD(NN-I+1)                                                 PRONYX   1923
 1010 CONTINUE                                                          PRONYX   1924
C                                                                       PRONYX   1925
C *** REMOVE THE ZEROS AT THE ORIGIN IF ANY                             PRONYX   1926
 1015 IF (P(NN).NE.0.D0) GO TO 1040                                     PRONYX   1927
      J = NDEG - N + 1                                                  PRONYX   1928
      ZEROR(J) = 0.D0                                                   PRONYX   1929
      ZEROI(J) = 0.D0                                                   PRONYX   1930
      NN = NN - 1                                                       PRONYX   1931
      N = N - 1                                                         PRONYX   1932
      GO TO 1015                                                        PRONYX   1933
C START THE ALGORITHM FOR ONE ZERO                                      PRONYX   1934
 1040 IF(N.GT.2) GO TO 1070                                             PRONYX   1935
      IF(N.LT.1) GO TO 1870                                             PRONYX   1936
C CALCULATE THE FINAL ZERO OR PAIR OF ZEROS                             PRONYX   1937
      IF(N.EQ.2) GO TO 1050                                             PRONYX   1938
      ZEROR(NDEG) = -P(2)/P(1)                                          PRONYX   1939
      ZEROI(NDEG) = 0.D0                                                PRONYX   1940
      GO TO 1870                                                        PRONYX   1941
 1050 MPASS = 0                                                         PRONYX   1942
      AQ = P(1)                                                         PRONYX   1943
      BQ = P(2)                                                         PRONYX   1944
      CQ = P(3)                                                         PRONYX   1945
      GO TO 5000                                                        PRONYX   1946
 1060 ZEROR(NDEG-1) = SZR                                               PRONYX   1947
      ZEROI(NDEG-1) = SZI                                               PRONYX   1948
      ZEROR(NDEG) = LZR                                                 PRONYX   1949
      ZEROI(NDEG) = LZI                                                 PRONYX   1950
      GO TO 1870                                                        PRONYX   1951
C                                                                       PRONYX   1952
C FIND LARGEST AND SMALLEST MODULI OF COEFFICIENTS.                     PRONYX   1953
 1070 RMAX = 0.D0                                                       PRONYX   1954
      RMIN = RINFIN                                                     PRONYX   1955
      DO 1090 I = 1, NN                                                 PRONYX   1956
      VAL = DABS(P(I))                                                  PRONYX   1957
      IF (VAL.GT.RMAX) RMAX = VAL                                       PRONYX   1958
      IF (VAL.NE.0.D0.AND.VAL.LT.RMIN) RMIN = VAL                       PRONYX   1959
 1090 CONTINUE                                                          PRONYX   1960
C                                                                       PRONYX   1961
C  IF THERE EXIST VERY LARGE OR SMALL COEFFICIENTS, COMPUTE A           PRONYX   1962
C  MULTIPLYING SCALE FACTOR TO AVOID OVERFLOW OR UNDETECTED UNDERFLOW.  PRONYX   1963
C  THE FACTOR WILL BE A POWER OF THE BASE ABOVE.                        PRONYX   1964
C      IF (1.EQ.1) THEN                                                 PRONYX   1965
      IF (I.EQ.1) THEN                                                  PRONYX   1966
        GO TO 1150                                                      PRONYX   1967
      END IF                                                            PRONYX   1968
      SC = RLO/RMIN                                                     PRONYX   1969
      IF (SC.GT.1.D0) GO TO 1110                                        PRONYX   1970
      IF(RMAX.LT.10.D0) GO TO 1150                                      PRONYX   1971
      IF (SC.EQ.0.D0) THEN                                              PRONYX   1972
        SC = SMALNO                                                     PRONYX   1973
        GO TO 1120                                                      PRONYX   1974
      END IF                                                            PRONYX   1975
 1110 IF (RINFIN/SC.LT.RMAX) GO TO 1150                                 PRONYX   1976
 1120 L = DLOG(SC)/DLOG(BASE)+0.5D0                                     PRONYX   1977
      FACTOR = (BASE*1.D0)**L                                           PRONYX   1978
      IF (FACTOR.EQ.1.D0) GO TO 1150                                    PRONYX   1979
      DO 1140 I = 1, NN                                                 PRONYX   1980
 1140 P(I) = FACTOR*P(I)                                                PRONYX   1981
C                                                                       PRONYX   1982
C COMPUTE LOWER BOUND ON MODULI OF ZEROS.                               PRONYX   1983
 1150 DO 1160 I = 1, NN                                                 PRONYX   1984
 1160 PT(I) = DABS(P(I))                                                PRONYX   1985
      PT(NN) = -PT(NN)                                                  PRONYX   1986
C COMPUTE UPPER ESTIMATE OF BOUND                                       PRONYX   1987
      X = DEXP((DLOG(-PT(NN))-DLOG(PT(1)))/REAL(N))                     PRONYX   1988
      IF(PT(N).EQ.0.D0) GO TO 1180                                      PRONYX   1989
C IF NEWTON STEP AT THE ORIGIN IS BETTER, USE IT.                       PRONYX   1990
      XM = -PT(NN)/PT(N)                                                PRONYX   1991
      IF(XM.LT.X) X = XM                                                PRONYX   1992
C CHOP THE INTERVAL (0,X) UNTIL FF .LE. 0                               PRONYX   1993
 1180 XM = X*0.1D0                                                      PRONYX   1994
      FF = PT(1)                                                        PRONYX   1995
      DO 1190 I = 2,NN                                                  PRONYX   1996
 1190 FF = FF*XM + PT(I)                                                PRONYX   1997
      IF (FF.GT.0.D0) THEN                                              PRONYX   1998
        X = XM                                                          PRONYX   1999
        GO TO 1180                                                      PRONYX   2000
      END IF                                                            PRONYX   2001
C                                                                       PRONYX   2002
C *** DO NEWTON ITERATION UNTIL X CONVERGES TO THREE DECIMAL PLACES     PRONYX   2003
      DX = X                                                            PRONYX   2004
 1210 IF(DABS(DX/X).LE.0.0005D0) GO TO 1230                             PRONYX   2005
      FF = PT(1)                                                        PRONYX   2006
      DF = FF                                                           PRONYX   2007
      DO 1220 I = 2, N                                                  PRONYX   2008
      FF = FF*X + PT(I)                                                 PRONYX   2009
 1220 DF = DF*X + FF                                                    PRONYX   2010
      FF = FF*X + PT(NN)                                                PRONYX   2011
      DX = FF/DF                                                        PRONYX   2012
      X = X - DX                                                        PRONYX   2013
      GO TO 1210                                                        PRONYX   2014
 1230 BND = X                                                           PRONYX   2015
C                                                                       PRONYX   2016
C COMPUTE THE DERIVATIVE AS THE INTIAL K POLYNOMIAL                     PRONYX   2017
C AND DO 5 STEPS WITH NO SHIFT                                          PRONYX   2018
      NM1 = N - 1                                                       PRONYX   2019
      DO 1240 I = 2, N                                                  PRONYX   2020
      K(I) = REAL(NN-I)*P(I)/REAL(N)                                    PRONYX   2021
 1240 CONTINUE                                                          PRONYX   2022
      K(1) = P(1)                                                       PRONYX   2023
      AA = P(NN)                                                        PRONYX   2024
      BB = P(N)                                                         PRONYX   2025
      ZEROK=K(N).EQ.0.D0                                                PRONYX   2026
      DO 1330 JJ = 1, 5                                                 PRONYX   2027
      CC = K(N)                                                         PRONYX   2028
      IF(ZEROK) GO TO 1290                                              PRONYX   2029
C USE SCALED FORM OF RECURRENCE IF VALUE OF K AT 0 IS NON-ZERO.         PRONYX   2030
      T = -AA/CC                                                        PRONYX   2031
      DO 1280 I = 1, NM1                                                PRONYX   2032
      J = NN - I                                                        PRONYX   2033
 1280 K(J) = T*K(J-1) + P(J)                                            PRONYX   2034
      K(1) = P(1)                                                       PRONYX   2035
      ZEROK = DABS(K(N)).LE.DABS(BB)*ETA*10.D0                          PRONYX   2036
      GO TO 1330                                                        PRONYX   2037
C USE UNSCALED FORM OF RECURRENCE                                       PRONYX   2038
 1290 DO 1300 I = 1, NM1                                                PRONYX   2039
      J = NN - I                                                        PRONYX   2040
 1300 K(J) = K(J-1)                                                     PRONYX   2041
      K(1) = 0.D0                                                       PRONYX   2042
      ZEROK=K(N).EQ.0.D0                                                PRONYX   2043
 1330 CONTINUE                                                          PRONYX   2044
C                                                                       PRONYX   2045
C SAVE K FOR RESTARTS WITH NEW SHIFTS                                   PRONYX   2046
      DO 1340 I = 1, N                                                  PRONYX   2047
 1340 TEMP(I) = K(I)                                                    PRONYX   2048
C                                                                       PRONYX   2049
C LOOP TO SELECT THE QUADRATIC CORRESPONDING TO EACH NEW SHIFT.         PRONYX   2050
      KOUNT = 0                                                         PRONYX   2051
 1350 KOUNT = KOUNT + 1                                                 PRONYX   2052
C  QUADRATIC CORRESPONDS TO A DOUBLE SHIFT TO A NON-REAL POINT          PRONYX   2053
C  AND ITS COMPLEX CONJUGATE.  THE POINT HAS MODULUS BND WITH AMPLITUDE PRONYX   2054
C  ROTATED BY 94 DEGREES FROM THE PREVIOUS SHIFT.                       PRONYX   2055
      XXX = COSR*XX - SINR*YY                                           PRONYX   2056
      YY = SINR*XX + COSR*YY                                            PRONYX   2057
      XX = XXX                                                          PRONYX   2058
      SR = BND*XX                                                       PRONYX   2059
      SI = BND*YY                                                       PRONYX   2060
      U = -2.D0*SR                                                      PRONYX   2061
      V=BND*BND                                                         PRONYX   2062
C SECOND STAGE CALCULATION, FIXED QUADRATIC                             PRONYX   2063
C                                                                       PRONYX   2064
C  ***** SUBROUTINE FXSHFR(L2, NZ) ******                               PRONYX   2065
      L2 = 20 * KOUNT                                                   PRONYX   2066
C  COMPUTES UP TO L2 FIXED SHIFT K-POLYNOMIALS, TESTING FOR CONVERGENCE PRONYX   2067
C  IN THE LINEAR OR QUADRATIC CASE.  INITIATES ONE OF THE VARIABLE      PRONYX   2068
C  SHIFT ITERATIONS AND RETURNS WITH THE NUMBER OF ZEROS FOUND.         PRONYX   2069
C L2 - LIMIT OF FIXED SHIFT STEPS                                       PRONYX   2070
C NZ - NUMBER OF ZEROS FOUND                                            PRONYX   2071
      NZ = 0                                                            PRONYX   2072
      BETAV = 0.25D0                                                    PRONYX   2073
      BETAS = 0.25D0                                                    PRONYX   2074
      OSS = SR                                                          PRONYX   2075
      OVV = V                                                           PRONYX   2076
C EVALUATE POLYNOMIAL BY SYNTHETIC DIVISION                             PRONYX   2077
      B = P(1)                                                          PRONYX   2078
      QP(1) = B                                                         PRONYX   2079
      A = P(2) - U*B                                                    PRONYX   2080
      QP(2) = A                                                         PRONYX   2081
      DO 1360 I = 3, NN                                                 PRONYX   2082
      CD = P(I) - U*A - V*B                                             PRONYX   2083
      QP(I) = CD                                                        PRONYX   2084
      B = A                                                             PRONYX   2085
 1360 A = CD                                                            PRONYX   2086
      JPASS = 1                                                         PRONYX   2087
C *** TRANSFER TO CALCSE                                                PRONYX   2088
      GO TO 2000                                                        PRONYX   2089
 1370 JF = 0                                                            PRONYX   2090
 1380 JF = JF + 1                                                       PRONYX   2091
C CALCULATE NEXT K POLYNOMIAL AND ESTIMATE V, CALLING NEXTK,            PRONYX   2092
C CALCSC AND NEWEST                                                     PRONYX   2093
      KPASS = -1                                                        PRONYX   2094
      JPASS = 2                                                         PRONYX   2095
      LPASS = 0                                                         PRONYX   2096
C ***  TRANSFER TO NEXTK                                                PRONYX   2097
      GO TO 3000                                                        PRONYX   2098
 1390 VV = VI                                                           PRONYX   2099
C ESTIMATE S                                                            PRONYX   2100
      SS = 0.D0                                                         PRONYX   2101
      IF (K(N).NE.0.D0) SS = -P(NN) / K(N)                              PRONYX   2102
      TV = 1.D0                                                         PRONYX   2103
      TS = 1.D0                                                         PRONYX   2104
      IF(JF.EQ.1 .OR. NTYPE.EQ.3) GO TO 1800                            PRONYX   2105
C COMPUTE RELATIVE MEASURES OF CONVERGENCE OF S AND V SEQUENCES.        PRONYX   2106
      IF (VV.NE.0.D0) TV = DABS((VV-OVV)/VV)                            PRONYX   2107
      IF (SS.NE.0.D0) TS = DABS((SS-OSS)/SS)                            PRONYX   2108
C IF DECREASING, MULTIPLY TWO MOST RECENT CONVERGENCE MEASURES.         PRONYX   2109
      TVV = 1.D0                                                        PRONYX   2110
      IF(TV.LT.OTV) TVV = TV*OTV                                        PRONYX   2111
      TSS = 1.D0                                                        PRONYX   2112
      IF(TS.LT.OTS) TSS = TS*OTS                                        PRONYX   2113
C COMPARE WITH CONVERGENCE CRITERIA                                     PRONYX   2114
      VPASS = TVV.LT.BETAV                                              PRONYX   2115
      SPASS = TSS.LT.BETAS                                              PRONYX   2116
      IF(.NOT.(SPASS .OR. VPASS)) GO TO 1800                            PRONYX   2117
C  AT LEAST ONE SEQUENCE PASSED CONVERGENCE TEST.  STORE VARIABLES.     PRONYX   2118
      SVU = U                                                           PRONYX   2119
      SVV = V                                                           PRONYX   2120
      DO 1470 I = 1, N                                                  PRONYX   2121
      SVK(I) = K(I)                                                     PRONYX   2122
 1470 CONTINUE                                                          PRONYX   2123
      S = SS                                                            PRONYX   2124
C CHOOSE ITERATION ACCORDING TO THE FASTEST CONVERGING SEQUENCE.        PRONYX   2125
      VTRY = .FALSE.                                                    PRONYX   2126
      STRY = .FALSE.                                                    PRONYX   2127
      IF(SPASS .AND. ((.NOT.VPASS) .OR. TSS.LT.TVV)) GO TO 1640         PRONYX   2128
C                                                                       PRONYX   2129
C  ***** SUBROUTINE QUADIT(UI,VI,NZ) *****                              PRONYX   2130
C  VARIABLE-SHIFT K-POLYNOMIAL ITERATION FOR A QUADRATIC FACTOR.        PRONYX   2131
C  CONVERGES ONLY IF THE ZEROS ARE EQUIMODULAR OR NEARLY SO.            PRONYX   2132
C UI,VI - COEFFICIENTS OF STARTING QUADRATIC                            PRONYX   2133
C NZ - NUMBER OF ZERO FOUND                                             PRONYX   2134
 1480 NZ = 0                                                            PRONYX   2135
      TRIED = .FALSE.                                                   PRONYX   2136
      U = UI                                                            PRONYX   2137
      V = VI                                                            PRONYX   2138
      JQ = 0                                                            PRONYX   2139
C MAIN LOOP                                                             PRONYX   2140
 1490 MPASS = 1                                                         PRONYX   2141
      AQ = 1.D0                                                         PRONYX   2142
      BQ = U                                                            PRONYX   2143
      CQ = V                                                            PRONYX   2144
      GO TO 5000                                                        PRONYX   2145
C RETURN IF ROOTS OF THE QUADRATIC ARE REAL AND NOT                     PRONYX   2146
C CLOSE TO MULTIPLE OR NEARLY EQUAL AND OF OPPOSITE SIGN.               PRONYX   2147
 1500 IF(DABS(DABS(SZR)-DABS(LZR)).GT..01D0*DABS(LZR)) GO TO 1620       PRONYX   2148
C EVALUATE POLYNOMIAL BY QUADRATIC SYNTHETIC DIVISION                   PRONYX   2149
      B = P(1)                                                          PRONYX   2150
      QP(1) = B                                                         PRONYX   2151
      A = P(2) - U*B                                                    PRONYX   2152
      QP(2) = A                                                         PRONYX   2153
      DO 1510 I = 3, NN                                                 PRONYX   2154
      CD = P(I) - U*A - V*B                                             PRONYX   2155
      QP(I) = CD                                                        PRONYX   2156
      B = A                                                             PRONYX   2157
 1510 A = CD                                                            PRONYX   2158
      RMP = DABS(A-SZR*B) + DABS(SZI*B)                                 PRONYX   2159
C COMPUTE A RIGOROUS BOUND ON THE ROUNDING ERROR IN EVALUATING P.       PRONYX   2160
      ZM = DSQRT(DABS(V))                                               PRONYX   2161
      EE = 2.D0*DABS(QP(1))                                             PRONYX   2162
      TQ = -SZR*B                                                       PRONYX   2163
      DO 1520 I = 2, N                                                  PRONYX   2164
 1520 EE = EE*ZM + DABS(QP(I))                                          PRONYX   2165
      EE = EE*ZM + DABS(A+TQ)                                           PRONYX   2166
      EE = (5.D0*XMRE+4.D0*ARE)*EE - (5.D0*XMRE+2.D0*ARE)*(DABS(A+TQ)+  PRONYX   2167
     * DABS(B)*ZM) + 2.D0*ARE*DABS(TQ)                                  PRONYX   2168
C ITERATION HAS CONVERGED SUFFICIENTLY IF THE                           PRONYX   2169
C POLYNOMIAL VALUE IS LESS THAN 20 TIMES THIS BOUND                     PRONYX   2170
      IF(RMP.GT.20.D0*EE) GO TO 1530                                    PRONYX   2171
      NZ = 2                                                            PRONYX   2172
      GO TO 1620                                                        PRONYX   2173
 1530 JQ = JQ + 1                                                       PRONYX   2174
C STOP ITERATION AFTER 20 STEPS                                         PRONYX   2175
      IF(JQ.GT.20) GO TO 1620                                           PRONYX   2176
      IF(JQ.LT.2) GO TO 1580                                            PRONYX   2177
      IF (RELSTP.GT.0.01D0.OR.RMP.LT.OMP.OR.TRIED) GO TO 1580           PRONYX   2178
C A CLUSTER APPEARS TO BE STALLING THE CONVERGENCE.  FIVE FIXED         PRONYX   2179
C SHIFT STEPS ARE TAKEN WITH A (U,V) CLOSE TO THE CLUSTER.              PRONYX   2180
      IF(RELSTP.LT.ETA) RELSTP = ETA                                    PRONYX   2181
      RELSTP = DSQRT(RELSTP)                                            PRONYX   2182
      U = U - U*RELSTP                                                  PRONYX   2183
      V = V + V*RELSTP                                                  PRONYX   2184
      B = P(1)                                                          PRONYX   2185
      QP(1) = B                                                         PRONYX   2186
      A = P(2) - U*B                                                    PRONYX   2187
      QP(2) = A                                                         PRONYX   2188
      DO 1540 I = 3, NN                                                 PRONYX   2189
      CD = P(I) - U*A - V*B                                             PRONYX   2190
      QP(I) = CD                                                        PRONYX   2191
      B = A                                                             PRONYX   2192
 1540 A = CD                                                            PRONYX   2193
C  LOOP TO CALL CALCSC AND NEXTK FIVE TIMES.                            PRONYX   2194
      IQ = 0                                                            PRONYX   2195
      JPASS = 3                                                         PRONYX   2196
      KPASS = 0                                                         PRONYX   2197
 1550 IQ = IQ + 1                                                       PRONYX   2198
      GO TO 2000                                                        PRONYX   2199
 1560 IF (IQ.LT.5) GO TO 1550                                           PRONYX   2200
      TRIED = .TRUE.                                                    PRONYX   2201
      JQ = 0                                                            PRONYX   2202
 1580 OMP = RMP                                                         PRONYX   2203
C  CALCULATE NEXT K POLYNOMIAL AND NEW U AND V                          PRONYX   2204
C  CALL CALCSC, NEXTK, CALCSC, AND NEWEST IN SEQUENCE.                  PRONYX   2205
      JPASS = 4                                                         PRONYX   2206
      KPASS = 1                                                         PRONYX   2207
      GO TO 2000                                                        PRONYX   2208
 1590 JPASS = 5                                                         PRONYX   2209
      LPASS = 1                                                         PRONYX   2210
      GO TO 2000                                                        PRONYX   2211
C IF VI IS ZERO THE ITERATION IS NOT CONVERGING                         PRONYX   2212
 1600 IF (VI.NE.0.D0) THEN                                              PRONYX   2213
        RELSTP = DABS((VI-V)/VI)                                        PRONYX   2214
        U = UI                                                          PRONYX   2215
        V = VI                                                          PRONYX   2216
        GO TO 1490                                                      PRONYX   2217
      END IF                                                            PRONYX   2218
 1620 IF(NZ.GT.0) GO TO 1820                                            PRONYX   2219
C QUADRATIC ITERATION HAS FAILED. FLAG THAT IT HAS                      PRONYX   2220
C BEEN TRIED AND DECREASE THE CONVERGENCE CRITERION.                    PRONYX   2221
      VTRY = .TRUE.                                                     PRONYX   2222
      BETAV = BETAV*0.25D0                                              PRONYX   2223
C  ***** END OF SUBROUTINE QUADIT *****                                 PRONYX   2224
C                                                                       PRONYX   2225
C TRY LINEAR ITERATION IF IT HAS NOT BEEN TRIED AND                     PRONYX   2226
C THE S SEQUENCE IS CONVERGING                                          PRONYX   2227
      IF(STRY .OR. (.NOT.SPASS)) GO TO 1770                             PRONYX   2228
      DO 1630 I = 1, N                                                  PRONYX   2229
 1630 K(I) = SVK(I)                                                     PRONYX   2230
C                                                                       PRONYX   2231
C ***** SUBROUTINE REALIT(S, NZ, IFLAG) *****                           PRONYX   2232
C VARIABLE-SHIFT H POLYNOMIAL ITERATION FOR A REAL ZERO.                PRONYX   2233
C S     - STARTING ITERATE                                              PRONYX   2234
C NZ    - NUMBER OF ZERO FOUND                                          PRONYX   2235
C IFLAG - FLAG TO INDICATE A PAIR OF ZEROS NEAR REAL AXIS               PRONYX   2236
 1640 NM1 = N - 1                                                       PRONYX   2237
      NZ = 0                                                            PRONYX   2238
      SSS = S                                                           PRONYX   2239
      IFLAG = 0                                                         PRONYX   2240
      JR = 0                                                            PRONYX   2241
C MAIN LOOP                                                             PRONYX   2242
 1650 PV = P(1)                                                         PRONYX   2243
C EVALUATE P AT SSS                                                     PRONYX   2244
      QP(1) = PV                                                        PRONYX   2245
      DO 1660 I = 2, NN                                                 PRONYX   2246
      PV = PV*SSS + P(I)                                                PRONYX   2247
 1660 QP(I) = PV                                                        PRONYX   2248
      RMP = DABS(PV)                                                    PRONYX   2249
C COMPUTE A RIGOROUS BOUND ON THE ERROR IN EVALUATING P                 PRONYX   2250
      RMS = DABS(SSS)                                                   PRONYX   2251
      EE = (XMRE/(ARE+XMRE))*DABS(QP(1))                                PRONYX   2252
      DO 1670 I = 2, NN                                                 PRONYX   2253
 1670 EE = EE*RMS + DABS(QP(I))                                         PRONYX   2254
C ITERATION HAS CONVERGED SUFFICIENTLY IF THE                           PRONYX   2255
C POLYNOMIAL VALUE IS LESS THAN 20 TIMES THIS BOUND                     PRONYX   2256
      IF(RMP.GT.20.D0*((ARE+XMRE)*EE-XMRE*RMP)) GO TO 1680              PRONYX   2257
      NZ = 1                                                            PRONYX   2258
      SZR = SSS                                                         PRONYX   2259
      SZI = 0.D0                                                        PRONYX   2260
      GO TO 1760                                                        PRONYX   2261
 1680 JR = JR + 1                                                       PRONYX   2262
C STOP ITERATION AFTER 10 STEPS                                         PRONYX   2263
      IF(JR.GT.10) GO TO 1760                                           PRONYX   2264
      IF(JR.LT.2) GO TO 1690                                            PRONYX   2265
      IF (DABS(T).GT.0.001D0*DABS(SSS-T).OR.RMP.LE.OMP) GO TO 1690      PRONYX   2266
C A CLUSTER OF ZEROS NEAR THE REAL AXIS HAS BEEN ENCOUNTERED.           PRONYX   2267
C RETURN WITH IFLAG SET TO INITIATE A QUADRATIC ITERATION.              PRONYX   2268
      IFLAG = 1                                                         PRONYX   2269
      S = SSS                                                           PRONYX   2270
      GO TO 1760                                                        PRONYX   2271
C RETURN IF THE POLYNOMIAL VALUE HAS INCREASED SIGNIFICANTLY            PRONYX   2272
 1690 OMP = RMP                                                         PRONYX   2273
C COMPUTE T, THE NEXT POLYNOMIAL, AND THE NEW ITERATE                   PRONYX   2274
      KV = K(1)                                                         PRONYX   2275
      QK(1) = KV                                                        PRONYX   2276
      DO 1700 I = 2, N                                                  PRONYX   2277
      KV = KV*SSS + K(I)                                                PRONYX   2278
 1700 QK(I) = KV                                                        PRONYX   2279
      IF(DABS(KV).LE.DABS(K(N))*10.D0*ETA) GO TO 1720                   PRONYX   2280
C USE THE SCALED FORM OF THE RECURRENCE IF THE VALUE OF K AT            PRONYX   2281
C SSS IS NONZERO                                                        PRONYX   2282
      T = -PV/KV                                                        PRONYX   2283
      K(1) = QP(1)                                                      PRONYX   2284
      DO 1710 I = 2, N                                                  PRONYX   2285
 1710 K(I) = T*QK(I-1) + QP(I)                                          PRONYX   2286
      GO TO 1740                                                        PRONYX   2287
C USE UNSCALED FORM                                                     PRONYX   2288
 1720 K(1) = 0.D0                                                       PRONYX   2289
      DO 1730 I = 2, N                                                  PRONYX   2290
 1730 K(I) = QK(I-1)                                                    PRONYX   2291
 1740 KV = K(1)                                                         PRONYX   2292
      DO 1750 I = 2, N                                                  PRONYX   2293
 1750 KV = KV*SSS + K(I)                                                PRONYX   2294
      T = 0.D0                                                          PRONYX   2295
      IF(DABS(KV).GT.DABS(K(N))*10.D0*ETA) T = -PV/KV                   PRONYX   2296
      SSS = SSS + T                                                     PRONYX   2297
      GO TO 1650                                                        PRONYX   2298
C ****** END OF SUBROUTINE REALIT  *****                                PRONYX   2299
C                                                                       PRONYX   2300
 1760 IF(NZ.GT.0) GO TO 1820                                            PRONYX   2301
C LINEAR ITERATION HAS FAILED. FLAG THAT IT HAS BEEN                    PRONYX   2302
C TRIED AND DECREASE THE CONVERGENCE CRITERION                          PRONYX   2303
      STRY = .TRUE.                                                     PRONYX   2304
      BETAS = BETAS*0.25D0                                              PRONYX   2305
      IF(IFLAG.EQ.0) GO TO 1770                                         PRONYX   2306
C IF LINEAR ITERATION SIGNALS AN ALMOST DOUBLE REAL                     PRONYX   2307
C ZERO ATTEMPT QUADRATIC INTERATION                                     PRONYX   2308
      UI = -(S+S)                                                       PRONYX   2309
      VI = S*S                                                          PRONYX   2310
      GO TO 1480                                                        PRONYX   2311
C RESTORE VARIABLES                                                     PRONYX   2312
 1770 U = SVU                                                           PRONYX   2313
      V = SVV                                                           PRONYX   2314
      DO 1780 I = 1, N                                                  PRONYX   2315
 1780 K(I) = SVK(I)                                                     PRONYX   2316
C TRY QUADRATIC ITERATION IF IT HAS NOT BEEN TRIED                      PRONYX   2317
C AND THE V SEQUENCE IS CONVERGING                                      PRONYX   2318
      IF(VPASS .AND. (.NOT.VTRY)) GO TO 1480                            PRONYX   2319
C RECOMPUTE QP AND SCALAR VALUES TO CONTINUE THE SECOND STAGE.          PRONYX   2320
      B = P(1)                                                          PRONYX   2321
      QP(1) = B                                                         PRONYX   2322
      A = P(2) - U*B                                                    PRONYX   2323
      QP(2) = A                                                         PRONYX   2324
      DO 1790 I = 3, NN                                                 PRONYX   2325
      CD = P(I) - U*A - V*B                                             PRONYX   2326
      QP(I) = CD                                                        PRONYX   2327
      B = A                                                             PRONYX   2328
 1790 A = CD                                                            PRONYX   2329
      JPASS = 6                                                         PRONYX   2330
      GO TO 2000                                                        PRONYX   2331
 1800 OVV = VV                                                          PRONYX   2332
      OSS = SS                                                          PRONYX   2333
      OTV = TV                                                          PRONYX   2334
      OTS = TS                                                          PRONYX   2335
      IF (JF.LT.L2) GO TO 1380                                          PRONYX   2336
C ***** END OF SUBROUTINE FXSHFR *****                                  PRONYX   2337
C                                                                       PRONYX   2338
      IF (NZ.EQ.0) GO TO 1840                                           PRONYX   2339
C                                                                       PRONYX   2340
C  THE SECOND STAGE JUMPS DIRECTLY TO ONE OF THE THIRD STAGE ITERATIONS PRONYX   2341
C  AND RETURNS HERE IF SUCCESSFUL.  DEFLATE THE POLYNOMIAL, STORE THE   PRONYX   2342
C  ZERO OR ZEROS, AND RETURN TO THE MAIN ALGORITHM.                     PRONYX   2343
 1820 J = NDEG - N + 1                                                  PRONYX   2344
      ZEROR(J) = SZR                                                    PRONYX   2345
      ZEROI(J) = SZI                                                    PRONYX   2346
      NN = NN - NZ                                                      PRONYX   2347
      N = NN - 1                                                        PRONYX   2348
      DO 1830 I = 1, NN                                                 PRONYX   2349
 1830 P(I) = QP(I)                                                      PRONYX   2350
      IF(NZ.EQ.1) GO TO 1040                                            PRONYX   2351
      ZEROR(J+1) = LZR                                                  PRONYX   2352
      ZEROI(J+1) = LZI                                                  PRONYX   2353
      GO TO 1040                                                        PRONYX   2354
C                                                                       PRONYX   2355
C IF THE ITERATION IS UNSUCCESSFUL ANOTHER QUADRATIC                    PRONYX   2356
C IS CHOSEN AFTER RESTORING K                                           PRONYX   2357
 1840 DO 1850 I = 1, N                                                  PRONYX   2358
 1850 K(I) = TEMP(I)                                                    PRONYX   2359
      IF (KOUNT.LT.20) GO TO 1350                                       PRONYX   2360
C RETURN WITH FAILURE IF NO CONVERGENCE WITH 20 SHIFTS                  PRONYX   2361
      IFAIL = 1                                                         PRONYX   2362
      NDEG = NDEG - N                                                   PRONYX   2363
C                                                                       PRONYX   2364
      GO TO 8000                                                        PRONYX   2365
C                                                                       PRONYX   2366
 1870 IF (NDEG.LE.2) GO TO 8000                                         PRONYX   2367
C                                                                       PRONYX   2368
C *** SECTION TO RECONSTRUCT COEFFICIENTS AND CHECK ACCURACY            PRONYX   2369
      C1 = ZEROR(1)                                                     PRONYX   2370
      C2 = ZEROI(1)                                                     PRONYX   2371
      B1 = ZEROR(2)                                                     PRONYX   2372
      B2 = ZEROI(2)                                                     PRONYX   2373
      K(1) = C1 * B1 - C2 * B2                                          PRONYX   2374
      P(1) = C1 * B2 + C2 * B1                                          PRONYX   2375
      K(2) = -C1 - B1                                                   PRONYX   2376
      P(2) = -C2 - B2                                                   PRONYX   2377
      DO 1880 I = 3, NDEG                                               PRONYX   2378
      B1 = ZEROR(I)                                                     PRONYX   2379
      B2 = ZEROI(I)                                                     PRONYX   2380
      K(I) = 1.D0                                                       PRONYX   2381
      P(I) = 0.D0                                                       PRONYX   2382
      C1=K(1)                                                           PRONYX   2383
      C2=P(1)                                                           PRONYX   2384
      K(1) = -C1 * B1 + C2 * B2                                         PRONYX   2385
      P(1) = -C2 * B1 - C1 * B2                                         PRONYX   2386
      DO 1880 J=2, I                                                    PRONYX   2387
      C3 = K(J)                                                         PRONYX   2388
      C4 = P(J)                                                         PRONYX   2389
      K(J) = C1 - B1 * C3 + B2 * C4                                     PRONYX   2390
      P(J) = C2 - B1 * C4 - B2 * C3                                     PRONYX   2391
      C1=C3                                                             PRONYX   2392
 1880 C2=C4                                                             PRONYX   2393
C                                                                       PRONYX   2394
C *** COMPUTE EUCLIDIAN NORMS OF ERROR VECTOR AND ORIGINAL COEFFICIENTS PRONYX   2395
      ANORM = 0.D0                                                      PRONYX   2396
      ENORM = 0.D0                                                      PRONYX   2397
      DO 1890 I=1,NDEG                                                  PRONYX   2398
      ENORM = ENORM + (AD(I) - K(I))**2                                 PRONYX   2399
 1890 ANORM = ANORM + AD(I)**2                                          PRONYX   2400
      ENORM = DSQRT(ENORM/ANORM)                                        PRONYX   2401
      IF (ENORM.GT.ERRTOL) THEN                                         PRONYX   2402
        IFAIL=-1                                                        PRONYX   2403
        WRITE (*,9020) ENORM,ERRTOL                                     PRONYX   2404
        WRITE (6,9020) ENORM,ERRTOL                                     PRONYX   2405
        WRITE (*,9040) IFAIL                                            PRONYX   2406
        WRITE (6,9040) IFAIL                                            PRONYX   2407
      END IF                                                            PRONYX   2408
C                                                                       PRONYX   2409
      GO TO 8000                                                        PRONYX   2410
C                                                                       PRONYX   2411
C ***** SUBROUTINE CALCSC(NTYPE) *****                                  PRONYX   2412
C  THIS ROUTINE CALCULATES SCALAR QUANTITIES USED TO COMPUTE THE NEXT   PRONYX   2413
C  K-POLYNOMIAL AND NEW ESTIMATES OF THE QUADRATIC COEFFICIENTS.        PRONYX   2414
C  NTYPE- INTEGER VARIABLE SET HERE INDICATING HOW THE CALCULATIONS     PRONYX   2415
C      ARE NORMALIZED TO AVOID OVERFLOW.                                PRONYX   2416
C                                                                       PRONYX   2417
C SYNTHETIC DIVISION OF K BY THE QUADRATIC (1,U,V)                      PRONYX   2418
 2000 D = K(1)                                                          PRONYX   2419
      QK(1) = D                                                         PRONYX   2420
      C = K(2) - U*D                                                    PRONYX   2421
      QK(2) = C                                                         PRONYX   2422
      DO 2010 I = 3, N                                                  PRONYX   2423
      CD = K(I) - U*C - V*D                                             PRONYX   2424
      QK(I) = CD                                                        PRONYX   2425
      D = C                                                             PRONYX   2426
 2010 C = CD                                                            PRONYX   2427
      IF(DABS(C).GT.DABS(K(N))*100.D0*ETA) GO TO 2020                   PRONYX   2428
      IF(DABS(D).GT.DABS(K(N-1))*100.D0*ETA) GO TO 2020                 PRONYX   2429
      NTYPE = 3                                                         PRONYX   2430
C NTYPE=3 INDICATES THE QUADRATIC IS ALMOST A FACTOR OF K.              PRONYX   2431
      GO TO (1370, 4000, 3000, 3000, 4000, 1800), JPASS                 PRONYX   2432
 2020 IF(DABS(D).LT.DABS(C)) GO TO 2030                                 PRONYX   2433
      NTYPE = 2                                                         PRONYX   2434
C NTYPE=2 INDICATES THAT ALL FORMULAS ARE DIVIDED BY D                  PRONYX   2435
      E = A/D                                                           PRONYX   2436
      F = C/D                                                           PRONYX   2437
      G = U*B                                                           PRONYX   2438
      H = V*B                                                           PRONYX   2439
      A3 = (A+G)*E + H*(B/D)                                            PRONYX   2440
      A1 = B*F - A                                                      PRONYX   2441
      A7 = (F+U)*A + H                                                  PRONYX   2442
      GO TO (1370, 4000, 3000, 3000, 4000, 1800), JPASS                 PRONYX   2443
 2030 NTYPE = 1                                                         PRONYX   2444
C NTYPE=1 INDICATES THAT ALL FORMULAS ARE DIVIDED BY C                  PRONYX   2445
      E = A/C                                                           PRONYX   2446
      F = D/C                                                           PRONYX   2447
      G = U*E                                                           PRONYX   2448
      H = V*B                                                           PRONYX   2449
      A3 = A*E + (H/C+G)*B                                              PRONYX   2450
      A1 = B - A*(D/C)                                                  PRONYX   2451
      A7 = A + G*D + H*F                                                PRONYX   2452
      GO TO (1370, 4000, 3000, 3000, 4000, 1800), JPASS                 PRONYX   2453
C  ***** END OF SUBROUTINE CALCSC *****                                 PRONYX   2454
C                                                                       PRONYX   2455
C  ***** SUBROUTINE NEXTK(NTYPE) ******                                 PRONYX   2456
C COMPUTES THE NEXT K POLYNOMIALS USING SCALARS COMPUTED IN CALCSC.     PRONYX   2457
 3000 IF(NTYPE.EQ.3) GO TO 3040                                         PRONYX   2458
      TNEX = A                                                          PRONYX   2459
      IF(NTYPE.EQ.1) TNEX = B                                           PRONYX   2460
      IF(DABS(A1).GT.DABS(TNEX)*ETA*10.D0) GO TO 3020                   PRONYX   2461
C IF A1 IS NEARLY ZERO THEN USE A SPECIAL FORM OF THE RECURRENCE        PRONYX   2462
      K(1) = 0.D0                                                       PRONYX   2463
      K(2) = -A7*QP(1)                                                  PRONYX   2464
      DO 3010 I = 3, N                                                  PRONYX   2465
 3010 K(I) = A3*QK(I-2) - A7*QP(I-1)                                    PRONYX   2466
      IF(KPASS) 2000, 1560, 1590                                        PRONYX   2467
C USE SCALED FORM OF THE RECURRENCE                                     PRONYX   2468
 3020 A7 = A7/A1                                                        PRONYX   2469
      A3 = A3/A1                                                        PRONYX   2470
      K(1) = QP(1)                                                      PRONYX   2471
      K(2) = QP(2) - A7*QP(1)                                           PRONYX   2472
      DO 3030 I = 3, N                                                  PRONYX   2473
 3030 K(I) = A3*QK(I-2) - A7*QP(I-1) + QP(I)                            PRONYX   2474
      IF(KPASS) 2000, 1560, 1590                                        PRONYX   2475
C USE UNSCALED FORM OF THE RECURRENCE IF NTYPE IS 3                     PRONYX   2476
 3040 K(1) = 0.D0                                                       PRONYX   2477
      K(2) = 0.D0                                                       PRONYX   2478
      DO 3050 I = 3, N                                                  PRONYX   2479
 3050 K(I) = QK(I-2)                                                    PRONYX   2480
C  ***** END OF SUBROUTINE NEXTK *****                                  PRONYX   2481
C                                                                       PRONYX   2482
      IF (KPASS.LT.0) GO TO 2000                                        PRONYX   2483
      IF (KPASS.EQ.0) GO TO 1560                                        PRONYX   2484
      GO TO 1590                                                        PRONYX   2485
C                                                                       PRONYX   2486
C ***** SUBROUTINE NEWEST(NTYPE, UI, VI) *****                          PRONYX   2487
C  COMPUTES NEW ESTIMATES OF THE QUADRATIC COEFFICIENTS USING SCALARS   PRONYX   2488
C  COMPUTED IN CALCSC.  USE THE FOMULAS APPROPRIATE TO SETTING OF NTYPE.PRONYX   2489
 4000 IF(NTYPE.EQ.3) GO TO 4040                                         PRONYX   2490
      IF(NTYPE.EQ.2) GO TO 4010                                         PRONYX   2491
      A4 = A + U*B + H*F                                                PRONYX   2492
      A5 = C + (U+V*F)*D                                                PRONYX   2493
      GO TO 4020                                                        PRONYX   2494
 4010 A4 = (A+G)*F + H                                                  PRONYX   2495
      A5 = (F+U)*C + V*D                                                PRONYX   2496
C EVALUATE NEW QUADRATIC COEFFICIENTS.                                  PRONYX   2497
 4020 B1 = -K(N)/P(NN)                                                  PRONYX   2498
      B2 = -(K(N-1)+B1*P(N))/P(NN)                                      PRONYX   2499
      C1 = V*B2*A1                                                      PRONYX   2500
      C2 = B1*A7                                                        PRONYX   2501
      C3 = B1*B1*A3                                                     PRONYX   2502
      C4 = C1 - C2 - C3                                                 PRONYX   2503
      TNEW = A5 + B1*A4 - C4                                            PRONYX   2504
      IF(TNEW) 4030, 4040, 4030                                         PRONYX   2505
 4030 UI = U - (U*(C3+C2)+V*(B1*A1+B2*A7))/TNEW                         PRONYX   2506
      VI = V*(1.D0+C4/TNEW)                                             PRONYX   2507
      IF(LPASS) 1390, 1390, 1600                                        PRONYX   2508
C IF NTYPE=3 THE QUADRATIC IS ZEROED                                    PRONYX   2509
 4040 UI = 0.D0                                                         PRONYX   2510
      VI = 0.D0                                                         PRONYX   2511
      IF (LPASS.LE.0) GO TO 1390                                        PRONYX   2512
      GO TO 1600                                                        PRONYX   2513
C  ***** END OF SUBROUTINE NEWEST *****                                 PRONYX   2514
C                                                                       PRONYX   2515
C  ***** SUBROUTINE QUAD(AQ, BQ, CQ, SZR, SZI, LZR, LZI)  ******        PRONYX   2516
C  CALCULATE THE ZEROS OF THE QUADRATIC     AQ*Z**2 + BQ*Z + CQ         PRONYX   2517
C  WITH MODIFICATIONS TO AVOID OVERFLOW.  ALSO USED TO FIND THE LARGER  PRONYX   2518
C  ZERO IF THE ZEROS ARE REAL AND BOTH ARE COMPLEX.  THE SMALLER REAL   PRONYX   2519
C  ZERO IS FOUND DIRECTLY FROM THE PRODUCT OF ZEROS CQ/AQ.              PRONYX   2520
 5000 IF (AQ.NE.0.D0) GO TO 5050                                        PRONYX   2521
      SZR = 0.D0                                                        PRONYX   2522
      IF (BQ.EQ.0.D0) GO TO 5030                                        PRONYX   2523
      SZR = - CQ/BQ                                                     PRONYX   2524
 5030 LZR = 0.D0                                                        PRONYX   2525
 5040 SZI = 0.D0                                                        PRONYX   2526
      LZI = 0.D0                                                        PRONYX   2527
      IF (MPASS.LE.0) GO TO 1060                                        PRONYX   2528
      GO TO 1500                                                        PRONYX   2529
 5050 IF (CQ.NE.0.D0) GO TO 5070                                        PRONYX   2530
      SZR = 0.D0                                                        PRONYX   2531
      LZR = -BQ/AQ                                                      PRONYX   2532
      GO TO 5040                                                        PRONYX   2533
C COMPUTE DISCRIMINANT AVOIDING OVERFLOW                                PRONYX   2534
 5070 QB = BQ/2.D0                                                      PRONYX   2535
      IF (DABS(QB).GE.DABS(CQ)) THEN                                    PRONYX   2536
        EQ = 1.D0 - (AQ/QB)*(CQ/QB)                                     PRONYX   2537
        DISC = DSQRT(DABS(EQ))*DABS(QB)                                 PRONYX   2538
        GO TO 5090                                                      PRONYX   2539
      END IF                                                            PRONYX   2540
      EQ = AQ                                                           PRONYX   2541
      IF(CQ.LT.0.D0) EQ = -AQ                                           PRONYX   2542
      EQ = QB*(QB/DABS(CQ)) - EQ                                        PRONYX   2543
      DISC = DSQRT(DABS(EQ))*DSQRT(DABS(CQ))                            PRONYX   2544
 5090 IF (EQ.LT.0.D0) GO TO 5110                                        PRONYX   2545
C REAL ZEROS                                                            PRONYX   2546
      IF(QB.GE.0.D0) DISC = -DISC                                       PRONYX   2547
      LZR = (-QB+DISC)/AQ                                               PRONYX   2548
      SZR = 0.D0                                                        PRONYX   2549
      IF (LZR.NE.0.D0) SZR=(CQ/LZR)/AQ                                  PRONYX   2550
      GO TO 5040                                                        PRONYX   2551
C COMPLEX CONJUGATE ZEROS                                               PRONYX   2552
 5110 SZR = -QB/AQ                                                      PRONYX   2553
      LZR = SZR                                                         PRONYX   2554
      SZI = DABS(DISC/AQ)                                               PRONYX   2555
      LZI = -SZI                                                        PRONYX   2556
      IF (MPASS.LE.0) GO TO 1060                                        PRONYX   2557
      GO TO 1500                                                        PRONYX   2558
C  ***** END OF SUBROUTINE QUAD *****                                   PRONYX   2559
C                                                                       PRONYX   2560
C *** LOAD QUAD-PRECISION RESULTS INTO DOUBLE-PRECISION RETURN          PRONYX   2561
 8000 TRACK=.FALSE.                                                     PRONYX   2562
      IF (TRACK) WRITE (*,9920)                                         PRONYX   2563
      IF (TRACK) WRITE (6,9920)                                         PRONYX   2564
      DO 8020 I=1,NDEG                                                  PRONYX   2565
      ZERORQ(I)=ZEROR(I)                                                PRONYX   2566
      ZEROIQ(I)=ZEROI(I)                                                PRONYX   2567
      IF (TRACK) WRITE (*,9922) ZEROR(I),ZERORQ(I),ZEROI(I),ZEROIQ(I)   PRONYX   2568
      IF (TRACK) WRITE (6,9922) ZEROR(I),ZERORQ(I),ZEROI(I),ZEROIQ(I)   PRONYX   2569
 8020 CONTINUE                                                          PRONYX   2570
C                                                                       PRONYX   2571
 9020 FORMAT (2X,'##IN DPOLY: ENORM =',E15.7,' ERRTOL =',E15.7)         PRONYX   2572
 9040 FORMAT (4X,'SETTING IFAIL =',I3)                                  PRONYX   2573
 9920 FORMAT (2X,'##IN DPOLY: DIAGNOSTICS')                             PRONYX   2574
 9922 FORMAT (5X,4D20.10)                                               PRONYX   2575
C                                                                       PRONYX   2576
*ENDIF                                                                  PRONYX   2577
      RETURN                                                            PRONYX   2578
      END                                                               PRONYX   2579
C                                                                       PRONYX   2580
        SUBROUTINE DSETTP (X,M,C,N,P,N1,FB,FLAG,LUOUT,BATCH,MRET)       PRONYX   2581
CD           This routine fills the Toeplitz M matrix according to the  PRONYX   2582
CD           linear prediction methods (Forword, backward or both),     PRONYX   2583
CD           (Corrn, covariance presse or post windowed).               PRONYX   2584
CD                                                                      PRONYX   2585
C Date : Nov 1987                                                       PRONYX   2586
CD                                                                      PRONYX   2587
CD For example, for the covariance forward-backward method, the system lPRONYX   2588
CD                                                                      PRONYX   2589
CD       |x(p-1) .......  x(0)|  | a(1) |     | x(p) |                  PRONYX   2590
CD       | x(p)  .......  x(1)|  | a(2) |     |x(p+1)|                  PRONYX   2591
CD       |  .    .......    . |  |  .   |     |  .   |                  PRONYX   2592
CD       |  .    .......    . | .|  .   | = - |  .   |                  PRONYX   2593
CD       |  .    .......    . |  |  .   |     |  .   |                  PRONYX   2594
CD       |  .    .......    . |  |a(p-1)|     |x(n-2)|                  PRONYX   2595
CD       |x(n-2) .....x(n-p-1)|  | a(p) |     |x(n-1)|                  PRONYX   2596
CD       ----------------------               --------                  PRONYX   2597
CD       | x(1)  .......  x(p)|               | x(0) |                  PRONYX   2598
CD       |  .    .......    . |               |  .   |                  PRONYX   2599
CD       |  .    .......    . |               |  .   |                  PRONYX   2600
CD       |  .    .......    . |               |  .   |                  PRONYX   2601
CD       |  .    .......    . |               |  .   |                  PRONYX   2602
CD       |x(n-p) .......x(n-1)|              |x(n-p-1)|                 PRONYX   2603
CD                                                                      PRONYX   2604
CD Calling Parameters :                                                 PRONYX   2605
CD                                                                      PRONYX   2606
CD   Input :     X(0:N-1): Data sequence.                               PRONYX   2607
CD               N       : Size of the vector                           PRONYX   2608
CD               P       : AR order                                     PRONYX   2609
CD               FLAG    : Method option                                PRONYX   2610
CD                               = 1   =>  Correlation method           PRONYX   2611
CD                               = 2   =>  Pre-windowed method          PRONYX   2612
CD                               = 3   =>  Covariance method            PRONYX   2613
CD                               = 4   =>  Post-windowed method         PRONYX   2614
CD               FB      : Forward-Backward option                      PRONYX   2615
CD                               = 1   =>  Forward                      PRONYX   2616
CD                               = 2   =>     "   -Backward             PRONYX   2617
CD                               = 3   =>  Backward                     PRONYX   2618
CD   Ouput :     M       : (Forward)-(Backward) matrix                  PRONYX   2619
CD               C       : Left side of the equation                    PRONYX   2620
CD               N1      : Number of equations in the system            PRONYX   2621
CD               LUOUT   : LOGICAL UNIT FOR DISPLAY                     PRONYX   2622
CD               BATCH   : FLAG FOR BATCH EXECUTION                     PRONYX   2623
CD               MRET    : RETURN MODE                                  PRONYX   2624
CD                                                                      PRONYX   2625
C *** LAST JFH MODIFICATION 07/05/89                                    PRONYX   2626
        INTEGER NSIZE,NTAIL                                             PRONYX   2627
        PARAMETER (NSIZE=128)                                           PRONYX   2628
        PARAMETER (NTAIL=1024)                                          PRONYX   2629
        DOUBLE PRECISION X(0:*),C(1:*),M(1:NTAIL,1:NSIZE)               PRONYX   2630
        INTEGER P,N,N1,N2,I,J,FB,FLAG,LUOUT,MRET,K,L,INDEX              PRONYX   2631
        LOGICAL BATCH                                                   PRONYX   2632
                                                                        PRONYX   2633
        IF (N.EQ.P) THEN                                                PRONYX   2634
          P = N - 1                                                     PRONYX   2635
          IF (BATCH) WRITE (*,9920) P                                   PRONYX   2636
          WRITE (LUOUT,9920) P                                          PRONYX   2637
        END IF                                                          PRONYX   2638
C                                                                       PRONYX   2639
C *** THE NUMBER OF EQUATIONS (ROWS IN M) IS 2*N1 FOR FORWARD-BACKWARD  PRONYX   2640
C     PREDICTION, N1 OTHERWISE                                          PRONYX   2641
        L = N + P - 2                                                   PRONYX   2642
        K = 0                                                           PRONYX   2643
        IF (FLAG.LE.2) K = - P                                          PRONYX   2644
        IF ((FLAG.EQ.2).OR.(FLAG.EQ.3)) L = N - 2                       PRONYX   2645
        N1 = L - K - P + 2                                              PRONYX   2646
C#      IF (BATCH) WRITE (*,9900) N1,N,P,FB,FLAG                        PRONYX   2647
C#        WRITE (LUOUT,9900) N1,N,P,FB,FLAG                             PRONYX   2648
      IF (N1.LE.0) THEN                                                 PRONYX   2649
        MRET=99                                                         PRONYX   2650
        RETURN                                                          PRONYX   2651
      END IF                                                            PRONYX   2652
C                                                                       PRONYX   2653
C *** FORWARD PREDICTION LOGIC                                          PRONYX   2654
        IF (FB.LE.2) THEN                                               PRONYX   2655
          DO 20 I=1,N1                                                  PRONYX   2656
            INDEX = K + P + I - 1                                       PRONYX   2657
            IF (INDEX.LT.0) C(I) = 0.                                   PRONYX   2658
            IF (INDEX.GE.0) C(I) = - X(INDEX)                           PRONYX   2659
            DO 10 J=1,P                                                 PRONYX   2660
              INDEX = K + P - J + I - 1                                 PRONYX   2661
              IF (INDEX.LT.0) M(I,J) = 0.                               PRONYX   2662
              IF (INDEX.GE.0) M(I,J) = X(INDEX)                         PRONYX   2663
 10         CONTINUE                                                    PRONYX   2664
 20       CONTINUE                                                      PRONYX   2665
        END IF                                                          PRONYX   2666
C                                                                       PRONYX   2667
C *** BACKWARD PREDICTION                                               PRONYX   2668
        IF (FB.GE.2) THEN                                               PRONYX   2669
          N2 = N1                                                       PRONYX   2670
          IF (FB.EQ.2) N2 = 2 * N2                                      PRONYX   2671
          DO 40 I=N2-N1+1,N2                                            PRONYX   2672
            DO 30 J=1,P                                                 PRONYX   2673
              INDEX = I + J - N2 - N1 + K - 1                           PRONYX   2674
              IF (INDEX.LT.0) M(I,J) = 0.                               PRONYX   2675
              IF (INDEX.GE.0) M(I,J) = X(INDEX)                         PRONYX   2676
 30         CONTINUE                                                    PRONYX   2677
            INDEX = I - N2 - N1 + K - 1                                 PRONYX   2678
            IF (INDEX.LT.0) C(I) = 0.                                   PRONYX   2679
            IF (INDEX.GE.0) C(I) = - X(INDEX)                           PRONYX   2680
 40       CONTINUE                                                      PRONYX   2681
        END IF                                                          PRONYX   2682
        IF (FB.EQ.2) N1 = 2 * N1                                        PRONYX   2683
C                                                                       PRONYX   2684
C                                                                       PRONYX   2685
C 9900 FORMAT (2X,'##IN DSET_TOEPLITZ: NO. OF EQUATIONS =',I5           PRONYX   2686
C     */8X,'N,P,FB,FLAG =',10I5)                                        PRONYX   2687
 9920 FORMAT (2X,'##IN DSET_TOEPLITZ: NOT ENOUGH DATA, OR ORDER TOO '   PRONYX   2688
     *,'LARGE'/4X,'SETTING P = N-1 = ',I5)                              PRONYX   2689
        RETURN                                                          PRONYX   2690
        END                                                             PRONYX   2691
      SUBROUTINE DSTVDM (M,MAG,PHA,N,Q,IND,IC,DEGREE,NOGO)              PRONYX   2692
CD           This routine fills the real version of the Van der Monde   PRONYX   2693
CD           matrix M with the set of modes Z(i) selected by the        PRONYX   2694
CD           binary vector IC.                                          PRONYX   2695
CD                                                                      PRONYX   2696
C Date : 01/26/89  [JFH]                                                PRONYX   2697
CD Calling Parameters:                                                  PRONYX   2698
CD                                                                      PRONYX   2699
CD   Input :     MAG     : Magnitude of the modes [1,...,Q]             PRONYX   2700
CD               PHA     : Phase of the modes [1,...,Q] (in degrees)    PRONYX   2701
CD               N       : Number of data points.                       PRONYX   2702
CD               Q       : Number of desired modes.                     PRONYX   2703
CD               IND     : Total number of available modes.             PRONYX   2704
CD               IC      : Binary array indicating the posn of the      PRONYX   2705
CD                         modes.                                       PRONYX   2706
CD   Ouput :     M       : Real version of the Van Der Monde Matrix     PRONYX   2707
CD               DEGREE  : Number of columns in M (N by DEGREE)         PRONYX   2708
CD                                                                      PRONYX   2709
      INTEGER NSIZE ,NTAIL                                              PRONYX   2710
      PARAMETER (NSIZE=128)                                             PRONYX   2711
      PARAMETER (NTAIL=1024)                                            PRONYX   2712
      DOUBLE PRECISION MAG(1:*) ,PHA(1:*) ,M(1:NTAIL,1:NSIZE)           PRONYX   2713
      DOUBLE PRECISION TEMP,TEMP1,TEMP2                                 PRONYX   2714
      DOUBLE PRECISION PI,DEGI,PHAJ,A,B,EPS                             PRONYX   2715
      INTEGER Q ,N ,I ,J ,DEGREE ,IND ,IC(1:*)                          PRONYX   2716
      LOGICAL COMP                                                      PRONYX   2717
      INTEGER NOGO                                                      PRONYX   2718
C                                                                       PRONYX   2719
      PI = 4.0D0 * ATAN(1.0D0)                                          PRONYX   2720
      DEGI=PI/180.D0                                                    PRONYX   2721
      EPS = 1.0D-12                                                     PRONYX   2722
      IF (Q.LE.0) THEN                                                  PRONYX   2723
        NOGO=-2                                                         PRONYX   2724
        WRITE(*,'(/10X,A)') ' ERROR ** Not enough modes in SET_VDM'     PRONYX   2725
        RETURN                                                          PRONYX   2726
      ENDIF                                                             PRONYX   2727
C                                                                       PRONYX   2728
C *** FILL ONE COLUMN FOR REAL MODES, TWO COLUMNS FOR COMPLEX           PRONYX   2729
      DEGREE=0                                                          PRONYX   2730
      DO 20 J=1,IND                                                     PRONYX   2731
        IF (IC(J).EQ.1) THEN                                            PRONYX   2732
          TEMP1 = 1.D0                                                  PRONYX   2733
          TEMP2 = 0.D0                                                  PRONYX   2734
          DEGREE = DEGREE + 1                                           PRONYX   2735
          PHAJ=PHA(J)*DEGI                                              PRONYX   2736
          A = DCOS(PHAJ)                                                PRONYX   2737
          B = DSIN(PHAJ)                                                PRONYX   2738
          COMP=DABS(MAG(J)*B).GT.EPS                                    PRONYX   2739
          DO 10 I=1,N                                                   PRONYX   2740
            M(I,DEGREE) = TEMP1                                         PRONYX   2741
            IF (COMP) M(I,DEGREE+1) = TEMP2                             PRONYX   2742
            TEMP = TEMP1                                                PRONYX   2743
            TEMP1 = MAG(J) * (TEMP1 * A - TEMP2 * B)                    PRONYX   2744
            TEMP2 = MAG(J) * (TEMP2 * A + TEMP * B)                     PRONYX   2745
            IF (DABS(TEMP1).GE.1.D30.OR.DABS(TEMP2).GE.1.D30) THEN      PRONYX   2746
              WRITE (*,9000) TEMP1,TEMP2                                PRONYX   2747
              WRITE (6,9000) TEMP1,TEMP2                                PRONYX   2748
            END IF                                                      PRONYX   2749
 10       CONTINUE                                                      PRONYX   2750
          IF (COMP) DEGREE = DEGREE + 1                                 PRONYX   2751
        END IF                                                          PRONYX   2752
 20   CONTINUE                                                          PRONYX   2753
C                                                                       PRONYX   2754
      RETURN                                                            PRONYX   2755
C                                                                       PRONYX   2756
 9000 FORMAT (2X,'##IN DSTVDM: TEMP1,TEMP2 = ',2E12.4)                  PRONYX   2757
C                                                                       PRONYX   2758
      END                                                               PRONYX   2759
C                                                                       PRONYX   2760
        DOUBLE PRECISION FUNCTION DERGY(X,N,M)                          PRONYX   2761
CD                                                                      PRONYX   2762
CD      This routine computes square mean of a vector                   PRONYX   2763
CD      (supposed centered)                                             PRONYX   2764
CD                                                                      PRONYX   2765
C Date : Oct 28 1987                                                    PRONYX   2766
CD Calling Parameters :                                                 PRONYX   2767
CD                                                                      PRONYX   2768
CD   Input :     X     : Vector (0,....N-1)                             PRONYX   2769
CD               N     : Number of values                               PRONYX   2770
CD  Output :    ENERGY : Energy of X = sum (square(X(i))) divided by N  PRONYX   2771
        DOUBLE PRECISION X(0:*)                                         PRONYX   2772
        INTEGER N ,I,M                                                  PRONYX   2773
        M=0                                                             PRONYX   2774
        DERGY = 0.D0                                                    PRONYX   2775
        IF (N.LE.0) THEN                                                PRONYX   2776
          WRITE(*,'(/10X,A,/)')                                         PRONYX   2777
     1    'ERROR Trying to compute ENERGY of zero length vector '       PRONYX   2778
          M=-2                                                          PRONYX   2779
          RETURN                                                        PRONYX   2780
        END IF                                                          PRONYX   2781
        DO 10 I=0,N-1                                                   PRONYX   2782
          DERGY = DERGY + X(I) * X(I)                                   PRONYX   2783
 10     CONTINUE                                                        PRONYX   2784
        DERGY = DERGY / REAL(N)                                         PRONYX   2785
        RETURN                                                          PRONYX   2786
        END                                                             PRONYX   2787
      SUBROUTINE HESORT (A,B,C,D,N,SRULE)                               PRONYX   2788
CD                                                                      PRONYX   2789
CD     This routine re-orders the incoming array in decreasing order    PRONYX   2790
CD                                                                      PRONYX   2791
CD     CALLING PARAMETERS :                                             PRONYX   2792
CD                                                                      PRONYX   2793
CD                     A  - ARRAY TO BE SORTED                          PRONYX   2794
CD                 B,C,D  - ARRAYS TO BE SWAPED ACCORDINGLY             PRONYX   2795
CD                     N  - NO. OF POINTS                               PRONYX   2796
CD                  SRULE  - CONTROL FLAG                               PRONYX   2797
CD                                                                      PRONYX   2798
      DOUBLE PRECISION A(*) ,B(*) ,C(*) ,D(*)                           PRONYX   2799
      INTEGER N,I,SRULE                                                 PRONYX   2800
C                                                                       PRONYX   2801
      IF (SRULE.GT.0) THEN                                              PRONYX   2802
        DO 10 I=1,N                                                     PRONYX   2803
          A(I)=-A(I)                                                    PRONYX   2804
 10     CONTINUE                                                        PRONYX   2805
      END IF                                                            PRONYX   2806
C                                                                       PRONYX   2807
C Initialize the heap properly : partially ordered tree                 PRONYX   2808
      DO 20 I=(N/2),1,-1                                                PRONYX   2809
        CALL DPUSHD (A,B,C,D,I,N)                                       PRONYX   2810
 20   CONTINUE                                                          PRONYX   2811
      DO 30 I=N,2,-1                                                    PRONYX   2812
C Remove minimum from top of tree and put it at the end                 PRONYX   2813
        CALL DEXCH (A(1),A(I),B(1),B(I),C(1),C(I),D(1),D(I))            PRONYX   2814
C Re-establish partially ordered tree property                          PRONYX   2815
        CALL DPUSHD (A,B,C,D,1,I-1)                                     PRONYX   2816
 30   CONTINUE                                                          PRONYX   2817
C                                                                       PRONYX   2818
      IF (SRULE.GT.0) THEN                                              PRONYX   2819
        DO 40 I=1,N                                                     PRONYX   2820
          A(I)=-A(I)                                                    PRONYX   2821
 40     CONTINUE                                                        PRONYX   2822
      END IF                                                            PRONYX   2823
      RETURN                                                            PRONYX   2824
C                                                                       PRONYX   2825
      END                                                               PRONYX   2826
C                                                                       PRONYX   2827
      SUBROUTINE DPUSHD (A,B,C,D,FIRST,LAST)                            PRONYX   2828
CD                                                                      PRONYX   2829
CD       This routine re-order the incoming array A as follows:         PRONYX   2830
CD                                                                      PRONYX   2831
CD         Assumes that A[FIRST],...,A[LAST] obeys partially ordered    PRONYX   2832
CD         tree property, except possibly for the children of A[FIRST]. PRONYX   2833
CD         The subroutine pushes down A[FIRST] down until the partially PRONYX   2834
CD         ordered tree property is restored.                           PRONYX   2835
CD                                                                      PRONYX   2836
CD  Calling Parameters:                                                 PRONYX   2837
CD                                                                      PRONYX   2838
CD     Input :  A            : Array to be sorted                       PRONYX   2839
CD              B,C,D        : Arrays to be swaped accordingly          PRONYX   2840
CD              FIRST,LAST   : Index range to be re-aranged             PRONYX   2841
CD    Output :  A,B,C,D      : Re-arranged arrays                       PRONYX   2842
CD                                                                      PRONYX   2843
      DOUBLE PRECISION A(*) ,B(*) ,C(*) ,D(*)                           PRONYX   2844
      INTEGER FIRST ,LAST ,R ,R2                                        PRONYX   2845
                                                                        PRONYX   2846
C Initialization                                                        PRONYX   2847
      R = FIRST                                                         PRONYX   2848
C#    DO WHILE (R .LE. (LAST/2))                                        PRONYX   2849
 100  IF (R .GT. (LAST/2)) GO TO 110                                    PRONYX   2850
        R2 = 2*R                                                        PRONYX   2851
        IF (LAST .EQ. R2) THEN                                          PRONYX   2852
C Then R has one child at R2                                            PRONYX   2853
          IF (A(R) .GT. A(R2)) CALL DEXCH (A(R),A(R2),B(R),B(R2),       PRONYX   2854
     *         C(R),C(R2),D(R),D(R2))                                   PRONYX   2855
          R = LAST                                                      PRONYX   2856
          R2 = 2*R                                                      PRONYX   2857
C Otherwise, R has two children, elements at R2 and R2+1                PRONYX   2858
        ELSE                                                            PRONYX   2859
          IF ((A(R) .GT. A(R2)) .AND. (A(R2) .LE. A(R2+1))) THEN        PRONYX   2860
C Swap R with left child                                                PRONYX   2861
            CALL DEXCH (A(R),A(R2),B(R),B(R2),C(R),C(R2),D(R),D(R2))    PRONYX   2862
            R = 2*R                                                     PRONYX   2863
            R2 = 2*R                                                    PRONYX   2864
          ELSE IF ((A(R).GT.A(R2+1)) .AND. (A(R2+1).LT.A(R2))) THEN     PRONYX   2865
C Swap R with right child                                               PRONYX   2866
         CALL DEXCH(A(R),A(R2+1),B(R),B(R2+1),C(R),C(R2+1),D(R),D(R2+1))PRONYX   2867
            R = 2*R + 1                                                 PRONYX   2868
            R2 = 2*R                                                    PRONYX   2869
          ELSE                                                          PRONYX   2870
C R does not violate partially ordered tree property, therefore break loPRONYX   2871
            R = LAST                                                    PRONYX   2872
          END IF                                                        PRONYX   2873
        END IF                                                          PRONYX   2874
      GOTO 100                                                          PRONYX   2875
C#    END DO                                                            PRONYX   2876
 110  RETURN                                                            PRONYX   2877
      END                                                               PRONYX   2878
      SUBROUTINE DEXCH (X,Y,A,B,C,D,E,F)                                PRONYX   2879
CD                                                                      PRONYX   2880
CD          This routine exchanges the values of X and Y using          PRONYX   2881
CD          a temporary variable TEMP. Do the same for                  PRONYX   2882
CD          A and B, C and D, and E and F.                              PRONYX   2883
CD                                                                      PRONYX   2884
CD    Calling Parameters :                                              PRONYX   2885
CD                                                                      PRONYX   2886
CD      Input and Outputs : X,Y      : The real values to be swaped     PRONYX   2887
CD                          A,B      : The real values to be swaped     PRONYX   2888
CD                          C,D      : The real values to be swaped     PRONYX   2889
CD                          E,F      : The real values to be swaped     PRONYX   2890
CD                                                                      PRONYX   2891
      DOUBLE PRECISION X,Y,A,B,C,D,E,F,TEMP                             PRONYX   2892
      TEMP = X                                                          PRONYX   2893
      X = Y                                                             PRONYX   2894
      Y = TEMP                                                          PRONYX   2895
      TEMP = A                                                          PRONYX   2896
      A = B                                                             PRONYX   2897
      B = TEMP                                                          PRONYX   2898
      TEMP = C                                                          PRONYX   2899
      C = D                                                             PRONYX   2900
      D = TEMP                                                          PRONYX   2901
      TEMP = E                                                          PRONYX   2902
      E = F                                                             PRONYX   2903
      F = TEMP                                                          PRONYX   2904
      RETURN                                                            PRONYX   2905
      END                                                               PRONYX   2906
      SUBROUTINE PRSPAK (XSIG,TSTEP,NFSAVE,NSTART                       PRONYX   2907
     *,QCON,SNR,SHAT,DAMP,FRQHZ,RELMAG,PHASE,AMPOPT,PHIOPT              PRONYX   2908
     *,LPOORG,LPMORG,FBSAVE,PIRORG,QSAVE,NHSAVE                         PRONYX   2909
     *,TRMLPM,FTRIMH,FTRIML                                             PRONYX   2910
     *,LUIN,LUOUT,LURES,MRETR,DBP,NOGO)                                 PRONYX   2911
CD                                                                      PRONYX   2912
CD       THIS ROUTINE IS THE MAIN ROUTINE PERFORMS THE SPECTRAL         PRONYX   2913
CD       ANALYSIS                                                       PRONYX   2914
CD                                                                      PRONYX   2915
CD       CALLING PARAMETERS :                                           PRONYX   2916
CD                                                                      PRONYX   2917
CD           XSIG   - ORIGINAL Y-AXIS ARRAY                             PRONYX   2918
CD           TSTEP  - TIME INTERVAL                                     PRONYX   2919
CD           NFSAVE - NUMBER OF DATA POINT                              PRONYX   2920
CD           NSTART - NOT USED                                          PRONYX   2921
CD           QCON   - NUMBER OF MODE                                    PRONYX   2922
CD           SNR    - SIGNAL TO NOISE RATIO                             PRONYX   2923
CD           SHAT   - CALCULATED Y-AXIS ARRAY AFTER PRONY ANALYSIS      PRONYX   2924
CD           DAMP   - THE DAMPING RATIO ARRAY OF THE MODES              PRONYX   2925
CD           FRQHZ  - THE NATURAL FREQUENCIES ARRAY OF THE MODES        PRONYX   2926
CD           RELMAG - THE RELATIVE MAGNITUDE ARRAY OF THE MODES         PRONYX   2927
CD           PHASE  - THE PHASE ANGLES ARRAY OF THE MODES               PRONYX   2928
CD           AMPOPT - AMPITUDE ARRAY                                    PRONYX   2929
CD           PHIOPT - PAHSE ARRAY                                       PRONYX   2930
CD           LPOORG - CONTROL FLAG FOR LINEAR PREDICTION ORDER          PRONYX   2931
CD           LPMORG - CONTROL FLAG FOR LINEAR PREDICTION METHOD         PRONYX   2932
CD           FBSAVE - CONTROL FLAG FOR BACKWARD/FORWARD COTROL          PRONYX   2933
CD           PIRORG -                                                   PRONYX   2934
CD           QSAVE  - NUMBER OF MODE                                    PRONYX   2935
CD           NHSAVE -                                                   PRONYX   2936
CD           TRMLPM - MODE TRIM IN %                                    PRONYX   2937
CD           FTRIMH - UPPER BOUND FOR THE MODE TRIM                     PRONYX   2938
CD           FTRIML - LOWER BOUND FOR THE MODE TRIM                     PRONYX   2939
CD     LUIN ,LUOUT  - POINTER FOR THE TERMINATION KEY                   PRONYX   2940
CD           LURES  - LOGICAL UNIT FOR DEBUG OUTPUT                     PRONYX   2941
CD           MRETR  - TERMINATION LOGIC                                 PRONYX   2942
CD           DBP    - LOGICAL FLAG FOR DEBUG                            PRONYX   2943
CD           NOGO   - I/O ERROR FLAG                                    PRONYX   2944
CD                                                                      PRONYX   2945
      INCLUDE 'CHDISP.INS'                                              PRONYX   2946
C                                                                       PRONYX   2947
      INTEGER NSIZE, NTAIL                                              PRONYX   2948
      PARAMETER (NSIZE=128)                                             PRONYX   2949
      PARAMETER (NTAIL=1024)                                            PRONYX   2950
      DOUBLE PRECISION XSIG(0:*),SHAT(0:*),SNR,SIGEN                    PRONYX   2951
      REAL TSTEP                                                        PRONYX   2952
      DOUBLE PRECISION DAMP(1:*),FRQHZ(1:*),RELMAG(1:*),PHASE(1:*),     PRONYX   2953
     +                 AMPREF                                           PRONYX   2954
      DOUBLE PRECISION FTRIMH,FTRIML,TRMLPM,ZTRIMH,ZTRIML,ZZTRIM        PRONYX   2955
      DOUBLE PRECISION FAC,PCDAMP                                       PRONYX   2956
      INTEGER NNHAT,MRETR,MRETC                                         PRONYX   2957
      INTEGER NFSAVE,NSTART,NHSAVE                                      PRONYX   2958
      INTEGER LPOORG,LPMORG,FBSAVE,PIRORG,QSAVE,DQPSET                  PRONYX   2959
      INTEGER LPOMAX                                                    PRONYX   2960
      LOGICAL BATCH,TRACK,PIRADJ                                        PRONYX   2961
      INTEGER ILAST,QTRIM,SRULE                                         PRONYX   2962
      INTEGER LUIN,LUOUT,LURES,PIRCON,I,NNFIT,FBCON                     PRONYX   2963
      INTEGER LPOCON,QCON,DMODES,DMREAL,DMCMPX,N1,J,LPMCON,K,NINIT,INFO PRONYX   2964
      INTEGER IC(1:NSIZE)                                               PRONYX   2965
      INTEGER ICOUNT,IQ,IFAIL                                           PRONYX   2966
      DOUBLE PRECISION SIGDIF(0:NTAIL-1)                                PRONYX   2967
      DOUBLE PRECISION MAT(1:NTAIL,1:NSIZE),WORK(1:NTAIL),SUM           PRONYX   2968
      DOUBLE PRECISION DERGY,ALPHA,PI,DEG,DEGI,EPS1,HZFAC               PRONYX   2969
      DOUBLE PRECISION U(1:NTAIL,1:NSIZE),V(1:NSIZE,1:NSIZE)            PRONYX   2970
      DOUBLE PRECISION A(0:NSIZE-1),B(1:NSIZE),C(1:NTAIL),              PRONYX   2971
     +                 DMATSV(1:NSIZE)                                  PRONYX   2972
      DOUBLE PRECISION OPR(1:NSIZE),ZEROR(1:NSIZE),ZEROI(1:NSIZE)       PRONYX   2973
      DOUBLE PRECISION ROOTR(1:NSIZE),ROOTI(1:NSIZE),MAG(1:NSIZE),      PRONYX   2974
     +                 PHA(1:NSIZE)                                     PRONYX   2975
      DOUBLE PRECISION AMP(1:NSIZE),PHI(1:NSIZE),DERR,DMNERR            PRONYX   2976
      DOUBLE PRECISION AMPOPT(1:NSIZE),PHIOPT(1:NSIZE)                  PRONYX   2977
      DOUBLE PRECISION MAGOPT(1:NSIZE),PHAOPT(1:NSIZE)                  PRONYX   2978
C                                                                       PRONYX   2979
C! FLF                                                                  PRONYX   2980
      CHARACTER*1 QQ                                                    PRONYX   2981
      INTEGER LINES1,NOGO                                               PRONYX   2982
      LOGICAL DBP                                                       PRONYX   2983
      DOUBLE PRECISION RPART,QPART,DFLF1,DFLF2                          PRONYX   2984
C                                                                       PRONYX   2985
      COMMON / PRSCOM / SIGDIF,WORK,U,V,A,B,C,DMATSV,OPR,ZEROR,ZEROI,   PRONYX   2986
     +                  ROOTR,ROOTI,PHA,AMP,PHI,PHAOPT,                 PRONYX   2987
     +                  IC,MAT,MAG,MAGOPT                               PRONYX   2988
C! FLF                                                                  PRONYX   2989
                                                                        PRONYX   2990
C *** MATHEMATICAL CONSTANTS                                            PRONYX   2991
      INTEGER IZERO,IONE                                                PRONYX   2992
C                                                                       PRONYX   2993
C *** SET PRECISION CONTROLS                                            PRONYX   2994
C     [LATER VERSIONS WILL ENTER THESE AS SPECIAL PROCESSING CONTROLS]  PRONYX   2995
C-----------------------------------------------C                       PRONYX   2996
C     FOR VAX, ROOTING VIA 'QPOLRT'             C                       PRONYX   2997
C     DQPSET=0                                                          PRONYX   2998
C     FOR COMPUTERS WITHOUT EXTENDED PRECISION, C                       PRONYX   2999
C     ROOTING VIA 'DPOLRT'                      C                       PRONYX   3000
C     DQPSET=1                                  C                       PRONYX   3001
*IF -VAX                                                                PRONYX   3002
      DQPSET=1                                                          PRONYX   3003
*ENDIF                                                                  PRONYX   3004
*IF VAX                                                                 PRONYX   3005
*ENDIF                                                                  PRONYX   3007
C-----------------------------------------------C                       PRONYX   3008
C                                                                       PRONYX   3009
C *** SET MATHEMATICAL CONSTANTS                                        PRONYX   3010
      IZERO=0                                                           PRONYX   3011
      IONE=1                                                            PRONYX   3012
      PI=4.D0*DATAN(1.D0)                                               PRONYX   3013
      DEG=180.D0/PI                                                     PRONYX   3014
      DEGI=1.D0/DEG                                                     PRONYX   3015
      HZFAC=1./(360.*TSTEP)                                             PRONYX   3016
      EPS1=1.D-8                                                        PRONYX   3017
C                                                                       PRONYX   3018
C *** CHECK STORAGE AND DISPLAY CONTROLS                                PRONYX   3019
C     [LATER VERSIONS WILL ENTER THESE AS SPECIAL PROCESSING CONTROLS]  PRONYX   3020
      BATCH=.TRUE.                                                      PRONYX   3021
      TRACK=.FALSE.                                                     PRONYX   3022
C     TRACK=.TRUE.                                                      PRONYX   3023
C                                                                       PRONYX   3024
C *** SET INTEGER CONTROLS                                              PRONYX   3025
      NNFIT=NFSAVE                                                      PRONYX   3026
      IF (NNFIT.LE.0) THEN                                              PRONYX   3027
        IF(DBP)WRITE (LURES,9940) NNFIT                                 PRONYX   3028
        WRITE (LUOUT,9940) NNFIT                                        PRONYX   3029
        GO TO 8000                                                      PRONYX   3030
      END IF                                                            PRONYX   3031
      NINIT=NSTART                                                      PRONYX   3032
      LPOCON=LPOORG                                                     PRONYX   3033
      LPOMAX=NNFIT/2                                                    PRONYX   3034
      IF (LPOCON.LT.1) THEN                                             PRONYX   3035
        LPOCON=11*LPOMAX/12                                             PRONYX   3036
        IF(DBP)WRITE (LURES,9060) LPOCON                                PRONYX   3037
      END IF                                                            PRONYX   3038
      IF (LPOCON.GT.LPOMAX) THEN                                        PRONYX   3039
        LPOCON=11*LPOMAX/12                                             PRONYX   3040
        IF(DBP)WRITE (LURES,9062) LPOMAX,LPOCON                         PRONYX   3041
      END IF                                                            PRONYX   3042
      IF (LPOCON.GT.NSIZE) THEN                                         PRONYX   3043
        LPOCON=NSIZE                                                    PRONYX   3044
        IF(DBP)WRITE (LURES,9064) NSIZE,LPOCON                          PRONYX   3045
      END IF                                                            PRONYX   3046
      IF (LPOMAX.GT.NSIZE) LPOMAX=NSIZE                                 PRONYX   3047
      LPMCON=LPMORG                                                     PRONYX   3048
      IF (LPMCON.LT.1.OR.LPMCON.GT.4) LPMCON=3                          PRONYX   3049
      FBCON=FBSAVE                                                      PRONYX   3050
      IF (FBCON.LT.1.OR.FBCON.GT.3) FBCON=3                             PRONYX   3051
      PIRCON=PIRORG                                                     PRONYX   3052
      IF (PIRCON.LT.1.OR.PIRCON.GT.LPOCON) PIRCON=15*LPOCON/16          PRONYX   3053
      QCON=QSAVE                                                        PRONYX   3054
      IF (QCON.LE.0.0) QCON=LPOCON                                      PRONYX   3055
      NNHAT=NHSAVE                                                      PRONYX   3056
      IF (NNHAT.LT.NNFIT) NNHAT=NNFIT                                   PRONYX   3057
      QTRIM=0                                                           PRONYX   3058
C                                                                       PRONYX   3059
C *** DIAGNOSTIC                                                        PRONYX   3060
      IF(DBP)WRITE (LURES,9910) NNFIT,NINIT,QCON,LPOCON,LPMCON          PRONYX   3061
     *,FBCON,PIRCON,MRETR                                               PRONYX   3062
      IF(DBP)WRITE (LURES,9912) TSTEP                                   PRONYX   3063
C                                                                       PRONYX   3064
C *** TEST FOR NULL SIGNAL                                              PRONYX   3065
      SIGEN=DERGY(XSIG,NNFIT,NOGO)                                      PRONYX   3066
      IF(NOGO.LT.0) RETURN                                              PRONYX   3067
      IF(DBP)WRITE (LURES,9914) SIGEN                                   PRONYX   3068
      IF (SIGEN.LT.1.D-12) THEN                                         PRONYX   3069
        IF(DBP)WRITE (LURES,9920) SIGEN                                 PRONYX   3070
        WRITE (LUOUT,9920) SIGEN                                        PRONYX   3071
        DMODES=0                                                        PRONYX   3072
        GO TO 8000                                                      PRONYX   3073
      END IF                                                            PRONYX   3074
C                                                                       PRONYX   3075
C *** SIGNAL DISPLAY                                                    PRONYX   3076
      IF (TRACK) THEN                                                   PRONYX   3077
        IF(DBP)THEN                                                     PRONYX   3078
           WRITE (LURES,9100)                                           PRONYX   3079
           DO 10 K=0,NNFIT-1                                            PRONYX   3080
             WRITE (LURES,9622) K,XSIG(K)                               PRONYX   3081
 10        CONTINUE                                                     PRONYX   3082
        ENDIF                                                           PRONYX   3083
      END IF                                                            PRONYX   3084
      IF(DBP)WRITE (LURES,9120) LPOCON                                  PRONYX   3085
      IF(DBP)WRITE (LURES,9130) LPMCON                                  PRONYX   3086
      LPOORG=LPOCON                                                     PRONYX   3087
      LPMORG=LPMCON                                                     PRONYX   3088
      IF(DBP)WRITE (LURES,9140) FBCON                                   PRONYX   3089
      FBSAVE=FBCON                                                      PRONYX   3090
C                                                                       PRONYX   3091
C *** FORM LEAST SQUARES TOEPLITZ SYSTEM FOR DENOMINATOR COEFFICIENTS   PRONYX   3092
      MRETC=0                                                           PRONYX   3093
      CALL DSETTP (XSIG,MAT,C,NNFIT,LPOCON,N1,FBCON,LPMCON,LURES        PRONYX   3094
     1,BATCH,MRETC)                                                     PRONYX   3095
      IF (MRETC.GT.0) GO TO 8000                                        PRONYX   3096
C                                                                       PRONYX   3097
C *** COMPUTE THE SVD OF THE MATRIX MAT=U*DIAG(DMATSV)*TRANS(V)         PRONYX   3098
      CALL DSVDC (MAT,NTAIL,N1,LPOCON,DMATSV,B,U,NTAIL,V                PRONYX   3099
     *,NSIZE,WORK,21,INFO)                                              PRONYX   3100
      IF(DBP)WRITE (LURES,9330) DMATSV(1),DMATSV(LPOCON)                PRONYX   3101
      IF(DBP)WRITE (LURES,9200)                                         PRONYX   3102
      IF(DBP) THEN                                                      PRONYX   3103
        DO 20 I=1,LPOCON                                                PRONYX   3104
          WRITE (LURES,9622) I,DMATSV(I)                                PRONYX   3105
 20     CONTINUE                                                        PRONYX   3106
      ENDIF                                                             PRONYX   3107
C                                                                       PRONYX   3108
      IF(DBP)WRITE (LURES,9320) LPOCON                                  PRONYX   3109
      IF(DBP)WRITE (LURES,9322) PIRCON                                  PRONYX   3110
C                                                                       PRONYX   3111
C *** SOLVE LEAST SQUARE SYSTEM TO GET THE VECTOR A(0:LPOCON), A(0)=1.  PRONYX   3112
  100 CALL DSQSVD (A,U,V,DMATSV,C,PIRCON,N1,LPOCON,LURES,BATCH,NOGO)    PRONYX   3113
C FLF                                                                   PRONYX   3114
      IF(NOGO.LT.0) RETURN                                              PRONYX   3115
CFLF                                                                    PRONYX   3116
      IF (TRACK) THEN                                                   PRONYX   3117
         IF(DBP)WRITE (LURES,9450) LPOCON,PIRCON                        PRONYX   3118
      ENDIF                                                             PRONYX   3119
C                                                                       PRONYX   3120
C *** DIAGNOSTIC DISPLAY OF DENOMINATOR Z-COEFFICIENTS                  PRONYX   3121
      IF (TRACK) THEN                                                   PRONYX   3122
        IF(DBP)THEN                                                     PRONYX   3123
           WRITE (LURES,9451)                                           PRONYX   3124
           DO 30 I=0,LPOCON                                             PRONYX   3125
            WRITE (LURES,9622) I,A(I)                                   PRONYX   3126
 30        CONTINUE                                                     PRONYX   3127
        ENDIF                                                           PRONYX   3128
      END IF                                                            PRONYX   3129
C                                                                       PRONYX   3130
C *** ROOT DENOMINATOR POLYNOMIAL FOR FREQUENCIES AND DAMPING FACTORS   PRONYX   3131
      DO 40 I=0,LPOCON                                                  PRONYX   3132
        OPR(I+1)=A(LPOCON-I)                                            PRONYX   3133
C                            note reversal                              PRONYX   3134
 40   CONTINUE                                                          PRONYX   3135
*IF VAX                                                                 PRONYX   3136
*ENDIF                                                                  PRONYX   3145
*IF -VAX                                                                PRONYX   3146
      IF (DQPSET.GE.1) THEN                                             PRONYX   3147
        CALL DPOLRT (OPR,LPOCON,ZEROR,ZEROI,IFAIL)                      PRONYX   3148
        IF (IFAIL.NE.0) THEN                                            PRONYX   3149
          WRITE (LURES,*) 'DPOLY FAILURE =',IFAIL                       PRONYX   3150
     1,   ' FOR DENOMINATOR'                                            PRONYX   3151
          WRITE (LUOUT,*) 'DPOLY FAILURE =',IFAIL,' FOR DENOMINATOR'    PRONYX   3152
        END IF                                                          PRONYX   3153
      END IF                                                            PRONYX   3154
*ENDIF                                                                  PRONYX   3155
C                                                                       PRONYX   3156
C *** DIAGNOSTIC DISPLAY OF DONOMINATOR Z-ROOTS                         PRONYX   3157
      IF (TRACK) THEN                                                   PRONYX   3158
        IF(DBP)THEN                                                     PRONYX   3159
           WRITE (LURES,9210)                                           PRONYX   3160
           DO 50 I=1,QCON                                               PRONYX   3161
              WRITE (LURES,9622) I,ZEROR(I),ZEROI(I)                    PRONYX   3162
 50        CONTINUE                                                     PRONYX   3163
        ENDIF                                                           PRONYX   3164
      END IF                                                            PRONYX   3165
C                                                                       PRONYX   3166
C *** SIFT OUT COMPLEX DENOMINATOR ROOTS IN LOWER HALF OF Z-PLANE       PRONYX   3167
C *** SIFT OUT OVERLY LARGE ROOTS, RESTARTING IF APPROPRIATE            PRONYX   3168
C *** COMPUTE MAGNITUDE AND PHASE OF RETAINED Z-ROOTS AND ORDER         PRONYX   3169
C     BY DECREASING MAGNITUDE                                           PRONYX   3170
      ZZTRIM=10.**(30.D0/NNFIT)                                         PRONYX   3171
      PIRADJ=.FALSE.                                                    PRONYX   3172
      DMODES=0                                                          PRONYX   3173
      DMREAL=0                                                          PRONYX   3174
      DMCMPX=0                                                          PRONYX   3175
      DO 120 I=1,LPOCON                                                 PRONYX   3176
      IF (ZEROI(I).LT.0.D0) GO TO 120                                   PRONYX   3177
      DMODES=DMODES+1                                                   PRONYX   3178
      ROOTR(DMODES) = ZEROR(I)                                          PRONYX   3179
      ROOTI(DMODES) = ZEROI(I)                                          PRONYX   3180
      IF (DABS(ROOTI(DMODES)).LT.EPS1) ROOTI(DMODES) = 0.D0             PRONYX   3181
      MAG(DMODES) = DSQRT(ZEROR(I)**2 + ZEROI(I)**2)                    PRONYX   3182
      PHA(DMODES) = DABS(DATAN2(ZEROI(I),ZEROR(I)))*DEG                 PRONYX   3183
      IF (PHA(DMODES).LT.EPS1) THEN                                     PRONYX   3184
        ZEROI(DMODES)=0.D0                                              PRONYX   3185
        PHA(DMODES)=0.D0                                                PRONYX   3186
      END IF                                                            PRONYX   3187
      IF (MAG(DMODES).LT.1.D0) GO TO 110                                PRONYX   3188
      IF(DBP)WRITE (LURES,9212) MAG(DMODES),PHA(DMODES)                 PRONYX   3189
      IF (MAG(DMODES).GT.ZZTRIM) THEN                                   PRONYX   3190
        WRITE (LURES,9214)                                              PRONYX   3191
        IF (PIRADJ.AND.PIRCON.GE.LPOCON-6) THEN                         PRONYX   3192
          PIRCON=PIRCON-2                                               PRONYX   3193
          IF(DBP)WRITE (LURES,9220) PIRCON                              PRONYX   3194
          GO TO 100                                                     PRONYX   3195
        ELSE                                                            PRONYX   3196
          IF(DBP)WRITE (LURES,9222)                                     PRONYX   3197
          DMODES=DMODES-1                                               PRONYX   3198
          GO TO 120                                                     PRONYX   3199
        END IF                                                          PRONYX   3200
      END IF                                                            PRONYX   3201
  110 IF (ZEROI(DMODES).EQ.0.D0) DMREAL=DMREAL+1                        PRONYX   3202
      IF (ZEROI(DMODES).NE.0.D0) DMCMPX=DMCMPX+1                        PRONYX   3203
  120 CONTINUE                                                          PRONYX   3204
      SRULE=-1                                                          PRONYX   3205
      CALL HESORT (MAG,PHA,ROOTR,ROOTI,DMODES,SRULE)                    PRONYX   3206
      IF(DBP)WRITE (LURES,9452) LPOCON,DMODES,DMREAL,DMCMPX             PRONYX   3207
      QCON=DMODES                                                       PRONYX   3208
C                                                                       PRONYX   3209
      IF(DBP)THEN                                                       PRONYX   3210
        WRITE (LURES,9210)                                              PRONYX   3211
        DO 60 I=1,QCON                                                  PRONYX   3212
          WRITE (LURES,9622) I,ROOTR(I),ROOTI(I)                        PRONYX   3213
 60     CONTINUE                                                        PRONYX   3214
      ENDIF                                                             PRONYX   3215
C                                                                       PRONYX   3216
C *** COMPUTE THE NUMBER OF SETS OF SIZE QCON OUT OF DMODES, WITH A LITTPRONYX   3217
C     TRICK TO AVOID INTEGER OVERFLOW                                   PRONYX   3218
      ICOUNT=1                                                          PRONYX   3219
      IQ=QCON                                                           PRONYX   3220
      IF (IQ.GT.(DMODES/2)) IQ=DMODES-QCON                              PRONYX   3221
      DO 70 I=0,IQ-1                                                    PRONYX   3222
        ICOUNT=ICOUNT*(DMODES-I)/(I+1)                                  PRONYX   3223
 70   CONTINUE                                                          PRONYX   3224
      IF(DBP)WRITE (LURES,*) 'ICOUNT =',ICOUNT                          PRONYX   3225
        DO 80 I=1,QCON                                                  PRONYX   3226
          IC(I)=1                                                       PRONYX   3227
 80     CONTINUE                                                        PRONYX   3228
        DO 90 I=QCON+1,DMODES                                           PRONYX   3229
          IC(I)=0                                                       PRONYX   3230
 90     CONTINUE                                                        PRONYX   3231
        DMNERR=DERGY(XSIG,NNFIT,NOGO)*REAL(NNFIT)                       PRONYX   3232
        IF(NOGO.LT.0)RETURN                                             PRONYX   3233
        DO 150 I=1,ICOUNT                                               PRONYX   3234
        CALL DSQVDM (XSIG,ROOTR,ROOTI,MAG,PHA,AMP,PHI,NNFIT,QCON        PRONYX   3235
     *  ,DMODES,IC,DERR,NOGO)                                           PRONYX   3236
        IF(NOGO.LT.0)RETURN                                             PRONYX   3237
        IF (DERR.LT.DMNERR) THEN                                        PRONYX   3238
          DMNERR=DERR                                                   PRONYX   3239
          J=1                                                           PRONYX   3240
          DO 140 K=1,DMODES                                             PRONYX   3241
            IF (IC(K).EQ.1) THEN                                        PRONYX   3242
              MAGOPT(J)=MAG(K)                                          PRONYX   3243
              PHAOPT(J)=PHA(K)                                          PRONYX   3244
              AMPOPT(J)=AMP(K)                                          PRONYX   3245
              PHIOPT(J)=PHI(K)                                          PRONYX   3246
              J=J+1                                                     PRONYX   3247
            END IF                                                      PRONYX   3248
 140      CONTINUE                                                      PRONYX   3249
        END IF                                                          PRONYX   3250
        CALL COMBIN (DMODES,QCON,IC)                                    PRONYX   3251
 150  CONTINUE                                                          PRONYX   3252
C                                                                       PRONYX   3253
C *** FORCE ALL MODE AMPLITUDES POSITIVE                                PRONYX   3254
      DO 160 I=1,QCON                                                   PRONYX   3255
      IF (AMPOPT(I).LT.0.0) THEN                                        PRONYX   3256
        AMPOPT(I)=-AMPOPT(I)                                            PRONYX   3257
        PHIOPT(I)=PHIOPT(I)+180.0                                       PRONYX   3258
      END IF                                                            PRONYX   3259
 160  CONTINUE                                                          PRONYX   3260
C                                                                       PRONYX   3261
C *** ORDER MODES BY DECREASING AMPLITUDE                               PRONYX   3262
      SRULE=-1                                                          PRONYX   3263
      CALL HESORT (AMPOPT,PHIOPT,MAGOPT,PHAOPT,QCON,SRULE)              PRONYX   3264
C                                                                       PRONYX   3265
C *** FIND STRONGEST MODE NOT AT SAMPLING FREQUENCY,                    PRONYX   3266
C     AND OUTSIDE INTERVAL [FTRIML,FTRIMH]                              PRONYX   3267
      ZTRIMH=FTRIMH/HZFAC                                               PRONYX   3268
      IF (ZTRIMH.LE.0.0) ZTRIMH=360.0                                   PRONYX   3269
      ZTRIML=FTRIML/HZFAC                                               PRONYX   3270
      AMPREF=AMPOPT(1)                                                  PRONYX   3271
      DO 300 I=1,QCON                                                   PRONYX   3272
      IF (PHAOPT(I).LT.ZTRIML.OR.PHAOPT(I).GT.ZTRIMH) GO TO 300         PRONYX   3273
      IF (DABS(PHAOPT(I)-180.0).GT.1.D-6) THEN                          PRONYX   3274
        AMPREF=AMPOPT(I)                                                PRONYX   3275
        GO TO 310                                                       PRONYX   3276
      END IF                                                            PRONYX   3277
  300 CONTINUE                                                          PRONYX   3278
  310 CONTINUE                                                          PRONYX   3279
C                                                                       PRONYX   3280
C *** TEMPORARY DIAGNOSTIC                                              PRONYX   3281
      IF(DBP)THEN                                                       PRONYX   3282
        WRITE (LURES,9160) QCON                                         PRONYX   3283
        WRITE (LURES,9618)                                              PRONYX   3284
        DO 320 I=1,QCON                                                 PRONYX   3285
           WRITE (LURES,9622) I,AMPOPT(I)                               PRONYX   3286
 320    CONTINUE                                                        PRONYX   3287
      ENDIF                                                             PRONYX   3288
C                                                                       PRONYX   3289
C      WRITE (LURES,*) 'QCON = ',QCON                                   PRONYX   3290
C FLF                                                                   PRONYX   3291
      WRITE (LURES,*) 'TOTAL NUMBER OF MODES FITTED : ',QCON            PRONYX   3292
C FLF                                                                   PRONYX   3293
      IF (TRMLPM.GT.0.0) THEN                                           PRONYX   3294
         IF(DBP)WRITE (LURES,9460) TRMLPM                               PRONYX   3295
      ENDIF                                                             PRONYX   3296
      IF(DBP)WRITE (LURES,9620)                                         PRONYX   3297
      DO 330 I=1,QCON                                                   PRONYX   3298
        RELMAG(I)=AMPOPT(I)/AMPREF                                      PRONYX   3299
        IF (RELMAG(I).LE.TRMLPM) GO TO 340                              PRONYX   3300
        ILAST=I                                                         PRONYX   3301
        IF(DBP)WRITE (LURES,9622) I,RELMAG(I)                           PRONYX   3302
 330  CONTINUE                                                          PRONYX   3303
 340  WRITE (LURES,9470) ILAST                                          PRONYX   3304
      QCON=ILAST                                                        PRONYX   3305
      QTRIM=QCON                                                        PRONYX   3306
C                                                                       PRONYX   3307
C *** COMPUTE THE RECONSTRUCTED AND PREDICTED SIGNAL SEGMENTS FROM THE  PRONYX   3308
C     QCON SELECTED MODES                                               PRONYX   3309
C *** COMPUTE THE FITTING ERROR AND SNR (OVER NNFIT SAMPLES ONLY)       PRONYX   3310
      DO 360 I=0,NNHAT-1                                                PRONYX   3311
        SUM=0.0                                                         PRONYX   3312
        DO 350 K=1,QCON                                                 PRONYX   3313
          ALPHA=(PHAOPT(K)*I+PHIOPT(K) )*DEGI                           PRONYX   3314
          SUM=SUM+AMPOPT(K)*DCOS(ALPHA)*MAGOPT(K)**I                    PRONYX   3315
 350    CONTINUE                                                        PRONYX   3316
        SHAT(I)=SUM                                                     PRONYX   3317
        SIGDIF(I)=XSIG(I)-SHAT(I)                                       PRONYX   3318
 360  CONTINUE                                                          PRONYX   3319
C FLF                                                                   PRONYX   3320
      DFLF1=DERGY(SHAT,NNFIT,NOGO)                                      PRONYX   3321
      IF(NOGO.LT.0) RETURN                                              PRONYX   3322
      DFLF2=DERGY(SIGDIF,NNFIT,NOGO)                                    PRONYX   3323
      IF(NOGO.LT.0) RETURN                                              PRONYX   3324
C      SNR=10.*DLOG10(DERGY(SHAT,NNFIT)/DERGY(SIGDIF,NNFIT))            PRONYX   3325
      SNR=10.*DLOG10(DFLF1/DFLF2)                                       PRONYX   3326
C FLF                                                                   PRONYX   3327
      IF(DBP)WRITE (LURES,9600) SNR,QCON                                PRONYX   3328
C                                                                       PRONYX   3329
C *** CONVERT Z-DOMAIN MODES TO S-DOMAIN                                PRONYX   3330
      DO 370 I=1,QCON                                                   PRONYX   3331
C!        DAMP(I)=-1.0/TSTEP*DLOG(MAGOPT(I))/(2.*PI)                    PRONYX   3332
        DAMP(I)=DLOG(MAGOPT(I))/TSTEP                                   PRONYX   3333
        FRQHZ(I)=PHAOPT(I)*HZFAC                                        PRONYX   3334
        PHASE(I)=PHIOPT(I)                                              PRONYX   3335
 370  CONTINUE                                                          PRONYX   3336
C                                                                       PRONYX   3337
C *** ORDER MODES BY INCREASING NATURAL FREQUENCY                       PRONYX   3338
      SRULE=1                                                           PRONYX   3339
      CALL HESORT (FRQHZ,DAMP,AMPOPT,PHASE,QCON,SRULE)                  PRONYX   3340
C                                                                       PRONYX   3341
C *** RECALCULATE RELATIVE MAGNITUDES                                   PRONYX   3342
      DO 400 I=1,QCON                                                   PRONYX   3343
  400 RELMAG(I)=AMPOPT(I)/AMPREF                                        PRONYX   3344
C                                                                       PRONYX   3345
C *** WRITE RESULTS                                                     PRONYX   3346
C                                                                       PRONYX   3347
      LINES1=1                                                          PRONYX   3348
      WRITE (LUOUT,9710)                                                PRONYX   3349
      WRITE (LURES,9710)                                                PRONYX   3350
      DO 410 I=1,QCON                                                   PRONYX   3351
        FAC=DSQRT((2.0*PI*FRQHZ(I))**2 + DAMP(I)**2)                    PRONYX   3352
        IF (FAC.GT.0.0) FAC=1./FAC                                      PRONYX   3353
        PCDAMP=-DAMP(I)*FAC                                             PRONYX   3354
C        WRITE (LUOUT,9714) I,DAMP(I),FRQHZ(I),RELMAG(I),               PRONYX   3355
C     +  PHASE(I),AMPOPT(I),PCDAMP                                      PRONYX   3356
C        WRITE (LURES,9714) I,DAMP(I),FRQHZ(I),RELMAG(I),               PRONYX   3357
C     +  PHASE(I),AMPOPT(I),PCDAMP                                      PRONYX   3358
C! FLF                                                                  PRONYX   3359
        RPART=DAMP(I)                                                   PRONYX   3360
        QPART=FRQHZ(I)*(2.*PI)                                          PRONYX   3361
        WRITE (LUOUT,9714) I,RPART,QPART,FRQHZ(I),PCDAMP,AMPOPT(I),     PRONYX   3362
     +  PHASE(I)                                                        PRONYX   3363
        WRITE (LURES,9714) I,RPART,QPART,FRQHZ(I),PCDAMP,AMPOPT(I),     PRONYX   3364
     +  PHASE(I)                                                        PRONYX   3365
C! FLF                                                                  PRONYX   3366
        LINES1=LINES1+1                                                 PRONYX   3367
        IF(LINES1.GE.22) THEN                                           PRONYX   3368
           LINES1=0                                                     PRONYX   3369
           IF(CFLRUN) GOTO 420                                          PRONYX   3370
           IF(QQ(LUIN,LUOUT).EQ.'Q') GOTO 420                           PRONYX   3371
        ENDIF                                                           PRONYX   3372
C! FLF                                                                  PRONYX   3373
 410  CONTINUE                                                          PRONYX   3374
      WRITE(LUOUT,9602)SNR                                              PRONYX   3375
        LINES1=LINES1+3                                                 PRONYX   3376
 420  CONTINUE                                                          PRONYX   3377
C                                                                       PRONYX   3378
C *** CASE TERMINATION LOGIC                                            PRONYX   3379
      MRETR=0                                                           PRONYX   3380
      GO TO 8800                                                        PRONYX   3381
C                                                                       PRONYX   3382
C *** GLOBAL TERMINATION LOGIC                                          PRONYX   3383
 8000 MRETR=99                                                          PRONYX   3384
C                                                                       PRONYX   3385
C *** FINAL DISPLAY                                                     PRONYX   3386
 8800 IF(DBP)WRITE (LURES,9998)                                         PRONYX   3387
C      WRITE (LUOUT,9998)                                               PRONYX   3388
C! FLF                                                                  PRONYX   3389
      IF(LINES1.GE.18) THEN                                             PRONYX   3390
           LINES1=0                                                     PRONYX   3391
           IF(QQ(LUIN,LUOUT).EQ.'Q') GOTO 430                           PRONYX   3392
      ENDIF                                                             PRONYX   3393
 430  CONTINUE                                                          PRONYX   3394
C! FLF                                                                  PRONYX   3395
                                                                        PRONYX   3396
C                                                                       PRONYX   3397
 9060 FORMAT (3X,'##LINEAR PREDICTION ORDER LESS THAN 1:'               PRONYX   3398
     1,' SETTING LPOCON =',I4)                                          PRONYX   3399
 9062 FORMAT (3X,'##LINEAR PREDICTION ORDER EXCEEDS NNFIT/2 =',I4,':'   PRONYX   3400
     2,' SETTING LPOCON =',I4)                                          PRONYX   3401
 9064 FORMAT (3X,'##LINEAR PREDICTION ORDER EXCEEDS NSIZE =',I4,':'     PRONYX   3402
     1,' SETTING LPOCON =',I4)                                          PRONYX   3403
 9100 FORMAT (6X,'K',12X,'XSIG(K)')                                     PRONYX   3404
 9120 FORMAT (3X,'##LINEAR PREDICTION ORDER  = LPOCON =',I4)            PRONYX   3405
 9130 FORMAT (3X,'##LINEAR PREDICTION METHOD = LPMCON =',I4)            PRONYX   3406
 9140 FORMAT (3X,'##FORWARD/BACKWARD CONTROL = FBCON  =',I4)            PRONYX   3407
 9160 FORMAT (3X,'##QCON =',I4)                                         PRONYX   3408
 9200 FORMAT (3X,'##DATA MATRIX SINGULAR VALUES:'                       PRONYX   3409
     */4X,'NUMBER',10X,'VALUE')                                         PRONYX   3410
 9210 FORMAT (3X,'##DENOMINATOR Z-ROOTS:'                               PRONYX   3411
     */5X,'MODE',8X,'REAL PART',3X,'IMAGINARY')                         PRONYX   3412
 9212 FORMAT (5X,'Z-ROOT ON OR OUTSIDE UNIT CIRCLE: ',                  PRONYX   3413
     *F6.2,' @ ',F8.4,' DEGREES')                                       PRONYX   3414
 9214 FORMAT (41X,'WARNING: VERY LARGE ROOT')                           PRONYX   3415
 9220 FORMAT (3X,'##RESTART WITH PIRCON REDUCED TO',I4//)               PRONYX   3416
 9222 FORMAT (41X,'(DELETED)')                                          PRONYX   3417
 9320 FORMAT (3X,'##MAXIMUM RANK OF THE DATA MATRIX   = LPOCON =',I4)   PRONYX   3418
 9322 FORMAT (3X,'##SETTING DATA RANK LIMIT           = PIRCON =',I4)   PRONYX   3419
 9330 FORMAT (3X,'##SINGULAR-VALUE RANGE:   SIGMA MAX = ',E12.4         PRONYX   3420
     */29X,'SIGMA MIN = ',E12.4)                                        PRONYX   3421
 9450 FORMAT (3X,'##AFTER DLEASTSVD:   LPOCON,PIRCON  =',2I4)           PRONYX   3422
 9451 FORMAT (3X,'##DONOMINATOR Z-COEFFICIENTS:')                       PRONYX   3423
 9452 FORMAT (3X,'##AFTER ROOTSIFT :   LPOCON,DMODES  =',2I4            PRONYX   3424
     1/24X,'DMREAL,DMCMPX  =',2I4)                                      PRONYX   3425
 9460 FORMAT (3X,'##LPM MODE TRIMMING: TRIM LEVEL     = ',D12.6)        PRONYX   3426
 9470 FORMAT (3X,'##NUMBER OF RETAINED MODES          =',I4)            PRONYX   3427
 9600 FORMAT ('0',2X,'##SIGNAL-TO-NOISE RATIO (SNR)   = ',F10.4         PRONYX   3428
     1,6X,'FOR QCON = ',I4)                                             PRONYX   3429
 9602 FORMAT (/10X,'SIGNAL-TO-NOISE RATIO (SNR)   = ',F10.4,/)          PRONYX   3430
 9618 FORMAT (3X,'##MODAL AMPLITUDES:'                                  PRONYX   3431
     */5X,'MODE',8X,'AMPLITUDE')                                        PRONYX   3432
 9620 FORMAT (3X,'##MODAL WEIGHTS (RELATIVE):'                          PRONYX   3433
     */5X,'MODE',8X,'REL. AMP.')                                        PRONYX   3434
 9622 FORMAT (3X,I5,5X,6D17.8)                                          PRONYX   3435
C 9710 FORMAT (/2X,'MODE',5X,'DAMPING',5X,'FRQ (HZ)',3X,'REL WT',3X     PRONYX   3436
C     1,'PHASE',3X,'AMPLITUDE',2X,'%DAMPING'/)                          PRONYX   3437
 9710 FORMAT(/,T12,'REAL',T24,'IMAGINARY',T48,'DAMPING',                PRONYX   3438
     1/2X,'MODE',T12,'PART',T26,'PART',T34,'FREQ(HZ)',T49,'RATIO',      PRONYX   3439
     1T61,'AMPLITUDE',T72,'PHASE')                                      PRONYX   3440
C 9714 FORMAT (2X,I3,2X,2F12.7,F8.3,F9.2,F12.5,F8.3)                    PRONYX   3441
 9714 FORMAT (2X,I3,2X,3F12.7,2F12.5,F12.2)                             PRONYX   3442
 9910 FORMAT (3X,'##IN PRSPAK: NNFIT,NINIT,QCON,LPOCON,LPMCON,FBCON     PRONYX   3443
     *,PIRCON,MRETR = ',10I4)                                           PRONYX   3444
 9912 FORMAT (4X,'TSTEP         = ',F8.6)                               PRONYX   3445
 9914 FORMAT (4X,'SIGNAL ENERGY = ',E12.6)                              PRONYX   3446
C                                                                       PRONYX   3447
 9920 FORMAT (3X,'##IN PRSPAK: NULL SIGNAL:   ENERGY = ',E12.4)         PRONYX   3448
 9940 FORMAT (3X,'##IN PRSPAK: NO DATA        NNFIT  =',I5)             PRONYX   3449
 9998 FORMAT (3X,'##PRSPAK DONE'//)                                     PRONYX   3450
C                                                                       PRONYX   3451
      END                                                               PRONYX   3452
      SUBROUTINE QPOLRT (ADDP,NDEG,ZEROR$,ZEROI$,IFAIL)                 PRONYX   3453
*IF VAX                                                                 PRONYX   3454
*ENDIF                                                                  PRONYX   4237
      RETURN                                                            PRONYX   4238
      END                                                               PRONYX   4239
