diff --git a/src/accuracy_constants.f b/src/accuracy_constants.f index c250ccf..7b90572 100644 --- a/src/accuracy_constants.f +++ b/src/accuracy_constants.f @@ -9,4 +9,5 @@ ! integer, parameter :: iacc = int16 !int*2 integer, parameter :: iacc = int32 !int*4 ! integer, parameter :: iacc = int64 !int*8 + real(dp),parameter:: pi = acos(-1.0_dp) end module diff --git a/src/data_module.f b/src/data_module.f index 9713923..e99ad53 100644 --- a/src/data_module.f +++ b/src/data_module.f @@ -8,14 +8,14 @@ double precision,protected, dimension(:,:), allocatable :: y_m double precision,protected, dimension(:,:), allocatable :: wt_m double precision,protected, dimension(:,:), allocatable :: ny_m - + double precision,protected, dimension(:,:,:), allocatable:: U_m contains !------------------------------ - subroutine init_data(numdatpt,q,x1,x2,y,wt,ny) + subroutine init_data(numdatpt,q,x1,x2,y,wt,ny,U) - use dim_parameter, only: qn, ntot + use dim_parameter, only: qn, ntot,ndiab implicit none @@ -26,6 +26,7 @@ double precision y(ntot,*) double precision wt(ntot,*) double precision ny(ntot,*) + double precision U(ndiab,ndiab,*) allocate(q_m(qn,numdatpt)) allocate(x1_m(qn,numdatpt)) @@ -33,6 +34,7 @@ allocate(y_m(ntot,numdatpt)) allocate(wt_m(ntot,numdatpt)) allocate(ny_m(ntot,numdatpt)) + allocate(U_m(ndiab,ndiab,numdatpt)) do i=1,numdatpt q_m(1:qn,i)=q(1:qn,i) @@ -41,6 +43,7 @@ y_m(1:ntot,i)=y(1:ntot,i) wt_m(1:ntot,i)=wt(1:ntot,i) ny_m(1:ntot,i)=ny(1:ntot,i) + U_m(1:ndiab,1:ndiab,i)=U(1:ndiab,1:ndiab,i) enddo end subroutine @@ -48,7 +51,7 @@ !------------------------------ subroutine dealloc_data() - deallocate(q_m,x1_m,x2_m,y_m,wt_m,ny_m) + deallocate(q_m,x1_m,x2_m,y_m,wt_m,ny_m,U_m) end subroutine end module data_module diff --git a/src/funcs.f b/src/funcs.f index 39583eb..e46d1b1 100644 --- a/src/funcs.f +++ b/src/funcs.f @@ -102,8 +102,8 @@ integer i, npar double precision p(npar), dp(npar) ! double precision, parameter :: d = 1.d-4 - double precision, parameter :: d = 1.d-6 !Standard -! double precision, parameter :: d = 1.d-8 +! double precision, parameter :: d = 1.d-6 !Standard + double precision, parameter :: d = 1.d-8 double precision, parameter :: thr = 1.d-12 do i=1,npar dp(i)=abs(p(i)*d) diff --git a/src/genetic.f b/src/genetic.f index 12838bb..3c6ef24 100644 --- a/src/genetic.f +++ b/src/genetic.f @@ -1,6 +1,6 @@ program genetic ! module for dimensioning parameters - use dim_parameter,only: qn,ntot,numdatpt,dealloc_dim + use dim_parameter,only: qn,ntot,numdatpt,dealloc_dim,ndiab !data module use data_module, only: init_data, dealloc_data ! parser module @@ -51,6 +51,8 @@ double precision rms,old character(len=80) filename character(len=80) chkpnt +! JP + double precision,allocatable :: U_pes(:,:,:) ! AtD matrix from pes ! ----------------------------- @@ -73,6 +75,7 @@ allocate(par(npar,nset),prange(2,npar)) allocate(x1_in(qn,numdatpt),x2_in(qn,numdatpt)) + allocate(U_pes(ndiab,ndiab,numdatpt)) ! JP allocate U call rinit(p,prange,p_spread,p_act,npar) par=0.d0 @@ -81,9 +84,10 @@ !------------------------------------------------- - call data_transform(q_in,x1_in,x2_in,y_in,wt_in,p,npar,p_act) + call data_transform(q_in,x1_in,x2_in,y_in,U_pes, + > wt_in,p,npar,p_act) !Fabian: Read data into module - call init_data(numdatpt,q_in,x1_in,x2_in,y_in,wt_in,y_in) + call init_data(numdatpt,q_in,x1_in,x2_in,y_in,wt_in,y_in,U_pes) !------------------------------------------------- #ifdef mpi_version @@ -147,7 +151,7 @@ #endif deallocate(q_in,x1_in,x2_in,y_in,wt_in, - $ p,par,p_act,p_spread,prange) + $ p,par,p_act,p_spread,prange,U_pes) call dealloc_data call dealloc_dim ! call dealloc_dw_ptr diff --git a/src/marq.f b/src/marq.f index cf3bf1a..027d7a8 100644 --- a/src/marq.f +++ b/src/marq.f @@ -94,7 +94,8 @@ !> Initialize skip, Marquardt parameter, error variables and parameter work arrays skip=.false. - ilamda=0.1d0 !Initial Marquardt parameter + !ilamda=0.1d0 !Initial Marquardt parameter + ilamda=0.05d0 !Initial Marquardt parameter alamda=ilamda !Current Marquardt parameter rms=1.d6 atry(1:npar)=par(1:npar) @@ -309,6 +310,7 @@ c############################################################### ! Internal variables integer i,j,k,l,m,n !< iteration variables integer nloop + integer active(mfit) !------------------------------------------------------------------------- @@ -323,11 +325,18 @@ c############################################################### !if(.not.hybrid) nloop=nstat ! commented out by JP on 08.10.2025 ! In fitting dipole nstat is meaningless, I always go with ntot ! and send to zero the unused parts of y_m and ymod + j=0 + do i=1,npar + if(ma(i).ge.1) then + j=j+1 + active(j)=i + endif + enddo do i=1,numdatpt - + call funcs(i,par,ymod,dyda,npar,ma,skip) - write(58,*) "ymod ",i,nloop,ymod(1:ntot) + !write(58,*) "ymod ",i,nloop,ymod(1:ntot) if (skip) return !Idea: Since the quantities dyda,dy and wt_m are rather small, one might consider scaling them @@ -335,41 +344,26 @@ c############################################################### !Idea: Scale dyda,dy and wt_m by 1D+5; final rescale of alpha,beta and chisq by 1D-10 do n=1,nloop - dy(n)=y_m(n,i)-ymod(n) - !write(58,*) "dy ",n,i,ymod(n),y_m(n,i),dy(n) - j=0 - do l=1,npar - ! Nicole: values of ma (active parameter) changed - if (ma(l).ge.1) then - j=j+1 - k=0 - do m=1,l -! Nicole: values of ma (active parameter) changed - if (ma(m).ge.1) then - k=k+1 - !(wt*J)^T*(wt*J) - alpha(j,k)=alpha(j,k)+ - $ (dyda(n,l)*dyda(n,m))*(wt_m(n,i)*wt_m(n,i)) - endif - enddo - !(wt*J)^T*(wt*delta_y) - beta(j)=beta(j)+ - $ (dy(n)*dyda(n,l))*(wt_m(n,i)*wt_m(n,i)) - endif - enddo - !(wt*delta_y)*(wt*delta_y) - chisq=chisq+ - $ (dy(n)*dy(n))*(wt_m(n,i)*wt_m(n,i)) - enddo - - enddo + dy(n)=(y_m(n,i)-ymod(n))*wt_m(n,i) + dyda(n,:) = dyda(n,:)*wt_m(n,i) + !write(58,"(A3,2I0,3ES18.9)")"dy ",n,i,ymod(n),y_m(n,i),dy(n) + do l=1,mfit + do m=1,l + alpha(m,l)=alpha(m,l)+ + & (dyda(n,active(l))*dyda(n,active(m)))!*wt_m(n,i) !JP + enddo + beta(l) = beta(l) + (dy(n) * dyda(n,active(l))) + enddo + chisq = chisq + (dy(n)**2) + enddo + enddo !------------------------------------------------------------------------- !Fill in missing parts of the symmetric matrix alpha do i=2,mfit do j=1,i-1 - alpha(j,i)=alpha(i,j) + alpha(i,j)=alpha(j,i) enddo enddo diff --git a/src/parser/io_parameters.f b/src/parser/io_parameters.f index 57e95f8..1d46e14 100644 --- a/src/parser/io_parameters.f +++ b/src/parser/io_parameters.f @@ -25,10 +25,12 @@ character(len=64) :: errcat(maxerrors) !< list of generic error Messages defined in errcat.incl ! parameter key declaration - integer, parameter :: maxpar_keys=400 ! key,maxpar_keys,parkeynum,parkeylen,ec_read,ec_dim + > key,maxpar_keys,parkeynum,parkeylen,ec_read,ec_dim, + > input_key_order use dim_parameter,only: pst,max_par use keys_mod, only: init_keys implicit none @@ -47,7 +48,6 @@ !Fabian: Include user specific keys call init_keys ! include 'keys.incl' - !Fabian: No need that this is within keys.incl since it is generic do j=1,maxpar_keys if (key(1,j)(1:1).eq.' ') then @@ -72,6 +72,9 @@ call write_oneline(fmt,std_out) endif +! Prepare empty list for tracking key appearances + allocate(input_key_order(0)) + ! reading cards for the number of parameters' @@ -87,9 +90,18 @@ if (infile(i)(1:key_end).eq.key(ktype,j)) then if(dbg) write(6,*) key(ktype,j),' read' read (infile(i)(key_end+1:llen),*) pst(2,j) + ! track of order of appearence of keys in the input file + input_key_order = [input_key_order,j] ! JP + endif enddo enddo + !Append all keys that are defined (in keys.f90) but are missing in the input file + do i=1, parkeynum + if (.not. any(input_key_order.eq.i)) then + input_key_order=[input_key_order,i] + endif + enddo !.. compute total number of parameters: do i=1, parkeynum