Imported Upstream version 4.8.1
[platform/upstream/gcc48.git] / 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