← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

Viewing changes to vendor/IREGI/src/ti_reduce.f90

  • Committer: olivier Mattelaer
  • Date: 2015-03-05 00:14:16 UTC
  • mfrom: (258.1.9 2.3)
  • mto: (258.8.1 2.3)
  • mto: This revision was merged to the branch mainline in revision 259.
  • Revision ID: olivier.mattelaer@uclouvain.be-20150305001416-y9mzeykfzwnl9t0j
partial merge

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/IREGI/src/ti_reduce.f90
25
25
    LOGICAL,INTENT(OUT)::TSTABLE
26
26
    INTEGER,DIMENSION(1,1)::sol11
27
27
    REAL(KIND(1d0)),DIMENSION(1)::factor1
28
 
    INTEGER::first=0,i,j,numzerp,n,r,nr,init,ntot
 
28
    INTEGER::first=0,i,j,jk,numzerp,n,r,nr,init,ntot
29
29
    INTEGER,DIMENSION(NLOOPLINE)::zerp
30
 
    REAL(KIND(1d0)),DIMENSION(xiarraymax2)::syfactor
31
 
    INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE)::sy
 
30
    ! f[i,n]=(i+n-1)!/(n-1)!/i!
 
31
    ! nmax=5,rmax=6,NLOOPLINE=nmax+1
 
32
    ! see syntensor in ti_reduce.f90
 
33
    ! xiarraymax2=(f[0,nmax])+(f[1,nmax])+(f[2,nmax]+f[0,nmax])
 
34
    ! +(f[3,nmax]+f[1,nmax])+(f[4,nmax]+f[2,nmax]+f[0,nmax])
 
35
    ! +(f[5,nmax]+f[3,nmax]+f[1,nmax])+(f[6,nmax]+f[4,nmax]+f[2,nmax]+f[0,nmax])
 
36
    ! when rmax=5,nmax=5 -> xiarraymax2=314
 
37
    ! when rmax=6,nmax=5 -> xiarraymax2=610
 
38
    INTEGER::xiarraymax2,xiarraymax,xiarraymax3
 
39
    !REAL(KIND(1d0)),DIMENSION(xiarraymax2)::syfactor
 
40
    REAL(KIND(1d0)),DIMENSION(:),ALLOCATABLE::syfactor 
 
41
    !INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE)::sy
 
42
    INTEGER,DIMENSION(:,:),ALLOCATABLE::sy 
32
43
!    REAL(KIND(1d0)),PARAMETER::EPS=1d-10
33
44
    REAL(KIND(1d0))::temp
34
45
    INTEGER::ctmode_local
35
46
    LOGICAL::do_shift
36
 
    SAVE first
 
47
    SAVE first,xiarraymax,xiarraymax2,xiarraymax3,syfactor,sy
37
48
    IF(PRESENT(CTMODE))THEN
38
49
       ctmode_local=CTMODE
39
50
    ELSE
42
53
    MU_R_IREGI=MU
43
54
    IF(first.EQ.0)THEN
44
55
       IF(.NOT.print_banner)THEN
45
 
       WRITE(*,*)"#####################################################################################"
46
 
       WRITE(*,*)"#                                                                                   #"
47
 
       WRITE(*,*)"#                           IREGI-alpha-1.0.0                                       #"
48
 
       WRITE(*,*)"#    package for one-loop tensor Integral REduction with General propagator Indices #"
49
 
       WRITE(*,*)"#    Author: Hua-Sheng Shao (huasheng.shao@cern.ch/erdissshaw@gmail.com)            #"
50
 
       WRITE(*,*)"#    a) Physics School, Peking University, Beijing, China                           #"
51
 
       WRITE(*,*)"#    b) PH Department, TH Unit, CERN, Geneva, Switzerland                           #"
52
 
       WRITE(*,*)"#                                                                                   #"
53
 
       WRITE(*,*)"#####################################################################################"
54
 
       print_banner=.TRUE.
 
56
          INCLUDE "banner.inc"
 
57
          print_banner=.TRUE.
 
58
       ENDIF
 
59
       xiarraymax2=0
 
