← Back to branch summary

~frederix/+junk/virtual_polyfit

« back to all changes in this revision

Viewing changes to Template/NLO/SubProcesses/polfit.f

  • Committer: Rikkert Frederix
  • Date: 2019-08-13 13:33:44 UTC
  • Revision ID: frederix@physik.uzh.ch-20190813133344-ulmd3cqo0rjqt02j
updated polfit an d pvalue subroutines to comply with more modern compilers

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/polfit.f
102
102
 
103
103
*DECK POLFIT      
104
104
      SUBROUTINE POLFIT (N, X, Y, W, MAXDEG, NDEG, EPS, R, IERR, A)
 
105
c
 
106
cRF: make double precision and replace DO-loops and arirthmic
 
107
c    IF-statements to comply more modern standards
 
108
c
105
109
C***BEGIN PROLOGUE  POLFIT
106
110
C***PURPOSE  Fit discrete data in a least squares sense by polynomials
107
111
C            in one variable.
225
229
C***  END PROLOGUE  POLFIT
226
230
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 
227
231
      DOUBLE PRECISION TEMD1,TEMD2
228
 
      DIMENSION X(*), Y(*), W(*), R(*), A(*)
 
232
      DIMENSION X(*), Y(*), W(*), R(*), A(*), YP(0)
229
233
      DIMENSION CO(4,3)
230
234
      SAVE CO
231
235
      DATA  CO(1,1), CO(2,1), CO(3,1), CO(4,1), CO(1,2), CO(2,2),
244
248
      XM = M
245
249
      ETST = EPS*EPS*XM
246
250
      IF (W(1) .LT. 0.0) GO TO 2
247
 
      DO 1 I = 1,M
 
251
      DO I = 1,M
248
252
        IF (W(I) .LE. 0.0) GO TO 30
249
 
 1      CONTINUE
 
253
      enddo
250
254
      GO TO 4
251
 
 2    DO 3 I = 1,M
252
 
 3      W(I) = 1.0
 
255
 2    DO I = 1,M
 
256
         W(I) = 1.0
 
257
      enddo
253
258
 4    IF (EPS .GE. 0.0) GO TO 8
254
259
C
255
260
C DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST FOR
272
277
      K3 = K2 + MAXDEG + 2
273
278
      K4 = K3 + M
274
279
      K5 = K4 + M
275
 
      DO 9 I = 2,K4
276
 
 9      A(I) = 0.0
 
280
      DO I = 2,K4
 
281
         A(I) = 0.0
 
282
      enddo
277
283
      W11 = 0.0
278
284
      IF (N .LT. 0) GO TO 11
279
285
C
280
286
C UNCONSTRAINED CASE
281
287
C
282
 
      DO 10 I = 1,M
 
288
      DO I = 1,M
283
289
        K4PI = K4 + I
284
290
        A(K4PI) = 1.0
285
 
 10     W11 = W11 + W(I)
 
291
        W11 = W11 + W(I)
 
292
      enddo
286
293
      GO TO 13
287
294
C
288
295
C CONSTRAINED CASE
289
296
C
290
 
 11   DO 12 I = 1,M
 
297
 11   DO I = 1,M
291
298
        K4PI = K4 + I
292
 
 12     W11 = W11 + W(I)*A(K4PI)**2
 
299
        W11 = W11 + W(I)*A(K4PI)**2
 
300
      enddo
293
301
C
294
302
C COMPUTE FIT OF DEGREE ZERO
295
303
C
296
304
 13   TEMD1 = 0.0D0
297
 
      DO 14 I = 1,M
 
305
      DO I = 1,M
298
306
        K4PI = K4 + I
299
307
        TEMD1 = TEMD1 + DBLE(W(I))*DBLE(Y(I))*DBLE(A(K4PI))
300
 
 14     CONTINUE
 
308
      enddo
301
309
      TEMD1 = TEMD1/DBLE(W11)
