Commit c2a95e94 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

switch to fortran90 version

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@95 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent 3c41cd2a
This diff is collapsed.
# #################################################################
# Makefile for MillePede II with possible input from C
# Makefile for MillePede II (Fortran90) with possible input from C
#
# Works on 64-bit SL5 with gcc 434. See comments about different
# gcc versions inline to get a hint about the necessary adjustments.
# Works on 64-bit SL5 with gcc version 4.4.4 (or higher).
# See comments about different gcc versions inline to get a
# hint about the necessary adjustments.
# #################################################################
# On 32-bit systems:
# LARGE_MEMORY_OPT =
# LARGE_SIZE=4
# #################################################################
# If compiler doesn't understand INTEGER(KIND=) (FORTRAN95) use
# LARGE_SIZE=4 and replace INTEGER(KIND=) by INTEGER*4 or
# LARGE_SIZE=8 and replace INTEGER(KIND=) by INTEGER*8
# e.g. change source files with sed:
# for name in pede.F dynal.F dynal.inc mpinds.inc; do sed 's/(KIND=LARGE)/*4/' $name >! tmp; mv tmp $name; done
# #################################################################
#
# compiler options
#
# To make it link with a BIG (>2 GB) static array in dynal.inc,
# we need '-mcmodel=medium', but does not work on 32 bit machines.
# (Seems not to be needed for recent gcc like gcc461,
# but until at least gcc434 it is needed.)
LARGE_MEMORY_OPT=-mcmodel=medium
# All but 'yes' disables support of reading C-binaries:
SUPPORT_READ_C = yes
# If yes (and if SUPPORT_READ_C is yes), uses rfio, i.e. shift library, for IO,
# requires shift library to be installed
# (valid only for default executables, pede_*rfio*-ones will always get it):
SUPPORT_C_RFIO =
# yes
# If yes (and if SUPPORT_READ_C is yes), use zlib to read gzipped binary files
SUPPORT_ZLIB = yes
#
# Define the number of words for matrix/vector storage for the default pede program:
NUMBER_OF_WORDS = 100000000 # =400 MB
# Define the size of LARGE integers (4: INTERGER*4, 8: INTEGER*8)
LARGE_SIZE=8
# ompP profiler (http://www.ompp-tool.com)
OMPP =
#kinst-ompp
OMPP = kinst-ompp
#
FCOMP = $(OMPP) gcc
F_FLAGS = -Wall -fautomatic -fno-backslash -O3 ${LARGE_MEMORY_OPT}
FCOMP = $(OMPP) gcc44
F_FLAGS = -Wall -fautomatic -fno-backslash -O3 -cpp -DLARGE_SIZE=${LARGE_SIZE}
#
CCOMP = $(OMPP) gcc
C_FLAGS = -Wall -O3 -Df2cFortran ${LARGE_MEMORY_OPT}
CCOMP = $(OMPP) gcc44
C_FLAGS = -Wall -O3 -Df2cFortran
C_INCLUDEDIRS = # e.g. -I .
# gcc3: C_LIBS = -lg2c -lfrtbegin
# gcc4 till gcc44:
C_LIBS = -lgfortran -lgfortranbegin
# gcc45, gcc46: C_LIBS = -lgfortran -lm
DEBUG = # e.g. -g
#
# Multithreading with OpenMP (TM) (needs gcc41 or above)
# Multithreading with OpenMP (TM)
C_LIBS += -lgomp
F_FLAGS += -fopenmp
#
LOADER = $(OMPP) gcc
L_FLAGS = -Wall -O3 ${LARGE_MEMORY_OPT}
LOADER = $(OMPP) gcc44
L_FLAGS = -Wall -O3
#
# objects for this project
#
USER_OBJ_PEDE_STATIC = mptest.o mptst2.o mille.o mpnum.o mptext.o mphistab.o \
minresblas.o minres.o randoms.o vertpr.o linesrch.o bandmatrix/Dbandmatrix.o
USER_OBJ_PEDE = ${USER_OBJ_PEDE_STATIC} pede.o dynal.o
USER_OBJ_PEDE = mpdef.o mpdalc.o mpmod.o mpbits.o mptest1.o mptest2.o mille.o mpnum.o mptext.o mphistab.o \
minresblas.o minres.o randoms.o vertpr.o linesrch.o Dbandmatrix.o pede.o
#
# Chose flags/object files for C-binary support:
#
......@@ -80,32 +60,20 @@ ifeq ($(SUPPORT_READ_C),yes)
endif
endif
endif
#
F_FLAGS += -DLARGE_SIZE=${LARGE_SIZE}
#
#
# Make the executables
# The specific ones (part of hack, see below) + the single:
EXECUTABLES = pede pede_1GB pede_2GB
# Or also those supporting rfio-reading:
#EXECUTABLES = pede pede_1GB pede_2GB pede_1GB_rfio pede_2GB_rfio
# If you need also the executables with more memory, you need a
# 64-bit environment and the '-mcmodel=medium' option (see above).
# pede_4GB pede_4GB_rfio pede_8GB pede_8GB_rfio
#EXECUTABLES = pede pede_1GB pede_2GB pede_1GB_rfio pede_2GB_rfio \
# pede_4GB pede_4GB_rfio pede_8GB pede_8GB_rfio
EXECUTABLES = pede
#
all: $(EXECUTABLES)
# The single special one:
pede : ${USER_OBJ_PEDE} Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE} $(C_LIBS)
#
#
-o $@ ${USER_OBJ_PEDE} $(C_LIBS)
#
clean:
rm -f *.o *~ */*.o */*~
rm -f *.o *~ */*.o */*~ *.mod */*.mod
#
clobber: clean
rm -f $(EXECUTABLES)
......@@ -116,102 +84,11 @@ install: $(EXECUTABLES) #clean
# Make the object files - depend on source and include file
#
%.o : %.F Makefile
%.o : %.f90 Makefile
${FCOMP} ${F_FLAGS} -c $< -o $@
# These two depend on the memory defined by NUMBER_OF_WORDS, LARGE_SIZE
# the second also on variable definitions:
dynal.o : dynal.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=${NUMBER_OF_WORDS} -c $< -o $@
pede.o : pede.F dynal.inc largeint.inc mpinds.inc localfit.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=${NUMBER_OF_WORDS} -c $< -o $@
#
%.o: %.c Makefile
$(CCOMP) -c $(C_FLAGS) $(DEFINES) $(C_INCLUDEDIRS) $(DEBUG) -o $@ $<
#
####################################################################
# Here we start the hack for the various executables...
%_1GB.o: %.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=250000000 -c $< -o $@
%_2GB.o: %.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=500000000 -c $< -o $@
%_4GB.o: %.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=1000000000 -c $< -o $@
%_8GB.o: %.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=2000000000 -c $< -o $@
%_16GB.o: %.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=4000000000_8 -c $< -o $@
%_24GB.o: %.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=6000000000_8 -c $< -o $@
%_32GB.o: %.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=8000000000_8 -c $< -o $@
%_48GB.o: %.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=12000000000_8 -c $< -o $@
%_64GB.o: %.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=16000000000_8 -c $< -o $@
%_96GB.o: %.F dynal.inc largeint.inc Makefile
${FCOMP} ${F_FLAGS} -DNUMBER_OF_WORDS=24000000000_8 -c $< -o $@
%_rfio.o: %.c Makefile
$(CCOMP) -c $(C_FLAGS) -DUSE_SHIFT_RFIO $(DEFINES) $(C_INCLUDEDIRS) \
$(DEBUG) -o $@ $<
#
#
#
pede_1GB: ${USER_OBJ_PEDE_STATIC} pede_1GB.o dynal_1GB.o readc.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_1GB.o dynal_1GB.o readc.o $(C_LIBS)
#
pede_2GB: ${USER_OBJ_PEDE_STATIC} pede_2GB.o dynal_2GB.o readc.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_2GB.o dynal_2GB.o readc.o $(C_LIBS)
#
pede_1GB_rfio: ${USER_OBJ_PEDE_STATIC} pede_1GB.o dynal_1GB.o readc_rfio.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_1GB.o dynal_1GB.o readc_rfio.o $(C_LIBS) -lshift
#
pede_2GB_rfio: ${USER_OBJ_PEDE_STATIC} pede_2GB.o dynal_2GB.o readc_rfio.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_2GB.o dynal_2GB.o readc_rfio.o $(C_LIBS) -lshift
#
pede_4GB: ${USER_OBJ_PEDE_STATIC} pede_4GB.o dynal_4GB.o readc.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_4GB.o dynal_4GB.o readc.o $(C_LIBS)
#
pede_4GB_rfio: ${USER_OBJ_PEDE_STATIC} pede_4GB.o dynal_4GB.o readc_rfio.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_4GB.o dynal_4GB.o readc_rfio.o $(C_LIBS) -lshift
#
pede_8GB: ${USER_OBJ_PEDE_STATIC} pede_8GB.o dynal_8GB.o readc.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_8GB.o dynal_8GB.o readc.o $(C_LIBS)
#
pede_8GB_rfio: ${USER_OBJ_PEDE_STATIC} pede_8GB.o dynal_8GB.o readc_rfio.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_8GB.o dynal_8GB.o readc_rfio.o $(C_LIBS) -lshift
#
pede_16GB: ${USER_OBJ_PEDE_STATIC} pede_16GB.o dynal_16GB.o readc.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_16GB.o dynal_16GB.o readc.o $(C_LIBS)
#
pede_24GB: ${USER_OBJ_PEDE_STATIC} pede_24GB.o dynal_24GB.o readc.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_24GB.o dynal_24GB.o readc.o $(C_LIBS)
pede_32GB: ${USER_OBJ_PEDE_STATIC} pede_32GB.o dynal_32GB.o readc.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_32GB.o dynal_32GB.o readc.o $(C_LIBS)
pede_48GB: ${USER_OBJ_PEDE_STATIC} pede_48GB.o dynal_48GB.o readc.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_48GB.o dynal_48GB.o readc.o $(C_LIBS)
pede_64GB: ${USER_OBJ_PEDE_STATIC} pede_64GB.o dynal_64GB.o readc.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_64GB.o dynal_64GB.o readc.o $(C_LIBS)
pede_96GB: ${USER_OBJ_PEDE_STATIC} pede_96GB.o dynal_96GB.o readc.o Makefile
$(LOADER) $(L_FLAGS) \
-o $@ ${USER_OBJ_PEDE_STATIC} pede_96GB.o dynal_96GB.o readc.o $(C_LIBS)
#
# End hack for the various executables...
####################################################################
#
# ##################################################################
# END
# ##################################################################
/**
* \file Mille.cc
* \file
* Create Millepede-II C-binary record.
*
* \author : Gero Flucke
* date : October 2006
* $Revision: 1.3 $
......@@ -14,12 +16,16 @@
//___________________________________________________________________________
/// Opens outFileName (by default as binary file).
/**
* \param[in] outFileName file name
* \param[in] asBinary flag for binary
* \param[in] writeZero flag for keeping of zeros
*/
Mille::Mille(const char *outFileName, bool asBinary, bool writeZero) :
myOutFile(outFileName, (asBinary ? (std::ios::binary | std::ios::out) : std::ios::out)),
myAsBinary(asBinary), myWriteZero(writeZero), myBufferPos(-1), myHasSpecial(false)
{
// opens outFileName, by default as binary file
// Instead myBufferPos(-1), myHasSpecial(false) and the following two lines
// we could call newSet() and kill()...
myBufferInt[0] = 0;
......@@ -32,15 +38,23 @@ Mille::Mille(const char *outFileName, bool asBinary, bool writeZero) :
}
//___________________________________________________________________________
/// Closes file.
Mille::~Mille()
{
// closes file
myOutFile.close();
}
//___________________________________________________________________________
/// Add measurement to buffer.
/**
* \param[in] NLC number of local derivatives
* \param[in] derLc local derivatives
* \param[in] NGL number of global derivatives
* \param[in] derGl global derivatives
* \param[in] label global labels
* \param[in] rMeas measurement (residuum)
* \param[in] sigma error
*/
void Mille::mille(int NLC, const float *derLc,
int NGL, const float *derGl, const int *label,
float rMeas, float sigma)
......@@ -68,7 +82,7 @@ void Mille::mille(int NLC, const float *derLc,
myBufferFloat[myBufferPos] = sigma;
myBufferInt [myBufferPos] = 0;
// store global derivatives and their lables
// store global derivatives and their labels
for (int i = 0; i < NGL; ++i) {
if (derGl[i] || myWriteZero) { // by default store only non-zero derivatives
if ((label[i] > 0 || myWriteZero) && label[i] <= myMaxLabel) { // and for valid labels
......@@ -84,6 +98,12 @@ void Mille::mille(int NLC, const float *derLc,
}
//___________________________________________________________________________
/// Add special data to buffer.
/**
* \param[in] nSpecial number of floats/ints
* \param[in] floatings floats
* \param[in] integers ints
*/
void Mille::special(int nSpecial, const float *floatings, const int *integers)
{
if (nSpecial == 0) return;
......@@ -118,18 +138,16 @@ void Mille::special(int nSpecial, const float *floatings, const int *integers)
}
//___________________________________________________________________________
/// Reset buffers, i.e. kill derivatives accumulated for current set.
void Mille::kill()
{
// reset buffers, i.e. kill derivatives accumulated for current set
myBufferPos = -1;
}
//___________________________________________________________________________
/// Write buffer (set of derivatives with same local parameters) to file.
void Mille::end()
{
// write set of derivatives with same local parameters to file
if (myBufferPos > 0) { // only if anything stored...
const int numWordsToWrite = (myBufferPos + 1)*2;
......@@ -157,10 +175,9 @@ void Mille::end()
}
//___________________________________________________________________________
/// Initialize for new set of locals, e.g. new track.
void Mille::newSet()
{
// initilise for new set of locals, e.g. new track
myBufferPos = 0;
myHasSpecial = false;
myBufferFloat[0] = 0.0;
......@@ -168,11 +185,14 @@ void Mille::newSet()
}
//___________________________________________________________________________
/// Enough space for next nLocal + nGlobal derivatives incl. measurement?
/**
* \param[in] nLocal number of local derivatives
* \param[in] nGlobal number of global derivatives
* \return true if sufficient space available (else false)
*/
bool Mille::checkBufferSize(int nLocal, int nGlobal)
{
// enough space for next nLocal + nGlobal derivatives incl. measurement?
if (myBufferPos + nLocal + nGlobal + 2 >= myBufferSize) {
++(myBufferInt[0]); // increase error count
std::cerr << "Mille::checkBufferSize: Buffer too short ("
......
......@@ -7,13 +7,14 @@
* \class Mille
*
* Class to write a C binary (cf. below) file of a given name and to fill it
* with information used as input to pede.
* Use its member functions mille(...), special(...), kill() and end() as you would
* use the fortran MILLE and its entry points MILLSP, KILLE and ENDLE.
* with information used as input to **pede**.
* Use its member functions \c mille(), \c special(), \c kill() and \c end()
* as you would use the fortran \ref mille.f90 "MILLE"
* and its entry points \c MILLSP, \c KILLE and \c ENDLE.
*
* For debugging purposes constructor flags enable switching to text output and/or
* to write also derivatives and lables which are ==0.
* But note that pede will not be able to read text output and has not been tested with
* to write also derivatives and labels which are ==0.
* But note that **pede** will not be able to read text output and has not been tested with
* derivatives/labels ==0.
*
* \author : Gero Flucke
......@@ -23,6 +24,7 @@
* (last update by $Author: flucke $)
*/
/// Class to write C binary file.
class Mille
{
public:
......@@ -39,16 +41,16 @@ class Mille
void newSet();
bool checkBufferSize(int nLocal, int nGlobal);
std::ofstream myOutFile; // C-binary for output
bool myAsBinary; // if false output as text
bool myWriteZero; // if true also write out derivatives/lables ==0
enum {myBufferSize = 5000};
int myBufferInt[myBufferSize]; // to collect labels etc.
float myBufferFloat[myBufferSize]; // to collect derivatives etc.
int myBufferPos;
bool myHasSpecial; // if true, special(..) already called for this record
enum {myMaxLabel = (0xFFFFFFFF - (1 << 31))}; // largest label allowed: 2^31 - 1
std::ofstream myOutFile; ///< C-binary for output
bool myAsBinary; ///< if false output as text
bool myWriteZero; ///< if true also write out derivatives/labels ==0
/// buffer size for ints and floats
enum {myBufferSize = 5000}; ///< buffer size for ints and floats
int myBufferInt[myBufferSize]; ///< to collect labels etc.
float myBufferFloat[myBufferSize]; ///< to collect derivatives etc.
int myBufferPos; ///< position in buffer
bool myHasSpecial; ///< if true, special(..) already called for this record
/// largest label allowed: 2^31 - 1
enum {myMaxLabel = (0xFFFFFFFF - (1 << 31))};
};
#endif
This diff is collapsed.
This diff is collapsed.
*
INTEGER MQ
INTEGER(KIND=LARGE) MEGA,MEGAH
PARAMETER (MEGA=NUMBER_OF_WORDS) ! make preprocessor decide
c PARAMETER (MEGA=220 000 000) ! 100 mio words: 400 MB
* -----------
COMMON/COMEGA/MQ(MEGA)
REAL QM(MEGA)
PARAMETER (MEGAH=MEGA/2)
DOUBLE PRECISION DQ(MEGAH)
INTEGER(KIND=LARGE) LQ(MEGA/LSCALE)
EQUIVALENCE (MQ(1),QM(1),DQ(1),LQ(1))
INTEGER(KIND=LARGE) MQU,INDICS,ISIZES,
+ IND0,IND1,IND2,IND3,IND4,IND5,IND6,IND7,IND8,IND9,
+ INDA,INDB,INDC,INDD,INDE,INDF,INDG,INDH,INDI,INDJ,INDK,
+ INDL,INDM,INDN,INDO,INDP,INDQ,INDR,INDS,INDT,INDU,INDV,
+ INDW,INDX,INDY,INDZ,IDOT1,IDOT2,IDOT3,IDOT4,
+ LND0,LND1,LND2,LND3,LND4,LND5,LND6,LND7,LND8,LND9,
+ LNDA,LNDB,LNDC,LNDD,LNDE,LNDF,LNDG,LNDH,LNDI,LNDJ,LNDK,
+ LNDL,LNDM,LNDN,LNDO,LNDP,LNDQ,LNDR,LNDS,LNDT,LNDU,LNDV,
+ LNDW,LNDX,LNDY,LNDZ,LDOT1,LDOT2,LDOT3,LDOT4
COMMON/CMINDC/MQU,INDICS(0:39),ISIZES(0:39)
EQUIVALENCE (INDICS( 0),IND0),(INDICS( 1),IND1),(INDICS( 2),IND2),
+ (INDICS( 3),IND3),(INDICS( 4),IND4),(INDICS( 5),IND5),
+ (INDICS( 6),IND6),(INDICS( 7),IND7),(INDICS( 8),IND8),
+ (INDICS( 9),IND9)
EQUIVALENCE (INDICS(10),INDA),(INDICS(11),INDB),(INDICS(12),INDC),
+ (INDICS(13),INDD),(INDICS(14),INDE),(INDICS(15),INDF),
+ (INDICS(16),INDG),(INDICS(17),INDH),(INDICS(18),INDI),
+ (INDICS(19),INDJ),(INDICS(20),INDK),(INDICS(21),INDL),
+ (INDICS(22),INDM),(INDICS(23),INDN),(INDICS(24),INDO),
+ (INDICS(25),INDP),(INDICS(26),INDQ),(INDICS(27),INDR),
+ (INDICS(28),INDS),(INDICS(29),INDT),(INDICS(30),INDU),
+ (INDICS(31),INDV),(INDICS(32),INDW),(INDICS(33),INDX),
+ (INDICS(34),INDY),(INDICS(35),INDZ)
EQUIVALENCE (INDICS(36),IDOT1),(INDICS(37),IDOT2),
+ (INDICS(38),IDOT3),(INDICS(39),IDOT4)
EQUIVALENCE (ISIZES( 0),LND0),(ISIZES( 1),LND1),(ISIZES( 2),LND2),
+ (ISIZES( 3),LND3),(ISIZES( 4),LND4),(ISIZES( 5),LND5),
+ (ISIZES( 6),LND6),(ISIZES( 7),LND7),(ISIZES( 8),LND8),
+ (ISIZES( 9),LND9)
EQUIVALENCE (ISIZES(10),LNDA),(ISIZES(11),LNDB),(ISIZES(12),LNDC),
+ (ISIZES(13),LNDD),(ISIZES(14),LNDE),(ISIZES(15),LNDF),
+ (ISIZES(16),LNDG),(ISIZES(17),LNDH),(ISIZES(18),LNDI),
+ (ISIZES(19),LNDJ),(ISIZES(20),LNDK),(ISIZES(21),LNDL),
+ (ISIZES(22),LNDM),(ISIZES(23),LNDN),(ISIZES(24),LNDO),
+ (ISIZES(25),LNDP),(ISIZES(26),LNDQ),(ISIZES(27),LNDR),
+ (ISIZES(28),LNDS),(ISIZES(29),LNDT),(ISIZES(30),LNDU),
+ (ISIZES(31),LNDV),(ISIZES(32),LNDW),(ISIZES(33),LNDX),
+ (ISIZES(34),LNDY),(ISIZES(35),LNDZ)
EQUIVALENCE (ISIZES(36),LDOT1),(ISIZES(37),LDOT2),
+ (ISIZES(38),LDOT3),(ISIZES(39),LDOT4)
INTEGER LUNLOG,LVLLOG
COMMON/CLULOG/LUNLOG,LVLLOG
* Large Integers are INTEGER(KIND=LARGE)
* if not understood by compiler replace (KIND=LARGE) by *4 or *8 with preprocessor
INTEGER LARGE,LSCALE,MAXI4
PARAMETER (LARGE=LARGE_SIZE) ! size of LARGE integers (*4 or *8) from preprocessor
PARAMETER (LSCALE=LARGE/4) ! number of (32bit) words
PARAMETER (MAXI4=2147483647) ! max. INTEGER*4
*
SUBROUTINE PTLINE(N,X,F,G,S,STEP, INFO) ! - 2 arguments
IMPLICIT NONE
* ------------------------------------------------------------------
* "Line search routine with sufficient decrease of slope"
*
* In many minimization problems the objective function is close to
* quadratic, except far from the solution. Close to the minimum the
* behaviour may be almost quadratic or, due to round-off errors,
* it may have a non-smooth behaviour, which often complicates any
* further progress and the recognition of convergence.
* Round-off errors affect the function value, which may be large and
* small parameter changes result in small relative changes of the
* function value. Close to the minimum the gradient becomes small
* and the behaviour is not so much affected by Round-off errors.
*
* ------------------------------------------------------------------
*
* CALL PTLDEF(0.0,0.0, 0,0) ! init line search
* N=...
* X(.)=...
* D(.)=...
* ALPHA=1.0D0
* 10 F(X)=...
* G(X)=...
* IF(.) S(X)=..
* CALL PTLINE(N,X,F,G,D,ALPHA,INFO)
* IF(INFO.LT.0) GOTO 10
*
* ------------------------------------------------------------------
* N = dimension of problem
* X(N) = current iterate
* F = associated function value
* G(N) = associated gradient
* S(N) = search vector
*
* ALPHA = step factor (initially = 1.0)
*
* INFO = information
* = -1 repeat function evaluation
* = 0 input error (e.g. gradient not negative)
* = 1 convergence reached
* = 2 convergence assumed, but round-off errors
* = 3 too many function calls
* = 4 ALPHA to small (ALPHA <= TOL)
*
* ------------------------------------------------------------------
INTEGER N,INFO,MINFE,MAXFE ! arguments
DOUBLE PRECISION X(N),F,G(N),S(N),STEP ! arguments
REAL GTOLE,STMAX ! arguments
* ------------------------------------------------------------------
C DOUBLE PRECISION AUX(1000)
INTEGER I,IM ! internal
DOUBLE PRECISION ALPHA ! internal
DOUBLE PRECISION DGINIT,DG,FSAVED,TOT,FP1,FP2 ! internal
* DGINIT = initial slope
* DG = current slope
* FSAVE = initla function value
* TOT = total step
* -----------------------------------------------------------------
INTEGER MSFD,NSFD,IDGL,IDGR,IDGM,MINF,MAXF,LSINFO
PARAMETER (MSFD=20)
DOUBLE PRECISION SFD,STMX,GTOL
COMMON/CLNSEA/SFD(4,MSFD),STMX,GTOL,NSFD,IDGL,IDGR,IDGM,MINF,MAXF
+ ,LSINFO
* -----------------------------------------------------------------
INTEGER I1,I2
C DOUBLE PRECISION CENTER,REVAL,COMP
* -----------------------------------------------------------------
* -----------------------------------------------------------------
INTEGER LUNLOG
COMMON/CLULOG/LUNLOG
* -----------------------------------------------------------------
SAVE
*
* initialization ---------------------------------------------------
*
INFO=0 ! reset INFO flag
DG=0.0D0
DO I=1,N !
DG=DG-G(I)*S(I) ! DG = scalar product: grad x search
END DO!
IF(NSFD.EQ.0) THEN ! initial call
DGINIT=DG ! DG = initial directional gradient
IF(DGINIT.GE.0.0D0) GOTO 100 ! error: step not decreasing
STEP=1.0D0 ! initial step factor is one
ALPHA=STEP ! get initial step factor
TOT=0.0D0 ! reset total step
IDGL=1 ! index of smallest negative slope
IDGR=0 ! index of smallest positive slope
FSAVED=F ! initial Function value
NSFD=1 ! starting point of iteration
SFD(1,1)=0.0 ! abscissa
SFD(2,1)=0.0 ! referencefunction value
SFD(3,1)=DGINIT ! slope
SFD(4,1)=0.0 ! predicted zero
IM=1 ! optimum
ELSE ! subsequent call
NSFD=NSFD+1
SFD(1,NSFD)=TOT ! abscissa
SFD(2,NSFD)=F-FSAVED ! function value difference to reference
SFD(3,NSFD)=DG ! slope
SFD(4,NSFD)=0.0 ! predicted zero (see below)
IF(DG.LT.SFD(3,IM)) THEN
IM=NSFD
C DO I=1,N
C AUX(I)=G(I)
C END DO
END IF
* define interval indices IDGL and IDGR
IF(DG.LE.0.0D0) THEN
IF(DG.GE.SFD(3,IDGL)) IDGL=NSFD
END IF
IF(DG.GE.0.0D0) THEN ! limit to the right
IF(IDGR.EQ.0) IDGR=NSFD
IF(DG.LE.SFD(3,IDGR)) IDGR=NSFD
END IF
IF(IDGR.EQ.0) THEN
I1=NSFD-1
I2=NSFD
ELSE
I1=IDGL
I2=IDGR
END IF
FP1=SFD(3,I1)
FP2=SFD(3,I2) ! interpolation
SFD(4,NSFD)=(SFD(1,I1)*FP2-SFD(1,I2)*FP1)/(FP2-FP1)
c CENTER=0.5D0*(SFD(1,I1)+SFD(1,I2))
c REVAL=ABS(2.0D0*MIN(ABS(FP1),ABS(FP2))/ABS(FP1-FP2)-1.0D0)
c IF(ABS(FP1).GT.ABS(FP2)) THEN