60
       DO i=0,MAXRANK_IREGI
 
61
          j=i/2+1
 
62
          DO jk=1,j
 
63
             xiarraymax2=xiarraymax2+&
 
64
                  xiarray_arg1(2*jk-2+MOD(i,2),MAXNLOOP_IREGI-1)
 
65
          ENDDO
 
66
       ENDDO
 
67
       xiarraymax=xiarray_arg1(MAXRANK_IREGI,MAXNLOOP_IREGI-1)
 
68
       xiarraymax3=xiarray_arg1(MAXRANK_IREGI,4)
 
69
       IF(.NOT.ALLOCATED(syfactor))THEN
 
70
          ALLOCATE(syfactor(xiarraymax2))
 
71
       ENDIF
 
72
       IF(.NOT.ALLOCATED(sy))THEN
 
73
          ALLOCATE(sy(xiarraymax2,-1:MAXNLOOP_IREGI))
55
74
       ENDIF
56
75
       ! initialization xiarray and metric
57
76
       CALL all_Integers(1,1,1,sol11,factor1)
98
117
       ENDIF
99
118
    ENDDO
100
119
    n=NLOOPLINE-numzerp
101
 
    IF(n.GE.6.OR.MAXRANK.GT.6)THEN
102
 
       WRITE(*,*)"ERROR: out of range of tensor_integral_reduce (N<7,R<7)"
 
120
    IF(n.GE.MAXNLOOP_IREGI.OR.MAXRANK.GT.MAXRANK_IREGI)THEN
 
121
       WRITE(*,100)"ERROR: out of range of tensor_integral_reduce (N<=",MAXNLOOP_IREGI,",R<=",MAXRANK_IREGI,")"
103
122
       STOP
104
123
    ENDIF
105
 
    CALL sytensor(n,MAXRANK,ntot,sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2))
106
 
    CALL symmetry(IMODE,NLOOPLINE,ntot,sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2),&
 
124
    CALL sytensor(n,MAXRANK,xiarraymax,xiarraymax2,ntot,&
 
125
         sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2))
 
126
    CALL symmetry(IMODE,NLOOPLINE,ntot,xiarraymax,xiarraymax2,xiarraymax3,&
 
127
         sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2),&
107
128
         numzerp,zerp,PDEN,M2L,NCOEFS,TICOEFS)
108
129
    TSTABLE=STABLE_IREGI
109
130
    IF(TSTABLE.AND.do_shift)THEN
113
134
       CALL SHIFT_MOM(PDEN_N,NCOEFS,TICOEFS_SAVE,TICOEFS)
114
135
    ENDIF
115
136
    RETURN
 
137
100 FORMAT(2X,A50,I2,A4,I2,A1)
116
138
  END SUBROUTINE tensor_integral_reduce
117
139
 
118
140
  SUBROUTINE tensor_integral_reduce2(IMODE,NLOOPLINE,MAXRANK,NCOEFS,PDEN,PijMatrix,M2L,MU,TICOEFS,TSTABLE)
129
151
    LOGICAL,INTENT(OUT)::TSTABLE
130
152
    INTEGER,DIMENSION(1,1)::sol11
131
153
    REAL(KIND(1d0)),DIMENSION(1)::factor1
132
 
    INTEGER::first=0,i,j,numzerp,n,r,nr,init,ntot
 
154
    INTEGER::first=0,i,j,jk,numzerp,n,r,nr,init,ntot
133
155
    INTEGER,DIMENSION(NLOOPLINE)::zerp
134
 
    REAL(KIND(1d0)),DIMENSION(xiarraymax2)::syfactor
135
 
    INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE)::sy
 
156
    ! f[i,n]=(i+n-1)!/(n-1)!/i!
 
157
    ! nmax=5,rmax=6,NLOOPLINE=nmax+1
 
158
    ! see syntensor in ti_reduce.f90 
 
159
    ! xiarraymax2=(f[0,nmax])+(f[1,nmax])+(f[2,nmax]+f[0,nmax])
 
160
    ! +(f[3,nmax]+f[1,nmax])+(f[4,nmax]+f[2,nmax]+f[0,nmax])
 
