SUBROUTINE RAT2D(U,X,Y,NXP,NYP,NX,NY,PP,A,IERR) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Computes for given functional values U(I,J) at points X(I) and C Y(J) the coefficients A(I,J,K,L) of a birational spline C interpolation. The boundary points are determined through a C simple difference scheme. C C INPUT C U Data matrix (NXP,NYP). C X,Y Coordinates of the data. C NXP,NYP Physical length of the above. C NX,NY Length of the above to be processed. C PP Parameter of the birational spline interpolation. C PP > -1. C PP = 0: Bicubic Spline interpolation. C PP --> + Infinity: Towards bilinear interpolation. C IERR Flag for matrix dimensions C = 1: wrong matrix dimension C C OUTPUT C A Coefficient matrix (NX,NY,4,4) of the birational C spline interpolation. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER (EPS=1.E-10,NNX=401,NNY=401,NB=4) C REAL U(NXP,NYP),X(NXP),Y(NYP),A(NXP,NYP,4,4) C ! REAL P(NNX,NNY),Q(NNX,NNY),R(NNX,NNY) ! REAL AX(NNX),BX(NNX),CX(NNX),RX(NNX) ! REAL AY(NNY),BY(NNY),CY(NNY),RY(NNY) ! REAL DX(NNX),DY(NNY) REAL P(NXP,NYP),Q(NXP,NYP),R(NXP,NYP) REAL AX(NXP),BX(NXP),CX(NXP),RX(NXP) REAL AY(NYP),BY(NYP),CY(NYP),RY(NYP) REAL DX(NXP),DY(NYP) C REAL B(NB,NB),C(NB,NB),D(NB,NB),E(NB,NB) !sjl: alternate matrix multiply approach: real Bt(NX-1,NB,NB) real Et(NB,NB) real Cj(NX-1,NB,NB) real Dj(NX-1,NB,NB) integer, parameter :: BLKSZ=192 ! pretty effective on westmere ierrcnt1 = 0 C C...PARAMETERS AND CHECK OF ARRAY BOUNDS C ! MMX=NNX ! MMY=NNY MMX=NXP MMY=NYP C IF ((NX.LT.2).OR.(NY.LT.2)) THEN PRINT*,' RAT2D: NX, NY SMALLER THAN TWO, FULL STOP' PRINT*,' NX, NY: ',NX,NY IERR=1 STOP ENDIF C C...PARAMETERS AND NONEQUIDISTANT GRID C NX1=NX-1 NX2=NX-2 NY1=NY-1 NY2=NY-2 II=MOD(NX1,4) DO 10 I=1,II DX(I)=1./(X(I+1)-X(I)) 10 CONTINUE DO 11 I=1+II,NX1,4 DX(I)=1./(X(I+1)-X(I)) DX(I+1)=1./(X(I+2)-X(I+1)) DX(I+2)=1./(X(I+3)-X(I+2)) DX(I+3)=1./(X(I+4)-X(I+3)) 11 CONTINUE JJ=MOD(NY1,4) DO 12 J=1,JJ DY(J)=1./(Y(J+1)-Y(J)) 12 CONTINUE DO 13 J=1+JJ,NY1,4 DY(J)=1./(Y(J+1)-Y(J)) DY(J+1)=1./(Y(J+2)-Y(J+1)) DY(J+2)=1./(Y(J+3)-Y(J+2)) DY(J+3)=1./(Y(J+4)-Y(J+3)) 13 CONTINUE P1=PP+1. P2=PP+2. P3=PP+3. DK=1./P1 B(1,1)=P2*DK B(2,1)=-DK B(4,3)=-DK B(4,1)=DK B(2,3)=B(1,1) B(1,3)=-DK B(3,1)=-DK B(3,3)=DK E(1,1)=B(1,1) E(2,1)=B(2,1) E(3,1)=B(3,1) E(4,1)=B(4,1) E(1,3)=B(1,3) E(2,3)=B(2,3) E(3,3)=B(3,3) E(4,3)=B(4,3) C C...COMPUTE PARTIAL DERIVATIVES AT THE BOUNDARIES C JJ=MOD(NY,4) DO 49 J=1,JJ P(1,J)=(U(2,J)-U(1,J))*DX(1) P(NX,J)=(U(NX,J)-U(NX1,J))*DX(NX1) 49 CONTINUE TEMP=DX(1) DO 50 J=1+JJ,NY,4 P(1,J)=(U(2,J)-U(1,J))*TEMP P(1,J+1)=(U(2,J+1)-U(1,J+1))*TEMP P(1,J+2)=(U(2,J+2)-U(1,J+2))*TEMP P(1,J+3)=(U(2,J+3)-U(1,J+3))*TEMP 50 CONTINUE TEMP=DX(NX1) DO 51 J=1+JJ,NY,4 P(NX,J)=(U(NX,J)-U(NX1,J))*TEMP P(NX,J+1)=(U(NX,J+1)-U(NX1,J+1))*TEMP P(NX,J+2)=(U(NX,J+2)-U(NX1,J+2))*TEMP P(NX,J+3)=(U(NX,J+3)-U(NX1,J+3))*TEMP 51 CONTINUE II=MOD(NX,4) DO 59 I=1,II Q(I,1)=(U(I,2)-U(I,1))*DY(1) Q(I,NY)=(U(I,NY)-U(I,NY1))*DY(NY1) 59 CONTINUE TEMP=DY(1) DO 60 I=1+II,NX,4 Q(I,1)=(U(I,2)-U(I,1))*TEMP Q(I+1,1)=(U(I+1,2)-U(I+1,1))*TEMP Q(I+2,1)=(U(I+2,2)-U(I+2,1))*TEMP Q(I+3,1)=(U(I+3,2)-U(I+3,1))*TEMP 60 CONTINUE TEMP=DY(NY1) DO 61 I=1+II,NX,4 Q(I,NY)=(U(I,NY)-U(I,NY1))*TEMP Q(I+1,NY)=(U(I+1,NY)-U(I+1,NY1))*TEMP Q(I+2,NY)=(U(I+2,NY)-U(I+2,NY1))*TEMP Q(I+3,NY)=(U(I+3,NY)-U(I+3,NY1))*TEMP 61 CONTINUE R(1,1)=.5*((P(1,2)-P(1,1))*DY(1)+ & (Q(2,1)-Q(1,1))*DX(1)) R(1,NY)=.5*((P(1,NY)-P(1,NY1))*DY(NY1)+ & (Q(2,NY)-Q(1,NY))*DX(1)) R(NX,1)=.5*((P(NX,2)-P(NX,1))*DY(1)+ & (Q(NX,1)-Q(NX1,1))*DX(NX1)) R(NX,NY)=.5*((P(NX,NY)-P(NX,NY1))*DY(NY1)+ & (Q(NX,NY)-Q(NX1,NY))*DX(NX1)) C C...PREPARE THE COEFFICIENT MATRIX FOR X C 70 IF ((NX.EQ.2).AND.(NY.EQ.2)) GOTO 130 IF (NX.EQ.2) GOTO 90 II=MOD(NX2-1,4) DO 79 I=1,II AX(I)=DX(I+1) BX(I)=P2*(DX(I+1)+DX(I)) IF (ABS(BX(I)).LT.EPS) GOTO 82 79 CONTINUE DO 80 I=1+II,NX2-1,4 AX(I)=DX(I+1) AX(I+1)=DX(I+2) AX(I+2)=DX(I+3) AX(I+3)=DX(I+4) 80 CONTINUE DO 81 I=1+II,NX2-1,4 BX(I)=P2*(DX(I+1)+DX(I)) BX(I+1)=P2*(DX(I+2)+DX(I+1)) BX(I+2)=P2*(DX(I+3)+DX(I+2)) BX(I+3)=P2*(DX(I+4)+DX(I+3)) IF ((ABS(BX(I)).LT.EPS).OR.(ABS(BX(I+1)).LT.EPS).OR. & (ABS(BX(I+2)).LT.EPS).OR.(ABS(BX(I+3)).LT.EPS)) GOTO 82 81 CONTINUE GOTO 83 82 PRINT*,' RAT2D: BX(',I,') TOO SMALL, FULL STOP' STOP 83 CONTINUE BX(NX2)=P2*(DX(NX2+1)+DX(NX2)) C C...LU- DECOMPOSITION FOR X C DO 84 K=2,NX2 KM1=K-1 H=AX(KM1)/BX(KM1) CX(KM1)=H BX(K)=BX(K)-H*AX(KM1) 84 CONTINUE IVJ=0 IPU=1 CALL RATPE(U,NXP,NYP,P,RX,MMX,MMY,IVJ,IPU,P3,AX,BX,CX,DX,NX,NY) 90 IF (NY.EQ.2) GOTO 110 C C...PREPARE THE COEFFICIENT MATRIX FOR Y C JJ=MOD(NY2-1,4) DO 99 J=1,JJ AY(J)=DY(J+1) BY(J)=P2*(DY(J+1)+DY(J)) IF (ABS(BY(J)).LT.EPS) GOTO 102 99 CONTINUE DO 100 J=1+JJ,NY2-1,4 AY(J)=DX(J+1) AY(J+1)=DY(J+2) AY(J+2)=DY(J+3) AY(J+3)=DY(J+4) 100 CONTINUE DO 101 J=1+JJ,NY2-1,4 BY(J)=P2*(DY(J+1)+DY(J)) BY(J+1)=P2*(DY(J+2)+DY(J+1)) BY(J+2)=P2*(DY(J+3)+DY(J+2)) BY(J+3)=P2*(DY(J+4)+DY(J+3)) IF ((ABS(BY(J)).LT.EPS).OR.(ABS(BY(J+1)).LT.EPS).OR. & (ABS(BY(J+2)).LT.EPS).OR.(ABS(BY(J+3)).LT.EPS)) GOTO 102 101 CONTINUE GOTO 103 102 PRINT*,' RAT2D: BY(',J,') TOO SMALL, FULL STOP' STOP 103 CONTINUE BY(NY2)=P2*(DY(NY2+1)+DY(NY2)) C C...LU- DECOMPOSITION FOR Y C DO 104 K=2,NY2 KM1=K-1 H=AY(KM1)/BY(KM1) CY(KM1)=H BY(K)=BY(K)-H*AY(KM1) 104 CONTINUE IVJ=1 IPU=1 CALL RATPE(U,NXP,NYP,Q,RY,MMX,MMY,IVJ,IPU,P3,AY,BY,CY,DY,NY,NX) 110 IF (NX.EQ.2) GOTO 120 IVJ=0 IPU=NY1 CALL RATPE(Q,MMX,MMY,R,RX,MMX,MMY,IVJ,IPU,P3,AX,BX,CX,DX,NX,NY) 120 IF (NY.EQ.2) GOTO 130 IVJ=1 IPU=1 CALL RATPE(P,MMX,MMY,R,RY,MMX,MMY,IVJ,IPU,P3,AY,BY,CY,DY,NY,NX) C C...COMPUTE THE COEFFICIENT MATRIX C ! sjl: We are computing for A(i,j,k1,k2), so make i the inner loop 130 do i=1,nx1 ! calc B at each i B(k1,k2) (also transpose) FA=(X(I+1)-X(I))/(PP+3.) Bt(i,1,1) = B(1,1) Bt(i,2,1) = B(1,1)*FA Bt(i,3,1) = B(1,3) Bt(i,1,2) = B(2,1) Bt(i,2,2) = B(2,1)*FA Bt(i,4,1) = -Bt(i,2,2) Bt(i,3,2) = B(2,3) Bt(i,4,2) = Bt(i,2,2)*(PP+2.) Bt(i,1,3) = B(3,1) Bt(i,2,3) = -Bt(i,2,1) Bt(i,3,3) = B(3,3) Bt(i,4,3) = Bt(i,2,2) Bt(i,1,4) = B(4,1) Bt(i,2,4) = -Bt(i,2,2) Bt(i,3,4) = B(4,3) Bt(i,4,4) = -Bt(i,4,2) end do do j=1,ny1 ! calc E at this j, and transpose it for better cache use FB=(Y(J+1)-Y(J))/(PP+3.) Et(1,1) = E(1,1) Et(1,2) = E(2,1) Et(1,3) = E(3,1) Et(1,4) = E(4,1) Et(2,1) = E(1,1)*FB Et(2,2) = E(2,1)*FB Et(2,3) = -Et(2,1) Et(2,4) = -Et(2,2) Et(3,1) = E(1,3) Et(3,2) = E(2,3) Et(3,3) = E(3,3) Et(3,4) = E(4,3) Et(4,1) = -Et(2,2) Et(4,2) = Et(2,2)*(PP+2.) Et(4,3) = Et(2,2) Et(4,4) = -Et(4,2) do i=1,nx1 ! calc C at each i for current j Cj(i,1,1:3:2) = U(I,J:J+1) Cj(i,1,2:4:2) = Q(I,J:J+1) Cj(i,2,1:3:2) = P(I,J:J+1) Cj(i,2,2:4:2) = R(I,J:J+1) Cj(i,3,1:3:2) = U(I+1,J:J+1) Cj(i,3,2:4:2) = Q(I+1,J:J+1) Cj(i,4,1:3:2) = P(I+1,J:J+1) Cj(i,4,2:4:2) = R(I+1,J:J+1) end do ! do k1=1,nb do k2=1,nb do i1=1,nx1,BLKSZ Dj(i1:min(nx1,i1+BLKSZ),k2,1) = 0.0 Dj(i1:min(nx1,i1+BLKSZ),k2,2) = 0.0 Dj(i1:min(nx1,i1+BLKSZ),k2,3) = 0.0 Dj(i1:min(nx1,i1+BLKSZ),k2,4) = 0.0 do k=1,nb ! block these loops to keep i in cache do i=i1, min(nx1,i1+BLKSZ) Dj(i,k2,1) = Dj(i,k2,1) + Bt(i,k,1)*Cj(i,k,k2) Dj(i,k2,2) = Dj(i,k2,2) + Bt(i,k,2)*Cj(i,k,k2) Dj(i,k2,3) = Dj(i,k2,3) + Bt(i,k,3)*Cj(i,k,k2) Dj(i,k2,4) = Dj(i,k2,4) + Bt(i,k,4)*Cj(i,k,k2) end do end do end do end do ! end do ! do k1=1,nb do k2=1,nb do i1=1,nx1,BLKSZ A(i1:min(nx1,i1+BLKSZ),j,1,k2) = 0.0 A(i1:min(nx1,i1+BLKSZ),j,2,k2) = 0.0 A(i1:min(nx1,i1+BLKSZ),j,3,k2) = 0.0 A(i1:min(nx1,i1+BLKSZ),j,4,k2) = 0.0 do k=1,nb ! block these loops to keep i in cache do i=i1, min(nx1,i1+BLKSZ) A(i,j,1,k2) = A(i,j,1,k2) + Dj(i,k,1)*Et(k,k2) A(i,j,2,k2) = A(i,j,2,k2) + Dj(i,k,2)*Et(k,k2) A(i,j,3,k2) = A(i,j,3,k2) + Dj(i,k,3)*Et(k,k2) A(i,j,4,k2) = A(i,j,4,k2) + Dj(i,k,4)*Et(k,k2) end do end do end do end do ! end do end do ! j RETURN END