Merge pull request #923 from MethodicalAcceleratorDesign/ld-dev

fix warnings in plot module for non-twiss tables with haxis=s.

Makefile_f90

@@ -64,6 +64,8 @@ ifneq ($(FCNAME),nagfor)
 Sg_sagan_wiggler.o: FFLAGS += $(if $(call eq,$(FCNAME),lf95),-dollar-ok,-fdollar-ok)
 endif
 
+# BOZ flag set at the end of this file, i.e. need gfortran version
+# util.o: FFLAGS += $(if $(call search,$(FCVER1),$(FCBOZV)),-fallow-invalid-boz,)
 
 #######################
 # Fortran dependencies (case not automatic!)
@@ -179,19 +181,26 @@ endif
 
 ifneq ($(filter $(BUILDGOALS),$(MAKECMDGOALS)),)
 ifeq ($(FCNAME),gfortran)
-FCVER  := $(strip $(shell $(FC) -dumpversion))
-ifeq ($(findstring GNU Fortran (GCC),$(FCVER)),GNU Fortran (GCC))
-FCVER  := $(word 4,$(FCVER)) # gfortran -dumpversion bug
-endif
-FCVERL := $(subst ., ,$(FCVER))
-FCVER1 := $(findstring $(word 1,$(FCVERL)),"456789")
-FCVER2 := 0
-  ifeq ($(FCVER1),4)
-  FCVER2 := $(findstring $(word 2,$(FCVERL)),"456789")
+  FCVER  := $(strip $(shell $(FC) -dumpversion))
+  ifeq ($(findstring GNU Fortran (GCC),$(FCVER)),GNU Fortran (GCC))
+    FCVER := $(word 4,$(FCVER)) # gfortran -dumpversion bug
+  endif
+  FCVERL := $(subst ., ,$(FCVER))
+  FCVER1 := $(word 1,$(FCVERL))
+  FCVER2 := $(word 2,$(FCVERL))
+  FCVER3 := $(word 3,$(FCVERL))
+  FCVERS := 4 5 6 7 8 9 10 11 12 13 14 15 16 17 19 20
+  FCISOK := $(call search,$(FCVER1),$(FCVERS))
+  ifeq ($(FCISOK),4)
+    FCISOK := $(call search,$(FCVER2),$(FCVERS))
   endif
-ifeq ($(and $(FCVER1),$(FCVER2)),)
-$(error MAD-X requires gfortran >= 4.4 ($(if $(FCVER),$(FCVER),none) detected))
-endif # version
+  ifeq ($(FCISOK),)
+    $(error MAD-X requires gfortran >= 4.4 ($(if $(FCVER),$(FCVER),none) detected))
+  endif # version
+
+  FCBOZV := 10 11 12 13 14 15 16 17 19 20
+  util.o: FFLAGS += $(if $(call search,$(FCVER1),$(FCBOZV)),-fallow-invalid-boz,)
+
 endif # gfortran
 endif # cmdgoal
 

libs/DISTlib/source/distinterface.h

@@ -31,4 +31,4 @@ void free_distribution(int i);
 void settasmatrix_element(double value, int row, int column);
 void addclosedorbit(double *clo);
 void getarraylength(int *totlength);
-void getrefpara(double *energy0, double *mass0, int *a0, int *z0);
+void getrefpara(double *energy0, double *mass0, int *a0, int *z0);

make/make.lib

@@ -18,7 +18,7 @@
 
 # some useful variables
 EMPTY :=
-SPACE := 
+SPACE :=
 SPACE +=
 COMMA := ,
 
@@ -32,14 +32,11 @@ first  = $(word 1,$1)
 last   = $(word $(words $1),$1)
 rest   = $(wordlist 2,$(words $1),$1)
 chop   = $(wordlist 2,$(words $1),_ $1)