161
    ! +(f[5,nmax]+f[3,nmax]+f[1,nmax])+(f[6,nmax]+f[4,nmax]+f[2,nmax]+f[0,nmax])
 
162
    ! when rmax=5,nmax=5 -> xiarraymax2=314
 
163
    ! when rmax=6,nmax=5 -> xiarraymax2=610 
 
164
    INTEGER::xiarraymax2,xiarraymax,xiarraymax3
 
165
    !REAL(KIND(1d0)),DIMENSION(xiarraymax2)::syfactor
 
166
    REAL(KIND(1d0)),DIMENSION(:),ALLOCATABLE::syfactor
 
167
    !INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE)::sy
 
168
    INTEGER,DIMENSION(:,:),ALLOCATABLE::sy
136
169
    REAL(KIND(1d0))::temp
137
 
    SAVE first
 
170
    SAVE first,xiarraymax,xiarraymax2,xiarraymax3,syfactor,sy
138
171
    MU_R_IREGI=MU
139
172
    IF(first.EQ.0)THEN
140
173
       IF(.NOT.print_banner)THEN
141
 
       WRITE(*,*)"#####################################################################################"
142
 
       WRITE(*,*)"#                                                                                   #"
143
 
       WRITE(*,*)"#                           IREGI-alpha-1.0.0                                       #"
144
 
       WRITE(*,*)"#    package for one-loop tensor Integral REduction with General propagator Indices #"
145
 
       WRITE(*,*)"#    Author: Hua-Sheng Shao (huasheng.shao@cern.ch/erdissshaw@gmail.com)            #"
146
 
       WRITE(*,*)"#    a) Physics School, Peking University, Beijing, China                           #"
147
 
       WRITE(*,*)"#    b) PH Department, TH Unit, CERN, Geneva, Switzerland                           #"
148
 
       WRITE(*,*)"#                                                                                   #"
149
 
       WRITE(*,*)"#####################################################################################"
150
 
       print_banner=.TRUE.
 
174
          INCLUDE "banner.inc"
 
175
          print_banner=.TRUE.
 
176
       ENDIF
 
177
       xiarraymax2=0
 
178
       DO i=0,MAXRANK_IREGI
 
179
          j=i/2+1
 
180
          DO jk=1,j
 
181
             xiarraymax2=xiarraymax2+&
 
182
                  xiarray_arg1(2*jk-2+MOD(i,2),MAXNLOOP_IREGI-1)
 
183
          ENDDO
 
184
       ENDDO
 
185
       xiarraymax=xiarray_arg1(MAXRANK_IREGI,MAXNLOOP_IREGI-1)
 
186
       xiarraymax3=xiarray_arg1(MAXRANK_IREGI,4)
 
187
       IF(.NOT.ALLOCATED(syfactor))THEN
 
188
          ALLOCATE(syfactor(xiarraymax2))
 
189
       ENDIF
 
190
       IF(.NOT.ALLOCATED(sy))THEN
 
191
          ALLOCATE(sy(xiarraymax2,-1:MAXNLOOP_IREGI))
151
192
       ENDIF
152
193
       ! initialization xiarray and metric 
153
194
       CALL all_Integers(1,1,1,sol11,factor1)
175
216
       ENDIF
176
217
    ENDDO
177
218
    n=NLOOPLINE-numzerp
178
 
    IF(n.GE.6.OR.MAXRANK.GT.6)THEN
179
 
       WRITE(*,*)"ERROR: out of range of tensor_integral_reduce2 (N<7,R<7)"
 
219
    IF(n.GE.MAXNLOOP_IREGI.OR.MAXRANK.GT.MAXRANK_IREGI)THEN
 
220
       WRITE(*,100)"ERROR: out of range of tensor_integral_reduce2 (N<=",MAXNLOOP_IREGI,",R<=",MAXRANK_IREGI,")"
180
221
       STOP
181
222
    ENDIF
182
 
    CALL sytensor(n,MAXRANK,ntot,sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2))
183
 
    CALL symmetry2(IMODE,NLOOPLINE,ntot,sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2),&
 