302
310
      A(K2+1) = TEMD1
303
311
      SIGJ = 0.0
304
 
      DO 15 I = 1,M
 
312
      DO I = 1,M
305
313
        K4PI = K4 + I
306
314
        K5PI = K5 + I
307
315
        TEMD2 = TEMD1*DBLE(A(K4PI))
308
316
        R(I) = TEMD2
309
317
        A(K5PI) = TEMD2 - DBLE(R(I))
310
 
 15     SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
 
318
        SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
 
319
      enddo
311
320
      J = 0
312
321
C
313
322
C SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES THE DEGREE SELECTION CRITERION
314
323
C
315
 
      IF (EPS) 24,26,27
 
324
      if (eps.lt.0d0) then
 
325
         goto 24
 
326
      elseif (eps.eq.0d0) then
 
327
         goto 26
 
328
      else
 
329
         goto 27
 
330
      endif
 
331
      
316
332
C
317
333
C INCREMENT DEGREE
318
334
C
329
345
C COMPUTE NEW A COEFFICIENT
330
346
C
331
347
      TEMD1 = 0.0D0
332
 
      DO 18 I = 1,M
 
348
      DO I = 1,M
333
349
        K4PI = K4 + I
334
350
        TEMD2 = A(K4PI)
335
351
        TEMD1 = TEMD1 + DBLE(X(I))*DBLE(W(I))*TEMD2*TEMD2
336
 
 18     CONTINUE
 
352
      enddo
337
353
      A(JP1) = TEMD1/DBLE(W11)
338
354
C
339
355
C EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS
340
356
C
341
357
      W1 = W11
342
358
      W11 = 0.0
343
 
      DO 19 I = 1,M
 
359
      DO I = 1,M
344
360
        K3PI = K3 + I
345
361
        K4PI = K4 + I
346
362
        TEMP = A(K3PI)
347
363
        A(K3PI) = A(K4PI)
348
364
        A(K4PI) = (X(I)-A(JP1))*A(K3PI) - A(K1PJ)*TEMP
349
 
 19     W11 = W11 + W(I)*A(K4PI)**2
 
365
        W11 = W11 + W(I)*A(K4PI)**2
 
366
      enddo
350
367
C
351
368
C GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE
352
369
C PRECISION
353
370
C
354
371
      TEMD1 = 0.0D0
355
 
      DO 20 I = 1,M
 
372
      DO I = 1,M
356
373
        K4PI = K4 + I
357
374
        K5PI = K5 + I
358
375
        TEMD2 = DBLE(W(I))*DBLE((Y(I)-R(I))-A(K5PI))*DBLE(A(K4PI))
359
 
 20     TEMD1 = TEMD1 + TEMD2
 
376
        TEMD1 = TEMD1 + TEMD2
 
377
      enddo
360
378
      TEMD1 = TEMD1/DBLE(W11)
361
379
      A(K2PJ+1) = TEMD1
362
380
C
367
385
C SIGNIFICANT BITS ARE IN  A(K5PI) .
368
386
C
369
387
      SIGJ = 0.0
370
 
      DO 21 I = 1,M
 
388
      DO I = 1,M
371
389
        K4PI = K4 + I
372
390
        K5PI = K5 + I
373
391
        TEMD2 = DBLE(R(I)) + DBLE(A(K5PI)) + TEMD1*DBLE(A(K4PI))
374
392
        R(I) = TEMD2
375
393
        A(K5PI) = TEMD2 - DBLE(R(I))
376
 
 21     SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
 
394
        SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
 
395
      enddo
377
396
C
378
397
C SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE
379
398
C MAXDEG  HAS BEEN REACHED
380
399
C
381
 
      IF (EPS) 23,26,27
 
400
      if (eps.lt.0d0) then
 
401
         goto 23
 
402
      elseif (eps.eq.0d0) then
 
403
         goto 26
 
404
      else
 
405
         goto 27
 