-uniq   = $(strip $(if $1,$(call uniq,$(call chop,$1)) \
-            $(if $(filter $(call last,$1),$(call chop,$1)),,$(call last,$1))))
-trans  = $(if $1,$(call trans,$(call rest,$1),$(call rest,$2)\
-                    ,$(patsubst $(firstword $1),$(firstword $2),$3)),$3) 
-find   = $(if $(and $1,$2),$(if $(findstring $1,$(firstword $2)),$2,\
-              $(call find,$1,$(call rest,$2))),)
-trunc  = $(if $(and $1,$2),$(if $(findstring $1,$(lastword $2)),$2,\
-              $(call trunc,$1,$(call chop,$2))),)
+uniq   = $(strip $(if $1,$(call uniq,$(call chop,$1)) $(if $(filter $(call last,$1),$(call chop,$1)),,$(call last,$1))))
+trans  = $(if $1,$(call trans,$(call rest,$1),$(call rest,$2),$(patsubst $(firstword $1),$(firstword $2),$3)),$3)
+find   = $(if $(and $1,$2),$(if $(findstring $1,$(firstword $2)),$2,$(call find,$1,$(call rest,$2))),)
+trunc  = $(if $(and $1,$2),$(if $(findstring $1,$(lastword $2)),$2,$(call trunc,$1,$(call chop,$2))),)
+search = $(if $(and $1,$2),$(if $(call eq,$1,$(firstword $2)),$1,$(call search,$1,$(call rest,$2))),)
 
 assert = $(if $1,,$(error $2))
 exists = $(if $(wildcard $1),,$(error $2))

src/mad_dict.c

@@ -309,6 +309,7 @@ const char *const_command_def =
 "update_from_parent = [l, false, true], "
 "keep_exp_move = [l, false, true], "
 "thin_cf = [l, false, true], "
+"fdstep = [r, 0], " /* ld 07.2020, finite difference step for JACOBIAN and LMDIF */
   /* BB and SPCH related options */
 "bborbit          = [l, false, true], " /* frs */
 "bb_ultra_relati  = [l, false, true], " /* frs */

src/mad_dist.h

@@ -5,4 +5,4 @@
 void pro_distribution(struct in_cmd* p);
 void setdistparameters(char *type, int cut_l, double* cuts, int index, int *dist_type, double *start_v, double *stop_v, double *coord_t);
 
-#endif // MAD_DIST_H
+#endif // MAD_DIST_H

src/mad_gvar.h

@@ -319,8 +319,8 @@ extern char
            tag_type[MAX_TAG][16],
            tag_code[MAX_TAG][16];
 
-extern time_t     last_time,
-           start_time;
+extern time_t last_time,
+              start_time;
 
 extern char filenames[100][500];
 extern int  currentline[100];

src/mad_plot.c