223
    CALL sytensor(n,MAXRANK,xiarraymax,xiarraymax2,ntot,&
 
224
         sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2))
 
225
    CALL symmetry2(IMODE,NLOOPLINE,ntot,xiarraymax,xiarraymax2,xiarraymax3,&
 
226
         sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2),&
184
227
         numzerp,zerp,PDEN,PijMatrix,M2L,NCOEFS,TICOEFS)
185
228
    TSTABLE=STABLE_IREGI
186
229
    RETURN
 
230
100 FORMAT(2X,A51,I2,A4,I2,A1)
187
231
  END SUBROUTINE tensor_integral_reduce2
188
232
 
189
 
  SUBROUTINE sytensor(n,rmax,ntot,sy,syfactor)
 
233
  SUBROUTINE sytensor(n,rmax,xiarraymax,xiarraymax2,ntot,sy,syfactor)
190
234
    IMPLICIT NONE
191
 
    INTEGER,INTENT(IN)::n,rmax
 
235
    INTEGER,INTENT(IN)::n,rmax,xiarraymax,xiarraymax2
192
236
    INTEGER,INTENT(OUT)::ntot
193
237
    INTEGER,DIMENSION(xiarraymax2,-1:n),INTENT(OUT)::sy
 
238
    !INTEGER,DIMENSION(*,*),INTENT(OUT)::sy
194
239
    INTEGER::r,i,i0,i0max,rm2x0,nntot,init
195
240
    INTEGER,DIMENSION(xiarraymax,n)::sol
196
241
    REAL(KIND(1d0)),DIMENSION(xiarraymax2),INTENT(OUT)::syfactor
233
278
    RETURN
234
279
  END SUBROUTINE sytensor
235
280
 
236
 
  SUBROUTINE symmetry2(IMODE,NLOOPLINE,ntot,sy,syfactor,numzerp,zerp,PDEN,PijMatrix,M2L,NCOEFS,coefs)
 
281
  SUBROUTINE symmetry2(IMODE,NLOOPLINE,ntot,xiarraymax,xiarraymax2,xiarraymax3,&
 
282
       sy,syfactor,numzerp,zerp,PDEN,PijMatrix,M2L,NCOEFS,coefs)
237
283
    ! IMODE=0, IBP reduction
238
284
    ! IMODE=1, PaVe reduction
239
285
    ! IMODE=2, PaVe reduction with stablility improved by IBP reduction
240
286
    IMPLICIT NONE
241
 
    INTEGER,INTENT(IN)::ntot,numzerp,NLOOPLINE,IMODE,NCOEFS
 
287
    INTEGER,INTENT(IN)::ntot,numzerp,NLOOPLINE,IMODE,NCOEFS,xiarraymax,xiarraymax2,xiarraymax3
242
288
    INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE-numzerp),INTENT(IN)::sy
 
289
    !INTEGER,DIMENSION(*,*),INTENT(IN)::sy
243
290
    REAL(KIND(1d0)),DIMENSION(xiarraymax2),INTENT(IN)::syfactor
244
291
    INTEGER,DIMENSION(NLOOPLINE),INTENT(IN)::zerp
245
292
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3),INTENT(IN)::PDEN
255
302
    INTEGER::i,j,k,init,r,nntot,nt
256
303
    LOGICAL::lzero
257
304
    INTEGER,DIMENSION(0:NLOOPLINE)::nj
258
 
    INTEGER,DIMENSION(0:NLOOPLINE,1:5)::jlist
 
305
    INTEGER,DIMENSION(0:NLOOPLINE,1:MAXRANK_IREGI)::jlist
259
306
    REAL(KIND(1d0))::res
260
307
    COMPLEX(KIND(1d0)),DIMENSION(1:4)::scalar
261
308
    REAL(KIND(1d0)),DIMENSION(xiarraymax3)::coco
262
 
    INTEGER,DIMENSION(5)::lor
 
309
    !REAL(KIND(1d0)),DIMENSION(xiarraymax)::coco
 
310
    INTEGER,DIMENSION(MAXRANK_IREGI)::lor
