pede.f90 306 KB
Newer Older
Claus Kleinwort's avatar
Claus Kleinwort committed
1
2
3
4

! Code converted using TO_F90 by Alan Miller
! Date: 2012-03-16  Time: 11:06:00

5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
!> \file
!! Millepede II program, subroutines.
!!
!! \author Volker Blobel, University Hamburg, 2005-2009 (initial Fortran77 version)
!! \author Gero Flucke, University Hamburg (support of C-type binary files)
!! \author Claus Kleinwort, DESY (maintenance and developement)
!!
!! \copyright
!! Copyright (c) 2009 - 2015 Deutsches Elektronen-Synchroton,
!! Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
!! This library is free software; you can redistribute it and/or modify
!! it under the terms of the GNU Library General Public License as
!! published by the Free Software Foundation; either version 2 of the
!! License, or (at your option) any later version. \n\n
!! This library is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU Library General Public License for more details. \n\n
!! You should have received a copy of the GNU Library General Public
!! License along with this program (see the file COPYING.LIB for more
!! details); if not, write to the Free Software Foundation, Inc.,
!! 675 Mass Ave, Cambridge, MA 02139, USA.
!!

Claus Kleinwort's avatar
Claus Kleinwort committed
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
!> \mainpage Overview
!!
!! \section intro_sec Introduction
!! In certain least squares fit problems with a very large number of parameters
!! the set of parameters can be divided into two classes, global and local parameters.
!! Local parameters are those parameters which are present only in subsets of the
!! data. Detector alignment and calibration based on track fits is one of the problems,
!! where the interest is only in optimal values of the global parameters, the
!! alignment parameters. The method, called Millepede, to solve the linear least
!! squares problem with a simultaneous fit of all global and local parameters,
!! irrespectively of the number of local parameters, is described in the draft manual.
!!
!! The Millepede method and the initial implementation has been
!! developed by [V. Blobel](http://www.desy.de/~blobel) from he University of Hamburg.
!! Meanwhile the code is maintained at DESY by the statistics tools group of the
!! analysis center of the Helmholtz Terascale alliance
!! ([www.terascale.de](https://www.wiki.terascale.de/index.php/Millepede_II)).
!!
47
48
49
!! The Millepede II software is provided by DESY under the terms of the
!! [LGPLv2 license](http://www.gnu.org/licenses/old-licenses/lgpl-2.0-standalone.html).
!!
Claus Kleinwort's avatar
Claus Kleinwort committed
50
51
52
53
54
!! \section install_sec Installation
!! To install **Millepede** (on a linux system):
!! 1. Download the software package from the DESY \c svn server to
!!    \a target directory, e.g.:
!!
55
!!         svn checkout http://svnsrv.desy.de/public/MillepedeII/tags/V04-03-09 target
Claus Kleinwort's avatar
Claus Kleinwort committed
56
57
58
59
60
61
62
63
64
65
66
!!
!! 2. Create **Pede** executable (in \a target directory):
!!
!!         make pede
!!
!! 3. Optionally check the installation by running the simple test case:
!!
!!         ./pede -t
!!
!!    This will create (and use) the necessary text and binary files.
!!
67
68
69
!! \section news_sec News
!! * 131008: New solution method \ref ch-minresqlp "MINRES-QLP"
!! [\ref ref_sec "ref 9"] implemented.
70
!! * 140226: Reading of C binary files containing *doubles* implemented.
71
!! * 141020: Storage of values read from text files as *doubles* implemented.
72
73
!! * 141125: Dynamic entries (from accepted local fits) check implemented.
!!   (Rejection of local fits may lead to the loss of degrees of freedom.)
74
75
!!   Printout of global parameter counters with new command \ref cmd-printcounts.
!! * 141126: Weighted constraints implemented (with new command \ref cmd-weightedcons).
76
77
78
!! * 150210: Solution by elimination for problems with linear equality constraints
!!   has been implemented (as default, new command \ref cmd-withelim) in addition to the
!!   Lagrange multiplier method (new command \ref cmd-withmult).
79
80
81
82
83
84
85
86
87
88
!! * 150218: Skipping *empty* constraints (without variable parameters).
!!   With new command \ref cmd-checkinput detailed check of input data (binary files,
!!   constraints) is performed, but no solution will be determined.
!!   Some input statistics is available in the output file <tt>millepede.res</tt>.
!! * 150226: Iteration of entries cut with new command \ref cmd-iterateentries.
!!   In the second iteration measurements with any parameters fixed by the 
!!   previous entries cut are skipped. Useful if parameters of measurements have
!!   different number of entries. 
!! * 150420: Skipping of empty constraints has to be enabled by new command \ref
!!   cmd-skipemptycons.
Claus Kleinwort's avatar
Claus Kleinwort committed
89
!! * 150901: Preconditioning for MINRES with skyline matrix (avoiding rank deficits 
90
91
92
93
!!   of band matrix) added (selected by second argument in \ref cmd-bandwidth >0).
!! * 150925: Monitoring of residuals per local fit cycle is selected by \ref cmd-monres.
!!   The normalized residuals are grouped by the first global label and the median 
!!   and the RMS (from the median of the absolute deviations) per group are
94
95
96
97
98
99
100
!!   written to <tt>millepede.mon</tt>.
!! * 170502: Monitoring of pulls per local fit cycle is selected by \ref cmd-monpull.
!!   The scaling of measurement errors is enabled by \ref cmd-scaleerrors.
!!   Pede will abort now for constraints with a singular QL decomposition
!!   of the constraints matrix (solution by elemination).
!!   This problem is usually caused by *empty* constraints (see \ref
!!   cmd-skipemptycons).
101
102
!! * 170831: More debug information for problems with reading Cfiles. Don't stop
!!   after read error for \ref cmd-checkinput mode.
103
104
!! * 180525: Some fixes: Proper handling of special (debug) data blocks in binary
!!   files, proper exit code (3) for 'function not decreasing'.
105
!!
Claus Kleinwort's avatar
Claus Kleinwort committed
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
!! \section tools_sec Tools
!! The subdirectory \c tools contains some useful scripts:
!! * \c readMilleBinary.py: Python script to read binary files and print
!!   records in text form.
!! * \c readPedeHists.C: ROOT script to read and convert the **Millepede**
!!   histogram file <tt>millepede.his</tt>.
!!
!! \section details_sec Details
!!
!! Detailed information is available at:
!!
!! \subpage draftman_page
!!
!! \subpage changes_page
!!
!! \subpage option_page
!!
Claus Kleinwort's avatar
Claus Kleinwort committed
123
124
!! \subpage exit_code_page
!!
Claus Kleinwort's avatar
Claus Kleinwort committed
125
126
127
128
129
130
131
132
133
134
135
136
!! \section Contact
!!
!! For information exchange the **Millepede** mailing list
!! anacentre-millepede2@desy.de should be used.
!!
!! \section ref_sec References
!!
!! 1. A New Method for the High-Precision Alignment of Track Detectors,
!!    Volker Blobel and Claus Kleinwort, Proceedings of the Conference on
!!    Adcanced Statistical Techniques in Particle Physics, Durham, 18 - 22 March 2002,
!!    Report DESY 02-077 (June 2002) and
!!    [hep-ex/0208021](http://arxiv.org/abs/hep-ex/0208021)
Claus Kleinwort's avatar
Claus Kleinwort committed
137
!! 2.  Alignment Algorithms, V. Blobel,
Claus Kleinwort's avatar
Claus Kleinwort committed
138
139
!!    [Proceedings](http://cdsweb.cern.ch/search?p=reportnumber%3ACERN-2007-004)
!!    of the LHC Detector Alignment Workshop, September 4 - 6 2006, CERN
Claus Kleinwort's avatar
Claus Kleinwort committed
140
!! 3. Software alignment for Tracking Detectors, V. Blobel,
Claus Kleinwort's avatar
Claus Kleinwort committed
141
142
!!    NIM A, 566 (2006), pp. 5-13,
!!    [doi:10.1016/j.nima.2006.05.157](http://dx.doi.org/10.1016/j.nima.2006.05.157)
Claus Kleinwort's avatar
Claus Kleinwort committed
143
!! 4. A new fast track-fit algorithm based on broken lines, V. Blobel,
Claus Kleinwort's avatar
Claus Kleinwort committed
144
145
!!    NIM A, 566 (2006), pp. 14-17,
!!    [doi:10.1016/j.nima.2006.05.156](http://dx.doi.org/10.1016/j.nima.2006.05.156)
Claus Kleinwort's avatar
Claus Kleinwort committed
146
!! 5. Millepede 2009, V. Blobel, [Contribution]
Claus Kleinwort's avatar
Claus Kleinwort committed
147
148
!!    (https://indico.cern.ch/conferenceOtherViews.py?view=standard&confId=50502)
!!    to the 3rd LHC Detector Alignment Workshop, June 15 - 16 2009, CERN
Claus Kleinwort's avatar
Claus Kleinwort committed
149
!! 6. General Broken Lines as advanced track fitting method, C. Kleinwort,
Claus Kleinwort's avatar
Claus Kleinwort committed
150
151
!!    NIM A, 673 (2012), pp. 107-110,
!!    [doi:10.1016/j.nima.2012.01.024](http://dx.doi.org/10.1016/j.nima.2012.01.024)
152
153
154
!! 7. Volker Blobel und Erich Lohrmann, Statistische und numerische Methoden der
!!    Datenanalyse, Teubner Studienb&uuml;cher, B.G. Teubner, Stuttgart, 1998.
!!    [Online-Ausgabe](http://www.desy.de/~blobel/eBuch.pdf).
155
!! 8. [Systems Optimization Laboratory](http://web.stanford.edu/group/SOL/software/minres),
156
157
158
159
!!    Stanford University;\n
!!    C. C. Paige and M. A. Saunders (1975),
!!    Solution of sparse indefinite systems of linear equations,
!!    SIAM J. Numer. Anal. 12(4), pp. 617-629.
160
!! 9. [Systems Optimization Laboratory](http://web.stanford.edu/group/SOL/software/minresqlp),
161
162
163
164
165
!!    Stanford University;\n
!!    Sou-Cheng Choi, Christopher Paige, and Michael Saunders,
!!    MINRES-QLP: A Krylov subspace method for indefinite or singular
!!    symmetric systems, SIAM Journal of Scientific Computing 33:4, 1810-1836, 2011,
!!    [doi:10.1137/100787921](http://dx.doi.org/10.1137/100787921)
Claus Kleinwort's avatar
Claus Kleinwort committed
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
!!

!> \page changes_page Major changes
!! Major changes with respect to the \ref draftman_page "draft manual".
!! \tableofcontents
!!
!! \section ch-methods Solution methods
!! The following methods to obtain the solution \f$\Vek{x}\f$ from a
!! linear equation system \f$\Vek{A}\cdot\Vek{x}=\Vek{b}\f$ are implemented:
!! \subsection ch-inv Inversion
!! The solution and the covariance matrix \f$\Vek{A}^{-1}\f$ are obtained by
!! \ref an-inv "inversion" of \f$\Vek{A}\f$.
!! Available are the value, error and global correlation for all global parameters.
!! The matrix inversion \ref sqminl "routine" has been \ref ch-openmp "parallelized"
!! and can be used for up to several 10000 parameters.
!! \subsection ch-diag Diagonalization
!! The solution and the covariance matrix \f$\Vek{A}^{-1}\f$ are obtained by
!! \ref an-diag "diagonalization" of \f$\Vek{A}\f$.
!! Available are the value, error, global correlation and
!! eigenvalue (and eigenvector) for all global parameters.
!! \subsection ch-minres Minimal Residual Method (MINRES)
!! The solution is obtained by minimizing \f$\Vert\Vek{A}\cdot\Vek{x}-\Vek{b}\Vert_2\f$
188
!! iteratively. \ref minresmodule::minres "MINRES"  [\ref ref_sec "ref 8"] is a special case of the
Claus Kleinwort's avatar
Claus Kleinwort committed
189
190
191
192
193
194
195
196
197
198
!! generalized minimal residual method (\ref an-gmres "GMRES") for symmetric matrices.
!! Preconditioning with a band matrix of zero or finite
!! \ref mpmod::mbandw "bandwidth" is possible.
!! Individual columns \f$\Vek{c_i}\f$ of the covariance matrix can be calculated by
!! solving \f$\Vek{A}\cdot\Vek{c}_i=\Vek{1}_i\f$ where \f$\Vek{1}_i\f$ is the i-th
!! column on the unit matrix.
!! The most time consuming part (\ref avprod "product" matrix times vector per iteration)
!! has been \ref ch-openmp "parallelized".
!! Available are the value for all (and optionally error, global correlation
!! for few) global parameters.
199
!! \subsection ch-minresqlp Advanced Minimal Residual Method (MINRES-QLP)
200
201
202
203
204
205
206
207
208
209
210
!! The \ref minresqlpmodule::minresqlp "MINRES-QLP" implementation [\ref ref_sec "ref 9"]
!! is a MINRES evolution with improved norm estimates and stopping conditions
!! (leading potentially to different numbers of internal iterations).
!! Internally it uses QLP instead of the QR factorization in
!! MINRES which should be numerically superior and allows to find for
!! singular systems the minimal length (pseudo-inverse) solution.
!!
!! The default behavior is to start (the internal iterations) with QR factorization
!! and to switch to QLP if the (estimated) matrix condition exceeds
!! \ref cmd-mrestranscond "mrtcnd". Pure QR or QLP factorization can be enforced
!! by \ref cmd-mresmode "mrmode".
Claus Kleinwort's avatar
Claus Kleinwort committed
211
!!
212
213
214
215
216
217
218
!! \subsection ch-elim-const Elimination of constraints
!! As alternative to the Lagrange multiplier method the solution by elimination
!! has been added for problems with linear equality constraints.
!! A \ref mpqldec::qldec "QL factorization" (with Householder reflections) of the
!! transposed constraints matrix is used to transform to an unconstrained problem.
!! For sparse matrix storage the sparsity of the global matrix is preserved.
!!
Claus Kleinwort's avatar
Claus Kleinwort committed
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
!! \section ch-regul Regularization
!! Optionally a term \f$\tau\cdot\Vert\Vek{x}\Vert\f$ can be added to the objective function
!! (to be minimized) where \f$\Vek{x}\f$ is the vector of global parameters
!! weighted with the inverse of their individual pre-sigma values.
!!
!! \section ch-locfit Local fit
!! In case the \ref par-locfitv "local fit" is a track fit with proper description of multiple
!! scattering in the detector material additional local parameters have to be introduced
!! for each scatterer and solution by *inversion* can get time consuming
!! (~ \f$n_{lp}^3\f$ for \f$n_{lp}\f$ local parameters). For trajectories based on
!! **broken lines** [\ref ref_sec "ref 4,6"] the corresponding matrix \f$\Vek{\Gamma}\f$
!! has a bordered band structure (\f$\Gamma_{ij}=0\f$ for \f$\min(i,j)>b\f$
!! (border size) and \f$|i-j|>m\f$ (bandwidth)). With
!! <i>root-free Cholesky decomposition</i> the time for the solution is linear
!! and for the calculation of \f$\Gamma^{-1}\f$
!! (needed for the construction of the global matrix) quadratic in \f$n_{lp}\f$.
!! For each local fit the structure of \f$\Vek{\Gamma}\f$ is checked and the faster
!! solution method selected automatically.
!!
!! \section ch-openmp Parallelization
!! The code has been largely parallelized using [OpenMP&tm;](www.openmp.org).
!! This includes the reading of binary files, the local fits, the construction of the
!! sparsity structure and filling of the global matrix and the global fit
!! (except by diagonalization). The number of threads is set by the command
!! \ref cmd-threads.
!!
!! \b Caching. The records are read in blocks into a *read cache* and processed from
!! there in parallel, each record by a single thread. For the filling of the global
!! matrix the (zero-compressed) update matrices (\f$\Vek{\D C}_1+\Vek{\D C}_2\f$ from
!! equations \ref eq-c1 "(15)", \ref eq-c2 "(16)")
!! produced by each local fit are collected in a
!! *write cache*. After processing the block of records this is used to update
!! the global matrix in parallel, each row by a single thread.
!! The total cache size can be changed by the command \ref cmd-cache.
!!
!! \section ch-compression Compressed sparse matrix
!! In sparse storage mode for each row the list of column indices (and values) for the
!! non-zero elements are stored. With compression regions of continous column indices
!! are represented by the first index and their number (packed into a single 32bit
!! integer). Compression is selected by the command \ref cmd-compress.
!! In addition rare elements can be neglected (,histogrammed) or stored in single instead
!! of double precision according to the \ref cmd-pairentries command.
!!
!! \section ch-gzip Gzipped C binary files
!! The [zlib](zlib.net) can be used to directly read *gzipped* C binary files.
!! In this case reading with multiple threads
!! (each file by single thread) can speed up the decompression.
!!
!! \section ch-transf Transformation from FORTRAN77 to Fortran90
!! The **Millepede** source code has been formally transformed from <i>fixed form</i>
!! FORTRAN77 to <i>free form</i> Fortran90 (using TO_F90 by Alan Miller)
!! and (most parts) modernized:
!! - <tt>IMPLICIT NONE</tt> everywhere. Unused variables removed.
!! - \c COMMON blocks replaced by \c MODULEs.
!! - Backward \c GOTOs replaced by proper \c DO loops.
!! - \c INTENT (input/output) of arguments described.
!! - Code documented with doxygen.
!!
!! Unused parts of the code (like the interactive mode) have been removed.
!! The reference compiler for the Fortran90 version is gcc-4.6.2 (gcc-4.4.4 works too).
!!
!! \section ch-memmanage Memory management
!! The memory management for dynamic data structures (matrices, vectors, ..)
!! has been changed from a \ref an-dynal "subdivided" *static* \c COMMON block to
!! *dynamic* (\c ALLOCATABLE) Fortran90 arrays. One **Pede** executable is now
!! sufficient for all application sizes.
!!
!! \section ch-readbuf Read buffer size
!! In the \ref sssec-loop1 "first loop" over all binary files a preset
!! \ref mpmod::ndimbuf "read buffer size" is used. Too large records are skipped,
!! but the maximal record length is still being updated. If any records had to be skipped
!! the read buffer size is afterwards adjusted according to the maximal record length
!! and the first loop is repeated.
!!
!! \section ch-numbin Number of binary files
!! The number of binary files has no hard-coded limit anymore, but is calculated from
!! the steering file and resources (file names, descriptors, ..)
!! are allocated dynamically. Some resources may be limited by the system.
!!

!> \page option_page List of options and commands
!!
!! \tableofcontents
!!
!! \section sec-opt Command line options:
!! \subsection opt-t1 -t
!! Create text and binary files for \ref mptest1.f90 "wire chamber" test case, set
!! \ref mpmod::ictest "ictest" to 1.
!! \subsection opt-t2 -t=track-model
!! Create text and binary files for \ref mptest2.f90 "silicon strip tracker" test case
!! using \a track-models with different accounting for multiple scattering, set
!! \ref mpmod::ictest "ictest" to 2..6.
!! \subsection opt-s -s
!! Solution is not iterated.
!! Automatically switched on in case of rank deficits for constraints.
!! \subsection opt-f -f
!! Force iterating of solution (in case of rank deficits for constraints).
316
317
!! \subsection opt-c -c
!! Check input (binary files, constraints). No solution is determined.
Claus Kleinwort's avatar
Claus Kleinwort committed
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
!!
!! \section sec-cmd Steering file commands:
!! In general the commands are defined by a single line:
!!
!!         keyword   number1  number2  ...
!!
!! For those specifying \ref sssec-parinf "properties" of the global parameters
!! (\a keyword = \c parameter, \c constraint or \c measurement)
!! for each involved global parameter (identified by a \ref an-glolab "label")
!! one additional line follows:
!!
!!         label     number1  number2  ...
!!
!! Default values for the numerical arguments are shown in
!! the command descriptions in '[]'. Missing arguments without default
!! values have no effect.
!!
!! \subsection cmd-bandwidth bandwidth
!! Set band width \ref mpmod::mbandw "mbandw" for
337
!! \ref minresmodule::minres "MINRES" preconditioner to \a number1 [0]
Claus Kleinwort's avatar
Claus Kleinwort committed
338
!! and additional flag \ref mpmod::lprecm "lprecm" to \a number2 [0].
Claus Kleinwort's avatar
Claus Kleinwort committed
339
340
341
342
343
!! \subsection cmd-cache cache
!! Set (read+write) cache size \ref mpmod::ncache "ncache" to \a number1.
!! Define cache size and average fill level.
!! \subsection cmd-cfiles Cfiles
!! Following binaries are C files.
344
345
346
!! \subsection cmd-checkinput checkinput
!! Set check input flag \ref mpmod::icheck "icheck" to 1 (true).
!! Same as \ref opt-c "-c".
Claus Kleinwort's avatar
Claus Kleinwort committed
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
!! \subsection cmd-chisqcut chisqcut
!! For local fit \ref an-chisq "setChi^2" cut \ref mpmod::chicut "chicut" to \a number1 [1.],
!! \ref mpmod::chirem "chirem" to \a number2 [1.].
!! \subsection cmd-compress compress
!! Set compression flag \ref mpmod::mcmprs "mcmprs" for \ref mpbits.f90 "sparse storage"
!! to 1 (true) (and \ref mpmod::msngpe "msngpe" to 1).
!! \subsection cmd-constraint constraint
!! Define \ref sssec_consinf "constraints" for global parameters.
!! \subsection cmd-debug debug
!! Set number of records with debug printout \ref mpmod::mdebug "mdebug" to
!! \a number1 [3], number of measurements with printout \ref mpmod::mdebg2 "mdebg2" to \a number2.
!! \subsection cmd-dwfractioncut dwfractioncut
!! Set \ref an-dwcut "down-weighting fraction" cut \ref mpmod::dwcut "dwcut"
!! to \a number1 (max. 0.5).
!! \subsection cmd-entries entries
362
!! Set \ref an-entries "entries" cuts for variable global parameter
363
364
365
!! \ref mpmod::mreqenf "mreqenf" to \a number1 [25],
!! \ref mpmod::mreqena "mreqena" to \a number2 [10] and
!! \ref mpmod::iteren "iteren" to the product of \a number1 and \a number3 [0].
Claus Kleinwort's avatar
Claus Kleinwort committed
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
!! \subsection cmd-errlabels errlabels
!! Define (up to 100 in total) global labels \a number1 .. \a numberN
!! for which the parameter errors are calculated for method MINRES too
!! (by \ref solglo "solving" \f$\Vek{C}\cdot\Vek{x}_i = \Vek{b}^i, b^i_j = \delta_{ij} \f$).
!! \subsection cmd-force force
!! Set force (iterations) flag \ref mpmod::iforce "iforce" to 1 (true).
!! Same as \ref opt-f "-f".
!! \subsection cmd-fortranfiles fortranfiles
!! Following binaries are Fortran files.
!! \subsection cmd-globalcorr globalcorr
!! Set flag \ref mpmod::igcorr "igcorr" for output of global correlations to 1 (true).
!! \subsection cmd-histprint histprint
!! Set flag \ref mpmod::nhistp "nhistp" for \ref an-histpr "histogram printout"
!! to 1 (true).
!! \subsection cmd-hugecut hugecut
!! For local fit set Chi^2 cut \ref mpmod::chhuge "chhuge"
!! for \ref sssec-outlierdeb "unreasonable data" to \a number1 [1.].
383
384
385
386
387
!! \subsection cmd-iterateentries iterateentries
!! Set maximum value \ref mpmod::iteren "iteren" for iteration of entries cut to
!! \a number1 [maxint]. Can alternatively be set by the \ref cmd-entries command.
!! For parameters with less entries the cut will be iterated ignoring measurements with
!! at least one parameter below \ref mpmod::mreqenf "mreqenf".
Claus Kleinwort's avatar
Claus Kleinwort committed
388
389
390
!! \subsection cmd-linesearch linesearch
!! The mode \ref mpmod::lsearch "lsearch" of the \ref par-linesearch "line search"
!! to improve the solution is set to \a number1.
Claus Kleinwort's avatar
Claus Kleinwort committed
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
!! \subsection cmd-localfit localfit
!! For local fit set number of iterations \ref mpmod::lfitnp "lfitnp"
!! with calculation of pulls to \a number1, flag \ref mpmod::lfitbb "lfitbb"
!! for auto-detection of bordered band matrices to \a number2.
!! \subsection cmd-matiter matiter
!! Set number of iterations \ref mpmod::matrit "matrit" with (re)calcuation of
!! global matrix to \a number1.
!! \subsection cmd-matmoni matmoni
!! Set record interval \ref mpmod::matmon "matmon" for monitoring of (sparse) matrix
!! construction to \a number1.
!! \subsection cmd-maxrecord maxrecord
!! Set record limit \ref mpmod::mxrec "mxrec" to \a number1.
!! \subsection cmd-measurement measurement
!! Define (additional) \ref sssec_gpm "measurements" for global parameters.
!! \subsection cmd-memorydebug memorydebug
!! Set debug flag \ref mpmod::memdbg "memdbg" for memory management
!! to \a number1 [1].
!! \subsection cmd-method method
!! Has special format:
!!
!!         method   name     number1  number2
!!
413
!! Set \ref ch-methods "solution method" \ref mpmod::metsol "metsol" and
Claus Kleinwort's avatar
Claus Kleinwort committed
414
415
!! storage mode \ref mpmod::matsto "matsto" according to \a name,
!! (\c inversion : (1,1), \c diagonalization : (2,1),
416
417
!! \c fullMINRES : (3,1) or \c sparseMINRES : (3,2),
!! \c fullMINRES-QLP : (4,1) or \c sparseMINRES-QLP : (4,2)),
Claus Kleinwort's avatar
Claus Kleinwort committed
418
419
!! (minimum) number of iterations \ref mpmod::mitera "mitera" to \a number1,
!! convergence limit \ref mpmod::dflim "dflim" to \a number2.
420
!! \subsection cmd-monres monitorresiduals
Claus Kleinwort's avatar
Claus Kleinwort committed
421
422
!! Set flag \ref mpmod::imonit "imonit" for monitoring of residuals to \a number1 [3]
!! and increase number of bins (of size 0.1) for internal storage to \a number2 [100].
423
424
425
426
427
!! Monitoring mode \ref mpmod::imonmd "imonmd" is 0.
!! \subsection cmd-monpull monitorpulls
!! Set flag \ref mpmod::imonit "imonit" for monitoring of pulls to \a number1 [3]
!! and increase number of bins (of size 0.1) for internal storage to \a number2 [100].
!! Monitoring mode \ref mpmod::imonmd "imonmd" is 1.
428
429
430
431
432
433
!! \subsection cmd-mresmode mresmode
!! Set \ref minresqlpmodule::minresqlp "MINRES-QLP" factorization mode
!!  \ref mpmod::mrmode "mrmode" to \a number1.
!! \subsection cmd-mrestranscond mrestranscond
!! Set \ref minresqlpmodule::minresqlp "MINRES-QLP" transition (matrix) condition
!!  \ref mpmod::mrtcnd "mrtcnd" to \a number1.
Claus Kleinwort's avatar
Claus Kleinwort committed
434
!! \subsection cmd-mrestol mrestol
435
!! Set tolerance criterion \ref mpmod::mrestl "mrestl" for \ref minresmodule::minres "MINRES"
Claus Kleinwort's avatar
Claus Kleinwort committed
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
!! to \a number1 (\f$10^{-10}\f$ .. \f$10^{-4}\f$).
!! \subsection cmd-nofeasiblestart nofeasiblestart
!! Set flag \ref mpmod::nofeas "nofeas" for \ref an-nofeas "skipping"
!! making parameters feasible to \a number1 [1].
!! \subsection cmd-outlierdownweighting outlierdownweighting
!! For local fit set number of \ref sssec-outlow "outlier"
!! \ref an-downw "down-weighting" iterations
!! \ref mpmod::lhuber "lhuber" to \a number1.
!! \subsection cmd-pairentries pairentries
!! Set entries cut for variable global parameter pairs \ref mpmod::mreqpe "mreqpe"
!! to \a number1, histogram upper bound \ref mpmod::mhispe "mhispe" for pairs
!! to \a number2 (<1: no histogramming), upper bound \ref mpmod::msngpe "msngpe"
!! for pair entries with single precision storage
!! to \a number3.
!! \subsection cmd-parameter parameter
!! Define \ref sssec-parinf "initial value, pre-sigma" for global parameters.
!! \subsection cmd-presigma presigma
!! Set default pre-sigma \ref mpmod::regpre "regpre" to \a number1 [1].
!! \subsection cmd-print print
!! Set print level \ref mpmod::mprint "mprint" to \a number1 [1].
456
!! \subsection cmd-printcounts printcounts
457
458
459
!! Set flag \ref mpmod::ipcntr "ipcntr" to \a number1 [1].
!! The counters for the global parameters from the accepted local fits (=1)
!! or from the binary files (>1) will be printed in the result file.
Claus Kleinwort's avatar
Claus Kleinwort committed
460
461
462
463
464
465
466
467
468
469
470
471
472
!! \subsection cmd-printrecord printrecord
!! \ref an-recpri "Record" numbers with printout.
!! \subsection cmd-pullrange pullrange
!! Set (symmetric) range \ref mpmod::prange "prange" for histograms
!! of pulls, normalized residuals to \a number1 (=0: auto-ranging).
!! \subsection cmd-regularisation regularisation
!! Set flag \ref mpmod::nregul "nregul" for regularization to 1 (true),
!! regularization parameter \ref mpmod::regula "regula" to \a number2,
!! default pre-sigma \ref mpmod::regpre "regpre" to \a number3.
!! \subsection cmd-regularization regularization
!! Set flag \ref mpmod::nregul "nregul" for regularization to 1 (true),
!! regularization parameter \ref mpmod::regula "regula" to \a number2,
!! default pre-sigma \ref mpmod::regpre "regpre" to \a number3.
473
474
475
476
477
!! \subsection cmd-scaleerrors scaleerrors
!! Set measurement scaling factors \ref mpmod::dscerr "dscerr"
!! to \a number1 [1.] and \a number2 [\a number1].
!! First value is for "global" measurements (with global derivatives),
!! second for "local" measurements (without global derivatives).
478
479
480
!! \subsection cmd-skipemptycons skipemptycons
!! Set flag \ref mpmod::iskpec "iskpec" to 1 (true).
!! Empty constraints (without variable parameters) will be skipped.
Claus Kleinwort's avatar
Claus Kleinwort committed
481
482
483
484
485
486
487
488
!! \subsection cmd-subito subito
!! Set subito (no iterations) flag \ref mpmod::isubit "isubit" to 1 (true).
!! Same as \ref opt-s "-s".
!! \subsection cmd-threads threads
!! Set number \ref mpmod::mthrd "mthrd" of OpenMP&tm; threads for processing
!! to \a number1,
!! number \ref mpmod::mthrdr "mthrdr" of threads for reading
!! binary files to \a number2 [\a number1].
489
490
491
!! \subsection cmd-weightedcons weightedcons
!! Set flag \ref mpmod::iwcons "iwcons" to \a number1 [1].
!! Implements \ref sssec_consinf "weighted constraints" for global parameters.
492
493
494
495
496
497
!! \subsection cmd-withelim withelimination
!! Set flag \ref mpmod::icelim "icelim" to 1 (true).
!! Selects solution by elimination for linear equality constraints.
!! \subsection cmd-withmult withmultipliers
!! Set flag \ref mpmod::icelim "icelim" to 0 (false).
!! Selects solution by Lagrange multipliers for linear equality constraints.
Claus Kleinwort's avatar
Claus Kleinwort committed
498
!! \subsection cmd-wolfe wolfe
Claus Kleinwort's avatar
Claus Kleinwort committed
499
!! For strong Wolfe condition in \ref par-linesearch "line search"
Claus Kleinwort's avatar
Claus Kleinwort committed
500
501
502
!! set parameter \ref mpmod::wolfc1 "wolfc1" to \a number1, \ref mpmod::wolfc2
!! "wolfc2" to \a number2.

Claus Kleinwort's avatar
Claus Kleinwort committed
503
504
505
506
507
508
!> \page exit_code_page List of exit codes
!! The exit code and message of the **Pede** executable can be found in the
!! file <tt>millepede.end</tt> :
!!    + <b>-1</b>   Still running or crashed
!!    + **00**   Ended normally
!!    + **01**   Ended with warnings (bad measurements)
509
510
!!    + **02**   Ended with severe warnings (insufficient measurements)
!!    + **03**   Ended with severe warnings (bad global matrix)
Claus Kleinwort's avatar
Claus Kleinwort committed
511
512
513
514
515
516
517
518
!!    + **10**   Aborted, no steering file
!!    + **11**   Aborted, open error for steering file
!!    + **12**   Aborted, second text file in command line
!!    + **13**   Aborted, unknown keywords in steering file
!!    + **14**   Aborted, no binary files
!!    + **15**   Aborted, open error(s) for binary files
!!    + **16**   Aborted, open error(s) for text files
!!    + **17**   Aborted, file name too long
519
!!    + **18**   Aborted, read error(s) for binary files
Claus Kleinwort's avatar
Claus Kleinwort committed
520
521
522
523
524
525
526
!!    + **20**   Aborted, bad binary records
!!    + **21**   Aborted, no labels/parameters defined
!!    + **22**   Aborted, no variable global parameters
!!    + **23**   Aborted, bad matrix index
!!    + **24**   Aborted, vector/matrix size mismatch
!!    + **25**   Aborted, result vector contains NaNs
!!    + **26**   Aborted, too many rejects
527
!!    + **27**   Aborted, singular QL decomposition of constraints matrix
Claus Kleinwort's avatar
Claus Kleinwort committed
528
529
530
531
!!    + **30**   Aborted, memory allocation failed
!!    + **31**   Aborted, memory deallocation failed
!!    + **32**   Aborted, iteration limit reached in diagonalization
!!    + **33**   Aborted, stack overflow in quicksort
532
!!    + **34**   Aborted, pattern string too long
Claus Kleinwort's avatar
Claus Kleinwort committed
533

Claus Kleinwort's avatar
Claus Kleinwort committed
534
535
536
537
538
539
540
541
!> Millepede II main program \ref sssec-stalone "Pede".
PROGRAM mptwo
    USE mpmod
    USE mpdalc
    USE mptest1, ONLY: nplan,del,dvd
    USE mptest2, ONLY: nlyr,nmx,nmy,sdevx,sdevy

    IMPLICIT NONE
542
543
544
545
546
547
548
549
    REAL(mps) :: andf
    REAL(mps) :: c2ndf
    REAL(mps) :: deltat
    REAL(mps) :: diff
    REAL(mps) :: err
    REAL(mps) :: gbu
    REAL(mps) :: gmati
    REAL(mps) :: rej
Claus Kleinwort's avatar
Claus Kleinwort committed
550
551
552
    REAL :: rloop1
    REAL :: rloop2
    REAL :: rstext
553
    REAL(mps) :: secnd
Claus Kleinwort's avatar
Claus Kleinwort committed
554
    REAL :: rst
555
    REAL :: rstp
Claus Kleinwort's avatar
Claus Kleinwort committed
556
    REAL, DIMENSION(2) :: ta
557
558
559
560
561
562
563
564
    INTEGER(mpi) :: i
    INTEGER(mpi) :: ii
    INTEGER(mpi) :: ix
    INTEGER(mpi) :: ixv
    INTEGER(mpi) :: iy
    INTEGER(mpi) :: k
    INTEGER(mpi) :: kfl
    INTEGER(mpi) :: lun
Claus Kleinwort's avatar
Claus Kleinwort committed
565
566
    INTEGER :: minut
    INTEGER :: nhour
567
568
569
570
571
    INTEGER(mpi) :: nmxy
    INTEGER(mpi) :: nrc
    INTEGER(mpi) :: nsecnd
    INTEGER(mpi) :: ntot
    INTEGER(mpi) :: ntsec
Claus Kleinwort's avatar
Claus Kleinwort committed
572
573
574
575

    CHARACTER (LEN=24) :: chdate
    CHARACTER (LEN=24) :: chost

576
577
    INTEGER(mpl) :: rows
    INTEGER(mpl) :: cols
Claus Kleinwort's avatar
Claus Kleinwort committed
578

579
580
581
582
    REAL(mpd) :: sums(9)
    !$    INTEGER(mpi) :: OMP_GET_NUM_PROCS,OMP_GET_MAX_THREADS
    !$    INTEGER(mpi) :: MXTHRD
    !$    INTEGER(mpi) :: NPROC
Claus Kleinwort's avatar
Claus Kleinwort committed
583
584
585
586
587
588

    SAVE
    !     ...
    CALL etime(ta,rstp)
    CALL fdate(chdate)

589
590
    !     millepede monitoring file
    lunmon=0
Claus Kleinwort's avatar
Claus Kleinwort committed
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
    !     millepede.log file
    lunlog=8
    lvllog=1
    CALL mvopen(lunlog,'millepede.log')
    CALL getenv('HOSTNAME',chost)
    IF (chost(1:1) == ' ') CALL getenv('HOST',chost)
    WRITE(*,*) '($Rev$)'
    !$    WRITE(*,*) 'using OpenMP (TM)'
#ifdef __GFORTRAN__
    WRITE(*,111)  __GNUC__ , __GNUC_MINOR__ , __GNUC_PATCHLEVEL__
111 FORMAT(' compiled with gcc ',i0,'.',i0,'.',i0)
#endif
    WRITE(*,*) ' '
    WRITE(*,*) '  <  Millepede II-P starting ... ',chdate
    WRITE(*,*) '                                 ',chost
    WRITE(*,*) ' '

    WRITE(8,*) ' '
    WRITE(8,*) 'Log-file Millepede II-P                        ', chdate
    WRITE(8,*) '                                               ', chost
Claus Kleinwort's avatar
Claus Kleinwort committed
611
    CALL peend(-1,'Still running or crashed')
Claus Kleinwort's avatar
Claus Kleinwort committed
612
613
614
615
    !     read command line and text files

    CALL filetc   ! command line and steering file analysis
    CALL filetx   ! read text files
616
617
618
619
    IF (icheck > 0) THEN
        WRITE(*,*) '!!!   Checking input only, no calculation of a solution   !!!'
        WRITE(8,*) '!!!   Checking input only, no calculation of a solution   !!!'
    END IF
Claus Kleinwort's avatar
Claus Kleinwort committed
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
    lvllog=mprint ! export print level
    IF (memdbg > 0) printflagalloc=1 ! debug memory management
    !$    WRITE(*,*)
    !$    NPROC=1
    !$    MXTHRD=1
    !$    NPROC=OMP_GET_NUM_PROCS()         ! number of processors available
    !$    CALL OMP_SET_NUM_THREADS(MTHRD)   ! set max number of threads to MTHRD
    !$    MXTHRD=OMP_GET_MAX_THREADS()      ! get max number of threads back
    !$    WRITE(*,*) 'Number of processors available:   ', NPROC
    !$    WRITE(*,*) 'Maximum number of OpenMP threads: ', MXTHRD
    !$    WRITE(*,*) 'Number of threads for processing: ', MTHRD
    !$    IF (MXREC.GT.0) MTHRDR=1          ! to get allways the same MXREC records
    !$    WRITE(*,*) 'Number of threads for reading:    ', MTHRDR
    !$POMP INST INIT                        ! start profiling with ompP
    IF (ncache < 0) THEN
        ncache=25000000*mthrd  ! default cache size (100 MB per thread)
    ENDIF
    rows=6; cols=mthrdr
    CALL mpalloc(readBufferInfo,rows,cols,'read buffer header')
    !     histogram file
    lun=7
    CALL mvopen(lun,'millepede.his')
    CALL hmplun(lun) ! unit for histograms
    CALL gmplun(lun) ! unit for xy data

    !     debugging
    IF(nrecpr /= 0.OR.nrecp2 /= 0) THEN
        CALL mvopen(1,'mpdebug.txt')
    END IF

    CALL etime(ta,rstext)
    times(0)=rstext-rstp ! time for text processing

    !     preparation of data sub-arrays

    CALL loop1
    CALL etime(ta,rloop1)
    times(1)=rloop1-rstext ! time for LOOP1

    CALL loop2
    IF(chicut /= 0.0) THEN
        WRITE(8,*) 'Chi square cut equiv 3 st.dev applied ...'
        WRITE(8,*) ' in  first iteration with factor',chicut
        WRITE(8,*) ' in second iteration with factor',chirem
        WRITE(8,*) ' (reduced by sqrt in next iterations)'
    END IF
666
    
Claus Kleinwort's avatar
Claus Kleinwort committed
667
668
669
670
671
672
673
674
    IF(lhuber /= 0) THEN
        WRITE(8,*) 'Down-weighting of outliers in', lhuber,' iterations'
        WRITE(8,*) 'Cut on downweight fraction',dwcut
    END IF

    CALL etime(ta,rloop2)
    times(2)=rloop2-rloop1 ! time for LOOP2

675
676
677
678
679
680
    IF(icheck > 0) THEN
        CALL prtstat
        CALL peend(0,'Ended normally')
        GOTO 99 ! only checking input
    END IF

Claus Kleinwort's avatar
Claus Kleinwort committed
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
    !     use different solution methods

    CALL mstart('Iteration')   ! Solution module starting

    CALL xloopn                ! all methods

    !     ------------------------------------------------------------------

    IF(nloopn > 2.AND.nhistp /= 0) THEN       ! last iteration
        CALL hmprnt(3)  ! scaled residual of single measurement (with global deriv.)
        CALL hmprnt(12) ! scaled residual of single measurement (no global deriv.)
        CALL hmprnt(4)  ! chi^2/Ndf
    END IF
    IF(nloopn > 2) THEN
        CALL hmpwrt(3)
        CALL hmpwrt(12)
        CALL hmpwrt(4)
        CALL gmpwrt(4) ! location, dispersion (res.) as a function of record nr
        IF (nloopn <= lfitnp) THEN
            CALL hmpwrt(13)
            CALL hmpwrt(14)
            CALL gmpwrt(5)
        END IF
    END IF
    IF(nhistp /= 0) THEN
        CALL gmprnt(1)
        CALL gmprnt(2)
    END IF
    CALL gmpwrt(1)             ! output of xy data
    CALL gmpwrt(2)             ! output of xy data
    !     'track quality' per binary file
    IF (nfilb > 1) THEN
        CALL gmpdef(6,1,'log10(#records) vs file number')
        CALL gmpdef(7,1,'final rejection fraction vs file number')
        CALL gmpdef(8,1,  &
716
            'final <Chi^2/Ndf> from accepted local fits vs file number')
Claus Kleinwort's avatar
Claus Kleinwort committed
717
718
719
720
721
722
        CALL gmpdef(9,1, '<Ndf> from accepted local fits vs file number')
  
        DO i=1,nfilb
            kfl=kfd(2,i)
            nrc=-kfd(1,i)
            IF (nrc > 0) THEN
723
724
725
                rej=REAL(nrc-jfd(kfl),mps)/REAL(nrc,mps)
                CALL gmpxy(6,REAL(kfl,mps),LOG10(REAL(nrc,mps))) ! log10(#records) vs file
                CALL gmpxy(7,REAL(kfl,mps),rej)    ! rejection fraction vs file
Claus Kleinwort's avatar
Claus Kleinwort committed
726
727
            END IF
            IF (jfd(kfl) > 0) THEN
728
729
730
731
                c2ndf=cfd(kfl)/REAL(jfd(kfl),mps)
                CALL gmpxy(8,REAL(kfl,mps),c2ndf)  ! <Chi2/NDF> vs file
                andf=REAL(dfd(kfl),mps)/REAL(jfd(kfl),mps)
                CALL gmpxy(9,REAL(kfl,mps),andf)  ! <NDF> vs file
Claus Kleinwort's avatar
Claus Kleinwort committed
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
            END IF
        END DO
        IF(nhistp /= 0) THEN
            CALL gmprnt(6)
            CALL gmprnt(7)
            CALL gmprnt(8)
            CALL gmprnt(9)
        END IF
        CALL gmpwrt(6)             ! output of xy data
        CALL gmpwrt(7)             ! output of xy data
        CALL gmpwrt(8)             ! output of xy data
        CALL gmpwrt(9)             ! output of xy data
    END IF

    IF(ictest == 1) THEN
        WRITE(*,*) ' '
        WRITE(*,*) 'Misalignment test wire chamber'
        WRITE(*,*) ' '
  
        CALL hmpdef( 9,-0.0015,+0.0015,'True - fitted displacement')
        CALL hmpdef(10,-0.0015,+0.0015,'True - fitted Vdrift')
        DO i=1,4
754
            sums(i)=0.0_mpd
Claus Kleinwort's avatar
Claus Kleinwort committed
755
756
        END DO
        DO i=1,nplan
757
            diff=REAL(-del(i)-globalParameter(i),mps)
Claus Kleinwort's avatar
Claus Kleinwort committed
758
759
            sums(1)=sums(1)+diff
            sums(2)=sums(2)+diff*diff
760
            diff=REAL(-dvd(i)-globalParameter(100+i),mps)
Claus Kleinwort's avatar
Claus Kleinwort committed
761
762
763
            sums(3)=sums(3)+diff
            sums(4)=sums(4)+diff*diff
        END DO
764
765
766
767
        sums(1)=0.01_mpd*sums(1)
        sums(2)=SQRT(0.01_mpd*sums(2))
        sums(3)=0.01_mpd*sums(3)
        sums(4)=SQRT(0.01_mpd*sums(4))
Claus Kleinwort's avatar
Claus Kleinwort committed
768
769
770
771
772
773
774
775
776
        WRITE(*,143) 'Parameters   1 - 100: mean =',sums(1), 'rms =',sums(2)
        WRITE(*,143) 'Parameters 101 - 200: mean =',sums(3), 'rms =',sums(4)
143     FORMAT(6X,a28,f9.6,3X,a5,f9.6)
        WRITE(*,*) ' '
        WRITE(*,*) ' '
        WRITE(*,*) '    I '
        WRITE(*,*) '   --- '
        DO i=1,100
            WRITE(*,102) i,-del(i),globalParameter(i),-del(i)-globalParameter(i),  &
777
                -dvd(i),globalParameter(100+i),-dvd(i)-globalParameter(100+i)
778
            diff=REAL(-del(i)-globalParameter(i),mps)
Claus Kleinwort's avatar
Claus Kleinwort committed
779
            CALL hmpent( 9,diff)
780
            diff=REAL(-dvd(i)-globalParameter(100+i),mps)
Claus Kleinwort's avatar
Claus Kleinwort committed
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
            CALL hmpent(10,diff)
        END DO
        IF(nhistp /= 0) THEN
            CALL hmprnt( 9)
            CALL hmprnt(10)
        END IF
        CALL hmpwrt( 9)
        CALL hmpwrt(10)
    END IF
    IF(ictest > 1) THEN
        WRITE(*,*) ' '
        WRITE(*,*) 'Misalignment test Si tracker'
        WRITE(*,*) ' '
  
        CALL hmpdef( 9,-0.0025,+0.0025,'True - fitted displacement X')
        CALL hmpdef(10,-0.025,+0.025,'True - fitted displacement Y')
        DO i=1,9
798
            sums(i)=0.0_mpd
Claus Kleinwort's avatar
Claus Kleinwort committed
799
800
801
802
803
804
805
        END DO
        nmxy=nmx*nmy
        ix=0
        iy=ntot
        DO i=1,nlyr
            DO k=1,nmxy
                ix=ix+1
806
807
                diff=REAL(-sdevx((i-1)*nmxy+k)-globalParameter(ix),mps)
                sums(1)=sums(1)+1.0_mpd
Claus Kleinwort's avatar
Claus Kleinwort committed
808
809
810
811
812
                sums(2)=sums(2)+diff
                sums(3)=sums(3)+diff*diff
                ixv=globalParLabelIndex(2,ix)
                IF (ixv > 0.AND.metsol == 1.OR.metsol == 2) THEN
                    ii=(ixv*ixv+ixv)/2
813
                    gmati=REAL(globalMatD(ii),mps)
Claus Kleinwort's avatar
Claus Kleinwort committed
814
815
                    ERR=SQRT(ABS(gmati))
                    diff=diff/ERR
816
                    sums(7)=sums(7)+1.0_mpd
Claus Kleinwort's avatar
Claus Kleinwort committed
817
818
819
820
821
822
823
                    sums(8)=sums(8)+diff
                    sums(9)=sums(9)+diff*diff
                END IF
            END DO
            IF (MOD(i,3) == 1) THEN
                DO k=1,nmxy
                    iy=iy+1
824
825
                    diff=-REAL(sdevy((i-1)*nmxy+k)-globalParameter(iy),mps)
                    sums(4)=sums(4)+1.0_mpd
Claus Kleinwort's avatar
Claus Kleinwort committed
826
827
828
829
830
                    sums(5)=sums(5)+diff
                    sums(6)=sums(6)+diff*diff
                    ixv=globalParLabelIndex(2,iy)
                    IF (ixv > 0.AND.metsol == 1.OR.metsol == 2) THEN
                        ii=(ixv*ixv+ixv)/2
831
                        gmati=REAL(globalMatD(ii),mps)
Claus Kleinwort's avatar
Claus Kleinwort committed
832
833
                        ERR=SQRT(ABS(gmati))
                        diff=diff/ERR
834
                        sums(7)=sums(7)+1.0_mpd
Claus Kleinwort's avatar
Claus Kleinwort committed
835
836
837
838
839
840
841
842
843
844
845
846
                        sums(8)=sums(8)+diff
                        sums(9)=sums(9)+diff*diff
                    END IF
                END DO
            END IF
        END DO
        sums(2)=sums(2)/sums(1)
        sums(3)=SQRT(sums(3)/sums(1))
        sums(5)=sums(5)/sums(4)
        sums(6)=SQRT(sums(6)/sums(4))
        WRITE(*,143) 'Parameters   1 - 500: mean =',sums(2), 'rms =',sums(3)
        WRITE(*,143) 'Parameters 501 - 700: mean =',sums(5), 'rms =',sums(6)
847
        IF (sums(7) > 0.5_mpd) THEN
Claus Kleinwort's avatar
Claus Kleinwort committed
848
849
850
851
852
853
854
855
856
857
858
859
860
            sums(8)=sums(8)/sums(7)
            sums(9)=SQRT(sums(9)/sums(7))
            WRITE(*,143) 'Parameter pulls, all: mean =',sums(8), 'rms =',sums(9)
        END IF
        WRITE(*,*) ' '
        WRITE(*,*) ' '
        WRITE(*,*) '    I '
        WRITE(*,*) '   --- '
        ix=0
        iy=ntot
        DO i=1,nlyr
            DO k=1,nmxy
                ix=ix+1
861
                diff=REAL(-sdevx((i-1)*nmxy+k)-globalParameter(ix),mps)
Claus Kleinwort's avatar
Claus Kleinwort committed
862
863
864
865
866
867
868
869
                CALL hmpent( 9,diff)
                WRITE(*,102) ix,-sdevx((i-1)*nmxy+k),globalParameter(ix),-diff
            END DO
        END DO
        DO i=1,nlyr
            IF (MOD(i,3) == 1) THEN
                DO k=1,nmxy
                    iy=iy+1
870
                    diff=REAL(-sdevy((i-1)*nmxy+k)-globalParameter(iy),mps)
Claus Kleinwort's avatar
Claus Kleinwort committed
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
                    CALL hmpent(10,diff)
                    WRITE(*,102) iy,-sdevy((i-1)*nmxy+k),globalParameter(iy),-diff
                END DO
            END IF
        END DO
        IF(nhistp /= 0) THEN
            CALL hmprnt( 9)
            CALL hmprnt(10)
        END IF
        CALL hmpwrt( 9)
        CALL hmpwrt(10)
    END IF

    IF(nrec1+nrec2 > 0) THEN
        WRITE(8,*) ' '
        IF(nrec1 > 0) THEN
            WRITE(8,*) 'Record',nrec1,' has largest residual:',value1
        END IF
        IF(nrec2 > 0) THEN
            WRITE(8,*) 'Record',nrec2,' has largest Chi^2/Ndf:',value2
        END IF
    END IF
893
    IF(nrec3 < huge(nrec3)) THEN
Claus Kleinwort's avatar
Claus Kleinwort committed
894
895
        WRITE(8,*) 'Record',nrec3, ' is first with error (rank deficit/NaN)'
    END IF
896
897
898
899
900
901
99  WRITE(8,*) ' '
    IF (iteren > mreqenf) THEN
        WRITE(8,*) 'In total 3 +',nloopn,' loops through the data files'
    ELSE
        WRITE(8,*) 'In total 2 +',nloopn,' loops through the data files'
    ENDIF   
Claus Kleinwort's avatar
Claus Kleinwort committed
902
903
904
905
906
907
    IF (mnrsit > 0) THEN
        WRITE(8,*) ' '
        WRITE(8,*) 'In total    ',mnrsit,' internal MINRES iterations'
    END IF

    WRITE(8,103) times(0),times(1),times(2),times(4),times(7),  &
908
        times(5),times(8),times(3),times(6)
909

Claus Kleinwort's avatar
Claus Kleinwort committed
910
911
    CALL etime(ta,rst)
    deltat=rst-rstp
912
    ntsec=nint(deltat,mpi)
Claus Kleinwort's avatar
Claus Kleinwort committed
913
    CALL sechms(deltat,nhour,minut,secnd)
914
    nsecnd=nint(secnd,mpi)  ! round
Claus Kleinwort's avatar
Claus Kleinwort committed
915
    WRITE(8,*) 'Total time =',ntsec,' seconds =',nhour,' h',minut,  &
916
        ' m',nsecnd,' seconds'
Claus Kleinwort's avatar
Claus Kleinwort committed
917
918
    CALL fdate(chdate)
    WRITE(8,*) 'end                                            ', chdate
919
    gbu=1.0E-9*REAL(maxwordsalloc*(BIT_SIZE(1_mpi)/8),mps)             ! GB used
Claus Kleinwort's avatar
Claus Kleinwort committed
920
921
922
923
924
925
926
927
928
    WRITE(8,*) ' '
    WRITE(8,105) gbu

    !     Rejects ----------------------------------------------------------

    IF(nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) /= 0) THEN
        WRITE(8,*) ' '
        WRITE(8,*) 'Data rejected in last iteration:   '
        WRITE(8,*) '   ',  &
929
930
            nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0)   ',  &
            nrejec(2), ' (huge)   ',nrejec(3),' (large)'
Claus Kleinwort's avatar
Claus Kleinwort committed
931
932
        WRITE(8,*) ' '
    END IF
933
    IF (icheck <= 0) CALL explfc(8)
Claus Kleinwort's avatar
Claus Kleinwort committed
934
935
936
937

    WRITE(*,*) ' '
    WRITE(*,*) '  <  Millepede II-P ending   ... ', chdate ! with exit code',ITEXIT,' >'
    WRITE(*,*) ' '
938
    gbu=1.0E-9*REAL(maxwordsalloc*(BIT_SIZE(1_mpi)/8),mps)             ! GB used
Claus Kleinwort's avatar
Claus Kleinwort committed
939
940
941
942
943
    WRITE(*,105) gbu
    WRITE(*,*) ' '

102 FORMAT(2X,i4,2X,3F10.5,2X,3F10.5)
103 FORMAT(' Times [in sec] for     text processing',f12.3/  &
944
945
946
947
948
        '                                  LOOP1',f12.3/  &
        '                                  LOOP2',f12.3/  &
        '   func. value                         ',f12.3,' *',f4.0/  &
        '   func. value, global matrix, solution',f12.3,' *',f4.0/  &
        '                           new solution',f12.3,' *',f4.0/)
Claus Kleinwort's avatar
Claus Kleinwort committed
949
950
951
105 FORMAT('      Peak dynamic memory allocation: ',f11.6,' GB')
END PROGRAM mptwo                              ! Mille

952
!> Error for single global parameter from \ref minresmodule::minres "MINRES".
Claus Kleinwort's avatar
Claus Kleinwort committed
953
954
955
956
957
958
959
960
!!
!! Calculate single row 'x_i' from inverse matrix by solving A*x_i=b
!! with b=0 except b_i=1.
!!
!! \param [in]  ivgbi index of variable parameter

SUBROUTINE solglo(ivgbi)
    USE mpmod
961
    USE minresModule, ONLY: minres
Claus Kleinwort's avatar
Claus Kleinwort committed
962
963

    IMPLICIT NONE
964
    REAL(mps) :: par
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
    REAL(mps) :: dpa
    REAL(mps) :: err
    REAL(mps) :: gcor2
    INTEGER(mpi) :: iph
    INTEGER(mpi) :: istop
    INTEGER(mpi) :: itgbi
    INTEGER(mpi) :: itgbl
    INTEGER(mpi) :: itn
    INTEGER(mpi) :: itnlim
    INTEGER(mpi) :: nout

    INTEGER(mpi), INTENT(IN)                      :: ivgbi

    REAL(mpd) :: shift
    REAL(mpd) :: rtol
    REAL(mpd) :: anorm
    REAL(mpd) :: acond
982
    REAL(mpd) :: arnorm
983
984
985
986
987
988
989
    REAL(mpd) :: rnorm
    REAL(mpd) :: ynorm
    REAL(mpd) :: gmati
    REAL(mpd) :: diag
    INTEGER(mpl) :: ijadd
    INTEGER(mpl) :: jk
    INTEGER(mpl) :: ii
Claus Kleinwort's avatar
Claus Kleinwort committed
990
991
992
993
994
995
996
997
998
999
1000
    LOGICAL :: checka
    EXTERNAL avprod, mcsolv, mvsolv
    SAVE
    DATA iph/0/
    !     ...
    IF(iph == 0) THEN
        iph=1
        WRITE(*,101)
    END IF
    itgbi=globalParVarToTotal(ivgbi)
    itgbl=globalParLabelIndex(1,itgbi)