← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to stdred/spec/libsrc/prd.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
stdred/spec/libsrc/prd.for
 
1
C @(#)prd.for   19.1 (ES0-DMD) 02/25/03 14:29:02
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Corresponding concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C
 
29
CC
 
30
CC  *RANGS* : SORTS A VECTOR OF LENGTH 'NCAS'.
 
31
CC
 
32
      SUBROUTINE RANGS(AM,LAME)
 
33
 
 
34
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
35
C      DIMENSION JLV(30),JRV(30),AM(LAME)
 
36
 
 
37
      INTEGER J,JNC,JNDL,JR,JSS,JTWE,LAME
 
38
      REAL AMM,XX
 
39
 
 
40
      INTEGER JLV(30),JRV(30)
 
41
      REAL AM(LAME)
 
42
 
 
43
      JSS=1
 
44
      JLV(1)=1
 
45
      JRV(1)=LAME
 
46
  10  JNDL=JLV(JSS)
 
47
      JR=JRV(JSS)
 
48
      JSS=JSS-1
 
49
  20  JNC=JNDL
 
50
      J=JR
 
51
      JTWE=(JNDL+JR)/2
 
52
      XX=AM(JTWE)
 
53
  30  IF(AM(JNC).GE.XX)GOTO 40
 
54
      JNC=JNC+1
 
55
      GOTO 30
 
56
  40  IF(XX.GE.AM(J))GOTO 50
 
57
      J=J-1
 
58
      GOTO 40
 
59
  50  IF(JNC.GT.J)GOTO 60
 
60
      AMM=AM(JNC)
 
61
      AM(JNC)=AM(J)
 
62
      AM(J)=AMM
 
63
      JNC=JNC+1
 
64
      J=J-1
 
65
  60  IF (JNC.LE.J) GOTO 30
 
66
      IF ((J-JNDL).LT.(JR-JNC)) GOTO 80
 
67
      IF (JNDL.GE.J) GOTO 70
 
68
      JSS=JSS+1
 
69
      JLV(JSS)=JNDL
 
70
      JRV(JSS)=J
 
71
  70  JNDL=JNC
 
72
      GOTO 100
 
73
  80  IF (JNC.GE.JR) GOTO 90
 
74
      JSS=JSS+1
 
75
      JLV(JSS)=JNC
 
76
      JRV(JSS)=JR
 
77
  90  JR=J
 
78
  100 IF (JNDL.LT.JR) GOTO 20
 
79
      IF (JSS.NE.0) GOTO 10
 
80
      RETURN
 
81
      END
 
82
CC
 
83
CC  *RDUAL* : CALCULATES THE RESIDUALS OF ALL CASES
 
84
CC
 
85
       SUBROUTINE RDUAL(AA,JKEUS,JDA,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,
 
86
     1 JREG,X,Y,RESDU,WEIGHTS,XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL,
 
87
     1 NDXX,NDXY,JHEAD,LAB)
 
88
 
 
89
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
90
 
 
91
C       DIMENSION AA(JDA),X(NDXX,NDXY),Y(NDXY),RESDU(NDXY),WEIGHTS(NDXY)
 
92
C       DIMENSION XMED(NDXX),XMAD(NDXX),AW(NDXY),NMVAL(NDXY)
 
93
       INTEGER J,JCST,JDA,JKEUS,JNC,JPL2,JPLT,JPRT,JREG
 
94
       INTEGER NCAS,NDXX,NDXY,NVAD,NVAR,NZWE,LUB
 
95
       REAL AL,AVW,BBB,PREC,YJ
 
96
 
 
97
       REAL AA(JDA),X(NDXX,NDXY),Y(NDXY),RESDU(NDXY),WEIGHTS(NDXY)
 
98
       REAL XMED(NDXX),XMAD(NDXX),AW(NDXY)
 
99
       INTEGER NMVAL(NDXY)
 
100
 
 
101
 
 
102
       CHARACTER*10 LAB(NDXX)
 
103
       CHARACTER*60 JHEAD
 
104
       JPL2=1
 
105
       AL=NCAS
 
106
       IF(JPRT.EQ.0) GOTO 10
 
107
       IF (NVAR.EQ.1.AND.JCST.EQ.1) THEN
 
108
          IF (JKEUS.NE.2) WRITE(LUB,8000) LAB(NVAD)
 
109
          IF (JKEUS.EQ.2) WRITE(LUB,8010) LAB(NVAD)
 
110
                                    ELSE
 
111
          IF (JKEUS.NE.2) WRITE(LUB,8020) LAB(NVAD),LAB(NVAD)
 
112
          IF (JKEUS.EQ.2) WRITE(LUB,8030) LAB(NVAD),LAB(NVAD)
 
113
       ENDIF
 
114
 10    DO 20 JNC=1,NCAS
 
115
       IF (.NOT.(NVAR.EQ.1.AND.JCST.EQ.1)) GOTO 40
 
116
       BBB=Y(JNC)-AA(1)
 
117
       GOTO 60
 
118
 40    AW(JNC)=0.0
 
119
       DO 30 J=1,NVAR
 
120
 30    AW(JNC)=AW(JNC)+AA(J)*(X(J,JNC)*XMAD(J)+XMED(J))
 
121
       YJ=Y(JNC)*XMAD(NVAD)+XMED(NVAD)
 
122
       BBB=YJ-AW(JNC)
 
123
 60    IF (AA(NVAD).GT.PREC) GOTO 70
 
124
       IF (ABS(BBB).LT.PREC) BBB=0.0
 
125
       RESDU(JNC)=BBB
 
126
       GOTO 50
 
127
 70    RESDU(JNC)=BBB/AA(NVAD)
 
128
 50    IF (JPRT.EQ.0) GOTO 20
 
129
       IF (NVAR.EQ.1.AND.JCST.EQ.1) GOTO 80
 
130
       IF (AA(NVAD).GT.PREC.AND.JKEUS.NE.2)
 
131
     1 WRITE(LUB,8040) YJ,AW(JNC),BBB,NMVAL(JNC),RESDU(JNC)
 
132
       IF (AA(NVAD).LE.PREC)
 
133
     1 WRITE(LUB,8050) YJ,AW(JNC),BBB,NMVAL(JNC)
 
134
       IF (AA(NVAD).GT.PREC.AND.JKEUS.EQ.2)
 
135
     1 WRITE(LUB,8040) YJ,AW(JNC),BBB,NMVAL(JNC),RESDU(JNC),
 
136
     1 WEIGHTS(JNC)
 
137
       GOTO 20
 
138
 80    IF (JKEUS.NE.2) WRITE(LUB,8060) Y(JNC),BBB,JNC,RESDU(JNC)
 
139
       IF (JKEUS.EQ.2) WRITE(LUB,8060) Y(JNC),BBB,JNC,RESDU(JNC),
 
140
     1 WEIGHTS(JNC)
 
141
 20    CONTINUE
 
142
       IF(JPLT.LT.1) GOTO 90
 
143
       IF (AA(NVAD).LE.PREC) JPL2=0
 
144
       CALL GRAF(AW,RESDU,NCAS,0,JPL2,JPLT,NMVAL,NDXY,LUB,JREG,
 
145
     1 JHEAD,NDXX,NVAD,LAB)
 
146
       IF (JPLT.LT.3) GOTO 90
 
147
       JPLT=2
 
148
       CALL GRAF(AW,RESDU,NCAS,0,JPL2,JPLT,NMVAL,NDXY,LUB,JREG,
 
149
     1 JHEAD,NDXX,NVAD,LAB)
 
150
       JPLT=3
 
151
 90    CONTINUE
 
152
 8000  FORMAT(//11X,'OBSERVED',12X,'RESIDUAL',
 
153
     1 '   NO','  RES/SC'/9X,A10/)
 
154
 8010  FORMAT(//11X,'OBSERVED',12X,'RESIDUAL',
 
155
     1 '   NO','  RES/SC',' WEIGHT'/9X,A10/)
 
156
 8020  FORMAT(//11X,'OBSERVED',11X,'ESTIMATED',12X,'RESIDUAL',
 
157
     1 3X,'NO','  RES/SC'/9X,A10,10X,A10/)
 
158
 8030  FORMAT(//11X,'OBSERVED',11X,'ESTIMATED',12X,'RESIDUAL',
 
159
     1 3X,'NO','  RES/SC',' WEIGHT'/9X,A10,10X,A10/)
 
160
 8040  FORMAT(1X,F18.5,2(1X,F19.5),I5,1X,F7.2,F6.1)
 
161
 8050  FORMAT(1X,F18.5,2(1X,F19.5),I5,'   *****')
 
162
 8060  FORMAT(1X,F18.5,1X,F19.5,I5,1X,F7.2,F6.1)
 
163
       RETURN
 
164
       END
 
165
CC
 
166
CC  *SHHLF* : CALCULATES THE LMS IN THE ONE-DIMENSIONAL CASE.
 
167
CC
 
168
       SUBROUTINE SHHLF(W,NCAS,JQU,SLUTN,BSTD,PREC)
 
169
 
 
170
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
171
C       DIMENSION W(NCAS)
 
172
       INTEGER JNC,JPK,JQU,JTEL,K,LNK,NCAS
 
173
       REAL ABSW,BSTD,PREC,SLUTN,TEL,VERSC
 
174
 
 
175
       REAL W(NCAS)
 
176
 
 
177
       K=JQU-1
 
178
       VERSC=W(K+1)-W(1)
 
179
       SLUTN=(W(1)+W(K+1))/2.0
 
180
       BSTD=((W(K+1)-W(1))/2.0)
 
181
       IF (NCAS.EQ.2) RETURN
 
182
       JTEL=1
 
183
       LNK=NCAS-K
 
184
       DO 10 JNC=2,LNK
 
185
       JPK=JNC+K
 
186
       ABSW=W(JPK)-W(JNC)
 
187
       IF ( ABS(ABSW-VERSC).GT.PREC) GOTO 20
 
188
       JTEL=JTEL+1
 
189
       SLUTN=SLUTN+(W(JNC)+W(JPK))/2.0
 
190
       GOTO 10
 
191
 20    IF (ABSW.LT.VERSC) THEN
 
192
       VERSC=ABSW
 
193
       JTEL=1
 
194
       BSTD=((W(JPK)-W(JNC))/2.0)
 
195
       SLUTN=(W(JNC)+W(JPK))/2.0
 
196
       ENDIF
 
197
 10    CONTINUE
 
198
       TEL=JTEL
 
199
       SLUTN=SLUTN/TEL
 
200
       RETURN
 
201
       END
 
202
CC
 
203
CC  *TRC* : RETRANSFORMATION OF THE VARIANCE-COVARIANCE MATRIX
 
204
CC
 
205
       SUBROUTINE TRC(H,DA,NMXV,NDXX,NDXY,NVAR,JCST,NVSB,NVAD,
 
206
     1 XMED,XMAD)
 
207
 
 
208
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
209
C       DIMENSION H(NMXV,NDXX),DA(NMXV)
 
210
C       DOUBLE PRECISION H,XMP2,HNN
 
211
C       DIMENSION XMED(NDXX),XMAD(NDXX)
 
212
       INTEGER  J,JCST,JU,K,K2,NDXX,NDXY,NMXV,NVAD,NVAR,NVSB
 
213
 
 
214
       REAL DA(NMXV)
 
215
       DOUBLE PRECISION H(NMXV,NDXX),XMP2,HNN
 
216
       REAL XMED(NDXX),XMAD(NDXX)
 
217
 
 
218
       XMP2=DBLE(XMAD(NVAD))*DBLE(XMAD(NVAD))
 
219
       IF (JCST.EQ.0) THEN
 
220
          DO 10 J=1,NVAR
 
221
          DO 20 K=1,J
 
222
          H(J,K)=H(J,K)*(XMP2/(DBLE(XMAD(J))*DBLE(XMAD(K))))
 
223
 20          CONTINUE
 
224
          DA(J)=DSQRT(H(J,J))
 
225
 10          CONTINUE
 
226
       ELSE
 
227
          DO 25 J=1,NVAR
 
228
          H(J,NVAD)=H(J,J)
 
229
 25          CONTINUE
 
230
          DO 30 J=1,NVAR
 
231
          DO 40 K=1,J
 
232
          H(J,K)=H(J,K)*XMP2/(DBLE(XMAD(J))*DBLE(XMAD(K)))
 
233
 40          CONTINUE
 
234
          DA(J)=DSQRT(H(J,J))
 
235
 30          CONTINUE
 
236
          DO 50 K=1,NVSB
 
237
          H(NVAR,K)=H(K,NVAR)*XMP2/DBLE(XMAD(K))
 
238
          DO 60 K2=1,NVAR
 
239
          IF (K.EQ.K2) THEN
 
240
            H(NVAR,K)=H(NVAR,K)-DBLE(XMED(K))*XMP2/
 
241
     1                      (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K2,NVAD)
 
242
            GOTO 60
 
243
          ENDIF
 
244
          IF (K.LT.K2) THEN
 
245
               H(NVAR,K)=H(NVAR,K)-(DBLE(XMED(K2))*XMP2)/
 
246
     1                      (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K,K2)
 
247
                        ELSE
 
248
               H(NVAR,K)=H(NVAR,K)-DBLE(XMED(K2))*XMP2/
 
249
     1                      (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K2,K)
 
250
          ENDIF
 
251
 60          CONTINUE
 
252
 50          CONTINUE
 
253
          H(NVAR,NVAR)=H(NVAR,NVAD)*XMP2
 
254
          DO 70 K=1,NVAR
 
255
          H(NVAR,NVAR)=H(NVAR,NVAR)+
 
256
     1            (DBLE(XMED(K))*DBLE(XMED(K)))*XMP2/
 
257
     1                         (DBLE(XMAD(K))*DBLE(XMAD(K)))*H(K,NVAD)
 
258
 70          CONTINUE
 
259
          DO 80 K=1,NVAR
 
260
          IF (K.NE.NVAR) THEN
 
261
             H(NVAR,NVAR)=H(NVAR,NVAR)-2.0D0*XMP2*DBLE(XMED(K))/
 
262
     1                                    (DBLE(XMAD(K)))*H(K,NVAR)
 
263
                         ELSE
 
264
             H(NVAR,NVAR)=H(NVAR,NVAR)-2.0D0*XMP2*DBLE(XMED(K))/
 
265
     1                                    (DBLE(XMAD(K)))*H(NVAR,NVAD)
 
266
          ENDIF
 
267
 80          CONTINUE
 
268
          DO 90 J=1,NVSB
 
269
          JU=J+1
 
270
          DO 90 K=JU,NVAR
 
271
          HNN=2.0D0*DBLE(XMED(J))*DBLE(XMED(K))*XMP2
 
272
          H(NVAR,NVAR)=H(NVAR,NVAR)+HNN/
 
273
     1                              (DBLE(XMAD(J))*DBLE(XMAD(K)))*H(J,K)
 
274
 90          CONTINUE
 
275
          DA(NVAR)=DSQRT(H(NVAR,NVAR))
 
276
       ENDIF
 
277
       RETURN
 
278
       END
 
279
CC
 
280
CC  *SCHCV* : PRINTS THE VARIANCE-COVARIANCE MATRIX OF THE
 
281
CC            REGRESSION COEFFICIENTS.
 
282
CC
 
283
       SUBROUTINE SCHCV(NVAR,NDXX,NMXV,H,LUB)
 
284
 
 
285
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
286
       INTEGER J,K,NDXX,NMXV,NVAR,LUB
 
287
 
 
288
       DOUBLE PRECISION H(NMXV,NDXX)
 
289
       WRITE(LUB,9000)
 
290
       DO 10 J=1,NVAR
 
291
       WRITE(LUB,9010) (H(J,K),K=1,J)
 
292
 10    CONTINUE
 
293
 9000  FORMAT(//' VARIANCE - COVARIANCE MATRIX = '/)
 
294
 9010  FORMAT(5(5X,D10.4))
 
295
       RETURN
 
296
       END
 
297
CC
 
298
CC  *RTRAN* : RETRANSFORMATION OF THE REGRESSION COEFFICIENTS.
 
299
CC
 
300
       SUBROUTINE RTRAN(NVAR,JCST,NVSB,NVAD,NDXX,XMED,XMAD,AA,JAL,FCKW)
 
301
 
 
302
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
303
C       DIMENSION AA(JAL),XMED(NDXX),XMAD(NDXX)
 
304
 
 
305
       INTEGER J,JAL,JCST,NDXX,NVAD,NVAR,NVSB
 
306
       REAL  FCKW
 
307
 
 
308
       REAL AA(JAL),XMED(NDXX),XMAD(NDXX)
 
309
 
 
310
       IF (NVAR.LE.1) THEN
 
311
         AA(1)=AA(1)*XMAD(NVAD)/XMAD(1)
 
312
         GOTO 30
 
313
       ENDIF
 
314
       DO 10 J=1,NVSB
 
315
 10    AA(J)=AA(J)*XMAD(NVAD)/XMAD(J)
 
316
       IF (JCST.EQ.0) THEN
 
317
         AA(NVAR)=AA(NVAR)*XMAD(NVAD)/XMAD(NVAR)
 
318
                      ELSE
 
319
         AA(NVAR)=AA(NVAR)*XMAD(NVAD)
 
320
         DO 20 J=1,NVSB
 
321
 20         AA(NVAR)=AA(NVAR)-AA(J)*XMED(J)
 
322
         AA(NVAR)=AA(NVAR)+XMED(NVAD)
 
323
       ENDIF
 
324
 30    FCKW=FCKW*(XMAD(NVAD)*XMAD(NVAD))
 
325
       RETURN
 
326
       END
 
327
CC
 
328
CC  *RANPN* : RANDOM NUMBER GENERATOR .
 
329
CC
 
330
       SUBROUTINE RANPN(NCAS,NVAR,JNDVC,NDXX,JEY,JLP,MAXRP)
 
331
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
332
C       DIMENSION JNDVC(NDXX)
 
333
       INTEGER J,JEY,JEYQT,JLP,JNC,JNCM,JRAN,MAXRP,NCAS,NDXX,NVAR
 
334
       REAL AL,RNULE,RY
 
335
 
 
336
       INTEGER JNDVC(NDXX)
 
337
 
 
338
CC     WE PROGRAMMED THIS RANDOM NUMBER GENERATOR OURSELVES BECAUSE
 
339
CC     WE WANTED IT TO BE MACHINE INDEPENDENT. IT SHOULD RUN ON MOST
 
340
CC     COMPUTERS BECAUSE THE LARGEST INTEGER USED IS LESS THAN 2**30.
 
341
CC     THE PERIOD IS 2**16=65536, WHICH IS SUFFICIENT FOR THIS
 
342
CC     APPLICATION.
 
343
       AL=NCAS
 
344
       JLP=JLP+1
 
345
       IF (JLP.LE.MAXRP) GOTO 10
 
346
       RETURN
 
347
 10    DO 20 JNC=1,NVAR
 
348
 30    JEY=JEY*5761+999
 
349
       JEYQT=JEY/65536
 
350
       JEY=JEY-JEYQT*65536
 
351
       RY=JEY
 
352
       RNULE=RY/65536.0
 
353
       JRAN=INT(RNULE*AL)+1
 
354
       IF (JNC.EQ.1) GOTO 50
 
355
       JNCM=JNC-1
 
356
       DO 40 J=1,JNCM
 
357
       IF (JNDVC(J).EQ.JRAN) GOTO 30
 
358
 40    CONTINUE
 
359
 50    JNDVC(JNC)=JRAN
 
360
 20    CONTINUE
 
361
       RETURN
 
362
       END
 
363
CC
 
364
CC  *GENPN* : GENERATES ALL COMBINATIONS OF 'NVAR' INTEGERS
 
365
CC            BETWEEN 1 AND 'NCAS', WHERE NVAR IS AT MOST FIVE.
 
366
CC
 
367
       SUBROUTINE GENPN(NCAS,NVAR,JNDVC,NDXX,JJJ)
 
368
 
 
369
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
370
C       DIMENSION JNDVC(NDXX)
 
371
       INTEGER J,JJJ,NDXX,NCAS,NVAR
 
372
       INTEGER JNDVC(NDXX)
 
373
 
 
374
       IF (JJJ.LE.1) THEN
 
375
         DO 5 J=1,NVAR
 
376
 5         JNDVC(J)=J
 
377
           RETURN
 
378
       ENDIF
 
379
       IF (JNDVC(NVAR).EQ.NCAS) GOTO 10
 
380
       JNDVC(NVAR)=JNDVC(NVAR)+1
 
381
       RETURN
 
382
 10    IF (JNDVC(NVAR-1).EQ.NCAS-1) GOTO 20
 
383
       JNDVC(NVAR-1)=JNDVC(NVAR-1)+1
 
384
       JNDVC(NVAR)=JNDVC(NVAR-1)+1
 
385
       RETURN
 
386
 20    IF(JNDVC(NVAR-2).EQ.NCAS-2) GOTO 30
 
387
       JNDVC(NVAR-2)=JNDVC(NVAR-2)+1
 
388
       JNDVC(NVAR-1)=JNDVC(NVAR-2)+1
 
389
       JNDVC(NVAR)=JNDVC(NVAR-1)+1
 
390
       RETURN
 
391
 30    IF(JNDVC(NVAR-3).EQ.NCAS-3) GOTO 40
 
392
       JNDVC(NVAR-3)=JNDVC(NVAR-3)+1
 
393
       JNDVC(NVAR-2)=JNDVC(NVAR-3)+1
 
394
       JNDVC(NVAR-1)=JNDVC(NVAR-2)+1
 
395
       JNDVC(NVAR)=JNDVC(NVAR-1)+1
 
396
       RETURN
 
397
 40    IF(JNDVC(NVAR-4).EQ.NCAS-4) GOTO 50
 
398
       JNDVC(NVAR-4)=JNDVC(NVAR-4)+1
 
399
       JNDVC(NVAR-3)=JNDVC(NVAR-4)+1
 
400
       JNDVC(NVAR-2)=JNDVC(NVAR-3)+1
 
401
       JNDVC(NVAR-1)=JNDVC(NVAR-2)+1
 
402
       JNDVC(NVAR)=JNDVC(NVAR-1)+1
 
403
       RETURN
 
404
 50    IF(JNDVC(NVAR-5).EQ.NCAS-5) RETURN
 
405
       JNDVC(NVAR-5)=JNDVC(NVAR-5)+1
 
406
       JNDVC(NVAR-4)=JNDVC(NVAR-5)+1
 
407
       JNDVC(NVAR-3)=JNDVC(NVAR-4)+1
 
408
       JNDVC(NVAR-2)=JNDVC(NVAR-3)+1
 
409
       JNDVC(NVAR-1)=JNDVC(NVAR-2)+1
 
410
       JNDVC(NVAR)=JNDVC(NVAR-1)+1
 
411
       RETURN
 
412
       END
 
413
CC
 
414
CC  *EQUAT* : SOLVES A SYSTEM OF LINEAR EQUATIONS.
 
415
CC
 
416
      SUBROUTINE EQUAT(AM,NMXV,NDXX,HVEC,JDMB,NA,NB,NERR,DTEM)
 
417
 
 
418
C       INCLUDE 'MID_REL_INCL:implicit.inc'
 
419
C      DIMENSION HVEC(JDMB),AM(NMXV,NDXX)
 
420
C      DOUBLE PRECISION HVEC,TURN,SWAP,DTEM,DETER
 
421
      INTEGER J,JBEGC,JBEGX,JDEL,JENDC,JENDX,JDM,JDMB,JHFD
 
422
      INTEGER JMAT,JNC,JNCB,JNCC,JNCD,JNCE,JNCF,JNK,JROW
 
423
      INTEGER LCLPL,LDEL,N,NA,NB,NC,NDXX,NEQA,NERR,NMXV,NZNDE
 
424
 
 
425
      REAL AM(NMXV,NDXX)
 
426
      DOUBLE PRECISION HVEC(JDMB),TURN,SWAP,DTEM,DETER
 
427
 
 
428
      JDM=NMXV
 
429
      DETER=1.0
 
430
      N=NA
 
431
      JMAT=N+NB
 
432
      JNK=0
 
433
      DO 10 J=1,JMAT
 
434
      JNK=(J-1)*NMXV
 
435
      DO 10 NC=1,NMXV
 
436
      JNK=JNK+1
 
437
      HVEC(JNK)=AM(NC,J)
 
438
 10   CONTINUE
 
439
      NZNDE=N-1
 
440
      LCLPL=-JDM
 
441
      DO 120 JHFD=1,N
 
442
      TURN=0.
 
443
      LCLPL=LCLPL+JDM+1
 
444
      JDEL=LCLPL+N-JHFD
 
445
      DO 40 JNCB=LCLPL,JDEL
 
446
      IF(DABS(HVEC(JNCB))-DABS(TURN)) 40,40,30
 
447
 30   TURN=HVEC(JNCB)
 
448
      LDEL=JNCB
 
449
 40   CONTINUE
 
450
      IF(TURN) 50,170,50
 
451
 50   IF(LDEL-LCLPL) 60,80,60
 
452
 60   DETER=-DETER
 
453
      LDEL=LDEL-JDM
 
454
      JNCB=LCLPL-JDM
 
455
      DO 70 JNCC=JHFD,JMAT
 
456
      LDEL=LDEL+JDM
 
457
      JNCB=JNCB+JDM
 
458
      SWAP=HVEC(JNCB)
 
459
      HVEC(JNCB)=HVEC(LDEL)
 
460
 70   HVEC(LDEL)=SWAP
 
461
 80   DETER=DETER*TURN
 
462
      IF (JHFD.EQ.N)  GOTO 120
 
463
      TURN=1./TURN
 
464
      JNCB=LCLPL+1
 
465
      DO 90 JNCC=JNCB,JDEL
 
466
 90   HVEC(JNCC)=HVEC(JNCC)*TURN
 
467
      JNCD=LCLPL
 
468
      JROW=JHFD+1
 
469
      DO 110 JNCB=JROW,N
 
470
      JNCD=JNCD+1
 
471
      JNCE=LCLPL
 
472
      JNCF=JNCD
 
473
      DO 100 JNCC=JROW,JMAT
 
474
      JNCE=JNCE+JDM
 
475
      JNCF=JNCF+JDM
 
476
 100  HVEC(JNCF)=HVEC(JNCF)-HVEC(JNCE)*HVEC(JNCD)
 
477
 110  CONTINUE
 
478
 120  CONTINUE
 
479
      DTEM=DETER
 
480
      NERR=0
 
481
      NEQA=N+1
 
482
      JBEGX=NZNDE*JDM+1
 
483
      DO 150 JNC=NEQA,JMAT
 
484
      JBEGX=JBEGX+JDM
 
485
      JENDX=JBEGX+N
 
486
      JBEGC=N*JDM+1
 
487
      JENDC=JBEGC+NZNDE
 
488
      DO 140 JNCB=1,NZNDE
 
489
      JENDX=JENDX-1
 
490
      JBEGC=JBEGC-JDM
 
491
      JENDC=JENDC-JDM-1
 
492
      HVEC(JENDX)=HVEC(JENDX)/HVEC(JENDC+1)
 
493
      SWAP=HVEC(JENDX)
 
494
      JNCD=JBEGX-1
 
495
      DO 130 JNCC=JBEGC,JENDC
 
496
      JNCD=JNCD+1
 
497
 130  HVEC(JNCD)=HVEC(JNCD)-HVEC(JNCC)*SWAP
 
498
 140  CONTINUE
 
499
      HVEC(JBEGX)=HVEC(JBEGX)/HVEC(1)
 
500
 150  CONTINUE
 
501
      JNC=-JDM
 
502
      JBEGX=NZNDE*JDM+1
 
503
      JENDX=JBEGX+NZNDE
 
504
      DO 160 JNCB=NEQA,JMAT
 
505
      JBEGX=JBEGX+JDM
 
506
      JENDX=JENDX+JDM
 
507
      JNC=JNC+JDM
 
508
      JNCD=JNC
 
509
      DO 160 JNCC=JBEGX,JENDX
 
510
      JNCD=JNCD+1
 
511
 160  HVEC(JNCD)=HVEC(JNCC)
 
512
      GOTO 180
 
513
 170  NERR=-1
 
514
      DTEM=DETER
 
515
 180  JNK=0
 
516
      DO 190 J=1,JMAT
 
517
      DO 190 NC=1,NMXV
 
518
      JNK=JNK+1
 
519
      AM(NC,J)=HVEC(JNK)
 
520
 190  CONTINUE
 
521
      RETURN
 
522
      END