263
311
    coefs(0:NCOEFS-1,1:4)=DCMPLX(0d0)
264
312
    j=1
265
313
    DO i=1,NLOOPLINE
295
343
             init=init+sol(j,4)
296
344
          ENDIF
297
345
          nj(0:nt)=0
298
 
          jlist(0:nt,1:5)=0
 
346
          jlist(0:nt,1:MAXRANK_IREGI)=0
299
347
          CALL recursive_symmetry(nt,1,r,PCLL(1:nt,0:3),sy(i,-1:nt),lor(1:r),&
300
 
               nj(0:nt),jlist(0:nt,1:5),res)
 
348
               nj(0:nt),jlist(0:nt,1:MAXRANK_IREGI),res)
301
349
          lzero=lzero.AND.(ABS(res).LT.EPS)
302
350
          IF(.NOT.ML5_CONVENTION)THEN
303
351
             coco(j)=res*factor(j)
344
392
                ! optimized PaVe reduction
345
393
                ! the stability has been improved by IBP reduction
346
394
                scalar(1:4)=pave_opt_reduce2(NLOOPLINE,paveindices,PijMatrix,M2L)                
347
 
             ENDIF
 
395
             ENDIF             
348
396
             IF(.NOT.STABLE_IREGI)THEN
349
397
                WRITE(*,*)"IREGI:WARNING, it detects unstable case, some integrals may set to be 0."
350
398
                RETURN
359
407
    RETURN
360
408
  END SUBROUTINE symmetry2
361
409
 
362
 
  SUBROUTINE symmetry(IMODE,NLOOPLINE,ntot,sy,syfactor,numzerp,zerp,PDEN,M2L,NCOEFS,coefs)
 
410
  SUBROUTINE symmetry(IMODE,NLOOPLINE,ntot,xiarraymax,xiarraymax2,xiarraymax3,&
 
411
       sy,syfactor,numzerp,zerp,PDEN,M2L,NCOEFS,coefs)
363
412
    ! IMODE=0, IBP reduction  
364
413
    ! IMODE=1, PaVe reduction
365
414
    ! IMODE=2, PaVe reduction with stablility improved by IBP reduction 
366
415
    IMPLICIT NONE
367
 
    INTEGER,INTENT(IN)::ntot,numzerp,NLOOPLINE,IMODE,NCOEFS
 
416
    INTEGER,INTENT(IN)::ntot,numzerp,NLOOPLINE,IMODE,NCOEFS,xiarraymax,xiarraymax2,xiarraymax3
368
417
    INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE-numzerp),INTENT(IN)::sy
 
418
    !INTEGER,DIMENSION(*,*),INTENT(IN)::sy 
369
419
    REAL(KIND(1d0)),DIMENSION(xiarraymax2),INTENT(IN)::syfactor
370
420
    INTEGER,DIMENSION(NLOOPLINE),INTENT(IN)::zerp
371
421
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3),INTENT(IN)::PDEN
380
430
    INTEGER::i,j,k,init,r,nntot,nt
381
431
    LOGICAL::lzero
382
432
    INTEGER,DIMENSION(0:NLOOPLINE)::nj
383
 
    INTEGER,DIMENSION(0:NLOOPLINE,1:5)::jlist
 
433
    INTEGER,DIMENSION(0:NLOOPLINE,1:MAXRANK_IREGI)::jlist
384
434
    REAL(KIND(1d0))::res
385
435
    COMPLEX(KIND(1d0)),DIMENSION(1:4)::scalar
386
436
    REAL(KIND(1d0)),DIMENSION(xiarraymax3)::coco
387
 
    INTEGER,DIMENSION(5)::lor
 
437
    !REAL(KIND(1d0)),DIMENSION(84)::coco
 
438
    INTEGER,DIMENSION(MAXRANK_IREGI)::lor
388
439
    coefs(0:NCOEFS-1,1:4)=DCMPLX(0d0)
389
440
    j=1
390
441
    DO i=1,NLOOPLINE
423
474
             init=init+sol(j,4)
424
475
          ENDIF
425
476
          nj(0:nt)=0
