Reduce stack usage for better compatibility with macOS

General cleanup of the code, to reduce array temporaries creation by gfortran
12.2. These changes should lead to lower stack usage, particularly important
on macOS and other systems with low stack limits. The effect on performance
on other systems should be neutral, or even positive.

Additionally, all unsed wrapper routines were removed from lapack.f90 glue
interface to LAPACK. The remaining routines (diagonalization of general
complex matrices) were rewritten to use the recommended call sequence for
allocating temporary storage.

There are still several crashes on macOS in the test set, but the present state
is a marked improvement.

Author: Serguei Patchkovskii

Makefile

@@ -19,12 +19,12 @@ ACT2 = -e 's/^!\*nm/    /' # Disable MPI statements
 #
 # System-specific overrides
 #
-  include vanilla.mak
+# include vanilla.mak
 # include configs/babel-gfortran_opt.mak
 # include configs/babel-ifort18_opt.mak
 # include configs/zen-gfortran-7_opt.mak
 # include configs/zen-gfortran-11_opt.mak
-# include configs/zen-gfortran-12_opt.mak
+  include configs/zen-gfortran-12_opt.mak
 # include configs/zen-oneapi_opt.mak
 # include configs/zen-oneapi_opt_mpi.mak
 # include configs/zen-aocc-1.1_opt.mak      # VERY SLOW CODE. DO NOT USE.
@@ -59,12 +59,6 @@ LIBS = $(LAPACK) $(LAPACK) $(LIBEXTRA)
 	$(ACT) $(ACT2) $< >preprocess/$<
 	$(F90) -c preprocess/$<
 
-dgefa.o:        dgefa.f
-	$(F90) -c dgefa.f
-
-dgedi.o:        dgedi.f
-	$(F90) -c dgedi.f
-
 #hacks.o:	hacks.f90 accuracy.o
 #	$(F90) -O0 -c hacks.f90
 
@@ -86,8 +80,6 @@ LIBSPHERICAL += composition_analysis.o
 LIBSPHERICAL += constants.o
 LIBSPHERICAL += coulomb_functions.o
 LIBSPHERICAL += cubic_spline.o
-LIBSPHERICAL += dgefa.o
-LIBSPHERICAL += dgedi.o
 LIBSPHERICAL += hacks.o
 LIBSPHERICAL += lapack.o
 LIBSPHERICAL += math.o

checkpoint_tools.f90

@@ -61,7 +61,7 @@ module checkpoint_tools
   integer(ik), save         :: ckpt_max_checkpoints = 3_ik      ! Maximum number of checkpoints to keep.
   integer(ik), save         :: ckpt_interval        = 10000_ik  ! Number of timesteps between checkpoints.
   !
-  character(len=clen), save :: rcsid_checkpoint_tools = "$Id: checkpoint_tools.f90,v 1.16 2021/04/26 15:44:44 ps Exp ps $"
+  character(len=clen), save :: rcsid_checkpoint_tools = "$Id: checkpoint_tools.f90,v 1.17 2022/10/08 17:24:26 ps Exp ps $"
   !
   type ckpt_data
     private                                                         ! This data is not for exterrrnal consumption
@@ -243,6 +243,9 @@ subroutine ckpt_load_checkpoint(ck,wfn_l,wfn_r,tsurf)
       action = 'verifying header'
       if (tag/='SCID CKPT' .or. version/=ckpt_version .or. any(options.neqv.options2)) then
         write (out,"('Checkpoint header mismatch')")
+        ! gfortran introduces an array temporary for the transfer() below. There
+        ! does not seem to be anything we could do about it: adding an explicit
+        ! temporary still creates a temporary array ...
         write (out,"('Expected: ',a,' V.',i0,' options = ',2l1,'(',2z8,')')") &
                'SCID CKPT', ckpt_version, options, transfer(options,1_ik,2)
         write (out,"('Received: ',a,' V.',i0,' options = ',2l1,'(',2z8,')')") &

composition_analysis.f90

@@ -37,7 +37,7 @@ module composition_analysis
   public ca_maxram
   public rcsid_composition_analysis
   !
