source: PSPA/madxPSPA/libs/ptc/src/b_da_arrays_all.f90 @ 430

Last change on this file since 430 was 430, checked in by touze, 11 years ago

import madx-5.01.00
File size: 14.3 KB
Line 
1!The Full Polymorphic Package
2!Copyright (C) Etienne Forest and Frank Schmidt
3! See file a_scratch_size
4
5module da_arrays
6  use precision_constants
7  use my_own_1D_TPSA
8  use scratch_size
9  implicit none
10  public
11  integer lda,lea,lia,lst
12  integer,private,parameter::nmax=400
13  real(dp),private,parameter::tiny=1e-20_dp
14  private ludcmp_nr,lubksb_nr
15  ! johan
16  ! integer,parameter::lno=2,lnv=6,lnomax=8,lnvmax=9,lstmax=800500,ldamax=16000,leamax=5000,liamax=50000
17  !
18  integer,parameter::lno=200,lnv=100,lnomax=8,lnvmax=9,lstmax=800500,ldamax=16000,leamax=5000,liamax=50000
19  logical(lp) :: reallocate = .true.
20  logical(lp) :: notallocated = .true.
21  logical(lp),parameter::etiennefix=.true.
22
23
24  real(dp), allocatable,dimension(:)::cc                               ! lst
25  integer, allocatable,dimension(:)::i_1,i_2                             ! lst
26  integer, allocatable,dimension(:)::ie1,ie2,ieo                       !lea
27  integer, allocatable,dimension(:)::ia1,ia2                           ! 0:lia
28  integer, allocatable,dimension(:)::idano,idanv,idapo,idalm,idall     ! lda
29  character(10), allocatable,dimension(:)::daname                      ! daname(lda)*10
30  logical(lp), allocatable,dimension(:)::allvec                            !lda
31
32  integer nda_dab
33  integer :: ndamaxi=0
34  real(dp),TARGET :: total_da_size = 1.D38  !c_300
35  integer nst0,nomax,nvmax,nmmax,nocut,lfi
36  real(dp) facint(0:lno)
37  integer nhole
38  integer,TARGET :: lda_used =15000
39
40contains
41
42  subroutine alloc_all(no,nv)
43    implicit none
44    integer no,nv
45    if(reallocate) then
46       call dealloc_all
47       call danum0(no,nv)
48       call alloc_
49       notallocated=.false.
50    endif
51    if(notallocated) then
52       lia=liamax
53       lst=lstmax
54       lda=ldamax
55       lea=leamax
56       notallocated=.false.
57       call alloc_
58    endif
59  end subroutine alloc_all
60
61  subroutine alloc_
62    implicit none
63    integer i
64    allocate(cc(lst))
65    allocate(i_1(lst));allocate(i_2(lst));
66    allocate(ie1(lea));allocate(ie2(lea));allocate(ieo(lea));
67    allocate(ia1(0:lia));allocate(ia2(0:lia));
68    allocate(idano(lda));allocate(idanv(lda));allocate(idapo(lda));
69    allocate(idalm(lda));allocate(idall(lda));
70    allocate(daname(lda));
71    allocate(allvec(lda));
72    do i=1,lst
73       cc(I)=0.0_dp  ! ADDED BY ETIENNE
74       i_1(i)=0
75       i_2(i)=0
76    enddo
77    do i=1,lea
78       ie1(i)=0
79       ie2(i)=0
80       ieo(i)=0
81    enddo
82    do i=0,lia
83       ia1(i)=0
84       ia2(i)=0
85    enddo
86    do i=1,lda
87       idano(i)=0
88       idanv(i)=0
89       idapo(i)=0
90       idalm(i)=0
91       idall(i)=0
92    enddo
93  end subroutine alloc_
94
95
96  subroutine dealloc_all
97    implicit none
98    integer error,ipause,mypauses
99    IF (ALLOCATED(cc)) THEN
100       DEALLOCATE (cc, STAT = error)
101       IF(ERROR==0) THEN
102       ELSE
103          ipause=mypauses(100," cc ARRAY not DEALLOCATED : PROBLEMS")
104       ENDIF
105    ENDIF
106    IF (ALLOCATED(i_1)) THEN
107       DEALLOCATE (i_1, STAT = error)
108       IF(ERROR==0) THEN
109       ELSE
110          ipause=mypauses(101," i_1 ARRAY not DEALLOCATED : PROBLEMS")
111       ENDIF
112    ENDIF
113    IF (ALLOCATED(i_2)) THEN
114       DEALLOCATE (i_2, STAT = error)
115       IF(ERROR==0) THEN
116       ELSE
117          ipause=mypauses(102," i_2 ARRAY not DEALLOCATED : PROBLEMS")
118       ENDIF
119    ENDIF
120    IF (ALLOCATED(ie1)) THEN
121       DEALLOCATE (ie1, STAT = error)
122       IF(ERROR==0) THEN
123       ELSE
124          ipause=mypauses(103," ie1 ARRAY not DEALLOCATED : PROBLEMS")
125       ENDIF
126    ENDIF
127    IF (ALLOCATED(ie2)) THEN
128       DEALLOCATE (ie2, STAT = error)
129       IF(ERROR==0) THEN
130       ELSE
131          ipause=mypauses(104," ie2 ARRAY not DEALLOCATED : PROBLEMS")
132       ENDIF
133    ENDIF
134    IF (ALLOCATED(ieo)) THEN
135       DEALLOCATE (ieo, STAT = error)
136       IF(ERROR==0) THEN
137       ELSE
138          ipause=mypauses(105," ieo ARRAY not DEALLOCATED : PROBLEMS")
139       ENDIF
140    ENDIF
141    IF (ALLOCATED(ia1)) THEN
142       DEALLOCATE (ia1, STAT = error)
143       IF(ERROR==0) THEN
144       ELSE
145          ipause=mypauses(106," ia1 ARRAY not DEALLOCATED : PROBLEMS")
146       ENDIF
147    ENDIF
148    IF (ALLOCATED(ia2)) THEN
149       DEALLOCATE (ia2, STAT = error)
150       IF(ERROR==0) THEN
151       ELSE
152          ipause=mypauses(107," ia2 ARRAY not DEALLOCATED : PROBLEMS")
153       ENDIF
154    ENDIF
155    IF (ALLOCATED(idano)) THEN
156       DEALLOCATE (idano, STAT = error)
157       IF(ERROR==0) THEN
158       ELSE
159          ipause=mypauses(108," idano ARRAY not DEALLOCATED : PROBLEMS")
160       ENDIF
161    ENDIF
162    IF (ALLOCATED(idanv)) THEN
163       DEALLOCATE (idanv, STAT = error)
164       IF(ERROR==0) THEN
165       ELSE
166          ipause=mypauses(109," idanv ARRAY not DEALLOCATED : PROBLEMS")
167       ENDIF
168    ENDIF
169    IF (ALLOCATED(idapo)) THEN
170       DEALLOCATE (idapo, STAT = error)
171       IF(ERROR==0) THEN
172       ELSE
173          ipause=mypauses(110," idapo ARRAY not DEALLOCATED : PROBLEMS")
174       ENDIF
175    ENDIF
176    IF (ALLOCATED(idalm)) THEN
177       DEALLOCATE (idalm, STAT = error)
178       IF(ERROR==0) THEN
179       ELSE
180          ipause=mypauses(111," idalm ARRAY not DEALLOCATED : PROBLEMS")
181       ENDIF
182    ENDIF
183    IF (ALLOCATED(idall)) THEN
184       DEALLOCATE (idall, STAT = error)
185       IF(ERROR==0) THEN
186       ELSE
187          ipause=mypauses(112," idall ARRAY not DEALLOCATED : PROBLEMS")
188       ENDIF
189    ENDIF
190    IF (ALLOCATED(daname)) THEN
191       DEALLOCATE (daname, STAT = error)
192       IF(ERROR==0) THEN
193       ELSE
194          ipause=mypauses(112," daname ARRAY not DEALLOCATED : PROBLEMS")
195       ENDIF
196    ENDIF
197    IF (ALLOCATED(allvec)) THEN
198       DEALLOCATE (allvec, STAT = error)
199       IF(ERROR==0) THEN
200       ELSE
201          ipause=mypauses(112," allvec ARRAY not DEALLOCATED : PROBLEMS")
202       ENDIF
203    ENDIF
204  end subroutine dealloc_all
205
206  subroutine danum(no,nv,numda)
207    implicit none
208    integer i,mm,no,numda,nv
209    !     *****************************
210    !
211    !     THIS SUBROUTINE COMPUTES THE NUMBER OF MONOMIALS OF
212    !     ORDER NO AND NUMBER OF VARIABLES NV
213    !
214    numda = 1
215    mm = max(nv,no)
216    !
217    do i=1,min(nv,no)
218       numda = (numda*(mm+i))/i
219    enddo
220  end  subroutine danum
221
222  subroutine danum0(no,nv)
223    use scratch_size
224    use file_handler
225    implicit none
226    integer i,mm,no,nv,ldamin ,nvt
227    integer mf
228    real(dp) size
229    !     *****************************
230    !
231    !     THIS SUBROUTINE COMPUTES THE NUMBER OF MONOMIALS OF
232    !     ORDER NO AND NUMBER OF VARIABLES NV
233    !
234    lea = 1
235    if(etiennefix) then
236       nvt=nv
237    else
238       if(nv==1) then   ! still unknown feature of daini which must be taken into account
239          nvt=2           ! perhaps related to Ying Wu modification in daini
240       else
241          nvt=nv
242       endif
243    endif
244    mm = max(nvt,no)
245    do i=1,min(nvt,no)
246       lea = (lea*(mm+i))/i
247    enddo
248    ldamin=2+no + 2+ 4  !+ nd2t*(2+ndum)      ! (2+ no)daini  + varf +varc + assignmap(tpsalie.f90)
249
250
251    ldamin=ldamin+2        !  2 + sum of ndumuser  in assign tpsa.f90
252
253
254    lia=(no+1)**((nv+mod(nv,2))/2)
255    lda=lda_used+ldamin
256    lst=lda*lea
257
258    if(lingyun_yang) then
259       size=REAL(lst,kind=DP)*(8.0_dp/1024.0_dp**2+2.0_dp*4.0_dp/1024.0_dp**2)
260    else
261       size=(3.0_dp*REAL(lea,kind=DP)+2.0_dp*(REAL(lia,kind=DP)+1.0_dp)+5.0_dp*REAL(lda,kind=DP))*4.0_dp/1024.0_dp**2
262       size=size+10.0_dp*REAL(lda,kind=DP)/1024.0_dp**2+REAL(lda,kind=DP)/8.0_dp/1024.0_dp**2
263       size=size+REAL(lst,kind=DP)*(8.0_dp/1024.0_dp**2+2.0_dp*4.0_dp/1024.0_dp**2)
264    endif
265    if(size>total_da_size.or.printdainfo) then
266       w_p=0
267       w_p%nc=13
268       w_p%fc='(12(1X,A72,/),(1X,A72))'
269       write(w_p%c(1),'(a10,1x,i4,1x,i4)') " no,nv  = ",no,nv
270       write(w_p%c(2),'(a10,1x,i8)') "    LEA = ",lea
271       write(w_p%c(3),'(a24,1x,i8)') " ldamin (with nd2=6)  = ",ldamin
272       write(w_p%c(4),'(a8,1x,i8)') " lia  = ",lia
273       write(w_p%c(5),'(a8,1x,i8)') " lda  = ",lda
274       write(w_p%c(6),'(a8,1x,i8)') " lst  = ",lst
275       write(w_p%c(7),'(a14,1x,i8)') " ndamaxi    = ",ndamaxi
276       write(w_p%c(8),'(a18,1x,g21.14)') " size in Mbytes = ",size
277       write(w_p%c(9),'(a25,1x,g21.14)') " Total_da_size Allowed = ",Total_da_size
278       w_p%c(10)=" "
279       w_p%c(11)="************************************"
280       w_p%c(12)="* Execution Continues Nevertheless *"
281       w_p%c(13)="************************************"
282       ! call !write_e(1000)
283       call kanalnummer(mf)
284       open(unit=mf,file='too_big_da.txt')
285       write(mf,*) "no,nv  = ",no,nv
286       write(mf,*) "LEA = ",lea
287       write(mf,*) "ldamin (with nd2=6)  = ",ldamin
288       write(mf,*) "lia  = ",lia
289       write(mf,*) "lda  = ",lda
290       write(mf,*) "lst  = ",lst
291       write(mf,*) "ndamaxi    = ",ndamaxi
292       write(mf,*) "size in Mbytes = ",size
293       write(mf,*) "Total_da_size Allowed = ",Total_da_size
294       write(6,*) "no,nv  = ",no,nv
295       write(6,*) "LEA = ",lea
296       write(6,*) "ldamin (with nd2=6)  = ",ldamin
297       write(6,*) "lia  = ",lia
298       write(6,*) "lda  = ",lda
299       write(6,*) "lst  = ",lst
300       write(6,*) "ndamaxi    = ",ndamaxi
301       write(6,*) "size in Mbytes = ",size
302       write(6,*) "Total_da_size Allowed = ",Total_da_size
303       close(mf)
304    endif
305
306
307  end  subroutine danum0
308
309!!!!! some da stuff
310
311  real(dp) function bran(xran)
312    implicit none
313    !     ************************************
314    !
315    !     VERY SIMPLE RANDOM NUMBER GENERATOR
316    !
317    !-----------------------------------------------------------------------------
318    !
319    real(dp) xran
320    !
321    xran = xran + 10.0_dp
322    if(xran.gt.1e4_dp) xran = xran - 9999.12345e0_dp
323    bran = abs(sin(xran))
324    bran = 10*bran
325    bran = bran - int(bran)
326    !      IF(BRAN.LT. c_0_1) BRAN = BRAN + c_0_1
327    !
328    return
329  end function bran
330
331  subroutine matinv(a,ai,n,nmx,ier)
332
333    implicit none
334
335    integer i,ier,j,n,nmx
336    integer,dimension(nmax)::indx
337    real(dp) d
338    real(dp),dimension(nmx,nmx)::a,ai
339    real(dp),dimension(nmax,nmax)::aw
340    !
341    !    if((.not.C_%STABLE_DA)) then
342    !       if(c_%watch_user) then
343    !          write(6,*) "big problem in dabnew ", sqrt(crash)
344    !       endif
345    !       return
346    !    endif
347
348    aw(1:n,1:n) = a(1:n,1:n)
349
350    call ludcmp_nr(aw,n,nmax,indx,d,ier)
351    if (ier .eq. 132) return
352
353    ai(1:n,1:n) = 0.0_dp
354    !    forall (i = 1:n) ai(i,i) = one
355    do i=1,n
356       ai(i,i) = 1.0_dp
357    enddo
358
359    do j=1,n
360       call lubksb_nr(aw,n,nmax,indx,ai(1,j),nmx)
361    enddo
362
363  end subroutine matinv
364
365
366  !
367  subroutine ludcmp_nr(a,n,np,indx,d,ier)
368    implicit none
369    !     ************************************
370    !
371    !     THIS SUBROUTINE DECOMPOSES A MATRIX INTO LU FORMAT
372    !     INPUT A: NXN MATRIX - WILL BE OVERWRITTEN BY THE LU DECOMP.
373    !           NP: PHYSICAL DIMENSION OF A
374    !           INDX: ROW PERMUTATION VECTOR
375    !           D: EVEN OR ODD ROW INTERCHANGES
376    !
377    !     REFERENCE: NUMERICAL RECIPIES BY PRESS ET AL (CAMBRIDGE) PG. 35
378    !
379    !-----------------------------------------------------------------------------
380    !
381    integer i,ier,imax,j,k,n,np
382    integer,dimension(np)::indx
383    real(dp) aamax,d,dum,sum
384    real(dp),dimension(np,np)::a
385    real(dp),dimension(nmax)::vv
386    !
387    !    if((.not.C_%STABLE_DA)) then
388    !       if(c_%watch_user) then
389    !          write(6,*) "big problem in dabnew ", sqrt(crash)
390    !       endif
391    !       return
392    !    endif
393    ier=0
394    d=1.0_dp
395    do i=1,n
396       aamax=0.0_dp
397       do j=1,n
398          if(abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
399       enddo
400       if(aamax.eq.0.0_dp) then
401          ier=132
402          return
403       endif
404       vv(i)=1.0_dp/aamax
405    enddo
406    do j=1,n
407       if(j.gt.1) then
408          do i=1,j-1
409             sum=a(i,j)
410             if(i.gt.1) then
411                do k=1,i-1
412                   sum=sum-a(i,k)*a(k,j)
413                enddo
414                a(i,j)=sum
415             endif
416          enddo
417       endif
418       aamax=0.0_dp
419       do i=j,n
420          sum=a(i,j)
421          if (j.gt.1) then
422             do k=1,j-1
423                sum=sum-a(i,k)*a(k,j)
424             enddo
425             a(i,j)=sum
426          endif
427          dum=vv(i)*abs(sum)
428          if(dum.ge.aamax) then
429             imax=i
430             aamax=dum
431          endif
432       enddo
433       if (j.ne.imax) then
434          do k=1,n
435             dum=a(imax,k)
436             a(imax,k)=a(j,k)
437             a(j,k)=dum
438          enddo
439          d=-d
440          vv(imax)=vv(j)
441       endif
442       indx(j)=imax
443       if(j.ne.n) then
444          if(a(j,j).eq.0.0_dp) a(j,j)=tiny
445          dum=1.0_dp/a(j,j)
446          do i=j+1,n
447             a(i,j)=a(i,j)*dum
448          enddo
449       endif
450    enddo
451    if(a(n,n).eq.0.0_dp) a(n,n)=tiny
452    return
453  end subroutine ludcmp_nr
454  !
455  subroutine lubksb_nr(a,n,np,indx,b,nmx)
456    implicit none
457    !     ************************************
458    !
459    !     THIS SUBROUTINE SOLVES SET OF LINEAR EQUATIONS AX=B,
460    !     INPUT A: NXN MATRIX IN lu FORM GIVEN BY ludcmp_nr
461    !           NP: PHYSICAL DIMENSION OF A
462    !           INDX: ROW PERMUTATION VECTOR
463    !           D: EVEN OR ODD ROW INTERCHANGES
464    !           B: RHS OF LINEAR EQUATION - WILL BE OVERWRITTEN BY X
465    !
466    !     REFERENCE: NUMERICAL RECIPIES BY PRESS ET AL (CAMBRIDGE) PG. 36
467    !
468    !-----------------------------------------------------------------------------
469    !
470    integer i,ii,j,ll,n,nmx,np
471    integer,dimension(np)::indx
472    real(dp) sum
473    real(dp),dimension(np,np)::a
474    real(dp),dimension(nmx)::b
475    !
476    !    if((.not.C_%STABLE_DA)) then
477    !       if(c_%watch_user) then
478    !          write(6,*) "big problem in dabnew ", sqrt(crash)
479    !       endif
480    !       return
481    !    endif
482    ii = 0
483    do i=1,n
484       ll = indx(i)
485       sum = b(ll)
486       b(ll) = b(i)
487       if(ii.ne.0) then
488          do j=ii,i-1
489             sum = sum-a(i,j)*b(j)
490          enddo
491       else if (sum.ne.0.0_dp) then
492          ii = i
493       endif
494       b(i)=sum
495    enddo
496    do i=n,1,-1
497       sum=b(i)
498       if(i.lt.n) then
499          do j=i+1,n
500             sum = sum-a(i,j)*b(j)
501          enddo
502       endif
503
504       b(i)=sum/a(i,i)
505
506    enddo
507    return
508  end subroutine lubksb_nr
509  !
510
511
512end  module da_arrays
Note: See TracBrowser for help on using the repository browser.