source: PSPA/madxPSPA/src/c_tpsa_interface.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: 29.3 KB
Line 
1!/*
2! * Copyright(C) 2008 by Lingyun Yang
3! * lingyun(.dot.]yang@gmail.com
4! * http://www.lingyunyang.com
5! *
6! * Please get permission from Lingyun Yang before you redistribute this file.
7! *
8! * Version: $Id: c_tpsa_interface.F90,v 1.13 2009-06-24 09:06:43 frs Exp $
9! */
10
11
12module dabnew
13  !$$$$  use da_arrays
14  use dabnew_b !$$$$
15  implicit none
16  public
17  ! integer,private,parameter:: lsw=1
18
19  ! integer,private,parameter::nmax=400,lsw=1
20  ! real(dp),private,parameter::tiny=c_1d_20
21  character(120),private :: line
22
23#ifdef _WIN32_DLL
24  !DEC$ ATTRIBUTES DLLIMPORT :: ad_tra, ad_shift, ad_print, ad_save_block, ad_read_block
25  !DEC$ ATTRIBUTES DLLIMPORT :: ad_fill_ran, ad_nvar, ad_length, ad_derivative
26  !DEC$ ATTRIBUTES DLLIMPORT :: ad_subst, ad_cos, ad_sin, ad_log, ad_exp, ad_sqrt, ad_abs
27  !DEC$ ATTRIBUTES DLLIMPORT :: ad_div_c, ad_c_div, ad_mult_const, ad_add_const
28  !DEC$ ATTRIBUTES DLLIMPORT :: ad_div, ad_mult, ad_sub, ad_reset, ad_pok, ad_pek
29  !DEC$ ATTRIBUTES DLLIMPORT :: ad_var, ad_truncate, ad_const, ad_count, ad_free, ad_add
30  !DEC$ ATTRIBUTES DLLIMPORT :: ad_copy, ad_clean, ad_alloc, ad_reserve, ad_init, ad_elem
31
32  DLL_IMPORT ad_tra, ad_shift, ad_print, ad_save_block, ad_read_block
33  DLL_IMPORT ad_fill_ran, ad_nvar, ad_length, ad_derivative
34  DLL_IMPORT ad_subst, ad_cos, ad_sin, ad_log, ad_exp, ad_sqrt, ad_abs
35  DLL_IMPORT ad_div_c, ad_c_div, ad_mult_const, ad_add_const
36  DLL_IMPORT ad_div, ad_mult, ad_sub, ad_reset, ad_pok, ad_pek
37  DLL_IMPORT ad_var, ad_truncate, ad_const, ad_count, ad_free, ad_add
38  DLL_IMPORT ad_copy, ad_clean, ad_alloc, ad_reserve, ad_init, ad_elem ,ad_resetvars
39#endif
40
41  private trx_cpp
42
43
44contains
45
46  !done
47  subroutine daini(nd,nv,k)
48    implicit none
49    !    integer, intent(in) :: nd,nv,nd2,k
50    integer nd,nv,k,last_nv
51    if(lingyun_yang) then !%%%%
52       if(last_tpsa==1.and.lielib_print(10)==0) then
53          call ad_nvar(last_nv)
54          call ad_resetvars(last_nv)
55       endif
56       call danum0(nd,nv)
57       call ad_init(nv, nd)
58       call ad_reserve(lda)
59       !     write(6,*) " LDA Lingyun ", lda
60       last_tpsa=1
61    else !%%%%
62       call daini_b(nd,nv,k)
63       !           last_tpsa=2   done in daini_b
64    endif !%%%%
65  end subroutine daini
66
67  !done
68  subroutine daall0(i)
69    implicit none
70    !    integer, intent(in) :: i
71    integer  i
72    if(lingyun_yang) then !%%%%
73       call ad_alloc(i)
74    else !%%%%
75       call daall0_b(i)
76    endif !%%%%
77  end subroutine daall0
78
79  !done
80  subroutine daall1(i,ccc,n0,nv)
81    implicit none
82    !    integer, intent(in) :: i
83    integer i,n0, nv
84    character(10) ccc
85
86    if(lingyun_yang) then !%%%%
87       call ad_alloc(i)
88    else !%%%%
89       call daall1_b(i,ccc,n0,nv)
90    endif !%%%%
91
92  end subroutine daall1
93
94  subroutine daclean(ina,value)
95    implicit none
96    integer ina
97    !
98    real(dp) value
99    !
100    !
101    if(lingyun_yang) then !%%%%
102       call ad_clean(ina, value)
103    else !%%%%
104       call daclean_b(ina, value)
105    endif !%%%%
106    return
107  end subroutine daclean
108
109  !done
110  subroutine daadd(i,j,k)
111    implicit none
112    !    integer, intent(in) :: i,j,k
113    integer  i,j,k
114    integer itmp
115    !write(*,*) k
116    ! if j!= k and i != k, then optimize
117    if(lingyun_yang) then !%%%%
118       call ad_alloc(itmp)
119       call ad_copy(i, itmp)
120       call ad_add(itmp, j)
121       call ad_copy(itmp, k)
122       call ad_free(itmp)
123    else !%%%%
124       call daadd_b(i,j,k)
125    endif !%%%%
126  end subroutine daadd
127
128  !done
129  subroutine dacon(i,r)
130    implicit none
131    !    integer, intent(in) :: i
132    !    real(dp), intent(in) :: r
133    integer  i
134    real(dp)  r
135
136    if(lingyun_yang) then !%%%%
137       call ad_const(i, r)
138    else !%%%%
139       call dacon_b(i,r)
140    endif !%%%%
141  end subroutine dacon
142
143  !done
144  subroutine dadal(idal,j)
145    implicit none
146    integer  j
147    !    integer, intent(in) :: j
148    integer,dimension(:)::idal
149    integer k
150    if(lingyun_yang) then !%%%%
151       do k=1,j
152          call ad_free(idal(k))
153       enddo
154    else !%%%%
155       call dadal_b(idal,j)
156    endif !%%%%
157  end subroutine dadal
158
159  !done
160  subroutine dadal1(idal)
161    implicit none
162    integer  idal
163    !    integer, intent(inout) :: idal
164    if(lingyun_yang) then !%%%%
165       call ad_free(idal)
166    else !%%%%
167       call dadal1_b(idal)
168    endif !%%%%
169  end subroutine dadal1
170
171  !done
172  subroutine count_da(idal)
173    implicit none
174    !    integer, intent(inout) :: idal
175    integer  idal
176    !call ad_all(i)
177    if(lingyun_yang) then !%%%%
178       call ad_count(idal)
179    else !%%%%
180       call count_da_b(idal)
181    endif !%%%%
182  end subroutine count_da
183
184  !done
185  subroutine davar(ina,ckon,i)
186    implicit none
187    !    integer, intent(in) :: ina,i
188    integer ina,i
189    real(dp), intent(in) :: ckon
190    if(lingyun_yang) then !%%%%
191       call ad_var(ina, ckon, i-1)
192    else !%%%%
193       call davar_b(ina,ckon,i)
194    endif !%%%%
195
196  end subroutine davar
197
198  ! done
199  subroutine danot(not)
200    implicit none
201    !    integer, intent(in) :: not
202    integer  not
203    if(lingyun_yang) then !%%%%
204       print *, 'ERROR: This is not used in new TPSA routines.'
205       STOP
206    else !%%%%
207       call danot_b(not)
208    endif !%%%%
209  end subroutine danot
210
211  !done
212  subroutine datrunc(ina, imd, inb)
213    implicit none
214    !    integer, intent(in) :: ina,imd,inb
215    integer  ina,imd,inb
216    if(lingyun_yang) then !%%%%
217       call ad_copy(ina, inb)
218       ! truncate imd and above.
219       call ad_truncate(inb, imd)
220    else !%%%%
221       call datrunc_b(ina, imd, inb)
222    endif !%%%%
223  end subroutine datrunc
224
225  ! done
226  subroutine daeps(deps)
227    implicit none
228    real(dp) deps
229    !    real(dp), intent(inout) :: deps
230
231    if(lingyun_yang) then !%%%%
232       print *, 'ERROR: We use machine dependent eps instead.'
233       if (deps<zero) then
234          deps = eps
235       else
236          STOP
237       endif
238    else !%%%%
239       call daeps_b(deps)
240    endif !%%%%
241
242  end subroutine daeps
243
244  ! done
245  subroutine dapek(ina,jv,cjj)
246    implicit none
247    !    integer,intent(in) :: ina
248    integer  ina
249    real(dp),intent(inout) :: cjj
250    integer,dimension(:)::jv     ! 2002.12.4
251    if(lingyun_yang) then !%%%%
252       call ad_pek(ina, jv, size(jv), cjj)
253    else !%%%%
254       call dapek_b(ina,jv,cjj)
255    endif !%%%%
256  end subroutine dapek
257
258  !done
259  subroutine dapok(ina,jv,cjj)
260    implicit none
261    integer  ina
262    !    integer,intent(in) :: ina
263    real(dp),intent(in) :: cjj
264    integer,dimension(:)::jv     ! 2002.12.4
265    if(lingyun_yang) then !%%%%
266       call ad_pok(ina, jv, size(jv), cjj)
267    else !%%%%
268       call dapok_b(ina,jv,cjj)
269    endif !%%%%
270  end subroutine dapok
271
272  !done
273  subroutine daclr(ina)
274    implicit none
275    !    integer,intent(in) :: ina
276    integer  ina
277    if(lingyun_yang) then !%%%%
278       call ad_reset(ina)
279    else !%%%%
280       call daclr_b(ina)
281    endif !%%%%
282  end subroutine daclr
283
284  !done
285  subroutine dacop(ina,inb)
286    implicit none
287    !    integer,intent(in) :: ina,inb
288    integer  ina,inb
289    if(lingyun_yang) then !%%%%
290       call ad_copy(ina, inb)
291    else !%%%%
292       call dacop_b(ina,inb)
293    endif !%%%%
294  end subroutine dacop
295
296  !done
297  subroutine dasub(ina,inb,inc)
298    implicit none
299    !    integer,intent(in) :: ina,inb,inc
300    integer  ina,inb,inc
301    integer i1,i2
302    if(lingyun_yang) then !%%%%
303       call ad_alloc(i1)
304       call ad_alloc(i2)
305       call ad_copy(ina, i1)
306       call ad_copy(inb, i2)
307       call ad_sub(i1, i2)
308       call ad_copy(i1, inc)
309       call ad_free(i1)
310       call ad_free(i2)
311    else !%%%%
312       call dasub_b(ina,inb,inc)
313    endif !%%%%
314  end subroutine dasub
315
316  !done
317  subroutine damul(ina,inb,inc)
318    implicit none
319    !    integer,intent(in) :: ina,inb,inc
320    integer  ina,inb,inc
321    integer i1, i2
322    if(lingyun_yang) then !%%%%
323       call ad_alloc(i1)
324       call ad_alloc(i2)
325       call ad_copy(ina, i1)
326       call ad_copy(inb, i2)
327       call ad_mult(i1, i2, inc)
328       call ad_free(i1)
329       call ad_free(i2)
330    else !%%%%
331       call damul_b(ina,inb,inc)
332    endif !%%%%
333  end subroutine damul
334
335  !done
336  subroutine dadiv(ina,inb,inc)
337    implicit none
338    !    integer,intent(in) :: ina,inb,inc
339    integer  ina,inb,inc
340    integer i1, i2
341
342    if(lingyun_yang) then !%%%%
343       call ad_alloc(i1)
344       call ad_alloc(i2)
345       call ad_copy(ina, i1)
346       call ad_copy(inb, i2)
347       call ad_div(i1, i2, inc)
348       call ad_free(i1)
349       call ad_free(i2)
350    else !%%%%
351       call dadiv_b(ina,inb,inc)
352    endif !%%%%
353
354  end subroutine dadiv
355
356  !done
357  subroutine dacad(ina,ckon,inc)
358    implicit none
359    integer ina,inc
360    !    integer,intent(in) :: ina,inc
361    real(dp), intent(in):: ckon
362    if(lingyun_yang) then !%%%%
363       call ad_copy(ina, inc)
364       call ad_add_const(inc, ckon)
365    else !%%%%
366       call dacad_b(ina,ckon,inc)
367    endif !%%%%
368  end subroutine dacad
369
370  !done
371  subroutine dacsu(ina,ckon,inc)
372    implicit none
373    !    integer,intent(in) :: ina,inc
374    integer  ina,inc
375    real(dp)  ckon
376    call dacad(ina, -ckon, inc)
377  end subroutine dacsu
378
379  !done
380  subroutine dasuc(ina,ckon,inc)
381    implicit none
382    integer  ina,inc
383    real(dp)  ckon
384    if(lingyun_yang) then !%%%%
385       call ad_copy(ina, inc)
386       call ad_mult_const(inc, -1.d0)
387       call dacad(inc, ckon, inc)
388    else !%%%%
389       call dasuc_b(ina,ckon,inc)
390    endif !%%%%
391  end subroutine dasuc
392
393  !done
394  subroutine dacmu(ina,ckon,inc)
395    implicit none
396    integer  ina,inc
397    integer it1
398    real(dp)  ckon
399    if(lingyun_yang) then !%%%%
400       call ad_alloc(it1)
401       call ad_copy(ina, it1)
402       call ad_mult_const(it1, ckon)
403       call ad_copy(it1, inc)
404       call ad_free(it1)
405    else !%%%%
406       call dacmu_b(ina,ckon,inc)
407    endif !%%%%
408  end subroutine dacmu
409
410  !done
411  subroutine dadic(ina,ckon,inc)
412    implicit none
413    integer ina,inc
414    real(dp)  ckon
415    integer i1
416
417    if(lingyun_yang) then !%%%%
418       call ad_alloc(i1)
419       call ad_c_div(ina, ckon, i1)
420       call ad_copy(i1, inc)
421       call ad_free(i1)
422    else !%%%%
423       call dadic_b(ina,ckon,inc)
424    endif !%%%%
425  end subroutine dadic
426
427  !done
428  subroutine dacdi(ina,ckon,inc)
429    implicit none
430    integer  ina,inc
431    real(dp)  ckon
432    integer i1
433
434    if(lingyun_yang) then !%%%%
435       call ad_alloc(i1)
436       call ad_copy(ina, i1)
437       call ad_div_c(i1, ckon)
438       call ad_copy(i1, inc)
439       call ad_free(i1)
440    else !%%%%
441       call dacdi_b(ina,ckon,inc)
442    endif !%%%%
443  end subroutine dacdi
444
445  !done
446  subroutine daabs(ina,r)
447    implicit none
448    integer  ina
449    real(dp)  r
450    if(lingyun_yang) then !%%%%
451       call ad_abs(ina, r)
452    else !%%%%
453       call daabs_b(ina,r)
454    endif !%%%%
455
456  end subroutine daabs
457
458  subroutine dalin(ina,afac,inb,bfac,inc)
459    implicit none
460    integer  ina,inb,inc
461    integer it1,it2
462    real(dp)  afac,bfac
463    !write(*,*) "Not implemented"
464    !call exit
465
466    if(lingyun_yang) then !%%%%
467       call ad_alloc(it1)
468       call ad_alloc(it2)
469       call ad_copy(ina,it1)
470       call ad_copy(inb,it2)
471       call ad_mult_const(it1,afac)
472       call ad_mult_const(it2,bfac)
473       call ad_copy(it1, inc)
474       call ad_add(inc, it2)
475       call ad_free(it1)
476       call ad_free(it2)
477    else !%%%%
478       call dalin_b(ina,afac,inb,bfac,inc)
479    endif !%%%%
480  end subroutine dalin
481
482  subroutine dafun(cf,ina,inc)
483    implicit none
484    integer  ina,inc
485    character(4)  cf
486    INTEGER IPAUSE,MYPAUSES
487    if(lingyun_yang) then !%%%%
488       select case (cf)
489       case('INV ')
490          !call exit
491          call ad_c_div(ina, 1.d0, inc)
492       case('SQRT')
493          call ad_sqrt(ina, inc)
494       case('EXP ')
495          call ad_exp(ina, inc)
496       case('LOG ')
497          call ad_log(ina, inc)
498       case('SIN ')
499          call ad_sin(ina, inc)
500       case('COS ')
501          call ad_cos(ina, inc)
502       case('SINH')
503          write(*,*) "BUG in sinh"
504          STOP
505       case('COSH')
506          write(*,*) "BUG in cosh"
507          STOP
508       case default
509          write(line,'(a28,1x,a4)')  'ERROR, UNSOPPORTED FUNCTION ',cf
510          ipause=mypauses(35,line)
511       end select
512       !
513    else !%%%%
514       call dafun_b(cf,ina,inc)
515    endif !%%%%
516  end subroutine dafun
517
518  subroutine dacct(ma,ia,mb,ib,mc,ic)
519    implicit none
520    integer,dimension(:)::ma,mb,mc
521    integer  ia,ib,ic
522    integer, allocatable :: c(:)
523    integer i
524
525    if(lingyun_yang) then !%%%%
526
527       if((size(ma)<ia).or.(size(mb)<ib).or.(size(mc)<ic))then
528          write(6,*) "Error caught in interface dacct 1"
529          stop
530       endif
531       if(ia/=ic)then
532          write(6,*) "Error caught in interface dacct 2"
533          stop
534       endif
535
536       !    if(ib/=nv)then
537       !     write(6,*) "Error caught in interface dacct 3"
538       !    stop
539       !    endif
540
541       allocate(c(ic))
542
543       do i=1,ic
544          call ad_alloc(c(i))
545       enddo
546
547       do i=1,ia
548          !write(*,*) ib
549          call ad_subst(ma(i), mb, ib, c(i))
550       enddo
551
552       do i=1,ic
553          call ad_copy(c(i),mc(i))
554          call ad_free(c(i))
555       enddo
556
557       deallocate(c)
558
559    else !%%%%
560       call dacct_b(ma,ia,mb,ib,mc,ic)
561    endif !%%%%
562
563
564  end subroutine dacct
565
566  subroutine mtree(mb,ib,mc,ic)
567    implicit none
568    integer,dimension(:)::mb,mc
569    integer  ib,ic
570    integer i
571    if(lingyun_yang) then !%%%%
572       if(ib/=ic) then
573          write(6,*) "In mtree not compatible to C++ TPSA "
574          stop 666
575       endif
576       do i=1,ib
577          call dacop(mb(i),mc(i))
578       enddo
579    else !%%%%
580       call  mtree_b(mb,ib,mc,ic)
581    endif !%%%%
582
583  end subroutine mtree
584
585  subroutine ppushstore(mc,nd2,coef,ml,mv)
586    implicit none
587    !
588    integer nd2
589    integer,dimension(:), intent(in)::mc
590    integer,dimension(:), intent(out)::ml,mv
591    real(dp),dimension(:),intent(out)::coef
592
593
594    if(lingyun_yang) then !%%%%
595       mv=0
596       ml=0
597       coef=0
598       write(*,*) "ppushstore should be called using the LBNL Version of Berz TPSA"
599       STOP 666
600    else !%%%%
601       call   ppushstore_b(mc,nd2,coef,ml,mv)
602    endif !%%%%
603  end subroutine ppushstore
604
605  subroutine ppushGETN(mc,ND2,ntot)
606    implicit none
607    !
608    integer ntot,ND2
609    integer,dimension(:), intent(inout)::mc
610    !
611    if(lingyun_yang) then !%%%%
612       write(*,*) "ppushGETN should be called using the LBNL Version of Berz TPSA"
613       STOP 666
614    else !%%%%
615       call   ppushGETN_b(mc,ND2,ntot)
616    endif !%%%%
617  end subroutine ppushGETN
618
619  subroutine ppush1(mc,xi,xf)
620    implicit none
621    integer mc
622    real(dp) xf
623    real(dp),dimension(:)::xi
624
625    if(lingyun_yang) then !%%%%
626       call trx_cpp(mc,xi,xf)
627    else !%%%%
628       call   ppush1_b(mc,xi,xf)
629    endif !%%%%
630  end subroutine ppush1
631
632  !  used for slow pushing
633  subroutine trx_cpp(h,xi,xf)  !,rh,y)
634    implicit none
635    real(dp) xf
636    real(dp),dimension(:)::xi
637    integer,dimension(1)::rh,hs
638    integer i,h,rh1
639    integer,dimension(lnv)::y
640    integer, allocatable :: jv(:)
641
642    if(.not.c_%stable_da) return
643    hs(1)=h
644    allocate(jv(c_%nv))
645    jv=0
646    call ad_alloc(rh1)
647    rh(1)=rh1
648
649
650    do i=1,c_%nv
651       call ad_alloc(y(i))
652       call ad_const(y(i), xi(i))
653    enddo
654
655    call dacct(hs,1,y,c_%nv,rh,1)
656    rh1=rh(1)
657
658    call ad_pek(rh1, jv, c_%nv, xf)
659
660
661    do i=1,c_%nv
662       call ad_free(y(i))
663    enddo
664    call ad_free(rh1)
665    deallocate(jv)
666
667    return
668  end subroutine trx_cpp
669
670  subroutine dainv(ma,ia,mc,ib)
671    implicit none
672    integer,dimension(:)::ma,mc
673    integer,allocatable::jj(:),ms(:),ml(:),mb(:)
674    real(dp),allocatable ::aa(:,:), ai(:,:)
675    real(dp),allocatable::ac(:)
676    integer i,j,k,ia,ib,ier
677    real(dp) amsjj
678
679    if(lingyun_yang) then !%%%%
680       allocate(jj(ib),ms(ib),ml(ib))
681       allocate(aa(ib,ib), ai(ib,ib))
682       allocate(ac(ib))
683       allocate(mb(ib))
684       jj=0
685       do i=1,ib
686          call ad_alloc(mb(i))
687          call ad_alloc(ms(i))
688          call ad_alloc(ml(i))
689       enddo
690
691       do i=1,ia
692          call dapek(ma(i), jj, ac(i))
693          call dapok(ma(i), jj, zero)
694       enddo
695
696       do i=1,ia
697          call dacon(mb(i),zero)
698          do j=1,ib
699             jj(j)=1
700             call dapek(ma(i), jj, aa(i,j))
701             call dapok(ma(i), jj, zero)
702             jj(j)=0
703          enddo
704          call dacmu(ma(i), -one, ma(i))
705       enddo
706
707
708       call matinv(aa,ai,ia,ia,ier)
709       if (ier.eq.132) then
710          write(*,*) "Can not inverse matrix:"
711          write(*,*) aa
712          STOP
713       endif
714       do i=1,ib
715          do j=1,ib
716             do k=1,ib
717                jj(k)=0
718             enddo
719             jj(j)=1
720             call dapok(mb(i),jj,ai(i,j))
721             call dapok(ml(i),jj,ai(i,j))
722          enddo
723       enddo
724       !
725       !
726       do i=2,c_%no
727          call dacct(ma,ia,mb,ib,ms,ia)
728          do j=1,ib
729             do k=1,ib
730                jj(k)=0
731             enddo
732             jj(j)=1
733             call dapek(ms(j),jj,amsjj)
734             call dapok(ms(j),jj,amsjj+one)
735          enddo
736          call dacct(ml,ia,ms,ia,mb,ib)
737       enddo
738       do i=1,ia
739          call dacmu(ma(i),-one,ma(i))
740          do j=1,ia
741             do k=1,ib
742                jj(k)=0
743             enddo
744             call dapok(ma(i),jj,ac(i))
745             jj(j)=1
746             call dapok(ma(i),jj,aa(i,j))
747          enddo
748       enddo
749
750       do i=1,ib
751          call ad_copy(mb(i), mc(i))
752          call ad_free(mb(i))
753          call ad_free(ml(i))
754          call ad_free(ms(i))
755       enddo
756
757       deallocate(mb)
758       deallocate(jj,ms,ml)
759       deallocate(aa, ai)
760       deallocate(ac)
761    else !%%%%
762       call   dainv_b(ma,ia,mc,ib)
763    endif !%%%%
764
765  end subroutine dainv
766
767  subroutine dapin(ma,ia,mb,ib,jx)
768    implicit none
769    integer ia,ib,i,nv
770    integer,dimension(:)::ma,mb,jx
771    !
772    integer,allocatable :: me(:), mn(:),mi(:),jj(:)
773
774    if(lingyun_yang) then !%%%%
775       allocate(me(ib),mn(ib),mi(ib),jj(ib))
776
777       ! me = const = zero, Identity matrix, linear part
778       do i=1,ib
779          jj(i) = 0
780       enddo
781
782       do i=1,ib
783          call ad_alloc(me(i))
784          call ad_alloc(mn(i))
785          call dapok(me(i), jj, zero)
786          jj(i)=1
787          call dapok(me(i),jj, one)
788          jj(i) = 0
789       enddo
790       ! if jx(i) == 0) mn(i)=me(i)
791       ! if(jx(i)/=0) mn(i)=ma(i)
792       if (jx(i).eq.0) then
793          call ad_copy(me(i), mn(i))
794       else
795          call ad_copy(ma(i), mn(i))
796       endif
797       ! dainv(mn,nv,mi,nv)
798       call dainv(mn,nv,mi,nv)
799       ! if jx(i) == 0) me(i)=ma(i)
800       !  dacct(me,nv,mi,nv,mb)
801       do i=1,ib
802          if (jx(i).eq.0) then
803             call ad_copy(ma(i),me(i))
804          endif
805       enddo
806       call dacct(me,ib,mi,ib,mb,ib)
807
808       deallocate(me, mn, mi,jj)
809    else !%%%%
810       call   dapin_b(ma,ia,mb,ib,jx)
811    endif !%%%%
812
813  end subroutine dapin
814
815  subroutine dader(idif,ina,inc)
816    implicit none
817    integer idif,ina,inc,itmp
818    if(lingyun_yang) then !%%%%
819       call ad_alloc(itmp)
820       call ad_derivative(ina,idif-1,itmp)
821       call ad_copy(itmp, inc)
822       call ad_free(itmp)
823    else !%%%%
824       call dader_b(idif,ina,inc)
825    endif !%%%%
826  end subroutine dader
827
828  subroutine dacfuR(ina,fun,inc)
829    implicit none
830    integer ina,inc
831    integer, allocatable :: j1d(:)
832    real(dp), allocatable :: v(:)
833    integer ns,nvar,i
834
835    interface
836       function fun(abc)
837         use precision_constants
838         implicit none
839         complex(dp) fun
840         integer,dimension(:)::abc
841       end function fun
842    end interface
843
844    if(lingyun_yang) then !%%%%
845       call ad_length(ina,ns)
846
847       call ad_nvar(nvar)
848       allocate(v(ns),j1d(ns*nvar))
849
850       call ad_read_block(ina, v, j1d, ns)
851
852       do i=1,ns
853          if(v(i)/=zero) then
854             v(i)=real(fun(j1d((i-1)*nvar+1: i*nvar)),kind=dp)*v(i)
855          endif
856       enddo
857
858       call ad_save_block(inc, v, j1d, ns)
859       deallocate(v,j1d)
860    else !%%%%
861       call   dacfuR_b(ina,fun,inc)
862    endif !%%%%
863
864
865
866  end subroutine dacfuR
867
868  subroutine dacfuI(ina,fun,inc)
869    implicit none
870    integer ina,inc
871    integer, allocatable :: j1d(:)
872    real(dp), allocatable :: v(:)
873    integer ns,nvar,i
874
875    interface
876       function fun(abc)
877         use precision_constants
878         implicit none
879         complex(dp) fun
880         integer,dimension(:)::abc
881       end function fun
882    end interface
883
884    if(lingyun_yang) then !%%%%
885       call ad_length(ina,ns)
886
887       call ad_nvar(nvar)
888       allocate(v(ns),j1d(ns*nvar))
889
890       call ad_read_block(ina, v, j1d, ns)
891
892       do i=1,ns
893          if(v(i)/=zero) then
894             v(i)=aimag(fun(j1d((i-1)*nvar+1: i*nvar)))*v(i)
895          endif
896       enddo
897
898       call ad_save_block(inc, v, j1d, ns)
899       deallocate(v,j1d)
900    else !%%%%
901       call dacfuI_b(ina,fun,inc)
902    endif !%%%%
903
904  end subroutine dacfuI
905
906  ! done
907  subroutine dacfu(ina,fun,inc)
908    implicit none
909    integer ina,inc
910    integer, allocatable :: j1d(:)
911    real(dp), allocatable :: v(:)
912    integer ns,nvar,i
913    !real(dp),external::fun
914    interface
915       function fun(abc)
916         use precision_constants
917         implicit none
918         real(dp) fun
919         integer,dimension(:)::abc
920       end function fun
921    end interface
922    if(lingyun_yang) then !%%%%
923       !pause 1
924       call ad_length(ina,ns)
925       ! write(6,*) ns
926
927       !pause 2
928       call ad_nvar(nvar)
929       ! write(6,*) nvar
930       !pause 3
931       allocate(v(ns),j1d(ns*nvar))
932       !call ad_print(ina)
933
934       call ad_read_block(ina, v, j1d, ns)
935
936       do i=1,ns
937          if(v(i)/=zero) then
938             v(i)=fun(j1d((i-1)*nvar+1: i*nvar))*v(i)
939          endif
940       enddo
941
942       call ad_save_block(inc, v, j1d, ns)
943       deallocate(v,j1d)
944    else !%%%%
945       call dacfu_b(ina,fun,inc)
946    endif !%%%%
947
948  end subroutine dacfu
949
950  ! done
951  !  subroutine GET_C_J(ina,I,C,J)
952  !    implicit none
953  !    INTEGER I,nv,ns,k,ina
954  !    integer, dimension(lnv)::j
955  !    integer, allocatable :: j1d(:)
956  !    real(dp), allocatable :: v(:)
957  !    real(dp) C
958  !    if(lingyun_yang) then !%%%%
959  !       call ad_length(ina, ns)
960  !       call ad_nvar(nv)
961  !       allocate(j1d(nv*ns))
962  !       call ad_read_block(ina, v, j1d, ns)
963  !       C = v(i)
964  !       do k=1,nv
965  !          j(k) = j1d((i-1)*nv+k)
966  !       enddo
967  !       deallocate(j1d)
968  !    else !%%%%
969  !       call GET_C_J_b(ina,I,C,J)
970  !    endif !%%%%
971  !  end  subroutine GET_C_J
972
973  subroutine dapri(ina,iunit)
974    implicit none
975    !    INTEGER ina,iunit
976    integer i,ina,ipresent,illa,iunit,iii,ioa,cn
977    real(dp) value
978    integer, allocatable :: j(:)
979    character(10) c10,k10
980
981    if(lingyun_yang) then !%%%%
982       ipresent=1
983       call dacycle(ina,ipresent,value,illa)
984
985       allocate(j(c_%nv))
986       cn=0
987       do i=1,illa
988          call dacycle(ina,i,value,illa,j)
989
990          if(abs(value)<=eps) then
991             cn=cn+1
992          endif
993       enddo
994
995       if(cn==illa) then
996          illa=0
997       endif
998
999
1000       j=0
1001       write(iunit,'(/1X,A,A,I5,A,I5,A,I5/1X,A/)') " Lingyun  ",', NO =',c_%no,', NV =',c_%nv,', INA =',ina,&
1002            '*********************************************'
1003       if(illa.ne.0) write(iunit,'(A)') '    I  COEFFICIENT          ORDER   EXPONENTS'
1004       if(illa.eq.0) write(iunit,'(A)') '   ALL COMPONENTS ZERO '
1005       c10='      NO ='
1006       k10='      NV ='
1007       i=0
1008
1009
1010       !    write(iunit,'(A10,I6,A10,I6)') c10,c_%no,k10,c_%nv
1011       cn=0
1012       do i=1,illa
1013          call dacycle(ina,i,value,illa,j)
1014
1015          ioa=0
1016          do iii=1,c_%nv
1017             ioa= j(iii)+ioa
1018          enddo
1019          if(abs(value)>eps) then
1020             cn=cn+1
1021             write(iunit,'(I6,2X,G21.14,I5,4X,18(2i2,1X))') cn,value,ioa,(j(iii),iii=1,c_%nv)
1022             !ETIENNE
1023             write(iunit,*) value
1024          endif
1025       enddo
1026       j=0
1027       write(iunit,'(A)') '                                      '
1028       deallocate(j)
1029    else !%%%%
1030       call dapri_b(ina,iunit)
1031    endif !%%%%
1032
1033
1034  end subroutine dapri
1035
1036  subroutine dapri77(ina,iunit)
1037    implicit none
1038    !    INTEGER ina,iunit
1039    integer i,ina,ipresent,illa,iunit,iii,ioa,cn
1040    real(dp) value
1041    integer, allocatable :: j(:)
1042    character(10) c10,k10
1043
1044    if(lingyun_yang) then !%%%%
1045       ipresent=1
1046       call dacycle(ina,ipresent,value,illa)
1047
1048       allocate(j(c_%nv))
1049       cn=0
1050       do i=1,illa
1051          call dacycle(ina,i,value,illa,j)
1052
1053          if(abs(value)<=eps) then
1054             cn=cn+1
1055          endif
1056       enddo
1057
1058       if(cn==illa) then
1059          illa=0
1060       endif
1061
1062
1063       j=0
1064       write(iunit,'(/1X,A,A,I5,A,I5,A,I5/1X,A/)') " Lingyun  ",', NO =',c_%no,', NV =',c_%nv,', INA =',ina,&
1065            '*********************************************'
1066       if(illa.ne.0) write(iunit,'(A)') '    I  COEFFICIENT          ORDER   EXPONENTS'
1067       if(illa.eq.0) write(iunit,'(A)') '   ALL COMPONENTS ZERO '
1068       c10='      NO ='
1069       k10='      NV ='
1070       i=0
1071
1072
1073       write(iunit,'(A10,I6,A10,I6)') c10,c_%no,k10,c_%nv
1074       cn=0
1075       do i=1,illa
1076          call dacycle(ina,i,value,illa,j)
1077
1078          ioa=0
1079          do iii=1,c_%nv
1080             ioa= j(iii)+ioa
1081          enddo
1082          if(abs(value)>eps) then
1083             cn=cn+1
1084             write(iunit,501) ioa,value,(j(iii),iii=1,c_%nv)
1085          endif
1086       enddo
1087501    format(' ', i3,1x,g23.16,1x,100(1x,i2))
1088502    format(' ', i5,1x,g23.16,1x,100(1x,i2))
1089       j=0
1090       write(iunit,502) -cn,zero,(j(iii),iii=1,c_%nv)
1091       deallocate(j)
1092    else !%%%%
1093       call dapri77_b(ina,iunit)
1094    endif !%%%%
1095
1096  end subroutine dapri77
1097
1098  subroutine dashift(ina,inc,ishift)
1099    implicit none
1100    integer ina,inc,ishift,itmp
1101    real(dp) eps
1102    if(lingyun_yang) then !%%%%
1103       call ad_alloc(itmp)
1104       eps = 1e-20
1105       call ad_shift(ina, ishift, itmp, eps)
1106       call ad_copy(itmp, inc)
1107       call ad_free(itmp)
1108    else !%%%%
1109       call dashift_b(ina,inc,ishift)
1110    endif !%%%%
1111  end  subroutine dashift
1112
1113  subroutine darea(ina,iunit)
1114    implicit none
1115    integer ina,iunit
1116    integer ii,iin,i,nno,io,iwarin,iwarnv,iwarno,io1
1117    character(10) c10
1118    integer,dimension(lnv)::j
1119    real(dp) c
1120
1121    if(lingyun_yang) then !%%%%
1122       read(iunit,'(A10)') c10
1123       read(iunit,'(18X,I4)') nno
1124       read(iunit,'(A10)') c10
1125       read(iunit,'(A10)') c10
1126       read(iunit,'(A10)') c10
1127
1128       iin = 0
1129       !
113010     continue
1131       iin = iin + 1
1132       read(iunit,'(I6,2X,G21.14,I5,4X,18(2i2,1X))') ii,c,io,(j(i),i=1,c_%nv)
1133       if(ii.eq.0) goto 20
1134       !ETIENNE
1135       read(iunit,*) c
1136       !ETIENNE
1137       if(ii.ne.iin) then
1138          iwarin = 1
1139       endif
1140       io1 = 0
1141       do i=1,c_%nv
1142          io1 = io1 + j(i)
1143       enddo
1144       !
1145       if(io1.ne.io) then
1146          iwarnv = 1
1147          goto 10
1148       endif
1149       if(io.gt.c_%no) then
1150          !        IF(IWARNO.EQ.0) PRINT*,'WARNING IN DAREA, FILE ',
1151          !    *              'CONTAINS HIGHER ORDERS THAN VECTOR '
1152          iwarno = 1
1153          goto 10
1154       endif
1155       !
1156       call ad_pok(ina, j, c_%nv, c)
1157       goto 10
1158       !
115920     continue
1160
1161       !
1162       return
1163    else !%%%%
1164       call darea_b(ina,iunit)
1165    endif !%%%%
1166
1167
1168  end subroutine darea
1169
1170
1171
1172  subroutine darea77(ina,iunit)
1173    implicit none
1174    integer ina,iunit
1175    integer nojoh,nvjoh,ii,iche,k,i
1176    character(10) c10,k10
1177    integer,dimension(lnv)::j
1178    real(dp) c
1179
1180    if(lingyun_yang) then !%%%%
1181       read(iunit,'(A10)') c10
1182       read(iunit,'(A10)') c10
1183       read(iunit,'(A10)') c10
1184       read(iunit,'(A10)') c10
1185       read(iunit,'(A10)') c10
1186       read(iunit,'(A10,I6,A10,I6)') c10,nojoh,k10,nvjoh
1187
118810     continue
1189       read(iunit,*) ii,c,(j(k),k=1,nvjoh)
1190       if(ii.lt.0) goto 20
1191
1192       do i=c_%nv+1,nvjoh
1193          if(j(i).ne.0) goto 10
1194       enddo
1195       iche=0
1196       do i=1,c_%nv
1197          iche=iche+j(i)
1198       enddo
1199       if(iche.gt.c_%no) goto 10
1200
1201       call ad_pok(ina, j, c_%nv, c)
1202
1203       goto 10
1204
120520     continue
1206
1207    else !%%%%
1208       call darea77_b(ina,iunit)
1209    endif !%%%%
1210
1211  end subroutine darea77
1212
1213  ! done
1214  !  subroutine dainf(inc,inoc,invc,ipoc,ilmc,illc)
1215  !    implicit none
1216  !    integer inc,inoc,invc,ipoc,ilmc,illc
1217  !    if(lingyun_yang) then !%%%%
1218  !      write(6,*) " FILL_UNI of Sagan only in Berz"
1219  !      STOP 666
1220  !    else !%%%%
1221  !     call dainf_b(inc,inoc,invc,ipoc,ilmc,illc)
1222  !    endif !%%%%
1223  !  end subroutine dainf
1224
1225  subroutine datra(idif,ina,inc)
1226    implicit none
1227    integer idif,ina,inc,itmp
1228    if(lingyun_yang) then !%%%%
1229       call ad_alloc(itmp)
1230       call ad_tra(ina,idif-1,itmp)
1231       call ad_copy(itmp, inc)
1232       call ad_free(itmp)
1233    else !%%%%
1234       call datra_b(idif,ina,inc)
1235    endif !%%%%
1236  end subroutine datra
1237
1238  ! done
1239  subroutine daran(ina,cm,xran)
1240    implicit none
1241    integer ina
1242    real(dp) cm,xran
1243
1244    if(lingyun_yang) then !%%%%
1245       call ad_fill_ran(ina,cm,xran)
1246    else !%%%%
1247       call daran_b(ina,cm,xran)
1248    endif !%%%%
1249
1250  end subroutine daran
1251
1252  subroutine dacycle(ina,ipresent,value,illa,j)
1253    implicit none
1254    integer illa,ina,ipresent
1255    integer,optional,dimension(:)::j
1256    real(dp) value
1257
1258    if(lingyun_yang) then !%%%%
1259       call ad_length(ina, illa)
1260       if(.not.present(j)) return
1261
1262       call ad_elem(ina, ipresent, j, value)
1263    else !%%%%
1264       call dacycle_b(ina,ipresent,value,illa,j)
1265    endif !%%%%
1266  end subroutine dacycle
1267
1268  subroutine dawritefile(ina)
1269    implicit none
1270    integer ina,n, nvar,i
1271    integer, allocatable :: j(:)
1272    real(dp),allocatable :: v(:)
1273    ! get real length of ina, without zeros in the tail
1274    call ad_length(ina, n)
1275    ! how many
1276    call ad_nvar(nvar)
1277    allocate(j(n*nvar),v(n))
1278    call ad_read_block(ina, v, j, n)
1279    do i=1,n
1280       write(6,'(1x,e15.8,5(1x,i4))') v(i),j((i-1)*nvar+1:i*nvar)
1281    enddo
1282
1283    ! for compare
1284    call ad_print(ina)
1285
1286    deallocate(j,v)
1287  end subroutine dawritefile
1288
1289  subroutine dareadfile(ina)
1290    implicit none
1291    integer ina
1292  end subroutine dareadfile
1293
1294end module dabnew !$$$$
Note: See TracBrowser for help on using the repository browser.