gcc/testsuite/gfortran.dg/pr41928.f90

   1 ! { dg-do compile }
   2 ! { dg-options "-O -fbounds-check -w" }
   3 MODULE kinds
   4   INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND ( 14, 200 )
   5   INTEGER, DIMENSION(:), ALLOCATABLE     :: nco,ncoset,nso,nsoset
   6   INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: co,coset
   7 END MODULE kinds
   8 MODULE ai_moments
   9   USE kinds
  10 CONTAINS
  11   SUBROUTINE cossin(la_max,npgfa,zeta,rpgfa,la_min,&
  12                     lb_max,npgfb,zetb,rpgfb,lb_min,&
  13                     rac,rbc,kvec,cosab,sinab)
  14     REAL(KIND=dp), DIMENSION(ncoset(la_max),&
  15       ncoset(lb_max))                        :: sc, ss
  16     DO ipgf=1,npgfa
  17       DO jpgf=1,npgfb
  18         IF (la_max > 0) THEN
  19           DO la=2,la_max
  20             DO ax=2,la
  21               DO ay=0,la-ax
  22                 sc(coset(ax,ay,az),1) = rap(1)*sc(coset(ax-1,ay,az),1) +&
  23                               f2 *          kvec(1)*ss(coset(ax-1,ay,az),1)
  24                 ss(coset(ax,ay,az),1) = rap(1)*ss(coset(ax-1,ay,az),1) +&
  25                               f2 *          kvec(1)*sc(coset(ax-1,ay,az),1)
  26               END DO
  27             END DO
  28           END DO
  29           IF (lb_max > 0) THEN
  30             DO lb=2,lb_max
  31               ss(1,coset(0,0,lb)) = rbp(3)*ss(1,coset(0,0,lb-1)) +&
  32                            f2 *         kvec(3)*sc(1,coset(0,0,lb-1))
  33               DO bx=2,lb
  34                 DO by=0,lb-bx
  35                   ss(1,coset(bx,by,bz)) = rbp(1)*ss(1,coset(bx-1,by,bz)) +&
  36                                f2 *           kvec(1)*sc(1,coset(bx-1,by,bz))
  37                 END DO
  38               END DO
  39             END DO
  40           END IF
  41         END IF
  42        DO j=ncoset(lb_min-1)+1,ncoset(lb_max)
  43         END DO
  44       END DO
  45     END DO
  46   END SUBROUTINE cossin
  47   SUBROUTINE moment(la_max,npgfa,zeta,rpgfa,la_min,&
  48                     lb_max,npgfb,zetb,rpgfb,&
  49                     lc_max,rac,rbc,mab)
  50     REAL(KIND=dp), DIMENSION(:), INTENT(IN)  :: zeta, rpgfa
  51     REAL(KIND=dp), DIMENSION(:), INTENT(IN)  :: zetb, rpgfb
  52     REAL(KIND=dp), DIMENSION(:, :, :), &
  53       INTENT(INOUT)                          :: mab
  54     REAL(KIND=dp), DIMENSION(3)              :: rab, rap, rbp, rpc
  55     REAL(KIND=dp), DIMENSION(ncoset(la_max),&
  56       ncoset(lb_max), ncoset(lc_max))        :: s
  57     DO ipgf=1,npgfa
  58       DO jpgf=1,npgfb
  59         IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
  60           DO k=1, ncoset(lc_max)-1
  61             DO j=nb+1,nb+ncoset(lb_max)
  62               DO i=na+1,na+ncoset(la_max)
  63                 mab(i,j,k) = 0.0_dp
  64               END DO
  65             END DO
  66           END DO
  67         END IF
  68         rpc = zetp*(zeta(ipgf)*rac+zetb(jpgf)*rbc)
  69         DO l=2, ncoset(lc_max)
  70           lx = indco(1,l)
  71           l2 = 0
  72           IF ( lz > 0 ) THEN
  73             IF ( lz > 1 ) l2 = coset(lx,ly,lz-2)
  74           ELSE IF ( ly > 0 ) THEN
  75             IF ( ly > 1 ) l2 = coset(lx,ly-2,lz)
  76             IF ( lx > 1 ) l2 = coset(lx-2,ly,lz)
  77           END IF
  78           s(1,1,l) = rpc(i)*s(1,1,l1)
  79           IF ( l2 > 0 ) s(1,1,l) = s(1,1,l) + f2*REAL(ni,dp)*s(1,1,l2)
  80         END DO
  81         DO l = 1, ncoset(lc_max)
  82           IF ( lx > 0 ) THEN
  83             lx1 = coset(lx-1,ly,lz)
  84           END IF
  85           IF ( ly > 0 ) THEN
  86             ly1 = coset(lx,ly-1,lz)
  87           END IF
  88           IF (la_max > 0) THEN
  89             DO la=2,la_max
  90               IF ( lz1 > 0 ) s(coset(0,0,la),1,l) = s(coset(0,0,la),1,l) + &
  91                              f2z*s(coset(0,0,la-1),1,lz1)
  92               IF ( ly1 > 0 ) s(coset(0,1,az),1,l) = s(coset(0,1,az),1,l) + &
  93                              f2y*s(coset(0,0,az),1,ly1)
  94               DO ay=2,la
  95                 s(coset(0,ay,az),1,l) = rap(2)*s(coset(0,ay-1,az),1,l) +&
  96                                        f2*REAL(ay-1,dp)*s(coset(0,ay-2,az),1,l)
  97                 IF ( ly1 > 0 ) s(coset(0,ay,az),1,l) = s(coset(0,ay,az),1,l) + &
  98                              f2y*s(coset(0,ay-1,az),1,ly1)
  99               END DO
 100               DO ay=0,la-1
 101                 IF ( lx1 > 0 ) s(coset(1,ay,az),1,l) = s(coset(1,ay,az),1,l) + &
 102                              f2x*s(coset(0,ay,az),1,lx1)
 103               END DO
 104               DO ax=2,la
 105                 DO ay=0,la-ax
 106                   s(coset(ax,ay,az),1,l) = rap(1)*s(coset(ax-1,ay,az),1,l) +&
 107                                           f3*s(coset(ax-2,ay,az),1,l)
 108                   IF ( lx1 > 0 ) s(coset(ax,ay,az),1,l) = s(coset(ax,ay,az),1,l) + &
 109                                  f2x*s(coset(ax-1,ay,az),1,lx1)
 110                 END DO
 111               END DO
 112             END DO
 113             IF (lb_max > 0) THEN
 114               DO j=2,ncoset(lb_max)
 115                 DO i=1,ncoset(la_max)
 116                   s(i,j,l) = 0.0_dp
 117                 END DO
 118               END DO
 119               DO la=la_start,la_max-1
 120                 DO ax=0,la
 121                   DO ay=0,la-ax
 122                     s(coset(ax,ay,az),2,l) = s(coset(ax+1,ay,az),1,l) -&
 123                                            rab(1)*s(coset(ax,ay,az),1,l)
 124                     s(coset(ax,ay,az),4,l) = s(coset(ax,ay,az+1),1,l) -&
 125                                            rab(3)*s(coset(ax,ay,az),1,l)
 126                   END DO
 127                 END DO
 128               END DO
 129               DO ax=0,la_max
 130                 DO ay=0,la_max-ax
 131                   IF (ax == 0) THEN
 132                     s(coset(ax,ay,az),2,l) = rbp(1)*s(coset(ax,ay,az),1,l)
 133                   ELSE
 134                     s(coset(ax,ay,az),2,l) = rbp(1)*s(coset(ax,ay,az),1,l) +&
 135                                             fx*s(coset(ax-1,ay,az),1,l)
 136                   END IF
 137                   IF (lx1 > 0) s(coset(ax,ay,az),2,l) = s(coset(ax,ay,az),2,l) +&
 138                         f2x*s(coset(ax,ay,az),1,lx1)
 139                   IF (ay == 0) THEN
 140                     s(coset(ax,ay,az),3,l) = rbp(2)*s(coset(ax,ay,az),1,l)
 141                   ELSE
 142                     s(coset(ax,ay,az),3,l) = rbp(2)*s(coset(ax,ay,az),1,l) +&
 143                                             fy*s(coset(ax,ay-1,az),1,l)
 144                   END IF
 145                   IF (ly1 > 0) s(coset(ax,ay,az),3,l) = s(coset(ax,ay,az),3,l) +&
 146                         f2y*s(coset(ax,ay,az),1,ly1)
 147                   IF (az == 0) THEN
 148                     s(coset(ax,ay,az),4,l) = rbp(3)*s(coset(ax,ay,az),1,l)
 149                   ELSE
 150                     s(coset(ax,ay,az),4,l) = rbp(3)*s(coset(ax,ay,az),1,l) +&
 151                                             fz*s(coset(ax,ay,az-1),1,l)
 152                   END IF
 153                   IF (lz1 > 0) s(coset(ax,ay,az),4,l) = s(coset(ax,ay,az),4,l) +&
 154                         f2z*s(coset(ax,ay,az),1,lz1)
 155                 END DO
 156               END DO
 157               DO lb=2,lb_max
 158                 DO la=la_start,la_max-1
 159                   DO ax=0,la
 160                     DO ay=0,la-ax
 161                       s(coset(ax,ay,az),coset(0,0,lb),l) =&
 162                         rab(3)*s(coset(ax,ay,az),coset(0,0,lb-1),l)
 163                       DO bx=1,lb
 164                         DO by=0,lb-bx
 165                           s(coset(ax,ay,az),coset(bx,by,bz),l) =&
 166                             rab(1)*s(coset(ax,ay,az),coset(bx-1,by,bz),l)
 167                         END DO
 168                       END DO
 169                     END DO
 170                   END DO
 171                 END DO
 172                 DO ax=0,la_max
 173                   DO ay=0,la_max-ax
 174                     IF (az == 0) THEN
 175                       s(coset(ax,ay,az),coset(0,0,lb),l) =&
 176                         rbp(3)*s(coset(ax,ay,az),coset(0,0,lb-1),l) +&
 177                         f3*s(coset(ax,ay,az),coset(0,0,lb-2),l)
 178                     END IF
 179                     IF (lz1 > 0) s(coset(ax,ay,az),coset(0,0,lb),l) =&
 180                                  f2z*s(coset(ax,ay,az),coset(0,0,lb-1),lz1)
 181                     IF (ay == 0) THEN
 182                       IF (ly1 > 0) s(coset(ax,ay,az),coset(0,1,bz),l) =&
 183                                  f2y*s(coset(ax,ay,az),coset(0,0,bz),ly1)
 184                       DO by=2,lb
 185                         s(coset(ax,ay,az),coset(0,by,bz),l) =&
 186                           f3*s(coset(ax,ay,az),coset(0,by-2,bz),l)
 187                         IF (ly1 > 0) s(coset(ax,ay,az),coset(0,by,bz),l) =&
 188                                  f2y*s(coset(ax,ay,az),coset(0,by-1,bz),ly1)
 189                       END DO
 190                       s(coset(ax,ay,az),coset(0,1,bz),l) =&
 191                         fy*s(coset(ax,ay-1,az),coset(0,0,bz),l)
 192                     END IF
 193                     IF (ax == 0) THEN
 194                       DO by=0,lb-1
 195                         IF (lx1 > 0) s(coset(ax,ay,az),coset(1,by,bz),l) =&
 196                                  f2x*s(coset(ax,ay,az),coset(0,by,bz),lx1)
 197                       END DO
 198                       DO bx=2,lb
 199                         DO by=0,lb-bx
 200                           s(coset(ax,ay,az),coset(bx,by,bz),l) =&
 201                             f3*s(coset(ax,ay,az),coset(bx-2,by,bz),l)
 202                           IF (lx1 > 0) s(coset(ax,ay,az),coset(bx,by,bz),l) =&
 203                                  f2x*s(coset(ax,ay,az),coset(bx-1,by,bz),lx1)
 204                         END DO
 205                       END DO
 206                       DO by=0,lb-1
 207                         IF (lx1 > 0) s(coset(ax,ay,az),coset(1,by,bz),l) =&
 208                                  f2x*s(coset(ax,ay,az),coset(0,by,bz),lx1)
 209                       END DO
 210                       DO bx=2,lb
 211                         DO by=0,lb-bx
 212                           s(coset(ax,ay,az),coset(bx,by,bz),l) =&
 213                             f3*s(coset(ax,ay,az),coset(bx-2,by,bz),l)
 214                           IF (lx1 > 0) s(coset(ax,ay,az),coset(bx,by,bz),l) =&
 215                                  f2x*s(coset(ax,ay,az),coset(bx-1,by,bz),lx1)
 216                         END DO
 217                       END DO
 218                     END IF
 219                   END DO
 220                 END DO
 221               END DO
 222             END IF
 223             IF (lb_max > 0) THEN
 224               DO lb=2,lb_max
 225                 IF (lz1 > 0) s(1,coset(0,0,lb),l) = s(1,coset(0,0,lb),l) +&
 226                              f2z*s(1,coset(0,0,lb-1),lz1)
 227                 IF (ly1 > 0) s(1,coset(0,1,bz),l) = s(1,coset(0,1,bz),l) +&
 228                              f2y*s(1,coset(0,0,bz),ly1)
 229               DO by=2,lb
 230                 s(1,coset(0,by,bz),l) = rbp(2)*s(1,coset(0,by-1,bz),l) +&
 231                                        f2*REAL(by-1,dp)*s(1,coset(0,by-2,bz),l)
 232                 IF (lx1 > 0) s(1,coset(1,by,bz),l) = s(1,coset(1,by,bz),l) +&
 233                              f2x*s(1,coset(0,by,bz),lx1)
 234               END DO
 235               DO bx=2,lb
 236                 DO by=0,lb-bx
 237                   IF (lx1 > 0) s(1,coset(bx,by,bz),l) = s(1,coset(bx,by,bz),l) +&
 238                                f2x*s(1,coset(bx-1,by,bz),lx1)
 239                 END DO
 240               END DO
 241             END DO
 242           END IF
 243         END IF
 244         END DO
 245         DO k=2,ncoset(lc_max)
 246           DO j=1,ncoset(lb_max)
 247           END DO
 248         END DO
 249       END DO
 250     END DO
 251   END SUBROUTINE moment
 252   SUBROUTINE diff_momop(la_max,npgfa,zeta,rpgfa,la_min,&
 253                     order,rac,rbc,difmab,mab_ext)
 254     REAL(KIND=dp), DIMENSION(:, :, :), &
 255       OPTIONAL, POINTER                      :: mab_ext
 256     REAL(KIND=dp), ALLOCATABLE, &
 257       DIMENSION(:, :, :)                     :: difmab_tmp
 258     DO imom = 1,ncoset(order)-1
 259       CALL adbdr(la_max,npgfa,rpgfa,la_min,&
 260                  difmab_tmp(:,:,2), difmab_tmp(:,:,3))
 261     END DO
 262   END SUBROUTINE diff_momop
 263 END MODULE ai_moments