-  character(len=clen) :: rcsid_composition_analysis = "$Id: composition_analysis.f90,v 1.17 2021/04/26 15:44:44 ps Exp ps $"
+  character(len=clen) :: rcsid_composition_analysis = "$Id: composition_analysis.f90,v 1.18 2022/10/08 17:24:26 ps Exp ps $"
   !
   real(rk), save :: ca_maxram = 0._rk ! Maximum amount of memory which can be used during the analysis step
                                       ! This limit does NOT include the memory needed to compute atomic 
@@ -167,6 +167,8 @@ subroutine ca_analyze(verbose,threshold,wfn_l,wfn_r,tsurf)
       call wt_atomic_solutions(verbose,lval,block_eval,block_evec)
       en(:,lval) = block_eval(:)
       scan_spin: do ispin=1,sd_nspin
+        ! gfortran is having trouble generating good code for this array assignment
+        ! there does not seem to be anything we could do about it though ...
         amp(:,1,ispin,mmin:mmax,lval) = matmul(transpose(block_evec(:,:,1)),wfn_r%wfn(:,ispin,lval,mmin:mmax))
         amp(:,2,ispin,mmin:mmax,lval) = matmul(transpose(block_evec(:,:,2)),wfn_l%wfn(:,ispin,lval,mmin:mmax))
       end do scan_spin

configs/macos_m1-gfortran_opt.mak

@@ -1,14 +1,16 @@
-BUILD_ID :="Optimized gfortran (-O1), built on $(shell hostname) at $(shell date)"
-ACT = sed -e 's/^!\*nq/    /' # Enable quad-math statements
-# WARNING: gfortran 11.2.0 appears to have many bugs in OpenMP code generation, especially at -O3
-# WARNING: Using -O1 reduces (but does not completely eliminate) these bugs.
-# WARNING: Please make sure that the ENTIRE test suite runs correctly before using in production.
+BUILD_ID :="Optimized gfortran, built on $(shell hostname) at $(shell date)"
+ACT = sed -e 's/^!\*nq/    /' # Disable quad-math statements
+# -Ofast is known to produce code which crashes on M1, at least with gfortran 12.2
 F90 = gfortran -I. \
-      -O1 -march=native -mtune=native -fopenmp \
-      -floop-block \
+      -flto -O1 -fprotect-parens -march=native -mtune=native -fopenmp \
       -ffast-math -fcx-fortran-rules \
+      -floop-block \
       -fno-realloc-lhs -fbacktrace -g \
       -cpp -D__BUILD_ID__='$(BUILD_ID)' -ffree-line-length-none
 F90L = $(F90) 
-LAPACK = -L/opt/homebrew/opt/openblas/lib -lopenblas
+#
+#  Both openblas and Accelerate are known to work, with openblas being somewhat faster
+#
+#LAPACK = -L/opt/homebrew/opt/openblas/lib -lopenblas
+LAPACK = -framework Accelerate
 LIBEXTRA = 

coulomb_functions.f90

@@ -72,7 +72,7 @@ module coulomb_functions
    public coulombFG, coulombF, coulombBound
    public rcsid_coulomb_functions
    !
-   character(len=clen), save :: rcsid_coulomb_functions = "$Id: coulomb_functions.f90,v 1.13 2021/04/26 15:44:44 ps Exp ps $"
+   character(len=clen), save :: rcsid_coulomb_functions = "$Id: coulomb_functions.f90,v 1.14 2022/10/08 17:24:26 ps Exp ps $"
    !
  ! integer, parameter      :: out         = 6
  ! integer, parameter      :: ik          = selected_int_kind(15)
@@ -142,7 +142,7 @@ subroutine coulombFG(lambdamin,nlambda,x,eta,fg,skip_g)
      !
      if (x<=0) stop 'coulomb_functions%coulombFG - bad x'
      if (nlambda<=0) stop 'coulomb_functions%coulombFG - bad nlambda'
-     if (any(ubound(fg)/=(/nlambda,4_ik/))) stop 'coulomb_functions%coulombFG - bad array sizes'
+     if (ubound(fg,1)/=nlambda .or. ubound(fg,2)/=4_ik) stop 'coulomb_functions%coulombFG - bad array sizes'
      !
      !  Start by filling coefficient tables: we'll need these repeatedly.
      !