Ignore:
Timestamp:
Feb 3, 2009, 10:50:31 AM (15 years ago)
Author:
yann meurdesoif
Message:

Modifications Othman Bouzi : optimisation du filtre (remplacement opération matrice/vecteurs par matrice/matrice/matrice - BLAS), l'allocation des tableaux du filtre se fait maintenant dynamiquement (plus d'intervention manuelle dans parafiltre.h)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/filtreg_p.F

    r985 r1086  
    22
    33      SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
    4      .                       ifiltre, iaire, griscal ,iter)
    5       USE Parallel, only : OMP_CHUNK 
     4     &     ifiltre, iaire, griscal ,iter)
     5      USE Parallel, only : OMP_CHUNK
    66      USE mod_filtre_fft
    77      USE timer_filtre
     8     
     9      USE filtre
     10     
    811      IMPLICIT NONE
    9 
     12     
    1013c=======================================================================
    1114c
    … …  
    5053#include "dimensions.h"
    5154#include "paramet.h"
    52 #include "parafilt.h"
    5355#include "coefils.h"
    5456c
    … …  
    5759      INTEGER iim2,immjm
    5860      INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
    59 
     61     
    6062      REAL  champ( iip1,nlat,nbniv)
    61       REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs
    62       COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus)
    63      ,             , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs)
    64      ,             ,  matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus)
    65 cym      REAL  eignq(iim), sdd1(iim),sdd2(iim)
    66 
    67       REAL  eignq(iim)
    68       REAL :: sdd1(iim),sdd2(iim)
    6963     
    7064      LOGICAL    griscal
    … …  
    7468      REAL :: champ_in(iip1,nlat,nbniv)
    7569     
    76       REAL,SAVE,TARGET :: sddu_loc(iim)
    77       REAL,SAVE,TARGET :: sddv_loc(iim)
    78       REAL,SAVE,TARGET :: unsddu_loc(iim)
    79       REAL,SAVE,TARGET :: unsddv_loc(iim)
    80 c$OMP THREADPRIVATE(sddu_loc,sddv_loc,unsddu_loc,unsddv_loc)
    8170      LOGICAL,SAVE     :: first=.TRUE.
    8271c$OMP THREADPRIVATE(first)
    8372
     73      REAL, DIMENSION(iip1,nlat,nbniv) :: champ_loc
     74      INTEGER :: ll_nb, nbniv_loc
     75      REAL, SAVE :: sdd12(iim,4)
     76c$OMP THREADPRIVATE(sdd12)
     77
     78      INTEGER, PARAMETER :: type_sddu=1
     79      INTEGER, PARAMETER :: type_sddv=2
     80      INTEGER, PARAMETER :: type_unsddu=3
     81      INTEGER, PARAMETER :: type_unsddv=4
     82
     83      INTEGER :: sdd1_type, sdd2_type
     84
    8485      IF (first) THEN
    85         sddu_loc(1:iim)=sddu(1:iim)
    86         sddv_loc(1:iim)=sddv(1:iim)
    87         unsddu_loc(1:iim)=unsddu(1:iim)
    88         unsddv_loc(1:iim)=unsddv(1:iim)
    89         CALL Init_timer
    90         first=.FALSE.
    91 c       PRINT *,"----> sddu_loc=",sddu_loc
    92 c       PRINT *,"----> sddv_loc=",sddv_loc
    93 c       PRINT *,"----> unsddu_loc=",unsddu_loc
    94 c       PRINT *,"----> unsddv_loc=",unsddv_loc
     86         sdd12(1:iim,type_sddu) = sddu(1:iim)
     87         sdd12(1:iim,type_sddv) = sddv(1:iim)
     88         sdd12(1:iim,type_unsddu) = unsddu(1:iim)
     89         sdd12(1:iim,type_unsddv) = unsddv(1:iim)
     90
     91         CALL Init_timer
     92         first=.FALSE.
    9593      ENDIF
    9694
    … …  
    9997c$OMP END MASTER
    10098
     99c-------------------------------------------------------c
     100
    101101      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
    102      *    STOP'Pas de transformee simple dans cette version'
    103 
     102     &     STOP'Pas de transformee simple dans cette version'
     103     
    104104      IF( iter.EQ. 2 )  THEN
    105        PRINT *,' Pas d iteration du filtre dans cette version !'
    106      * , ' Utiliser old_filtreg et repasser !'
    107            STOP
     105         PRINT *,' Pas d iteration du filtre dans cette version !'
     106     &        , ' Utiliser old_filtreg et repasser !'
     107         STOP
    108108      ENDIF
    109109
    110110      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
    111        PRINT *,' Cette routine ne calcule le filtre inverse que ',
    112      * ' sur la grille des scalaires !'
    113            STOP
     111         PRINT *,' Cette routine ne calcule le filtre inverse que '
     112     &        , ' sur la grille des scalaires !'
     113         STOP
    114114      ENDIF
    115115
    116116      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
    117        PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
    118      *,' corriger et repasser !'
    119            STOP
     117         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
     118     &        , ' corriger et repasser !'
     119         STOP
    120120      ENDIF
    121121c
    … …  
    127127      IF( griscal )   THEN
    128128         IF( nlat. NE. jjp1 )  THEN
    129              PRINT  1111
    130              STOP
     129            PRINT  1111
     130            STOP
    131131         ELSE
    132 c
    133              IF( iaire.EQ.1 )  THEN
    134 cym                CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
    135 cym                CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
    136 cym               sdd1=>sddv_loc
    137 cym               sdd2=>unsddv_loc
    138                sdd1(1:iim)=sddv_loc(1:iim)
    139                sdd2(1:iim)=unsddv_loc(1:iim)
    140              ELSE
    141 cym                CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
    142 cym                CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
    143                sdd1(1:iim)=unsddv_loc(1:iim)
    144                sdd2(1:iim)=sddv_loc(1:iim)
    145              END IF
    146 c
    147              jdfil1 = 2
    148              jffil1 = jfiltnu
    149              jdfil2 = jfiltsu
    150              jffil2 = jjm
    151           END IF
     132c     
     133            IF( iaire.EQ.1 )  THEN
     134               sdd1_type = type_sddv
     135               sdd2_type = type_unsddv
     136            ELSE
     137               sdd1_type = type_unsddv
     138               sdd2_type = type_sddv
     139            ENDIF
     140c
     141            jdfil1 = 2
     142            jffil1 = jfiltnu
     143            jdfil2 = jfiltsu
     144            jffil2 = jjm
     145         ENDIF
    152146      ELSE
    153           IF( nlat.NE.jjm )  THEN
    154              PRINT  2222
    155              STOP
    156           ELSE
    157 c
    158              IF( iaire.EQ.1 )  THEN
    159 cym                CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
    160 cym                CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
    161 cym                sdd1=>sddu_loc
    162 cym                sdd2=>unsddu_loc
    163                 sdd1(1:iim)=sddu_loc(1:iim)
    164                 sdd2(1:iim)=unsddu_loc(1:iim)
    165 
    166              ELSE
    167 cym                CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
    168 cym                CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
    169 cym               sdd1=>unsddu_loc
    170 cym               sdd2=>sddu_loc
    171                sdd1(1:iim)=unsddu_loc(1:iim)
    172                sdd2(1:iim)=sddu_loc(1:iim)
    173              END IF
    174 c
    175              jdfil1 = 1
    176              jffil1 = jfiltnv
    177              jdfil2 = jfiltsv
    178              jffil2 = jjm
    179           END IF
    180       END IF
    181 
    182 c      PRINT *,"APPEL a filtreg --> sdd1=",sdd1
    183 c      PRINT *,"APPEL a filtreg --> sdd2=",sdd2
    184 c      PRINT *,"----> sddu_loc=",sddu_loc
    185 c       PRINT *,"----> sddv_loc=",sddv_loc
    186 c       PRINT *,"----> unsddu_loc=",unsddu_loc
    187 c       PRINT *,"----> unsddv_loc=",unsddv_loc
    188  
    189 c
    190 c
    191       DO 100  hemisph = 1, 2
    192 c
    193       IF ( hemisph.EQ.1 )  THEN
    194 c ym
    195           jdfil = max(jdfil1,ibeg)
    196           jffil = min(jffil1,iend)
    197       ELSE
    198 c ym
    199           jdfil = max(jdfil2,ibeg)
    200           jffil = min(jffil2,iend)
    201       END IF
     147         IF( nlat.NE.jjm )  THEN
     148            PRINT  2222
     149            STOP
     150         ELSE
     151c
     152            IF( iaire.EQ.1 )  THEN
     153               sdd1_type = type_sddu
     154               sdd2_type = type_unsddu
     155            ELSE
     156               sdd1_type = type_unsddu
     157               sdd2_type = type_sddu
     158            ENDIF
     159c     
     160            jdfil1 = 1
     161            jffil1 = jfiltnv
     162            jdfil2 = jfiltsv
     163            jffil2 = jjm
     164         ENDIF
     165      ENDIF
     166c     
     167      DO hemisph = 1, 2
     168c     
     169         IF ( hemisph.EQ.1 )  THEN
     170cym
     171            jdfil = max(jdfil1,ibeg)
     172            jffil = min(jffil1,iend)
     173         ELSE
     174cym
     175            jdfil = max(jdfil2,ibeg)
     176            jffil = min(jffil2,iend)
     177         ENDIF
    202178
    203179
    … …  
    206182cccccccccccccccccccccccccccccccccccccccccccc
    207183
    208       IF (.NOT. use_filtre_fft) THEN
    209      
    210 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
    211       DO 50  l = 1, nbniv
    212         DO 30  j = jdfil,jffil
    213  
    214  
    215           DO  5  i = 1, iim
    216             champ(i,j,l) = champ(i,j,l) * sdd1(i)
    217    5      CONTINUE
    218 c
    219 
    220           IF( hemisph. EQ. 1 )      THEN
    221 
    222             IF( ifiltre. EQ. -2 )   THEN
    223 
    224 
    225               CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
    226      .                     champ(1,j,l), 1, 0.0, eignq, 1)
    227 
    228 
    229             ELSE IF ( griscal )     THEN
    230 
    231               CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
    232      .                    champ(1,j,l), 1, 0.0, eignq, 1)
    233 
    234             ELSE
    235 
    236               CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
    237      .                   champ(1,j,l), 1, 0.0, eignq, 1)
    238             ENDIF
    239 
    240           ELSE
    241 
    242             IF( ifiltre. EQ. -2 )   THEN
    243      
    244               CALL SGEMV("N",iim,iim,1.0, matrinvs(1,1,j-jfiltsu+1),iim,
    245      .                   champ(1,j,l), 1, 0.0, eignq, 1)
    246      
    247             ELSE IF ( griscal )     THEN
    248      
    249               CALL SGEMV("N",iim,iim,1.0,matriceus(1,1,j-jfiltsu+1),iim,
    250      .                   champ(1,j,l), 1, 0.0, eignq, 1)
    251             ELSE
    252          
    253               CALL SGEMV("N",iim,iim,1.0,matricevs(1,1,j-jfiltsv+1),iim,
    254      .                    champ(1,j,l), 1, 0.0, eignq, 1)
    255             ENDIF
    256 
    257           ENDIF
    258 
    259 
    260 c
    261           IF( ifiltre.EQ. 2 )  THEN
    262          
    263             DO 15 i = 1, iim
    264               champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
    265   15        CONTINUE
    266          
    267           ELSE
    268        
    269             DO 16 i=1,iim
    270                champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
    271 16          CONTINUE
    272          
    273           ENDIF
    274 c
    275           champ( iip1,j,l ) = champ( 1,j,l )
    276 c
    277   30    CONTINUE
    278 c
    279   50  CONTINUE
     184         IF (.NOT. use_filtre_fft) THEN
     185     
     186c     !---------------------------------!
     187c     ! Agregation des niveau verticaux !
     188c     ! uniquement necessaire pour une  !
     189c     ! execution OpenMP                !
     190c     !---------------------------------!
     191            ll_nb = 0
     192c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     193            DO l = 1, nbniv
     194               ll_nb = ll_nb+1
     195               DO j = jdfil,jffil
     196                  DO i = 1, iim
     197                     champ_loc(i,j,ll_nb) =
     198     &                    champ(i,j,l) * sdd12(i,sdd1_type)
     199                  ENDDO
     200               ENDDO
     201            ENDDO
    280202c$OMP END DO NOWAIT
    281203
     204            nbniv_loc = ll_nb
     205
     206            IF( hemisph.EQ.1 )      THEN
     207               
     208               IF( ifiltre.EQ.-2 )   THEN
     209                  DO j = jdfil,jffil
     210                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
     211     &                    matrinvn(1,1,j), iim,
     212     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
     213     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     214                  ENDDO
     215                 
     216               ELSE IF ( griscal )     THEN
     217                  DO j = jdfil,jffil
     218                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
     219     &                    matriceun(1,1,j), iim,
     220     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
     221     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     222                  ENDDO
     223                 
     224               ELSE
     225                  DO j = jdfil,jffil
     226                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
     227     &                    matricevn(1,1,j), iim,
     228     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
     229     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     230                  ENDDO
     231                 
     232               ENDIF
     233               
     234            ELSE
     235               
     236               IF( ifiltre.EQ.-2 )   THEN
     237                  DO j = jdfil,jffil
     238                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
     239     &                    matrinvs(1,1,j-jfiltsu+1), iim,
     240     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
     241     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     242                  ENDDO
     243                 
     244               ELSE IF ( griscal )     THEN
     245                 
     246                  DO j = jdfil,jffil
     247                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
     248     &                    matriceus(1,1,j-jfiltsu+1), iim,
     249     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
     250     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     251                  ENDDO
     252                 
     253               ELSE
     254                 
     255                  DO j = jdfil,jffil
     256                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
     257     &                    matricevs(1,1,j-jfiltsv+1), iim,
     258     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
     259     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     260                  ENDDO
     261                 
     262               ENDIF
     263               
     264            ENDIF
     265!     c     
     266            IF( ifiltre.EQ.2 )  THEN
     267               
     268c     !-------------------------------------!
     269c     ! Dés-agregation des niveau verticaux !
     270c     ! uniquement necessaire pour une      !
     271c     ! execution OpenMP                    !
     272c     !-------------------------------------!
     273               ll_nb = 0
     274c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     275               DO l = 1, nbniv
     276                  ll_nb = ll_nb + 1
     277                  DO j = jdfil,jffil
     278                     DO i = 1, iim
     279                        champ( i,j,l ) = (champ_loc(i,j,ll_nb)
     280     &                       + champ_fft(i,j-jdfil+1,ll_nb))
     281     &                       * sdd12(i,sdd2_type)
     282                     ENDDO
     283                  ENDDO
     284               ENDDO
     285c$OMP END DO NOWAIT
     286               
     287            ELSE
     288               
     289               ll_nb = 0
     290c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     291               DO l = 1, nbniv_loc
     292                  ll_nb = ll_nb + 1
     293                  DO j = jdfil,jffil
     294                     DO i = 1, iim
     295                        champ( i,j,l ) = (champ_loc(i,j,ll_nb)
     296     &                       - champ_fft(i,j-jdfil+1,ll_nb))
     297     &                       * sdd12(i,sdd2_type)
     298                     ENDDO
     299                  ENDDO
     300               ENDDO
     301c$OMP END DO NOWAIT
     302               
     303            ENDIF
     304           
     305c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     306            DO l = 1, nbniv
     307               DO j = jdfil,jffil
     308                  champ( iip1,j,l ) = champ( 1,j,l )
     309               ENDDO
     310            ENDDO
     311c$OMP END DO NOWAIT
     312           
    282313ccccccccccccccccccccccccccccccccccccccccccccc
    283314c Utilisation du filtre FFT
    284315ccccccccccccccccccccccccccccccccccccccccccccc
    285316       
    286        ELSE
     317         ELSE
    287318       
    288319c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    289           DO l=1,nbniv
    290             DO j=jdfil,jffil
    291               DO  i = 1, iim
    292                 champ( i,j,l)= champ(i,j,l)*sdd1(i)
    293                 champ_fft( i,j,l) = champ(i,j,l)
    294               ENDDO
     320            DO l=1,nbniv
     321               DO j=jdfil,jffil
     322                  DO  i = 1, iim
     323                     champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
     324                     champ_fft( i,j,l) = champ(i,j,l)
     325                  ENDDO
     326               ENDDO
    295327            ENDDO
    296           ENDDO
    297328c$OMP END DO NOWAIT
    298329
    299       IF (jdfil<=jffil) THEN
    300         IF( ifiltre. EQ. -2 )   THEN
    301           CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
    302         ELSE IF ( griscal )     THEN
    303           CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
    304         ELSE
    305           CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
    306         ENDIF
    307       ENDIF
    308 
    309 
    310         IF( ifiltre.EQ. 2 )  THEN
     330            IF (jdfil<=jffil) THEN
     331               IF( ifiltre. EQ. -2 )   THEN
     332                  CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
     333               ELSE IF ( griscal )     THEN
     334                  CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
     335               ELSE
     336                  CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
     337               ENDIF
     338            ENDIF
     339
     340
     341            IF( ifiltre.EQ. 2 )  THEN
    311342c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
    312           DO l=1,nbniv
    313             DO j=jdfil,jffil
    314               DO  i = 1, iim
    315                 champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
    316      .                             *sdd2(i)
    317               ENDDO
    318             ENDDO
    319           ENDDO
     343               DO l=1,nbniv
     344                  DO j=jdfil,jffil
     345                     DO  i = 1, iim
     346                        champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
     347     &                       *sdd12(i,sdd2_type)
     348                     ENDDO
     349                  ENDDO
     350               ENDDO
    320351c$OMP END DO NOWAIT       
    321         ELSE
     352            ELSE
    322353       
    323354c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
    324           DO l=1,nbniv
    325             DO j=jdfil,jffil
    326               DO  i = 1, iim
    327                 champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
    328      .                            *sdd2(i)
    329               ENDDO
     355               DO l=1,nbniv
     356                  DO j=jdfil,jffil
     357                     DO  i = 1, iim
     358                        champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
     359     &                       *sdd12(i,sdd2_type)
     360                     ENDDO
     361                  ENDDO
     362               ENDDO
     363c$OMP END DO NOWAIT         
     364            ENDIF
     365c
     366c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     367            DO l=1,nbniv
     368               DO j=jdfil,jffil
     369!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
     370                  champ( iip1,j,l ) = champ( 1,j,l )
     371               ENDDO
    330372            ENDDO
    331           ENDDO
    332 c$OMP END DO NOWAIT         
    333         ENDIF
    334 c
    335 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    336         DO l=1,nbniv
    337           DO j=jdfil,jffil
    338 !            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
    339             champ( iip1,j,l ) = champ( 1,j,l )
    340           ENDDO
    341         ENDDO
    342373c$OMP END DO NOWAIT             
    343       ENDIF
     374         ENDIF
    344375c Fin de la zone de filtrage
    345376
    346377       
    347  100  CONTINUE
     378      ENDDO
    348379
    349380!      DO j=1,nlat
    … …  
    359390     
    360391c
    361 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
    362      *filtrer, sur la grille des scalaires'/)
    363 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
    364      *ltrer, sur la grille de V ou de Z'/)
     392 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
     393     &     filtrer, sur la grille des scalaires'/)
     394 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
     395     &     ltrer, sur la grille de V ou de Z'/)
    365396c$OMP MASTER     
    366397      CALL stop_timer
Note: See TracChangeset for help on using the changeset viewer.