426
 
          jlist(0:nt,1:5)=0
 
477
          jlist(0:nt,1:MAXRANK_IREGI)=0
427
478
          CALL recursive_symmetry(nt,1,r,PCLL(1:nt,0:3),sy(i,-1:nt),lor(1:r),&
428
 
               nj(0:nt),jlist(0:nt,1:5),res)
 
479
               nj(0:nt),jlist(0:nt,1:MAXRANK_IREGI),res)
429
480
          lzero=lzero.AND.(ABS(res).LT.EPS)
430
481
          IF(.NOT.ML5_CONVENTION)THEN
431
482
             coco(j)=res*factor(j)
495
546
    REAL(KIND(1d0)),INTENT(INOUT)::res
496
547
    REAL(KIND(1d0))::temp,temp2
497
548
    INTEGER,DIMENSION(0:n),INTENT(INOUT)::njj
498
 
    INTEGER,DIMENSION(0:n,5),INTENT(INOUT)::jjlist
 
549
    INTEGER,DIMENSION(0:n,MAXRANK_IREGI),INTENT(INOUT)::jjlist
499
550
    REAL(KIND(1d0)),DIMENSION(n,0:3),INTENT(IN)::PCL
500
551
    INTEGER,DIMENSION(-1:n),INTENT(IN)::sy
501
552
    INTEGER,DIMENSION(r),INTENT(IN)::lor
575
626
    INTEGER,DIMENSION(nn),INTENT(IN)::lor2
576
627
    INTEGER,DIMENSION(nnj),INTENT(IN)::jlist
577
628
    INTEGER::ifirst=0,i,j,k,l,t
578
 
    INTEGER,PARAMETER::maxrank=2,maxrank2=16 ! 2*maxrank*4
579
 
    REAL(KIND(1d0)),DIMENSION(0:maxrank,0:maxrank,0:maxrank,0:maxrank)::gfsave
580
 
    INTEGER,DIMENSION(maxrank2)::lor
 
629
    INTEGER::maxgrank,maxgrank2 ! maxrank2=2*maxrank*4
 
630
    REAL(KIND(1d0)),DIMENSION(:,:,:,:),ALLOCATABLE::gfsave
 
631
    INTEGER,DIMENSION(:),ALLOCATABLE::lor
581
632
    INTEGER::n
582
633
    INTEGER,DIMENSION(0:3)::indexnum
583
634
    REAL(KIND(1d0))::gfunction
584
635
    SAVE ifirst,gfsave
585
636
    IF(ifirst.EQ.0)THEN
586
 
       DO i=0,maxrank
587
 
          DO j=0,maxrank
588
 
             DO k=0, maxrank
589
 
                DO l=0, maxrank
 
637
       maxgrank=MAXRANK_IREGI/2
 
638
       maxgrank2=2*maxgrank*4
 
639
       IF(.NOT.ALLOCATED(gfsave))THEN
 
640
          ALLOCATE(gfsave(0:maxgrank,0:maxgrank,0:maxgrank,0:maxgrank))
 
641
       ENDIF
 
642
       IF(.NOT.ALLOCATED(lor))THEN
 
643
          ALLOCATE(lor(maxgrank2))
 
644
       ENDIF
 
645
       DO i=0,maxgrank
 
646
          DO j=0,maxgrank
 
647
             DO k=0, maxgrank
 
648
                DO l=0, maxgrank
590
649
                   n=2*i+2*j+2*k+2*l
591
650
                   IF(n.EQ.0)THEN
592
651
                      gfsave(i,j,k,l)=1d0
644
703
    IMPLICIT NONE
645
704
    INTEGER,DIMENSION(0:3),INTENT(IN)::pos
646
705
    INTEGER::calc_pos
647
 
    INTEGER,PARAMETER::maxrank=5
 
706
    INTEGER,PARAMETER::maxrank=MAXRANK_IREGI
648
707
    INTEGER,DIMENSION(0:maxrank,0:maxrank,0:maxrank,0:maxrank)::possave
649
708
    INTEGER::ifirst=0,i,j,k,l,rk,init
650
709
    SAVE possave,ifirst