diff --git a/src/trans/CMakeLists.txt b/src/trans/CMakeLists.txt index 9eeae222c..f12c632ae 100644 --- a/src/trans/CMakeLists.txt +++ b/src/trans/CMakeLists.txt @@ -34,6 +34,7 @@ endif() ## Sources which are precision independent can go into a common library list( APPEND ectrans_common_src + algor/ectrans_blas_mod.F90 sharedmem/sharedmem_mod.F90 sharedmem/sharedmem.c internal/abort_trans_mod.F90 @@ -76,6 +77,7 @@ ecbuild_add_library( LINKER_LANGUAGE Fortran SOURCES ${ectrans_common_src} PUBLIC_LIBS fiat + PRIVATE_LIBS ${LAPACK_LIBRARIES} ) ectrans_target_fortran_module_directory( TARGET ectrans_common diff --git a/src/trans/algor/butterfly_alg_mod.F90 b/src/trans/algor/butterfly_alg_mod.F90 index 71bd7f4d8..364bd411e 100644 --- a/src/trans/algor/butterfly_alg_mod.F90 +++ b/src/trans/algor/butterfly_alg_mod.F90 @@ -12,8 +12,7 @@ MODULE BUTTERFLY_ALG_MOD USE PARKIND1, ONLY : JPRD, JPRM, JPIM, JPRB, JPIB USE INTERPOL_DECOMP_MOD, ONLY : COMPUTE_ID USE SHAREDMEM_MOD, ONLY : SHAREDMEM, SHAREDMEM_ASSOCIATE - -use, intrinsic :: ieee_exceptions +USE ECTRANS_BLAS_MOD, ONLY : GEMM, GEMV IMPLICIT NONE @@ -64,12 +63,6 @@ MODULE BUTTERFLY_ALG_MOD REAL(KIND=JPRB) , ALLOCATABLE :: COMMSBUF(:) ! for communicating packed bufferfly_structs END TYPE CLONE ! between MPI tasks -#ifdef WITH_IEEE_HALT -LOGICAL, PARAMETER :: LL_IEEE_HALT = .TRUE. -#else -LOGICAL, PARAMETER :: LL_IEEE_HALT = .FALSE. -#endif - LOGICAL, PARAMETER :: LLDOUBLE = (JPRB == JPRD) CONTAINS @@ -610,15 +603,9 @@ SUBROUTINE MULT_BUTV(CDTRANS,YD_STRUCT,PVECIN,PVECOUT) IFR = YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%IFROW ILR = YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%ILROW IROWS=YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%IROWS - IF (LLDOUBLE) THEN - CALL DGEMV('T',IROWS,YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%IRANK,& - & 1.0_JPRD,YNODE%B,IROWS,PVECIN(IFR:ILR),1,& - & 0.0_JPRD,ZBETA(IBTST:IBTEN,IBETALV),1) - ELSE - CALL SGEMV('T',IROWS,YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%IRANK,& - & 1.0_JPRM,YNODE%B,IROWS,PVECIN(IFR:ILR),1,& - & 0.0_JPRM,ZBETA(IBTST:IBTEN,IBETALV),1) - ENDIF + CALL GEMV('T',IROWS,YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%IRANK,& + & 1.0_JPRB,YNODE%B,IROWS,PVECIN(IFR:ILR),1,& + & 0.0_JPRB,ZBETA(IBTST:IBTEN,IBETALV),1) ENDIF ILM1 = JL-1 IBETALVM1=MOD(ILM1,2) @@ -689,15 +676,9 @@ SUBROUTINE MULT_BUTV(CDTRANS,YD_STRUCT,PVECIN,PVECOUT) IFR = YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%IFROW ILR = YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%ILROW IROWS = YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%IROWS - IF (LLDOUBLE) THEN - CALL DGEMV('N',IROWS,YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%IRANK,& - & 1.0_JPRD,YNODE%B,IROWS,ZBETA(IBTST:IBTEN,IBETALV),1,& - & 0.0_JPRD,PVECOUT(IFR:ILR),1) - ELSE - CALL SGEMV('N',IROWS,YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%IRANK,& - & 1.0_JPRM,YNODE%B,IROWS,ZBETA(IBTST:IBTEN,IBETALV),1,& - & 0.0_JPRM,PVECOUT(IFR:ILR),1) - ENDIF + CALL GEMV('N',IROWS,YD_STRUCT%SLEV(ILEVS)%NODE(JJ,JK)%IRANK,& + & 1.0_JPRB,YNODE%B,IROWS,ZBETA(IBTST:IBTEN,IBETALV),1,& + & 0.0_JPRB,PVECOUT(IFR:ILR),1) ENDIF ENDDO ENDDO @@ -724,7 +705,6 @@ SUBROUTINE MULT_BUTM(CDTRANS,YD_STRUCT,KF,PVECIN,PVECOUT,KWV) REAL(KIND=JPRB) :: ZVECIN(YD_STRUCT%N_ORDER,KF),ZVECOUT(YD_STRUCT%N_ORDER,KF) REAL(KIND=JPRB),ALLOCATABLE :: ZBETA(:,:,:) LOGICAL :: LLTRANSPOSE -LOGICAL :: LL_HALT_INVALID ! IKWV==0 only, LLTRANSPOSE = true only REAL(KIND=JPRD),ALLOCATABLE :: ZPNONIM_D(:,:) @@ -770,7 +750,7 @@ SUBROUTINE MULT_BUTM(CDTRANS,YD_STRUCT,KF,PVECIN,PVECOUT,KWV) IF( IM <=0 ) CALL ABOR1('mult_butm: IM<=0 not allowed') IF(IN>0) THEN IF (LLDOUBLE.OR.(IKWV == 0)) THEN - IF(.not.LLDOUBLE) THEN + IF(.NOT.LLDOUBLE) THEN ALLOCATE(ZPNONIM_D(IM,IN)) II=0 DO JN=1,IN @@ -780,25 +760,20 @@ SUBROUTINE MULT_BUTM(CDTRANS,YD_STRUCT,KF,PVECIN,PVECOUT,KWV) ENDDO ENDDO ZBETA_D(1:IM,1:KF)=REAL(ZBETA(IBTST:IBTST+IM-1,1:KF,IBETALV),JPRD) - CALL DGEMM('T','N',IN,KF,IM,1.0_JPRD,& - & ZPNONIM_D,IM,ZBETA_D,ILBETA,0.0_JPRD,& - & ZOUT_D,YD_STRUCT%N_ORDER) + CALL GEMM('T','N',IN,KF,IM,1.0_JPRD,& + & ZPNONIM_D(1,1),IM,ZBETA_D(1,1),ILBETA,0.0_JPRD,& + & ZOUT_D(1,1),YD_STRUCT%N_ORDER) ZVECOUT(YNODE%IRANK+1:YNODE%IRANK+IN,1:KF) = REAL(ZOUT_D(1:IN,1:KF),JPRM) DEALLOCATE(ZPNONIM_D) ELSE - CALL DGEMM('T','N',IN,KF,IM,1.0_JPRD,& + CALL GEMM('T','N',IN,KF,IM,1.0_JPRD,& & YNODE%PNONIM(1),IM,ZBETA(IBTST,1,IBETALV),ILBETA,0.0_JPRD,& & ZVECOUT(YNODE%IRANK+1,1),YD_STRUCT%N_ORDER) ENDIF ELSE - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('T','N',IN,KF,IM,1.0_JPRM,& + CALL GEMM('T','N',IN,KF,IM,1.0_JPRM,& & YNODE%PNONIM(1),IM,ZBETA(IBTST,1,IBETALV),ILBETA,0.0_JPRM,& & ZVECOUT(YNODE%IRANK+1,1),YD_STRUCT%N_ORDER) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) ENDIF ENDIF DO JF=1,KF @@ -818,12 +793,12 @@ SUBROUTINE MULT_BUTM(CDTRANS,YD_STRUCT,KF,PVECIN,PVECOUT,KWV) IROWS =YNODE%IROWS IRANK = YNODE%IRANK IF (LLDOUBLE.OR.(IKWV == 0)) THEN - IF(.not.LLDOUBLE) THEN + IF(.NOT.LLDOUBLE) THEN ALLOCATE(ZB_D(IROWS,IRANK)) ZB_D(1:IROWS,1:IRANK) = REAL(YNODE%B(1:IROWS,1:IRANK),JPRD) ZIN_D(1:ILR-IFR+1,1:KF) = REAL(PVECIN(IFR:ILR,1:KF),JPRD) - CALL DGEMM('T','N',IRANK,KF,IROWS,1.0_JPRD,& + CALL GEMM('T','N',IRANK,KF,IROWS,1.0_JPRD,& & ZB_D,IROWS,ZIN_D,IRIN,0.0_JPRD,& & ZBETA_D,ILBETA) @@ -831,19 +806,14 @@ SUBROUTINE MULT_BUTM(CDTRANS,YD_STRUCT,KF,PVECIN,PVECOUT,KWV) DEALLOCATE(ZB_D) ELSE - CALL DGEMM('T','N',IRANK,KF,IROWS,1.0_JPRD,& - & YNODE%B,IROWS,PVECIN(IFR,1),IRIN,0.0_JPRD,& + CALL GEMM('T','N',IRANK,KF,IROWS,1.0_JPRD,& + & YNODE%B(1,1),IROWS,PVECIN(IFR,1),IRIN,0.0_JPRD,& & ZBETA(IBTST,1,IBETALV),ILBETA) END IF ELSE - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('T','N',IRANK,KF,IROWS,1.0_JPRM,& - & YNODE%B,IROWS,PVECIN(IFR,1),IRIN,0.0_JPRM,& + CALL GEMM('T','N',IRANK,KF,IROWS,1.0_JPRM,& + & YNODE%B(1,1),IROWS,PVECIN(IFR,1),IRIN,0.0_JPRM,& & ZBETA(IBTST,1,IBETALV),ILBETA) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) ENDIF ENDIF ILM1 = JL-1 @@ -867,7 +837,7 @@ SUBROUTINE MULT_BUTM(CDTRANS,YD_STRUCT,KF,PVECIN,PVECOUT,KWV) IF( IM <=0 ) CALL ABOR1('mult_butm: IM<=0 not allowed') IF(IN>0) THEN IF (LLDOUBLE.OR.(IKWV == 0)) THEN - IF(.not.LLDOUBLE) THEN + IF(.NOT.LLDOUBLE) THEN ALLOCATE(ZPNONIM_D(IM,IN)) II=0 DO JN=1,IN @@ -878,26 +848,21 @@ SUBROUTINE MULT_BUTM(CDTRANS,YD_STRUCT,KF,PVECIN,PVECOUT,KWV) ENDDO ZBETA_D(1:IM,1:KF)=REAL(ZBETA(IBTST:IBTST+IM-1,1:KF,IBETALV),JPRD) - CALL DGEMM('T','N',IN,KF,IM,1.0_JPRD,& + CALL GEMM('T','N',IN,KF,IM,1.0_JPRD,& & ZPNONIM_D,IM,ZBETA_D,ILBETA,0.0_JPRD,& & ZOUT_D,YD_STRUCT%N_ORDER) ZVECOUT(YNODE%IRANK+1:YNODE%IRANK+IN,1:KF) = REAL(ZOUT_D(1:IN,1:KF),JPRM) DEALLOCATE(ZPNONIM_D) ELSE - CALL DGEMM('T','N',IN,KF,IM,1.0_JPRD,& + CALL GEMM('T','N',IN,KF,IM,1.0_JPRD,& & YNODE%PNONIM(1),IM,ZBETA(IBTST,1,IBETALV),ILBETA,0.0_JPRD,& & ZVECOUT(YNODE%IRANK+1,1),YD_STRUCT%N_ORDER) ENDIF ELSE - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('T','N',IN,KF,IM,1.0_JPRM,& + CALL GEMM('T','N',IN,KF,IM,1.0_JPRM,& & YNODE%PNONIM(1),IM,ZBETA(IBTST,1,IBETALV),ILBETA,0.0_JPRM,& & ZVECOUT(YNODE%IRANK+1,1),YD_STRUCT%N_ORDER) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) ENDIF ENDIF DO JF=1,KF @@ -963,20 +928,9 @@ SUBROUTINE MULT_BUTM(CDTRANS,YD_STRUCT,KF,PVECIN,PVECOUT,KWV) ENDDO IF( IRANK <=0 ) CALL ABOR1('mult_butm: IRANK<=0 not allowed') IF(YNODE%ICOLS > IRANK) THEN - IF (LLDOUBLE) THEN - CALL DGEMM('N','N',IRANK,KF,IN,1.0_JPRD,& - & YNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1,1),YD_STRUCT%N_ORDER,1.0_JPRD,& + CALL GEMM('N','N',IRANK,KF,IN,1.0_JPRB,& + & YNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1,1),YD_STRUCT%N_ORDER,1.0_JPRB,& & ZBETA(IBTST,1,IBETALV),ILBETA) - ELSE - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('N','N',IRANK,KF,IN,1.0_JPRM,& - & YNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1,1),YD_STRUCT%N_ORDER,1.0_JPRM,& - & ZBETA(IBTST,1,IBETALV),ILBETA) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) - ENDIF ENDIF ELSE ILM1 = JL-1 @@ -1011,20 +965,9 @@ SUBROUTINE MULT_BUTM(CDTRANS,YD_STRUCT,KF,PVECIN,PVECOUT,KWV) ENDDO IF( IRANK <=0 ) CALL ABOR1('mult_butm: IRANK<=0 not allowed') IF(YNODE%ICOLS > IRANK) THEN - IF (LLDOUBLE) THEN - CALL DGEMM('N','N',IRANK,KF,IN,1.0_JPRD,& - & YNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1,1),YD_STRUCT%N_ORDER,1.0_JPRD,& + CALL GEMM('N','N',IRANK,KF,IN,1.0_JPRB,& + & YNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1,1),YD_STRUCT%N_ORDER,1.0_JPRB,& & ZBETA(IBTST,1,IBETALV),ILBETA) - ELSE - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('N','N',IRANK,KF,IN,1.0_JPRM,& - & YNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1,1),YD_STRUCT%N_ORDER,1.0_JPRM,& - & ZBETA(IBTST,1,IBETALV),ILBETA) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) - ENDIF ENDIF ENDIF IF( IRANK <=0 ) CALL ABOR1('mult_butm: IRANK<=0 not allowed') @@ -1032,20 +975,9 @@ SUBROUTINE MULT_BUTM(CDTRANS,YD_STRUCT,KF,PVECIN,PVECOUT,KWV) IFR = YNODE%IFROW ILR = YNODE%ILROW IROWS = YNODE%IROWS - IF (LLDOUBLE) THEN - CALL DGEMM('N','N',IROWS,KF,YNODE%IRANK,1.0_JPRD,& - & YNODE%B,IROWS,ZBETA(IBTST,1,IBETALV),YD_STRUCT%IBETALEN_MAX,0.0_JPRD,& + CALL GEMM('N','N',IROWS,KF,YNODE%IRANK,1.0_JPRB,& + & YNODE%B(1,1),IROWS,ZBETA(IBTST,1,IBETALV),YD_STRUCT%IBETALEN_MAX,0.0_JPRB,& & PVECOUT(IFR,1),IROUT) - ELSE - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('N','N',IROWS,KF,YNODE%IRANK,1.0_JPRM,& - & YNODE%B,IROWS,ZBETA(IBTST,1,IBETALV),YD_STRUCT%IBETALEN_MAX,0.0_JPRM,& - & PVECOUT(IFR,1),IROUT) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) - ENDIF ENDIF ENDDO ENDDO @@ -1079,10 +1011,10 @@ SUBROUTINE MULT_P(YDNODE,PVECIN,PVECOUT) IM = IRANK IN = YDNODE%ICOLS-IRANK IF (JPRB == JPRD) THEN - CALL DGEMV('N',IM,IN,1.0_JPRD,YDNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1),1,1.0_JPRD,ZVECOUT,1) + CALL GEMV('N',IM,IN,1.0_JPRD,YDNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1),1,1.0_JPRD,ZVECOUT(1),1) PVECOUT(:)=ZVECOUT(:) ELSE - CALL SGEMV('N',IM,IN,1.0_JPRM,YDNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1),1,1.0_JPRM,PVECOUT,1) + CALL GEMV('N',IM,IN,1.0_JPRM,YDNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1),1,1.0_JPRM,PVECOUT(1),1) ENDIF ENDIF @@ -1100,7 +1032,6 @@ SUBROUTINE MULT_PM(YDNODE,KF,KLBETA,PVECIN,PVECOUT) REAL(KIND=JPRB) :: ZVECIN(YDNODE%ICOLS,KF), ZVECOUT(SIZE(PVECOUT(:,1)),KF) INTEGER(KIND=JPIM) :: JN,IDX,IRANK,IM,IN,JF -LOGICAL :: LL_HALT_INVALID !--------------------------------------------------------- IRANK = YDNODE%IRANK @@ -1117,20 +1048,9 @@ SUBROUTINE MULT_PM(YDNODE,KF,KLBETA,PVECIN,PVECOUT) ENDDO ENDDO IF(YDNODE%ICOLS > IRANK) THEN - IF (JPRB == JPRD) THEN - CALL DGEMM('N','N',IRANK,KF,IN,1.0_JPRD,& - & YDNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1,1),YDNODE%ICOLS,1.0_JPRD,& - & PVECOUT,IRANK) - ELSE - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('N','N',IRANK,KF,IN,1.0_JPRM,& - & YDNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1,1),YDNODE%ICOLS,1.0_JPRM,& - & PVECOUT,IRANK) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) - ENDIF + CALL GEMM('N','N',IRANK,KF,IN,1.0_JPRB,& + & YDNODE%PNONIM(1),IRANK,ZVECIN(IRANK+1,1),YDNODE%ICOLS,1.0_JPRB,& + & PVECOUT(1,1),IRANK) ENDIF END SUBROUTINE MULT_PM !================================================================== @@ -1151,9 +1071,9 @@ SUBROUTINE MULT_P_TR(YDNODE,PVECIN,PVECOUT) IM = IRANK IF (JPRB == JPRD) THEN ZVECIN(:) = PVECIN(:) - CALL DGEMV('T',IM,IN,1.0_JPRD,YDNODE%PNONIM,IRANK,ZVECIN,1,0.0_JPRD,ZVECOUT(IRANK+1),1) + CALL GEMV('T',IM,IN,1.0_JPRD,YDNODE%PNONIM(1),IRANK,ZVECIN(1),1,0.0_JPRD,ZVECOUT(IRANK+1),1) ELSE - CALL SGEMV('T',IM,IN,1.0_JPRM,YDNODE%PNONIM,IRANK,PVECIN,1,0.0_JPRM,ZVECOUT(IRANK+1),1) + CALL GEMV('T',IM,IN,1.0_JPRM,YDNODE%PNONIM(1),IRANK,PVECIN(1),1,0.0_JPRM,ZVECOUT(IRANK+1),1) ENDIF ENDIF DO JK=1,IRANK @@ -1178,8 +1098,6 @@ SUBROUTINE MULT_P_TRM(YDNODE,KF,PVECIN,PVECOUT) REAL(KIND=JPRB) :: ZVECOUT(YDNODE%ICOLS,KF), ZVECIN(SIZE(PVECIN(:,1)),KF) INTEGER(KIND=JPIM) :: JK,JN,IDX,IM,IN,JF -LOGICAL :: LL_HALT_INVALID - !------------------------------------------------------------------ IN = YDNODE%ICOLS-YDNODE%IRANK @@ -1187,18 +1105,9 @@ SUBROUTINE MULT_P_TRM(YDNODE,KF,PVECIN,PVECOUT) IF(IN>0) THEN IF (JPRB == JPRD) THEN ZVECIN(:,:) = PVECIN(:,:) - CALL DGEMM('T','N',IN,KF,IM,1.0_JPRD,& - & YDNODE%PNONIM(1),IM,ZVECIN,IM,0.0_JPRD,& - & ZVECOUT(YDNODE%IRANK+1,1),YDNODE%ICOLS) + CALL GEMM('T','N',IN,KF,IM,1.0_JPRD,YDNODE%PNONIM(1),IM,ZVECIN(1,1),IM,0.0_JPRD,ZVECOUT(YDNODE%IRANK+1,1),YDNODE%ICOLS) ELSE - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('T','N',IN,KF,IM,1.0_JPRM,& - & YDNODE%PNONIM(1),IM,PVECIN,IM,0.0_JPRM,& - & ZVECOUT(YDNODE%IRANK+1,1),YDNODE%ICOLS) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) + CALL GEMM('T','N',IN,KF,IM,1.0_JPRM,YDNODE%PNONIM(1),IM,PVECIN(1,1),IM,0.0_JPRM,ZVECOUT(YDNODE%IRANK+1,1),YDNODE%ICOLS) ENDIF ENDIF DO JF=1,KF diff --git a/src/trans/algor/ectrans_blas_mod.F90 b/src/trans/algor/ectrans_blas_mod.F90 new file mode 100644 index 000000000..256836dff --- /dev/null +++ b/src/trans/algor/ectrans_blas_mod.F90 @@ -0,0 +1,310 @@ +! (C) Copyright 2024- ECMWF. +! +! This software is licensed under the terms of the Apache Licence Version 2.0 +! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +! In applying this licence, ECMWF does not waive the privileges and immunities +! granted to it by virtue of its status as an intergovernmental organisation +! nor does it submit to any jurisdiction. +! + +!==================================================================== +MODULE ECTRANS_BLAS_MOD +!==================================================================== +! Author: Willem Deconinck (ECMWF) +! +! This module provides interfaces for BLAS routines such as +! DGEMM/SGEMM and DGEMV/SGEMV +! The correct overload is used depending on the precision of the arguments +!==================================================================== + + +USE EC_PARKIND, ONLY : JPRD, JPRM, JPIM + +IMPLICIT NONE + +PRIVATE + +PUBLIC :: GEMM, GEMV + +!--------------------------------------------------------------------- + +INTERFACE GEMM +! GEMM performs one of the matrix-matrix operations +! +! C := alpha*op( A )*op( B ) + beta*C, +! +! where op( X ) is one of +! +! op( X ) = X or op( X ) = X**T, +! +! alpha and beta are scalars, and A, B and C are matrices, with op( A ) +! an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. + + ! SGEMM + MODULE PROCEDURE GEMM_SP ! Matrix arguments as array, (alpha,beta) in SP + MODULE PROCEDURE GEMM_SP_DP ! Matrix arguments as array, (alpha,beta) in DP + MODULE PROCEDURE GEMM_SCAL_SP ! Matrix arguments as scalar (address), (alpha,beta) in SP + MODULE PROCEDURE GEMM_SCAL_SP_DP ! Matrix arguments as scalar (address), (alpha,beta) in DP + + ! DGEMM + MODULE PROCEDURE GEMM_DP ! Matrix arguments as array, (alpha,beta) in DP + MODULE PROCEDURE GEMM_DP_SP ! Matrix arguments as array, (alpha,beta) in SP + MODULE PROCEDURE GEMM_SCAL_DP ! Matrix arguments as scalar (address), (alpha,beta) in DP + MODULE PROCEDURE GEMM_SCAL_DP_SP ! Matrix arguments as scalar (address), (alpha,beta) in SP +END INTERFACE + +!--------------------------------------------------------------------- + +INTERFACE GEMV +! GEMV performs one of the matrix-vector operations +! +! y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y, +! +! where alpha and beta are scalars, x and y are vectors and A is an +! m by n matrix. + + ! SGEMV + MODULE PROCEDURE GEMV_SP ! Matrix/Vector arguments as array, (alpha,beta) in SP + MODULE PROCEDURE GEMV_SP_DP ! Matrix/Vector arguments as array, (alpha,beta) in DP + MODULE PROCEDURE GEMV_SCAL_SP ! Matrix/Vector arguments as scalar (address), (alpha,beta) in SP + MODULE PROCEDURE GEMV_SCAL_SP_DP ! Matrix/Vector arguments as scalar (address), (alpha,beta) in DP + + ! DGEMV + MODULE PROCEDURE GEMV_DP ! Matrix/Vector arguments as array, (alpha,beta) in DP + MODULE PROCEDURE GEMV_DP_SP ! Matrix/Vector arguments as array, (alpha,beta) in SP + MODULE PROCEDURE GEMV_SCAL_DP ! Matrix/Vector arguments as scalar (address), (alpha,beta) in DP + MODULE PROCEDURE GEMV_SCAL_DP_SP ! Matrix/Vector arguments as scalar (address), (alpha,beta) in SP +END INTERFACE + +!--------------------------------------------------------------------- + +!==================================================================== +CONTAINS +!==================================================================== + +SUBROUTINE GEMM_SCAL_DP(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) + REAL(KIND=JPRD) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: K, LDA, LDB, LDC, M, N + CHARACTER ,INTENT(IN) :: TRANSA, TRANSB + REAL(KIND=JPRD) ,INTENT(IN) :: A, B + REAL(KIND=JPRD) ,INTENT(INOUT) :: C + + CALL DGEMM(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) + +END SUBROUTINE GEMM_SCAL_DP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMM_SCAL_DP_SP(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) + REAL(KIND=JPRM) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: K, LDA, LDB, LDC, M, N + CHARACTER ,INTENT(IN) :: TRANSA, TRANSB + REAL(KIND=JPRD) ,INTENT(IN) :: A, B + REAL(KIND=JPRD) ,INTENT(INOUT) :: C + + CALL GEMM_SCAL_DP(TRANSA, TRANSB, M, N, K, REAL(ALPHA,JPRD), A, LDA, B, LDB, REAL(BETA,JPRD), C, LDC) + +END SUBROUTINE GEMM_SCAL_DP_SP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMM_DP(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) + REAL(KIND=JPRD) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: K, LDA, LDB, LDC, M, N + CHARACTER ,INTENT(IN) :: TRANSA, TRANSB + REAL(KIND=JPRD) ,INTENT(IN) :: A(LDA,*), B(LDB,*) + REAL(KIND=JPRD) ,INTENT(INOUT) :: C(LDC,*) + + CALL GEMM_SCAL_DP(TRANSA, TRANSB, M, N, K, ALPHA, A(1,1), LDA, B(1,1), LDB, BETA, C(1,1), LDC) + +END SUBROUTINE GEMM_DP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMM_DP_SP(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) + REAL(KIND=JPRM) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: K, LDA, LDB, LDC, M, N + CHARACTER ,INTENT(IN) :: TRANSA, TRANSB + REAL(KIND=JPRD) ,INTENT(IN) :: A(LDA,*), B(LDB,*) + REAL(KIND=JPRD) ,INTENT(INOUT) :: C(LDC,*) + + CALL GEMM_SCAL_DP(TRANSA, TRANSB, M, N, K, REAL(ALPHA,JPRD), A(1,1), LDA, B(1,1), LDB, REAL(BETA,JPRD), C(1,1), LDC) + +END SUBROUTINE GEMM_DP_SP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMM_SCAL_SP(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) + USE, INTRINSIC :: IEEE_EXCEPTIONS, ONLY : IEEE_GET_HALTING_MODE, IEEE_SET_HALTING_MODE, IEEE_INVALID + + REAL(KIND=JPRM) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: K, LDA, LDB, LDC, M, N + CHARACTER ,INTENT(IN) :: TRANSA, TRANSB + REAL(KIND=JPRM) ,INTENT(IN) :: A, B + REAL(KIND=JPRM) ,INTENT(INOUT) :: C + +#ifdef WITH_IEEE_HALT + LOGICAL, PARAMETER :: LL_IEEE_HALT = .TRUE. +#else + LOGICAL, PARAMETER :: LL_IEEE_HALT = .FALSE. +#endif + LOGICAL :: LL_HALT_INVALID = .FALSE. + + IF (LL_IEEE_HALT) THEN + CALL IEEE_GET_HALTING_MODE(IEEE_INVALID,LL_HALT_INVALID) + IF (LL_HALT_INVALID) CALL IEEE_SET_HALTING_MODE(IEEE_INVALID, .FALSE.) + ENDIF + + CALL SGEMM(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) + + IF (LL_IEEE_HALT .AND. LL_HALT_INVALID) CALL IEEE_SET_HALTING_MODE(IEEE_INVALID, .TRUE.) + +END SUBROUTINE GEMM_SCAL_SP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMM_SCAL_SP_DP(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) + REAL(KIND=JPRD) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: K, LDA, LDB, LDC, M, N + CHARACTER ,INTENT(IN) :: TRANSA, TRANSB + REAL(KIND=JPRM) ,INTENT(IN) :: A, B + REAL(KIND=JPRM) ,INTENT(INOUT) :: C + + CALL GEMM_SCAL_SP(TRANSA, TRANSB, M, N, K, REAL(ALPHA,JPRM), A, LDA, B, LDB, REAL(BETA,JPRM), C, LDC) + +END SUBROUTINE GEMM_SCAL_SP_DP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMM_SP(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) + REAL(KIND=JPRM) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: K, LDA, LDB, LDC, M, N + CHARACTER ,INTENT(IN) :: TRANSA, TRANSB + REAL(KIND=JPRM) ,INTENT(IN) :: A(LDA,*), B(LDB,*) + REAL(KIND=JPRM) ,INTENT(INOUT) :: C(LDC,*) + + CALL GEMM_SCAL_SP(TRANSA, TRANSB, M, N, K, ALPHA, A(1,1), LDA, B(1,1), LDB, BETA, C(1,1), LDC) + +END SUBROUTINE GEMM_SP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMM_SP_DP(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) + REAL(KIND=JPRD) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: K, LDA, LDB, LDC, M, N + CHARACTER ,INTENT(IN) :: TRANSA, TRANSB + REAL(KIND=JPRM) ,INTENT(IN) :: A(LDA,*), B(LDB,*) + REAL(KIND=JPRM) ,INTENT(INOUT) :: C(LDC,*) + + CALL GEMM_SCAL_SP(TRANSA, TRANSB, M, N, K, REAL(ALPHA,JPRM), A(1,1), LDA, B(1,1), LDB, REAL(BETA,JPRM), C(1,1), LDC) + +END SUBROUTINE GEMM_SP_DP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMV_SCAL_SP(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY) + REAL(KIND=JPRM) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: LDA, INCX, INCY, M, N + CHARACTER ,INTENT(IN) :: TRANS + REAL(KIND=JPRM) ,INTENT(IN) :: A, X + REAL(KIND=JPRM) ,INTENT(INOUT) :: Y + + CALL SGEMV(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY) + +END SUBROUTINE GEMV_SCAL_SP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMV_SCAL_SP_DP(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY) + REAL(KIND=JPRD) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: LDA, INCX, INCY, M, N + CHARACTER ,INTENT(IN) :: TRANS + REAL(KIND=JPRM) ,INTENT(IN) :: A, X + REAL(KIND=JPRM) ,INTENT(INOUT) :: Y + + CALL GEMV_SCAL_SP(TRANS, M, N, REAL(ALPHA,JPRM), A, LDA, X, INCX, REAL(BETA,JPRM), Y, INCY) + +END SUBROUTINE GEMV_SCAL_SP_DP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMV_SP(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY) + REAL(KIND=JPRM) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: LDA, INCX, INCY, M, N + CHARACTER ,INTENT(IN) :: TRANS + REAL(KIND=JPRM) ,INTENT(IN) :: A(:,:), X(:) + REAL(KIND=JPRM) ,INTENT(INOUT) :: Y(:) + + CALL GEMV_SCAL_SP(TRANS, M, N, ALPHA, A(1,1), LDA, X(1), INCX, BETA, Y(1), INCY) + +END SUBROUTINE GEMV_SP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMV_SP_DP(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY) + REAL(KIND=JPRD) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: LDA, INCX, INCY, M, N + CHARACTER ,INTENT(IN) :: TRANS + REAL(KIND=JPRM) ,INTENT(IN) :: A(:,:), X(:) + REAL(KIND=JPRM) ,INTENT(INOUT) :: Y(:) + + CALL GEMV_SCAL_SP(TRANS, M, N, REAL(ALPHA,JPRM), A(1,1), LDA, X(1), INCX, REAL(BETA,JPRM), Y(1), INCY) + +END SUBROUTINE GEMV_SP_DP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMV_SCAL_DP(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY) + REAL(KIND=JPRD) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: LDA, INCX, INCY, M, N + CHARACTER ,INTENT(IN) :: TRANS + REAL(KIND=JPRD) ,INTENT(IN) :: A, X + REAL(KIND=JPRD) ,INTENT(INOUT) :: Y + + CALL DGEMV(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY) + +END SUBROUTINE GEMV_SCAL_DP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMV_SCAL_DP_SP(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY) + REAL(KIND=JPRM) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: LDA, INCX, INCY, M, N + CHARACTER ,INTENT(IN) :: TRANS + REAL(KIND=JPRD) ,INTENT(IN) :: A, X + REAL(KIND=JPRD) ,INTENT(INOUT) :: Y + + CALL GEMV_SCAL_DP(TRANS, M, N, REAL(ALPHA,JPRD), A, LDA, X, INCX, REAL(BETA,JPRD), Y, INCY) + +END SUBROUTINE GEMV_SCAL_DP_SP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMV_DP(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY) + REAL(KIND=JPRD) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: LDA, INCX, INCY, M, N + CHARACTER ,INTENT(IN) :: TRANS + REAL(KIND=JPRD) ,INTENT(IN) :: A(:,:), X(:) + REAL(KIND=JPRD) ,INTENT(INOUT) :: Y(:) + + CALL GEMV_SCAL_DP(TRANS, M, N, ALPHA, A(1,1), LDA, X(1), INCX, BETA, Y(1), INCY) + +END SUBROUTINE GEMV_DP + +!--------------------------------------------------------------------- + +SUBROUTINE GEMV_DP_SP(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY) + REAL(KIND=JPRM) ,INTENT(IN) :: ALPHA, BETA + INTEGER(KIND=JPIM) ,INTENT(IN) :: LDA, INCX, INCY, M, N + CHARACTER ,INTENT(IN) :: TRANS + REAL(KIND=JPRD) ,INTENT(IN) :: A(:,:), X(:) + REAL(KIND=JPRD) ,INTENT(INOUT) :: Y(:) + + CALL GEMV_SCAL_DP(TRANS, M, N, REAL(ALPHA,JPRD), A(1,1), LDA, X(1), INCX, REAL(BETA,JPRD), Y(1), INCY) + +END SUBROUTINE GEMV_DP_SP + +!==================================================================== + +END MODULE ECTRANS_BLAS_MOD + diff --git a/src/trans/internal/ledir_mod.F90 b/src/trans/internal/ledir_mod.F90 index 2fa0ee5ea..9c0e4e009 100644 --- a/src/trans/internal/ledir_mod.F90 +++ b/src/trans/internal/ledir_mod.F90 @@ -61,9 +61,7 @@ SUBROUTINE LEDIR(KM,KMLOC,KFC,KIFC,KSL,KDGLU,KLED2,PAIA,PSIA,POA1,PW) USE TPM_DIM ,ONLY : R USE TPM_FLT ,ONLY : S USE BUTTERFLY_ALG_MOD, ONLY : MULT_BUTM - -use, intrinsic :: ieee_exceptions - +USE ECTRANS_BLAS_MOD, ONLY : GEMM IMPLICIT NONE @@ -86,12 +84,6 @@ SUBROUTINE LEDIR(KM,KMLOC,KFC,KIFC,KSL,KDGLU,KLED2,PAIA,PSIA,POA1,PW) INTEGER(KIND=JPIM) :: ITHRESHOLD REAL(KIND=JPRB) :: ZB(KDGLU,KIFC), ZCA((R%NTMAX-KM+2)/2,KIFC), ZCS((R%NTMAX-KM+3)/2,KIFC) REAL(KIND=JPRD), allocatable :: ZB_D(:,:), ZCA_D(:,:), ZCS_D(:,:),ZRPNMA(:,:), ZRPNMS(:,:) -LOGICAL :: LL_HALT_INVALID -#ifdef WITH_IEEE_HALT -LOGICAL, PARAMETER :: LL_IEEE_HALT = .TRUE. -#else -LOGICAL, PARAMETER :: LL_IEEE_HALT = .FALSE. -#endif LOGICAL, PARAMETER :: LLDOUBLE = (JPRB == JPRD) CHARACTER(LEN=1) :: CLX REAL(KIND=JPHOOK) :: ZHOOK_HANDLE @@ -136,17 +128,12 @@ SUBROUTINE LEDIR(KM,KMLOC,KFC,KIFC,KSL,KDGLU,KLED2,PAIA,PSIA,POA1,PW) IF (LHOOK) CALL DR_HOOK('LEDIR_'//CLX//'GEMM_1',0,ZHOOK_HANDLE) IF (LLDOUBLE) THEN - CALL DGEMM('T','N',ILA,KIFC,KDGLU,1.0_JPRD,S%FA(KMLOC)%RPNMA,KDGLU,& + CALL GEMM('T','N',ILA,KIFC,KDGLU,1.0_JPRD,S%FA(KMLOC)%RPNMA,KDGLU,& &ZB,KDGLU,0._JPRD,ZCA,ILA) ELSE IF(KM>=1)THEN ! DGEM for the mean to improve mass conservation - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('T','N',ILA,KIFC,KDGLU,1.0_JPRM,S%FA(KMLOC)%RPNMA,KDGLU,& + CALL GEMM('T','N',ILA,KIFC,KDGLU,1.0_JPRM,S%FA(KMLOC)%RPNMA,KDGLU,& &ZB,KDGLU,0._JPRM,ZCA,ILA) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) ELSE I1 = size(S%FA(KMLOC)%RPNMA(:,1)) I2 = size(S%FA(KMLOC)%RPNMA(1,:)) @@ -165,7 +152,7 @@ SUBROUTINE LEDIR(KM,KMLOC,KFC,KIFC,KSL,KDGLU,KLED2,PAIA,PSIA,POA1,PW) ZRPNMA(I3,I4) = S%FA(KMLOC)%RPNMA(I3,I4) END DO END DO - CALL DGEMM('T','N',ILA,KIFC,KDGLU,1.0_JPRD,ZRPNMA,KDGLU,& + CALL GEMM('T','N',ILA,KIFC,KDGLU,1.0_JPRD,ZRPNMA,KDGLU,& &ZB_D,KDGLU,0._JPRD,ZCA_D,ILA) IFLD=0 DO JK=1,KFC,ISKIP @@ -210,17 +197,12 @@ SUBROUTINE LEDIR(KM,KMLOC,KFC,KIFC,KSL,KDGLU,KLED2,PAIA,PSIA,POA1,PW) IF (LHOOK) CALL DR_HOOK('LEDIR_'//CLX//'GEMM_2',0,ZHOOK_HANDLE) IF (LLDOUBLE) THEN - CALL DGEMM('T','N',ILS,KIFC,KDGLU,1.0_JPRD,S%FA(KMLOC)%RPNMS,KDGLU,& + CALL GEMM('T','N',ILS,KIFC,KDGLU,1.0_JPRD,S%FA(KMLOC)%RPNMS,KDGLU,& &ZB,KDGLU,0._JPRD,ZCS,ILS) ELSE IF(KM>=1)THEN ! DGEM for the mean to improve mass conservation - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('T','N',ILS,KIFC,KDGLU,1.0_JPRM,S%FA(KMLOC)%RPNMS,KDGLU,& + CALL GEMM('T','N',ILS,KIFC,KDGLU,1.0_JPRM,S%FA(KMLOC)%RPNMS,KDGLU,& &ZB,KDGLU,0._JPRM,ZCS,ILS) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) ELSE I1 = size(S%FA(KMLOC)%RPNMS(:,1)) I2 = size(S%FA(KMLOC)%RPNMS(1,:)) @@ -239,7 +221,7 @@ SUBROUTINE LEDIR(KM,KMLOC,KFC,KIFC,KSL,KDGLU,KLED2,PAIA,PSIA,POA1,PW) ZRPNMS(I3,I4) = S%FA(KMLOC)%RPNMS(I3,I4) END DO END DO - CALL DGEMM('T','N',ILS,KIFC,KDGLU,1.0_JPRD,ZRPNMS,KDGLU,& + CALL GEMM('T','N',ILS,KIFC,KDGLU,1.0_JPRD,ZRPNMS,KDGLU,& &ZB_D,KDGLU,0._JPRD,ZCS_D,ILS) IFLD=0 DO JK=1,KFC,ISKIP diff --git a/src/trans/internal/ledirad_mod.F90 b/src/trans/internal/ledirad_mod.F90 index fa8a3b955..96decfe87 100644 --- a/src/trans/internal/ledirad_mod.F90 +++ b/src/trans/internal/ledirad_mod.F90 @@ -68,6 +68,7 @@ SUBROUTINE LEDIRAD(KM,KMLOC,KFC,KIFC,KDGLU,KLED2,PAIA,PSIA,POA1) USE TPM_FLT ,ONLY : S USE TPM_FIELDS ,ONLY : F USE BUTTERFLY_ALG_MOD ,ONLY: MULT_BUTM +USE ECTRANS_BLAS_MOD, ONLY : GEMM IMPLICIT NONE @@ -86,7 +87,6 @@ SUBROUTINE LEDIRAD(KM,KMLOC,KFC,KIFC,KDGLU,KLED2,PAIA,PSIA,POA1) INTEGER(KIND=JPIM) :: IA, ILA, ILS, IS, ISKIP, ISL, J, JK,JGL,J1 INTEGER(KIND=JPIM) :: IFLD,ITHRESHOLD REAL(KIND=JPRB) :: ZB(KDGLU,KIFC), ZCA((R%NTMAX-KM+2)/2,KIFC), ZCS((R%NTMAX-KM+3)/2,KIFC) -LOGICAL, PARAMETER :: LLDOUBLE = (JPRD == JPRB) CHARACTER(LEN=1) :: CLX REAL(KIND=JPHOOK) :: ZHOOK_HANDLE @@ -98,7 +98,7 @@ SUBROUTINE LEDIRAD(KM,KMLOC,KFC,KIFC,KDGLU,KLED2,PAIA,PSIA,POA1) !* 1.1 PREPARATIONS. CLX = 'S' -IF (LLDOUBLE) CLX = 'D' +IF (JPRD == JPRB) CLX = 'D' IA = 1+MOD(R%NTMAX-KM+2,2) IS = 1+MOD(R%NTMAX-KM+1,2) @@ -134,15 +134,9 @@ SUBROUTINE LEDIRAD(KM,KMLOC,KFC,KIFC,KDGLU,KLED2,PAIA,PSIA,POA1) ENDDO IF(ILA <= ITHRESHOLD .OR. .NOT.S%LUSEFLT) THEN - IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_1',0,ZHOOK_HANDLE) - IF(LLDOUBLE)THEN - CALL DGEMM('N','N',KDGLU,KIFC,ILA,1.0_JPRD,S%FA(KMLOC)%RPNMA,KDGLU,& - &ZCA,ILA,0._JPRD,ZB,KDGLU) - ELSE - CALL SGEMM('N','N',KDGLU,KIFC,ILA,1.0_JPRM,S%FA(KMLOC)%RPNMA,KDGLU,& - &ZCA,ILA,0._JPRM,ZB,KDGLU) - END IF - IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_1',1,ZHOOK_HANDLE) + IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_1',0,ZHOOK_HANDLE) + CALL GEMM('N','N',KDGLU,KIFC,ILA,1.0_JPRB,S%FA(KMLOC)%RPNMA,KDGLU,ZCA,ILA,0._JPRB,ZB,KDGLU) + IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_1',1,ZHOOK_HANDLE) ELSE @@ -173,14 +167,7 @@ SUBROUTINE LEDIRAD(KM,KMLOC,KFC,KIFC,KDGLU,KLED2,PAIA,PSIA,POA1) IF(ILS <= ITHRESHOLD .OR. .NOT.S%LUSEFLT) THEN IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_2',0,ZHOOK_HANDLE) - IF(LLDOUBLE)THEN - CALL DGEMM('N','N',KDGLU,KIFC,ILS,1.0_JPRD,S%FA(KMLOC)%RPNMS,KDGLU,& - &ZCS,ILS,0._JPRD,ZB,KDGLU) - ELSE - CALL SGEMM('N','N',KDGLU,KIFC,ILS,1.0_JPRM,S%FA(KMLOC)%RPNMS,KDGLU,& - &ZCS,ILS,0._JPRM,ZB,KDGLU) - - END IF + CALL GEMM('N','N',KDGLU,KIFC,ILS,1.0_JPRD,S%FA(KMLOC)%RPNMS,KDGLU,ZCS,ILS,0._JPRD,ZB,KDGLU) IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_2',1,ZHOOK_HANDLE) ELSE diff --git a/src/trans/internal/leinv_mod.F90 b/src/trans/internal/leinv_mod.F90 index a753c96c8..5d6a33080 100644 --- a/src/trans/internal/leinv_mod.F90 +++ b/src/trans/internal/leinv_mod.F90 @@ -50,7 +50,7 @@ SUBROUTINE LEINV(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KSL,KDGLU,PIA,PAOA1,PSOA1) ! ! Modifications. ! -------------- -! J.Hague : Oct 2012 DR_HOOK round calls to DGEMM: +! J.Hague : Oct 2012 DR_HOOK round calls to GEMM: ! F. Vana 05-Mar-2015 Support for single precision ! ------------------------------------------------------------------ @@ -60,8 +60,7 @@ SUBROUTINE LEINV(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KSL,KDGLU,PIA,PAOA1,PSOA1) USE TPM_DIM ,ONLY : R USE TPM_FLT ,ONLY : S USE BUTTERFLY_ALG_MOD, ONLY : MULT_BUTM - -use, intrinsic :: ieee_exceptions +USE ECTRANS_BLAS_MOD, ONLY : GEMM IMPLICIT NONE @@ -80,13 +79,6 @@ SUBROUTINE LEINV(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KSL,KDGLU,PIA,PAOA1,PSOA1) INTEGER(KIND=JPIM) :: IA, ILA, ILS, IS, ISKIP, ISL, J1, IFLD, JGL,JK, J,JI, IEND INTEGER(KIND=JPIM) :: ITHRESHOLD REAL(KIND=JPRB) :: ZBA((R%NSMAX-KM+2)/2,KIFC), ZBS((R%NSMAX-KM+3)/2,KIFC), ZC(KDGLU,KIFC) -LOGICAL :: LL_HALT_INVALID -#ifdef WITH_IEEE_HALT -LOGICAL, PARAMETER :: LL_IEEE_HALT = .TRUE. -#else -LOGICAL, PARAMETER :: LL_IEEE_HALT = .FALSE. -#endif -LOGICAL, PARAMETER :: LLDOUBLE = (JPRB == JPRD) CHARACTER(LEN=1) :: CLX REAL(KIND=JPHOOK) :: ZHOOK_HANDLE @@ -98,7 +90,7 @@ SUBROUTINE LEINV(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KSL,KDGLU,PIA,PAOA1,PSOA1) !* 1.1 PREPARATIONS. CLX = 'S' -IF (LLDOUBLE) CLX = 'D' +IF (JPRB == JPRD) CLX = 'D' !ISL = MAX(R%NDGNH-G%NDGLU(KM)+1,1) ISL = KSL @@ -138,20 +130,9 @@ SUBROUTINE LEINV(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KSL,KDGLU,PIA,PAOA1,PSOA1) IF(ILA <= ITHRESHOLD .OR. .NOT.S%LUSEFLT) THEN IF (LHOOK) CALL DR_HOOK('LEINV_'//CLX//'GEMM_1',0,ZHOOK_HANDLE) - IF (LLDOUBLE) THEN - CALL DGEMM('N','N',KDGLU,KIFC,ILA,1.0_JPRD,S%FA(KMLOC)%RPNMA,KDGLU,& - &ZBA,ILA,0._JPRD,ZC,KDGLU) - ELSE - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('N','N',KDGLU,KIFC,ILA,1.0_JPRM,S%FA(KMLOC)%RPNMA,KDGLU,& - &ZBA,ILA,0._JPRM,ZC,KDGLU) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) - ENDIF + CALL GEMM('N','N',KDGLU,KIFC,ILA,1.0_JPRB,S%FA(KMLOC)%RPNMA,KDGLU,ZBA,ILA,0._JPRB,ZC,KDGLU) IF (LHOOK) CALL DR_HOOK('LEINV_'//CLX//'GEMM_1',1,ZHOOK_HANDLE) - + ELSE IF (LHOOK) CALL DR_HOOK('LEINV_'//CLX//'BUTM_1',0,ZHOOK_HANDLE) @@ -182,18 +163,7 @@ SUBROUTINE LEINV(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KSL,KDGLU,PIA,PAOA1,PSOA1) IF(ILS <= ITHRESHOLD .OR. .NOT.S%LUSEFLT ) THEN IF (LHOOK) CALL DR_HOOK('LEINV_'//CLX//'GEMM_2',0,ZHOOK_HANDLE) - IF (LLDOUBLE) THEN - CALL DGEMM('N','N',KDGLU,KIFC,ILS,1.0_JPRD,S%FA(KMLOC)%RPNMS,KDGLU,& - &ZBS,ILS,0._JPRD,ZC,KDGLU) - ELSE - IF (LL_IEEE_HALT) THEN - call ieee_get_halting_mode(ieee_invalid,LL_HALT_INVALID) - if (LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.false.) - ENDIF - CALL SGEMM('N','N',KDGLU,KIFC,ILS,1.0_JPRM,S%FA(KMLOC)%RPNMS,KDGLU,& - &ZBS,ILS,0._JPRM,ZC,KDGLU) - if (LL_IEEE_HALT .and. LL_HALT_INVALID) call ieee_set_halting_mode(ieee_invalid,.true.) - ENDIF + CALL GEMM('N','N',KDGLU,KIFC,ILS,1.0_JPRB,S%FA(KMLOC)%RPNMS,KDGLU,ZBS,ILS,0._JPRB,ZC,KDGLU) IF (LHOOK) CALL DR_HOOK('LEINV_'//CLX//'GEMM_2',1,ZHOOK_HANDLE) ELSE diff --git a/src/trans/internal/leinvad_mod.F90 b/src/trans/internal/leinvad_mod.F90 index b6f37c01a..9eb8f95a9 100644 --- a/src/trans/internal/leinvad_mod.F90 +++ b/src/trans/internal/leinvad_mod.F90 @@ -62,6 +62,7 @@ SUBROUTINE LEINVAD(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KDGLU,PIA,PAOA1,PSOA1) ! USE TPM_FLT ,ONLY : S USE BUTTERFLY_ALG_MOD,ONLY : MULT_BUTM +USE ECTRANS_BLAS_MOD, ONLY : GEMM IMPLICIT NONE @@ -79,7 +80,6 @@ SUBROUTINE LEINVAD(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KDGLU,PIA,PAOA1,PSOA1) INTEGER(KIND=JPIM) :: IA, ILA, ILS, IS, ISKIP, ISL, IOAD1, JK,JI INTEGER(KIND=JPIM) :: IFLD,ITHRESHOLD REAL(KIND=JPRB) :: ZBA((R%NSMAX-KM+2)/2,KIFC), ZBS((R%NSMAX-KM+3)/2,KIFC), ZC(KDGLU,KIFC) -LOGICAL, PARAMETER :: LLDOUBLE = (JPRD == JPRB) CHARACTER(LEN=1) :: CLX REAL(KIND=JPHOOK) :: ZHOOK_HANDLE @@ -91,7 +91,7 @@ SUBROUTINE LEINVAD(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KDGLU,PIA,PAOA1,PSOA1) !* 1.1 PREPARATIONS. CLX = 'S' -IF (LLDOUBLE) CLX = 'D' +IF (JPRD == JPRB) CLX = 'D' IA = 1+MOD(R%NSMAX-KM+2,2) IS = 1+MOD(R%NSMAX-KM+1,2) @@ -123,15 +123,9 @@ SUBROUTINE LEINVAD(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KDGLU,PIA,PAOA1,PSOA1) ENDDO IF(ILA <= ITHRESHOLD .OR. .NOT.S%LUSEFLT) THEN - IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_1',0,ZHOOK_HANDLE) - IF(LLDOUBLE)THEN - CALL DGEMM('T','N',ILA,KIFC,KDGLU,1.0_JPRD,S%FA(KMLOC)%RPNMA,KDGLU,& - &ZC,KDGLU,0._JPRD,ZBA,ILA) - ELSE - CALL SGEMM('T','N',ILA,KIFC,KDGLU,1.0_JPRM,S%FA(KMLOC)%RPNMA,KDGLU,& - &ZC,KDGLU,0._JPRM,ZBA,ILA) - END IF - IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_1',1,ZHOOK_HANDLE) + IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_1',0,ZHOOK_HANDLE) + CALL GEMM('T','N',ILA,KIFC,KDGLU,1.0_JPRB,S%FA(KMLOC)%RPNMA,KDGLU,ZC,KDGLU,0._JPRB,ZBA,ILA) + IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_1',1,ZHOOK_HANDLE) ELSE @@ -161,13 +155,7 @@ SUBROUTINE LEINVAD(KM,KMLOC,KFC,KIFC,KF_OUT_LT,KDGLU,PIA,PAOA1,PSOA1) IF(ILS <= ITHRESHOLD .OR. .NOT.S%LUSEFLT ) THEN IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_2',0,ZHOOK_HANDLE) - IF(LLDOUBLE)THEN - CALL DGEMM('T','N',ILS,KIFC,KDGLU,1.0_JPRD,S%FA(KMLOC)%RPNMS,KDGLU,& - &ZC,KDGLU,0._JPRD,ZBS,ILS) - ELSE - CALL SGEMM('T','N',ILS,KIFC,KDGLU,1.0_JPRM,S%FA(KMLOC)%RPNMS,KDGLU,& - &ZC,KDGLU,0._JPRM,ZBS,ILS) - END IF + CALL GEMM('T','N',ILS,KIFC,KDGLU,1.0_JPRB,S%FA(KMLOC)%RPNMS,KDGLU,ZC,KDGLU,0._JPRB,ZBS,ILS) IF (LHOOK) CALL DR_HOOK('LE_'//CLX//'GEMM_2',1,ZHOOK_HANDLE) ELSE