Commit b04085de authored by Albrecht Gebhardt's avatar Albrecht Gebhardt
Browse files

fix Fortran compiler warnings

parent 86f1668b
Package: akima
Version: 0.6-3.3
Date: 2022-03-30
Date: 2022-04-18
Title: Interpolation of Irregularly and Regularly Spaced Data
Authors@R: c(person("Hiroshi", "Akima", role = c("aut", "cph"),
comment = "Fortran code (TOMS 760, 761, 697 and 433)"),
......@@ -26,7 +26,7 @@ License: ACM | file LICENSE
Depends: R (>= 2.0.0)
Imports: sp
NeedsCompilation: yes
Packaged: 2020-05-30 09:17:28 UTC; ripley
Packaged: 2022-04-18 14:34:55 UTC; agebhard
Author: Hiroshi Akima [aut, cph] (Fortran code (TOMS 760, 761, 697 and 433)),
Albrecht Gebhardt [aut, cre, cph] (R port (interp*, bicubic*
functions), bilinear code),
......@@ -36,4 +36,4 @@ Author: Hiroshi Akima [aut, cph] (Fortran code (TOMS 760, 761, 697 and 433)),
TOMS 760, 761, 697 and 433)
License_restricts_use: yes
Repository: CRAN
Date/Publication: 2020-05-30 09:50:55 UTC
Date/Publication: 2022-04-18 14:34:55 UTC
......@@ -36,7 +36,7 @@ C DECLARATION STATEMENTS
EQUIVALENCE (UK,DX),(IMN,X2,A1,M1),(IMX,X5,A5,M5),
1 (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3)
C PRELIMINARY PROCESSING
10 L0=L
L0=L
LM1=L0-1
LM2=LM1-1
LP1=L0+1
......@@ -56,14 +56,14 @@ C MAIN DO-LOOP
DO 80 K=1,N0
UK=U(K)
C ROUTINE TO LOCATE THE DESIRED POINT
20 IF(LM2.EQ.0) GO TO 27
IF(LM2.EQ.0) GO TO 27
IF(UK.GE.X(L0)) GO TO 26
IF(UK.LT.X(1)) GO TO 25
IMN=2
IMX=L0
21 I=(IMN+IMX)/2
IF(UK.GE.X(I)) GO TO 23
22 IMX=I
IMX=I
GO TO 24
23 IMN=I+1
24 IF(IMX.GT.IMN) GO TO 21
......@@ -79,7 +79,7 @@ C CHECK IF I = IPV
IPV=I
C ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND
C TO ESTIMATE THEM IF NECESSARY
40 J=I
J=I
IF(J.EQ.1) J=2
IF(J.EQ.LP1) J=L0
X3=X(J-1)
......@@ -167,6 +167,4 @@ C ERROR EXIT
RETURN
96 ERR = 7
RETURN
99 ERR = 10
RETURN
END
......@@ -86,7 +86,7 @@ C Initialization added by ThPe
ID2=0
ID3=0
C Error check
10 IF (ND.LE.1) GO TO 90
IF (ND.LE.1) GO TO 90
IF (NI.LE.0) GO TO 91
DO 11 ID=2,ND
IF (XD(ID).LE.XD(ID-1)) GO TO 92
......@@ -95,13 +95,13 @@ C Branches off special cases.
IF (ND.LE.4) GO TO 50
C General case -- Five data points of more
C Calculates some local variables.
20 NP0=MAX(3,NP)
NP0=MAX(3,NP)
NPM1=NP0-1
RENPM1=NPM1
RENNM2=NP0*(NP0-2)
C Main calculation for the general case
C First (outermost) DO-loop with respect to the desired points
30 DO 39 II=1,NI
DO 39 II=1,NI
IF (II.EQ.1) IINTPV=-1
XII=XI(II)
C Locates the interval that includes the desired point by binary
......@@ -396,6 +396,4 @@ C Error exit
RETURN
92 ERR=3
RETURN
99 ERR=10
RETURN
END
......@@ -158,17 +158,7 @@ C WRITE (*,FMT=9040)
90 CONTINUE
C WRITE (*,FMT=9050) NXD,NYD,NIP
RETURN
* Format statements for error messages
9000 FORMAT (1X,/,'*** RGBI3P Error 1: NXD = 1 or less')
9010 FORMAT (1X,/,'*** RGBI3P Error 2: NYD = 1 or less')
9020 FORMAT (1X,/,'*** RGBI3P Error 3: Identical XD values or',
+ ' XD values out of sequence',/,' IX =',I6,', XD(IX) =',
+ E11.3)
9030 FORMAT (1X,/,'*** RGBI3P Error 4: Identical YD values or',
+ ' YD values out of sequence',/,' IY =',I6,', YD(IY) =',
+ E11.3)
9040 FORMAT (1X,/,'*** RGBI3P Error 5: NIP = 0 or less')
9050 FORMAT (' NXD =',I5,', NYD =',I5,', NIP =',I5,/)
END
......@@ -349,19 +339,6 @@ C WRITE (*,FMT=9050)
120 CONTINUE
C WRITE (*,FMT=9060) NXD,NYD,NXI,NYI
RETURN
* Format statements for error messages
9000 FORMAT (1X,/,'*** RGSF3P Error 1: NXD = 1 or less')
9010 FORMAT (1X,/,'*** RGSF3P Error 2: NYD = 1 or less')
9020 FORMAT (1X,/,'*** RGSF3P Error 3: Identical XD values or',
+ ' XD values out of sequence',/,' IX =',I6,', XD(IX) =',
+ E11.3)
9030 FORMAT (1X,/,'*** RGSF3P Error 4: Identical YD values or',
+ ' YD values out of sequence',/,' IY =',I6,', YD(IY) =',
+ E11.3)
9040 FORMAT (1X,/,'*** RGSF3P Error 5: NXI = 0 or less')
9050 FORMAT (1X,/,'*** RGSF3P Error 6: NYI = 0 or less')
9060 FORMAT (' NXD =',I5,', NYD =',I5,', NXI =',I5,', NYI =',I5,
+ /)
END
......@@ -443,6 +420,27 @@ C WRITE (*,FMT=9060) NXD,NYD,NXI,NYI
+ (ZZ1-ZZ0)* (XX3-XX2)/XX1)*
+ (XX3/ (XX2-XX1)) + ZZ0
* ..
* initialize some variables to silence compiler warnings
Z00=0.0D0
Z01=0.0D0
Z02=0.0D0
Z03=0.0D0
Z10=0.0D0
Z11=0.0D0
Z12=0.0D0
Z13=0.0D0
Z20=0.0D0
Z21=0.0D0
Z22=0.0D0
Z23=0.0D0
Z30=0.0D0
Z31=0.0D0
Z32=0.0D0
Z33=0.0D0
X2=0.0D0
X3=0.0D0
Y2=0.0D0
Y3=0.0D0
* Calculation
* Initial setting of some local variables
NX0 = MAX(4,NXD)
......@@ -905,6 +903,9 @@ C WRITE (*,FMT=9060) NXD,NYD,NXI,NYI
DOUBLE PRECISION XII,YII
INTEGER IIP,IMD,IMN,IMX,IXD,IYD,NINTX,NINTY
* ..
* initialize some variables to silence compiler warnings
IXD=0
IYD=0
* DO-loop with respect to IIP, which is the point number of the
* output point
DO 30 IIP = 1,NIP
......@@ -1064,6 +1065,33 @@ C WRITE (*,FMT=9060) NXD,NYD,NXI,NYI
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* initialize some variables to silence compiler warnings
X0=0.0D0
Y0=0.0D0
Z00=0.0D0
Z11=0.0D0
ZII=0.0D0
ZX00=0.0D0
ZY00=0.0D0
ZXY00=0.0D0
P00=0.0D0
P01=0.0D0
P02=0.0D0
P03=0.0D0
P10=0.0D0
P11=0.0D0
P12=0.0D0
P13=0.0D0
P20=0.0D0
P21=0.0D0
P22=0.0D0
P23=0.0D0
P30=0.0D0
P31=0.0D0
P32=0.0D0
P33=0.0D0
IYD0=0
IXD0=0
* Calculation
* Outermost DO-loop with respect to the output point
DO 10 IIP = 1,NIP
......
......@@ -180,7 +180,7 @@
* CALL SDLCTN(NDP,XD,YD,NT,IPT,NL,IPL,
* 1 NIP,XI,YI, KTLI,ITLI)
IF (LINEAR) THEN
CALL SDLIPL(NDP,XD,YD,ZD,NT,IWK(1,1),NL,IWK(1,7),
CALL SDLIPL(NDP,XD,YD,ZD,NT,IWK(1,1),
+ NIPI,XI(IIP),YI(IIP),KTLI,ITLI, ZI(IIP),
+ EXTRPI(IIP))
ELSE
......@@ -219,14 +219,6 @@ c 60 WRITE (*,FMT=9030)
c triangle removal fails
IER = 11
RETURN
* Format statement for error message
9000 FORMAT (' ',/,'*** SDBI3P Error 1: NDP = 9 or less',/,' MD =',
+ I5,', NDP =',I5,/)
9010 FORMAT (' ',/,'*** SDBI3P Error 2: NDP not equal to NDPPV',/,
+ ' MD =',I5,', NDP =',I5,', NDPPV =',I5,/)
9020 FORMAT (' ',/,'*** SDBI3P Error 3: NIP = 0 or less',/,' MD =',
+ I5,', NDP =',I5,', NIP =',I5,/)
9030 FORMAT (' Error detected in SDTRAN called by SDBI3P',/)
END
......@@ -428,7 +420,7 @@ c triangle removal fails
* 1 NIP,XI,YI, KTLI,ITLI)
* agebhard: add linear interpolation:
IF (LINEAR) THEN
CALL SDLIPL(NDP,XD,YD,ZD,NT,IWK(1,1),NL,IWK(1,7),
CALL SDLIPL(NDP,XD,YD,ZD,NT,IWK(1,1),
+ NIPI,XI(IXI),YII,KTLI,ITLI, ZI(IXI,IYI),
+ EXTRPI(IXI,IYI))
ELSE
......@@ -472,16 +464,6 @@ c 90 WRITE (*,FMT=9040)
c triangle removal fails
IER = 11
RETURN
* Format statement for error message
9000 FORMAT (' ',/,'*** SDSF3P Error 1: NDP = 9 or less',/,' MD =',
+ I5,', NDP =',I5,/)
9010 FORMAT (' ',/,'*** SDSF3P Error 2: NDP not equal to NDPPV',/,
+ ' MD =',I5,', NDP =',I5,' NDPPV =',I5,/)
9020 FORMAT (' ',/,'*** SDSF3P Error 3: NXI = 0 or less',/,' MD =',
+ I5,', NDP =',I5,' NXI =',I5,', NYI =',I5,/)
9030 FORMAT (' ',/,'*** SDSF3P Error 4: NYI = 0 or less',/,' MD =',
+ I5,', NDP =',I5,' NXI =',I5,', NYI =',I5,/)
9040 FORMAT (' Error detected in SDTRAN called by SDSF3P',/)
END
......@@ -595,11 +577,6 @@ c WRITE (*,FMT=9010)
IERT = 6
END IF
RETURN
* Format statements
9000 FORMAT (' ',/,'*** SDTRAN Error 4: NDP outside its valid',
+ ' range',/,' NDP =',I5)
9010 FORMAT (' ',/,'*** SDTRAN Error 5: ',
+ 'Invalid data structure (LIST,LPTR,LEND)',/)
END
......@@ -2097,6 +2074,34 @@ c WRITE (*,FMT=9010)
* .. Intrinsic Functions ..
INTRINSIC MOD
* ..
* initialize some variables to silence compiler warnings
P00=0.0D0
P01=0.0D0
P02=0.0D0
P03=0.0D0
P04=0.0D0
P05=0.0D0
P10=0.0D0
P11=0.0D0
P12=0.0D0
P13=0.0D0
P14=0.0D0
P20=0.0D0
P21=0.0D0
P22=0.0D0
P23=0.0D0
P30=0.0D0
P31=0.0D0
P32=0.0D0
P40=0.0D0
P41=0.0D0
P50=0.0D0
Y0=0.0D0
X0=0.0D0
AP=0.0D0
BP=0.0D0
CP=0.0D0
DP=0.0D0
* Outermost DO-loop with respect to the output point
DO 120 IIP = 1,NIP
XII = XI(IIP)
......@@ -2653,6 +2658,9 @@ C
* entering into the least squares fit
* XK,YK,ZK = Coordinates and data value associated with K
*
* initialize some variables to silence compiler warnings
RS=0.0D0
KK = K
*
* Test for errors and initialize LMIN and LMAX.
......@@ -3065,7 +3073,7 @@ C
* agebhard: add a linear interpolator, along the lines of sdplnl
SUBROUTINE SDLIPL(NDP,XD,YD,ZD,NT,IPT,NL,IPL,NIP,XI,YI,KTLI,
SUBROUTINE SDLIPL(NDP,XD,YD,ZD,NT,IPT,NIP,XI,YI,KTLI,
+ ITLI, ZI, EXTRPI)
*
* A. Gebhardt:
......@@ -3083,10 +3091,6 @@ C
* IPT = two-dimensional integer array of dimension 3*NT
* containing the point numbers of the vertexes of
* the triangles,
* NL = number of border line segments,
* IPL = two-dimensional integer array of dimension 2*NL
* containing the point numbers of the end points of
* the border line segments,
* NIP = number of output points at which interpolation is
* to be performed,
* XI = array of dimension NIP containing the x
......@@ -3119,7 +3123,7 @@ C
*
* Specification statements
* .. Scalar Arguments ..
INTEGER NDP,NIP,NL,NT
INTEGER NDP,NIP,NT
* ..
* .. Array Arguments ..
DOUBLE PRECISION XD(NDP),XI(NIP),YD(NDP),
......@@ -3132,8 +3136,7 @@ C
INTEGER I,IDP,IIP,ITLII,ITLIPV,KTLII,KTLIPV
* ..
* .. Local Arrays ..
DOUBLE PRECISION X(3),Y(3),Z(3),
+ ZUV(3),ZV(3),ZVV(3)
DOUBLE PRECISION X(3),Y(3),Z(3)
* ..
* .. Intrinsic Functions ..
INTRINSIC DABS,MOD
......
......@@ -1582,10 +1582,7 @@ C
C Error flag returned by OPTIM.
C
25 IER = 5
c WRITE (6,100) NIT, IERR
RETURN
100 FORMAT (//5X,'*** Error in OPTIM: NIT = ',I4,
. ', IER = ',I1,' ***'/)
END
SUBROUTINE EDGE (IN1,IN2,X,Y, LWK,IWK,LIST,LPTR,
. LEND, IER)
......@@ -2078,18 +2075,11 @@ C Invalid triangulation data structure or collinear nodes
C on convex hull boundary.
C
33 IER = 3
c WRITE (6,130) IN1, IN2
130 FORMAT (//5X,'*** Error in EDGE: Invalid triangula',
. 'tion or null triangles on boundary'/
. 9X,'IN1 =',I4,', IN2=',I4/)
RETURN
C
C Error flag returned by OPTIM.
C
34 IER = 4
c WRITE (6,140) NIT, IERR
140 FORMAT (//5X,'*** Error in OPTIM: NIT = ',I4,
. ', IER = ',I1,' ***'/)
RETURN
END
SUBROUTINE GETNP (NCC,LCC,N,X,Y,LIST,LPTR,LEND,
......@@ -2202,7 +2192,11 @@ C
LOGICAL ISW, VIS, NCF, NJF, SKIP, SKSAV, LFT1, LFT2,
. LFT12
DOUBLE PRECISION DC, DL, X1, XC, XJ, XK, Y1, YC, YJ, YK
C
C initialize some vars
NL=0
ILAST=0
C Store parameters in local variables and test for errors.
C LCC1 indexes the first constraint node.
C
......@@ -4518,180 +4512,6 @@ C
13 NT = 0
IER = 2
RETURN
END
SUBROUTINE TRLPRT (NCC,LCT,N,X,Y,NROW,NT,LTRI,LOUT,
. PRNTX)
INTEGER NCC, LCT(*), N, NROW, NT, LTRI(NROW,NT), LOUT
LOGICAL PRNTX
DOUBLE PRECISION X(N), Y(N)
C
C***********************************************************
C
C From TRLPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C (817) 565-2767
C 08/22/91
C
C Given a triangulation of a set of points in the plane,
C this subroutine prints the triangle list created by
C subroutine TRLIST and, optionally, the nodal coordinates
C on logical unit LOUT. The numbers of boundary nodes,
C triangles, and arcs, and the constraint region triangle
C indexes, if any, are also printed.
C
C All parameters other than LOUT and PRNTX should be
C unaltered from their values on output from TRLIST.
C
C
C On input:
C
C NCC = Number of constraints.
C
C LCT = List of constraint triangle starting indexes
C (or dummy array of length 1 if NCC = 0).
C
C N = Number of nodes in the triangulation.
C 3 .LE. N .LE. 9999.
C
C X,Y = Arrays of length N containing the coordinates
C of the nodes in the triangulation -- not used
C unless PRNTX = TRUE.
C
C NROW = Number of rows (entries per triangle) re-
C served for the triangle list LTRI. The value
C must be 6 if only the vertex indexes and
C neighboring triangle indexes are stored, or 9
C if arc indexes are also stored.
C
C NT = Number of triangles in the triangulation.
C 1 .LE. NT .LE. 9999.
C
C LTRI = NROW by NT array whose J-th column contains
C the vertex nodal indexes (first three rows),
C neighboring triangle indexes (second three
C rows), and, if NROW = 9, arc indexes (last
C three rows) associated with triangle J for
C J = 1,...,NT.
C
C LOUT = Logical unit number for output. 0 .LE. LOUT
C .LE. 99. Output is printed on unit 6 if LOUT
C is outside its valid range on input.
C
C PRNTX = Logical variable with value TRUE if and only
C if X and Y are to be printed (to 6 decimal
C places).
C
C None of the parameters are altered by this routine.
C
C Modules required by TRLPRT: None
C
C***********************************************************
C
INTEGER I, K, LUN, NA, NB, NL, NLMAX, NMAX
DATA NMAX/9999/, NLMAX/60/
C
C Local parameters:
C
C I = DO-loop, nodal index, and row index for LTRI
C K = DO-loop and triangle index
C LUN = Logical unit number for output
C NA = Number of triangulation arcs
C NB = Number of boundary nodes
C NL = Number of lines printed on the current page
C NLMAX = Maximum number of print lines per page
C NMAX = Maximum value of N and NT (4-digit format)
C
LUN = LOUT
IF (LUN .LT. 0 .OR. LUN .GT. 99) LUN = 6
C
C Print a heading and test for invalid input.
C
c WRITE (LUN,100)
NL = 1
IF (N .LT. 3 .OR. N .GT. NMAX .OR.
. (NROW .NE. 6 .AND. NROW .NE. 9) .OR.
. NT .LT. 1 .OR. NT .GT. NMAX) THEN
C
C Print an error message and bypass the loops.
C
c WRITE (LUN,110) N, NROW, NT
GO TO 3
ENDIF
IF (PRNTX) THEN
C
C Print X and Y.
C
c WRITE (LUN,101)
NL = 6
DO 1 I = 1,N
IF (NL .GE. NLMAX) THEN
c WRITE (LUN,106)
NL = 0
ENDIF
c WRITE (LUN,102) I, X(I), Y(I)
NL = NL + 1
1 CONTINUE
ENDIF
C
C Print the triangulation LTRI.
C
IF (NL .GT. NLMAX/2) THEN
c WRITE (LUN,106)
NL = 0
ENDIF
IF (NROW .EQ. 6) THEN
c WRITE (LUN,103)
ELSE
c WRITE (LUN,104)
ENDIF
NL = NL + 5
DO 2 K = 1,NT
IF (NL .GE. NLMAX) THEN
c WRITE (LUN,106)
NL = 0
ENDIF
c WRITE (LUN,105) K, (LTRI(I,K), I = 1,NROW)
NL = NL + 1
2 CONTINUE
C
C Print NB, NA, and NT (boundary nodes, arcs, and
C triangles).
C
NB = 2*N - NT - 2
NA = NT + N - 1
c IF (NL .GT. NLMAX-6) WRITE (LUN,106)
c WRITE (LUN,107) NB, NA, NT
C
C Print NCC and LCT.
C
3 CONTINUE
c WRITE (LUN,108) NCC
c IF (NCC .GT. 0) WRITE (LUN,109) (LCT(I), I = 1,NCC)
RETURN
C
C Print formats:
C
100 FORMAT ('1',24X,'TRIPACK (TRLIST) OUTPUT')
101 FORMAT (//16X,'NODE',7X,'X(NODE)',10X,'Y(NODE)'//)
102 FORMAT (16X,I4,2E17.6)
103 FORMAT (//1X,'TRIANGLE',8X,'VERTICES',12X,'NEIGHBORS'/
. 4X,'KT',7X,'N1',5X,'N2',5X,'N3',4X,'KT1',4X,
. 'KT2',4X,'KT3'/)
104 FORMAT (//1X,'TRIANGLE',8X,'VERTICES',12X,'NEIGHBORS',
. 14X,'ARCS'/
. 4X,'KT',7X,'N1',5X,'N2',5X,'N3',4X,'KT1',4X,
. 'KT2',4X,'KT3',4X,'KA1',4X,'KA2',4X,'KA3'/)
105 FORMAT (2X,I4,2X,6(3X,I4),3(2X,I5))
106 FORMAT ('1')
107 FORMAT (/1X,'NB = ',I4,' BOUNDARY NODES',5X,
. 'NA = ',I5,' ARCS',5X,'NT = ',I5,
. ' TRIANGLES')
108 FORMAT (/1X,'NCC =',I3,' CONSTRAINT CURVES')
109 FORMAT (1X,9X,14I5)
110 FORMAT (//1X,10X,'*** INVALID PARAMETER: N =',I5,
. ', NROW =',I5,', NT =',I5,' ***')
END
SUBROUTINE TRMESH (N,X,Y, LIST,LPTR,LEND,LNEW,IER)
INTEGER N, LIST(*), LPTR(*), LEND(N), LNEW, IER
......@@ -5397,194 +5217,3 @@ C
12 IER = K
RETURN
END
SUBROUTINE TRPRNT (NCC,LCC,N,X,Y,LIST,LPTR,LEND,LOUT,
. PRNTX)
INTEGER NCC, LCC(*), N, LIST(*), LPTR(*), LEND(N),
. LOUT
LOGICAL PRNTX
DOUBLE PRECISION X(N), Y(N)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C (817) 565-2767
C 08/22/91
C
C Given a triangulation of a set of points in the plane,
C this subroutine prints the adjacency lists and, option-
C ally, the nodal coordinates on logical unit LOUT. The
C list of neighbors of a boundary node is followed by index
C 0. The numbers of boundary nodes, triangles, and arcs,
C and the constraint curve starting indexes, if any, are
C also printed.
C
C
C On input:
C
C NCC = Number of constraints.
C
C LCC = List of constraint curve starting indexes (or
C dummy array of length 1 if NCC = 0).
C
C N = Number of nodes in the triangulation.
C 3 .LE. N .LE. 9999.
C
C X,Y = Arrays of length N containing the coordinates
C of the nodes in the triangulation -- not used
C unless PRNTX = TRUE.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to subroutine
C TRMESH.
C
C LOUT = Logical unit number for output. 0 .LE. LOUT
C .LE. 99. Output is printed on unit 6 if LOUT
C is outside its valid range on input.