406
      endif
382
407
C
383
408
C COMPUTE F STATISTICS  (INPUT EPS .LT. 0.)
384
409
C
448
473
C
449
474
      IF(EPS .GE. 0.0  .OR.  NDEG .EQ. MAXDEG) GO TO 36
450
475
      NDER = 0
451
 
      DO 35 I = 1,M
 
476
      DO I = 1,M
452
477
        CALL PVALUE (NDEG,NDER,X(I),R(I),YP,A)
453
 
 35     CONTINUE
 
478
      enddo
454
479
 36   EPS = SQRT(SIG/XM)
455
480
 37   RETURN
456
481
      END
458
483
 
459
484
*DECK PVALUE
460
485
      SUBROUTINE PVALUE (L, NDER, X, YFIT, YP, A)
 
486
c
 
487
cRF: make double precision and replace DO-loops to comply with more
 
488
c    modern standards
 
489
c
461
490
C***BEGIN PROLOGUE  PVALUE
462
491
C***PURPOSE  Use the coefficients generated by POLFIT to evaluate the
463
492
C            polynomial fit of degree L, along with the first NDER of
528
557
      IF (L .GT. NORD) GO TO 11
529
558
      K4 = K3 + L + 1
530
559
      IF (NDER .LT. 1) GO TO 2
531
 
      DO 1 I = 1,NDER
532
 
 1      YP(I) = 0.0
 
560
      DO I = 1,NDER
 
561
         YP(I) = 0.0
 
562
      enddo
533
563
 2    IF (L .GE. 2) GO TO 4
534
564
      IF (L .EQ. 1) GO TO 3
535
565
C
554
584
      LM1 = L - 1
555
585
      ILO = K3 + 3
556
586
      IUP = K4 + NDP1
557
 
      DO 5 I = ILO,IUP
558
 
 5      A(I) = 0.0
 
587
      DO I = ILO,IUP
 
588
         A(I) = 0.0
 
589
      enddo
559
590
      DIF = X - A(LP1)
560
591
      KC = K2 + LP1
561
592
      A(K4P1) = A(KC)
564
595
C
565
596
C EVALUATE RECURRENCE RELATIONS FOR FUNCTION VALUE AND DERIVATIVES
566
597
C
567
 
      DO 9 I = 1,LM1
 
598
      DO I = 1,LM1
568
599
        IN = L - I
569
600
        INP1 = IN + 1
570
601
        K1I = K1 + INP1
572
603
        DIF = X - A(INP1)
573
604
        VAL = A(IC) + DIF*A(K3P1) - A(K1I)*A(K4P1)
574
605
        IF (NDO .LE. 0) GO TO 8
575
 
        DO 6 N = 1,NDO
 
606
        DO N = 1,NDO
576
607
          K3PN = K3P1 + N
577
608
          K4PN = K4P1 + N
578
 
 6        YP(N) = DIF*A(K3PN) + N*A(K3PN-1) - A(K1I)*A(K4PN)
 
609
          YP(N) = DIF*A(K3PN) + N*A(K3PN-1) - A(K1I)*A(K4PN)
 
610
       enddo
579
611
C
580
612
C SAVE VALUES NEEDED FOR NEXT EVALUATION OF RECURRENCE RELATIONS
581
613
C
582
 
        DO 7 N = 1,NDO
 
614
        DO N = 1,NDO
583
615
          K3PN = K3P1 + N
584
616
          K4PN = K4P1 + N
585
617
          A(K4PN) = A(K3PN)
586
 
 7        A(K3PN) = YP(N)
 
618
          A(K3PN) = YP(N)
 
619
        enddo
587
620
 8      A(K4P1) = A(K3P1)
588
 
 9      A(K3P1) = VAL
 
621
        A(K3P1) = VAL
 
622
      enddo
589
623
C
590
624
C NORMAL RETURN OR ABORT DUE TO ERROR
591
625
C