From 2a4aac357b8c18493ce4138b8560ee087a6c55f9 Mon Sep 17 00:00:00 2001 From: David Sagan Date: Thu, 28 Sep 2023 14:57:30 -0400 Subject: [PATCH] Update from Etienne. --- forest/code/Ci_tpsa.f90 | 201 ++- forest/code/Se_status.f90 | 16 +- forest/code/Si_def_element.f90 | 68 +- forest/code/Sr_spin.f90 | 4 + forest/code/Su_duan_zhe_map.f90 | 12 +- forest/code/a_scratch_size.f90 | 13 +- forest/code/b_da_arrays_all.f90 | 127 +- forest/code/c_dabnew.f90 | 2651 +++++++++++++++++++++++++----- forest/code/cb_da_arrays_all.f90 | 33 +- forest/code/cc_dabnew.f90 | 2492 +++++++++++++++++++++++++--- 10 files changed, 4854 insertions(+), 763 deletions(-) diff --git a/forest/code/Ci_tpsa.f90 b/forest/code/Ci_tpsa.f90 index 5caeb9158d..7d58025e91 100644 --- a/forest/code/Ci_tpsa.f90 +++ b/forest/code/Ci_tpsa.f90 @@ -148,7 +148,7 @@ MODULE c_TPSA private compute_lattice_functions_1,compute_lattice_functions_2 !private c_clean_vector,c_clean_matrix,c_clean_vector_complex,c_clean_matrix_complex private A_opt_c_vector,K_OPT_c_vector,r_field_for_demin,clean_c_universal_taylor -logical :: old_phase_calculation=.false.,secondorder=.true. +logical :: old_phase_calculation=.false. ! These routines computes lattice functions a la Ripken-Forest-Wolski @@ -3463,7 +3463,7 @@ FUNCTION from_phasor(k) enddo endif - + from_phasor=from_phasor**(k1) call kill(from_phasori) @@ -5831,13 +5831,17 @@ SUBROUTINE print_ql(S2,imaginary,quaternion_only,mf) if(.not.q_only) then write(mff,*) " Orbital Matrix " -if(present(imaginary) )write(mff,*) "Real part " +if(present(imaginary) ) then + + if(imaginary) then + write(mff,*) "Real part " do i=1,6 write(mff,'(6(1x,G21.14))') real(s2%mat(i,1:min(6,nd2))) enddo -if(present(imaginary) ) then +endif + if(imaginary) then write(mff,*) "Imaginary part " do i=1,6 @@ -5846,8 +5850,14 @@ SUBROUTINE print_ql(S2,imaginary,quaternion_only,mf) enddo endif + +else + do i=1,6 + + write(mff,'(12(1x,G21.14))') (s2%mat(i,1:min(6,nd2))) + + enddo endif - endif write(mff,*) " Quaternion Matrix " @@ -8741,6 +8751,8 @@ subroutine c_init(NO1,NV1,np1,ndpt1,AC_rf,ptc) !,spin integer, optional :: np1,ndpt1,AC_RF logical(lp), optional :: ptc !spin, integer ndpt_ptc,i,i1,i2,i3,i4,i5,i6,noo,j + real tim(0:10) + logical present_newtpsa ! if(use_quaternion) spin_def_tune=1 !NP=np1 !NO=no1 @@ -8748,6 +8760,10 @@ subroutine c_init(NO1,NV1,np1,ndpt1,AC_rf,ptc) !,spin !ND2=2*nd !!!! total dimension of phase space !nv=nd2+np !!!! total number of Taylor variables + + nphere=0 +if(use_np) nphere=np1 + ip_mat=0; jp_mat=0; jt_mat=0; do i=1,3 @@ -8767,10 +8783,14 @@ subroutine c_init(NO1,NV1,np1,ndpt1,AC_rf,ptc) !,spin formatlf(6)="(5(25x),1(1x,g23.16,1x))" ! order_gofix=no1 + if(associated(dz_c)) then + present_newtpsa=newtpsa + newtpsa=previous_newtpsa call kill(dz_c) deallocate(dz_c) nullify(dz_c) + newtpsa=present_newtpsa endif call set_da_pointers() @@ -8779,7 +8799,10 @@ subroutine c_init(NO1,NV1,np1,ndpt1,AC_rf,ptc) !,spin read77=.true. print77=.true. if(c_last_tpsa/=0) then + present_newtpsa=newtpsa + newtpsa=previous_newtpsa call c_DEASSIGN + newtpsa=present_newtpsa c_last_tpsa=1 endif @@ -8852,13 +8875,15 @@ subroutine c_init(NO1,NV1,np1,ndpt1,AC_rf,ptc) !,spin !write(6,*) "ndc2t,nd2t,nd2harm,nd2" !write(6,*) ndc2t,nd2t,nd2harm,nd2 - +!if(no==2.and.newtpsa) call dacctt2datest call c_daini(no,nv,0) +if(no==2.and.newtpsa) call c_dacctt2datest(ND2) c_master=0 ! master=1 2002.12.2 - + CALL c_ASSIGN + allocate(dz_c(nv)) call alloc(dz_c) @@ -8873,9 +8898,13 @@ subroutine c_init(NO1,NV1,np1,ndpt1,AC_rf,ptc) !,spin sj(2*i,2*i-1)=-1 enddo if(present(np1)) then - q_phasor=c_phasor() + q_phasor=c_phasor() + qi_phasor=ci_phasor() -endif + + + + endif c_%rf=>rf c_%nd2t=>nd2t c_%nd2harm=>nd2harm @@ -8897,6 +8926,7 @@ subroutine c_init(NO1,NV1,np1,ndpt1,AC_rf,ptc) !,spin if(i/=0) c_%pos_of_delta=nd2harm+1 endif + @@ -8932,7 +8962,7 @@ subroutine c_init(NO1,NV1,np1,ndpt1,AC_rf,ptc) !,spin jp_mat(3,5,6)=0 endif - + end subroutine c_init @@ -9022,7 +9052,12 @@ subroutine c_etcct(x,n1,y,n2,z) do i=1,n2 iv(i)=y(i) enddo + call c_dacct(x,n1,iv,nv,z,n1) + + + + if(nt.gt.0) then do k=1,nt call c_DADAL1(ie(k)) @@ -9038,8 +9073,9 @@ subroutine c_etinv(x,y) integer i type(c_damap) ie1,ie2 type(c_damap), intent(inout):: x,y + logical rad1 + if(.not.c_stable_da) return - ie1%n=nv;ie2%n=nv; call alloc(ie1) call alloc(ie2) @@ -9049,15 +9085,22 @@ subroutine c_etinv(x,y) do i=x%n+1,nv ie1%v(i)=1.0_dp.cmono.i enddo - + + call c_dainv(ie1%v(1:nv)%i,nv,ie2%v(1:nv)%i,nv) do i=1,x%n y%v(i)=ie2%v(i) enddo - + if(use_quaternion) THEN + rad1=assume_da_map + ! assume_da_map=.false. + y%q=x%q*y + y%q=y%q**(-1) + assume_da_map=rad1 + else y%s=x%s*y call c_inv_as(y%s,y%s) @@ -9179,23 +9222,37 @@ FUNCTION c_concat(S1,S2) enddo ! v1=s1 ! change oct 2004.10 - + ! if(old) then - + + call c_etcct(t1%v%i,t1%n,t2%v%i,t2%n,tempnew%v%i) - + + ! call c_dacct(t1%v%i,t1%n,t2%v%i,t2%n,tempnew%v%i,tempnew%n) + if(add_constant_part_concat) then do i=1,t1%n tempnew%v(i)=tempnew%v(i)+(s1%v(i).sub.'0') enddo endif - - t1%s = t1%s*t2 + + if(newspin) then + rad1=assume_da_map + t1%q = t1%q*t2 - - tempnew%s=t1%s*t2%s + tempnew%q=t1%q*t2%q + + assume_da_map=rad1 + else + t1%s = t1%s*t2 + tempnew%s=t1%s*t2%s + endif + ! t1%q = t1%q*t2 + ! t1%s = t1%s*t2 + ! tempnew%s=t1%s*t2%s + ! tempnew%q=t1%q*t2%q if(.not.c_similarity) then call c_check_rad(t1%e_ij,rad1) call c_check_rad(t2%e_ij,rad2) @@ -10472,13 +10529,17 @@ FUNCTION POWMAP( S1, R2 ) s11%n=s1%n call alloc(s11) - s11=1 + + s11=1 + R22=IABS(R2) DO I=1,R22 s11=s1*s11 ENDDO - + + + IF(R2.LT.0) THEN ! CALL c_etinv1(S11%v%i,S11%v%i,s11%n) @@ -10486,7 +10547,7 @@ FUNCTION POWMAP( S1, R2 ) ENDIF - + powmap=s11 if(complex_extra_order==1.and.special_extra_order_1) powmap=powmap.cut.no @@ -10509,6 +10570,8 @@ FUNCTION POWMAP_INV( S1, R2 ) localmaster=c_master +stop 2333 + do i=1,lnv jn(i)=0 enddo @@ -11655,11 +11718,15 @@ subroutine c_gofix(xy,a1) if(mod(i,2)==0) ndloc=i/2 ! degree of freedom of coasting plane endif enddo - + v=v.cut.(order_gofix+1) + + w=v**(-1) ! W= (Map-1)**-1 + + x=0 x%s=1 ! spin part is identity x%q=1.0_dp @@ -11667,13 +11734,16 @@ subroutine c_gofix(xy,a1) x%v(i)=1.0e0_dp.cmono.i ! Identity in all the parameter planes enddo if(ndpt/=0) x%v(ndpt)=1.0e0_dp.cmono.ndpt ! If coasting, then energy is a parameter + + a1=v + v=w*x ! v contains the fixed point, for example v(1)= eta_x * x_pt + ... - + a1=v - - ! however a1 is not a transformation, we must add the identity (done at the end) + + ! however a1 is not a transformation, we must add the identity (done at the end) ! also we must add some stuff to time to make it symplectic if coasting ! because the Lie operator which produces v(1)= eta_x * x_pt + ... ! is, within a sign, eta_x * x_pt * p_x-eta_x_prime * x_pt * x-... which affects time @@ -11767,13 +11837,12 @@ subroutine c_gofix(xy,a1) endif else ! npdt=0 - + do i=1,nd2 a1%v(i)=(1e0_dp.cmono.i)+a1%v(i) ! ndpt is already identity enddo - + endif - call kill(t1); @@ -11876,6 +11945,7 @@ subroutine c_full_canonise(at,a_cs,as,a0,a1,a2,rotation,phase,nu_spin,irot) call alloc(ri) ! at= (a,s) = (a,I) o (I,s) ! call c_full_norm_spin(at%s,kspin,norm) + call c_full_norm_spin_map(at,kspin,norm) ! storing the spin because the code is careless (sets it to zero) @@ -11902,13 +11972,13 @@ subroutine c_full_canonise(at,a_cs,as,a0,a1,a2,rotation,phase,nu_spin,irot) phi=1 -!call print(phi%v(6),6) -!pause 1 + if(ir/=1) then call extract_a1(att,a1t) else call extract_a1(att,a1t,phi) endif + ! if(present(phase)) ar=ar*phi !call print(phi%v(6),6) !pause 2 @@ -12070,18 +12140,10 @@ subroutine c_full_factorise(at,as,a0,a1,a2,dir) call extract_only_a0(att,a0t) -!call print(phi%v(6),6) -!pause 1 + call extract_only_a1(att,a1t) -! if(present(phase)) ar=ar*phi -!call print(phi%v(6),6) -!pause 2 -! if(no>1) - ! call extract_a2(att,ar) - -!call print(phi%v(6),6) -!pause 3 + a2t=att @@ -12132,17 +12194,18 @@ subroutine c_normal_spin_linear_quaternion(m_in,m_out,as,alpha) type(quaternion) q0,q1,e_y,q3,qs real(dp) alpha,cosalpha,sinalpha,tone + q0=m_in%q.sub.0 -!call print(q0) -!pause 3 + as=1 q1=q0 q1%x(0)=0.0_dp qs=1.0_dp/sqrt(q1%x(1)**2+q1%x(2)**2+q1%x(3)**2) + + q1=q1*qs ! q1=n -!call print(q1) -!pause 4 + e_y=0.0_dp e_y%x(2)=1.0_dp !call print(e_y) @@ -16800,20 +16863,24 @@ subroutine extract_a1(a,a1,phi1) endif - + a=b1**(-1)*a a1=b1 ! imposed Teng-Edward A_12=0 or, for fun, Anti-Teng-Edwards A_21=0 + if(present(phi1)) then phi1=1 + do i=1,nd2/2 if((i<=ndt).or.(i>nd-rf)) then if(courant_snyder_teng_edwards) then + t=sqrt((b1%v(2*i-1).d.(2*i))**2+(b1%v(2*i-1).d.(2*i-1))**2) + cphi=(b1%v(2*i-1).d.(2*i-1))/t sphi=(b1%v(2*i-1).d.(2*i))/t else @@ -16835,7 +16902,7 @@ subroutine extract_a1(a,a1,phi1) endif enddo - + a1=a1*phi1**(-1) a=phi1*a*phi1**(-1) @@ -18286,6 +18353,7 @@ subroutine get_c_yu_w(n,yu,a0,a1,a2,ugiven) u=ugiven else call c_canonise(n%a_t,u,a0=b0,a1=b1,a2=bn) +n%a_t=u ui=ci_phasor()*bn*c_phasor() u=ui**(-1) f=n%ker @@ -18804,34 +18872,40 @@ subroutine c_normal(xyso3,n,dospin,no_used,rot,phase,nu_spin,canonize) ! but energy is constant. (Momentum compaction, phase slip etc.. falls from there) ! etienne - if(c_skip_gofix) then + if(c_skip_gofix) then a1=1 else call c_gofix(m1,a1) endif - m1=c_simil(a1,m1,-1) + + m1=c_simil(a1,m1,-1) + + + ! Does the the diagonalisation into a rotation call c_linear_a(m1,a2) + + !!! Now the linear map is normalised m1=c_simil(a2,m1,-1) + - !!! We go into the phasors' basis ri=from_phasor(-1) - m1=c_simil(ri,m1,1) - +!stop 999 ri=(m1.sub.-1)**(-1) + ri%s=1 ! make spin identity ri%q=1.0_dp ! make spin identity - + !!! The tunes are stored for the nonlinear normal form recursive algorithm do k=1,xy%n @@ -18858,6 +18932,7 @@ subroutine c_normal(xyso3,n,dospin,no_used,rot,phase,nu_spin,canonize) nonl= exp_inv(n%ker,nonl) nonl=nonl.sub.i + do k=1,xy%n if(lielib_print(13)/=0) then @@ -18925,14 +19000,18 @@ subroutine c_normal(xyso3,n,dospin,no_used,rot,phase,nu_spin,canonize) enddo ! over vector index m1=c_simil(n%g%f(i),m1,-1) - + enddo + + + n%a_t=a1*a2*from_phasor()*texp(n%g)*from_phasor(-1) - + ! n%a1=a1 n%a2=a2 + !!!!! here we put the normalised linear part into the factored vector field !!!!! not necessary but useful do k=1,xy%n @@ -19001,7 +19080,9 @@ subroutine c_normal(xyso3,n,dospin,no_used,rot,phase,nu_spin,canonize) if(use_quaternion)then + call c_normal_spin_linear_quaternion(m1,m1,n%AS,alpha) + n%quaternion_angle=alpha/2.0_dp ri=1 ; ri%q=m1%q.sub.0 ; ! exp(theta_0 L_y) (2) ! sx=sqrt(ri%q%x(1)**2+ri%q%x(2)**2+ri%q%x(3)**2) @@ -19052,6 +19133,8 @@ subroutine c_normal(xyso3,n,dospin,no_used,rot,phase,nu_spin,canonize) !!! tune is taken from egspin(1) or egspin(3) spin_def_tune= +/- 1 n%spin_tune=aimag(log(egspin(2+spin_def_tune))/twopi) + + ! because exp(a L_y) x = x- a z + O(a**2) ri=ri**(-1) ! exp(-alpha_0 L_y) (3) @@ -19072,7 +19155,8 @@ subroutine c_normal(xyso3,n,dospin,no_used,rot,phase,nu_spin,canonize) write(mkers,*) "Order ",i endif - + + mt=m1*ri ! S*exp(-theta_0 L_y) (5) @@ -19147,6 +19231,9 @@ subroutine c_normal(xyso3,n,dospin,no_used,rot,phase,nu_spin,canonize) call c_nr_to_n0(nr,nr) ! (10) +!write(6,*) " i spin ",i +!call print(nr,6) + if(use_quaternion)then qnr=nr diff --git a/forest/code/Se_status.f90 b/forest/code/Se_status.f90 index 085e67c48c..05f4d82ea6 100644 --- a/forest/code/Se_status.f90 +++ b/forest/code/Se_status.f90 @@ -144,7 +144,7 @@ module S_status private track_TREE_G_complexr,track_TREE_G_complexp,track_TREE_probe_complexr,track_TREE_probe_complexp_new real(dp),TARGET ::INITIAL_CHARGE=1 logical :: mcmillan=.false. - real(dp) :: radfac=1 ! to fudge radiation (lower it) + real(dp) :: radfac=1,flucfac=1 ! to fudge radiation (lower it) TYPE B_CYL integer firsttime integer, POINTER :: nmul,n_mono !,nmul_e,n_mono_e @@ -264,7 +264,7 @@ real(dp) function cflucf(p) ! if(junk_e) then ! cflucf=cfluc*p%p0c**5 ! else - cflucf=cfluc0*twopii/p%GAMMA0I**5/p%MASS**2 + cflucf=flucfac*cfluc0*twopii/p%GAMMA0I**5/p%MASS**2 ! endif end function cflucf @@ -1529,6 +1529,8 @@ subroutine S_init(STATE,NO1,NP1,pack,ND2,NPARA,number_of_clocks) INTEGER,optional :: ND2,NPARA,number_of_clocks INTEGER ND2l,NPARAl,n_acc,no1c,nv,i LOGICAL(lp) package + + do_damping=.false. do_spin=.false. if(state%radiation) do_damping=.true. @@ -1589,7 +1591,7 @@ subroutine S_init(STATE,NO1,NP1,pack,ND2,NPARA,number_of_clocks) ! endif ! write(6,*) NO1,ND1,NP1,NDEL,NDPT1 !pause 678 - + CALL INIT(NO1,ND1,NP1+NDEL+2*n_acc,NDPT1,PACKAGE) nv=2*nd1+NP1+NDEL+2*n_acc @@ -1611,12 +1613,14 @@ subroutine S_init(STATE,NO1,NP1,pack,ND2,NPARA,number_of_clocks) ND1=ND1+n_acc if(use_complex_in_ptc) call c_init(NO1c,nd1,np1+ndel,ndpt1,n_acc,ptc=my_false) ! PTC false because we will not use the real FPP for acc modulation n_rf=n_acc - + previous_newtpsa=newtpsa END subroutine S_init subroutine kill_map_cp() implicit none - + logical present_newtpsa + present_newtpsa=newtpsa + newtpsa=previous_newtpsa if(associated(dz_8)) then call kill(dz_8) deallocate(dz_8) @@ -1635,7 +1639,7 @@ subroutine kill_map_cp() deallocate(dz_c) nullify(dz_c) endif - + newtpsa=present_newtpsa end subroutine kill_map_cp diff --git a/forest/code/Si_def_element.f90 b/forest/code/Si_def_element.f90 index cc6b1a0dfa..585482c887 100644 --- a/forest/code/Si_def_element.f90 +++ b/forest/code/Si_def_element.f90 @@ -2031,68 +2031,68 @@ SUBROUTINE transfer_ANBN(EL,ELP,VR,DVR,VP,DVP,T) if(EL%KIND==kind1) return - if(associated(EL%ramp)) then - - if(EL%KIND/=kind15) then - do n=1,EL%P%NMUL - EL%BN(N)= EL%ramp%table(0)%bn(n) - EL%AN(N)= EL%ramp%table(0)%an(n) - ELP%BN(N)= ELP%ramp%table(0)%bn(n) - ELP%AN(N)= ELP%ramp%table(0)%an(n) - enddo - else - EL%VOLT=EL%ramp%table(0)%bn(1)*COS(twopi*EL%ramp%table(0)%an(1)*T/clight+EL%ramp%table(0)%bn(2))+EL%ramp%table(0)%an(2) - ELP%VOLT=EL%ramp%table(0)%bn(1)*COS(twopi*EL%ramp%table(0)%an(1)*T/clight+EL%ramp%table(0)%bn(2))+EL%ramp%table(0)%an(2) - write(6,*) " volt ",el%volt,EL%ramp%table(0)%bn(1) - endif - - if(EL%ramp%table(0)%b_t/=0.0_dp) then - if(EL%parent_fibre%PATCH%TIME==0) EL%parent_fibre%PATCH%TIME=2 - if(EL%parent_fibre%PATCH%TIME==1) EL%parent_fibre%PATCH%TIME=3 - EL%parent_fibre%PATCH%b_T=EL%ramp%table(0)%b_t - else - if(EL%parent_fibre%PATCH%TIME==2) EL%parent_fibre%PATCH%TIME=0 - if(EL%parent_fibre%PATCH%TIME==3) EL%parent_fibre%PATCH%TIME=1 - EL%parent_fibre%PATCH%b_T=0.0_dp - endif +! if(associated(EL%ramp)) then +! +! if(EL%KIND/=kind15) then +! do n=1,EL%P%NMUL +! EL%BN(N)= EL%ramp%table(0)%bn(n) +! EL%AN(N)= EL%ramp%table(0)%an(n) +! ELP%BN(N)= ELP%ramp%table(0)%bn(n) +! ELP%AN(N)= ELP%ramp%table(0)%an(n) +! enddo +! else +! EL%VOLT=EL%ramp%table(0)%bn(1)*COS(twopi*EL%ramp%table(0)%an(1)*T/clight+EL%ramp%table(0)%bn(2))+EL%ramp%table(0)%an(2) +! ELP%VOLT=EL%ramp%table(0)%bn(1)*COS(twopi*EL%ramp%table(0)%an(1)*T/clight+EL%ramp%table(0)%bn(2))+EL%ramp%table(0)%an(2) +! write(6,*) " volt ",el%volt,EL%ramp%table(0)%bn(1) +! endif +! +! if(EL%ramp%table(0)%b_t/=0.0_dp) then +! if(EL%parent_fibre%PATCH%TIME==0) EL%parent_fibre%PATCH%TIME=2 +! if(EL%parent_fibre%PATCH%TIME==1) EL%parent_fibre%PATCH%TIME=3 +! EL%parent_fibre%PATCH%b_T=EL%ramp%table(0)%b_t +! else +! if(EL%parent_fibre%PATCH%TIME==2) EL%parent_fibre%PATCH%TIME=0 +! if(EL%parent_fibre%PATCH%TIME==3) EL%parent_fibre%PATCH%TIME=1 +! EL%parent_fibre%PATCH%b_T=0.0_dp + ! endif - else + ! else IF(EL%P%NMUL>=1.and.associated(EL%D_BN)) THEN if(present(VR))then do n=1,EL%P%NMUL EL%BN(N)= vR*EL%D0_BN(N)+DVR*EL%D_BN(N) EL%AN(N)= vR*EL%D0_AN(N)+DVR*EL%D_AN(N) - ELP%BN(N)= vR*EL%D0_BN(N)+DVR*EL%D_BN(N) - ELP%AN(N)= vR*EL%D0_AN(N)+DVR*EL%D_AN(N) + ELP%BN(N)= vR*ELp%D0_BN(N)+DVR*ELp%D_BN(N) + ELP%AN(N)= vR*ELp%D0_AN(N)+DVR*ELp%D_AN(N) enddo else do n=1,EL%P%NMUL EL%BN(N)= vp*EL%D0_BN(N)+DVp*EL%D_BN(N) EL%AN(N)= vp*EL%D0_AN(N)+DVp*EL%D_AN(N) - ELP%BN(N)= vp*EL%D0_BN(N)+DVp*EL%D_BN(N) - ELP%AN(N)= vp*EL%D0_AN(N)+DVp*EL%D_AN(N) + ELP%BN(N)= vp*ELp%D0_BN(N)+DVp*ELp%D_BN(N) + ELP%AN(N)= vp*ELp%D0_AN(N)+DVp*ELp%D_AN(N) enddo endif - endif + ! endif if(associated(el%volt)) then if(present(VR))then EL%volt= vR*EL%D0_Volt+DVR*EL%D_Volt - ELP%volt= vR*EL%D0_Volt+DVR*EL%D_Volt + ELP%volt= vR*ELp%D0_Volt+DVR*ELp%D_Volt else EL%volt= vp*EL%D0_Volt+DVp*EL%D_Volt - ELP%volt= vp*EL%D0_Volt+DVp*EL%D_Volt + ELP%volt= vp*ELp%D0_Volt+DVp*ELp%D_Volt endif endif if(associated(el%phas)) then if(present(VR))then EL%phas= vR*EL%D0_phas+DVR*EL%D_phas - ELP%phas= vR*EL%D0_phas+DVR*EL%D_phas + ELP%phas= vR*ELp%D0_phas+DVR*ELp%D_phas else EL%phas= vp*EL%D0_phas+DVp*EL%D_phas - ELP%phas= vp*EL%D0_phas+DVp*EL%D_phas + ELP%phas= vp*ELp%D0_phas+DVp*ELp%D_phas endif endif diff --git a/forest/code/Sr_spin.f90 b/forest/code/Sr_spin.f90 index 0e237b2ada..af444b0ff9 100644 --- a/forest/code/Sr_spin.f90 +++ b/forest/code/Sr_spin.f90 @@ -1572,8 +1572,10 @@ SUBROUTINE TRACK_NODE_FLAG_probe_QUAP(C,XS,K) ! endif IF(K%MODULATION.and.xs%nac/=0) then + knob=.true. if(c%parent_fibre%mag%slow_ac/=0) CALL MODULATE(C,XS,K) !modulate CALL TRACK_MODULATION(C,XS,K) !modulate + knob=.false. ENDIF !modulate @@ -1939,8 +1941,10 @@ SUBROUTINE TRACK_NODE_FLAG_probe_P(C,XS,K) ! endif IF(K%MODULATION.and.xs%nac/=0) then + knob=.true. if(c%parent_fibre%mag%slow_ac/=0) CALL MODULATE(C,XS,K) !modulate CALL TRACK_MODULATION(C,XS,K) !modulate + knob=.false. ENDIF !modulate diff --git a/forest/code/Su_duan_zhe_map.f90 b/forest/code/Su_duan_zhe_map.f90 index 1bb3d4bcf9..b17e4a6bdf 100644 --- a/forest/code/Su_duan_zhe_map.f90 +++ b/forest/code/Su_duan_zhe_map.f90 @@ -1073,18 +1073,19 @@ end SUBROUTINE orthonormalise !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! new zhe tracking !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - SUBROUTINE track_TREE_probe_complex_zhe(T,xs,spin,rad,stoch,linear,slim) + SUBROUTINE track_TREE_probe_complex_zhe(T,xs,spin,rad,stoch,linear,slim,flucfac) ! use da_arrays IMPLICIT NONE TYPE(TREE_ELEMENT),target, INTENT(INout) :: T(3) type(probe) xs real(dp) x(size_tree),x0(size_tree),s0(3,3),r(3,3),dx6,beta,q(3),p(3),qg(3),qf(3) - real(dp) normb,norm,x0_begin(size_tree),xr(6),normbb + real(dp) normb,norm,x0_begin(size_tree),xr(6),normbb,flucfac1 integer i,j,k,ier,is logical, optional :: spin,stoch,rad,linear,slim logical spin0,stoch0,rad0,doit,as_is0,slim0 integer no1 + real(dp), optional :: flucfac type(quaternion) qu as_is0=t(1)%usenonsymp spin0=.true. @@ -1099,7 +1100,8 @@ SUBROUTINE track_TREE_probe_complex_zhe(T,xs,spin,rad,stoch,linear,slim) if(present(spin)) spin0=spin if(present(stoch)) stoch0=stoch if(present(rad)) rad0=rad - + flucfac1=1 + if(present(flucfac)) flucfac1=flucfac if(as_is0) rad0=.true. doit=rad0.or.stoch0 x=0 @@ -1117,11 +1119,11 @@ SUBROUTINE track_TREE_probe_complex_zhe(T,xs,spin,rad,stoch,linear,slim) xr=0.0_dp if(use_gaussian_zhe) then do i=1,6 - xr(i)=GRNF_zhe_gaussian()*t(2)%fix0(i) + xr(i)=flucfac1*GRNF_zhe_gaussian()*t(2)%fix0(i) enddo else do i=1,6 - xr(i)=GRNF_zhe()*t(2)%fix0(i) + xr(i)=flucfac1*GRNF_zhe()*t(2)%fix0(i) enddo endif xr =matmul(t(2)%rad,xr) diff --git a/forest/code/a_scratch_size.f90 b/forest/code/a_scratch_size.f90 index 019ec769d6..76bc250d0b 100644 --- a/forest/code/a_scratch_size.f90 +++ b/forest/code/a_scratch_size.f90 @@ -226,6 +226,8 @@ module precision_constants character(18) :: format8="(8(1x,g23.16,1x))" character(18) :: format9="(9(1x,g23.16,1x))" character(19) :: format10="(10(1x,g23.16,1x))" + character(19) :: format11="(11(1x,g23.16,1x))" + character(19) :: format12="(12(1x,g23.16,1x))" ! logical(lp) :: fixed_found ! lielib_print(1)=1 lieinit prints info ! lielib_print(2)=1 expflo warning if no convergence @@ -246,7 +248,16 @@ module precision_constants ! lielib_print(17)=1 print magnets with excessive cutting integer , target :: old_integrator =1 ! before making spin high order character*255 :: preffile="pref.txt" - + logical :: newtpsa=.false., assume_da_map = .true.,previous_newtpsa=.false.,use_np=.true.,newspin=.true. + integer :: with_para=2 + integer :: nphere=0 + integer(4), pointer :: inds(:,:),ind1(:),ind2(:) => null() + integer(4), pointer :: nind1(:),nind2(:) => null() + integer(4) combien, poscombien,ninds + real(dp), pointer:: reel(:),finds(:,:,:) => null(),finds1(:) => null() + complex(dp), pointer:: reelc(:) => null() + integer :: ipric=0,ipric2=0 + INTERFACE read MODULE PROCEDURE read_d MODULE PROCEDURE read_int diff --git a/forest/code/b_da_arrays_all.f90 b/forest/code/b_da_arrays_all.f90 index 69acd74218..382174939f 100644 --- a/forest/code/b_da_arrays_all.f90 +++ b/forest/code/b_da_arrays_all.f90 @@ -9,13 +9,14 @@ module da_arrays implicit none public integer lda,lea,lia,lst + integer, private :: nohere,nvhere integer,private,parameter::nmax=400 real(dp),private,parameter::tiny=1e-20_dp ! johan ! integer,parameter::lno=2,lnv=6,lnomax=8,lnvmax=9,lstmax=800500,ldamax=16000,leamax=5000,liamax=50000 ! - integer,parameter :: lnv=100,lno=200 + integer,parameter :: lnv=200,lno=200 integer :: lstmax=800500,ldamax=16000,leamax=5000,liamax=50000 !integer,parameter::lno=200,lnv=100,lnomax=8,lnvmax=9,lstmax=800500,ldamax=16000,leamax=5000,liamax=50000 logical(lp) :: reallocate = .true. @@ -45,7 +46,15 @@ module da_arrays subroutine alloc_all(no,nv) implicit none integer no,nv +nohere=no +nvhere=nv if(reallocate) then + if(associated(ind1))deallocate(ind1) + if(associated(ind2))deallocate(ind2) + if(associated(inds))deallocate(inds) + if(associated(reel))deallocate(reel) + if(associated(reelc))deallocate(reelc) +! if(Allocated(comp))deallocate(comp) call dealloc_all call danum0(no,nv) call alloc_ @@ -63,10 +72,18 @@ end subroutine alloc_all subroutine alloc_ implicit none +integer i,j,k allocate(cc(lst)) -allocate(i_1(lst));allocate(i_2(lst)); -allocate(ie1(lea));allocate(ie2(lea));allocate(ieo(lea)); -allocate(ia1(0:lia));allocate(ia2(0:lia)); +if(nohere>2.or.(.not.newtpsa)) then + allocate(i_1(lst));allocate(i_2(lst)); + allocate(ie1(lea));allocate(ie2(lea));allocate(ieo(lea)); + allocate(ia1(0:lia));allocate(ia2(0:lia)); +else + allocate(i_1(1));allocate(i_2(1)); + allocate(ie1(1));allocate(ie2(1));allocate(ieo(1)); + allocate(ia1(0:1));allocate(ia2(0:1)); + +endif allocate(idano(lda));allocate(idanv(lda));allocate(idapo(lda)); allocate(idalm(lda));allocate(idall(lda)); allocate(daname(lda)); @@ -85,6 +102,108 @@ subroutine alloc_ idapo=0 idalm=0 idall=0 +combien=1 + +if(nohere==1.and.newtpsa) then + if(associated(ind1)) then + DEALLOCATE(ind1,ind2,inds,reel,reelc); + endif + if(associated(nind2)) then + DEALLOCATE(nind1,nind2); + ninds=0 + endif +combien=nvhere+1 +allocate(ind1(combien),ind2(combien),inds(0:nvhere,0:nvhere)) +allocate(reel(combien),reelc(combien) ) !,comp(combien)) +reel=0.0_dp +reelc=0.0_dp + +!comp=0.0_dp +k=nvhere +inds(0,0)=1 +ind1(1)=0 +ind2(1)=0 +k=1 +do i=1,nvhere + k=k+1 + ind1(k)=i + ind2(k)=0 + inds(i,0)=k + inds(0,i)=k +enddo +poscombien=0 + + endif + + +if(nohere==2.and.newtpsa) then + if(associated(ind1)) then + DEALLOCATE(ind1,ind2,inds,reel,reelc,finds); + endif + if(associated(finds)) then + DEALLOCATE(finds); + endif + if(associated(finds1)) then + DEALLOCATE(finds1); + endif + if(associated(nind2)) then + DEALLOCATE(nind1,nind2); + ninds=0 + endif + +combien=(nvhere+2)*(nvhere+1)/2 +allocate(ind1(combien),ind2(combien),inds(0:nvhere,0:nvhere)) +allocate(reel(combien),reelc(combien) ) !,comp(combien)) +reel=0.0_dp +reelc=0.0_dp + +!comp=0.0_dp +k=nvhere +inds(0,0)=1 +ind1(1)=0 +ind2(1)=0 +k=1 +do i=1,nvhere + k=k+1 + ind1(k)=i + ind2(k)=0 + inds(i,0)=k + inds(0,i)=k +enddo +poscombien=k+1 + +if(with_para<2) then + allocate(finds(nvhere,nvhere,poscombien:combien )) + finds=1 +endif + +do i=1,nvhere +do j=i,nvhere + k=k+1 + ind1(k)=i + ind2(k)=j + inds(i,j)=k + inds(j,i)=k +enddo +enddo + +if(with_para<2) then + + do i=1,nvhere + do k=1,nvhere + do j=poscombien,combien + if(ind1(j)/=ind2(j)) finds(i,k,j)=2 + enddo + enddo + enddo + +endif + !write(6,*) k,poscombien,combien +!pause 765 +!do k=1,combien !+ finds(ind1(jk),ind2(jk),ab)* & +!write(6,*) ind1(k),ind2(k),inds(ind1(k),ind2(k)),inds(ind2(k),ind1(k)) +!enddo + endif end subroutine alloc_ diff --git a/forest/code/c_dabnew.f90 b/forest/code/c_dabnew.f90 index 8dab373738..d0bcb05847 100644 --- a/forest/code/c_dabnew.f90 +++ b/forest/code/c_dabnew.f90 @@ -9,8 +9,8 @@ module dabnew ! private public public ldamax, lstmax,leamax,liamax - private daallno1,daall,damult,dasqrt,dacmut,dacma,DALINt,dafunt,dacctt - private dainvt,dapint,dadert,dacfuRt,dacfuIt,dacfut + private daallno1,daall,damult,dasqrt,dacmut,dacma,DALINt,dafunt,dacctt,dacctt2da,dacctt2tpsa + private dainvt,dapint,dadert,dacfut,dacctt1,dainvt1,dainvt2 private dadeb,dapac,dachk,damch,dadcd,dancd,hash,dehash ! public ppushstore,ppushprint,ppushGETN ! public DEALLOC_ALL,dainf @@ -23,10 +23,15 @@ module dabnew ! public count_da ! write_ removed integer,private,parameter:: lsw=1 - integer :: lda_max_used=0 + integer :: lda_max_used=0 + ! integer,private,parameter::nmax=400,lsw=1 ! real(dp),private,parameter::tiny=c_1d_20 character(120),private :: line + + + + contains !****************************************************************************** ! * @@ -212,6 +217,7 @@ subroutine change_package(i) write(6,*) " INPUT IGNORED " endif end subroutine change_package + subroutine daini(no,nv,iunit) implicit none ! ***************************** @@ -235,6 +241,61 @@ subroutine daini(no,nv,iunit) integer,dimension(0:lnv)::k integer,dimension(lnv)::j,jj character(10) aa + + nomax = no + nvmax = nv +last_tpsa=2 + +if(newtpsa) then + if(nomax.eq.1) then + call alloc_all(no,nv) + ndamaxi=0 + do i=1, lda + allvec(i) = .false. + enddo + nhole=0 + nda_dab = 0 + nst0 = 0 + nomax = no + nvmax = nv + nocut=no + call danum(no,nv,nmmax) + + do i=1,lda + idapo(i)=nmmax*(i-1)+1 + enddo + idano=1 + idanv=nv + idalm=nmmax + idall=nmmax + return + endif + + if(nomax==2) then + call alloc_all(no,nv) + ndamaxi=0 + do i=1, lda + allvec(i) = .false. + enddo + nhole=0 + nda_dab = 0 + nst0 = 0 + nomax = no + nvmax = nv + nocut=no + call danum(no,nv,nmmax) + + do i=1,lda + idapo(i)=nmmax*(i-1)+1 + enddo + idano=2 + idanv=nv + idalm=nmmax + idall=nmmax + return + endif + +endif ! !frs if(eps.le.zero) eps=c_1d_38 ! if(EPS.le.zero) eps=c_1d_90 @@ -245,9 +306,10 @@ subroutine daini(no,nv,iunit) endif return endif - last_tpsa=2 +! last_tpsa=2 if(nv.eq.0) return call alloc_all(no,nv) + ndamaxi=0 ! do i=1, lda @@ -263,6 +325,7 @@ subroutine daini(no,nv,iunit) nomax = no nvmax = nv call danum(no,nv,nmmax) + nocut = no lfi = 0 ! @@ -569,7 +632,73 @@ subroutine daallno1(ic,ccc) ! endif ! return ! endif - ! +if(newtpsa) then + if(nomax.eq.1) then + + if(nhole.gt.0) then + ind=nda_dab +200 if (allvec(ind)) then + ind = ind - 1 + goto 200 + endif + incnda = .false. + nhole=nhole-1 + else + incnda = .true. + nda_dab = nda_dab + 1 + ind=nda_dab + if(nda_dab.gt.lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 200 + endif + endif + if(ind>lda_max_used) lda_max_used=ind + if(ind>lda) then + write(6,*) "ind>lda ",lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA + stop 201 + endif + allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + daname(ind) = ccc + ic = ind + return + endif + + if(nomax==2) then + + if(nhole.gt.0) then + ind=nda_dab +300 if (allvec(ind)) then + ind = ind - 1 + goto 300 + endif + incnda = .false. + nhole=nhole-1 + else + incnda = .true. + nda_dab = nda_dab + 1 + ind=nda_dab + if(nda_dab.gt.lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 300 + endif + endif + if(ind>lda_max_used) lda_max_used=ind + if(ind>lda) then + write(6,*) "ind>lda ",lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA + stop 301 + endif + allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + daname(ind) = ccc + ic = ind + return + endif +endif + + no=nomax nv=nvmax ind = 1 @@ -651,6 +780,7 @@ subroutine daallno1(ic,ccc) if(nda_dab.gt.ndamaxi) ndamaxi=nda_dab return end subroutine daallno1 + subroutine daall(ic,l,ccc,no,nv) implicit none ! ******************************** @@ -669,7 +799,76 @@ subroutine daall(ic,l,ccc,no,nv) ! write(6,*) "big problem in dabnew ", sqrt(crash) ! endif ! return -! endif +! +if(newtpsa) then + if(nomax.eq.1) then +do i=1,l + if(nhole.gt.0) then + ind=nda_dab +200 if (allvec(ind)) then + ind = ind - 1 + goto 200 + endif + incnda = .false. + nhole=nhole-1 + else + incnda = .true. + nda_dab = nda_dab + 1 + ind=nda_dab + if(nda_dab.gt.lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 200 + endif + endif + if(ind>lda_max_used) lda_max_used=ind + if(ind>lda) then + write(6,*) "ind>lda ",lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA + stop 201 + endif + allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + daname(ind) = ccc + ic(i) = ind + enddo + return + endif + + + if(nomax==2) then + do i=1,l + if(nhole.gt.0) then + ind=nda_dab +300 if (allvec(ind)) then + ind = ind - 1 + goto 300 + endif + incnda = .false. + nhole=nhole-1 + else + incnda = .true. + nda_dab = nda_dab + 1 + ind=nda_dab + if(nda_dab.gt.lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 300 + endif + endif + if(ind>lda_max_used) lda_max_used=ind + if(ind>lda) then + write(6,*) "ind>lda ",lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA + stop 301 + endif + allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + daname(ind) = ccc + ic(i) = ind + enddo + return + endif +endif + ! ind = 1 do i=1,l @@ -755,6 +954,7 @@ subroutine daall(ic,l,ccc,no,nv) if(nda_dab.gt.ndamaxi) ndamaxi=nda_dab return end subroutine daall + subroutine daall1(ic,ccc,no,nv) implicit none ! ******************************** @@ -767,6 +967,73 @@ subroutine daall1(ic,ccc,no,nv) logical(lp) incnda integer ic,ind,ndanum,no,nv,ipause,mypauses character(10) c,ccc +if(newtpsa) then + if(nomax.eq.1) then + + if(nhole.gt.0) then + ind=nda_dab +200 if (allvec(ind)) then + ind = ind - 1 + goto 200 + endif + incnda = .false. + nhole=nhole-1 + else + incnda = .true. + nda_dab = nda_dab + 1 + ind=nda_dab + if(nda_dab.gt.lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 200 + endif + endif + if(ind>lda_max_used) lda_max_used=ind + if(ind>lda) then + write(6,*) "ind>lda ",lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA + stop 201 + endif + allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + daname(ind) = ccc + ic = ind + return + endif + + if(nomax==2) then + + if(nhole.gt.0) then + ind=nda_dab +300 if (allvec(ind)) then + ind = ind - 1 + goto 300 + endif + incnda = .false. + nhole=nhole-1 + else + incnda = .true. + nda_dab = nda_dab + 1 + ind=nda_dab + if(nda_dab.gt.lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 300 + endif + endif + if(ind>lda_max_used) lda_max_used=ind + if(ind>lda) then + write(6,*) "ind>lda ",lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA + stop 301 + endif + allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + daname(ind) = ccc + ic = ind + return + endif +endif + + ! if((.not.C_%STABLE_DA)) then ! if(c_%watch_user) then ! write(6,*) "big problem in dabnew ", sqrt(crash) @@ -854,6 +1121,7 @@ subroutine daall1(ic,ccc,no,nv) if(nda_dab.gt.ndamaxi) ndamaxi=nda_dab return end subroutine daall1 + subroutine daall0(ic) implicit none ! ******************************** @@ -866,6 +1134,73 @@ subroutine daall0(ic) logical(lp) incnda integer ic,ind,ndanum,no,nv,ipause,mypauses character(10) c,ccc + +if(newtpsa) then + if(nomax.eq.1) then + + if(nhole.gt.0) then + ind=nda_dab +200 if (allvec(ind)) then + ind = ind - 1 + goto 200 + endif + incnda = .false. + nhole=nhole-1 + else + incnda = .true. + nda_dab = nda_dab + 1 + ind=nda_dab + if(nda_dab.gt.lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 200 + endif + endif + if(ind>lda_max_used) lda_max_used=ind + if(ind>lda) then + write(6,*) "ind>lda ",lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA + stop 201 + endif + allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + daname(ind) = ccc + ic = ind + return + endif + + if(nomax==2) then + + if(nhole.gt.0) then + ind=nda_dab +300 if (allvec(ind)) then + ind = ind - 1 + goto 300 + endif + incnda = .false. + nhole=nhole-1 + else + incnda = .true. + nda_dab = nda_dab + 1 + ind=nda_dab + if(nda_dab.gt.lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 300 + endif + endif + if(ind>lda_max_used) lda_max_used=ind + if(ind>lda) then + write(6,*) "ind>lda ",lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',LDA + stop 301 + endif + allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + daname(ind) = ccc + ic = ind + return + endif +endif + ccc=' ' no=nomax nv=nvmax @@ -968,23 +1303,55 @@ subroutine dadal(idal,l) integer i,l,ipause,mypauses integer,dimension(:)::idal ! -! if((.not.C_%STABLE_DA)) then -! if(c_%watch_user) then -! write(6,*) "big problem in dabnew ", sqrt(crash) -! endif -! return -! endif +if(newtpsa) then + if(nomax.eq.1) then do i=l,1,-1 - if(idal(i).le.nomax+2.or.idal(i).gt.nda_dab) then - write(line,'(a38,i8,1x,i8)') 'ERROR IN ROUTINE DADAL, IDAL(I),NDA = ',idal(i),nda_dab -! ipause=mypauses(13,line) - C_%STABLE_DA = .false. - l = 1 - return - call dadeb !(31,'ERR DADAL ',1) - endif - if(idal(i).eq.nda_dab) then - ! deallocate + if(idal(i).eq.nda_dab) then + ! deallocate + nda_dab = nda_dab - 1 + else + nhole=nhole+1 + endif + cc(idapo(idal(i)):idapo(idal(i))+nmmax-1)=0 + allvec(idal(i)) = .false. + idal(i) = 0 + enddo + return + endif + + if(nomax==2) then + do i=l,1,-1 + if(idal(i).eq.nda_dab) then + ! deallocate + nda_dab = nda_dab - 1 + else + nhole=nhole+1 + endif + cc(idapo(idal(i)):idapo(idal(i))+nmmax-1)=0 + allvec(idal(i)) = .false. + idal(i) = 0 + enddo + return + endif +endif + +! if((.not.C_%STABLE_DA)) then +! if(c_%watch_user) then +! write(6,*) "big problem in dabnew ", sqrt(crash) +! endif +! return +! endif + do i=l,1,-1 + if(idal(i).le.nomax+2.or.idal(i).gt.nda_dab) then + write(line,'(a38,i8,1x,i8)') 'ERROR IN ROUTINE DADAL, IDAL(I),NDA = ',idal(i),nda_dab +! ipause=mypauses(13,line) + C_%STABLE_DA = .false. + l = 1 + return + call dadeb !(31,'ERR DADAL ',1) + endif + if(idal(i).eq.nda_dab) then + ! deallocate nst0 = idapo(nda_dab) - 1 nda_dab = nda_dab - 1 else @@ -1000,6 +1367,7 @@ subroutine dadal(idal,l) enddo return end subroutine dadal + subroutine dadal1(idal) implicit none ! ************************ @@ -1010,6 +1378,34 @@ subroutine dadal1(idal) ! integer idal,ipause,mypauses ! +if(newtpsa) then + if(nomax.eq.1) then + if(idal.eq.nda_dab) then + ! deallocate + nda_dab = nda_dab - 1 + else + nhole=nhole+1 + endif + cc(idapo(idal):idapo(idal)+nmmax-1)=0 + allvec(idal) = .false. + idal= 0 + return + endif + + if(nomax==2) then + if(idal.eq.nda_dab) then + ! deallocate + nda_dab = nda_dab - 1 + else + nhole=nhole+1 + endif + cc(idapo(idal):idapo(idal)+nmmax-1)=0 + allvec(idal) = .false. + idal= 0 + return + endif +endif + ! if((.not.C_%STABLE_DA)) then ! if(c_%watch_user) then ! write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1040,6 +1436,7 @@ subroutine dadal1(idal) idal = 0 return end subroutine dadal1 + subroutine count_da(n) implicit none ! ************************ @@ -1056,6 +1453,7 @@ subroutine count_da(n) enddo return end subroutine count_da + subroutine davar(ina,ckon,i) implicit none ! **************************** @@ -1068,6 +1466,25 @@ subroutine davar(ina,ckon,i) integer i,ibase,ic1,ic2,illa,ilma,ina,inoa,inva,ipoa,ipause,mypauses real(dp) ckon ! + if(newtpsa) then + if(nomax.eq.1) then + cc(idapo(ina):idapo(ina)+nmmax-1) =0 + cc(idapo(ina)) = ckon + cc(idapo(ina)+i) = one + return + endif + + if(nomax==2) then + + cc(idapo(ina):idapo(ina)+nmmax-1) =0 + cc(idapo(ina)) = ckon + cc(idapo(ina)+i) = one + return + endif +endif + + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1138,7 +1555,20 @@ subroutine dacon(ina,ckon) ! integer illa,ilma,ina,inoa,inva,ipoa real(dp) ckon - ! + + if(newtpsa) then + if(nomax==1) then + call daclr(ina) + cc(idapo(ina)) = ckon + return + endif + if(nomax==2) then + call daclr(ina) + cc(idapo(ina)) = ckon + return + endif +endif + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1248,6 +1678,88 @@ subroutine dapek(ina,jv,cjj) integer,dimension(lnv)::jj real(dp) cjj ! + + icu=0 + do i=1,nvmax + icu=jv(i)+icu + enddo + ic1=0 + do i=nvmax+1,size(jv) + ic1=jv(i)+ic1 + enddo +if(icu>nomax.or.ic1>0) then + write(6,*) "in c_dabnew Crash 1, 0 continue ",icu,ic1 + read(5,*) i + if(i==1) then + cjj=sqrt(-cjj) + stop + endif + +! return +endif + if(newtpsa) then +ipoa=idapo(ina) + + if(nomax==1) then + + icu = 0 +jj1=0 + jj=0 + do i=1,size(jv) + jj(i)=jv(i) + enddo + do i=1,nvmax + icu = icu + jj(i) + if(jj(i).eq.1) jj1=i + enddo + + if(icu>1) then + cjj=0 + return + endif + + ipek = ipoa + jj1 + cjj = cc(ipek) + return + endif + + if(nomax==2) then + jj=0 + do i=1,size(jv) + jj(i)=jv(i) + enddo + ic1=0 + cjj=0 + do i=1,nvmax + ic1=ic1+jj(i) + enddo + if(ic1>nomax) return + ic2=0 + ic1=0 + cjj=0 + do i=1,nvmax + if(jj(i)==2) then + ic1=i + ic2=i + exit + endif + + if(jj(i)==1) then + if(ic1==0) then + ic1=i + else + ic2=i + exit + endif + endif + enddo + + ipek = ipoa + inds(ic1,ic2) - 1 + cjj = cc(ipek) + return + endif +endif + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1377,6 +1889,87 @@ subroutine dapok(ina,jv,cjj) integer,dimension(lnv)::jj real(dp) cjj ! + icu=0 + do i=1,size(jv) + icu=jv(i)+icu + enddo + ic1=0 + do i=nvmax+1,size(jv) + ic1=jv(i)+ic1 + enddo +if(icu>nomax.or.ic1>0) then + write(6,*) "in c_dabnew Crash 1, 0 continue ",icu,ic1 + read(5,*) i + if(i==1) then + cjj=sqrt(-cjj) + stop + endif + +! return +endif + if(newtpsa) then +ipoa=idapo(ina) + + if(nomax==1) then + icu = 0 + jj1=0 + jj=0 + do i=1,size(jv) + jj(i)=jv(i) + enddo + do i=1,nvmax + icu = icu + jj(i) + if(jj(i).eq.1) jj1=i + enddo + + if(icu>1) then + return + endif + + ipok = ipoa + jj1 + cc(ipok) = cjj + + return + endif + + if(nomax==2) then + jj=0 + do i=1,size(jv) + jj(i)=jv(i) + enddo + ic1=0 + do i=1,nvmax + ic1=ic1+jj(i) + enddo + if(ic1>nomax) return + ic2=0 + ic1=0 + do i=1,nvmax + if(jj(i)==2) then + ic1=i + ic2=i + exit + endif + + if(jj(i)==1) then + if(ic1==0) then + ic1=i + else + ic2=i + exit + endif + endif + enddo + + ipok = ipoa + inds(ic1,ic2) - 1 + cc(ipok)=cjj + + return + endif +endif + + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1537,6 +2130,15 @@ subroutine daclr(inc) ! integer i,illc,ilmc,inc,inoc,invc,ipoc ! + + if(newtpsa) then + + cc(idapo(inc):idapo(inc)+nmmax-1)=0 + return + +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1576,6 +2178,13 @@ subroutine dacop(ina,inb) ! integer ia,ib,illa,ina,inb,ipoa,ipob ! + + if(newtpsa) then + cc(idapo(inb):idapo(inb)+nmmax-1)=cc(idapo(ina):idapo(ina)+nmmax-1) + return +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1605,6 +2214,7 @@ subroutine dacop(ina,inb) idall(inb) = ib - ipob + 1 return end subroutine dacop + subroutine daadd(ina,inb,inc) implicit none ! ***************************** @@ -1618,6 +2228,17 @@ subroutine daadd(ina,inb,inc) integer idaadd,inb,inc,ipoc integer ipob ! + + if(newtpsa) then + ipoc = idapo(inc) + ipoa = idapo(ina) + ipob = idapo(inb) + cc(ipoc:ipoc+nmmax-1) = cc(ipoa:ipoa+nmmax-1) + cc(ipob:ipob+nmmax-1) + return + +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1649,7 +2270,40 @@ end subroutine daadd ! subroutine datrunc(ina,io,inb) implicit none - integer ina,io,inb,nt,ipoca,ipocb,i + integer ina,io,inb,nt,ipoca,ipocb,i,iot + + if(newtpsa) then + ipoca=idapo(ina) + ipocb=idapo(inb) + if(nomax==1) then + + do i=1,nvmax + cc(ipocb+i) =0.0_dp + enddo + + if(io==1) then + cc(ipocb) =cc(ipoca) + elseif(io>1) then + cc(ipocb:ipocb+nmmax-1)=cc(ipoca:ipoca+nmmax-1) + endif + return + endif + + if(nomax==2) then + cc(ipocb+1:ipocb+combien-1) = 0.0_dp + if(io==1) then + cc(ipocb) =cc(ipoca) + elseif(io==2) then + cc(ipocb:ipocb+nvmax) =cc(ipoca:ipoca+nvmax) + cc(ipocb+nvmax+1:ipocb+combien-1)=0.0_dp + elseif(io>2) then + cc(ipocb:ipocb+nmmax-1)=cc(ipoca:ipoca+nmmax-1) + endif + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1687,6 +2341,16 @@ subroutine dasub(ina,inb,inc) integer inc,ipoc,inb integer ipob ! + + if(newtpsa) then + ipoc = idapo(inc) + ipoa = idapo(ina) + ipob = idapo(inb) + cc(ipoc:ipoc+nmmax-1) = cc(ipoa:ipoa+nmmax-1) - cc(ipob:ipob+nmmax-1) + return +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1716,6 +2380,7 @@ subroutine dasub(ina,inb,inc) ! return end subroutine dasub + subroutine damul(ina,inb,inc) implicit none ! ***************************** @@ -1726,9 +2391,47 @@ subroutine damul(ina,inb,inc) ! !----------------------------------------------------------------------------- ! - integer ina,inb,inc,incc,ipoc,ipoa,ipob,i + integer ina,inb,inc,incc,ipoc,ipoa,ipob,i,j real(dp) ccipoa,ccipob ! + + if(newtpsa) then + ipoa=idapo(ina) + ipob=idapo(inb) + ipoc=idapo(inc) + if(nomax==1) then + ccipoa = cc(ipoa) + ccipob = cc(ipob) + cc(ipoc) = ccipoa*ccipob + do i=1,nvmax + cc(ipoc+i) = ccipoa*cc(ipob+i) + ccipob*cc(ipoa+i) + enddo + + return + endif + + if(nomax==2) then + ccipoa = cc(ipoa) + ccipob = cc(ipob) + reel=0 + reel(1) = ccipoa*ccipob + do i=1,nvmax + reel(i+1) = ccipoa*cc(ipob+i) + ccipob*cc(ipoa+i)+reel(i+1) + enddo + do i=poscombien,combien + reel(i) = ccipoa*cc(ipob+i-1) + ccipob*cc(ipoa+i-1) +reel(i) + enddo + do i=1,nvmax + do j=1,nvmax + reel(inds(i,j))=reel(inds(i,j))+cc(ipoa+i)*cc(ipob+j) + enddo + enddo + cc(ipoc:ipoc+combien-1)=reel + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1762,6 +2465,7 @@ subroutine damul(ina,inb,inc) endif return end subroutine damul + subroutine damult(ina,inb,inc) implicit none ! ***************************** @@ -1786,6 +2490,20 @@ subroutine damult(ina,inb,inc) integer,dimension(0:lno)::ipno,noff real(dp) ccia,ccipoa,ccipob ! + + if(newtpsa) then + if(nomax==1) then + stop 100 + return + endif + + if(nomax==2) then + stop 101 + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -1896,6 +2614,34 @@ subroutine dadiv(ina,inb,inc) ! integer idadiv,inb,ina,inc,ipoc,ipoa,ipob,i real(dp) ck,ck1 + + if(newtpsa) then + ipoa = idapo(ina) + ipob = idapo(inb) + ipoc = idapo(inc) + + if(nomax==1) then + ck=1.0_dp/cc(ipob) + ck1=cc(ipoa)*ck + do i=1,nvmax + cc(ipoc+i) = (cc(ipoa+i)-cc(ipob+i)*ck1)*ck + enddo + cc(ipoc)=ck1 + return + endif + + if(nomax==2) then + idadiv = 0 + ! call dainf(inc,inoc,invc,ipoc,ilmc,illc) + call daall1(idadiv,'$$DADIV $$',nomax,nvmax) + call dafun('INV ',inb,idadiv) + call damul(ina,idadiv,inc) + call dadal1(idadiv) + return + endif +endif + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -1936,23 +2682,45 @@ subroutine dasqr(ina,inc) ! integer ina,inc,incc,ipoc,i,ipoa real(dp) ccipoa - ! - ! CASE OF FIRST ORDER ONLY - ! ************************ - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - if(nomax.eq.1) then + + + if(newtpsa) then + if(nomax==1) then ipoc = idapo(inc) ipoa = idapo(ina) ! minv = min(inva,invc) ccipoa = cc(ipoa) cc(ipoc) = ccipoa*ccipoa do i=1,nvmax - cc(ipoc+i) = two*ccipoa*cc(ipoa+i) + cc(ipoc+i) = 2.0_dp*ccipoa*cc(ipoa+i) + enddo + return + endif + + if(nomax==2) then + call damul(ina,ina,inc) + return + endif +endif + + + ! + ! CASE OF FIRST ORDER ONLY + ! ************************ + if((.not.C_%STABLE_DA)) then + if(c_%watch_user) then + write(6,*) "big problem in dabnew ", sqrt(crash) + endif + return + endif + if(nomax.eq.1) then + ipoc = idapo(inc) + ipoa = idapo(ina) + ! minv = min(inva,invc) + ccipoa = cc(ipoa) + cc(ipoc) = ccipoa*ccipoa + do i=1,nvmax + cc(ipoc+i) = two*ccipoa*cc(ipoa+i) enddo return endif @@ -1968,6 +2736,7 @@ subroutine dasqr(ina,inc) endif return end subroutine dasqr + subroutine dasqrt(ina,inc) implicit none ! ************************* @@ -1989,6 +2758,21 @@ subroutine dasqrt(ina,inc) inva,invc,ioffa,ioffb,ipoa,ipoc,ipos,noia,noib,nom integer,dimension(0:lno)::ipno,noff real(dp) ccia,ccipoa + + + if(newtpsa) then + if(nomax==1) then +stop 350 + return + endif + + if(nomax==2) then +stop 351 + return + endif +endif + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2103,6 +2887,17 @@ subroutine dacad(ina,ckon,inb) integer,parameter,dimension(lnv)::jjx=0 real(dp) ckon,const ! + + if(newtpsa) then + call dacop(ina,inb) + + cc(idapo(inb)) = cc(idapo(inb)) + ckon + + return + +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -2132,6 +2927,18 @@ subroutine dacsu(ina,ckon,inb) integer ina,inb integer,parameter,dimension(lnv)::jjx=0 real(dp) ckon,const + + + if(newtpsa) then + call dacop(ina,inb) + + cc(idapo(inb)) = cc(idapo(inb)) - ckon + + return + +endif + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2166,6 +2973,15 @@ subroutine dasuc(ina,ckon,inb) ! integer i,ina,inb,ipoa,ipob real(dp) ckon + + if(newtpsa) then + ipob=idapo(inb) + ipoa=idapo(ina) + cc(ipob) = ckon - cc(ipoa) + cc(ipob+1:ipob+nmmax-1) =-cc(ipoa+1:ipoa+nmmax-1) + return + endif + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2200,6 +3016,16 @@ subroutine dacmu(ina,ckon,inc) ! integer ipoa,i,ina,inc,incc,ipoc real(dp) ckon + + + if(newtpsa) then + ipoc=idapo(inc) + ipoa=idapo(ina) + cc(ipoc:ipoc+nmmax-1) = ckon*cc(ipoa:ipoa+nmmax-1) + return + endif + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2228,6 +3054,7 @@ subroutine dacmu(ina,ckon,inc) endif return end subroutine dacmu + subroutine dacmut(ina,ckon,inb) implicit none ! ****************************** @@ -2246,6 +3073,22 @@ subroutine dacmut(ina,ckon,inb) integer i,ia,ib,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,ipoa,& ipob,ipause,mypauses real(dp) ckon + + + if(newtpsa) then + if(nomax==1) then + stop 900 + return + endif + + if(nomax==2) then + stop 901 + + return + endif +endif + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2308,6 +3151,17 @@ subroutine dacdi(ina,ckon,inb) ! integer i,ina,inb,ipoa,ipob,ipause,mypauses real(dp) ckon + + + if(newtpsa) then + ipob=idapo(inb) + ipoa=idapo(ina) + cc(ipob:ipob+nmmax-1) = cc(ipoa:ipoa+nmmax-1)/ckon + return + endif + + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2351,6 +3205,31 @@ subroutine dadic(ina,ckon,inc) ! integer i,idadic,ina,inc,ipoa,ipoc real(dp) ckon,ck + + + if(newtpsa) then + if(nomax==1) then + ipoa = idapo(ina) + ipoc = idapo(inc) + cc(ipoc)=ckon/cc(ipoa) + ck=cc(ipoc)/cc(ipoa) + do i=1,nvmax + cc(ipoc+i) = -cc(ipoa+i)* ck + enddo + return + endif + + if(nomax==2) then + idadic = 0 + call daall1(idadic,'$$DADIC $$',nomax,nvmax) + call dafun('INV ',ina,idadic) + call dacmu(idadic,ckon,inc) + call dadal1(idadic) + return + endif +endif + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2406,6 +3285,18 @@ subroutine dacma(ina,inb,bfac,inc) ! integer idacma,ina,inb,inc,ipoc,ipob,ipoa,i real(dp) bfac + + if(newtpsa) then + ipoc = idapo(inc) + ipoa = idapo(ina) + ipob = idapo(inb) + + cc(ipoc:ipoc+nmmax-1) = cc(ipoa:ipoa+nmmax-1) + cc(ipob:ipob+nmmax-1) * bfac + return + +endif + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2447,6 +3338,20 @@ subroutine dalin(ina,afac,inb,bfac,inc) !----------------------------------------------------------------------------- ! integer ipob,ipoa,i + + + if(newtpsa) then + ipoc = idapo(inc) + ipoa = idapo(ina) + ipob = idapo(inb) + + cc(ipoc:ipoc+nmmax-1) = cc(ipoa:ipoa+nmmax-1) * afac + cc(ipob:ipob+nmmax-1) * bfac + + return + endif + + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2478,6 +3383,7 @@ subroutine dalin(ina,afac,inb,bfac,inc) endif return end subroutine dalin + subroutine dalint(ina,afac,inb,bfac,inc) implicit none ! *************************************** @@ -2495,6 +3401,21 @@ subroutine dalint(ina,afac,inb,bfac,inc) inob,inoc,inva,invb,invc,ipoa,ipob,ipoc,is,ismax,ismin,ja,jb,mchk,& ipause,mypauses real(dp) afac,bfac,ccc,copf + + + if(newtpsa) then + if(nomax==1) then + stop 400 + return + endif + + if(nomax==2) then + stop 401 + return + endif +endif + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2668,6 +3589,9 @@ subroutine dafun(cf,ina,inc) ! integer ina,inc,incc character(4) cf + + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2687,6 +3611,7 @@ subroutine dafun(cf,ina,inc) endif return end subroutine dafun + subroutine dafunt(cf,ina,inc) implicit none ! **************************** @@ -2708,6 +3633,9 @@ subroutine dafunt(cf,ina,inc) data abcs /'abcdefghijklmnopqrstuvwxyz'/ data abcc /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/ ! + + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -2946,12 +3874,18 @@ subroutine dafunt(cf,ina,inc) call dacop(ina,inon) call dapok(inon,jjx,zero) call dacon(ipow,one) - ! + do i=1,min(no,nocut) + ! call damul(inon,ipow,iscr) + + call dacop(iscr,ipow) + + call dacma(inc,ipow,xf(i),inc) + ! enddo ! @@ -2974,6 +3908,22 @@ subroutine daabs(ina,anorm) ! integer i,illa,ilma,ina,inoa,inva,ipoa real(dp) anorm + + + if(newtpsa) then + ipoa=idapo(ina) + anorm = 0.0_dp + if(nomax==1) illa=nmmax+1 + if(nomax==2) illa=combien + do i=ipoa,ipoa+illa-1 + anorm = anorm + abs(cc(i)) + enddo + return + endif + + + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -2997,8 +3947,8 @@ subroutine daabs(ina,anorm) return end subroutine daabs ! - ! - subroutine dacct(ma,ia,mb,ib,mc,ic) + + subroutine dacctt1(mb,ib,mc,ic,ma,ia) implicit none ! *********************************** ! @@ -3008,39 +3958,66 @@ subroutine dacct(ma,ia,mb,ib,mc,ic) ! !----------------------------------------------------------------------------- ! - integer i,ia,ib,ic,ij,illc,ilmc,inoc,invc,ipoc - integer,dimension(lnv)::monx - integer,dimension(:)::ma,mb,mc - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif + ! INTEGER MON(c_lno+1),Ic_cc(c_lnv) + ! INTEGER,dimension(:)::MB,MC,MA + !ETIENNE + !----------------------------------------------------------------------------- ! - if(ma(1).eq.mc(1).or.mb(1).eq.mc(1)) then - call dainf(mc(1),inoc,invc,ipoc,ilmc,illc) - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - do ij=1,ic - monx(ij)=0 - enddo - call daall(monx,ic,'$$DAJUNK$$',inoc,invc) - call dacctt(ma,ia,mb,ib,monx,ic) - do i=1,ic - call dacop(monx(i),mc(i)) - enddo - call dadal(monx,ic) - else - call dacctt(ma,ia,mb,ib,mc,ic) - endif + integer i0,i,j,k,ia,ib,ic,illa,illb,illc,ilma,ilmb,ilmc,inoa,inob,inoc,inva,invb,invc + + ! integer,dimension(c_lno)::icc + integer,dimension(:)::ma,mb,mc + integer ipoa(lnv),ipob(lnv),ipoc(lnv) + ! + !ETIENNE + ! + ! CONSISTENCY CHECKS + ! ****************** + ! + + + if(ib/=ic) then + write(6,*) " stopped in dacctt1 " + stop 880 + endif + +call dainf(ma(1),inoa,i0,ipoa(1),ilma,illa) + + + do i=1,i0 + call dainf(ma(i),inoa,inva,ipoa(i),ilma,illa) + call dainf(mb(i),inob,invb,ipob(i),ilmb,illb) + call dainf(mc(i),inoc,invc,ipoc(i),ilmc,illc) + + cc(ipoa(i):ipoa(i)+inva)=0.0_dp + enddo + + + do i=1,inva + cc(ipoa(i))=cc(ipob(i)) + enddo + + do i=1,inva + do j=1,inva + cc(ipoa(i))=cc(ipoa(i))+cc(ipob(i)+j)*cc(ipoc(j)) + enddo + + enddo + do i=1,inva + do j=1,inva + do k=1,inva + cc(ipoa(i)+k)=cc(ipoa(i)+k)+cc(ipob(i)+j)*cc(ipoc(j)+k) + enddo + enddo + enddo + + return - end subroutine dacct - subroutine dacctt(mb,ib,mc,ic,ma,ia) + end subroutine dacctt1 + + + + subroutine dacctt2da(mb,ib,mc,ic,ma,ia) implicit none ! *********************************** ! @@ -3050,49 +4027,429 @@ subroutine dacctt(mb,ib,mc,ic,ma,ia) ! !----------------------------------------------------------------------------- ! - ! INTEGER MON(LNO+1),ICC(LNV) + ! INTEGER MON(c_lno+1),Icc(lnv) ! INTEGER,dimension(:)::MB,MC,MA !ETIENNE !----------------------------------------------------------------------------- ! - integer i,ia,ib,ic,iia,iib,iic,illa,illb,illc,ilma,ilmb,ilmc,inoa,inob,inoc,inva,invb,invc,& - ipoa,ipob,ipoc,iv,jl,jv,ipause,mypauses - integer,dimension(lno+1)::mon - ! integer,dimension(lno)::icc - integer,dimension(lnv)::icc !johan 2008 March - integer,dimension(:)::ma,mb,mc - real(dp) ccf - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif + integer i,j,k,ia,ib,ic,illa,illb,illc,ilma,ilmb,ilmc,inoa,inob,inoc,inva,invb,invc + integer ab,jk,ik + + ! integer,dimension(c_lno)::icc + integer,dimension(:)::ma,mb,mc + integer ipoa(lnv),ipob(lnv),ipoc(lnv) + real(dp) xx ! !ETIENNE ! ! CONSISTENCY CHECKS ! ****************** ! - iia = ma(1) - iib = mb(1) - iic = mc(1) - call dainf(iia,inoa,inva,ipoa,ilma,illa) - call dainf(iib,inob,invb,ipob,ilmb,illb) - call dainf(iic,inoc,invc,ipoc,ilmc,illc) - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - ! - call damch(ma,ia) - call damch(mb,ib) - ! - if(ia.ne.ib) then - write(line,'(a26)') 'ERROR IN DACCT, IA .NE. IB' - ipause=mypauses(35,line) + + + if(ib/=ic) then + write(6,*) " stopped in dacctt1 " + stop 880 + endif + + + + do i=1,ia + call dainf(ma(i),inoa,inva,ipoa(i),ilma,illa) + cc(ipoa(i):ipoa(i)+combien-1)=0.0_dp + enddo + + do i=1,ib + call dainf(mb(i),inob,invb,ipob(i),ilmb,illb) + enddo + + do i=1,ic + call dainf(mc(i),inoc,invc,ipoc(i),ilmc,illc) + enddo + + + do i=invc-nphere+1,invc + cc(ipoa(i)+i)=1.0_dp + enddo +!cc(ipoa(i)+inva-nphere+1:ipoa(i)+inva)=1.0_dp + + do i=1,invc-nphere + do j=1,invc-nphere !? + do k=1,invc + cc(ipoa(i)+k)=cc(ipoa(i)+k)+cc(ipob(i)+j)*cc(ipoc(j)+k) + enddo + enddo + + do j=invc-nphere+1,invc +! do k=1,inva + cc(ipoa(i)+j)=cc(ipoa(i)+j)+cc(ipob(i)+j)*cc(ipoc(j)+j) +! enddo + enddo + + enddo + + + do i=1,invc-nphere + do j=1,invc -nphere !? + do jk=poscombien ,combien + cc(ipoa(i)+jk-1 )=cc(ipoa(i)+jk-1 ) + cc(ipob(i)+j)*cc(ipoc(j)+jk-1) + enddo + enddo + enddo + +if(with_para==0) then + do i=1,invc-nphere + do jk=poscombien,combien + do ab=poscombien,combien + cc(ipoa(i)+ab-1 )=cc(ipoa(i)+ab-1 ) + finds(ind1(jk),ind2(jk),ab)* & + cc(ipob(i)+jk-1 )*(cc(ipoc(ind1(jk))+ind1(ab))*cc(ipoc(ind2(jk)) & + +ind2(ab))+cc(ipoc(ind1(jk))+ind2(ab))*cc(ipoc(ind2(jk))+ind1(ab)))/2.0_dp + enddo + enddo + enddo + elseif(with_para==2) then + do i=1,invc-nphere + do ik=1,ninds + ab=nind1(ik) + jk=nind2(ik) + cc(ipoa(i)+ab-1 )=cc(ipoa(i)+ab-1 ) + finds1(ik)* & + cc(ipob(i)+jk-1 )*(cc(ipoc(ind1(jk))+ind1(ab))*cc(ipoc(ind2(jk))+ind2(ab)) & + +cc(ipoc(ind1(jk))+ind2(ab))*cc(ipoc(ind2(jk))+ind1(ab)))/2.0_dp + + enddo + enddo +elseif(with_para==1) then + do i=1,invc-nphere + do ik=1,ninds + ab=nind1(ik) + jk=nind2(ik) + cc(ipoa(i)+ab-1 )=cc(ipoa(i)+ab-1 ) + finds(ind1(jk),ind2(jk),ab)* & + cc(ipob(i)+jk-1 )*(cc(ipoc(ind1(jk))+ind1(ab))*cc(ipoc(ind2(jk))+ind2(ab)) & + +cc(ipoc(ind1(jk))+ind2(ab))*cc(ipoc(ind2(jk))+ind1(ab)))/2.0_dp + + enddo + enddo +endif + +do i=1,inva + cc(ipoa(i))=cc(ipob(i)) +enddo + + return + end subroutine dacctt2da + + + + subroutine dacctt2tpsa(mb,ib,mc,ic,ma,ia) + implicit none + ! *********************************** + ! + ! THIS SUBROUTINE PERFORMS A CONCATENATION MA = MB o MC + ! WHERE MA, MB AND MC ARE MATRICES CONSISTING OF IA, IB AND IC + ! DA VECTORS, RESPECTIVELY. + ! + !----------------------------------------------------------------------------- + ! + ! INTEGER MON(c_lno+1),Ic_cc(c_lnv) + ! INTEGER,dimension(:)::MB,MC,MA + !ETIENNE + !----------------------------------------------------------------------------- + ! + integer i0,i,j,k,ia,ib,ic,illa,illb,illc,ilma,ilmb,ilmc,inoa,inob,inoc,inva,invb,invc,m,n + real(dp) fac + + ! integer,dimension(c_lno)::icc + integer,dimension(:)::ma,mb,mc + integer ipoa(lnv),ipob(lnv),ipoc(lnv) + real(dp), allocatable :: g0(:),h0(:), G(:,:),H(:,:), Ga(:,:,:), He(:,:,:) + real(dp), allocatable :: f0(:),F(:,:),Fi(:,:,:) + + ! + !ETIENNE + ! + ! CONSISTENCY CHECKS + ! ****************** + ! +call dainf(ma(1),inoa,i0,ipoa(1),ilma,illa) + +allocate(g0(i0),h0(i0),f0(i0), G(i0,i0),H(i0,i0),F(i0,i0), Ga(i0,i0,i0), He(i0,i0,i0), fi(i0,i0,i0)) + g0=0 + h0=0 + f0=0 + g=0 + h=0 + f=0 + ga=0 + he=0 + fi=0 + + do i=1,i0 + call dainf(ma(i),inoa,inva,ipoa(i),ilma,illa) + call dainf(mb(i),inob,invb,ipob(i),ilmb,illb) + call dainf(mc(i),inoc,invc,ipoc(i),ilmc,illc) + if(inva/=i0) stop 1234 + if(invb/=i0) stop 1235 + if(invc/=i0) stop 1236 + cc(ipoa(i):ipoa(i)+nmmax-1)=0.0_dp + + g0(i)=cc(ipob(i)) + h0(i)=cc(ipoc(i)) + do j=1,i0 + g(i,j)=cc(ipob(i)+j) + h(i,j)=cc(ipoc(i)+j) + do k=1,i0 + fac=2.0_dp + if(j==k) fac=1.0_dp + ga(i,j,k)=cc(ipob(i)+inds(j,k)-1)/fac + he(i,j,k)=cc(ipoc(i)+inds(j,k)-1)/fac + enddo + enddo + enddo + + f0=g0 + + do i=1,i0 + do j=1,i0 + f0(i)=g(i,j)*h0(j) +f0(i) + + do k=1,i0 + f0(i)=ga(i,j,k)*h0(j)*h0(k) + f0(i) + f(i,k)=f(i,k)+g(i,j)*H(j,k) + do m=1,i0 + f(i,k)=f(i,k)+ga(i,j,m)*H(j,k)*h0(m)+ga(i,j,m)*H(m,k)*h0(j) + enddo + enddo + enddo + enddo + + + do i=1,i0 + do j=1,i0 + do n=1,i0 + do m=1,i0 + fi(i,m,n)=fi(i,m,n)+ g(i,j)*he(j,m,n) + do k=1,i0 + fi(i,m,n)=fi(i,m,n)+ ga(i,j,k)*he(k,m,n)*h0(j)+ ga(i,j,k)*he(j,m,n)*h0(k) + fi(i,m,n)=fi(i,m,n)+ ga(i,j,k)*h(j,m)*H(k,n) + enddo + enddo + enddo + enddo + enddo + + + do i=1,i0 + cc(ipoa(i))=f0(i) + cc(ipoa(i)+1:ipoa(i)+i0)=f(i,1:i0) + enddo + + + do k=1,i0 + do j=poscombien, combien + fac=2.0_dp + if(ind1(j)==ind2(j)) fac=1.0_dp + cc(ipoa(k)+j-1)=fi(k,ind1(j),ind2(j))*fac + enddo + enddo + + + +deallocate(g0,h0, G,H, Ga, He,f0,f,fi) + + return + end subroutine dacctt2tpsa + + subroutine dacct(ma,ia,mb,ib,mc,ic) + implicit none + ! *********************************** + ! + ! THIS SUBROUTINE PERFORMS A CONCATENATION MA = MB o MC + ! WHERE MA, MB AND MC ARE MATRICES CONSISTING OF IA, IB AND IC + ! DA VECTORS, RESPECTIVELY. + ! + !----------------------------------------------------------------------------- + ! + integer i,ia,ib,ic,ij,illc,ilmc,inoc,invc,ipoc + integer,dimension(lnv)::monx,maa,mcc + integer,dimension(:)::ma,mb,mc + integer iaa,icc + + if(newtpsa) then + if(ib/=nvmax) then + write(6,*) "problem in dacct " + stop 9 + endif + if(ia/=ic) then + write(6,*) "problem in dacct " + stop 10 + endif + maa=0 + mcc=0 + do i=1,ia + maa(i)=ma(i) + enddo + + if((mc(1)==ma(1)).or.(mc(1)==mb(1))) then + do i=1,ic + call daall0(mcc(i)) + enddo + else + do i=1,ic + mcc(i)=mc(i) + enddo + endif + + + do i=ia+1,nvmax + call daall0(maa(i)) + enddo + do i=ic+1,nvmax + call daall0(mcc(i)) + enddo + + + + + if(nomax==1) then + call dacctt1(maa,nvmax,mb,nvmax,mcc,nvmax) + + if((mc(1)==ma(1)).or.(mc(1)==mb(1))) then + do i=1,ic + call dacop(mcc(i),mc(i)) + call dadal1(mcc(i)) + enddo + endif + + do i=ia+1,nvmax + call dadal1(maa(i)) + enddo + do i=ic+1,nvmax + call dadal1(mcc(i)) + enddo + + return + endif + + if(nomax==2) then + if(assume_da_map) then + call dacctt2da(maa,nvmax,mb,nvmax,mcc,nvmax) + else + call dacctt2tpsa(maa,nvmax,mb,nvmax,mcc,nvmax) + endif + if((mc(1)==ma(1)).or.(mc(1)==mb(1))) then + do i=1,ic + call dacop(mcc(i),mc(i)) + call dadal1(mcc(i)) + enddo + endif + do i=ia+1,nvmax + call dadal1(maa(i)) + enddo + do i=ic+1,nvmax + call dadal1(mcc(i)) + enddo + return + endif +endif + + + if((.not.C_%STABLE_DA)) then + if(c_%watch_user) then + write(6,*) "big problem in dabnew ", sqrt(crash) + endif + return + endif + ! + if(ma(1).eq.mc(1).or.mb(1).eq.mc(1)) then + call dainf(mc(1),inoc,invc,ipoc,ilmc,illc) + if((.not.C_%STABLE_DA)) then + if(c_%watch_user) then + write(6,*) "big problem in dabnew ", sqrt(crash) + endif + return + endif + do ij=1,ic + monx(ij)=0 + enddo + call daall(monx,ic,'$$DAJUNK$$',inoc,invc) + call dacctt(ma,ia,mb,ib,monx,ic) + do i=1,ic + call dacop(monx(i),mc(i)) + enddo + call dadal(monx,ic) + else + call dacctt(ma,ia,mb,ib,mc,ic) + endif + return + end subroutine dacct + + subroutine dacctt(mb,ib,mc,ic,ma,ia) + implicit none + ! *********************************** + ! + ! THIS SUBROUTINE PERFORMS A CONCATENATION MA = MB o MC + ! WHERE MA, MB AND MC ARE MATRICES CONSISTING OF IA, IB AND IC + ! DA VECTORS, RESPECTIVELY. + ! + !----------------------------------------------------------------------------- + ! + ! INTEGER MON(LNO+1),ICC(LNV) + ! INTEGER,dimension(:)::MB,MC,MA + !ETIENNE + !----------------------------------------------------------------------------- + ! + integer i,ia,ib,ic,iia,iib,iic,illa,illb,illc,ilma,ilmb,ilmc,inoa,inob,inoc,inva,invb,invc,& + ipoa,ipob,ipoc,iv,jl,jv,ipause,mypauses + integer,dimension(lno+1)::mon + ! integer,dimension(lno)::icc + integer,dimension(lnv)::icc !johan 2008 March + integer,dimension(:)::ma,mb,mc + real(dp) ccf + + if(newtpsa) then + if(nomax==1) then + stop 750 + return + endif + + if(nomax==2) then + stop 751 + return + endif +endif + + + if((.not.C_%STABLE_DA)) then + if(c_%watch_user) then + write(6,*) "big problem in dabnew ", sqrt(crash) + endif + return + endif + ! + !ETIENNE + ! + ! CONSISTENCY CHECKS + ! ****************** + ! + iia = ma(1) + iib = mb(1) + iic = mc(1) + call dainf(iia,inoa,inva,ipoa,ilma,illa) + call dainf(iib,inob,invb,ipob,ilmb,illb) + call dainf(iic,inoc,invc,ipoc,ilmc,illc) + if((.not.C_%STABLE_DA)) then + if(c_%watch_user) then + write(6,*) "big problem in dabnew ", sqrt(crash) + endif + return + endif + ! + call damch(ma,ia) + call damch(mb,ib) + ! + if(ia.ne.ib) then + write(line,'(a26)') 'ERROR IN DACCT, IA .NE. IB' + ipause=mypauses(35,line) call dadeb !(31,'ERR DACCT1',1) elseif(ic.ne.invb) then write(line,'(a26)') 'ERROR IN DACCT, IC.NE.INVB' @@ -3163,6 +4520,21 @@ subroutine mtree(mb,ib,mc,ic) integer,dimension(0:lno)::jv integer,dimension(:)::mb,mc real(dp) apek,bbijj,chkjj + + + if(newtpsa) then + if(nomax==1) then + stop 888 + return + endif + + if(nomax==2) then + stop 889 + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -3370,6 +4742,20 @@ subroutine ppushprint(mc,ic,mf,jc,line) integer i,ic,iv,jc,jl,jv,mf integer,dimension(:)::mc character(20) line + + if(newtpsa) then + if(nomax==1) then + stop 788 + return + endif + + if(nomax==2) then + stop 789 + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -3407,6 +4793,20 @@ subroutine ppushstore(mc,nd2,coef,ml,mv) integer,dimension(:), intent(out)::ml,mv real(dp),dimension(:),intent(out)::coef ! + + if(newtpsa) then + if(nomax==1) then + stop 818 + return + endif + + if(nomax==2) then + stop 819 + return + endif +endif + + jc=0 jl=0;jv=0; ic=nd2 @@ -3434,6 +4834,7 @@ subroutine ppushstore(mc,nd2,coef,ml,mv) enddo return end subroutine ppushstore + subroutine ppushGETN(mc,ND2,ntot) implicit none ! @@ -3458,6 +4859,20 @@ subroutine ppush(mc,ic,xi,xf) real(dp),dimension(lno+1)::xm real(dp),dimension(lno)::xt real(dp),dimension(:)::xf,xi + + if(newtpsa) then + if(nomax==1) then + stop 882 + return + endif + + if(nomax==2) then + stop 883 + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -3494,6 +4909,7 @@ subroutine ppush(mc,ic,xi,xf) ! return end subroutine ppush + subroutine ppush1(mc,xi,xf) implicit none ! ***************************** @@ -3508,6 +4924,21 @@ subroutine ppush1(mc,xi,xf) real(dp),dimension(:)::xi real(dp),dimension(lno)::xt real(dp),dimension(lno+1)::xm + + + if(newtpsa) then + if(nomax==1) then + stop 838 + return + endif + + if(nomax==2) then + stop 839 + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -3538,24 +4969,106 @@ subroutine ppush1(mc,xi,xf) endif return end subroutine ppush1 - subroutine dainv(ma,ia,mb,ib) + + + +subroutine dainvt1(mb,ib,mc,ic) implicit none - ! ***************************** - ! - ! THIS SUBROUTINE INVERTS THE MATRIX MA WITH IA DA VECTORS AND - ! STORES THE RESULT IN MI + ! *********************************** ! + ! THIS SUBROUTINE PERFORMS A INVERSION first or second order + !ETIENNE !----------------------------------------------------------------------------- ! - integer i,ia,ib,ij,illb,ilmb,inob,invb,ipob - integer,dimension(lnv)::jj,ml - integer,dimension(:)::ma,mb - real(dp),dimension(lnv)::x - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return + integer i0,i,j,k,ib,ic,illa,illb,illc,ilmb,ilmc,inob,inoc,invb,invc,m,n + real(dp) fac + + ! integer,dimension(lno)::icc + integer,dimension(:):: mb,mc + integer ipob(lnv),ipoc(lnv),ier + real(dp), allocatable :: g0(:),h0(:), G(:,:),H(:,:) + + integer, allocatable :: icc(:),ipo(:) + ! + !ETIENNE + ! + ! CONSISTENCY CHECKS + ! ****************** + ! + + +call dainf(mb(1),inob,i0,ipob(1),ilmb,illb) + + +allocate(g0(i0),h0(i0) , G(i0,i0),H(i0,i0) ) + g0=0 + h0=0 + g=0 + h=0 + + + do i=1,i0 + call dainf(mc(i),inoc,invc,ipoc(i),ilmc,illc) + call dainf(mb(i),inob,invb,ipob(i),ilmb,illb) + if(invb/=i0) stop 1235 + if(invc/=i0) stop 1236 + + g0(i)=cc(ipob(i)) + do j=1,i0 + g(i,j)=cc(ipob(i)+j) + + enddo + enddo + + call matinv(g,h,i0,i0,ier) + + h0=-matmul(h,g0) + + + do i=1,i0 + cc(ipoc(i))=h0(i) + cc(ipoc(i)+1:ipoc(i)+i0)=H(i,1:i0) + enddo + + + +end subroutine dainvt1 + + subroutine dainv(ma,ia,mb,ib) + implicit none + ! ***************************** + ! + ! THIS SUBROUTINE INVERTS THE MATRIX MA WITH IA DA VECTORS AND + ! STORES THE RESULT IN MI + ! + !----------------------------------------------------------------------------- + ! + integer i,ia,ib,ij,illb,ilmb,inob,invb,ipob + integer,dimension(lnv)::jj,ml + integer,dimension(:)::ma,mb + real(dp),dimension(lnv)::x + + + + if(newtpsa) then + if(nomax==1) then + call dainvt1(ma,ia,mb,ib) + return + endif + + if(nomax==2) then + call dainvt2(ma,ia,mb,ib) + return + endif +endif + + + + if((.not.C_%STABLE_DA)) then + if(c_%watch_user) then + write(6,*) "big problem in dabnew ", sqrt(crash) + endif + return endif ! do i=1,lnv @@ -3593,6 +5106,135 @@ subroutine dainv(ma,ia,mb,ib) endif return end subroutine dainv + + +subroutine dainvt2(mb,ib,mc,ic) + implicit none + ! *********************************** + ! + ! THIS SUBROUTINE PERFORMS A INVERSION first or second order + !ETIENNE + !----------------------------------------------------------------------------- + ! + integer i0,i,j,k,ib,ic,illa,illb,illc,ilmb,ilmc,inob,inoc,invb,invc,m,n + real(dp) fac + + ! integer,dimension(lno)::mcc + integer,dimension(:):: mb,mc + integer ipob(lnv),ipoc(lnv),ier + real(dp), allocatable :: g0(:),h0(:), G(:,:),H(:,:), Ga(:,:,:), He(:,:,:),r0(:), r1(:,:) + + integer, allocatable :: mcc(:),ipo(:),mccc(:),ipoo(:) + ! + !ETIENNE + ! + ! CONSISTENCY CHECKS + ! ****************** + ! + +call dainf(mb(1),inob,i0,ipob(1),ilmb,illb) + + +allocate(r0(i0),g0(i0),h0(i0) ,r1(i0,i0), G(i0,i0),H(i0,i0) , Ga(i0,i0,i0), He(i0,i0,i0) ) + g0=0 + h0=0 + g=0 + h=0 + ga=0 + he=0 + + do i=1,i0 + call dainf(mc(i),inoc,invc,ipoc(i),ilmc,illc) + call dainf(mb(i),inob,invb,ipob(i),ilmb,illb) + if(invb/=i0) stop 1235 + if(invc/=i0) stop 1236 + + g0(i)=cc(ipob(i)) + do j=1,i0 + g(i,j)=cc(ipob(i)+j) + if(inob==2) then + do k=1,i0 + fac=2.0_dp + if(j==k) fac=1.0_dp + ga(i,j,k)=cc(ipob(i)+inds(j,k)-1)/fac + enddo + endif + enddo + enddo + + call matinv(g,h,i0,i0,ier) + + ! h0=-matmul(h,g0) + + + + allocate(mcc(i0),ipo(i0)) + allocate(mccc(i0),ipoo(i0)) + + call daall(mcc,i0,'$$DACCT $$',nomax,nvmax) + call daall(mccc,i0,'$$DACCT $$',nomax,nvmax) + + do i=1,i0 + call dainf(mcc(i),inoc,invc,ipo(i),ilmc,illc) + call dainf(mccc(i),inoc,invc,ipoo(i),ilmc,illc) + cc(ipo(i):ipo(i)+combien-1)=0 + cc(ipo(i))=h0(i) + cc(ipo(i)+1:ipo(i)+i0)=H(i,1:i0) + enddo + + if(assume_da_map) then + call dacctt2da(mb,i0,mcc,i0,mccc,i0) + else + call dacctt2tpsa(mb,i0,mcc,i0,mccc,i0) + endif + + do i=1,i0 + cc(ipoo(i)+poscombien-1:ipoo(i)+combien-1)=-cc(ipoo(i)+poscombien-1:ipoo(i)+combien-1) + enddo + + + if(assume_da_map) then + call dacctt2da(mcc,i0,mccc,i0,mc,ic) + else + call dacctt2tpsa(mcc,i0,mccc,i0,mc,ic) + endif + + do i=1,i0 + cc(idapo(mcc(i)):idapo(mcc(i))+nmmax-1 )=0 + cc(idapo(mccc(i)):idapo(mccc(i))+nmmax-1 )=0 + enddo + + do i=1,i0 + cc(idapo(mc(i)):idapo(mc(i)) ) = 0.0_dp + cc(idapo(mcc(i)):idapo(mcc(i)) ) = -g0(i) + cc(idapo(mcc(i))+i)=1.0_dp + enddo + + if(assume_da_map) then + call dacctt2da(mc,ic,mcc,i0,mccc,i0) + else + call dacctt2tpsa(mc,ic,mcc,i0,mccc,i0) + endif + + do i=1,i0 + cc(idapo(mc(i)):idapo(mc(i)) + nmmax-1 ) = cc(idapo(mccc(i)):idapo(mccc(i)) + nmmax-1 ) + enddo + + call dadal(mcc,i0) + call dadal(mccc,i0) + deallocate(mcc) + deallocate(mccc) + deallocate(ipo) + deallocate(ipoo) + +deallocate(g0,h0, G,H, Ga, He) + + + +end subroutine dainvt2 + + + subroutine dainvt(ma,ia,mb,ib) implicit none ! ***************************** @@ -3608,6 +5250,21 @@ subroutine dainvt(ma,ia,mb,ib) integer,dimension(:)::ma,mb real(dp),dimension(lnv,lnv)::aa,ai real(dp) amjj,amsjj,prod + + + if(newtpsa) then + if(nomax==1) then + stop 780 + return + endif + + if(nomax==2) then + stop 781 + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -3661,7 +5318,8 @@ subroutine dainvt(ma,ia,mb,ib) ! do i=1,ib do j=1,ib - do k=1,ib +! do k=1,ib + do k=1,size(jj) jj(k) = 0 enddo jj(j) = 1 @@ -3789,6 +5447,7 @@ subroutine dapin(ma,ia,mb,ib,jx) integer,dimension(:)::ma,mb,jx real(dp),dimension(lnv)::x ! + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -3830,6 +5489,7 @@ subroutine dapin(ma,ia,mb,ib,jx) endif return end subroutine dapin + subroutine dapint(ma,ia,mb,ib,jind) implicit none ! ********************************** @@ -3843,6 +5503,10 @@ subroutine dapint(ma,ia,mb,ib,jind) integer,dimension(lnv)::jj,mn,mi,me integer,dimension(:)::ma,mb,jind ! + + + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -3867,8 +5531,9 @@ subroutine dapint(ma,ia,mb,ib,jind) call daall(mi,ia,'$$PIN2 $$',inoa,inva) call daall(me,ia,'$$PIN3 $$',inoa,inva) ! + do i=1,ia - do k=1,nvmax + do k=1,size(jj) jj(k) = 0 enddo jj(i) = 1 @@ -3905,8 +5570,56 @@ subroutine dader(idif,ina,inc) ! !----------------------------------------------------------------------------- ! - integer idif,illc,ilmc,ina,inc,incc,inoc,invc,ipoc - ! + integer idif,illc,ilmc,ina,inc,incc,inoc,invc,ipoc,inoa,inva,ipoa,ilma,illa + integer i,jd(lnv),ic1,ic2 + real(dp) rr + + if(newtpsa) then + call dainf(ina,inoa,inva,ipoa,ilma,illa) + call dainf(inc,inoc,invc,ipoc,ilmc,illc) + if(nomax.eq.1) then + ! PRINT*,'ERROR, DADER CALLED WITH c_nomax = 1' + ! call dadeb !(31,'ERR DADER1',1) + ! stop 666 + do i=1,lnv + jd(i)=0 + enddo + jd(idif)=1 + call dapek(ina,jd,rr) + call dacon(inc,rr) + return + endif + if(nomax.eq.2) then + reel=0 + + + do i=1,combien + + ic1=ind1(i) + ic2=ind2(i) + + if(ic1==ic2) then + if(ic1==idif) then + ic1=0 + reel(inds(ic1,ic2))=2.0_dp*cc(ipoa+i-1) + endif + elseif(ic1==idif) then + ic1=0 + reel(inds(ic1,ic2))=cc(ipoa+i-1) + elseif(ic2==idif) then + ic2=0 + reel(inds(ic1,ic2))=cc(ipoa+i-1) + endif + + enddo + + cc(ipoc:ipoc+combien-1) = reel + + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -3931,6 +5644,7 @@ subroutine dader(idif,ina,inc) endif return end subroutine dader + subroutine dadert(idif,ina,inc) implicit none ! ****************************** @@ -3944,6 +5658,20 @@ subroutine dadert(idif,ina,inc) ilma,ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj,ipause,mypauses integer,dimension(lnv)::jd real(dp) rr,x,xdivi + + if(newtpsa) then + if(nomax==1) then + stop 310 + return + endif + + if(nomax==2) then + stop 311 + return + endif +endif + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -4019,318 +5747,105 @@ subroutine dadert(idif,ina,inc) cc(ic) = cc(i)*ifac i_1(ic) = i_1(i) - ider1s i_2(ic) = i_2(i) - ider2s - ! -100 continue - enddo - ! - idall(inc) = ic - ipoc + 1 - if(idall(inc).gt.idalm(inc)) then - write(line,'(a15)') 'ERROR IN DADER ' - ipause=mypauses(35,line) - call dadeb !(31,'ERR DADER2',1) - endif - ! - return - end subroutine dadert - ! - ! - subroutine dacfuR(ina,fun,inc) - implicit none - ! ***************************** - ! - ! THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION - ! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE - ! RESULT IN C - ! - !----------------------------------------------------------------------------- - ! - integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc - complex(dp),external::fun - ! - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - if(ina.eq.inc) then - call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - incc=0 - call daall1(incc,'$$DAJUNK$$',inoc,invc) - call dacfuRt(ina,fun,incc) - call dacop(incc,inc) - call dadal1(incc) - else - call dacfuRt(ina,fun,inc) - endif - return - end subroutine dacfuR - subroutine dacfuRt(ina,fun,inc) - implicit none - ! external fun - ! ***************************** - ! - ! THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION - ! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE - ! RESULT IN C - ! - !----------------------------------------------------------------------------- - ! - integer i,ia,ic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,& - invc,ipoa,ipoc,ipause,mypauses - integer,dimension(lnv)::j - real(dp) cfac,rr - ! - interface - ! complex(kind(one)) function fun(abc) - function fun(abc) - use precision_constants - implicit none - complex(dp) fun - integer,dimension(:)::abc - end function fun - end interface - ! - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - call dainf(ina,inoa,inva,ipoa,ilma,illa) - call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - ! - ! - if(nomax.eq.1) then - do i=1,lnv - j(i)=0 - enddo - call dapek(ina,j,rr) - cfac = REAL(fun(j),kind=DP) - rr=cfac*rr - call dapok(inc,j,rr) - do i=1,lnv - j(i)=1 - call dapek(ina,j,rr) - cfac = REAL(fun(j),kind=DP) - rr=cfac*rr - call dapok(inc,j,rr) - j(i)=0 - enddo - ! PRINT*,'ERROR, DACFU CALLED WITH NOMAX = 1' - ! call dadeb !(31,'ERR DACFU ',1) - ! stop 667 - return - endif - ! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INC,INOC,INVC) - ! - ic = ipoc - 1 - ! - do ia=ipoa,ipoa+illa-1 - ! - call dancd(i_1(ia),i_2(ia),j) - cfac = REAL(fun(j),kind=DP) - ! IF(abs(CFAC).LT.EPS) GOTO 100 - ! IF(abs(CFAC*CC(IA)).LT.EPS) GOTO 100 - if(abs(cfac*cc(ia)).lt.eps_da.or.abs(cc(ia)).lt.eps_da) goto 100 - ! - ic = ic + 1 - cc(ic) = cc(ia)*cfac - i_1(ic) = i_1(ia) - i_2(ic) = i_2(ia) - ! -100 continue - enddo - ! - idall(inc) = ic - ipoc + 1 - if(idall(inc).gt.idalm(inc)) then - write(line,'(a15)') 'ERROR IN DACFU ' - ipause=mypauses(36,line) - call dadeb !(31,'ERR DACFU ',1) - endif - ! - return - end subroutine dacfuRt - ! - subroutine dacfu(ina,fun,inc) - implicit none - ! ***************************** - ! - ! THIS SUBROUTINE APPLIES THE EXTERNAL real(dp) FUNCTION - ! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE - ! RESULT IN C - ! - !----------------------------------------------------------------------------- - ! - integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc - real(dp),external::fun - ! - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - if(ina.eq.inc) then - call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - incc=0 - call daall1(incc,'$$DAJUNK$$',inoc,invc) - call dacfut(ina,fun,incc) - call dacop(incc,inc) - call dadal1(incc) - else - call dacfut(ina,fun,inc) - endif - return - end subroutine dacfu - subroutine dacfuI(ina,fun,inc) - implicit none - ! ***************************** - ! - ! THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION - ! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE - ! RESULT IN C - ! - !----------------------------------------------------------------------------- - ! - integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc - complex(dp),external::fun - ! - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - if(ina.eq.inc) then - call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif - incc=0 - call daall1(incc,'$$DAJUNK$$',inoc,invc) - call dacfuIt(ina,fun,incc) - call dacop(incc,inc) - call dadal1(incc) - else - call dacfuIt(ina,fun,inc) + ! +100 continue + enddo + ! + idall(inc) = ic - ipoc + 1 + if(idall(inc).gt.idalm(inc)) then + write(line,'(a15)') 'ERROR IN DADER ' + ipause=mypauses(35,line) + call dadeb !(31,'ERR DADER2',1) endif + ! return - end subroutine dacfuI - subroutine dacfuIt(ina,fun,inc) + end subroutine dadert + ! + + subroutine dacfu(ina,fun,inc) implicit none - ! external fun ! ***************************** ! - ! THIS SUBROUTINE APPLIES THE EXTERNAL DOUBLE COMPLEX FUNCTION + ! THIS SUBROUTINE APPLIES THE EXTERNAL real(dp) FUNCTION ! OF THE EXPONENTS FUN TO EACH COEFFICIENT OF A AND STORES THE ! RESULT IN C ! !----------------------------------------------------------------------------- ! - integer i,ia,ic,illa,illc,ilma,ilmc,ina,inc,inoa,inoc,inva,& - invc,ipoa,ipoc,ipause,mypauses + integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc + integer inoa,inva,ipoa,ilma,illa integer,dimension(lnv)::j - real(dp) cfac,rr - ! + integer i interface - ! complex(kind(one)) function fun(abc) + ! real(kind(one)) function fun(abc) function fun(abc) use precision_constants implicit none - complex(dp) fun + real(dp) fun integer,dimension(:)::abc end function fun end interface ! - if((.not.C_%STABLE_DA)) then - if(c_%watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return - endif + + if(newtpsa) then call dainf(ina,inoa,inva,ipoa,ilma,illa) call dainf(inc,inoc,invc,ipoc,ilmc,illc) + if(nomax==1) then + do i=1,lnv + j(i)=0 + enddo + + do i=0,nvmax + if(i/=0) j(i)=1 + cc(ipoc+i)=fun(j)*cc(ipoa+i) + if(i/=0) j(i)=0 + enddo + return + endif + + if(nomax==2) then + do i=0,nvmax + if(i/=0) j(i)=1 + cc(ipoc+i)=fun(j)*cc(ipoa+i) + if(i/=0) j(i)=0 + enddo + do i=poscombien,combien + j=0 + j(ind1(i))=j(ind1(i))+1 + j(ind2(i))=j(ind2(i))+1 + cc(ipoc+i-1)=fun(j)*cc(ipoa+i-1) + j(ind1(i)) = 0 + j(ind2(i)) = 0 + enddo + return + endif +endif + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif return endif - ! - ! - if(nomax.eq.1) then - do i=1,lnv - j(i)=0 - enddo - call dapek(ina,j,rr) - cfac = AIMAG(fun(j)) - rr=cfac*rr - call dapok(inc,j,rr) - do i=1,lnv - j(i)=1 - call dapek(ina,j,rr) - cfac = AIMAG(fun(j)) - rr=cfac*rr - call dapok(inc,j,rr) - j(i)=0 - enddo - ! PRINT*,'ERROR, DACFU CALLED WITH NOMAX = 1' - ! call dadeb !(31,'ERR DACFU ',1) - ! stop 667 - return - endif - ! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INC,INOC,INVC) - ! - ic = ipoc - 1 - ! - do ia=ipoa,ipoa+illa-1 - ! - call dancd(i_1(ia),i_2(ia),j) - cfac = aimag(fun(j)) - ! IF(abs(CFAC).LT.EPS) GOTO 100 - ! IF(abs(CFAC*CC(IA)).LT.EPS) GOTO 100 - if(abs(cfac*cc(ia)).lt.eps_da.or.abs(cc(ia)).lt.eps_da) goto 100 - ! - ic = ic + 1 - cc(ic) = cc(ia)*cfac - i_1(ic) = i_1(ia) - i_2(ic) = i_2(ia) - ! -100 continue - enddo - ! - idall(inc) = ic - ipoc + 1 - if(idall(inc).gt.idalm(inc)) then - write(line,'(a15)') 'ERROR IN DACFU ' - ipause=mypauses(37,line) - call dadeb !(31,'ERR DACFU ',1) + if(ina.eq.inc) then + call dainf(inc,inoc,invc,ipoc,ilmc,illc) + if((.not.C_%STABLE_DA)) then + if(c_%watch_user) then + write(6,*) "big problem in dabnew ", sqrt(crash) + endif + return + endif + incc=0 + call daall1(incc,'$$DAJUNK$$',inoc,invc) + call dacfut(ina,fun,incc) + call dacop(incc,inc) + call dadal1(incc) + else + call dacfut(ina,fun,inc) endif - ! return - end subroutine dacfuIt - ! + end subroutine dacfu + subroutine dacfut(ina,fun,inc) implicit none ! external fun @@ -4356,6 +5871,11 @@ function fun(abc) integer,dimension(:)::abc end function fun end interface + + + + + ! if((.not.C_%STABLE_DA)) then if(c_%watch_user) then @@ -4377,21 +5897,13 @@ end function fun do i=1,lnv j(i)=0 enddo - call dapek(ina,j,rr) - cfac = fun(j) - rr=cfac*rr - call dapok(inc,j,rr) - do i=1,lnv - j(i)=1 - call dapek(ina,j,rr) - cfac = fun(j) - rr=cfac*rr - call dapok(inc,j,rr) - j(i)=0 - enddo - ! PRINT*,'ERROR, DACFU CALLED WITH NOMAX = 1' - ! call dadeb !(31,'ERR DACFU ',1) - ! stop 667 + + do i=0,nvmax + if(i/=0) j(i)=1 + cc(ipoc+i)=fun(j)*cc(ipoa+i) + if(i/=0) j(i)=0 + enddo + return endif ! CALL DACHK(INA,INOA,INVA, ' ',-1,-1,INC,INOC,INVC) @@ -4423,6 +5935,7 @@ end function fun ! return end subroutine dacfut + ! subroutine GET_C_J(ina,I,C,J) ! subroutine GET_C_J(ina,I,C,J) ! implicit none @@ -4446,6 +5959,26 @@ subroutine dapri(ina,iunit) integer i,ii,iii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,iunit,k integer,dimension(lnv)::j ! + +if(newtpsa) then + if(nomax.eq.1) then + write(iunit,*) "1st order polynomial ", ina,nvmax + do i=1,nvmax+1 + if(cc(idapo(ina)+i-1).ne.0.0_dp) write(iunit,*) i-1, cc(idapo(ina)+i-1) + enddo + + return + endif + + if(nomax==2) then + write(iunit,*) "2nd order polynomial ", ina,nvmax + do i=1,combien + if(cc(idapo(ina)+i-1).ne.0.0_dp) write(iunit,*) ind1(i),ind2(i),cc(idapo(ina)+i-1) + enddo + return + endif +endif + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -4514,6 +6047,7 @@ subroutine dapri(ina,iunit) ! return end subroutine dapri + subroutine dapri77(ina,iunit,ind) implicit none ! *************************** @@ -4527,6 +6061,57 @@ subroutine dapri77(ina,iunit,ind) integer,dimension(lnv)::j character(10) c10,k10 logical some + integer(2) xi(3),count + + some=.false. + +if(newtpsa) then + if(nomax.eq.1) then + count=0 + write(iunit,*) "1st order polynomial ", ina,nvmax + do i=1,nvmax+1 + if(cc(idapo(ina)+i-1).ne.0.0_dp) then + xi(1)=i-1 + xi(2)=ind1(i) + xi(3)=ind2(i) + if(xi(2)>nvmax-nphere) then + xi(2)=-(xi(2)-(nvmax-nphere)) + endif + if(xi(3)>nvmax-nphere) then + xi(3)=-(xi(3)-(nvmax-nphere)) + endif + write(iunit,*) xi, cc(idapo(ina)+i-1) + count=count+1 + endif + enddo + write(iunit,*) count, " monomial(s) printed " + return + endif + + if(nomax==2) then + count=0 + + write(iunit,*) "2nd order polynomial ", ina,nvmax + do i=1,combien + if(cc(idapo(ina)+i-1).ne.0.0_dp) then + xi(1)=i-1 + xi(2)=ind1(i) + xi(3)=ind2(i) + if(xi(2)>nvmax-nphere) then + xi(2)=-(xi(2)-(nvmax-nphere)) + endif + if(xi(3)>nvmax-nphere) then + xi(3)=-(xi(3)-(nvmax-nphere)) + endif + write(iunit,*) xi, cc(idapo(ina)+i-1) + count=count+1 + endif + enddo + write(iunit,*) count, " monomial(s) printed " + return + endif +endif + some=.false. ! if(iunit.eq.0) return @@ -4630,6 +6215,7 @@ subroutine dapri77(ina,iunit,ind) ! return end subroutine dapri77 + subroutine dashift(ina,inc,ishift) implicit none ! real(dp) c @@ -4638,9 +6224,77 @@ subroutine dashift(ina,inc,ishift) !----------------------------------------------------------------------------- ! integer i,ii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,inb,ishift,ich,& - ik,inc,k + ik,inc,k,ipob,ic1,ic2 integer,dimension(lnv)::j,jd ! + + if(newtpsa) then + inva = idanv(ina) + ipoa = idapo(ina) + ipob = idapo(ina) + + if(nomax==1) then + call daall1(inb,'$$DAJUNK$$',inoa,inva) + + cc(ipob:ipob+inva)=0 + + do ii=1,inva + if(cc(ipoa+ii)/=0.0_dp.and.ii<=ishift) then + write(6,*) " trouble in dashift " + stop 886 + else + cc(ipob+ii-ishift)=cc(ipoa+ii) + endif + enddo + cc(ipob)=cc(ipoa) + call dacop(inb,inc) + call dadal1(inb) + return + endif + + if(nomax.eq.2) then + call daall1(inb,'$$DAJUNK$$',inoa,inva) + ipob = idapo(inb) + + cc(ipob:ipob+combien-1)=0 + + do ii=1,inva !ipoa,ipoa+illa-1 + if(cc(ipoa+ii)/=0.0_dp.and.ii<=ishift) then + write(6,*) " trouble in dashift " + stop 887 + else + cc(ipob+ii-ishift)=cc(ipoa+ii) + endif + enddo + + do ii=poscombien,combien !ipoa,ipoa+illa-1 + + + if( ((ind1(ii)<=ishift).or.(ind2(ii)<=ishift)).and.cc(ipoa+ii-1)/=0.0_dp) then + write(6,*) " trouble in dashift " + stop 889 + else + if(((ind1(ii)>ishift).or.(ind1(ii)==0)).and.((ind2(ii)>ishift).or.(ind2(ii)==0))) then + ic1=0;ic2=0; + if(ind1(ii)>ishift) ic1 = ind1(ii)-ishift + if(ind2(ii)>ishift) ic2 = ind2(ii)-ishift + + cc(ipob+inds(ic1,ic2)-1)=cc(ipoa+ii-1) + + endif + endif + enddo + + cc(ipob)=cc(ipoa) + call dacop(inb,inc) + call dadal1(inb) + return + endif + + endif + + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -4741,6 +6395,7 @@ subroutine dashift(ina,inc,ishift) ! return end subroutine dashift + subroutine darea(ina,iunit) implicit none ! Frank @@ -4751,6 +6406,18 @@ subroutine darea(ina,iunit) integer,dimension(lnv)::j real(dp) c character(10) c10 + + + if(newtpsa) then + if(nomax==1) then + return + endif + + if(nomax==2) then + return + endif +endif + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -4859,6 +6526,18 @@ subroutine darea77(ina,iunit) real(dp) c character(10) c10,k10 ! + + if(newtpsa) then + if(nomax==1) then + return + endif + + if(nomax==2) then + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -4959,6 +6638,10 @@ subroutine dainf(inc,inoc,invc,ipoc,ilmc,illc) ! integer illc,ilmc,inc,inoc,invc,ipoc,ipause,mypauses ! + + + + if(inc.ge.1.and.inc.le.nda_dab) then inoc = idano(inc) invc = idanv(inc) @@ -5192,9 +6875,48 @@ subroutine datra(idif,ina,inc) !----------------------------------------------------------------------------- ! integer i,ibase,ic,ider1,ider1s,ider2,ider2s,idif,iee,ifac,illa,illc,ilma,& - ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj,ipause,mypauses + ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj,ipause,mypauses,ic1,ic2 real(dp) x,xdivi ! + + if(newtpsa) then + if(nomax==1) then + call dader(idif,ina,inc) + return + endif + + if(nomax==2) then + reel=0 + ipoa=idapo(ina) + ipoc=idapo(inc) + + do i=1,combien + + ic1=ind1(i) + ic2=ind2(i) + + if(ic1==ic2) then + if(ic1==idif) then + ic1=0 + reel(inds(ic1,ic2))= cc(ipoa+i-1) + endif + elseif(ic1==idif) then + ic1=0 + reel(inds(ic1,ic2))=cc(ipoa+i-1) + elseif(ic2==idif) then + ic2=0 + reel(inds(ic1,ic2))=cc(ipoa+i-1) + endif + + enddo + + cc(ipoc:ipoc+combien-1) = reel + + return + endif +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -5277,6 +6999,7 @@ subroutine datra(idif,ina,inc) ! return end subroutine datra + subroutine hash(no1,nv1,jj,ic1,ic2) implicit none ! **************************** @@ -5336,6 +7059,7 @@ subroutine dehash(no1,nv1,ic1,ic2,jj) ! return end subroutine dehash + subroutine daran(ina,cm,xran) implicit none ! ************************ @@ -5350,6 +7074,30 @@ subroutine daran(ina,cm,xran) integer i,illa,ilma,ina,inoa,inva,ipoa,ipause,mypauses real(dp) cm,xran ! + + if(newtpsa) then + +ipoa=idapo(ina) + + do i=ipoa,ipoa+nmmax-1 + if(cm.gt.zero) then + cc(i) = bran(xran) + if(cc(i).gt.cm) cc(i) = zero + elseif(cm.lt.zero) then + cc(i) = int(1+10*bran(xran)) + if(cc(i).gt.-ten*cm) cc(i) = zero + else + line='ERROR IN ROUTINE DARAN' + ipause=mypauses(31,line) + call dadeb !(31,'ERR DARAN2',1) + endif + enddo + + return + +endif + + if((.not.C_%STABLE_DA)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) @@ -5414,10 +7162,35 @@ subroutine dacycle(ina,ipresent,value,illa,j) ! !----------------------------------------------------------------------------- ! - integer ii,illa,ilma,ina,inoa,inva,iout,ipoa,ipresent + integer ii,illa,ilma,ina,inoa,inva,iout,ipoa,ipresent,i integer,optional,dimension(:)::j real(dp) value ! + + if(newtpsa) then + + ipoa=idapo(ina) + if(.not.present(j)) then + illa=nmmax + ! illa0$=0 + + ! do i=1,nmmax + ! if(abs(cc(ipoa+i-1))/=0) illa=illa+1 + ! enddo + else + + value=cc(ipoa+ipresent-1) + j=0 + + if(ind1(ipresent)/=0) j(ind1(ipresent))=1+j(ind1(ipresent)) + if(ind2(ipresent)/=0) j(ind2(ipresent))=1+j(ind2(ipresent)) + + + return + endif +endif + + if(ina.lt.1.or.ina.gt.nda_dab) then write(line,'(a22,i8)') 'ERROR IN dacycle, INA = ',ina ipause=mypauses(39,line) @@ -5446,35 +7219,7 @@ subroutine dacycle(ina,ipresent,value,illa,j) ! ipresent=1+ipresent return end subroutine dacycle -!!!! new stuff lingyun - subroutine daclean(ina,value) - implicit none - integer ipause, mypauses - ! *************************** - ! - ! THIS SUBROUTINE PRINTS THE DA VECTOR INA TO UNIT IUNIT. - ! - !----------------------------------------------------------------------------- - ! - integer ii,illa,ilma,ina,inoa,inva,iout,ipoa - real(dp) value - ! - if(ina.lt.1.or.ina.gt.nda_dab) then - write(line,'(a22,i8)') 'ERROR IN dacycle, INA = ',ina - ipause=mypauses(39,line) - stop - endif - ! - inoa = idano(ina) - inva = idanv(ina) - ipoa = idapo(ina) - ilma = idalm(ina) - illa = idall(ina) - iout = 0 - do ii=ipoa,ipoa+ilma-1 - if(abs(cc(ii))2.or.(.not.newtpsa)) then + allocate(c_i_1(c_lst));allocate(c_i_2(c_lst)); + allocate(c_ie1(c_lea));allocate(c_ie2(c_lea));allocate(c_ieo(c_lea)); + allocate(c_ia1(0:c_lia));allocate(c_ia2(0:c_lia)); +else + allocate(c_i_1(1));allocate(c_i_2(1)); + allocate(c_ie1(1));allocate(c_ie2(1));allocate(c_ieo(1)); + allocate(c_ia1(0:1));allocate(c_ia2(0:1)); + +endif +allocate(c_idano(c_lda));allocate(c_idanv(c_lda));allocate(c_idapo(c_lda)); +allocate(c_idalm(c_lda));allocate(c_idall(c_lda)); +allocate(c_daname(c_lda)); +allocate(c_allvec(c_lda)); + + ! do i=1,c_lst ! c_cc(I)=0.0_dp ! ADDED BY ETIENNE ! c_i_1(i)=0 diff --git a/forest/code/cc_dabnew.f90 b/forest/code/cc_dabnew.f90 index 0c9a69837c..841ea829ed 100644 --- a/forest/code/cc_dabnew.f90 +++ b/forest/code/cc_dabnew.f90 @@ -8,12 +8,13 @@ module c_dabnew implicit none private ! public - public c_ldamax, c_lstmax,c_leamax,c_liamax - private daalc_lno1,daall,damult,dasqrt,dacmut,dacma,DALINt,dacctt - private dainvt,dadert,dacfut + public c_ldamax, c_lstmax,c_leamax,c_liamax,c_dacctt2datest + + private daalc_lno1,daall,damult,dasqrt,dacmut,dacma,DALINt,dacctt,dacctt1,dacctt2tpsa,dacctt2da + private dainvt,dadert,dacfut,dainvt1 private dadeb,dapac,dachk,damch,dadcd,dancd,hash,dehash - public C_STABLE_DA,C_watch_user,c_daini,c_daall0,c_davar,c_damul,c_dacycle + public c_stable_da,C_watch_user,c_daini,c_daall0,c_davar,c_damul,c_dacycle public c_dapok,c_dapek,c_etall1,c_DAshift,c_dacop,c_dacon,c_daadd,c_dacad,c_dasub public c_dacsu,c_dasuc,c_dacmu,c_dadal1,c_dapri,c_dapri77,c_darea,c_darea77,c_daeps public c_dacdi,c_dadic,c_count_da,c_mtree,c_dafun,c_DAABS,c_dadiv,c_take,c_datrunc,c_dader,c_datra @@ -24,11 +25,43 @@ module c_dabnew ! integer,private,parameter::nmax=400,lsw=1 ! real(dp),private,parameter::tiny=c_1d_20 character(120),private :: line - logical(lp) C_STABLE_DA,C_watch_user,C_check_stable + logical(lp) c_stable_da,C_watch_user,C_check_stable real(dp), private :: eps=1.d-38,epsprint=1.d-38 !real(dp),public :: eps_clean=0 public c_print_c_nda_dab_c_lda complex(dp), private :: i_=(0.0_dp,1.0_dp) + + + + INTERFACE c_daall1$ + MODULE PROCEDURE daall1 + end INTERFACE c_daall1$ + + + INTERFACE c_daall$ + MODULE PROCEDURE daall + end INTERFACE c_daall$ + + INTERFACE c_daall0$ + MODULE PROCEDURE c_daall0 + end INTERFACE c_daall0$ + + INTERFACE c_dadal$ + MODULE PROCEDURE dadal + end INTERFACE c_dadal$ + + INTERFACE c_dadal1$ + MODULE PROCEDURE c_dadal1 + end INTERFACE c_dadal1$ + + INTERFACE c_davar$ + MODULE PROCEDURE c_davar + end INTERFACE c_davar$ + + INTERFACE c_dacma$ + MODULE PROCEDURE dacma + end INTERFACE c_dacma$ + contains !****************************************************************************** ! * @@ -227,16 +260,80 @@ subroutine c_daini(no,nv,iunit) integer,dimension(c_lnv)::j,jj character(10) aa ! + c_nomax = no + c_nvmax = nv + c_last_tpsa=1 +if(newtpsa) then + + if(c_nomax.eq.1) then + call alloc_all_c(no,nv) + c_ndamaxi=0 + do i=1, c_lda + c_allvec(i) = .false. + enddo + c_nhole=0 + + c_nda_dab = 0 + c_nst0 = 0 + c_nomax = no + c_nvmax = nv + c_nocut=no + call danum_c(no,nv,c_nmmax) + + do i=1,c_lda + c_idapo(i)=c_nmmax*(i-1)+1 + enddo + c_idano=1 + c_idanv=nv + c_idalm=c_nmmax + c_idall=c_nmmax + return + endif + + if(c_nomax==2) then + + call alloc_all_c(no,nv) + + c_ndamaxi=0 + do i=1, c_lda + c_allvec(i) = .false. + enddo + c_nhole=0 + + c_nda_dab = 0 + c_nst0 = 0 + c_nomax = no + c_nvmax = nv + c_nocut=no + + + call danum_c(no,nv,c_nmmax) + + + do i=1,c_lda + c_idapo(i)=c_nmmax*(i-1)+1 + enddo + c_idano=2 + c_idanv=nv + c_idalm=c_nmmax + c_idall=c_nmmax + + return + endif + +endif + + !frs if(eps.le.zero) eps=c_1d_38 ! if(EPS.le.zero) eps=c_1d_90 !frs epsmac=c_1d_7 - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif return endif - c_last_tpsa=1 + ! c_last_tpsa=1 if(nv.eq.0) return call alloc_all_c(no,nv) c_ndamaxi=0 @@ -452,6 +549,7 @@ subroutine c_daini(no,nv,iunit) enddo ! return + end subroutine c_daini ! subroutine daexter ! implicit none @@ -553,13 +651,80 @@ subroutine daalc_lno1(ic,ccc) logical(lp) incnda integer ind,ndanum,no,nv,ic,ipause,mypauses character(10) c,ccc - ! if((.not.C_STABLE_DA)) then + ! if((.not.c_stable_da)) then ! if(C_watch_user) then ! write(6,*) "big problem in dabnew ", sqrt(crash) ! endif ! return ! endif ! + +if(newtpsa) then + if(c_nomax.eq.1) then + + if(c_nhole.gt.0) then + ind=c_nda_dab +200 if (c_allvec(ind)) then + ind = ind - 1 + goto 200 + endif + incnda = .false. + c_nhole=c_nhole-1 + else + incnda = .true. + c_nda_dab = c_nda_dab + 1 + ind=c_nda_dab + if(c_nda_dab.gt.c_lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 200 + endif + endif + if(ind>c_lda_max_used) c_lda_max_used=ind + if(ind>c_lda) then + write(6,*) "ind>lda ",c_lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',c_LDA + stop 201 + endif + c_allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + c_daname(ind) = ccc + ic = ind + return + endif + + if(c_nomax==2) then + + if(c_nhole.gt.0) then + ind=c_nda_dab +300 if (c_allvec(ind)) then + ind = ind - 1 + goto 300 + endif + incnda = .false. + c_nhole=c_nhole-1 + else + incnda = .true. + c_nda_dab = c_nda_dab + 1 + ind=c_nda_dab + if(c_nda_dab.gt.c_lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 3001 + endif + endif + if(ind>c_lda_max_used) c_lda_max_used=ind + if(ind>c_lda) then + write(6,*) "ind>lda ",c_lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',c_LDA + stop 301 + endif + c_allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + c_daname(ind) = ccc + ic = ind + return + endif +endif + no=c_nomax nv=c_nvmax ind = 1 @@ -654,7 +819,78 @@ subroutine daall(ic,l,ccc,no,nv) integer i,ind,l,ndanum,no,nv,ipause,mypauses integer,dimension(:)::ic character(10) c,ccc - ! if((.not.C_STABLE_DA)) then + + +if(newtpsa) then + if(c_nomax.eq.1) then +do i=1,l + if(c_nhole.gt.0) then + ind=c_nda_dab +200 if (c_allvec(ind)) then + ind = ind - 1 + goto 200 + endif + incnda = .false. + c_nhole=c_nhole-1 + else + incnda = .true. + c_nda_dab = c_nda_dab + 1 + ind=c_nda_dab + if(c_nda_dab.gt.c_lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 200 + endif + endif + if(ind>c_lda_max_used) c_lda_max_used=ind + if(ind>c_lda) then + write(6,*) "ind>c_lda ",c_lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: c_lda = ',c_lda + stop 201 + endif + c_allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + c_daname(ind) = ccc + ic(i) = ind + enddo + return + endif + + + if(c_nomax==2) then + do i=1,l + if(c_nhole.gt.0) then + ind=c_nda_dab +300 if (c_allvec(ind)) then + ind = ind - 1 + goto 300 + endif + incnda = .false. + c_nhole=c_nhole-1 + else + incnda = .true. + c_nda_dab = c_nda_dab + 1 + ind=c_nda_dab + if(c_nda_dab.gt.c_lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 3002 + endif + endif + if(ind>c_lda_max_used) c_lda_max_used=ind + if(ind>c_lda) then + write(6,*) "ind>c_lda ",c_lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: c_lda = ',c_lda + stop 301 + endif + c_allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + c_daname(ind) = ccc + ic(i) = ind + enddo + return + endif +endif + + ! if((.not.c_stable_da)) then ! if(C_watch_user) then ! write(6,*) "big problem in dabnew ", sqrt(crash) ! endif @@ -745,6 +981,7 @@ subroutine daall(ic,l,ccc,no,nv) if(c_nda_dab.gt.c_ndamaxi) c_ndamaxi=c_nda_dab return end subroutine daall + subroutine daall1(ic,ccc,no,nv) implicit none ! ******************************** @@ -757,7 +994,75 @@ subroutine daall1(ic,ccc,no,nv) logical(lp) incnda integer ic,ind,ndanum,no,nv,ipause,mypauses character(10) c,ccc - ! if((.not.C_STABLE_DA)) then + + +if(newtpsa) then + if(c_nomax.eq.1) then + + if(c_nhole.gt.0) then + ind=c_nda_dab +200 if (c_allvec(ind)) then + ind = ind - 1 + goto 200 + endif + incnda = .false. + c_nhole=c_nhole-1 + else + incnda = .true. + c_nda_dab = c_nda_dab + 1 + ind=c_nda_dab + if(c_nda_dab.gt.c_lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 200 + endif + endif + if(ind>c_lda_max_used) c_lda_max_used=ind + if(ind>c_lda) then + write(6,*) "ind>lda ",c_lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',c_LDA + stop 201 + endif + c_allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + c_daname(ind) = ccc + ic = ind + return + endif + + if(c_nomax==2) then + + if(c_nhole.gt.0) then + ind=c_nda_dab +300 if (c_allvec(ind)) then + ind = ind - 1 + goto 300 + endif + incnda = .false. + c_nhole=c_nhole-1 + else + incnda = .true. + c_nda_dab = c_nda_dab + 1 + ind=c_nda_dab + if(c_nda_dab.gt.c_lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 3003 + endif + endif + if(ind>c_lda_max_used) c_lda_max_used=ind + if(ind>c_lda) then + write(6,*) "ind>lda ",c_lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',c_LDA + stop 301 + endif + c_allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + c_daname(ind) = ccc + ic = ind + return + endif +endif + + ! if((.not.c_stable_da)) then ! if(C_watch_user) then ! write(6,*) "big problem in dabnew ", sqrt(crash) ! endif @@ -871,10 +1176,77 @@ subroutine c_daall0(ic) logical(lp) incnda integer ic,ind,ndanum,no,nv,ipause,mypauses character(10) c,ccc + +if(newtpsa) then + if(c_nomax.eq.1) then + + if(c_nhole.gt.0) then + ind=c_nda_dab +200 if (c_allvec(ind)) then + ind = ind - 1 + goto 200 + endif + incnda = .false. + c_nhole=c_nhole-1 + else + incnda = .true. + c_nda_dab = c_nda_dab + 1 + ind=c_nda_dab + if(c_nda_dab.gt.c_lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 200 + endif + endif + if(ind>c_lda_max_used) c_lda_max_used=ind + if(ind>c_lda) then + write(6,*) "ind>lda ",c_lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',c_LDA + stop 201 + endif + c_allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + c_daname(ind) = ccc + ic = ind + return + endif + + if(c_nomax==2) then + + if(c_nhole.gt.0) then + ind=c_nda_dab +300 if (c_allvec(ind)) then + ind = ind - 1 + goto 300 + endif + incnda = .false. + c_nhole=c_nhole-1 + else + incnda = .true. + c_nda_dab = c_nda_dab + 1 + ind=c_nda_dab + if(c_nda_dab.gt.c_lda) then + write(line,'(a52)') 'ERROR IN DAALL, MAX NUMBER OF DA VECTORS EXHAUSTED,1' + stop 3004 + endif + endif + if(ind>c_lda_max_used) c_lda_max_used=ind + if(ind>c_lda) then + write(6,*) "ind>lda ",c_lda,ind + print*, 'ERROR IN DAALLNO1, MAX NUMBER OF DA VECTORS EXHAUSTED: LDA = ',c_LDA + stop 301 + endif + c_allvec(ind) = .true. + !idapo(ind)= ind*(nmmax-1)+1 + c_daname(ind) = ccc + ic = ind + return + endif +endif + ccc=' ' no=c_nomax nv=c_nvmax - ! if((.not.C_STABLE_DA)) then + ! if((.not.c_stable_da)) then ! if(C_watch_user) then ! write(6,*) "big problem in dabnew ", sqrt(crash) ! endif @@ -972,8 +1344,41 @@ subroutine dadal(idal,l) ! integer i,l,ipause,mypauses integer,dimension(:)::idal + +if(newtpsa) then + if(c_nomax.eq.1) then + do i=l,1,-1 + if(idal(i).eq.c_nda_dab) then + ! deallocate + c_nda_dab = c_nda_dab - 1 + else + c_nhole=c_nhole+1 + endif + c_cc(c_idapo(idal(i)):c_idapo(idal(i))+c_nmmax-1)=0 + c_allvec(idal(i)) = .false. + idal(i) = 0 + enddo + return + endif + + if(c_nomax==2) then + do i=l,1,-1 + if(idal(i).eq.c_nda_dab) then + ! deallocate + c_nda_dab = c_nda_dab - 1 + else + c_nhole=c_nhole+1 + endif + c_cc(c_idapo(idal(i)):c_idapo(idal(i))+c_nmmax-1)=0 + c_allvec(idal(i)) = .false. + idal(i) = 0 + enddo + return + endif +endif + ! - ! if((.not.C_STABLE_DA)) then + ! if((.not.c_stable_da)) then ! if(C_watch_user) then ! write(6,*) "big problem in dabnew ", sqrt(crash) ! endif @@ -983,7 +1388,7 @@ subroutine dadal(idal,l) if(idal(i).le.c_nomax+2.or.idal(i).gt.c_nda_dab) then write(line,'(a38,i8,1x,i8)') 'ERROR IN ROUTINE DADAL, IDAL(I),NDA = ',idal(i),c_nda_dab !ipause=mypauses(13,line) - C_%STABLE_DA = .false. + c_stable_da = .false. l = 1 return call dadeb !(31,'ERR DADAL ',1) @@ -1005,7 +1410,7 @@ subroutine dadal(idal,l) enddo return end subroutine dadal - + subroutine c_dadal1(idal) implicit none ! ************************ @@ -1015,13 +1420,42 @@ subroutine c_dadal1(idal) !----------------------------------------------------------------------------- ! integer idal - ! - ! if((.not.C_STABLE_DA)) then - ! if(C_watch_user) then - ! write(6,*) "big problem in dabnew ", sqrt(crash) - ! endif - ! return - ! endif + +if(newtpsa) then + if(c_nomax.eq.1) then + if(idal.eq.c_nda_dab) then + ! deallocate + c_nda_dab = c_nda_dab - 1 + else + c_nhole=c_nhole+1 + endif + c_cc(c_idapo(idal):c_idapo(idal)+c_nmmax-1)=0 + c_allvec(idal) = .false. + idal= 0 + return + endif + + if(c_nomax==2) then + if(idal.eq.c_nda_dab) then + ! deallocate + c_nda_dab = c_nda_dab - 1 + else + c_nhole=c_nhole+1 + endif + c_cc(c_idapo(idal):c_idapo(idal)+c_nmmax-1)=0 + c_allvec(idal) = .false. + idal= 0 + return + endif +endif + + ! + ! if((.not.c_stable_da)) then + ! if(C_watch_user) then + ! write(6,*) "big problem in dabnew ", sqrt(crash) + ! endif + ! return + ! endif if(idal.le.c_nomax+2.or.idal.gt.c_nda_dab) then write(6,'(a35,i8,1x,i8)') 'ERROR IN ROUTINE DADAL, IDAL,NDA = ',idal,c_nda_dab call dadeb !(31,'ERR DADAL ',1) @@ -1072,14 +1506,34 @@ subroutine c_davar(ina,ckon,i) integer i,ibase,ic1,ic2,illa,ilma,ina,inoa,inva,ipoa,ipause,mypauses complex(dp) ckon ! - if((.not.C_STABLE_DA)) then + + + if(newtpsa) then + if(c_nomax.eq.1) then + c_cc(c_idapo(ina):c_idapo(ina)+c_nmmax-1) =0 + c_cc(c_idapo(ina)) = ckon + c_cc(c_idapo(ina)+i) = one + return + endif + + if(c_nomax==2) then + + c_cc(c_idapo(ina):c_idapo(ina)+c_nmmax-1) =0 + c_cc(c_idapo(ina)) = ckon + c_cc(c_idapo(ina)+i) = one + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif return endif call dainf(ina,inoa,inva,ipoa,ilma,illa) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1143,14 +1597,30 @@ subroutine c_dacon(ina,ckon) integer illa,ilma,ina,inoa,inva,ipoa complex(dp) ckon ! - if((.not.C_STABLE_DA)) then + + + if(newtpsa) then + if(c_nomax==1) then + call daclr(ina) + c_cc(c_idapo(ina)) = ckon + return + endif + if(c_nomax==2) then + call daclr(ina) + c_cc(c_idapo(ina)) = ckon + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif return endif call dainf(ina,inoa,inva,ipoa,ilma,illa) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1182,7 +1652,7 @@ subroutine c_danot(not) ! integer not,ipause,mypauses ! - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1273,8 +1743,92 @@ subroutine c_dapek(ina,jv,cjj) integer,dimension(c_lnv)::jj complex(dp) cjj ! + + icu=0 + do i=1,c_nvmax + icu=jv(i)+icu + enddo + ic1=0 + do i=c_nvmax+1,size(jv) + ic1=jv(i)+ic1 + enddo +if(icu>c_nomax.or.ic1>0) then + write(6,*) "in cc_dabnew Crash 1, 0 continue ",icu,ic1 + read(5,*) i + if(i==1) then + cjj=sqrt(-cjj) + stop + endif + +! return +endif + + if(newtpsa) then +ipoa=c_idapo(ina) + + if(c_nomax==1) then + + icu = 0 +jj1=0 + jj=0 + do i=1,size(jv) + jj(i)=jv(i) + enddo + do i=1,c_nvmax + icu = icu + jj(i) + if(jj(i).eq.1) jj1=i + enddo + + if(icu>1) then + cjj=0 + return + endif + + ipek = ipoa + jj1 + cjj = c_cc(ipek) + return + endif + + if(c_nomax==2) then + jj=0 + do i=1,size(jv) + jj(i)=jv(i) + enddo + ic1=0 + cjj=0 + do i=1,c_nvmax + ic1=ic1+jj(i) + enddo + if(ic1>c_nomax) return + ic2=0 + ic1=0 + cjj=0 + do i=1,c_nvmax + if(jj(i)==2) then + ic1=i + ic2=i + exit + endif + + if(jj(i)==1) then + if(ic1==0) then + ic1=i + else + ic2=i + exit + endif + endif + enddo + + ipek = ipoa + inds(ic1,ic2) - 1 + cjj = c_cc(ipek) + return + endif +endif + + cjj=0 - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1285,7 +1839,7 @@ subroutine c_dapek(ina,jv,cjj) jj(i)=jv(i) enddo call dainf(ina,inoa,inva,ipoa,ilma,illa) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1405,8 +1959,90 @@ subroutine c_dapok(ina,jv,cjj) integer,dimension(:)::jv ! 2002.12.4 integer,dimension(c_lnv)::jj complex(dp) cjj + +icu=0 + do i=1,size(jv) + icu=jv(i)+icu + enddo + ic1=0 + do i=c_nvmax+1,size(jv) + ic1=jv(i)+ic1 + enddo +if(icu>c_nomax.or.ic1>0) then + write(6,*) "in cc_dabnew Crash 1, 0 continue ",icu,ic1 + read(5,*) i + if(i==1) then + cjj=sqrt(-cjj) + stop + endif + +! return +endif + if(newtpsa) then +ipoa=c_idapo(ina) + + if(c_nomax==1) then + icu = 0 + jj1=0 + jj=0 + do i=1,size(jv) + jj(i)=jv(i) + enddo + do i=1,c_nvmax + icu = icu + jj(i) + if(jj(i).eq.1) jj1=i + enddo + + if(icu>1) then + return + endif + + ipok = ipoa + jj1 + c_cc(ipok) = cjj + + return + endif + + if(c_nomax==2) then + jj=0 + do i=1,size(jv) + jj(i)=jv(i) + enddo + ic1=0 + do i=1,c_nvmax + ic1=ic1+jj(i) + enddo + if(ic1>c_nomax) return + ic2=0 + ic1=0 + do i=1,c_nvmax + if(jj(i)==2) then + ic1=i + ic2=i + exit + endif + + if(jj(i)==1) then + if(ic1==0) then + ic1=i + else + ic2=i + exit + endif + endif + enddo + + ipok = ipoa + inds(ic1,ic2) - 1 + c_cc(ipok)=cjj + + return + endif +endif + + + ! - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1418,7 +2054,7 @@ subroutine c_dapok(ina,jv,cjj) enddo ! call dainf(ina,inoa,inva,ipoa,ilma,illa) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1566,14 +2202,19 @@ subroutine daclr(inc) ! integer i,illc,ilmc,inc,inoc,invc,ipoc ! - if((.not.C_STABLE_DA)) then + if(newtpsa) then + c_cc(c_idapo(inc):c_idapo(inc)+c_nmmax-1)=0 + return +endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif return endif call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1604,8 +2245,14 @@ subroutine c_dacop(ina,inb) !----------------------------------------------------------------------------- ! integer ia,ib,illa,ina,inb,ipoa,ipob - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + c_cc(c_idapo(inb):c_idapo(inb)+c_nmmax-1)=c_cc(c_idapo(ina):c_idapo(ina)+c_nmmax-1) + return +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1634,6 +2281,7 @@ subroutine c_dacop(ina,inb) c_idall(inb) = ib - ipob + 1 return end subroutine c_dacop + subroutine c_daadd(ina,inb,inc) implicit none ! ***************************** @@ -1646,8 +2294,16 @@ subroutine c_daadd(ina,inb,inc) integer i,ina,ipoa integer idaadd,inb,inc,ipoc integer ipob - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + ipoc = c_idapo(inc) + ipoa = c_idapo(ina) + ipob = c_idapo(inb) + c_cc(ipoc:ipoc+c_nmmax-1) = c_cc(ipoa:ipoa+c_nmmax-1) + c_cc(ipob:ipob+c_nmmax-1) + return + +endif + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1679,7 +2335,40 @@ end subroutine c_daadd subroutine c_datrunc(ina,io,inb) implicit none integer ina,io,inb,nt,ipoca,ipocb,i - if((.not.C_%STABLE_DA)) then + + if(newtpsa) then + ipoca=c_idapo(ina) + ipocb=c_idapo(inb) + if(c_nomax==1) then + + do i=1,c_nvmax + c_cc(ipocb+i) =0.0_dp + enddo + + if(io==1) then + c_cc(ipocb) =c_cc(ipoca) + elseif(io>1) then + c_cc(ipocb:ipocb+c_nmmax-1)=c_cc(ipoca:ipoca+c_nmmax-1) + endif + return + endif + + if(c_nomax==2) then + c_cc(ipocb+1:ipocb+combien-1) = 0.0_dp + if(io==1) then + c_cc(ipocb) =c_cc(ipoca) + elseif(io==2) then + c_cc(ipocb:ipocb+c_nvmax) =c_cc(ipoca:ipoca+c_nvmax) + c_cc(ipocb+c_nvmax+1:ipocb+combien-1)=0.0_dp + elseif(io>2) then + c_cc(ipocb:ipocb+c_nmmax-1)=c_cc(ipoca:ipoca+c_nmmax-1) + endif + return + endif +endif + + + if((.not.c_stable_da)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1716,8 +2405,17 @@ subroutine c_dasub(ina,inb,inc) integer i,ina,ipoa integer inc,ipoc,inb integer ipob - ! - if((.not.C_STABLE_DA)) then + + + if(newtpsa) then + ipoc = c_idapo(inc) + ipoa = c_idapo(ina) + ipob = c_idapo(inb) + c_cc(ipoc:ipoc+c_nmmax-1) = c_cc(ipoa:ipoa+c_nmmax-1) - c_cc(ipob:ipob+c_nmmax-1) + return +endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1757,10 +2455,47 @@ subroutine c_damul(ina,inb,inc) ! !----------------------------------------------------------------------------- ! - integer ina,inb,inc,incc,ipoc,ipoa,ipob,i + integer ina,inb,inc,incc,ipoc,ipoa,ipob,i,j complex(dp) ccipoa,ccipob - ! - if((.not.C_STABLE_DA)) then + +if(newtpsa) then + ipoa=c_idapo(ina) + ipob=c_idapo(inb) + ipoc=c_idapo(inc) + if(c_nomax==1) then + ccipoa = c_cc(ipoa) + ccipob = c_cc(ipob) + c_cc(ipoc) = ccipoa*ccipob + do i=1,c_nvmax + c_cc(ipoc+i) = ccipoa*c_cc(ipob+i) + ccipob*c_cc(ipoa+i) + enddo + + return + endif + + if(c_nomax==2) then + ccipoa = c_cc(ipoa) + ccipob = c_cc(ipob) + reelc=0 + reelc(1) = ccipoa*ccipob + do i=1,c_nvmax + reelc(i+1) = ccipoa*c_cc(ipob+i) + ccipob*c_cc(ipoa+i)+reelc(i+1) + enddo + do i=poscombien,combien + reelc(i) = ccipoa*c_cc(ipob+i-1) + ccipob*c_cc(ipoa+i-1) +reelc(i) + enddo + do i=1,c_nvmax + do j=1,c_nvmax + reelc(inds(i,j))=reelc(inds(i,j))+c_cc(ipoa+i)*c_cc(ipob+j) + enddo + enddo + c_cc(ipoc:ipoc+combien-1)=reelc + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1817,8 +2552,20 @@ subroutine damult(ina,inb,inc) inva,invb,invc,ioffb,ipoa,ipob,ipoc,ipos,noib,nom integer,dimension(0:c_lno)::ipno,noff complex(dp) ccia,ccipoa,ccipob - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + if(c_nomax.eq.1) then + stop 1000 + return + endif + + if(c_nomax==2) then + stop 1001 + return + endif + endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1842,7 +2589,7 @@ subroutine damult(ina,inb,inc) call dainf(ina,inoa,inva,ipoa,ilma,illa) call dainf(inb,inob,invb,ipob,ilmb,illb) call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1928,8 +2675,35 @@ subroutine c_dadiv(ina,inb,inc) ! integer idadiv,inb,ina,inc,ipoc,ipoa,ipob,i complex(dp) ck,ck1 - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + ipoa = c_idapo(ina) + ipob = c_idapo(inb) + ipoc = c_idapo(inc) + + if(c_nomax==1) then + ck=1.0_dp/c_cc(ipob) + ck1=c_cc(ipoa)*ck + do i=1,c_nvmax + c_cc(ipoc+i) = (c_cc(ipoa+i)-c_cc(ipob+i)*ck1)*ck + enddo + c_cc(ipoc)=ck1 + return + endif + + if(c_nomax==2) then + idadiv = 0 + ! call dainf(inc,inoc,invc,ipoc,ilmc,illc) + call daall1(idadiv,'$$DADIV $$',c_nomax,c_nvmax) + call c_dafun('INV ',inb,idadiv) + call c_damul(ina,idadiv,inc) + call c_dadal1(idadiv) + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -1968,10 +2742,30 @@ subroutine dasqr(ina,inc) ! integer ina,inc,incc,ipoc,i,ipoa complex(dp) ccipoa - ! + + + if(newtpsa) then + if(c_nomax==1) then + ipoc = c_idapo(inc) + ipoa = c_idapo(ina) + ! minv = min(inva,invc) + ccipoa = c_cc(ipoa) + c_cc(ipoc) = ccipoa*ccipoa + do i=1,c_nvmax + c_cc(ipoc+i) = 2.0_dp*ccipoa*c_cc(ipoa+i) + enddo + return + endif + + if(c_nomax==2) then + call c_damul(ina,ina,inc) + return + endif +endif + ! CASE OF FIRST ORDER ONLY ! ************************ - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2000,6 +2794,7 @@ subroutine dasqr(ina,inc) endif return end subroutine dasqr + subroutine dasqrt(ina,inc) implicit none ! ************************* @@ -2021,8 +2816,20 @@ subroutine dasqrt(ina,inc) inva,invc,ioffa,ioffb,ipoa,ipoc,ipos,noia,noib,nom integer,dimension(0:c_lno)::ipno,noff complex(dp) ccia,ccipoa - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + if(c_nomax==1) then +stop 3501 + return + endif + + if(c_nomax==2) then +stop 3511 + return + endif +endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2041,7 +2848,7 @@ subroutine dasqrt(ina,inc) endif call dainf(ina,inoa,inva,ipoa,ilma,illa) call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2134,8 +2941,17 @@ subroutine c_dacad(ina,ckon,inb) integer ina,inb integer,parameter,dimension(c_lnv)::jjx=0 complex(dp) ckon,const - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + call c_dacop(ina,inb) + + c_cc(c_idapo(inb)) = c_cc(c_idapo(inb)) + ckon + + return + +endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2164,13 +2980,22 @@ subroutine c_dacsu(ina,ckon,inb) integer ina,inb integer,parameter,dimension(c_lnv)::jjx=0 complex(dp) ckon,const - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + call c_dacop(ina,inb) + + c_cc(c_idapo(inb)) = c_cc(c_idapo(inb)) - ckon + + return +endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif return endif + call c_dacop(ina,inb) ! if(c_nomax.eq.1) then @@ -2198,8 +3023,18 @@ subroutine c_dasuc(ina,ckon,inb) ! integer i,ina,inb,ipoa,ipob complex(dp) ckon - ! - if((.not.C_STABLE_DA)) then + + + if(newtpsa) then + ipob=c_idapo(inb) + ipoa=c_idapo(ina) + c_cc(ipob) = ckon - c_cc(ipoa) + c_cc(ipob+1:ipob+c_nmmax-1) =-c_cc(ipoa+1:ipoa+c_nmmax-1) + return + endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2232,9 +3067,17 @@ subroutine c_dacmu(ina,ckon,inc) ! integer ipoa,i,ina,inc,incc,ipoc complex(dp) ckon - ! - if((.not.C_STABLE_DA)) then - if(C_watch_user) then + + if(newtpsa) then + ipoc=c_idapo(inc) + ipoa=c_idapo(ina) + c_cc(ipoc:ipoc+c_nmmax-1) = ckon*c_cc(ipoa:ipoa+c_nmmax-1) + return + endif + + + if((.not.c_stable_da)) then + if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif return @@ -2279,8 +3122,22 @@ subroutine dacmut(ina,ckon,inb) integer i,ia,ib,illa,illb,ilma,ilmb,ina,inb,inoa,inob,inva,invb,ipoa,& ipob,ipause,mypauses complex(dp) ckon - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + if(c_nomax==1) then + stop 9001 + return + endif + + if(c_nomax==2) then + stop 9011 + + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2297,7 +3154,7 @@ subroutine dacmut(ina,ckon,inb) endif call dainf(ina,inoa,inva,ipoa,ilma,illa) call dainf(inb,inob,invb,ipob,ilmb,illb) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2341,8 +3198,16 @@ subroutine c_dacdi(ina,ckon,inb) ! integer i,ina,inb,ipoa,ipob,ipause,mypauses complex(dp) ckon - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + ipob=c_idapo(inb) + ipoa=c_idapo(ina) + c_cc(ipob:ipob+c_nmmax-1) = c_cc(ipoa:ipoa+c_nmmax-1)/ckon + return + endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2350,7 +3215,7 @@ subroutine c_dacdi(ina,ckon,inb) endif if(abs(ckon)==0.0_dp) then if(check_da) then - C_STABLE_DA=.false. + c_stable_da=.false. messagelost='constant part 0.0_dp in dacdi' return else @@ -2384,8 +3249,31 @@ subroutine c_dadic(ina,ckon,inc) ! integer i,idadic,ina,inc,ipoa,ipoc complex(dp) ckon,ck - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + if(c_nomax==1) then + ipoa = c_idapo(ina) + ipoc = c_idapo(inc) + c_cc(ipoc)=ckon/c_cc(ipoa) + ck=c_cc(ipoc)/c_cc(ipoa) + do i=1,c_nvmax + c_cc(ipoc+i) = -c_cc(ipoa+i)* ck + enddo + return + endif + + if(c_nomax==2) then + idadic = 0 + call daall1(idadic,'$$DADIC $$',c_nomax,c_nvmax) + call c_dafun('INV ',ina,idadic) + call c_dacmu(idadic,ckon,inc) + call c_dadal1(idadic) + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2393,7 +3281,7 @@ subroutine c_dadic(ina,ckon,inc) endif ipoa = c_idapo(ina) if(abs(c_cc(ipoa))==0.0_dp) then - if(check_da) C_STABLE_DA=.false. + if(check_da) c_stable_da=.false. messagelost='constant part 0.0_dp in dadic' return endif @@ -2439,8 +3327,19 @@ subroutine dacma(ina,inb,bfac,inc) ! integer idacma,ina,inb,inc,ipoc,ipob,ipoa,i complex(dp) bfac - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + ipoc = c_idapo(inc) + ipoa = c_idapo(ina) + ipob = c_idapo(inb) + + c_cc(ipoc:ipoc+c_nmmax-1) = c_cc(ipoa:ipoa+c_nmmax-1) + c_cc(ipob:ipob+c_nmmax-1) * bfac + return + +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2480,8 +3379,20 @@ subroutine dalin(ina,afac,inb,bfac,inc) !----------------------------------------------------------------------------- ! integer ipob,ipoa,i - ! - if((.not.C_STABLE_DA)) then + + + if(newtpsa) then + ipoc = c_idapo(inc) + ipoa = c_idapo(ina) + ipob = c_idapo(inb) + + c_cc(ipoc:ipoc+c_nmmax-1) = c_cc(ipoa:ipoa+c_nmmax-1) * afac + c_cc(ipob:ipob+c_nmmax-1) * bfac + + return + endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2529,8 +3440,21 @@ subroutine dalint(ina,afac,inb,bfac,inc) inob,inoc,inva,invb,invc,ipoa,ipob,ipoc,is,ismax,ismin,ja,jb,mchk,& ipause,mypauses complex(dp) afac,bfac,ccc,copf - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + if(c_nomax==1) then + stop 4001 + return + endif + + if(c_nomax==2) then + stop 4011 + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2551,7 +3475,7 @@ subroutine dalint(ina,afac,inb,bfac,inc) call dainf(ina,inoa,inva,ipoa,ilma,illa) call dainf(inb,inob,invb,ipob,ilmb,illb) call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2701,15 +3625,28 @@ subroutine c_daabs(ina,anorm) ! integer i,illa,ilma,ina,inoa,inva,ipoa real(dp) anorm - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + ipoa=c_idapo(ina) + anorm = 0.0_dp + if(c_nomax==1) illa=c_nmmax+1 + if(c_nomax==2) illa=combien + do i=ipoa,ipoa+illa-1 + anorm = anorm + abs(c_cc(i)) + enddo + return + endif + + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif return endif call dainf(ina,inoa,inva,ipoa,ilma,illa) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2725,6 +3662,428 @@ subroutine c_daabs(ina,anorm) end subroutine c_DAABS ! ! + subroutine dacctt1(mb,ib,mc,ic,ma,ia) + implicit none + ! *********************************** + ! + ! THIS SUBROUTINE PERFORMS A CONCATENATION MA = MB o MC + ! WHERE MA, MB AND MC ARE MATRICES CONSISTING OF IA, IB AND IC + ! DA VECTORS, RESPECTIVELY. + ! + !----------------------------------------------------------------------------- + ! + ! INTEGER MON(c_lno+1),Ic_cc(c_lnv) + ! INTEGER,dimension(:)::MB,MC,MA + !ETIENNE + !----------------------------------------------------------------------------- + ! + integer i0,i,j,k,ia,ib,ic,illa,illb,illc,ilma,ilmb,ilmc,inoa,inob,inoc,inva,invb,invc + + ! integer,dimension(c_lno)::icc + integer,dimension(:)::ma,mb,mc + integer ipoa(c_lnv),ipob(c_lnv),ipoc(c_lnv) + ! + !ETIENNE + ! + ! CONSISTENCY CHECKS + ! ****************** + ! + + + if(ib/=ic) then + write(6,*) " stopped in dacctt1 " + stop 880 + endif + +call dainf(ma(1),inoa,i0,ipoa(1),ilma,illa) + + + do i=1,i0 + call dainf(ma(i),inoa,inva,ipoa(i),ilma,illa) + call dainf(mb(i),inob,invb,ipob(i),ilmb,illb) + call dainf(mc(i),inoc,invc,ipoc(i),ilmc,illc) + + c_cc(ipoa(i):ipoa(i)+inva)=0.0_dp + enddo + + + do i=1,inva + c_cc(ipoa(i))=c_cc(ipob(i)) + enddo + + do i=1,inva + do j=1,inva + c_cc(ipoa(i))=c_cc(ipoa(i))+c_cc(ipob(i)+j)*c_cc(ipoc(j)) + enddo + + enddo + do i=1,inva + do j=1,inva + do k=1,inva + c_cc(ipoa(i)+k)=c_cc(ipoa(i)+k)+c_cc(ipob(i)+j)*c_cc(ipoc(j)+k) + enddo + enddo + enddo + + + return + end subroutine dacctt1 + + subroutine dacctt2tpsa(mb,ib,mc,ic,ma,ia) + implicit none + ! *********************************** + ! + ! THIS SUBROUTINE PERFORMS A CONCATENATION MA = MB o MC + ! WHERE MA, MB AND MC ARE MATRICES CONSISTING OF IA, IB AND IC + ! DA VECTORS, RESPECTIVELY. + ! + !----------------------------------------------------------------------------- + ! + ! INTEGER MON(c_lno+1),Ic_cc(c_lnv) + ! INTEGER,dimension(:)::MB,MC,MA + !ETIENNE + !----------------------------------------------------------------------------- + ! + integer i0,i,j,k,ia,ib,ic,illa,illb,illc,ilma,ilmb,ilmc,inoa,inob,inoc,inva,invb,invc,m,n + real(dp) fac + + ! integer,dimension(c_lno)::icc + integer,dimension(:)::ma,mb,mc + integer ipoa(c_lnv),ipob(c_lnv),ipoc(c_lnv) + complex(dp), allocatable :: g0(:),h0(:), G(:,:),H(:,:), Ga(:,:,:), He(:,:,:) + complex(dp), allocatable :: f0(:),F(:,:),Fi(:,:,:) + + ! + !ETIENNE + ! + ! CONSISTENCY CHECKS + ! ****************** + ! +call dainf(ma(1),inoa,i0,ipoa(1),ilma,illa) + +allocate(g0(i0),h0(i0),f0(i0), G(i0,i0),H(i0,i0),F(i0,i0), Ga(i0,i0,i0), He(i0,i0,i0), fi(i0,i0,i0)) + g0=0 + h0=0 + f0=0 + g=0 + h=0 + f=0 + ga=0 + he=0 + fi=0 + + do i=1,i0 + call dainf(ma(i),inoa,inva,ipoa(i),ilma,illa) + call dainf(mb(i),inob,invb,ipob(i),ilmb,illb) + call dainf(mc(i),inoc,invc,ipoc(i),ilmc,illc) + if(inva/=i0) stop 1234 + if(invb/=i0) stop 1235 + if(invc/=i0) stop 1236 + c_cc(ipoa(i):ipoa(i)+c_nmmax-1)=0.0_dp + + g0(i)=c_cc(ipob(i)) + h0(i)=c_cc(ipoc(i)) + do j=1,i0 + g(i,j)=c_cc(ipob(i)+j) + h(i,j)=c_cc(ipoc(i)+j) + do k=1,i0 + fac=2.0_dp + if(j==k) fac=1.0_dp + ga(i,j,k)=c_cc(ipob(i)+inds(j,k)-1)/fac + he(i,j,k)=c_cc(ipoc(i)+inds(j,k)-1)/fac + enddo + enddo + enddo + + f0=g0 + + do i=1,i0 + do j=1,i0 + f0(i)=g(i,j)*h0(j) +f0(i) + + do k=1,i0 + f0(i)=ga(i,j,k)*h0(j)*h0(k) + f0(i) + f(i,k)=f(i,k)+g(i,j)*H(j,k) + do m=1,i0 + f(i,k)=f(i,k)+ga(i,j,m)*H(j,k)*h0(m)+ga(i,j,m)*H(m,k)*h0(j) + enddo + enddo + enddo + enddo + + + do i=1,i0 + do j=1,i0 + do n=1,i0 + do m=1,i0 + fi(i,m,n)=fi(i,m,n)+ g(i,j)*he(j,m,n) + do k=1,i0 + fi(i,m,n)=fi(i,m,n)+ ga(i,j,k)*he(k,m,n)*h0(j)+ ga(i,j,k)*he(j,m,n)*h0(k) + fi(i,m,n)=fi(i,m,n)+ ga(i,j,k)*h(j,m)*H(k,n) + enddo + enddo + enddo + enddo + enddo + + + do i=1,i0 + c_cc(ipoa(i))=f0(i) + c_cc(ipoa(i)+1:ipoa(i)+i0)=f(i,1:i0) + enddo + + + do k=1,i0 + do j=poscombien, combien + fac=2.0_dp + if(ind1(j)==ind2(j)) fac=1.0_dp + c_cc(ipoa(k)+j-1)=fi(k,ind1(j),ind2(j))*fac + enddo + enddo + + + +deallocate(g0,h0, G,H, Ga, He,f0,f,fi) + + return + end subroutine dacctt2tpsa + + subroutine dacctt2da(mb,ib,mc,ic,ma,ia) + implicit none + ! *********************************** + ! + ! THIS SUBROUTINE PERFORMS A CONCATENATION MA = MB o MC + ! WHERE MA, MB AND MC ARE MATRICES CONSISTING OF IA, IB AND IC + ! DA VECTORS, RESPECTIVELY. + ! + !----------------------------------------------------------------------------- + ! + ! INTEGER MON(c_lno+1),Ic_cc(c_lnv) + ! INTEGER,dimension(:)::MB,MC,MA + !ETIENNE + !----------------------------------------------------------------------------- + ! + integer i,j,k,ia,ib,ic,illa,illb,illc,ilma,ilmb,ilmc,inoa,inob,inoc,inva,invb,invc + integer ab,jk,ik + + ! integer,dimension(c_lno)::icc + integer,dimension(:)::ma,mb,mc + integer ipoa(c_lnv),ipob(c_lnv),ipoc(c_lnv) + complex(dp)xx + ! + !ETIENNE + ! + ! CONSISTENCY CHECKS + ! ****************** + ! + + + if(ib/=ic) then + write(6,*) " stopped in dacctt1 " + stop 880 + endif + + + + do i=1,ia + call dainf(ma(i),inoa,inva,ipoa(i),ilma,illa) + c_cc(ipoa(i):ipoa(i)+combien-1)=0.0_dp + enddo + + do i=1,ib + call dainf(mb(i),inob,invb,ipob(i),ilmb,illb) + enddo + + do i=1,ic + call dainf(mc(i),inoc,invc,ipoc(i),ilmc,illc) + enddo + + + do i=invc-nphere+1,invc + c_cc(ipoa(i)+i)=1.0_dp + enddo +!c_cc(ipoa(i)+inva-nphere+1:ipoa(i)+inva)=1.0_dp + + do i=1,invc-nphere + do j=1,invc-nphere !? + do k=1,invc + c_cc(ipoa(i)+k)=c_cc(ipoa(i)+k)+c_cc(ipob(i)+j)*c_cc(ipoc(j)+k) + enddo + enddo + + do j=invc-nphere+1,invc +! do k=1,inva + c_cc(ipoa(i)+j)=c_cc(ipoa(i)+j)+c_cc(ipob(i)+j)*c_cc(ipoc(j)+j) +! enddo + enddo + + enddo + + + do i=1,invc-nphere + do j=1,invc -nphere !? + do jk=poscombien ,combien + c_cc(ipoa(i)+jk-1 )=c_cc(ipoa(i)+jk-1 ) + c_cc(ipob(i)+j)*c_cc(ipoc(j)+jk-1) + enddo + enddo + enddo + +if(with_para==0) then + do i=1,invc-nphere + do jk=poscombien,combien + do ab=poscombien,combien + c_cc(ipoa(i)+ab-1 )=c_cc(ipoa(i)+ab-1 ) + finds(ind1(jk),ind2(jk),ab)* & + c_cc(ipob(i)+jk-1 )*(c_cc(ipoc(ind1(jk))+ind1(ab))*c_cc(ipoc(ind2(jk)) & + +ind2(ab))+c_cc(ipoc(ind1(jk))+ind2(ab))*c_cc(ipoc(ind2(jk))+ind1(ab)))/2.0_dp + enddo + enddo + enddo + elseif(with_para==2) then + + + do i=1,invc-nphere + do ik=1,ninds + ab=nind1(ik) + jk=nind2(ik) + c_cc(ipoa(i)+ab-1 )=c_cc(ipoa(i)+ab-1 ) + finds1(ik)* & + c_cc(ipob(i)+jk-1 )*(c_cc(ipoc(ind1(jk))+ind1(ab))*c_cc(ipoc(ind2(jk))+ind2(ab)) & + +c_cc(ipoc(ind1(jk))+ind2(ab))*c_cc(ipoc(ind2(jk))+ind1(ab)))/2.0_dp + + enddo + enddo +elseif(with_para==1) then + do i=1,invc-nphere + do ik=1,ninds + ab=nind1(ik) + jk=nind2(ik) + c_cc(ipoa(i)+ab-1 )=c_cc(ipoa(i)+ab-1 ) + finds(ind1(jk),ind2(jk),ab)* & + c_cc(ipob(i)+jk-1 )*(c_cc(ipoc(ind1(jk))+ind1(ab))*c_cc(ipoc(ind2(jk))+ind2(ab)) & + +c_cc(ipoc(ind1(jk))+ind2(ab))*c_cc(ipoc(ind2(jk))+ind1(ab)))/2.0_dp + + enddo + enddo +endif + +do i=1,inva + c_cc(ipoa(i))=c_cc(ipob(i)) +enddo + + return + end subroutine dacctt2da + + + subroutine c_dacctt2datest(cnd2) + implicit none + + integer i,k,kp,knp,cnd2 + integer ab,jk + logical trash,trash1,trash2,trash3,trash4 + + +i=0 +k=0 +kp=0 +knp=0 + do jk=poscombien,combien + do ab=poscombien,combien +i=i+1 + +if(ind1(jk)>cnd2.or.ind2(jk)>cnd2) then +trash=.false. +k=k+1 +!write(6,*) i,jk,ab +!write(6,*) ind1(jk),ind2(jk),ind1(ab),ind2(ab) +if(ind1(jk)>cnd2) then + trash1=(ind1(jk)/=ind1(ab)) + trash3=(ind1(jk)/=ind2(ab)) +endif +if(ind2(jk)>cnd2) then + trash2=(ind2(jk)/=ind2(ab)) + trash4=(ind2(jk)/=ind1(ab)) +endif +trash=(trash1.or.trash2).and.(trash3.or.trash4) +!pause +endif +if(trash) then + kp=kp+1 +else + knp=knp+1 +endif +! c_cc(ipoa(i)+ab-1 )=c_cc(ipoa(i)+ab-1 ) + finds(ind1(jk),ind2(jk),ab)* & + ! c_cc(ipob(i)+jk-1 )*(c_cc(ipoc(ind1(jk))+ind1(ab))*c_cc(ipoc(ind2(jk)) & +! +ind2(ab))+c_cc(ipoc(ind1(jk))+ind2(ab))*c_cc(ipoc(ind2(jk))+ind1(ab)))/2.0_dp + enddo + enddo + +i=0 +k=0 +kp=0 +ninds=knp +knp=0 +allocate(nind1(ninds),nind2(ninds)) + +nind1=0 +nind2=0 + + + do jk=poscombien,combien + do ab=poscombien,combien +i=i+1 + +if(ind1(jk)>cnd2.or.ind2(jk)>cnd2) then +trash=.false. +k=k+1 +!write(6,*) i,jk,ab +!write(6,*) ind1(jk),ind2(jk),ind1(ab),ind2(ab) +if(ind1(jk)>cnd2) then + trash1=(ind1(jk)/=ind1(ab)) + trash3=(ind1(jk)/=ind2(ab)) +endif +if(ind2(jk)>cnd2) then + trash2=(ind2(jk)/=ind2(ab)) + trash4=(ind2(jk)/=ind1(ab)) +endif +trash=(trash1.or.trash2).and.(trash3.or.trash4) +!pause +endif +if(trash) then + kp=kp+1 +else + knp=knp+1 + nind1(knp)=ab + nind2(knp)=jk +endif + enddo + enddo + + + +!allocate(finds(nvhere,nvhere,poscombien:combien )) +if(with_para==2) then +write(6,*) combien,ninds,size(nind1), size(nind2) +allocate(finds1(ninds)) +finds1=1 + do i=1,ninds + if(ind1(nind1(i))/=ind2(nind1(i))) finds1(i)=2 + enddo +endif +!write(6,*) size(finds1),size(ind1), size(ind2) + +!pause 777 + + ! do ik=1,ninds +! ab=nind1(ik) +! jk=nind2(ik) +!finds1(i) + ! c_cc(ipoa(i)+ab-1 )=c_cc(ipoa(i)+ab-1 ) + finds(ind1(jk),ind2(jk),ab)* & + ! c_cc(ipoa(i)+ab-1 )=c_cc(ipoa(i)+ab-1 ) + finds1(ik)* & + + + return + end subroutine c_dacctt2datest + + + subroutine c_dacct(ma,ia,mb,ib,mc,ic) implicit none ! *********************************** @@ -2736,9 +4095,95 @@ subroutine c_dacct(ma,ia,mb,ib,mc,ic) !----------------------------------------------------------------------------- ! integer i,ia,ib,ic,ij,illc,ilmc,inoc,invc,ipoc - integer,dimension(c_lnv)::monx + integer,dimension(c_lnv)::monx,maa,mcc integer,dimension(:)::ma,mb,mc - if((.not.C_STABLE_DA)) then + + +if(newtpsa) then + if(ib/=c_nvmax) then + write(6,*) "problem in dacct " + stop 9 + endif + if(ia/=ic) then + write(6,*) "problem in dacct " + stop 10 + endif + maa=0 + mcc=0 + do i=1,ia + maa(i)=ma(i) + enddo + + if((mc(1)==ma(1)).or.(mc(1)==mb(1))) then + do i=1,ic + call c_daall0(mcc(i)) + enddo + else + do i=1,ic + mcc(i)=mc(i) + enddo + endif + + do i=ia+1,c_nvmax + call c_daall0(maa(i)) + enddo + do i=ic+1,c_nvmax + call c_daall0(mcc(i)) + enddo + + + + + if(c_nomax==1) then + call dacctt1(maa,c_nvmax,mb,c_nvmax,mcc,c_nvmax) + + if((mc(1)==ma(1)).or.(mc(1)==mb(1))) then + do i=1,ic + call c_dacop(mcc(i),mc(i)) + call c_dadal1(mcc(i)) + enddo + endif + + do i=ia+1,c_nvmax + call c_dadal1(maa(i)) + enddo + do i=ic+1,c_nvmax + call c_dadal1(mcc(i)) + enddo + + return + endif + + if(c_nomax==2) then + if(assume_da_map) then + + + call dacctt2da(maa,c_nvmax,mb,c_nvmax,mcc,c_nvmax) + + + + else + call dacctt2tpsa(maa,c_nvmax,mb,c_nvmax,mcc,c_nvmax) + endif + if((mc(1)==ma(1)).or.(mc(1)==mb(1))) then + do i=1,ic + call c_dacop(mcc(i),mc(i)) + call c_dadal1(mcc(i)) + enddo + endif + do i=ia+1,c_nvmax + call c_dadal1(maa(i)) + enddo + do i=ic+1,c_nvmax + call c_dadal1(mcc(i)) + enddo + return + endif +endif + + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2747,7 +4192,7 @@ subroutine c_dacct(ma,ia,mb,ib,mc,ic) ! if(ma(1).eq.mc(1).or.mb(1).eq.mc(1)) then call dainf(mc(1),inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2790,12 +4235,24 @@ subroutine dacctt(mb,ib,mc,ic,ma,ia) integer,dimension(c_lnv)::icc !johan 2008 March integer,dimension(:)::ma,mb,mc complex(dp) ccf - if((.not.C_STABLE_DA)) then - if(C_watch_user) then - write(6,*) "big problem in dabnew ", sqrt(crash) - endif - return + + + if(c_nomax.eq.1.and.newtpsa) then + return + endif + + if(newtpsa) then + if(c_nomax==1) then + stop 7501 + return + endif + + if(c_nomax==2) then + stop 7511 + return endif +endif + ! !ETIENNE ! @@ -2808,7 +4265,7 @@ subroutine dacctt(mb,ib,mc,ic,ma,ia) call dainf(iia,inoa,inva,ipoa,ilma,illa) call dainf(iib,inob,invb,ipob,ilmb,illb) call dainf(iic,inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2891,7 +4348,20 @@ subroutine c_mtree(mb,ib,mc,ic) integer,dimension(0:c_lno)::jv integer,dimension(:)::mb,mc complex(dp) apek,bbijj,chkjj - if((.not.C_%STABLE_DA)) then + + if(newtpsa) then + if(c_nomax==1) then + stop 8881 + return + endif + + if(c_nomax==2) then + stop 8891 + return + endif +endif + + if((.not.c_stable_da)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -2905,7 +4375,7 @@ subroutine c_mtree(mb,ib,mc,ic) iic = mc(1) call dainf(iib,inob,invb,ipob,ilmb,illb) call dainf(iic,inoc,invc,ipoc,ilmc,illc) - if((.not.C_%STABLE_DA)) then + if((.not.c_stable_da)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3099,7 +4569,20 @@ subroutine ppushprint(mc,ic,mf,jc,line) integer i,ic,iv,jc,jl,jv,mf integer,dimension(:)::mc character(20) line - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + if(c_nomax==1) then + stop 7881 + return + endif + + if(c_nomax==2) then + stop 7891 + return + endif +endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3135,6 +4618,19 @@ subroutine ppushstore(mc,nd2,coef,ml,mv) integer,dimension(:), intent(in)::mc integer,dimension(:), intent(out)::ml,mv complex(dp),dimension(:),intent(out)::coef + + if(newtpsa) then + if(c_nomax==1) then + stop 8183 + return + endif + + if(c_nomax==2) then + stop 8194 + return + endif +endif + ! jc=0 jl=0;jv=0; @@ -3163,6 +4659,7 @@ subroutine ppushstore(mc,nd2,coef,ml,mv) enddo return end subroutine ppushstore + subroutine ppushGETN(mc,ND2,ntot) implicit none ! @@ -3187,7 +4684,21 @@ subroutine ppush(mc,ic,xi,xf) complex(dp),dimension(c_lno+1)::xm complex(dp),dimension(c_lno)::xt complex(dp),dimension(:)::xf,xi - if((.not.C_STABLE_DA)) then + + + if(newtpsa) then + if(c_nomax==1) then + stop 8828 + return + endif + + if(c_nomax==2) then + stop 8839 + return + endif +endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3216,13 +4727,14 @@ subroutine ppush(mc,ic,xi,xf) enddo do i=1,c_nvmax if(abs(xf(i))>da_absolute_aperture.and.check_da) then - C_STABLE_DA=.false. + c_stable_da=.false. write(6,*) "unstable in ppush ",i,xf(i) endif enddo ! return end subroutine ppush + subroutine ppush1(mc,xi,xf) implicit none ! ***************************** @@ -3237,7 +4749,20 @@ subroutine ppush1(mc,xi,xf) complex(dp),dimension(:)::xi complex(dp),dimension(c_lno)::xt complex(dp),dimension(c_lno+1)::xm - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + if(c_nomax==1) then + stop 8389 + return + endif + + if(c_nomax==2) then + stop 8390 + return + endif +endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3261,13 +4786,203 @@ subroutine ppush1(mc,xi,xf) xf = xf + c_cc(c_idapo(mc)+i) * xx enddo if(abs(xf)>da_absolute_aperture.and.check_da) then - C_STABLE_DA=.false. + c_stable_da=.false. write(6,*) "unstable in ppush1 ", xf write(6,*) xi(1:c_nvmax) endif return end subroutine ppush1 + +subroutine dainvt1(mb,ib,mc,ic) + implicit none + ! *********************************** + ! + ! THIS SUBROUTINE PERFORMS A INVERSION first or second order + !ETIENNE + !----------------------------------------------------------------------------- + ! + integer i0,i,j,k,ib,ic,illa,illb,illc,ilmb,ilmc,inob,inoc,invb,invc,m,n + real(dp) fac + + ! integer,dimension(lno)::icc + integer,dimension(:):: mb,mc + integer ipob(c_lnv),ipoc(c_lnv),ier + complex(dp), allocatable :: g0(:),h0(:), G(:,:),H(:,:) + + integer, allocatable :: icc(:),ipo(:) + ! + !ETIENNE + ! + ! CONSISTENCY CHECKS + ! ****************** + ! + + +call dainf(mb(1),inob,i0,ipob(1),ilmb,illb) + + +allocate(g0(i0),h0(i0) , G(i0,i0),H(i0,i0) ) + g0=0 + h0=0 + g=0 + h=0 + + + do i=1,i0 + call dainf(mc(i),inoc,invc,ipoc(i),ilmc,illc) + call dainf(mb(i),inob,invb,ipob(i),ilmb,illb) + if(invb/=i0) stop 1235 + if(invc/=i0) stop 1236 + + g0(i)=c_cc(ipob(i)) + do j=1,i0 + g(i,j)=c_cc(ipob(i)+j) + + enddo + enddo + + call c_matinv(g,h,i0,i0,ier) + + h0=-matmul(h,g0) + + + do i=1,i0 + c_cc(ipoc(i))=h0(i) + c_cc(ipoc(i)+1:ipoc(i)+i0)=H(i,1:i0) + enddo + + + +end subroutine dainvt1 + + +subroutine dainvt2(mb,ib,mc,ic) + implicit none + ! *********************************** + ! + ! THIS SUBROUTINE PERFORMS A INVERSION first or second order + !ETIENNE + !----------------------------------------------------------------------------- + ! + integer i0,i,j,k,ib,ic,illa,illb,illc,ilmb,ilmc,inob,inoc,invb,invc,m,n + real(dp) fac + + ! integer,dimension(lno)::mcc + integer,dimension(:):: mb,mc + integer ipob(c_lnv),ipoc(c_lnv),ier + complex(dp), allocatable :: g0(:),h0(:), G(:,:),H(:,:), Ga(:,:,:), He(:,:,:),r0(:), r1(:,:) + + integer, allocatable :: mcc(:),ipo(:),mccc(:),ipoo(:) + ! + !ETIENNE + ! + ! CONSISTENCY CHECKS + ! ****************** + ! + +call dainf(mb(1),inob,i0,ipob(1),ilmb,illb) + + +allocate(r0(i0),g0(i0),h0(i0) ,r1(i0,i0), G(i0,i0),H(i0,i0) , Ga(i0,i0,i0), He(i0,i0,i0) ) + g0=0 + h0=0 + g=0 + h=0 + ga=0 + he=0 + + do i=1,i0 + call dainf(mc(i),inoc,invc,ipoc(i),ilmc,illc) + call dainf(mb(i),inob,invb,ipob(i),ilmb,illb) + if(invb/=i0) stop 1235 + if(invc/=i0) stop 1236 + + g0(i)=c_cc(ipob(i)) + do j=1,i0 + g(i,j)=c_cc(ipob(i)+j) + if(inob==2) then + do k=1,i0 + fac=2.0_dp + if(j==k) fac=1.0_dp + ga(i,j,k)=c_cc(ipob(i)+inds(j,k)-1)/fac + enddo + endif + enddo + enddo + + call c_matinv(g,h,i0,i0,ier) + + ! h0=-matmul(h,g0) + + + + allocate(mcc(i0),ipo(i0)) + allocate(mccc(i0),ipoo(i0)) + + call daall(mcc,i0,'$$DACCT $$',c_nomax,c_nvmax) + call daall(mccc,i0,'$$DACCT $$',c_nomax,c_nvmax) + + do i=1,i0 + call dainf(mcc(i),inoc,invc,ipo(i),ilmc,illc) + call dainf(mccc(i),inoc,invc,ipoo(i),ilmc,illc) + c_cc(ipo(i):ipo(i)+combien-1)=0 + c_cc(ipo(i))=h0(i) + c_cc(ipo(i)+1:ipo(i)+i0)=H(i,1:i0) + enddo + + + if(assume_da_map) then + call dacctt2da(mb,i0,mcc,i0,mccc,i0) + else + call dacctt2tpsa(mb,i0,mcc,i0,mccc,i0) + endif + + + do i=1,i0 + c_cc(ipoo(i)+poscombien-1:ipoo(i)+combien-1)=-c_cc(ipoo(i)+poscombien-1:ipoo(i)+combien-1) + enddo + + if(assume_da_map) then + call dacctt2da(mcc,i0,mccc,i0,mc,ic) + else + call dacctt2tpsa(mcc,i0,mccc,i0,mc,ic) + endif + + + do i=1,i0 + c_cc(c_idapo(mcc(i)):c_idapo(mcc(i))+c_nmmax-1 )=0 + c_cc(c_idapo(mccc(i)):c_idapo(mccc(i))+c_nmmax-1 )=0 + enddo + + do i=1,i0 + c_cc(c_idapo(mc(i)):c_idapo(mc(i)) ) = 0.0_dp + c_cc(c_idapo(mcc(i)):c_idapo(mcc(i)) ) = -g0(i) + c_cc(c_idapo(mcc(i))+i)=1.0_dp + enddo + + if(assume_da_map) then + call dacctt2da(mc,ic,mcc,i0,mccc,i0) + else + call dacctt2tpsa(mc,ic,mcc,i0,mccc,i0) + endif + + + do i=1,i0 + c_cc(c_idapo(mc(i)):c_idapo(mc(i)) + c_nmmax-1 ) = c_cc(c_idapo(mccc(i)):c_idapo(mccc(i)) + c_nmmax-1 ) + enddo + + call dadal(mcc,i0) + call dadal(mccc,i0) + deallocate(mcc) + deallocate(mccc) + deallocate(ipo) + deallocate(ipoo) + +deallocate(g0,h0, G,H, Ga, He) + +end subroutine dainvt2 + subroutine c_dainv(ma,ia,mb,ib) implicit none ! ***************************** @@ -3281,7 +4996,23 @@ subroutine c_dainv(ma,ia,mb,ib) integer,dimension(c_lnv)::jj,ml integer,dimension(:)::ma,mb complex(dp),dimension(c_lnv)::x - if((.not.C_STABLE_DA)) then + + + if(newtpsa) then + if(c_nomax==1) then + call dainvt1(ma,ia,mb,ib) + return + endif + + if(c_nomax==2) then + + call dainvt2(ma,ia,mb,ib) + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3293,7 +5024,7 @@ subroutine c_dainv(ma,ia,mb,ib) enddo if(ma(1).eq.mb(1)) then call dainf(mb(1),inob,invb,ipob,ilmb,illb) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3306,7 +5037,9 @@ subroutine c_dainv(ma,ia,mb,ib) ml(ij)=0 enddo call daall(ml,ib,'$$DAJUNK$$',inob,invb) + call dainvt(ma,ia,ml,ib) + do i=1,ib call c_dacop(ml(i),mb(i)) enddo @@ -3316,7 +5049,9 @@ subroutine c_dainv(ma,ia,mb,ib) call c_dapek(ma(i),jj,x(i)) call c_dapok(ma(i),jj,(0.0_dp,0.0_dp)) enddo + call dainvt(ma,ia,mb,ib) + do i=1,ia call c_dapok(ma(i),jj,x(i)) enddo @@ -3341,7 +5076,21 @@ subroutine dainvt(ma,ia,mb,ib,success) complex(dp) amjj,amsjj,prod logical, optional :: success if(present(success)) success=.true. - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + if(c_nomax==1) then + stop 1780 + return + endif + + if(c_nomax==2) then + stop 1781 + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3350,7 +5099,7 @@ subroutine dainvt(ma,ia,mb,ib,success) ! call dainf(ma(1),inoa,inva,ipoa,ilma,illa) call dainf(mb(1),inob,invb,ipob,ilmb,illb) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3396,9 +5145,12 @@ subroutine dainvt(ma,ia,mb,ib,success) ! EXTRACTING LINEAR MATRIX, GENERATING NONLINEAR PART OF A ! ******************************************************** ! + + do i=1,ib do j=1,ib - do k=1,ib +! do k=1,ib + do k=1,size(jj) jj(k) = 0 enddo jj(j) = 1 @@ -3408,6 +5160,9 @@ subroutine dainvt(ma,ia,mb,ib,success) enddo call c_dacmu(ma(i),-(1.0_dp,0.0_dp),ma(i)) enddo + + + ! ! INVERTING LINEAR MATRIX, CHECKING RESULT AND STORING IN ML ! ********************************************************** @@ -3418,7 +5173,7 @@ subroutine dainvt(ma,ia,mb,ib,success) if(present(success)) success=.false. if(check_da) then - C_STABLE_DA=.false. + c_stable_da=.false. C_check_stable=.false. messagelost='ERROR IN ROUTINE DAINV, ier=132 in matinv' write(6,*) messagelost @@ -3446,7 +5201,7 @@ subroutine dainvt(ma,ia,mb,ib,success) write(6,*) " abs(prod) > 100.0_dp*epsmac in dainvt",abs(prod), 100.0_dp*epsmac if(check_da) then - C_STABLE_DA=.false. + c_stable_da=.false. messagelost='ERROR IN ROUTINE DAINV, abs(prod).gt.100.0_dp*epsmac ' call dadal(ml,ia) call dadal(ms,ia) @@ -3529,8 +5284,9 @@ subroutine dapin(ma,ia,mb,ib,jx) integer,dimension(c_lnv)::jj,ml integer,dimension(:)::ma,mb,jx complex(dp),dimension(c_lnv)::x - ! - if((.not.C_STABLE_DA)) then + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3541,7 +5297,7 @@ subroutine dapin(ma,ia,mb,ib,jx) enddo if(ma(1).eq.mb(1)) then call dainf(mb(1),inob,invb,ipob,ilmb,illb) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3585,14 +5341,17 @@ subroutine c_dapint(ma,ia,mb,ib,jind) integer,dimension(c_lnv)::jj,mn,mi,me integer,dimension(:)::ma,mb,jind ! - if((.not.C_STABLE_DA)) then + + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif return endif call dainf(ma(1),inoa,inva,ipoa,ilma,illa) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3609,8 +5368,9 @@ subroutine c_dapint(ma,ia,mb,ib,jind) call daall(mi,ia,'$$PIN2 $$',inoa,inva) call daall(me,ia,'$$PIN3 $$',inoa,inva) ! + do i=1,ia - do k=1,c_nvmax + do k=1,size(jj) jj(k) = 0 enddo jj(i) = 1 @@ -3647,9 +5407,57 @@ subroutine c_dader(idif,ina,inc) ! !----------------------------------------------------------------------------- ! - integer idif,illc,ilmc,ina,inc,incc,inoc,invc,ipoc - ! - if((.not.C_STABLE_DA)) then + integer idif,illc,ilmc,ina,inc,incc,inoc,invc,ipoc,inoa,inva,ipoa,ilma,illa + integer i,jd(c_lnv),ic1,ic2 + complex(dp) rr + +if(newtpsa) then + call dainf(ina,inoa,inva,ipoa,ilma,illa) + call dainf(inc,inoc,invc,ipoc,ilmc,illc) + if(c_nomax.eq.1) then + ! PRINT*,'ERROR, DADER CALLED WITH c_nomax = 1' + ! call dadeb !(31,'ERR DADER1',1) + ! stop 666 + do i=1,c_lnv + jd(i)=0 + enddo + jd(idif)=1 + call c_dapek(ina,jd,rr) + call c_dacon(inc,rr) + return + endif + if(c_nomax.eq.2) then + reelc=0 + + + do i=1,combien + + ic1=ind1(i) + ic2=ind2(i) + + if(ic1==ic2) then + if(ic1==idif) then + ic1=0 + reelc(inds(ic1,ic2))=2.0_dp*c_cc(ipoa+i-1) + endif + elseif(ic1==idif) then + ic1=0 + reelc(inds(ic1,ic2))=c_cc(ipoa+i-1) + elseif(ic2==idif) then + ic2=0 + reelc(inds(ic1,ic2))=c_cc(ipoa+i-1) + endif + + enddo + + c_cc(ipoc:ipoc+combien-1) = reelc + + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3657,7 +5465,7 @@ subroutine c_dader(idif,ina,inc) endif if(ina.eq.inc) then call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3673,6 +5481,7 @@ subroutine c_dader(idif,ina,inc) endif return end subroutine c_dader + subroutine dadert(idif,ina,inc) implicit none ! ****************************** @@ -3686,8 +5495,22 @@ subroutine dadert(idif,ina,inc) ilma,ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj,ipause,mypauses integer,dimension(c_lnv)::jd complex(dp) rr,x,xdivi - ! - if((.not.C_STABLE_DA)) then + + + if(newtpsa) then + if(c_nomax==1) then + stop 1310 + return + endif + + if(c_nomax==2) then + stop 1311 + return + endif +endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3695,7 +5518,7 @@ subroutine dadert(idif,ina,inc) endif call dainf(ina,inoa,inva,ipoa,ilma,illa) call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3787,10 +5610,56 @@ subroutine c_dacfu(ina,fun,inc) ! !----------------------------------------------------------------------------- ! - integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc - complex(dp),external::fun + integer illc,ilmc,ina,inc,incc,inoc,invc,ipoc,inoa,inva,ipoa,ilma,illa,i,j(c_lnv) + + ! real(dp) rr,x,xdivi ! - if((.not.C_STABLE_DA)) then + interface + ! real(kind(one)) function fun(abc) + function fun(abc) + use precision_constants + implicit none + complex(dp) fun + integer,dimension(:)::abc + end function fun + end interface + + if(newtpsa) then + call dainf(ina,inoa,inva,ipoa,ilma,illa) + call dainf(inc,inoc,invc,ipoc,ilmc,illc) + + if(c_nomax==1) then + do i=1,c_lnv + j(i)=0 + enddo + + do i=0,c_nvmax + if(i/=0) j(i)=1 + c_cc(ipoc+i)=fun(j)*c_cc(ipoa+i) + if(i/=0) j(i)=0 + enddo + return + endif + + if(c_nomax==2) then + do i=0,c_nvmax + if(i/=0) j(i)=1 + c_cc(ipoc+i)=fun(j)*c_cc(ipoa+i) + if(i/=0) j(i)=0 + enddo + do i=poscombien,combien + j=0 + j(ind1(i))=j(ind1(i))+1 + j(ind2(i))=j(ind2(i))+1 + c_cc(ipoc+i-1)=fun(j)*c_cc(ipoa+i-1) + j(ind1(i)) = 0 + j(ind2(i)) = 0 + enddo + return + endif +endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3798,7 +5667,7 @@ subroutine c_dacfu(ina,fun,inc) endif if(ina.eq.inc) then call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3841,8 +5710,17 @@ function fun(abc) integer,dimension(:)::abc end function fun end interface - ! - if((.not.C_STABLE_DA)) then + + if(c_nomax.eq.1.and.newtpsa) then + return + endif + + if(c_nomax==2.and.newtpsa) then + return + endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3850,7 +5728,7 @@ end function fun endif call dainf(ina,inoa,inva,ipoa,ilma,illa) call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3920,6 +5798,7 @@ end subroutine dacfut ! call dancd(c_i_1(I),c_i_2(I),J)! ! END subroutine GET_C_J ! END subroutine GET_C_J + subroutine c_dapri(ina,iunit) implicit none ! *************************** @@ -3935,7 +5814,31 @@ subroutine c_dapri(ina,iunit) logical long long=longprint if(iunit/=6) longprint=.true. - if((.not.C_STABLE_DA)) then + + + +if(newtpsa) then + if(c_nomax.eq.1) then + write(iunit,*) "1st order polynomial ", ina,c_nvmax + do i=1,c_nvmax+1 + if(c_cc(c_idapo(ina)+i-1).ne.0.0_dp) write(iunit,*) i-1, c_cc(c_idapo(ina)+i-1) + enddo + + return + endif + + if(c_nomax==2) then + write(iunit,*) "2nd order polynomial ", ina,c_nvmax + do i=1,combien + if(c_cc(c_idapo(ina)+i-1).ne.0.0_dp) write(iunit,*) ind1(i),ind2(i),c_cc(c_idapo(ina)+i-1) + enddo + return + endif +endif + + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -3943,7 +5846,7 @@ subroutine c_dapri(ina,iunit) endif if(ina.lt.1.or.ina.gt.c_nda_dab) then print*,'ERROR IN c_dapri, INA = ',ina - C_STABLE_DA=.false. + c_stable_da=.false. endif ! inoa = c_idano(ina) @@ -4033,19 +5936,69 @@ subroutine c_dapri77(ina,iunit) ! !----------------------------------------------------------------------------- ! - integer i,ii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,iunit + integer i,ii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,iunit,count integer,dimension(c_lnv)::j character(10) c10,k10 real(dp) a,b complex(dp) ccc logical some,imprime logical long + integer(2) xi(3) long=longprint if(iunit/=6) longprint=.true. some=.false. - ! + +if(newtpsa) then + if(c_nomax.eq.1) then + count=0 + write(iunit,*) "1st order polynomial ", ina,c_nvmax + do i=1,c_nvmax+1 + if(c_cc(c_idapo(ina)+i-1).ne.0.0_dp) then + xi(1)=i-1 + xi(2)=ind1(i) + xi(3)=ind2(i) + if(xi(2)>c_nvmax-nphere) then + xi(2)=-(xi(2)-(c_nvmax-nphere)) + endif + if(xi(3)>c_nvmax-nphere) then + xi(3)=-(xi(3)-(c_nvmax-nphere)) + endif + write(iunit,*) xi, c_cc(c_idapo(ina)+i-1) + count=count+1 + endif + enddo + write(iunit,*) count, " monomial(s) printed " + return + endif + + if(c_nomax==2) then + count=0 + + write(iunit,*) "2nd order polynomial ", ina,c_nvmax + do i=1,combien + if(c_cc(c_idapo(ina)+i-1).ne.0.0_dp) then + xi(1)=i-1 + xi(2)=ind1(i) + xi(3)=ind2(i) + if(xi(2)>c_nvmax-nphere) then + xi(2)=-(xi(2)-(c_nvmax-nphere)) + endif + if(xi(3)>c_nvmax-nphere) then + xi(3)=-(xi(3)-(c_nvmax-nphere)) + endif + write(iunit,*) xi, c_cc(c_idapo(ina)+i-1) + count=count+1 + endif + enddo + write(iunit,*) count, " monomial(s) printed " + return + endif +endif + + + if(iunit.eq.0) return - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -4053,7 +6006,7 @@ subroutine c_dapri77(ina,iunit) endif if(ina.lt.1.or.ina.gt.c_nda_dab) then write(line,'(a22,i8)') 'ERROR IN c_dapri, INA = ',ina - C_STABLE_DA=.false. + c_stable_da=.false. endif ! inoa = c_idano(ina) @@ -4151,10 +6104,76 @@ subroutine c_dashift(ina,inc,ishift) !----------------------------------------------------------------------------- ! integer i,ii,illa,ilma,ina,inoa,inva,ioa,iout,ipoa,inb,ishift,ich,& - ik,inc,k + ik,inc,k,ipob,ic1,ic2 integer,dimension(c_lnv)::j,jd - ! - if((.not.C_STABLE_DA)) then + + if(newtpsa) then + inva = c_idanv(ina) + ipoa = c_idapo(ina) + + if(c_nomax==1) then + call daall1(inb,'$$DAJUNK$$',inoa,inva) + ipob = c_idapo(inb) + + c_cc(ipob:ipob+inva)=0 + + do ii=1,inva + if(c_cc(ipoa+ii)/=0.0_dp.and.ii<=ishift) then + write(6,*) " trouble in dashift " + stop 886 + else + c_cc(ipob+ii-ishift)=c_cc(ipoa+ii) + endif + enddo + c_cc(ipob)=c_cc(ipoa) + call c_dacop(inb,inc) + call c_dadal1(inb) + return + endif + + if(c_nomax==2) then + call daall1(inb,'$$DAJUNK$$',inoa,inva) + ipob = c_idapo(inb) + + c_cc(ipob:ipob+combien-1)=0 + + do ii=1,inva !ipoa,ipoa+illa-1 + if(c_cc(ipoa+ii)/=0.0_dp.and.ii<=ishift) then + write(6,*) " trouble in dashift " + stop 887 + else + c_cc(ipob+ii-ishift)=c_cc(ipoa+ii) + endif + enddo + + do ii=poscombien,combien !ipoa,ipoa+illa-1 + + + if( ((ind1(ii)<=ishift).or.(ind2(ii)<=ishift)).and.c_cc(ipoa+ii-1)/=0.0_dp) then + write(6,*) " trouble in dashift " + stop 889 + else + if(((ind1(ii)>ishift).or.(ind1(ii)==0)).and.((ind2(ii)>ishift).or.(ind2(ii)==0))) then + ic1=0;ic2=0; + if(ind1(ii)>ishift) ic1 = ind1(ii)-ishift + if(ind2(ii)>ishift) ic2 = ind2(ii)-ishift + + c_cc(ipob+inds(ic1,ic2)-1)=c_cc(ipoa+ii-1) + + endif + endif + enddo + + c_cc(ipob)=c_cc(ipoa) + call c_dacop(inb,inc) + call c_dadal1(inb) + return + endif + + endif + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -4163,7 +6182,7 @@ subroutine c_dashift(ina,inc,ishift) inb=0 if(ina.lt.1.or.ina.gt.c_nda_dab) then write(line,'(a22,i8)') 'ERROR IN dashift, INA = ',ina - C_STABLE_DA=.false. + c_stable_da=.false. endif ! inoa = c_idano(ina) @@ -4265,8 +6284,17 @@ subroutine c_darea(ina,iunit) integer,dimension(c_lnv)::j complex(dp) c character(10) c10 + + if(c_nomax.eq.1.and.newtpsa) then + return + endif + + if(c_nomax==2.and.newtpsa) then + return + endif - if((.not.C_STABLE_DA)) then + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -4274,7 +6302,7 @@ subroutine c_darea(ina,iunit) endif ! if(ina.lt.1.or.ina.gt.c_nda_dab) then - C_STABLE_DA=.false. + c_stable_da=.false. print*,'ERROR IN DAREA, INA = ',ina ! X = SQRT(-ONE) ! PRINT*,X @@ -4376,10 +6404,19 @@ subroutine c_darea77(ina,iunit) real(dp) cr,ci character(10) c10,k10 complex ik + + if(c_nomax.eq.1.and.newtpsa) then + return + endif + + if(c_nomax==2.and.newtpsa) then + return + endif + ik=( 0.0_dp,1.0_dp ) ! - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -4387,7 +6424,7 @@ subroutine c_darea77(ina,iunit) endif if(ina.lt.1.or.ina.gt.c_nda_dab) then write(line,'(a22,i8)') 'ERROR IN darea77, INA = ',ina - C_STABLE_DA=.false. + c_stable_da=.false. endif ! inoa = c_idano(ina) @@ -4460,7 +6497,7 @@ subroutine dadeb !(iunit,c,istop) ! !etienne ! if(check_da) then - C_STABLE_DA=.false. + c_stable_da=.false. write(6,*) "big problem in complex dadeb ", sqrt(crash) return @@ -4610,14 +6647,14 @@ subroutine damch(iaa,ia) integer i,ia,illa,ilma,ino1,inoi,inv1,invi,ipoa,ipause,mypauses integer,dimension(:)::iaa ! - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", "big problem in dabnew ", sqrt(crash) endif return endif call dainf(iaa(1),ino1,inv1,ipoa,ilma,illa) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", "big problem in dabnew ", sqrt(crash) endif @@ -4626,7 +6663,7 @@ subroutine damch(iaa,ia) ! do i=2,ia call dainf(iaa(i),inoi,invi,ipoa,ilma,illa) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -4717,10 +6754,50 @@ subroutine c_datra(idif,ina,inc) !----------------------------------------------------------------------------- ! integer i,ibase,ic,ider1,ider1s,ider2,ider2s,idif,iee,ifac,illa,illc,ilma,& - ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj,ipause,mypauses + ilmc,ina,inc,inoa,inoc,inva,invc,ipoa,ipoc,jj,ipause,mypauses,ic1,ic2 complex(dp) x,xdivi - ! - if((.not.C_STABLE_DA)) then + + + if(newtpsa) then + if(c_nomax==1) then + call c_dader(idif,ina,inc) + return + endif + + if(c_nomax==2) then + reelc=0 + ipoa=c_idapo(ina) + ipoc=c_idapo(inc) + + do i=1,combien + + ic1=ind1(i) + ic2=ind2(i) + + if(ic1==ic2) then + if(ic1==idif) then + ic1=0 + reelc(inds(ic1,ic2))= c_cc(ipoa+i-1) + endif + elseif(ic1==idif) then + ic1=0 + reelc(inds(ic1,ic2))=c_cc(ipoa+i-1) + elseif(ic2==idif) then + ic2=0 + reelc(inds(ic1,ic2))=c_cc(ipoa+i-1) + endif + + enddo + + c_cc(ipoc:ipoc+combien-1) = reelc + + return + endif +endif + + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -4728,7 +6805,7 @@ subroutine c_datra(idif,ina,inc) endif call dainf(ina,inoa,inva,ipoa,ilma,illa) call dainf(inc,inoc,invc,ipoc,ilmc,illc) - if((.not.C_STABLE_DA)) then + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -4877,7 +6954,33 @@ subroutine c_dacycle(ina,ipresent,value,illa,j) integer ii,illa,ilma,ina,inoa,inva,iout,ipoa,ipresent integer,optional,dimension(:)::j complex(dp) value - ! + + + if(newtpsa) then + + ipoa=c_idapo(ina) + if(.not.present(j)) then + illa=c_nmmax + ! illa0$=0 + + ! do i=1,c_nmmax + ! if(abs(c_cc(ipoa+i-1))/=0) illa=illa+1 + ! enddo + else + + value=c_cc(ipoa+ipresent-1) + j=0 + + if(ind1(ipresent)/=0) j(ind1(ipresent))=1+j(ind1(ipresent)) + if(ind2(ipresent)/=0) j(ind2(ipresent))=1+j(ind2(ipresent)) + + + return + endif +endif + + + if(ina.lt.1.or.ina.gt.c_nda_dab) then write(line,'(a22,i8)') 'ERROR IN dacycle, INA = ',ina ipause=mypauses(39,line) @@ -4908,45 +7011,8 @@ subroutine c_dacycle(ina,ipresent,value,illa,j) ! ipresent=1+ipresent return end subroutine c_dacycle -!!!! new stuff lingyun - subroutine dacc_lean(ina) - implicit none - integer ipause, mypauses - ! *************************** - ! - ! THIS SUBROUTINE PRINTS THE DA VECTOR INA TO UNIT IUNIT. - ! - !----------------------------------------------------------------------------- - ! - integer ii,illa,ilma,ina,inoa,inva,iout,ipoa - real(dp) a,b - if(ina.lt.1.or.ina.gt.c_nda_dab) then - write(line,'(a22,i8)') 'ERROR IN dacycle, INA = ',ina - ipause=mypauses(39,line) - stop - endif - ! - inoa = c_idano(ina) - inva = c_idanv(ina) - ipoa = c_idapo(ina) - ilma = c_idalm(ina) - illa = c_idall(ina) - iout = 0 - do ii=ipoa,ipoa+ilma-1 - ! if(abs(c_cc(ii)) eps) a=real(c_cc(ii)) - if(abs(aimag(c_cc(ii)))> eps) b=aimag(c_cc(ii)) - c_cc(ii)=a+(0.0_dp,1.0_dp)*b - ! - enddo - ! call dapac(ina) - return - end subroutine dacc_lean - - subroutine c_dafun(cf,ina,inc) implicit none ! **************************** @@ -4960,8 +7026,10 @@ subroutine c_dafun(cf,ina,inc) ! integer ina,inc,incc character(4) cf - ! - if((.not.C_STABLE_DA)) then + + + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -5001,8 +7069,16 @@ subroutine dafunt(cf,ina,inc) ! data abcs /'abcdefghijklmnopqrstuvwxyz'/ data abcc /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/ - ! - if((.not.C_STABLE_DA)) then + + ! if(c_nomax.eq.1.and.newtpsa) then + ! return + ! endif + + ! if(c_nomax==2.and.newtpsa) then + ! return + ! endif + + if((.not.c_stable_da)) then if(C_watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -5060,7 +7136,7 @@ subroutine dafunt(cf,ina,inc) if(abs(a0).eq.0) then if(check_da) then messagelost="a0.eq.0 for INV in dafun" - C_STABLE_DA=.false. + c_stable_da=.false. C_check_stable=.false. call c_dadal1(iscr) call c_dadal1(inon) @@ -5121,7 +7197,7 @@ subroutine dafunt(cf,ina,inc) if(abs(a0)>hyperbolic_aperture) then if(check_da) then messagelost="a0>hyperbolic_aperture for EXP in dafun" - C_STABLE_DA=.false. + c_stable_da=.false. C_check_stable=.false. call c_dadal1(iscr) call c_dadal1(inon) @@ -5189,6 +7265,8 @@ subroutine c_take(h,m,ht) integer i,m,h,ht,b1,b2,b3,no,nv integer,dimension(c_lnv)::j complex(dp) r + + if(.not.c_stable_da) return no = c_nomax @@ -5254,15 +7332,38 @@ subroutine c_daran(ina,cm,xran) ! integer i,illa,ilma,ina,inoa,inva,ipoa,ipause,mypauses real(dp) cm,xran - ! - if((.not.C_%STABLE_DA)) then + + if(newtpsa) then + +ipoa=c_idapo(ina) + + do i=ipoa,ipoa+c_nmmax-1 + if(cm.gt.zero) then + c_cc(i) = bran(xran)+(0.0_dp,1.0_dp)*bran(xran) + c_cc(i) =c_cc(i)/abs(c_cc(i)) + if(abs(c_cc(i)).gt.cm) c_cc(i) = zero + elseif(cm.lt.zero) then + c_cc(i) = int(1+10*bran(xran)) + if(abs(c_cc(i)).gt.-ten*cm) c_cc(i) = zero + else + line='ERROR IN ROUTINE DARAN' + ipause=mypauses(31,line) + call dadeb !(31,'ERR DARAN2',1) + endif + enddo + + return + +endif + + if((.not.c_stable_da)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif return endif call dainf(ina,inoa,inva,ipoa,ilma,illa) - if((.not.C_%STABLE_DA)) then + if((.not.c_stable_da)) then if(c_%watch_user) then write(6,*) "big problem in dabnew ", sqrt(crash) endif @@ -5320,7 +7421,9 @@ subroutine c_dapek0(v,x,jj) integer,dimension(:)::v complex(dp),dimension(:)::x if(.not.c_stable_da) return + + do i=1,c_lnv jd(i)=0 enddo @@ -5337,7 +7440,8 @@ subroutine c_dapok0(v,x,jj) integer,dimension(:)::v complex(dp),dimension(:)::x if(.not.c_stable_da) return - + + do i=1,c_lnv jd(i)=0 enddo @@ -5354,6 +7458,8 @@ subroutine c_etcom(x,y,h,nd2) integer,dimension(:)::h,x,y integer,dimension(c_lnv)::t3 if(.not.c_stable_da) return + + call c_etall1(t1) call c_etall1(t2)