@@ -122,8 +122,8 @@ plot_option(char* name)
 void
 exec_plot(struct in_cmd* cmd)
 {
-  int i, j, k, ierr, notitle = strcmp(title,"no-title") == 0 ? 1 : 0;
-  int nointerp = 0, multiple = 0, noversion = 0, nolegend = 0, s_haxis = 1, track_flag = 0;
+  int i, j, k, ierr, notitle = !strcmp(title,"no-title");
+  int nointerp = 0, multiple = 0, noversion = 0, nolegend = 0, s_haxis = 0, track_flag = 0;
   int tsm1 = TITLE_SIZE - 1, tsm2 = TITLE_SIZE - 2;
   int part_idx[100], curr, track_cols_length, haxis_idx = 0, vaxis_idx = 0;
   int size_plot_title = tsm1, size_version = tsm1;
@@ -162,7 +162,7 @@ exec_plot(struct in_cmd* cmd)
   /* get haxis_name & s_haxis flag */
   haxis_name = command_par_string_user("haxis", this_cmd->clone);
   if (haxis_name) {
-    s_haxis = strcmp(haxis_name,"s");
+    s_haxis = !strcmp(haxis_name,"s");
   }
 
   /* get table_name & track_flag */
@@ -321,7 +321,7 @@ exec_plot(struct in_cmd* cmd)
   else { /* normal plot */
 
     /* if haxis is "s" and no interpolation, check if table name is the same of the last twiss call */
-    if (nointerp == 0 && s_haxis == 0) {
+    if (nointerp == 0 && s_haxis == 1) {
 
       if (!current_sequ->tw_table) {
 	warning("PLOT - no TWISS table present", "PLOT command ignored");

src/match.f90

@@ -709,17 +709,23 @@ subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
 !----------------------------------------------------------------------*
       integer i,iflag,j,ldfjac,m,n
       double precision eps,epsfcn,fjac(ldfjac,n),fvec(m),h,temp,wa(m),  &
-     &x(n),zero,epsmch
+     &x(n),zero,epsmch,fdstep
       parameter(zero=0d0,epsmch=1d-16)
       external fcn
+      integer, external :: get_option
 
+      fdstep = get_option('fdstep ')
       eps = sqrt(max(epsfcn,epsmch))
       iflag = 0
 
       do j = 1, n
         temp = x(j)
-        h = eps*abs(temp)
-        if (h .eq. zero) h = eps
+        if (fdstep .ne. zero) then
+          h = fdstep
+        else
+          h = eps*abs(temp)
+          if (h .eq. zero) h = eps
+        endif
         x(j) = temp + h
         call fcn(m,n,x,wa,iflag)
         x(j) = temp

src/matchjc.f90

@@ -98,7 +98,7 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
 !----------------------------------------------------------------------*
 
       integer izero
-      integer i,iflag,info,j,level,m,n,calls,call_lim
+      integer i,iflag,info,j,level,m,n,calls,call_lim,ecalls
       integer ireset,strategy,bisec,match_mode
       double precision fmin_old,epsfcn
       double precision ftol,gtol
@@ -121,7 +121,7 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
       integer, external :: get_option
 
       debug = get_option('debug ')
-
+      ecalls = 0
       ireset = 0
       izero = 0
       info = 0
@@ -188,11 +188,24 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
 !---- Start main loop
   20  continue
 
+      if (debug .ne. 0) then
+        write(*,*) "Current variables values"
+        do i=1,N
+          write(*,'(i5,e20.12)') i, x(i)
+        enddo
+      endif
+
 !---- Reset ireset
       ireset=0
 
 !---- Calculate the jacobian
       call fdjac2(fcn,m,n,x,fvec,fjac,m,iflag,xtol,wa4)
+      ecalls=ecalls+n
+      if (iflag .ne. 0) then
+        call fort_warn('JACOBIAN', ' stopped, possibly unstable')
+        info = - 1
+        go to 300
+      endif
 
       if (strategy.eq.2) then
 !        call DGESDD('A',M,N,fjac,M,SV,U,M,VT,N,                         &
@@ -246,6 +259,8 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
         endif
       end do
 
+!-- BUG, effjac and effsol can be empty in case jacobian is e.g. full of zero...
+
  33   continue
 
 
@@ -366,7 +381,7 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
         enddo
 
         if (ireset.eq.20) then
-          write(*,*) "Too loops in system resizing, set strategy=1"
+          write(*,*) "Too many loops in system resizing, set strategy=1"
           strategy=1
 !----     Exit var cancel loop
           goto 34
@@ -376,7 +391,7 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
 
         if ((effvar.lt.effcon+1).or.(condnum.gt.1E10)) then
 !----     Impossible to get  better
-          write(*,*) "Too var to exclude, set strategy=1"
+          write(*,*) "Too many variables to exclude, set strategy=1"
           strategy=1
 !----     Exit var cancel loop
           goto 34
@@ -398,6 +413,12 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
 
 !----- Calculate the penalty function
       call FCN(M,N,x,fvec,IFLAG)
+      ecalls=ecalls+1
+      if (iflag .ne. 0) then
+        call fort_warn('JACOBIAN', ' stopped, possibly unstable')
+        info = - 1
+        go to 300
+      endif
 
       fmin_old=fmin
       fmin = vdot(m, fvec,fvec)
@@ -414,12 +435,18 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
         ! go back to average solution
         do i=1,N
           x(i)=(x(i)+xold(i))*.5
-          dx(i)=(x(i)-xold(i))*.5
+          dx(i)=(x(i)-xold(i))*.5                   ! LD: never used
         enddo
         j=j+1
         call FCN(M,N,x,fvec,IFLAG)
+        ecalls=ecalls+1
+        if (iflag .ne. 0) then
+          call fort_warn('JACOBIAN', ' stopped, possibly unstable')
+          info = - 1
+          go to 300
+        endif
         fmin = vdot(m,fvec,fvec)
-        dxnorm=sqrt(vdot(N, dx, dx)/vdot(N, x, x))
+        dxnorm=sqrt(vdot(N, dx, dx)/vdot(N, x, x))  ! LD: never used
         if (fmin.gt.fmin_old2) then
           goto 37
         endif
@@ -436,9 +463,13 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
 !---- calculate the target function and set new values
       call FCN(M,N,x,fvec,IFLAG)
       calls=calls+1
+      if (iflag .ne. 0) then
+        call fort_warn('JACOBIAN', ' stopped, possibly unstable')
+        info = - 1
+        go to 300
+      endif
 
       fmin = vdot(m,fvec,fvec)
-
       xnorm=DNRM2(N, x, 1)
       dxnorm=DNRM2(effvar, effsol, 1)/xnorm
       write(*,"('call: ',I5,' Dx = ', e16.8,                            &
@@ -464,7 +495,7 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
         goto 300
       endif
       ! Iteration do not converge
-      ! For the first calls it can happen due to not fisical penalty
+      ! For the first calls it can happen due to not physical penalty
       if((calls.ge.3).and.(abs(1-fmin/fmin_old).lt.1D-15) ) then
 !        print *, 'dbg',fmin,fmin_old,abs(1-fmin/fmin_old)
         print *, '++++++++++ JACOBIAN ended: minimum found'
@@ -497,12 +528,21 @@ subroutine jacob(fcn,m,n,strategy,calls,call_lim,                 &
       enddo
       ! set values to the previous solution
       call FCN(M,N,X,fvec,IFLAG)
+      ecalls=ecalls+1
+      if (iflag .ne. 0) then
+        call fort_warn('JACOBIAN', ' stopped, possibly unstable')
+        info = - 1
+      endif
 
 !---- Store the final distance
       do i = 1, n
         dx(i)=(x(i)-xstart(i))
       enddo
 
+      if (debug .ne. 0) then
+        write(*,*) 'Effective number of calls:', calls+ecalls
+      endif
+
 !      xnorm=DNRM2(N, x, 1)
       dxnorm=sqrt(DNRM2(N, dx, 1))
       write(*,*) 'Final difference norm:',dxnorm

src/plot.f90

@@ -310,7 +310,7 @@ subroutine pefill(ierr)
   !----------------------------------------------------------------------*
   integer, intent(OUT) :: ierr
 
-  integer :: i, j, k, l, nint, new1, new2, currtyp, crow, pltyp, pos_flag
+  integer :: i, j, k, l, nint, new1, new2, currtyp, crow, pltyp, pos_flag, warn
   integer :: ilist(mtype), p(mxplot), proc_n(2, mxcurv)
   double precision :: fact, d_val, d_val1
   double precision :: currpos, currleng, currtilt
@@ -335,6 +335,9 @@ subroutine pefill(ierr)
   !--- Initialize return code
   ierr = 0
 
+  !-- Retrieve warn value
+  warn = get_option('warn ')
+
   !--- Initialize marker_plot logical
   marker_plot = get_value('plot ','marker_plot ') .ne. zero
   range_plot  = get_value('plot ','range_plot ') .ne. zero
@@ -460,12 +463,14 @@ subroutine pefill(ierr)
         ! for codes 1 to 8 only and elements with length
         if (currleng.gt.0.d0 .and. currtyp.gt.1 .and. currtyp.lt.8) then
            currtilt = node_value('tilt ')
+           call set_option('warn ', 0) !-- LD: hugly fix as this code always expect a twiss table...
            k = double_from_table_row(tabname, 'k1l ' , j, currk1l);  currk1l  = currk1l /currleng
            k = double_from_table_row(tabname, 'k1sl ', j, currk1sl); currk1sl = currk1sl/currleng
            k = double_from_table_row(tabname, 'k2l ' , j, currk2l);  currk2l  = currk2l /currleng
            k = double_from_table_row(tabname, 'k2sl ', j, currk2sl); currk2sl = currk2sl/currleng
            k = double_from_table_row(tabname, 'k3l ' , j, currk3l);  currk3l  = currk3l /currleng
            k = double_from_table_row(tabname, 'k3sl ', j, currk3sl); currk3sl = currk3sl/currleng
+           call set_option('warn ', warn)
         endif
 
         !--- sbend, rbend