Loading...
Searching...
No Matches
pede.f90 File Reference
Millepede II program, subroutines. More...
Go to the source code of this file.
Functions/Subroutines | |
| program | mptwo |
| Millepede II main program Pede. | |
| subroutine | solglo (ivgbi) |
| Error for single global parameter from MINRES. | |
| subroutine | solgloqlp (ivgbi) |
| Error for single global parameter from MINRES-QLP. | |
| subroutine | addcst |
| Add constraint information to matrix and vector. | |
| subroutine | prpcon |
| Prepare constraints. | |
| subroutine | feasma |
| Matrix for feasible solution. | |
| subroutine | feasib (concut, iact) |
| Make parameters feasible. | |
| subroutine | peread (more) |
| Read (block of) records from binary files. | |
| subroutine | peprep (mode) |
| Prepare records. | |
| subroutine | pechk (ibuf, nerr) |
| Check Millepede record. | |
| subroutine | isjajb (nst, is, ja, jb, jsp) |
| Decode Millepede record. | |
| subroutine | loopn |
| Loop with fits and sums. | |
| subroutine | ploopa (lunp) |
| Print title for iteration. | |
| subroutine | ploopb (lunp) |
| Print iteration line. | |
| subroutine | ploopc (lunp) |
| Print sub-iteration line. | |
| subroutine | ploopd (lunp) |
| Print solution line. | |
| subroutine | explfc (lunit) |
| Print explanation of iteration table. | |
| subroutine | mupdat (i, j, add) |
| Update element of global matrix. | |
| subroutine | loopbf (nrej, ndfs, sndf, dchi2s, numfil, naccf, chi2f, ndff) |
| Loop over records in read buffer (block), fits and sums. | |
| subroutine | prtglo |
| Print final log file. | |
| subroutine | prtstat |
| Print input statistic. | |
| subroutine | avprd0 (n, x, b) |
| Product symmetric matrix times vector. | |
| subroutine | avprod (n, x, b) |
| Product symmetric matrix times vector. | |
| integer(mpl) function | ijadd (itema, itemb) |
| Index for sparse storage. | |
| subroutine | mhalf2 |
| Fill 2nd half of matrix for extended storage. | |
| subroutine | sechms (deltat, nhour, minut, secnd) |
| Time conversion. | |
| integer(mpi) function | inone (item) |
| Translate labels to indices (for global parameters). | |
| subroutine | upone |
| Update, redefine hash indices. | |
| integer(mpi) function | iprime (n) |
| largest prime number < N. | |
| subroutine | loop1 |
| First data loop (get global labels). | |
| subroutine | loop1i |
| Iteration of first data loop. | |
| subroutine | loop2 |
| Second data loop (number of derivatives, global label pairs). | |
| subroutine | monres |
| Monitor input residuals. | |
| subroutine | vmprep (msize) |
| Prepare storage for vectors and matrices. | |
| subroutine | minver |
| Solution by matrix inversion. | |
| subroutine | mdiags |
| Solution by diagonalization. | |
| subroutine | zdiags |
| Covariance matrix for diagonalization (,correction of eigenvectors). | |
| subroutine | mminrs |
| Solution with MINRES. | |
| subroutine | mminrsqlp |
| Solution with MINRES-QLP. | |
| subroutine | mcsolv (n, x, y) |
| Solution for zero band width preconditioner. | |
| subroutine | mvsolv (n, x, y) |
| Solution for finite band width preconditioner. | |
| subroutine | xloopn |
| Standard solution algorithm. | |
| subroutine | filetc |
| Interprete command line option, steering file. | |
| subroutine | filetx |
| Interprete text files. | |
| integer(mpi) function | nufile (fname) |
| Inquire on file. | |
| subroutine | intext (text, nline) |
| Interprete text. | |
| subroutine | additem (length, list, label, value) |
| add item to list | |
| subroutine | mstart (text) |
| Start of 'module' printout. | |
| subroutine | mend |
| End of 'module' printout. | |
| subroutine | mvopen (lun, fname) |
| Open file. | |
| subroutine | petime |
| Print times. | |
| subroutine | peend (icode, cmessage) |
| Print exit code. | |
| subroutine | binopn (kfile, ithr, ierr) |
| Open binary file. | |
| subroutine | bincls (kfile, ithr) |
| Close binary file. | |
| subroutine | binrwd (kfile) |
| Rewind binary file. | |
| subroutine | addsum (add) |
| Accurate summation. | |
| subroutine | getsum (asum) |
| Get accurate sum. | |
Detailed Description
Millepede II program, subroutines.
- Copyright
- Copyright (c) 2009 - 2018 Deutsches Elektronen-Synchroton, Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY
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.
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.
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.
Definition in file pede.f90.
Function/Subroutine Documentation
◆ addcst()
| subroutine addcst |
Add constraint information to matrix and vector.
Definition at line 1169 of file pede.f90.
1171
1172 IMPLICIT NONE
1173 REAL(mpd) :: climit
1174 REAL(mpd) :: factr
1175 REAL(mpd) :: sgm
1176
1177 INTEGER(mpi) :: i
1178 INTEGER(mpi) :: icgb
1179 INTEGER(mpi) :: irhs
1180 INTEGER(mpi) :: itgbi
1181 INTEGER(mpi) :: ivgb
1182 INTEGER(mpi) :: j
1183 INTEGER(mpi) :: jcgb
1184 INTEGER(mpi) :: l
1185 INTEGER(mpi) :: label
1186 INTEGER(mpi) :: nop
1187 INTEGER(mpi) :: inone
1188
1189 REAL(mpd) :: rhs
1190 REAL(mpd) :: drhs(4)
1191 INTEGER(mpi) :: idrh (4)
1192 SAVE
1193 ! ...
1194 nop=0
1196 climit=1.0e-5 ! limit for printout
1197 irhs=0 ! number of values in DRHS(.), to be printed
1198
1201 i=vecconsstart(icgb)
1205 label=listconstraints(j)%label
1206 factr=listconstraints(j)%value
1209
1212 END IF
1213
1215 END DO
1216 IF(abs(rhs) > climit) THEN
1217 irhs=irhs+1
1218 idrh(irhs)=jcgb
1219 drhs(irhs)=rhs
1220 nop=1
1221 IF(irhs == 4) THEN
1222 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1223 irhs=0
1224 END IF
1225 END IF
1226 vecconsresiduals(jcgb)=rhs
1228 END DO
1229
1230 IF(irhs /= 0) THEN
1231 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1232 END IF
1233 IF(nop == 0) RETURN
1234 WRITE(*,102) ' Constraints: only equation values >', climit,' are printed'
1235101 FORMAT(' ',4(i4,g11.3))
1236102 FORMAT(a,g11.2,a)
real(mpd), dimension(:), allocatable globalparameter
global parameters (start values + sum(x_i))
Definition mpmod.f90:174
real(mpd), dimension(:), allocatable vecconsresiduals
residuals of constraints
Definition mpmod.f90:203
integer(mpi), dimension(:,:), allocatable matconssort
keys and index for sorting
Definition mpmod.f90:207
type(listitem), dimension(:), allocatable listconstraints
list of constraints from steering file
Definition mpmod.f90:276
integer(mpi), dimension(:,:), allocatable globalparlabelindex
global parameters label, total -> var. index
Definition mpmod.f90:214
integer(mpi), dimension(:), allocatable vecconsstart
start of constraint in listConstraints (unsorted input)
Definition mpmod.f90:206
real(mpd), dimension(:), allocatable globalvector
global vector 'x' (in A*x=b)
Definition mpmod.f90:185
integer(mpi) function inone(item)
Translate labels to indices (for global parameters).
Definition pede.f90:4663
◆ additem()
| subroutine additem | ( | integer(mpi), intent(inout) | length, |
| type(listitem), dimension(:), intent(inout), allocatable | list, | ||
| integer(mpi), intent(in) | label, | ||
| real(mpd), intent(in) | value | ||
| ) |
add item to list
- Parameters
-
[in,out] length length of list [in,out] list list of items [in] label item label [in] value item value
Definition at line 8713 of file pede.f90.
8716
8717 INTEGER(mpi), INTENT(IN OUT) :: length
8718 TYPE(listItem), DIMENSION(:), INTENT(IN OUT), ALLOCATABLE :: list
8719 INTEGER(mpi), INTENT(IN) :: label
8720 REAL(mpd), INTENT(IN) :: value
8721
8722 INTEGER(mpl) :: newSize
8723 INTEGER(mpl) :: oldSize
8724 TYPE(listItem), DIMENSION(:), ALLOCATABLE :: tempList
8725
8726 IF (label > 0.AND.value == 0.) RETURN ! skip zero for valid labels
8727 IF (length == 0 ) THEN ! initial list with size = 100
8728 newsize = 100
8730 ENDIF
8732 IF (length >= oldsize) THEN ! increase sizeby 20% + 100
8733 newsize = oldsize + oldsize/5 + 100
8735 templist=list
8738 list(1:oldsize)=templist(1:oldsize)
8740 ENDIF
8741 ! add to end of list
8742 length=length+1
8743 list(length)%label=label
8744 list(length)%value=value
8745
◆ addsum()
| subroutine addsum | ( | real(mpd) | add | ) |
◆ avprd0()
| subroutine avprd0 | ( | integer(mpi), intent(in) | n, |
| real(mpd), dimension(n), intent(in) | x, | ||
| real(mpd), dimension(n), intent(out) | b | ||
| ) |
Product symmetric matrix times vector.
A(sym) * X => B. Used by MINRES method (Is most CPU intensive part). The matrix A is the global matrix in full symmetric or (compressed) sparse storage.
- Parameters
-
[in] n size of matrix [in] x vector X [in] b result vector B
Definition at line 4226 of file pede.f90.
4228
4229 IMPLICIT NONE
4230 INTEGER(mpi) :: i
4231 INTEGER(mpi) :: iencdb
4232 INTEGER(mpi) :: iencdm
4233 INTEGER(mpi) :: iproc
4234 INTEGER(mpi) :: ir
4235 INTEGER(mpi) :: j
4236 INTEGER(mpi) :: jc
4237 INTEGER(mpi) :: jj
4238 INTEGER(mpi) :: jn
4239
4240 INTEGER(mpi), INTENT(IN) :: n
4241 REAL(mpd), INTENT(IN) :: x(n)
4242 REAL(mpd), INTENT(OUT) :: b(n)
4243 INTEGER(mpl) :: k
4244 INTEGER(mpl) :: kk
4245 INTEGER(mpl) :: kl
4246 INTEGER(mpl) :: ku
4247 INTEGER(mpl) :: ll
4248 INTEGER(mpl) :: lj
4249 INTEGER(mpl) :: indij
4250 INTEGER(mpl) :: indid
4251 INTEGER(mpl) :: ij
4252 INTEGER(mpi) :: ichunk
4253 !$ INTEGER(mpi) OMP_GET_THREAD_NUM
4254 SAVE
4255 ! ...
4256 !$ DO i=1,n
4257 !$ b(i)=0.0_mpd ! reset 'global' B()
4258 !$ END DO
4261 ! full symmetric matrix
4262 ! parallelize row loop
4263 ! private copy of B(N) for each thread, combined at end, init with 0.
4264 ! slot of 1024 'I' for next idle thread
4265 !$OMP PARALLEL DO &
4266 !$OMP PRIVATE(J,IJ) &
4267 !$OMP REDUCTION(+:B) &
4268 !$OMP SCHEDULE(DYNAMIC,ichunk)
4269 DO i=1,n
4270 ij=i
4271 ij=(ij*ij-ij)/2
4272 b(i)=globalmatd(ij+i)*x(i)
4273 DO j=1,i-1
4274 b(j)=b(j)+globalmatd(ij+j)*x(i)
4275 b(i)=b(i)+globalmatd(ij+j)*x(j)
4276 END DO
4277 END DO
4278 !$OMP END PARALLEL DO
4279 ELSE
4280 ! sparse, compressed matrix
4283 stop 'AVPRD0: mismatched vector and matrix'
4284 END IF
4285 iencdb=nencdb
4286 iencdm=ishft(1,iencdb)-1
4287 ! parallelize row loop
4288 ! slot of 1024 'I' for next idle thread
4289 !$OMP PARALLEL DO &
4290 !$OMP PRIVATE(IR,K,KK,LL,KL,KU,INDID,INDIJ,J,JC,JN,LJ,JJ) &
4291 !$OMP REDUCTION(+:B) &
4292 !$OMP SCHEDULE(DYNAMIC,ichunk)
4293 DO i=1,n
4294 iproc=0
4295 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
4297 ! ! off-diagonals double precision
4298 ir=i
4301 kl=0
4302 ku=sparsematrixoffsets(1,ir+1)-1-kk
4303 indid=kk
4304 indij=ll
4306 DO k=kl,ku
4307 j=sparsematrixcolumns(indid+k)
4308 b(j)=b(j)+globalmatd(indij+k)*x(i)
4309 b(i)=b(i)+globalmatd(indij+k)*x(j)
4310 END DO
4311 ELSE
4312 lj=0
4313 ku=((ku+1)*8)/9-1 ! number of regions (-1)
4314 indid=indid+ku/8+1 ! skip group offsets
4316 ! extended storage
4317 DO kl=0,ku
4318 jc=sparsematrixcolumns(indid+kl)
4319 j=ishft(jc,-iencdb)
4320 jn=iand(jc, iencdm)
4321 b(i)=b(i)+dot_product(globalmatd(indij+lj:indij+lj+jn-1),x(j:j+jn-1))
4322 lj=lj+jn
4323 END DO
4324 ELSE
4325 DO kl=0,ku
4326 jc=sparsematrixcolumns(indid+kl)
4327 j=ishft(jc,-iencdb)
4328 jn=iand(jc, iencdm)
4329 DO jj=1,jn
4330 b(j)=b(j)+globalmatd(indij+lj)*x(i)
4331 b(i)=b(i)+globalmatd(indij+lj)*x(j)
4332 j=j+1
4333 lj=lj+1
4334 END DO
4335 END DO
4336 END IF
4337 END IF
4338
4340 ir=i+n+1 ! off-diagonals single precision
4343 kl=0
4344 ku=sparsematrixoffsets(1,ir+1)-1-kk
4345 indid=kk
4346 indij=ll
4348 DO k=kl,ku
4349 j=sparsematrixcolumns(indid+k)
4350 b(j)=b(j)+real(globalmatf(indij+k),mpd)*x(i)
4351 b(i)=b(i)+real(globalmatf(indij+k),mpd)*x(j)
4352 END DO
4353 ELSE
4354 lj=0
4355 ku=((ku+1)*8)/9-1 ! number of regions (-1)
4356 indid=indid+ku/8+1 ! skip group offsets
4357 DO kl=0,ku
4358 jc=sparsematrixcolumns(indid+kl)
4359 j=ishft(jc,-iencdb)
4360 jn=iand(jc, iencdm)
4362 ! extended storage
4363 DO jj=1,jn
4364 b(i)=b(i)+real(globalmatf(indij+lj),mpd)*x(j)
4365 j=j+1
4366 lj=lj+1
4367 END DO
4368 ELSE
4369 DO jj=1,jn
4370 b(j)=b(j)+real(globalmatf(indij+lj),mpd)*x(i)
4371 b(i)=b(i)+real(globalmatf(indij+lj),mpd)*x(j)
4372 j=j+1
4373 lj=lj+1
4374 END DO
4375 END IF
4376 END DO
4377 END IF
4378 END IF
4379 END DO
4380 ENDIF
4381
integer(mpi), dimension(:), allocatable sparsematrixcolumns
(compressed) list of columns for sparse matrix
Definition mpmod.f90:230
integer(mpl), dimension(:,:), allocatable sparsematrixoffsets
row offsets for column list, sparse matrix elements
Definition mpmod.f90:231
real(mpd), dimension(:), allocatable globalmatd
global matrix 'A' (double, full or sparse)
Definition mpmod.f90:183
real(mps), dimension(:), allocatable globalmatf
global matrix 'A' (float part for compressed sparse)
Definition mpmod.f90:184
integer(mpi) mextnd
flag for extended storage (both 'halves' of sym. mat. for improved access patterns)
Definition mpmod.f90:83
◆ avprod()
| subroutine avprod | ( | integer(mpi), intent(in) | n, |
| real(mpd), dimension(n), intent(in) | x, | ||
| real(mpd), dimension(n), intent(out) | b | ||
| ) |
Product symmetric matrix times vector.
A(sym) * X => B. Used by MINRES method (Is most CPU intensive part). The matrix A is the global matrix in full symmetric or (compressed) sparse storage. Allows for size of X and smaller than size of matrix in case of solution with constriants by elimination.
- Parameters
-
[in] n size of matrix ( <= size of global matrix) [in] x vector X [in] b result vector B
Definition at line 4394 of file pede.f90.
4396
4397 IMPLICIT NONE
4398
4399 INTEGER(mpi), INTENT(IN) :: n
4400 REAL(mpd), INTENT(IN) :: x(n)
4401 REAL(mpd), INTENT(OUT) :: b(n)
4402
4403 SAVE
4404 ! ...
4407 stop 'AVPROD: mismatched vector and matrix'
4408 END IF
4409 ! input to AVPRD0
4410 vecxav(1:n)=x
4412 !use elimination for constraints ?
4414 ! calclulate vecBav=globalMat*vecXav
4416 !use elimination for constraints ?
4418 ! output from AVPRD0
4419 b=vecbav(1:n)
4420
◆ bincls()
| subroutine bincls | ( | integer(mpi), intent(in) | kfile, |
| integer(mpi), intent(in) | ithr | ||
| ) |
Close binary file.
- Parameters
-
[in] kfile file number [in] ithr thread number ([1..nthrd] - close and reopen) for C files
Definition at line 8987 of file pede.f90.
8989
8990 IMPLICIT NONE
8991 INTEGER(mpi), INTENT(IN) :: kfile
8992 INTEGER(mpi), INTENT(IN) :: ithr
8993
8994 INTEGER(mpi) :: lun
8995
8996 lun=ithr
8997 !print *, " closing binary ", kfile, ithr
8999 lun=kfile+10
9000 CLOSE(lun)
9001#ifdef READ_C_FILES
9002 ELSE ! C file
9003 CALL closec(lun)
9004#endif
9005 END IF
9006
◆ binopn()
| subroutine binopn | ( | integer(mpi), intent(in) | kfile, |
| integer(mpi), intent(in) | ithr, | ||
| integer(mpi), intent(out) | ierr | ||
| ) |
Open binary file.
- Parameters
-
[in] kfile file number [in] ithr thread number ([1..nthrd] - close and reopen) or 0 (next file - keep open) for C files [out] ierr error flag
Definition at line 8909 of file pede.f90.
8911
8912 IMPLICIT NONE
8913 INTEGER(mpi), INTENT(IN) :: kfile
8914 INTEGER(mpi), INTENT(IN) :: ithr
8915 INTEGER(mpi), INTENT(OUT) :: ierr
8916
8917 INTEGER(mpi), DIMENSION(13) :: ibuff
8918 INTEGER(mpi) :: ioff
8919 INTEGER(mpi) :: ios
8920 INTEGER(mpi) :: k
8921 INTEGER(mpi) :: lfn
8922 INTEGER(mpi) :: lun
8923 INTEGER(mpi) :: moddate
8924 CHARACTER (LEN=1024) :: fname
8925
8926 ierr=0
8927 lun=ithr
8928 ! modification date (=0: open for first time, >0: reopen, <0: unknown )
8929 moddate=yfd(kfile)
8930 ! file name
8931 ioff=sfd(1,kfile)
8932 lfn=sfd(2,kfile)
8934 !print *, " opening binary ", kfile, ithr, moddate, " : ", fname(1:lfn)
8935 ! open
8936 ios=0
8938 ! Fortran file
8939 lun=kfile+10
8940 OPEN(lun,file=fname(1:lfn),iostat=ios, form='UNFORMATTED')
8941 print *, ' lun ', lun, ios
8942#ifdef READ_C_FILES
8943 ELSE
8944 ! C file
8945 CALL openc(fname(1:lfn),lun,ios)
8946#else
8947 WRITE(*,*) 'Opening of C-files not supported.'
8948 ierr=1
8949 RETURN
8950#endif
8951 END IF
8952 IF(ios /= 0) THEN
8953 ierr=1
8954 WRITE(*,*) 'Open error for file ',fname(1:lfn), ios
8955 IF (moddate /= 0) THEN
8957 stop 'PEREAD: open error '
8958 ENDIF
8959 RETURN
8960 END IF
8961 ! get status
8962 CALL stat(fname(1:lfn),ibuff,ios)
8963 !print *, ' STAT ', ios, ibuff(10), moddate
8964 IF(ios /= 0) THEN
8965 ierr=1
8966 WRITE(*,*) 'STAT error for file ',fname(1:lfn), ios
8967 ibuff(10)=-1
8968 END IF
8969 ! check/store modification date
8970 IF (moddate /= 0) THEN
8971 IF (ibuff(10) /= moddate) THEN
8973 stop 'PEREAD: file modified '
8974 END IF
8975 ELSE
8976 yfd(kfile)=ibuff(10)
8977 END IF
8978 RETURN
8979
integer(mpi), dimension(:,:), allocatable sfd
offset (1,..), length (2,..) of binary file name in tfd
Definition mpmod.f90:293
◆ binrwd()
| subroutine binrwd | ( | integer(mpi), intent(in) | kfile | ) |
Rewind binary file.
- Parameters
-
[in] kfile file number
Definition at line 9013 of file pede.f90.
9015
9016 IMPLICIT NONE
9017 INTEGER(mpi), INTENT(IN) :: kfile
9018
9019 INTEGER(mpi) :: lun
9020
9021 !print *, " rewinding binary ", kfile
9023 lun=kfile+10
9024 rewind lun
9025#ifdef READ_C_FILES
9026 ELSE
9027 lun=kfile-nfilf
9028 CALL resetc(lun)
9029#endif
9030 END IF
9031
◆ explfc()
| subroutine explfc | ( | integer(mpi) | lunit | ) |
Print explanation of iteration table.
Definition at line 2886 of file pede.f90.
2889
2890 IMPLICIT NONE
2891 INTEGER(mpi) :: lunit
2892 WRITE(lunit,*) ' '
2893 WRITE(lunit,102) 'Explanation of iteration table'
2894 WRITE(lunit,102) '=============================='
2895 WRITE(lunit,101) 'it', &
2896 'iteration number. Global parameters are improved for it > 0.'
2897 WRITE(lunit,102) 'First function evaluation is called iteraton 0.'
2898 WRITE(lunit,101) 'fc', 'number of function evaluations.'
2899 WRITE(lunit,101) 'fcn_value', 'value of 2 x Likelihood function (LF).'
2900 WRITE(lunit,102) 'The final value is the chi^2 value of the fit and should'
2901 WRITE(lunit,102) 'be about equal to the NDF (see below).'
2902 WRITE(lunit,101) 'dfcn_exp', &
2903 'expected reduction of the value of the Likelihood function (LF)'
2904 WRITE(lunit,101) 'slpr', 'ratio of the actual slope to inital slope.'
2905 WRITE(lunit,101) 'costh', &
2906 'cosine of angle between search direction and -gradient'
2908 WRITE(lunit,101) 'iit', &
2909 'number of internal iterations in MINRES algorithm'
2910 WRITE(lunit,101) 'st', 'stop code of MINRES algorithm'
2911 WRITE(lunit,102) '< 0: rhs is very special, with beta2 = 0'
2912 WRITE(lunit,102) '= 0: rhs b = 0, i.e. the exact solution is x = 0'
2913 WRITE(lunit,102) '= 1 requested accuracy achieved, as determined by rtol'
2914 WRITE(lunit,102) '= 2 reasonable accuracy achieved, given eps'
2915 WRITE(lunit,102) '= 3 x has converged to an eigenvector'
2916 WRITE(lunit,102) '= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
2917 WRITE(lunit,102) '= 5 the iteration limit was reached'
2918 WRITE(lunit,102) '= 6 Matrix x vector does not define a symmetric matrix'
2919 WRITE(lunit,102) '= 7 Preconditioner does not define a symmetric matrix'
2921 WRITE(lunit,101) 'iit', &
2922 'number of internal iterations in MINRES-QLP algorithm'
2923 WRITE(lunit,101) 'st', 'stop code of MINRES-QLP algorithm'
2924 WRITE(lunit,102) '= 1: beta_{k+1} < eps, iteration k is the final Lanczos step.'
2925 WRITE(lunit,102) '= 2: beta2 = 0. If M = I, b and x are eigenvectors of A.'
2926 WRITE(lunit,102) '= 3: beta1 = 0. The exact solution is x = 0.'
2927 WRITE(lunit,102) '= 4: A solution to (poss. singular) Ax = b found, given rtol.'
2928 WRITE(lunit,102) '= 5: A solution to (poss. singular) Ax = b found, given eps.'
2929 WRITE(lunit,102) '= 6: Pseudoinverse solution for singular LS problem, given rtol.'
2930 WRITE(lunit,102) '= 7: Pseudoinverse solution for singular LS problem, given eps.'
2931 WRITE(lunit,102) '= 8: The iteration limit was reached.'
2932 WRITE(lunit,102) '= 9: The operator defined by Aprod appears to be unsymmetric.'
2933 WRITE(lunit,102) '=10: The operator defined by Msolve appears to be unsymmetric.'
2934 WRITE(lunit,102) '=11: The operator defined by Msolve appears to be indefinite.'
2935 WRITE(lunit,102) '=12: xnorm has exceeded maxxnorm or will exceed it next iteration.'
2936 WRITE(lunit,102) '=13: Acond has exceeded Acondlim or 0.1/eps.'
2937 WRITE(lunit,102) '=14: Least-squares problem but no converged solution yet.'
2938 WRITE(lunit,102) '=15: A null vector obtained, given rtol.'
2939 ENDIF
2940 WRITE(lunit,101) 'ls', 'line search info'
2941 WRITE(lunit,102) '< 0 recalculate function'
2942 WRITE(lunit,102) '= 0: N or STP lt 0 or step not descending'
2943 WRITE(lunit,102) '= 1: Linesearch convergence conditions reached'
2944 WRITE(lunit,102) '= 2: interval of uncertainty at lower limit'
2945 WRITE(lunit,102) '= 3: max nr of line search calls reached'
2946 WRITE(lunit,102) '= 4: step at the lower bound'
2947 WRITE(lunit,102) '= 5: step at the upper bound'
2948 WRITE(lunit,102) '= 6: rounding error limitation'
2949 WRITE(lunit,101) 'step', &
2950 'the factor for the Newton step during the line search. Usually'
2951 WRITE(lunit,102) &
2952 'a value of 1 gives a sufficient reduction of the LF. Oherwise'
2953 WRITE(lunit,102) 'other step values are tried.'
2954 WRITE(lunit,101) 'cutf', &
2955 'cut factor. Local fits are rejected, if their chi^2 value'
2956 WRITE(lunit,102) &
2957 'is larger than the 3-sigma chi^2 value times the cut factor.'
2958 WRITE(lunit,102) 'A cut factor of 1 is used finally, but initially a larger'
2959 WRITE(lunit,102) 'factor may be used. A value of 0.0 means no cut.'
2960 WRITE(lunit,101) 'rejects', 'total number of rejected local fits.'
2961 WRITE(lunit,101) 'hmmsec', 'the time in hours (h), minutes (mm) and seconds.'
2962 WRITE(lunit,101) 'FMS', 'calculation of Function value, Matrix, Solution.'
2963 WRITE(lunit,*) ' '
2964
2965101 FORMAT(a9,' = ',a)
2966102 FORMAT(13x,a)
integer(mpi) metsol
solution method (1: inversion, 2: diagonalization, 3: MINRES-QLP)
Definition mpmod.f90:34
◆ feasib()
| subroutine feasib | ( | real(mps), intent(in) | concut, |
| integer(mpi), intent(out) | iact | ||
| ) |
Make parameters feasible.
Correct for constraint equation discrepancies.
- Parameters
-
[in] concut cut for discrepancies [out] iact =1 if correction needed, else =0
Definition at line 1576 of file pede.f90.
1579
1580 IMPLICIT NONE
1581 REAL(mpd) :: factr
1582 REAL(mpd) :: sgm
1583 INTEGER(mpi) :: i
1584 INTEGER(mpi) :: icgb
1585 INTEGER(mpi) :: iter
1586 INTEGER(mpi) :: itgbi
1587 INTEGER(mpi) :: ivgb
1588 INTEGER(mpi) :: iblck
1589 INTEGER(mpi) :: ieblck
1590 INTEGER(mpi) :: isblck
1591 INTEGER(mpi) :: ifirst
1592 INTEGER(mpi) :: ilast
1593 INTEGER(mpi) :: j
1594 INTEGER(mpi) :: jcgb
1595 INTEGER(mpi) :: label
1596 INTEGER(mpi) :: inone
1597 INTEGER(mpi) :: ncon
1598
1599 REAL(mps), INTENT(IN) :: concut
1600 INTEGER(mpi), INTENT(OUT) :: iact
1601
1602 REAL(mpd) :: rhs
1603 REAL(mpd) ::sum1
1604 REAL(mpd) ::sum2
1605 REAL(mpd) ::sum3
1606
1607 REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecCorrections
1608 SAVE
1609
1610 iact=0
1612
1613 DO iter=1,2
1614 vecconsresiduals=0.0_mpd
1615
1616 ! calculate right constraint equation discrepancies
1619 i=vecconsstart(icgb)
1623 label=listconstraints(j)%label
1624 factr=listconstraints(j)%value
1627 ENDDO
1629 END DO
1630
1631 ! constraint equation discrepancies -------------------------------
1632
1633 sum1=0.0_mpd
1634 sum2=0.0_mpd
1635 sum3=0.0_mpd
1637 sum1=sum1+vecconsresiduals(icgb)**2
1638 sum2=sum2+abs(vecconsresiduals(icgb))
1639 sum3=max(sum3,abs(vecconsresiduals(icgb)))
1640 END DO
1641 sum1=sqrt(sum1/real(ncgb,mpd))
1642 sum2=sum2/real(ncgb,mpd)
1643
1644 IF(iter == 1.AND.sum1 < concut) RETURN ! do nothing if correction small
1645
1647 WRITE(*,*) ' '
1648 WRITE(*,*) 'Constraint equation discrepancies:'
1650101 FORMAT(4x,4(i5,g12.4))
1651 WRITE(*,103) concut
1652103 FORMAT(10x,' Cut on rms value is',g8.1)
1653 END IF
1654
1655 IF(iact == 0) THEN
1656 WRITE(*,*) ' '
1657 WRITE(*,*) 'Improve constraints'
1658 END IF
1659 iact=1
1660
1661 WRITE(*,102) iter,sum1,sum2,sum3
1662102 FORMAT(i6,' rms',g12.4,' avrg_abs',g12.4,' max_abs',g12.4)
1663
1665 veccorrections=0.0_mpd
1666
1667 ! multiply (block-wise) inverse matrix and constraint vector
1668 isblck=0
1672 ncon=ilast+1-ifirst
1673 ieblck=isblck+(ncon*(ncon+1))/2
1674 CALL dbsvx(matconsproduct(isblck+1:ieblck),vecconsresiduals(ifirst:ilast),vecconssolution(ifirst:ilast),ncon)
1675 isblck=ieblck
1676 END DO
1677
1680 i=vecconsstart(icgb)
1684 label=listconstraints(j)%label
1685 factr=listconstraints(j)%value
1688 IF(ivgb > 0) THEN
1689 veccorrections(ivgb)=veccorrections(ivgb)+vecconssolution(jcgb)*factr
1690 END IF
1691 ENDDO
1692 END DO
1693
1695 itgbi=globalparvartototal(i)
1697 END DO
1698
1700
1701 END DO ! iteration 1 and 2
1702
real(mpd), dimension(:), allocatable matconsproduct
product matrix of constraints
Definition mpmod.f90:202
real(mpd), dimension(:), allocatable vecconssolution
solution for constraint elimination
Definition mpmod.f90:204
integer(mpi), dimension(:), allocatable globalparvartototal
global parameters variable -> total index
Definition mpmod.f90:216
integer(mpi), dimension(:,:), allocatable matconsblocks
start of constraint blocks, parameter range
Definition mpmod.f90:208
◆ feasma()
| subroutine feasma |
Matrix for feasible solution.
Check rank of product matrix of constraints.
Definition at line 1427 of file pede.f90.
1430
1431 IMPLICIT NONE
1432 REAL(mpd) :: factr
1433 REAL(mpd) :: sgm
1434 INTEGER(mpi) :: i
1435 INTEGER(mpi) :: iblck
1436 INTEGER(mpi) :: icgb
1437 INTEGER(mpi) :: ij
1438 INTEGER(mpi) :: ifirst
1439 INTEGER(mpi) :: ilast
1440 INTEGER(mpi) :: ioffc
1441 INTEGER(mpi) :: ioffp
1442 INTEGER(mpi) :: irank
1443 INTEGER(mpi) :: ipar0
1444 INTEGER(mpi) :: itgbi
1445 INTEGER(mpi) :: ivgb
1446 INTEGER(mpi) :: j
1447 INTEGER(mpi) :: jcgb
1448 INTEGER(mpi) :: l
1449 INTEGER(mpi) :: label
1450 INTEGER(mpi) :: ncon
1451 INTEGER(mpi) :: npar
1452 INTEGER(mpi) :: nrank
1453 INTEGER(mpi) :: inone
1454
1455 REAL(mpd):: rhs
1456 REAL(mpd):: evmax
1457 REAL(mpd):: evmin
1458 INTEGER(mpl):: length
1459 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matConstraintsT
1460 REAL(mpd), DIMENSION(:), ALLOCATABLE :: auxVectorD
1461 INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: auxVectorI
1462 SAVE
1463 ! ...
1464
1466
1467 ! QL decomposition
1469 ! product matrix A A^T (A is stored as transposed)
1470 length=mszprd
1472 matconsproduct=0.0_mpd
1473 length=ncgb
1478 ! constraint matrix A (A is stored as transposed)
1479 length = mszcon
1481 matconstraintst=0.0_mpd
1482
1483 ! loop over sorted constraints, fill matrices, get rank, inverted product matrix (in blocks)
1484 ioffc=0 ! block offset in constraint matrix
1485 ioffp=0 ! block offset in product matrix
1486 nrank=0
1490 ncon=ilast+1-ifirst
1493 DO jcgb=ifirst,ilast
1494 ! index in list
1495 icgb=matconssort(3,jcgb)
1496 ! fill constraint matrix
1497 i=vecconsstart(icgb)
1501 label=listconstraints(j)%label
1502 factr=listconstraints(j)%value
1505 IF(ivgb > 0) matconstraintst(ivgb-ipar0+ioffc+(jcgb-ifirst)*npar)=factr ! matrix element
1508 END DO
1510 END DO
1511
1512 ! get rank of blocks
1513 DO l=ioffc+1,ioffc+npar
1514 ij=ioffp
1515 DO i=1,ncon
1516 DO j=1,i
1517 ij=ij+1
1518 matconsproduct(ij)=matconsproduct(ij)+matconstraintst((i-1)*npar+l)*matconstraintst((j-1)*npar+l)
1519 END DO
1520 END DO
1521 END DO
1522 ! inversion of product matrix of constraints
1523 CALL sqminv(matconsproduct(ioffp+1:ij),vecconsresiduals(ifirst:ilast),ncon,irank, auxvectord, auxvectori)
1525 nrank=nrank+irank
1526 ioffc=ioffc+npar*ncon
1527 ioffp=ij
1528 END DO
1529
1531
1532 WRITE(*,*) ' '
1533 WRITE(*,*) 'Rank of product matrix of constraints is',nrank, &
1535 WRITE(8,*) 'Rank of product matrix of constraints is',nrank, &
1538 WRITE(*,*) 'Warning: insufficient constraint equations!'
1539 WRITE(8,*) 'Warning: insufficient constraint equations!'
1541 isubit=1
1542 WRITE(*,*) ' --> enforcing SUBITO mode'
1543 WRITE(8,*) ' --> enforcing SUBITO mode'
1544 END IF
1545 END IF
1546
1547 ! QL decomposition
1549 print *
1550 print *, 'QL decomposition of constraints matrix'
1552 ! check eignevalues of L
1554 print *, ' largest |eigenvalue| of L: ', evmax
1555 print *, ' smallest |eigenvalue| of L: ', evmin
1556 IF (evmin == 0.0_mpd) THEN
1558 stop 'FEASMA: stopping due to singular QL decomposition of constraints matrix'
1559 END IF
1560 END IF
1561
1565
1566 RETURN
integer(mpi) mszprd
(integrated block) matrix size for (constraint) product matrix
Definition mpmod.f90:125
integer(mpi), dimension(:), allocatable globalparcons
global parameters (number of) constraints
Definition mpmod.f90:181
◆ filetc()
| subroutine filetc |
Interprete command line option, steering file.
Fetch and interprete command line options, if steering file specified, check file existence (calling NUFILE)
If no steering file specified, check default steering file.
Create test files for command line option '-t'.
Read steering file, print some lines, detect names of text and binary files, check file existence, store all file names.
Open all binary files.
Definition at line 7283 of file pede.f90.
7286
7287 IMPLICIT NONE
7288 INTEGER(mpi) :: i
7289 INTEGER(mpi) :: ia
7290 INTEGER(mpi) :: iargc
7291 INTEGER(mpi) :: ib
7292 INTEGER(mpi) :: ie
7293 INTEGER(mpi) :: ierrf
7294 INTEGER(mpi) :: ieq
7295 INTEGER(mpi) :: ifilb
7296 INTEGER(mpi) :: ioff
7297 INTEGER(mpi) :: iopt
7298 INTEGER(mpi) :: ios
7299 INTEGER(mpi) :: iosum
7300 INTEGER(mpi) :: it
7301 INTEGER(mpi) :: k
7302 INTEGER(mpi) :: mat
7303 INTEGER(mpi) :: nab
7304 INTEGER(mpi) :: nline
7305 INTEGER(mpi) :: npat
7306 INTEGER(mpi) :: ntext
7307 INTEGER(mpi) :: nu
7308 INTEGER(mpi) :: nuf
7309 INTEGER(mpi) :: nums
7310 INTEGER(mpi) :: nufile
7311 INTEGER(mpi) :: lenfileInfo
7312 INTEGER(mpi) :: lenFileNames
7313 INTEGER(mpi) :: matint
7314 INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: vecfileInfo
7315 INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: tempArray
7316 INTEGER(mpl) :: rows
7317 INTEGER(mpl) :: cols
7318 INTEGER(mpl) :: newcols
7319 INTEGER(mpl) :: length
7320
7321 CHARACTER (LEN=1024) :: text
7322 CHARACTER (LEN=1024) :: fname
7323 CHARACTER (LEN=14) :: bite(3)
7324 CHARACTER (LEN=32) :: keystx
7325 REAL(mpd) :: dnum(100)
7326 SAVE
7327 DATA bite/'C_binary','text ','Fortran_binary'/
7328 ! ...
7330
7331 nuf=1 ! C binary is default
7332 DO i=1,8
7333 times(i)=0.0
7334 END DO
7335
7336 ! read command line options ----------------------------------------
7337
7339 DO i=1,iargc()
7340 IF(i == 1) THEN
7341 WRITE(*,*) ' '
7342 WRITE(*,*) 'Command line options: '
7343 WRITE(*,*) '--------------------- '
7344 END IF
7345 CALL getarg(i,text) ! get I.th text from command line
7347 WRITE(*,101) i,text(1:nab) ! echo print
7348 IF(text(ia:ia) /= '-') THEN
7350 IF(nu == 2) THEN ! existing text file
7352 WRITE(*,*) 'Second text file in command line - stop'
7354 stop
7355 ELSE
7356 filnam=text
7357 END IF
7358 ELSE
7359 WRITE(*,*) 'Open error for file:',text(ia:ib),' - stop'
7361 stop
7362 END IF
7363 ELSE
7364 IF(index(text(ia:ib),'b') /= 0) THEN
7366 WRITE(*,*) 'Debugging requested'
7367 END IF
7368 it=index(text(ia:ib),'t')
7369 IF(it /= 0) THEN
7371 ieq=index(text(ia+it:ib),'=')+it
7372 IF (it /= ieq) THEN
7378 END IF
7379 END IF
7384 END IF
7385 IF(i == iargc()) WRITE(*,*) '--------------------- '
7386 END DO
7387
7388
7389 ! create test files for option -t ----------------------------------
7390
7392 WRITE(*,*) ' '
7395 ELSE
7397 END IF
7399 WRITE(*,*) ' '
7400 END IF
7401
7402 ! check default steering file with file-name "steerfile" -----------
7403
7405 text='steerfile'
7408 IF(nu > 0) THEN
7409 filnam=text
7410 ELSE
7412 stop 'in FILETC: no steering file. .'
7413 END IF
7414 END IF
7415
7416
7417 ! open, read steering file:
7418 ! end
7419 ! fortranfiles
7420 ! cfiles
7421
7422
7424 WRITE(*,*) ' '
7426 WRITE(*,*) '-------------------------'
7428 IF(ios /= 0) THEN
7429 WRITE(*,*) 'Open error for steering file - stop'
7431 stop
7432 END IF
7433 ifile =0
7434 nfiles=0
7435
7436 lenfileinfo=2
7437 lenfilenames=0
7438 rows=6; cols=lenfileinfo
7440 nline=0
7441 DO
7442 READ(10,102,iostat=ierrf) text ! read steering file
7443 IF (ierrf < 0) EXIT ! eof
7445 nline=nline+1
7446 IF(nline <= 50) THEN ! print up to 50 lines
7447 WRITE(*,101) nline,text(1:nab)
7448 IF(nline == 50) WRITE(*,*) ' ...'
7449 END IF
7450
7452 IF(ib == ia+2) THEN
7454 IF(mat == 3) THEN
7455 text=' '
7457 WRITE(*,*) ' end-statement after',nline,' text lines'
7458 EXIT
7459 END IF
7460 END IF
7461
7462 keystx='fortranfiles'
7463 mat=matint(text(ia:ib),keystx,npat,ntext)
7464 IF(mat == ntext) THEN ! exact matching
7465 nuf=3
7466 ! WRITE(*,*) 'Fortran files'
7467 cycle
7468 END IF
7469
7470 keystx='Cfiles'
7471 mat=matint(text(ia:ib),keystx,npat,ntext)
7472 IF(mat == ntext) THEN ! exact matching
7473 nuf=1
7474 ! WRITE(*,*) 'Cfiles'
7475 cycle
7476 END IF
7477
7478 keystx='closeandreopen' ! don't keep binary files open
7479 mat=matint(text(ia:ib),keystx,npat,ntext)
7480 IF(mat == ntext) THEN ! exact matching
7481 keepopen=0
7482 cycle
7483 END IF
7484
7485 ! file names
7486 ! check for file options (' -- ')
7487 ie=ib
7488 iopt=index(text(ia:ib),' -- ')
7489 IF (iopt > 0) ie=iopt-1
7490
7491 IF(nab == 0) cycle
7493 IF(nu > 0) THEN ! existing file
7496 temparray=vecfileinfo
7498 lenfileinfo=lenfileinfo*2
7499 newcols=lenfileinfo
7501 vecfileinfo(:,1:cols)=temparray(:,1:cols)
7503 cols=newcols
7504 ENDIF
7506 IF(nu == 1) nu=nuf !
7507 lenfilenames=lenfilenames+ie-ia+1 ! total length of file names
7514 ELSE
7515 ! WRITE(*,*) 'Open error for file ',TEXT(IA:IB)
7516 ! STOP
7517 END IF
7518 END DO
7519 rewind 10
7520 ! read again to fill dynamic arrays with file info
7521 length=nfiles
7526 length=lenfilenames
7528 nline=0
7529 i=1
7530 ioff=0
7531 DO
7532 READ(10,102,iostat=ierrf) text ! read steering file
7533 IF (ierrf < 0) EXIT ! eof
7534 nline=nline+1
7535 IF (nline == vecfileinfo(1,i)) THEN
7536 nfd(i)=vecfileinfo(1,i)
7537 mfd(i)=vecfileinfo(2,i)
7538 ia=vecfileinfo(3,i)-1
7541 ! tfd(i)=text(vecFileInfo(3,i):vecFileInfo(4,i)) ! file name
7542 ioff=ioff+lfd(i)
7544 IF (vecfileinfo(5,i) > 0) THEN
7545 CALL ratext(text(vecfileinfo(5,i)+4:vecfileinfo(6,i)),nums,dnum) ! translate text to DP numbers
7547 END IF
7548 i=i+1
7550 ENDIF
7551 ENDDO
7553 rewind 10
7554 ! additional info for binary files
7555 length=nfiles; rows=2
7565 yfd=0
7566 !
7567 WRITE(*,*) '-------------------------'
7568 WRITE(*,*) ' '
7569
7570 ! print table of files ---------------------------------------------
7571
7573 WRITE(*,*) 'Table of files:'
7574 WRITE(*,*) '---------------'
7575 END IF
7576 WRITE(8,*) ' '
7577 WRITE(8,*) 'Text and data files:'
7578 ioff=0
7581 ! fname=tfd(i)(1:lfd(i))
7584 ioff=ioff+lfd(i)
7585 END DO
7587 WRITE(*,*) '---------------'
7588 WRITE(*,*) ' '
7589 END IF
7590
7591 ! open the binary Fortran (data) files on unit 11, 12, ...
7592
7593 iosum=0
7594 nfilf=0
7595 nfilb=0
7596 nfilw=0
7597 ioff=0
7598 ifilb=0
7604 ! next file name
7608 IF(ios == 0) THEN
7611 ELSE ! failure
7612 iosum=iosum+1
7615 END IF
7616 END IF
7617 ioff=ioff+lfd(i)
7618 END DO
7619
7620 ! open the binary C files
7621
7622 nfilc=-1
7623 ioff=0
7626#ifdef READ_C_FILES
7629 nfilc=0
7630 END IF
7633 ! next file name
7637 IF(ios == 0) THEN
7640 ELSE ! failure
7641 iosum=iosum+1
7644 END IF
7645#else
7646 WRITE(*,*) 'Opening of C-files not supported.'
7647 ! GF add
7648 iosum=iosum+1
7649 ! GF add end
7650#endif
7651 END IF
7652 ioff=ioff+lfd(i)
7653 END DO
7654
7660 END DO
7661
7662 IF(iosum /= 0) THEN
7664 stop 'FILETC: open error '
7665 END IF
7668 stop 'FILETC: no binary files '
7669 END IF
7672 ELSE
7674 END IF
7675101 FORMAT(i3,2x,a)
7676102 FORMAT(a)
7677103 FORMAT(i3,2x,a14,3x,a)
7678 ! CALL mend
7679 RETURN
integer(mpi) function matint(pat, text, npat, ntext)
Approximate string matching.
Definition mptext.f90:290
integer(mpi), dimension(:,:), allocatable kfd
(1,.)= number of records in file, (2,..)= file order
Definition mpmod.f90:285
integer(mpi), dimension(:), allocatable ifd
file: integrated record numbers (=offset)
Definition mpmod.f90:286
integer(mpi), dimension(:), allocatable mfd
file mode: cbinary =1, text =2, fbinary=3
Definition mpmod.f90:282
◆ filetx()
| subroutine filetx |
Interprete text files.
Reset flags and read steering and all other text files. Print some lines from each file.
Store parameter values, constraints and measurements.
Check flags METSOL (method of solution) and MATSTO (matrix storage mode). Set default values for flags, which are undefined.
Parameter values, format:
1 label 2 (initial) parameter value 3 pre-sigma 4 label 5 (initial) parameter value 6 pre-sigma 7 label ... ... (number of words is multiple of 3)
Constraint data, format:
1 0 ! constraint header of four words: 2 right-hand-side ! 0 and -1 ... 3 -1; -2 ! ... indicate (weighting) ... 4 sigma ! ... header 5 label 6 factor 7 label 8 factor 9 ... ... ... (number of words is multiple of 2, at least 6)
Measured data, format:
1 0 ! constraint header of four words: 2 right-hand-side ! 0 and -1 ... 3 -1 ! ... indicate ... 4 sigma ! ... header 5 label 6 factor 7 label 8 factor 9 ... ... ... (number of words is multiple of 2, at least 6)
Definition at line 7733 of file pede.f90.
7735
7736 IMPLICIT NONE
7737 INTEGER(mpi) :: i
7738 INTEGER(mpi) :: ia
7739 INTEGER(mpi) :: ib
7740 INTEGER(mpi) :: ierrf
7741 INTEGER(mpi) :: ioff
7742 INTEGER(mpi) :: ios
7743 INTEGER(mpi) :: iosum
7744 INTEGER(mpi) :: k
7745 INTEGER(mpi) :: mat
7746 INTEGER(mpi) :: nab
7747 INTEGER(mpi) :: nfiln
7748 INTEGER(mpi) :: nline
7749 INTEGER(mpi) :: nlinmx
7750 INTEGER(mpi) :: npat
7751 INTEGER(mpi) :: ntext
7752 INTEGER(mpi) :: matint
7753
7754 ! CALL MSTART('FILETX')
7755
7756 CHARACTER (LEN=1024) :: text
7757 CHARACTER (LEN=1024) :: fname
7758
7759 WRITE(*,*) ' '
7760 WRITE(*,*) 'Processing text files ...'
7761 WRITE(*,*) ' '
7762
7763 iosum=0
7764 ioff=0
7766 IF(i == 0) THEN
7768 nlinmx=100
7769 ELSE
7770 nlinmx=10
7771 ia=ioff
7772 ioff=ioff+lfd(i)
7778 IF(ios /= 0) THEN
7780 iosum=iosum+1
7781 cycle
7782 END IF
7783 END IF
7784
7785 nline=0
7786 nfiln=1
7787 ! read text file
7788 DO
7789 READ(10,102,iostat=ierrf) text
7790 IF (ierrf < 0) THEN
7791 text=' '
7793 WRITE(*,*) ' end-of-file after',nline,' text lines'
7794 EXIT ! eof
7795 ENDIF
7796 nline=nline+1
7799 nab=max(1,nab)
7800 WRITE(*,101) nline,text(1:nab)
7801 IF(nline == nlinmx) WRITE(*,*) ' ...'
7802 END IF
7803
7805 IF(ib == ia+2) THEN
7807 IF(mat == 3) THEN
7808 text=' '
7810 WRITE(*,*) ' end-statement after',nline,' text lines'
7811 EXIT
7812 END IF
7813 END IF
7814
7815 IF(i == 0) THEN ! first text file - exclude lines with file names
7817 nfiln=nfiln+1
7818 text=' '
7819 ! WRITE(*,*) 'line is excluded ',TEXT(1:10)
7820 END IF
7821 END IF
7822 ! WRITE(*,*) TEXT(1:40),' < interprete text'
7824 END DO
7825 WRITE(*,*) ' '
7826 rewind 10
7827 CLOSE(unit=10)
7828 END DO
7829
7830 IF(iosum /= 0) THEN
7832 stop 'FILETX: open error(s) in text files '
7833 END IF
7834
7835 WRITE(*,*) '... end of text file processing.'
7836 WRITE(*,*) ' '
7837
7839 WRITE(*,*) ' '
7841 'or file non-existing,'
7842 WRITE(*,*) ' see above!'
7843 WRITE(*,*) '------------> stop'
7844 WRITE(*,*) ' '
7846 stop
7847 END IF
7848
7849 ! check methods
7850
7853 ! METSOL=1 ! default is matrix inversion
7854 ! MATSTO=1 ! default is symmetric matrix
7859 END IF
7863 matsto=1
7865 ! MATSTO=2 or 1
7867 ! MATSTO=2 or 1
7869 ! MATSTO=2 or 1
7870 ELSE
7871 WRITE(*,*) 'MINRES forced with sparse matrix!'
7872 WRITE(*,*) ' '
7873 WRITE(*,*) 'MINRES forced with sparse matrix!'
7874 WRITE(*,*) ' '
7875 WRITE(*,*) 'MINRES forced with sparse matrix!'
7878 END IF
7880 WRITE(*,*) 'MINRES forced with sparse matrix!'
7881 WRITE(*,*) ' '
7882 WRITE(*,*) 'MINRES forced with sparse matrix!'
7883 WRITE(*,*) ' '
7884 WRITE(*,*) 'MINRES forced with sparse matrix!'
7887 END IF
7888
7889 ! print information about methods and matrix storage modes
7890
7891 WRITE(*,*) ' '
7892 WRITE(*,*) 'Solution method and matrix-storage mode:'
7894 WRITE(*,*) ' METSOL = 1: matrix inversion'
7896 WRITE(*,*) ' METSOL = 2: diagonalization'
7898 WRITE(*,*) ' METSOL = 3: MINRES'
7900 WRITE(*,*) ' METSOL = 4: MINRES-QLP'
7902 WRITE(*,*) ' METSOL = 5: GMRES (-> MINRES)'
7903
7904 END IF
7905
7907
7909 WRITE(*,*) ' MATSTO = 1: symmetric matrix, ', '(n*n+n)/2 elements'
7911 WRITE(*,*) ' MATSTO = 2: sparse matrix'
7912 END IF
7915 END IF
7916
7918 WRITE(*,*) 'Chi square cut equiv 3 st.dev applied ...'
7921 WRITE(*,*) ' (reduced by sqrt in next iterations)'
7922 END IF
7923
7927 END IF
7928
7929 WRITE(*,*) 'Iterations (solutions) with line search:'
7931 WRITE(*,*) ' All'
7933 WRITE(*,*) ' Last'
7935 WRITE(*,*) ' None'
7936 ELSE
7938 WRITE(*,*) ' All with Chi square cut scaling factor <= 1.'
7939 ELSE
7940 WRITE(*,*) ' All'
7941 ENDIF
7942 ENDIF
7943
7945 WRITE(*,*)
7947 ENDIF
7948
7950
7951101 FORMAT(i3,2x,a)
7952102 FORMAT(a)
real(mps) chirem
cut in terms of 3-sigma cut, other iterations, approaching 1.
Definition mpmod.f90:48
integer(mpi) nummeasurements
number of (external) measurements from steering file
Definition mpmod.f90:277
integer(mpi) lsearch
iterations (solutions) with line search: >2: all, =2: all with (next) Chi2 cut scaling factor =1....
Definition mpmod.f90:98
◆ getsum()
| subroutine getsum | ( | real(mpd), intent(out) | asum | ) |
Get accurate sum.
- Parameters
-
[out] asum accurate sum
Definition at line 9066 of file pede.f90.
9068
9069 IMPLICIT NONE
9070 REAL(mpd), INTENT(OUT) ::asum
9072 accuratedsum=0.0_mpd
9073 accuratensum=0
9074 accuratenexp=0
9075 RETURN
◆ ijadd()
| integer(mpl) function ijadd | ( | integer(mpi), intent(in) | itema, |
| integer(mpi), intent(in) | itemb | ||
| ) |
Index for sparse storage.
In case of (compressed) sparse storage calculate index for off-diagonal matrix element.
- Parameters
-
[in] itema row number [in] itemb column number
- Returns
- index (>(<) 0: double(single) precision element, =0: not existing)
Definition at line 4432 of file pede.f90.
4434
4435 IMPLICIT NONE
4436 INTEGER(mpi) :: iencdb
4437 INTEGER(mpi) :: iencdm
4438 INTEGER(mpi) :: isgn
4439 INTEGER(mpi) :: ispc
4440 INTEGER(mpi) :: item2
4441 INTEGER(mpi) :: jtem
4442 INTEGER(mpi) :: jtemc
4443 INTEGER(mpi) :: jtemn
4444
4445 INTEGER(mpi), INTENT(IN) :: itema
4446 INTEGER(mpi), INTENT(IN) :: itemb
4447
4448 INTEGER(mpl) :: ijadd
4449 INTEGER(mpl) :: k
4450 INTEGER(mpl) :: kk
4451 INTEGER(mpl) :: kl
4452 INTEGER(mpl) :: ku
4453 INTEGER(mpl) :: indid
4454 INTEGER(mpl) :: nd
4455 INTEGER(mpl) :: ll
4456 INTEGER(mpl) :: k8
4457 INTEGER(mpl) :: item1
4458 ! ...
4459 ijadd=0
4461 item1=max(itema,itemb) ! larger index
4462 item2=min(itema,itemb) ! smaller index
4463 IF(item2 <= 0.OR.item1 > nd) RETURN
4464 IF(item1 == item2) THEN ! diagonal element
4465 ijadd=item1
4466 RETURN
4467 END IF
4468 ! ! off-diagonal element
4470 iencdm=ishft(1,iencdb)-1
4471 isgn=-1
4475 kl=0
4476 ku=sparsematrixoffsets(1,item1+1)-1-kk
4477 indid=kk
4478 item1=item1+nd+1
4479 isgn=-isgn
4481
4482 ku=((ku+1)*8)/9-1 ! number of regions (-1)
4483 indid=indid+ku/8+1 ! skip group offsets
4484 kl=0
4485 IF(ku < kl) cycle outer ! not found
4486 DO
4487 k=(kl+ku)/2 ! binary search
4489 jtem =ishft(jtemc,-iencdb) ! first column of region
4490 jtemn=jtem+iand(jtemc,iencdm) ! first column after region
4491 IF(item2 >= jtem.AND.item2 < jtemn) EXIT ! found
4492 IF(item2 < jtem) THEN
4493 ku=k-1
4494 ELSE IF(item2 >= jtemn) THEN
4495 kl=k+1
4496 END IF
4497 IF(kl <= ku) cycle
4498 cycle outer ! not found
4499 END DO
4500 k8=k/8 ! region group (-1)
4502 DO kl=k8*8,k-1
4504 END DO
4505 ijadd=ll+item2-jtem
4506
4507 ELSE
4508
4509 IF(ku < kl) cycle outer ! not found
4510 DO
4511 k=(kl+ku)/2 ! binary search
4512 jtem=sparsematrixcolumns(indid+k)
4513 jtemn=jtem
4514 IF(item2 == jtem) EXIT ! found
4515 IF(item2 < jtem) THEN
4516 ku=k-1
4517 ELSE IF(item2 > jtem) THEN
4518 kl=k+1
4519 END IF
4520 IF(kl <= ku) cycle
4521 cycle outer ! not found
4522 END DO
4523 ijadd=ll+k
4524
4525 END IF
4527 RETURN
4528 END DO outer
4529
◆ inone()
| integer(mpi) function inone | ( | integer(mpi), intent(in) | item | ) |
Translate labels to indices (for global parameters).
Functions INONE and subroutine UPONE are used to collect items, i.e. labels, and to order and translate them.
In the first phase items are collected and stored by calling IRES=INONE(ITEM).
At the first entry the two sub-arrays "a" (globalParLabelIndex) and "b" (globalParHashTable) of length 2N are generated with a start length for N=128 entries. In array "a" two words are reserved for each item: (ITEM, count). The function INONE(ITEM) returns the number of the item. At each entry the argument is compared with the already stored items, new items are stored. Search for entries is done using hash-indices, stored in sub-array "b". The initial hash-index is
j = 1 + mod(ITEM, n_prime) + N
where n_prime is the largest prime number less than N. At each entry the count is increased by one. If N items are stored, the size of the sub-arrays is increased by calling CALL UPONE.
- Parameters
-
[in] item label
- Returns
- index
Definition at line 4662 of file pede.f90.
4665
4666 IMPLICIT NONE
4667 INTEGER(mpi), INTENT(IN) :: item
4668 INTEGER(mpi) :: j
4669 INTEGER(mpi) :: k
4670 INTEGER(mpi) :: iprime
4671 INTEGER(mpl) :: length
4672 INTEGER(mpl), PARAMETER :: two = 2
4673
4674 inone=0
4675 IF(item <= 0) RETURN
4677 length=128 ! initial number
4680 globalparhashtable = 0
4688 END IF
4689 outer: DO
4691 inner: DO ! normal case: find item
4692 k=j
4693 j=globalparhashtable(k)
4694 IF(j == 0) EXIT inner ! unused hash code
4696 END DO inner
4697 ! not found
4700 j=0
4701 RETURN
4702 END IF
4705 j=globalparheader(-1)
4710 ! update with larger dimension and redefine index
4715 END DO outer
4716
4720 END IF
4721 inone=j
integer(mpi), dimension(-7:0) globalparheader
global parameters (mapping) header
Definition mpmod.f90:217
integer(mpi), dimension(:), allocatable globalparhashtable
global parameters hash table
Definition mpmod.f90:215
◆ intext()
| subroutine intext | ( | character (len=*), intent(in) | text, |
| integer(mpi), intent(in) | nline | ||
| ) |
Interprete text.
Look for keywords and argument in text line.
- Parameters
-
[in] text text [in] nline line number
Definition at line 8019 of file pede.f90.
8022
8023 IMPLICIT NONE
8024 INTEGER(mpi) :: i
8025 INTEGER(mpi) :: ia
8026 INTEGER(mpi) :: ib
8027 INTEGER(mpi) :: ier
8028 INTEGER(mpi) :: iomp
8029 INTEGER(mpi) :: k
8030 INTEGER(mpi) :: kkey
8031 INTEGER(mpi) :: label
8032 INTEGER(mpi) :: lkey
8033 INTEGER(mpi) :: mat
8034 INTEGER(mpi) :: miter
8035 INTEGER(mpi) :: nab
8036 INTEGER(mpi) :: nkey
8037 INTEGER(mpi) :: nkeys
8038 INTEGER(mpi) :: nl
8039 INTEGER(mpi) :: nmeth
8040 INTEGER(mpi) :: npat
8041 INTEGER(mpi) :: ntext
8042 INTEGER(mpi) :: nums
8043 INTEGER(mpi) :: matint
8044
8045 CHARACTER (LEN=*), INTENT(IN) :: text
8046 INTEGER(mpi), INTENT(IN) :: nline
8047
8048 parameter(nkeys=5,nmeth=6)
8049 CHARACTER (LEN=16) :: methxt(nmeth)
8050 CHARACTER (LEN=16) :: keylst(nkeys)
8051 CHARACTER (LEN=32) :: keywrd
8052 CHARACTER (LEN=32) :: keystx
8053 REAL(mpd) :: dnum(100)
8054 INTEGER(mpi) :: lpvs ! ... integer
8055 REAL(mpd) :: plvs ! ... float
8056
8057 INTERFACE
8060 INTEGER(mpi), INTENT(IN OUT) :: length
8061 TYPE(listItem), DIMENSION(:), INTENT(IN OUT), ALLOCATABLE :: list
8062 INTEGER(mpi), INTENT(IN) :: label
8063 REAL(mpd), INTENT(IN) :: value
8065 END INTERFACE
8066
8067 DATA keylst/'unknown','parameter','constraint','measurement','method'/
8068
8069 SAVE
8070 DATA methxt/'diagonalization','inversion','fullMINRES', 'sparseMINRES', &
8071 'fullMINRES-QLP', 'sparseMINRES-QLP'/
8072 DATA lkey/-1/ ! last keyword
8073
8074 ! ...
8075 nkey=-1 ! new keyword
8077 IF(nab == 0) GOTO 10
8079
8080 IF(nums /= 0) nkey=0
8083 ! WRITE(*,*) 'Keyword is ',KEYWRD
8084
8085 ! compare keywords
8086
8087 DO nkey=2,nkeys ! loop over all pede keywords
8088 keystx=keylst(nkey) ! copy NKEY.th pede keyword
8090 IF(mat >= ntext-ntext/5) GO TO 10
8091 END DO
8092
8093 ! more comparisons
8094
8095 keystx='print'
8097 ! WRITE(*,*) KEYSTX,MAT,NTEXT
8098 ! IF(MAT.GE.NTEXT) THEN
8099 IF(mat >= (npat-npat/5)) THEN
8100 ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN
8101 mprint=1
8102 IF(nums > 0) mprint=nint(dnum(1),mpi)
8103 RETURN
8104 END IF
8105
8106 keystx='debug'
8108 IF(mat >= (npat-npat/5)) THEN
8109 ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN
8110 mdebug=3
8111 ! GF IF(NUMS.GT.0) MPRINT=DNUM(1)
8112 IF(nums > 0) mdebug=nint(dnum(1),mpi)
8113 IF(nums > 1) mdebg2=nint(dnum(2),mpi)
8114 RETURN
8115 END IF
8116
8117 keystx='entries'
8119 IF(mat >= (npat-npat/5)) THEN
8120 ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN
8121 IF(nums > 0 .AND. dnum(1) > 0.5) mreqenf=nint(dnum(1),mpi)
8122 IF(nums > 1 .AND. dnum(2) > 0.5) mreqena=nint(dnum(2),mpi)
8123 IF(nums > 2 .AND. dnum(3) > 0.5) iteren=nint(dnum(1)*dnum(3),mpi)
8124 RETURN
8125 END IF
8126
8127 keystx='printrecord'
8129 IF(mat >= (npat-npat/5)) THEN
8130 IF(nums > 0) nrecpr=nint(dnum(1),mpi)
8131 IF(nums > 1) nrecp2=nint(dnum(2),mpi)
8132 RETURN
8133 END IF
8134
8135 keystx='maxrecord'
8137 IF(mat >= (npat-npat/5)) THEN
8138 IF (nums > 0.AND.dnum(1) > 0.) mxrec=nint(dnum(1),mpi)
8139 RETURN
8140 END IF
8141
8142 keystx='cache'
8144 IF(mat >= (npat-npat/5)) THEN
8145 IF (nums > 0.AND.dnum(1) >= 0.) ncache=nint(dnum(1),mpi) ! cache size, <0 keeps default
8146 IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) & ! read cache fill level
8147 fcache(1)=real(dnum(2),mps)
8148 IF (nums >= 4) THEN ! explicit cache splitting
8149 DO k=1,3
8150 fcache(k)=real(dnum(k+1),mps)
8151 END DO
8152 END IF
8153 RETURN
8154 END IF
8155
8156 keystx='chisqcut'
8158 IF(mat >= (npat-npat/5)) THEN
8159 IF(nums == 0) THEN ! always 3-sigma cut
8160 chicut=1.0
8161 chirem=1.0
8162 ELSE
8163 chicut=real(dnum(1),mps)
8164 IF(chicut < 1.0) chicut=-1.0
8165 IF(nums == 1) THEN
8166 chirem=1.0 ! 3-sigma cut, if not specified
8167 ELSE
8168 chirem=real(dnum(2),mps)
8169 IF(chirem < 1.0) chirem=1.0
8170 IF(chicut >= 1.0) chirem=min(chirem,chicut)
8171 END IF
8172 END IF
8173 RETURN
8174 END IF
8175
8176 ! GF added:
8177 keystx='hugecut'
8179 IF(mat >= (npat-npat/5)) THEN
8180 IF(nums > 0) chhuge=real(dnum(1),mps)
8181 IF(chhuge < 1.0) chhuge=1.0 ! at least (!!) 3-sigma
8182 RETURN
8183 END IF
8184 ! GF added end
8185
8186 keystx='linesearch'
8188 IF(mat >= (npat-npat/5)) THEN
8189 IF(nums > 0) lsearch=nint(dnum(1),mpi)
8190 RETURN
8191 END IF
8192
8193 keystx='localfit'
8195 IF(mat >= (npat-npat/5)) THEN
8196 IF(nums > 0) lfitnp=nint(dnum(1),mpi)
8197 IF(nums > 1) lfitbb=nint(dnum(2),mpi)
8198 RETURN
8199 END IF
8200
8201 keystx='regularization'
8203 IF(mat >= (npat-npat/5)) THEN
8204 nregul=1
8205 regula=real(dnum(1),mps)
8206 IF(nums >= 2) regpre=real(dnum(2),mps)
8207 RETURN
8208 END IF
8209
8210 keystx='regularisation'
8212 IF(mat >= (npat-npat/5)) THEN
8213 nregul=1
8214 regula=real(dnum(1),mps)
8215 IF(nums >= 2) regpre=real(dnum(2),mps)
8216 RETURN
8217 END IF
8218
8219 keystx='presigma'
8221 IF(mat >= (npat-npat/5)) THEN
8222 regpre=real(dnum(1),mps)
8223 RETURN
8224 END IF
8225
8226 keystx='matiter'
8228 IF(mat >= (npat-npat/5)) THEN
8229 matrit=nint(dnum(1),mpi)
8230 RETURN
8231 END IF
8232
8233 keystx='matmoni'
8235 IF(mat >= (npat-npat/5)) THEN
8236 matmon=-1
8237 IF (nums > 0.AND.dnum(1) > 0.) matmon=nint(dnum(1),mpi)
8238 RETURN
8239 END IF
8240
8241 keystx='bandwidth'
8243 ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN
8244 IF(mat >= (npat-npat/5)) THEN
8245 IF(nums > 0) mbandw=nint(dnum(1),mpi)
8246 IF(mbandw < 0) mbandw=-1
8247 IF(nums > 1) lprecm=nint(dnum(2),mpi)
8248 RETURN
8249 END IF
8250
8251 ! KEYSTX='outlierrejection'
8252 ! MAT=MATINT(TEXT(KEYA:KEYB),KEYSTX,NPAT,NTEXT) ! comparison
8253 ! WRITE(*,*) KEYSTX,MAT,(NTEXT+NTEXT)/3
8254 ! IF(MAT.GE.(NTEXT+NTEXT+NTEXT-2)/3) THEN
8255 ! IF(MAT.GE.(NPAT-NPAT/5)) THEN
8256 ! CHDFRJ=DNUM(1)
8257 ! IF(CHDFRJ.LT.3.0) CHDFRJ=100.0
8258 ! RETURN
8259 ! END IF
8260
8261 ! KEYSTX='outliersuppression'
8262 ! MAT=MATINT(TEXT(KEYA:KEYB),KEYSTX,NPAT,NTEXT) ! comparison
8263 ! WRITE(*,*) KEYSTX,MAT,(NTEXT+NTEXT)/3
8264 ! IF(MAT.GE.(NTEXT+NTEXT+NTEXT-2)/3) THEN
8265 ! IF(MAT.GE.(NPAT-NPAT/5)) THEN
8266 ! LHUBER=DNUM(1)
8267 ! IF(LHUBER.LE.2) LHUBER=2 ! at least 2 Huber iterations
8268 ! RETURN
8269 ! END IF
8270
8271 keystx='outlierdownweighting'
8273 ! WRITE(*,*) KEYSTX,MAT,(NTEXT+NTEXT)/3
8274 ! IF(MAT.GE.(NTEXT+NTEXT+NTEXT-2)/3) THEN
8275 IF(mat >= (npat-npat/5)) THEN
8276 lhuber=nint(dnum(1),mpi)
8277 IF(lhuber > 0.AND.lhuber <= 2) lhuber=2 ! at least 2 Huber iterations (if any)
8278 RETURN
8279 END IF
8280
8281 keystx='dwfractioncut'
8283 IF(mat >= (npat-npat/5)) THEN
8284 dwcut=real(dnum(1),mps)
8285 IF(dwcut > 0.5) dwcut=0.5
8286 RETURN
8287 END IF
8288
8289 keystx='pullrange'
8291 IF(mat >= (npat-npat/5)) THEN
8292 prange=abs(real(dnum(1),mps))
8293 RETURN
8294 END IF
8295
8296 keystx='subito'
8298 IF(mat >= (npat-npat/5)) THEN
8299 isubit=1
8300 RETURN
8301 END IF
8302
8303 keystx='force'
8305 IF(mat >= (npat-npat/5)) THEN
8306 iforce=1
8307 RETURN
8308 END IF
8309
8310 keystx='memorydebug'
8312 IF(mat >= (npat-npat/5)) THEN
8313 memdbg=1
8314 IF (nums > 0.AND.dnum(1) > 0.0) memdbg=nint(dnum(1),mpi)
8315 RETURN
8316 END IF
8317
8318 keystx='globalcorr'
8320 IF(mat >= (npat-npat/5)) THEN
8321 igcorr=1
8322 RETURN
8323 END IF
8324
8325 keystx='printcounts'
8327 IF(mat >= (npat-npat/5)) THEN
8328 ipcntr=1
8329 IF (nums > 0.AND.dnum(1) > 0.0) ipcntr=nint(dnum(1),mpi)
8330 RETURN
8331 END IF
8332
8333 keystx='weightedcons'
8335 IF(mat >= (npat-npat/5)) THEN
8336 iwcons=1
8337 IF (nums > 0) iwcons=nint(dnum(1),mpi)
8338 RETURN
8339 END IF
8340
8341 keystx='skipemptycons'
8343 IF(mat >= (npat-npat/5)) THEN
8344 iskpec=1
8345 RETURN
8346 END IF
8347
8348 keystx='withelimination'
8350 IF(mat >= (npat-npat/5)) THEN
8351 icelim=1
8352 RETURN
8353 END IF
8354
8355 keystx='withmultipliers'
8357 IF(mat >= (npat-npat/5)) THEN
8358 icelim=0
8359 RETURN
8360 END IF
8361
8362 keystx='checkinput'
8364 IF(mat >= (npat-npat/5)) THEN
8365 icheck=1
8366 IF (nums > 0) icheck=nint(dnum(1),mpi)
8367 RETURN
8368 END IF
8369
8370 keystx='monitorresiduals'
8372 IF(mat >= (npat-npat/5)) THEN
8373 imonit=3
8374 IF (nums > 0) imonit=nint(dnum(1),mpi)
8375 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
8376 RETURN
8377 END IF
8378
8379 keystx='monitorpulls'
8381 IF(mat >= (npat-npat/5)) THEN
8382 imonit=3
8383 imonmd=1
8384 IF (nums > 0) imonit=nint(dnum(1),mpi)
8385 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
8386 RETURN
8387 END IF
8388
8389 keystx='scaleerrors'
8391 IF(mat >= (npat-npat/5)) THEN
8392 iscerr=1
8393 IF (nums > 0) dscerr(1:2)=dnum(1)
8394 IF (nums > 1) dscerr(2)=dnum(2)
8395 RETURN
8396 END IF
8397
8398 keystx='iterateentries'
8400 IF(mat >= (npat-npat/5)) THEN
8401 iteren=huge(iteren)
8402 IF (nums > 0) iteren=nint(dnum(1),mpi)
8403 RETURN
8404 END IF
8405
8406 keystx='threads'
8408 IF(mat >= (npat-npat/5)) THEN
8409 iomp=0
8410 !$ IOMP=1
8411 !$ IF (IOMP.GT.0) THEN
8412 !$ IF (NUMS.GE.1.AND.DNUM(1).GT.0.) MTHRD =NINT(dnum(1),mpi)
8413 !$ IF (NUMS.GE.2) THEN
8414 !$ MTHRDR=MTHRD
8415 !$ IF (DNUM(2).GT.0.) MTHRDR=NINT(dnum(2),mpi)
8416 !$ ENDIF
8417 !$ ELSE
8418 WRITE(*,*) 'WARNING: multithreading not available'
8419 !$ ENDIF
8420 RETURN
8421 END IF
8422
8423 keystx='compress'
8425 IF(mat >= (npat-npat/5)) THEN
8426 mcmprs=1
8427 RETURN
8428 END IF
8429
8430 ! still experimental
8431 !keystx='extendedStorage'
8432 !mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
8433 !IF(mat >= (npat-npat/5)) THEN
8434 ! mextnd=1
8435 ! ! compression enforced for extended storage (in mpbits)
8436 ! RETURN
8437 !END IF
8438
8439 keystx='errlabels'
8441 IF(mat >= (npat-npat/5).AND.mnrsel < 100) THEN
8442 nl=min(nums,100-mnrsel)
8443 DO k=1,nl
8444 lbmnrs(mnrsel+k)=nint(dnum(k),mpi)
8445 END DO
8446 mnrsel=mnrsel+nl
8447 RETURN
8448 END IF
8449
8450 keystx='pairentries'
8452 IF(mat >= (npat-npat/5)) THEN
8453 ! This option could be implemented to get rid of parameter pairs
8454 ! that have very few entries - to save matrix memory size.
8455 IF (nums > 0.AND.dnum(1) > 0.0) THEN
8456 mreqpe=nint(dnum(1),mpi)
8457 IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=nint(dnum(2),mpi)
8458 IF (nums >= 3.AND.dnum(3) >= dnum(1)) msngpe=nint(dnum(3),mpi)
8459 END IF
8460 RETURN
8461 END IF
8462
8463 keystx='wolfe'
8465 IF(mat >= (npat-npat/5)) THEN
8466 wolfc1=real(dnum(1),mps)
8467 wolfc2=real(dnum(2),mps)
8468 RETURN
8469 END IF
8470
8471 ! GF added:
8472 ! convergence tolerance for minres:
8473 keystx='mrestol'
8475 IF(mat >= (npat-npat/5)) THEN
8476 IF(nums > 0) THEN
8477 IF (dnum(1) < 1.0e-10_mpd.OR.dnum(1) > 1.0e-04_mpd) THEN
8478 WRITE(*,*) 'ERROR: need 1.0D-10 <= MRESTL ', &
8479 '<= 1.0D-04, but get ', dnum(1)
8480 ELSE
8481 mrestl=dnum(1)
8482 END IF
8483 END IF
8484 RETURN
8485 END IF
8486 ! GF added end
8487
8488 keystx='mrestranscond'
8490 IF(mat >= (npat-npat/5)) THEN
8491 IF(nums > 0) THEN
8492 mrtcnd = dnum(1)
8493 END IF
8494 RETURN
8495 END IF
8496
8497 keystx='mresmode'
8499 IF(mat >= (npat-npat/5)) THEN
8500 IF(nums > 0) THEN
8501 mrmode = int(dnum(1),mpi)
8502 END IF
8503 RETURN
8504 END IF
8505
8506 keystx='nofeasiblestart'
8508 IF(mat >= (npat-npat/5)) THEN
8509 nofeas=1 ! do not make parameters feasible at start
8510 RETURN
8511 END IF
8512
8513 keystx='histprint'
8515 IF(mat >= ntext-ntext/10) THEN
8516 ! IF(MAT.GE.(NPAT-NPAT/5)) THEN
8517 nhistp=1 ! print histograms
8518 RETURN
8519 END IF
8520
8521 keystx='fortranfiles'
8523 IF(mat >= ntext-ntext/10) RETURN
8524
8525 keystx='Cfiles'
8527 IF(mat >= ntext-ntext/10) RETURN
8528
8529 keystx='closeandreopen'
8531 IF(mat >= ntext-ntext/10) RETURN
8532
8533 keystx=keylst(1)
8534 nkey=1 ! unknown keyword
8535 IF(nums /= 0) nkey=0
8536
8537 WRITE(*,*) ' '
8538 WRITE(*,*) '**************************************************'
8539 WRITE(*,*) ' '
8540 WRITE(*,*) 'Unknown keyword(s): ',text(1:min(nab,50))
8541 WRITE(*,*) ' '
8542 WRITE(*,*) '**************************************************'
8543 WRITE(*,*) ' '
8544 lunkno=lunkno+1
8545
8546 END IF
8547 ! result: NKEY = -1 blank
8548 ! NKEY = 0 numerical data, no text keyword or unknown
8549 ! NKEY > 0 keyword NKEY from list, keyword = KEYSTX
8550
8551
8552 ! content/lastcontent
8553 ! -------------------
8554 ! blank -1
8555 ! data 0
8556 ! keyword
8557 ! unknown 1
8558 ! parameter 2
8559 ! constraint 3
8560 ! measurement 4
8561 ! method 5
8562
8563
856410 IF(nkey > 0) THEN ! new keyword
8565 lkey=nkey
8566 IF(lkey == 2) THEN ! parameter
8567 IF(nums == 3) THEN
8568 lpvs=nint(dnum(1),mpi) ! label
8569 IF(lpvs /= 0) THEN
8572 ELSE
8573 WRITE(*,*) 'Line',nline,' error, label=',lpvs
8574 END IF
8575 ELSE IF(nums /= 0) THEN
8576 kkey=1 ! switch to "unknown" ?
8577 WRITE(*,*) 'Wrong text in line',nline
8578 WRITE(*,*) 'Status: new parameter'
8579 WRITE(*,*) '> ',text(1:nab)
8580 END IF
8581 ELSE IF(lkey == 3) THEN ! constraint
8582 ! WRITE(*,*) 'Keyword is constraint!',NUMS,' numerical data'
8583 IF(nums >= 1.AND.nums <= 2) THEN ! start constraint
8584 lpvs=0 ! r = r.h.s. value
8586 lpvs=-1 ! constraint
8587 IF(iwcons > 0) lpvs=-2 ! weighted constraint
8588 plvs=0.0
8589 IF(nums == 2) plvs=dnum(2) ! sigma
8591 ELSE
8592 kkey=1 ! switch to "unknown"
8593 WRITE(*,*) 'Wrong text in line',nline
8594 WRITE(*,*) 'Status: new keyword constraint'
8595 WRITE(*,*) '> ',text(1:nab)
8596 END IF
8597 ELSE IF(lkey == 4) THEN ! measurement
8598 IF(nums == 2) THEN ! start measurement
8599 nummeasurements=nummeasurements+1
8600 lpvs=0 ! r = r.h.s. value
8602 lpvs=-1 ! sigma
8604 ELSE
8605 kkey=1 ! switch to "unknown"
8606 WRITE(*,*) 'Wrong text in line',nline
8607 WRITE(*,*) 'Status: new keyword measurement'
8608 WRITE(*,*) '> ',text(1:nab)
8609 END IF
8610
8611 ELSE IF(lkey == 5) THEN ! method
8612 miter=mitera
8613 IF(nums >= 1) miter=nint(dnum(1),mpi)
8614 IF(miter >= 1) mitera=miter
8615 dflim=real(dnum(2),mps)
8616 lkey=0
8617 DO i=1,nmeth
8618 keystx=methxt(i)
8620 IF(mat >= ntext-ntext/5) THEN
8621 IF(i == 1) THEN ! diagonalization
8622 metsol=2
8623 matsto=1
8624 ELSE IF(i == 2) THEN ! inversion
8625 metsol=1
8626 matsto=1
8627 ELSE IF(i == 3) THEN ! fullMINRES
8628 metsol=3
8629 matsto=1
8630 ELSE IF(i == 4) THEN ! sparseMINRES
8631 metsol=3
8632 matsto=2
8633 ELSE IF(i == 5) THEN ! fullMINRES-QLP
8634 metsol=4
8635 matsto=1
8636 ELSE IF(i == 6) THEN ! sparseMINRES-QLP
8637 metsol=4
8638 matsto=2
8639 END IF
8640 END IF
8641 END DO
8642 END IF
8643 ELSE IF(nkey == 0) THEN ! data for continuation
8644 IF(lkey == 2) THEN ! parameter
8645 IF(nums >= 3) THEN ! store data from this line
8646 lpvs=nint(dnum(1),mpi) ! label
8647 IF(lpvs /= 0) THEN
8650 ELSE
8651 WRITE(*,*) 'Line',nline,' error, label=',lpvs
8652 END IF
8653 ELSE IF(nums > 1.AND.nums < 3) THEN
8654 kkey=1 ! switch to "unknown" ?
8655 WRITE(*,*) 'Wrong text in line',nline
8656 WRITE(*,*) 'Status continuation parameter'
8657 WRITE(*,*) '> ',text(1:nab)
8658 END IF
8659
8660 ELSE IF(lkey == 3) THEN ! constraint
8661 ier=0
8662 DO i=1,nums,2
8663 label=nint(dnum(i),mpi)
8664 IF(label <= 0) ier=1
8665 END DO
8666 IF(mod(nums,2) /= 0) ier=1 ! reject odd number
8667 IF(ier == 0) THEN
8668 DO i=1,nums,2
8669 lpvs=nint(dnum(i),mpi) ! label
8670 plvs=dnum(i+1) ! factor
8672 END DO
8673 ELSE
8674 kkey=0
8675 WRITE(*,*) 'Wrong text in line',nline
8676 WRITE(*,*) 'Status continuation constraint'
8677 WRITE(*,*) '> ',text(1:nab)
8678 END IF
8679
8680 ELSE IF(lkey == 4) THEN ! measurement
8681 ! WRITE(*,*) 'continuation < ',NUMS
8682 ier=0
8683 DO i=1,nums,2
8684 label=nint(dnum(i),mpi)
8685 IF(label <= 0) ier=1
8686 END DO
8687 IF(mod(nums,2) /= 0) ier=1 ! reject odd number
8688 ! WRITE(*,*) 'IER NUMS ',IER,NUMS
8689 IF(ier == 0) THEN
8690 DO i=1,nums,2
8691 lpvs=nint(dnum(i),mpi) ! label
8692 plvs=dnum(i+1) ! factor
8694 END DO
8695 ELSE
8696 kkey=0
8697 WRITE(*,*) 'Wrong text in line',nline
8698 WRITE(*,*) 'Status continuation measurement'
8699 WRITE(*,*) '> ',text(1:nab)
8700 END IF
8701
8702 END IF
8703 END IF
◆ iprime()
| integer(mpi) function iprime | ( | integer(mpi), intent(in) | n | ) |
largest prime number < N.
- Parameters
-
[in] n N
- Returns
- largest prime number < N
Definition at line 4789 of file pede.f90.
4791
4792 IMPLICIT NONE
4793 INTEGER(mpi), INTENT(IN) :: n
4794 INTEGER(mpi) :: nprime
4795 INTEGER(mpi) :: nsqrt
4796 INTEGER(mpi) :: i
4797 ! ...
4798 SAVE
4799 nprime=n ! max number
4800 IF(mod(nprime,2) == 0) nprime=nprime+1 ! ... odd number
4801 outer: DO
4802 nprime=nprime-2 ! next lower odd number
4804 DO i=3,nsqrt,2 !
4805 IF(i*(nprime/i) == nprime) cycle outer ! test prime number
4806 END DO
4807 EXIT outer ! found
4808 END DO outer
4809 iprime=nprime
◆ isjajb()
| subroutine isjajb | ( | integer(mpi), intent(in) | nst, |
| integer(mpi), intent(inout) | is, | ||
| integer(mpi), intent(out) | ja, | ||
| integer(mpi), intent(out) | jb, | ||
| integer(mpi), intent(out) | jsp | ||
| ) |
Decode Millepede record.
Get indices JA, JB, IS for next measurement within record:
- Measurement is:
GLDER(JA) - Local derivatives are:
(INDER(JA+J),GLDER(JA+J),J=1,JB-JA-1)i.e. JB-JA-1 derivatives - Standard deviation is:
GLDER(JB) - Global derivatives are:
(INDER(JB+J),GLDER(JB+J),J=1,IS-JB)i.e. IST-JB derivatives
End_of_data is indicated by returned values JA=0 and JB=0 Special data are ignored. At end_of_data the info to the special data is returned: IS = pointer to special data; number of words is NSP=-GLDER(IS).
- Parameters
-
[in] nst index of last word of record [in,out] is index of last global derivative (index of first word of record at the first call) [out] ja index of measured value (=0 at end), = pointer to local derivatives [out] jb index of standard deviation (=0 at end), = pointer to global derivatives [out] jsp index to special data
Definition at line 2259 of file pede.f90.
2261
2262 IMPLICIT NONE
2263 REAL(mpr8) :: glder
2264 INTEGER(mpi) :: i
2265 INTEGER(mpi) :: inder
2266
2267 INTEGER(mpi), INTENT(IN) :: nst
2268 INTEGER(mpi), INTENT(IN OUT) :: is
2269 INTEGER(mpi), INTENT(OUT) :: ja
2270 INTEGER(mpi), INTENT(OUT) :: jb
2271 INTEGER(mpi), INTENT(OUT) :: jsp
2272 SAVE
2273 ! ...
2274 inder(i)=readbufferdatai(i)
2275 glder(i)=readbufferdatad(i)
2276
2277 jsp=0
2278 DO
2279 ja=0
2280 jb=0
2281 IF(is >= nst) RETURN
2282 DO
2283 is=is+1
2284 IF(inder(is) == 0) EXIT
2285 END DO
2286 ja=is
2287 DO
2288 is=is+1
2289 IF(inder(is) == 0) EXIT
2290 END DO
2291 jb=is
2292 IF(ja+1 == jb.AND.glder(jb) < 0.0_mpr8) THEN
2293 ! special data
2294 jsp=jb ! pointer to special data
2295 is=is+nint(-glder(jb),mpi) ! skip NSP words
2296 cycle
2297 END IF
2298 DO WHILE(inder(is+1) /= 0.AND.is < nst)
2299 is=is+1
2300 END DO
2301 EXIT
2302 END DO
2303
◆ loop1()
| subroutine loop1 |
First data loop (get global labels).
Read all data files and add all labels to global labels table, add labels from parameters, constraints and measurements (from text files).
Define variable and fixed global parameters (depending on entries and pre-sigma).
Iterate if records had been skipped due to too small read buffer size.
Definition at line 4821 of file pede.f90.
4824
4825 IMPLICIT NONE
4826 INTEGER(mpi) :: i
4827 INTEGER(mpi) :: idum
4828 INTEGER(mpi) :: in
4829 INTEGER(mpi) :: indab
4830 INTEGER(mpi) :: itgbi
4831 INTEGER(mpi) :: itgbl
4832 INTEGER(mpi) :: ivgbi
4833 INTEGER(mpi) :: j
4834 INTEGER(mpi) :: mqi
4835 INTEGER(mpi) :: nc31
4836 INTEGER(mpi) :: nr
4837 INTEGER(mpi) :: nwrd
4838 INTEGER(mpi) :: inone
4839 REAL(mpd) :: param
4840 REAL(mpd) :: presg
4841 REAL(mpd) :: prewt
4842
4843 INTEGER(mpl) :: length
4844 SAVE
4845 ! ...
4849 ! add labels from parameter, constraints, measurements -------------
4852 END DO
4855 END DO
4858 END DO
4861 END DO
4862
4865 END IF
4867
4868 DO
4871 END DO
4872
4873 ! read all data files and add all labels to global labels table ----
4874
4876 WRITE(*,*) 'Read all binary data files:'
4877 END IF
4878 CALL hmpldf(1,'Number of words/record in binary file')
4880 ! define read buffer
4882 nwrd=nc31+1
4883 length=nwrd*mthrdr
4885 nwrd=nc31*10+2+ndimbuf
4886 length=nwrd*mthrdr
4889 ! to read (old) float binary files
4892
4893 DO
4896 IF(nr <= 0) EXIT ! end of data?
4897 END DO
4898 ! release read buffer
4904 EXIT
4905 ELSE
4907 END IF
4908 END DO
4909
4911 CALL hmprnt(1)
4912 CALL hmprnt(8)
4913 END IF
4914 CALL hmpwrt(1)
4915 CALL hmpwrt(8)
4919 stop 'LOOP1: no labels/parameters defined'
4920 END IF
4923 ' is total number NTGB of labels/parameters'
4924 ! histogram number of entries per label ----------------------------
4925 CALL hmpldf(2,'Number of entries per label')
4928 END DO
4930 CALL hmpwrt(2) ! write to his file
4931
4932 ! three subarrays for all global parameters ------------------------
4933 length=ntgb
4935 globalparameter=0.0_mpd
4937 globalparpresigma=0.
4939 globalparstart=0.
4943 globalparcons=0
4944
4946 param=listparameters(i)%value
4948 IF(in /= 0) THEN
4949 globalparameter(in)=param
4950 globalparstart(in)=param
4951 ENDIF
4952 END DO
4953
4954 npresg=0
4956 presg=listpresigmas(i)%value
4958 IF(in /= 0) THEN
4961 END IF
4962 END DO
4966
4967 ! determine flag variable (active) or fixed (inactive) -------------
4968
4969 indab=0
4976 ELSE
4977 indab=indab+1
4979 END IF
4980 END DO
4985
4986 ! translation table of length NVGB of total global indices ---------
4987 length=nvgb
4989 indab=0
4992 indab=indab+1
4993 globalparvartototal(indab)=i
4994 END IF
4995 END DO
4996
4997 ! regularization ---------------------------------------------------
5000 ' (if no individual pre-sigma defined)'
5002
5004 WRITE(*,*) 'No regularization will be done'
5005 ELSE
5007 END IF
5008112 FORMAT(a,e9.2,a)
5011 stop '... no variable global parameters'
5012 ENDIF
5013
5017 prewt=0.0 ! pre-weight
5018 IF(presg > 0.0) THEN
5019 prewt=1.0/presg**2 ! 1/presigma^2
5022 END IF
5024 END DO
5025
5026 ! WRITE(*,*) 'GlPa_index GlPa_label array1 array6'
5028 itgbl=globalparlabelindex(1,i)
5029 ivgbi=globalparlabelindex(2,i)
5030 IF(ivgbi > 0) THEN
5031 ! WRITE(*,111) I,ITGBL,QM(IND1+I),QM(IND6+IVGBI)
5032 ELSE
5033 ! WRITE(*,111) I,ITGBL,QM(IND1+I)
5034 END IF
5035 END DO
5036 ! 111 FORMAT(I5,I10,F10.5,E12.4)
5039
5040 ! print overview over important numbers ----------------------------
5041
5044 WRITE(*,*) ' '
5058 WRITE(*,*) ' '
5059 WRITE(*,*) 'Global parameter labels:'
5060 mqi=ntgb
5061 IF(mqi <= 100) THEN
5063 ELSE
5065 WRITE(*,*) ' ...'
5066 mqi=((mqi-20)/20)*20+1
5068 END IF
5069 END IF
5070 WRITE(*,*) ' '
5071 WRITE(*,*) ' '
5072 END IF
5073 WRITE(8,*) ' '
5080
5084
5085101 FORMAT(1x,a8,' =',i10,' = ',a)
type(listitem), dimension(:), allocatable listparameters
list of parameters from steering file
Definition mpmod.f90:272
integer(mpi) iteren
entries cut is iterated for parameters with less entries (if > mreqenf)
Definition mpmod.f90:104
type(listitem), dimension(:), allocatable listmeasurements
list of (external) measurements from steering file
Definition mpmod.f90:279
integer(mpi) ndimbuf
default read buffer size (I/F words, half record length)
Definition mpmod.f90:304
real(mps) regula
regularization parameter, add regula * norm(global par.) to objective function
Definition mpmod.f90:70
real(mpd), dimension(:), allocatable globalparpreweight
weight from pre-sigma
Definition mpmod.f90:179
real(mpd), dimension(:), allocatable globalparpresigma
pre-sigma for global parameters
Definition mpmod.f90:178
integer(mpi) mreqenf
required number of entries (for variable global parameter from binary Files)
Definition mpmod.f90:39
integer(mpi) lenmeasurements
length of list of (external) measurements from steering file
Definition mpmod.f90:278
integer(mpi), dimension(:), allocatable globalparcounts
global parameters counts (from binary files)
Definition mpmod.f90:180
real(mpd), dimension(:), allocatable globalparstart
start value for global parameters
Definition mpmod.f90:177
integer(mpi) mreqena
required number of entries (for variable global parameter from Accepted local fits)
Definition mpmod.f90:40
integer(mpi), dimension(:), allocatable readbufferpointer
pointer to used buffers
Definition mpmod.f90:234
type(listitem), dimension(:), allocatable listpresigmas
list of pre-sgmas from steering file
Definition mpmod.f90:274
◆ loop1i()
| subroutine loop1i |
Iteration of first data loop.
Read all data files again skipping measurements with any parameter below the entries cut to update the number of entries.
Redefine variable and fixed global parameters (depending on updated entries).
Definition at line 5095 of file pede.f90.
5098
5099 IMPLICIT NONE
5100 INTEGER(mpi) :: i
5101 INTEGER(mpi) :: ibuf
5102 INTEGER(mpi) :: ij
5103 INTEGER(mpi) :: indab
5104 INTEGER(mpi) :: inder
5105 INTEGER(mpi) :: isfrst
5106 INTEGER(mpi) :: islast
5107 INTEGER(mpi) :: ist
5108 INTEGER(mpi) :: j
5109 INTEGER(mpi) :: ja
5110 INTEGER(mpi) :: jb
5111 INTEGER(mpi) :: jsp
5112 INTEGER(mpi) :: nc31
5113 INTEGER(mpi) :: nr
5114 INTEGER(mpi) :: nlow
5115 INTEGER(mpi) :: nst
5116 INTEGER(mpi) :: nwrd
5117 REAL(mpr8) :: glder
5118
5119 INTEGER(mpl) :: length
5120 INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: newCounter
5121 SAVE
5122
5123 isfrst(ibuf)=readbufferpointer(ibuf)+1
5125 inder(i)=readbufferdatai(i)
5126 glder(i)=readbufferdatad(i)
5127 ! ...
5130 WRITE(*,*) ' '
5131 WRITE(*,*) 'LOOP1: iterating'
5132
5133 length=ntgb
5135 newcounter=0
5136
5137 ! define read buffer
5139 nwrd=nc31+1
5140 length=nwrd*mthrdr
5142 nwrd=nc31*10+2+ndimbuf
5143 length=nwrd*mthrdr
5146 ! to read (old) float binary files
5149
5150 DO
5154 ist=isfrst(ibuf)
5155 nst=islast(ibuf)
5156 nwrd=nst-ist+1
5157 DO ! loop over measurements
5159 IF(ja == 0.AND.jb == 0) EXIT
5160 IF(ja /= 0) THEN
5161 nlow=0
5162 DO j=1,ist-jb
5163 ij=inder(jb+j) ! index of global parameter
5165 IF(ij == -2) nlow=nlow+1 ! fixed by entries cut
5166 END DO
5167 IF(nlow == 0) THEN
5168 DO j=1,ist-jb
5169 ij=inder(jb+j) ! index of global parameter
5170 newcounter(ij)=newcounter(ij)+1 ! count again
5171 END DO
5172 ENDIF
5173 END IF
5174 END DO
5175 ! end-of-event
5176 END DO
5177 IF(nr <= 0) EXIT ! end of data?
5178 END DO
5179
5180 ! release read buffer
5185
5186 indab=0
5190 indab=indab+1
5192 ELSE
5194 END IF
5195 END IF
5196 END DO
5201
◆ loop2()
| subroutine loop2 |
Second data loop (number of derivatives, global label pairs).
Calculate maximum number of local, global derivatives and equations per record.
For sparse storage count index pairs with bit (field) counters to construct sparsity structure (row offsets, (compressed) column lists).
Determine read/write cache splitting from average record values (length, global par. vector/matrix).
Check constraints for rank deficit.
Definition at line 5215 of file pede.f90.
5218
5219 IMPLICIT NONE
5220 REAL(mps) :: chin2
5221 REAL(mps) :: chin3
5222 REAL(mps) :: cpr
5223 REAL(mps) :: fsum
5224 REAL(mps) :: gbc
5225 REAL(mps) :: gbu
5226 REAL(mpr8) :: glder
5227 INTEGER(mpi) :: i
5228 INTEGER(mpi) :: ia
5229 INTEGER(mpi) :: ib
5230 INTEGER(mpi) :: ibuf
5231 INTEGER(mpi) :: icgb
5232 INTEGER(mpi) :: iext
5233 INTEGER(mpi) :: ihis
5234 INTEGER(mpi) :: ij
5235 INTEGER(mpi) :: ijn
5236 INTEGER(mpi) :: inder
5237 INTEGER(mpi) :: ioff
5238 INTEGER(mpi) :: iproc
5239 INTEGER(mpi) :: irecmm
5240 INTEGER(mpi) :: isfrst
5241 INTEGER(mpi) :: islast
5242 INTEGER(mpi) :: ist
5243 INTEGER(mpi) :: itgbi
5244 INTEGER(mpi) :: itgbij
5245 INTEGER(mpi) :: itgbik
5246 INTEGER(mpi) :: ivgbij
5247 INTEGER(mpi) :: ivgbik
5248 INTEGER(mpi) :: j
5249 INTEGER(mpi) :: ja
5250 INTEGER(mpi) :: jb
5251 INTEGER(mpi) :: jext
5252 INTEGER(mpi) :: jcgb
5253 INTEGER(mpi) :: jsp
5254 INTEGER(mpi) :: joff
5255 INTEGER(mpi) :: k
5256 INTEGER(mpi) :: kfile
5257 INTEGER(mpi) :: l
5258 INTEGER(mpi) :: label
5259 INTEGER(mpi) :: lu
5260 INTEGER(mpi) :: lun
5261 INTEGER(mpi) :: maeqnf
5262 INTEGER(mpi) :: naeqna
5263 INTEGER(mpi) :: naeqnf
5264 INTEGER(mpi) :: naeqng
5265 INTEGER(mpi) :: nc31
5266 INTEGER(mpi) :: ncachd
5267 INTEGER(mpi) :: ncachi
5268 INTEGER(mpi) :: ncachr
5269 INTEGER(mpi) :: nda
5270 INTEGER(mpi) :: ndf
5271 INTEGER(mpi) :: ndfmax
5272 INTEGER(mpi) :: nfixed
5273 INTEGER(mpi) :: nggd
5274 INTEGER(mpi) :: nggi
5275 INTEGER(mpi) :: nmatmo
5276 INTEGER(mpi) :: noff
5277 INTEGER(mpi) :: nr
5278 INTEGER(mpi) :: nrecf
5279 INTEGER(mpi) :: nrecmm
5280 INTEGER(mpi) :: nst
5281 INTEGER(mpi) :: nwrd
5282 INTEGER(mpi) :: inone
5283 INTEGER(mpi) :: inc
5284 REAL(mps) :: wgh
5285 REAL(mps) :: wolfc3
5286 REAL(mps) :: wrec
5287 REAL(mps) :: chindl
5288
5289 REAL(mpd)::dstat(3)
5290 REAL(mpd)::rerr
5291 INTEGER(mpl):: noff8
5292 INTEGER(mpl):: ndimbi
5293 INTEGER(mpl):: ndimsa(4)
5294 INTEGER(mpl):: ndgn
5295 INTEGER(mpl):: matsiz(2)
5296 INTEGER(mpl):: matwords
5297 INTEGER(mpl):: length
5298 INTEGER(mpl):: rows
5299 INTEGER(mpl):: cols
5300 INTEGER(mpl), PARAMETER :: two=2
5301 INTEGER(mpi) :: maxGlobalPar = 0
5302 INTEGER(mpi) :: maxLocalPar = 0
5303 INTEGER(mpi) :: maxEquations = 0
5304
5305 INTERFACE ! needed for assumed-shape dummy arguments
5308 INTEGER(mpl), DIMENSION(4), INTENT(OUT) :: ndims
5309 INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: ncmprs
5310 INTEGER(mpl), DIMENSION(:,:), INTENT(OUT) :: nsparr
5311 INTEGER(mpi), INTENT(IN) :: ihst
5315 INTEGER(mpl), DIMENSION(:,:), INTENT(IN) :: nsparr
5316 INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: nsparc
5317 INTEGER(mpi), DIMENSION(:), INTENT(IN) :: ncmprs
5321 INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: npair
5323 END INTERFACE
5324
5325 SAVE
5326
5327 !$ INTEGER(mpi) :: OMP_GET_THREAD_NUM
5328
5329 isfrst(ibuf)=readbufferpointer(ibuf)+1
5331 inder(i)=readbufferdatai(i)
5332 glder(i)=readbufferdatad(i)
5333 ! ...
5337
5338 ! two subarrays to get the global parameter indices, used in an event
5339 length=nvgb
5342 backindexusage=0
5343
5344 ! prepare constraints - determine number of constraints NCGB
5345 ! - sort and split into blocks
5347
5353 ELSE ! Lagrange multipliers
5358 ENDIF
5360
5361 ! read all data files and add all variable index pairs -------------
5362
5364
5366 CALL clbits(nagb,mreqpe,mhispe,msngpe,mcmprs,mextnd,ndimbi,nencdb,nspc) ! get dimension for bit storage, encoding, precision info
5367 END IF
5368
5370 length=ntgb
5372 measindex=0
5374 measres=0.0_mps
5375 lunmon=9
5377 ENDIF
5378
5379 ! reading events===reading events===reading events===reading events=
5380 nrecf =0 ! records with fixed global parameters
5381 naeqng=0 ! count number of equations (with global der.)
5382 naeqnf=0 ! count number of equations ( " , fixed)
5383 naeqna=0 ! all
5385 ! monitoring for sparse matrix?
5386 irecmm=0
5388 nmatmo=0
5390 nrecmm=matmon
5391 ELSE
5392 nrecmm=1
5393 END IF
5394 END IF
5395 DO k=1,3
5396 dstat(k)=0.0_mpd
5397 END DO
5398 ! define read buffer
5400 nwrd=nc31+1
5401 length=nwrd*mthrdr
5403 nwrd=nc31*10+2+ndimbuf
5404 length=nwrd*mthrdr
5407 ! to read (old) float binary files
5410
5411 ! for checking appearance
5413 length=5*ntgb
5415 appearancecounter=0
5416 length=ntgb
5418 paircounter=0
5419 END IF
5420
5421 DO
5424 ioff=0
5427 ! Printout for DEBUG
5429 nda=0
5432 WRITE(*,*) ' '
5434 IF (wgh /= 1.0) WRITE(*,*) ' weight ',wrec
5435 ist=isfrst(ibuf)
5436 nst=islast(ibuf)
5437 DO ! loop over measurements
5439 IF(ja == 0) EXIT
5440 nda=nda+1
5443 cycle
5444 END IF
5445 WRITE(*,*) ' '
5446 WRITE(*,*) nda, ' Measured value =',glder(ja),' +- ',glder(jb)
5447 WRITE(*,*) 'Local derivatives:'
5448 WRITE(*,107) (inder(ja+j),glder(ja+j),j=1,jb-ja-1)
5449107 FORMAT(6(i3,g12.4))
5450 IF (jb < ist) THEN
5451 WRITE(*,*) 'Global derivatives:'
5453 globalparlabelindex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb)
5454108 FORMAT(3i11,g12.4)
5455 END IF
5456 IF(nda == 1) THEN
5457 WRITE(*,*) 'total_par_label __label__ var_par_index derivative'
5458 END IF
5459 END DO
5460 WRITE(*,*) ' '
5461 END IF
5462
5466 maeqnf=naeqnf
5467 ist=isfrst(ibuf)
5468 nst=islast(ibuf)
5469 nwrd=nst-ist+1
5470 DO ! loop over measurements
5472 IF(ja == 0.AND.jb == 0) EXIT
5474 naeqna=naeqna+1
5475 IF(ja /= 0) THEN
5476 IF (ist > jb) THEN
5477 naeqng=naeqng+1
5478 ! monitoring, group measurements, sum up entries and errors
5480 rerr =real(glder(jb),mpd) ! the error
5481 ij=inder(jb+1) ! index of first global parameter, used to group measurements
5484 END IF
5485 END IF
5486 nfixed=0
5487 DO j=1,ist-jb
5488 ij=inder(jb+j) ! index of global parameter
5489 ! check appearance
5491 joff = 5*(ij-1)
5494 appearancecounter(joff+1) = kfile
5496 END IF
5498 appearancecounter(joff+3) = kfile
5500 ! count pairs
5501 DO k=1,j
5503 END DO
5504 END IF
5505
5507 IF(ij > 0) THEN
5509 IF(ijn == 0) THEN ! not yet included
5513 END IF
5514 ELSE
5515 nfixed=nfixed+1
5516 END IF
5517 END DO
5518 IF (nfixed > 0) naeqnf=naeqnf+1
5519 END IF
5520
5521 IF(ja /= 0.AND.jb /= 0) THEN
5522 DO j=1,jb-ja-1 ! local parameters
5523 ij=inder(ja+j)
5525 END DO
5526 END IF
5527 END DO
5528
5529 ! end-of-event
5530 IF (naeqnf > maeqnf) nrecf=nrecf+1
5531 irecmm=irecmm+1
5532 ! end-of-event-end-of-event-end-of-event-end-of-event-end-of-event-e
5533
5537
5538 ! sample statistics for caching
5539 dstat(1)=dstat(1)+real((nwrd+2)*2,mpd) ! record size
5542
5544 ! overwrite read buffer with lists of global labels
5545 ioff=ioff+1
5546 readbufferpointer(ibuf)=ioff
5549 iext=globalindexusage(i)
5550 backindexusage(iext)=0
5551 readbufferdatai(ioff+i)=iext
5552 END DO
5553 ioff=ioff+nagbn
5554
5555 END DO
5556 ioff=0
5557
5559 !$OMP PARALLEL &
5560 !$OMP DEFAULT(PRIVATE) &
5561 !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI,MTHRD)
5562 iproc=0
5563 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
5565 ist=isfrst(ibuf)
5566 nst=islast(ibuf)
5567 DO i=ist,nst ! store all combinations
5569 !$ IF (MOD(IEXT,MTHRD).EQ.IPROC) THEN ! distinct rows per thread
5570 DO l=ist,i
5571 jext=readbufferdatai(l)
5573 END DO
5574 !$ ENDIF
5575 END DO
5576 END DO
5577 !$OMP END PARALLEL
5578 ! monitoring
5581 IF (nmatmo == 0) THEN
5582 WRITE(*,*)
5583 WRITE(*,*) 'Monitoring of sparse matrix construction'
5584 WRITE(*,*) ' records ........ off-diagonal elements ', &
5585 '....... compression memory'
5586 WRITE(*,*) ' non-zero used(double) used', &
5587 '(float) [%] [GB]'
5588 END IF
5589 nmatmo=nmatmo+1
5591 gbc=1.0e-9*real((mpi*ndimsa(2)+mpd*ndimsa(3)+mps*ndimsa(4))/mpi*(bit_size(1_mpi)/8),mps) ! GB compressed
5592 gbu=1.0e-9*real(((mpi+mpd)*(ndimsa(3)+ndimsa(4)))/mpi*(bit_size(1_mpi)/8),mps) ! GB uncompressed
5593 cpr=100.0*gbc/gbu
5594 WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc
55951177 FORMAT(i9,3i13,f10.2,f11.6)
5596 DO WHILE(irecmm >= nrecmm)
5598 nrecmm=nrecmm+matmon
5599 ELSE
5600 nrecmm=nrecmm*2
5601 END IF
5602 END DO
5603 END IF
5604
5605 END IF
5606
5607 IF (nr <= 0) EXIT ! next block of events ?
5608 END DO
5609 ! release read buffer
5614
5616 DO k=1,3
5617 dstat(k)=dstat(k)/real(nrec,mpd)
5618 END DO
5619 ! end=of=data=end=of=data=end=of=data=end=of=data=end=of=data=end=of
5620
5623 END IF
5624
5626
5627 ! constraints and index pairs with Lagrange multiplier
5628
5629
5630 ! constraints - determine number of constraints NCGB and index-pairs
5631 ! Lagrange multiplier and global parameters
5632
5633
5635
5636 ! loop over (sorted) constraints
5640 label=listconstraints(i)%label
5641 itgbi=inone(label)
5645 END IF
5646 END DO
5647 END DO
5648
5649 ! measurements - determine index-pairs
5650
5651
5652 i=1
5654 i=i+2
5655 ! loop over label/factor pairs
5656 ia=i
5657 DO
5658 i=i+1
5661 END DO
5662 ib=i-1
5663
5664 DO j=ia,ib
5666 ! first index
5667 ivgbij=0
5669 DO k=ia,j
5671 ! second index
5672 ivgbik=0
5674 IF(ivgbij > 0.AND.ivgbik > 0) THEN
5677 END IF
5678 END DO
5679 END DO
5680
5681 END DO
5682 END IF
5683
5691 END IF
5692 END DO
5695 END IF
5696 ! print numbers ----------------------------------------------------
5697
5699 noff=int(noff8/1000,mpi)
5700 ELSE
5701 noff=int(noff8,mpi)
5702 END IF
5703 ndgn=0
5704 matwords=0
5706 ihis=0
5708 ihis=15
5710 END IF
5713 sparsematrixcompression=0
5717 ndgn=ndimsa(3)+ndimsa(4) ! actual number of off-diagonal elements
5718 matwords=ndimsa(2)+length ! size of sparsity structure
5719
5722 CALL hmpwrt(ihis)
5723 END IF
5724 END IF
5725
5731 ! matrices for event matrices
5732 ! split up cache
5735 fcache(2)=real(dstat(2),mps)
5736 fcache(3)=real(dstat(3),mps)
5737 END IF
5739 DO k=1,3
5741 END DO
5743 ! define read buffer
5745 nwrd=nc31+1
5746 length=nwrd*mthrdr
5748 nwrd=nc31*10+2+ndimbuf
5749 length=nwrd*mthrdr
5752 ! to read (old) float binary files
5755
5764 backindexusage=0
5767 length=nggd*mthrd
5771 length=nggi*mthrd
5773 rows=7; cols=mthrd
5775 rows=2; cols=mthrd
5779
5780 ! print all relevant dimension parameters
5781
5782 DO lu=6,8,2 ! unit 6 and 8
5783
5784 WRITE(*,*) ' '
5786 WRITE(lu,102) '(all parameters, appearing in binary files)'
5788 WRITE(lu,102) '(appearing in fit matrix/vectors)'
5790 WRITE(lu,102) '(including Lagrange multiplier or reduced)'
5793 WRITE(lu,102) '(if =0, no band matrix)'
5795 WRITE(lu,101) 'NOFF/K',noff,'max number of off-diagonal elements'
5796 ELSE
5797 WRITE(lu,101) 'NOFF',noff,'max number of off-diagonal elements'
5798 END IF
5799 IF(ndgn /= 0) THEN
5801 WRITE(lu,101) 'NDGN/K',ndgn/1000,'actual number of off-diagonal elements'
5802 ELSE
5803 WRITE(lu,101) 'NDGN',ndgn,'actual number of off-diagonal elements'
5804 ENDIF
5805 ENDIF
5811 WRITE(lu,101) 'NAEQNA',naeqna,'number of equations'
5812 WRITE(lu,101) 'NAEQNG',naeqng, &
5813 'number of equations with global derivatives'
5814 WRITE(lu,101) 'NAEQNF',naeqnf, &
5815 'number of equations with fixed global derivatives'
5816 WRITE(lu,101) 'NRECF',nrecf, &
5817 'number of records with fixed global derivatives'
5818 END IF
5822111 FORMAT(22x,'cache splitting ',3(f6.1,' %'))
5823 END IF
5824 WRITE(lu,*) ' '
5825
5826 WRITE(lu,*) ' '
5827 WRITE(lu,*) 'Solution method and matrix-storage mode:'
5829 WRITE(lu,*) ' METSOL = 1: matrix inversion'
5831 WRITE(lu,*) ' METSOL = 2: diagonalization'
5837 WRITE(lu,*) ' METSOL = 5: GMRES'
5838 END IF
5841 WRITE(lu,*) ' MATSTO = 1: symmetric matrix, ', '(n*n+n)/2 elements'
5843 WRITE(lu,*) ' MATSTO = 2: sparse matrix'
5844 END IF
5848 END IF
5851 WRITE(lu,*) 'Constraints handled by elimination'
5852 ELSE
5853 WRITE(lu,*) 'Constraints handled by Lagrange multipliers'
5854 ENDIF
5855 END IF
5856
5857 END DO ! print loop
5858
5859 ! Wolfe conditions
5860
5868 wolfc3=wolfc1
5870 wolfc2=wolfc3
5871 ELSE
5872 wolfc1=1.0e-4
5873 wolfc2=0.9
5874 END IF
5877105 FORMAT(' Constants C1, C2 for Wolfe conditions:',g12.4,', ',g12.4)
5878
5879 ! prepare matrix and gradient storage ------------------------------
5880 !32 CONTINUE
5882 matsiz(2)=0 ! number of words for single precision storage '.'
5884 matsiz(1)=ndimsa(3)+nagb
5885 matsiz(2)=ndimsa(4)
5888 END IF
5889 matwords=matwords+matsiz(1)*2+matsiz(2) ! #words for matrix storage
5890
5892
5894 WRITE(*,*) ' '
5895 IF (matwords < 250000) THEN
5896 WRITE(*,*) 'Size of global matrix: < 1 MB'
5897 ELSE
5898 WRITE(*,*) 'Size of global matrix:',int(real(matwords,mps)*4.0e-6,mpi),' MB'
5899 ENDIF
5900 ! print chi^2 cut tables
5901
5902 ndfmax=naeqn-1
5907 ' Chi^2/Ndf(3) Chi^2(3)'
5908 ndf=0
5909 DO
5911 IF(ndf < 10) THEN
5912 ndf=ndf+1
5913 ELSE IF(ndf < 20) THEN
5914 ndf=ndf+2
5915 ELSE IF(ndf < 100) THEN
5916 ndf=ndf+5
5917 ELSE IF(ndf < 200) THEN
5918 ndf=ndf+10
5919 ELSE
5920 EXIT
5921 END IF
5922 chin2=chindl(2,ndf)
5923 chin3=chindl(3,ndf)
5925 END DO
5926
5930101 FORMAT(1x,a8,' =',i10,' = ',a)
5931102 FORMAT(22x,a)
5932103 FORMAT(1x,a,g12.4)
5933106 FORMAT(i6,2(3x,f9.3,f12.1,3x))
subroutine clbits(in, jreqpe, jhispe, jsngpe, jcmprs, jextnd, idimb, iencdb, ispc)
Calculate bit (field) array size, encoding.
Definition mpbits.f90:148
integer(mpi), dimension(:), allocatable globalindexusage
indices of global par in record
Definition mpmod.f90:239
real(mpd), dimension(:), allocatable localglobalmatrix
matrix correlating local and global par
Definition mpmod.f90:256
integer(mpi) matmon
record interval for monitoring of (sparse) matrix construction
Definition mpmod.f90:86
integer(mpi), dimension(:), allocatable paircounter
number of paired parameters (in equations)
Definition mpmod.f90:242
integer(mpi), dimension(:), allocatable meashists
measurement histograms (100 bins per thread)
Definition mpmod.f90:211
integer(mpi), dimension(2) nprecond
number of constraints, matrix size for preconditioner
Definition mpmod.f90:126
real(mps), dimension(3) fcache
read cache, average fill level; write cache; dynamic size
Definition mpmod.f90:91
integer(mpi) icelim
flag for using elimination (instead of multipliers) for constraints
Definition mpmod.f90:102
integer(mpi), dimension(:), allocatable writebufferindices
write buffer for indices
Definition mpmod.f90:260
integer(mpi), dimension(:), allocatable backindexusage
list of global par in record
Definition mpmod.f90:240
real(mpd), dimension(:), allocatable writebufferupdates
write buffer for update matrices
Definition mpmod.f90:261
integer(mpi), dimension(:), allocatable measindex
mapping of 1. global label to measurement index
Definition mpmod.f90:210
integer(mpi), dimension(-6:6) writebufferheader
write buffer header (-6..-1: updates, 1..6: indices)
Definition mpmod.f90:262
integer(mpi), dimension(:), allocatable appearancecounter
appearance statistics for global par (first/last file,record)
Definition mpmod.f90:241
real(mps), dimension(:,:), allocatable writebufferdata
write buffer data (largest residual, Chi2/ndf, per thread)
Definition mpmod.f90:259
integer(mpi), dimension(:,:), allocatable writebufferinfo
write buffer management (per thread)
Definition mpmod.f90:258
integer(mpi), dimension(:), allocatable sparsematrixcompression
compression info (per row)
Definition mpmod.f90:229
integer(mpi) imonit
flag for monitoring residuals per local fit cycle (=0: none, <0: all, bit 0: first,...
Definition mpmod.f90:106
◆ loopbf()
| subroutine loopbf | ( | integer(mpi), dimension(0:3), intent(inout) | nrej, |
| integer(mpi), intent(inout) | ndfs, | ||
| real(mpd), intent(inout) | sndf, | ||
| real(mpd), intent(inout) | dchi2s, | ||
| integer(mpi), intent(in) | numfil, | ||
| integer(mpi), dimension(numfil), intent(inout) | naccf, | ||
| real(mps), dimension(numfil), intent(inout) | chi2f, | ||
| integer(mpi), dimension(numfil), intent(inout) | ndff | ||
| ) |
Loop over records in read buffer (block), fits and sums.
Loop over records in current read buffer block (with multiple threads). Perform local fits (optionally with outlier downweigthing) to calculate Chi2, ndf and r.h.s. 'b' of linear equation system A*x=b. In first iteration(s) fill global matrix A.
For the filling of the global matrix each thread creates from his share of local fits (small) udpdate matrices ( 
The matrices of the local fits are checked for bordered band structure. For border size b and band width m all elements (i,j) are zero for min(i,j)>b and abs(i-j)>m. For sufficient small (b,m) a solution by root free Cholesky decomposition and forward/backward substitution of the band part is much faster compared to inversion (see broken lines in references). Based on the expected computing cost the faster solution method is selected.
- Parameters
-
[in,out] nrej number of rejected records [in,out] ndfs sum(ndf) [in,out] sndf sum(weighted ndf) [in,out] dchi2s sum(weighted chi2) [in] numfil number of binary files [in,out] naccf number of accepted records per binary file [in,out] chi2f sum(chi2/ndf) per binary file [in,out] ndff sum(ndf) per binary file
Definition at line 3063 of file pede.f90.
3065
3066 IMPLICIT NONE
3067 REAL(mpd) :: cauchy
3068 REAL(mps) :: chichi
3069 REAL(mps) :: chlimt
3070 REAL(mps) :: chndf
3071 REAL(mpd) :: chuber
3072 REAL(mpd) :: down
3073 REAL(mpr8) :: glder
3074 REAL(mpd) :: pull
3075 REAL(mpd) :: r1
3076 REAL(mpd) :: r2
3077 REAL(mps) :: rec
3078 REAL(mpd) :: rerr
3079 REAL(mpd) :: resid
3080 REAL(mps) :: resing
3081 REAL(mpd) :: resmax
3082 REAL(mpd) :: rmeas
3083 REAL(mpd) :: rmloc
3084 REAL(mpd) :: suwt
3085 REAL(mps) :: used
3086 REAL(mpd) :: wght
3087 REAL(mps) :: chindl
3088 INTEGER(mpi) :: i
3089 INTEGER(mpi) :: ia
3090 INTEGER(mpi) :: ib
3091 INTEGER(mpi) :: ibuf
3092 INTEGER(mpi) :: ichunk
3093 INTEGER(mpl) :: icmn
3094 INTEGER(mpl) :: icost
3095 INTEGER(mpi) :: id
3096 INTEGER(mpi) :: idiag
3097 INTEGER(mpi) :: iext
3098 INTEGER(mpi) :: ij
3099 INTEGER(mpi) :: ije
3100 INTEGER(mpi) :: ijn
3101 INTEGER(mpi) :: ijsym
3102 INTEGER(mpi) :: ik
3103 INTEGER(mpi) :: ike
3104 INTEGER(mpi) :: im
3105 INTEGER(mpi) :: imeas
3106 INTEGER(mpi) :: in
3107 INTEGER(mpi) :: inder
3108 INTEGER(mpi) :: inv
3109 INTEGER(mpi) :: ioffb
3110 INTEGER(mpi) :: ioffc
3111 INTEGER(mpi) :: ioffd
3112 INTEGER(mpi) :: ioffe
3113 INTEGER(mpi) :: ioffi
3114 INTEGER(mpi) :: iprdbg
3115 INTEGER(mpi) :: iproc
3116 INTEGER(mpi) :: irbin
3117 INTEGER(mpi) :: irow
3118 INTEGER(mpi) :: isfrst
3119 INTEGER(mpi) :: islast
3120 INTEGER(mpi) :: ist
3121 INTEGER(mpi) :: iter
3122 INTEGER(mpi) :: itgbi
3123 INTEGER(mpi) :: ivgbj
3124 INTEGER(mpi) :: ivgbk
3125 INTEGER(mpi) :: j
3126 INTEGER(mpi) :: ja
3127 INTEGER(mpi) :: jb
3128 INTEGER(mpi) :: jk
3129 INTEGER(mpi) :: jn
3130 INTEGER(mpi) :: joffd
3131 INTEGER(mpi) :: joffi
3132 INTEGER(mpi) :: jproc
3133 INTEGER(mpi) :: jsp
3134 INTEGER(mpi) :: k
3135 INTEGER(mpi) :: kbdr
3136 INTEGER(mpi) :: kbdrx
3137 INTEGER(mpi) :: kbnd
3138 INTEGER(mpi) :: kfl
3139 INTEGER(mpi) :: kx
3140 INTEGER(mpi) :: mbdr
3141 INTEGER(mpi) :: mbnd
3142 INTEGER(mpi) :: mside
3143 INTEGER(mpi) :: nalc
3144 INTEGER(mpi) :: nalg
3145 INTEGER(mpi) :: nan
3146 INTEGER(mpi) :: nb
3147 INTEGER(mpi) :: ndf
3148 INTEGER(mpi) :: ndfsd
3149 INTEGER(mpi) :: ndown
3150 INTEGER(mpi) :: neq
3151 INTEGER(mpi) :: nfred
3152 INTEGER(mpi) :: nfrei
3153 INTEGER(mpi) :: ngg
3154 INTEGER(mpi) :: nprdbg
3155 INTEGER(mpi) :: nrank
3156 INTEGER(mpi) :: nrc
3157 INTEGER(mpi) :: nst
3158 INTEGER(mpi) :: nter
3159 INTEGER(mpi) :: nweig
3160
3161 INTEGER(mpi), INTENT(IN OUT) :: nrej(0:3)
3162 INTEGER(mpi), INTENT(IN OUT) :: ndfs
3163 REAL(mpd), INTENT(IN OUT) :: sndf
3164 REAL(mpd), INTENT(IN OUT) :: dchi2s
3165 INTEGER(mpi), INTENT(IN) :: numfil
3166 INTEGER(mpi), INTENT(IN OUT) :: naccf(numfil)
3167 REAL(mps), INTENT(IN OUT) :: chi2f(numfil)
3168 INTEGER(mpi), INTENT(IN OUT) :: ndff(numfil)
3169
3170 REAL(mpd):: dchi2
3171 REAL(mpd)::dvar
3172 REAL(mpd):: dw1
3173 REAL(mpd)::dw2
3174 REAL(mpd)::summ
3175
3176 !$ INTEGER(mpi) OMP_GET_THREAD_NUM
3177
3178 LOGICAL:: lprnt
3179 LOGICAL::lhist
3180 CHARACTER (LEN=3):: chast
3181 DATA chuber/1.345_mpd/ ! constant for Huber down-weighting
3182 DATA cauchy/2.3849_mpd/ ! constant for Cauchy down-weighting
3183 SAVE chuber,cauchy
3184 ! ...
3185 ijsym(i,j)=min(i,j)+(max(i,j)*max(i,j)-max(i,j))/2
3186 isfrst(ibuf)=readbufferpointer(ibuf)+1
3188 inder(i)=readbufferdatai(i)
3189 glder(i)=readbufferdatad(i)
3190
3192 ! reset header, 3 words per thread:
3193 ! number of entries, offset to data, indices
3194 writebufferinfo=0
3195 writebufferdata=0.
3196 nprdbg=0
3197 iprdbg=-1
3198
3199 ! parallelize record loop
3200 ! private copy of NDFS,.. for each thread, combined at end, init with 0.
3201 !$OMP PARALLEL DO &
3202 !$OMP DEFAULT(PRIVATE) &
3203 !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI, &
3204 !$OMP readBufferDataD,writeBufferHeader,writeBufferInfo, &
3205 !$OMP writeBufferData,writeBufferIndices,writeBufferUpdates,globalVector,globalCounter, &
3206 !$OMP globalParameter,globalParLabelIndex,globalIndexUsage,backIndexUsage, &
3207 !$OMP measBins,numMeas,measIndex,measRes,measHists, &
3208 !$OMP NAGB,NVGB,NAGBN,ICALCM,ICHUNK,NLOOPN,NRECER,NPRDBG,IPRDBG, &
3209 !$OMP NEWITE,CHICUT,LHUBER,CHUBER,ITERAT,NRECPR,MTHRD, &
3210 !$OMP DWCUT,CHHUGE,NRECP2,CAUCHY,LFITNP,LFITBB,IMONIT,IMONMD) &
3211 !$OMP REDUCTION(+:NDFS,SNDF,DCHI2S,NREJ,NBNDR,NACCF,CHI2F,NDFF) &
3212 !$OMP REDUCTION(MAX:NBNDX,NBDRX) &
3213 !$OMP REDUCTION(MIN:NREC3) &
3214 !$OMP SCHEDULE(DYNAMIC,ICHUNK)
3219 dw2=sqrt(dw1)
3220
3221 iproc=0
3222 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
3228 ! ----- reset ------------------------------------------------------
3229 lprnt=.false.
3230 lhist=(iproc == 0)
3231 rec=nrc ! floating point value
3233 WRITE(*,*) 'Record',nrc,' ... still reading'
3234 END IF
3235
3236 ! printout/debug only for one thread at a time
3237
3238
3239 ! flag for record printout -----------------------------------------
3240
3241 lprnt=.false.
3246 END IF
3247 IF (lprnt)THEN
3248 !$OMP ATOMIC
3249 nprdbg=nprdbg+1 ! number of threads with debug
3250 IF (nprdbg == 1) iprdbg=iproc ! first thread with debug
3251 IF (iproc /= iprdbg) lprnt=.false.
3252 ! print *, ' LPRNT ', NRC, NPRDBG, IPRDBG, IPROC, LPRNT
3253 END IF
3254 IF(lprnt) THEN
3255 WRITE(1,*) ' '
3257 ': Printout for record',nrc,iproc
3258 WRITE(1,*) ' '
3259 END IF
3260
3261 ! ----- print data -------------------------------------------------
3262
3263 IF(lprnt) THEN
3264 imeas=0 ! local derivatives
3265 ist=isfrst(ibuf)
3266 nst=islast(ibuf)
3267 DO ! loop over measurements
3269 IF(ja == 0) EXIT
3270 IF(imeas == 0) WRITE(1,1121)
3271 imeas=imeas+1
3272 WRITE(1,1122) imeas,glder(ja),glder(jb), &
3273 (inder(ja+j),glder(ja+j),j=1,jb-ja-1)
3274 END DO
32751121 FORMAT(/'Measured value and local derivatives'/ &
3276 ' i measured std_dev index...derivative ...')
32771122 FORMAT(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
3278
3279 imeas=0 ! global derivatives
3280 ist=isfrst(ibuf)
3281 nst=islast(ibuf)
3282 DO ! loop over measurements
3284 IF(ja == 0) EXIT
3285 IF(imeas == 0) WRITE(1,1123)
3286 imeas=imeas+1
3287 IF (jb < ist) THEN
3288 IF(ist-jb > 2) THEN
3290 globalparlabelindex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb)
3291 ELSE
3293 globalparlabelindex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb)
3294 END IF
3295 END IF
3296 END DO
32971123 FORMAT(/'Global derivatives'/ &
3298 ' i label gindex vindex derivative ...')
32991124 FORMAT(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
33001125 FORMAT(i3,2(i9,i7,i7,g12.4))
3301 END IF
3302
3303 ! ----- first loop -------------------------------------------------
3304 ! ------ prepare local fit ------
3305 ! count local and global derivates
3306 ! subtract actual alignment parameters from the measured data
3307
3308 IF(lprnt) THEN
3309 WRITE(1,*) ' '
3310 WRITE(1,*) 'Data corrections using values of global parameters'
3311 WRITE(1,*) '=================================================='
3312 WRITE(1,101)
3313 END IF
3314 nalg=0 ! count number of global derivatives
3315 nalc=0 ! count number of local derivatives
3316 neq=0 ! count number of equations
3317 ist=isfrst(ibuf)
3318 nst=islast(ibuf)
3319 DO ! loop over measurements
3321 IF(ja == 0) EXIT
3322 rmeas=real(glder(ja),mpd) ! data
3323 neq=neq+1 ! count equation
3324 ! subtract global ... from measured value
3325 DO j=1,ist-jb ! global parameter loop
3326 itgbi=inder(jb+j) ! global parameter label
3330 IF(ij > 0) THEN
3332 IF(ijn == 0) THEN ! not yet included
3333 nalg=nalg+1 ! count
3336 END IF
3337 END IF
3338 END IF
3339 END DO
3340 IF(lprnt) THEN
3341 IF (jb < ist) WRITE(1,102) neq,glder(ja),rmeas,glder(jb)
3342 END IF
3344 DO j=1,jb-ja-1 ! local parameter loop
3345 ij=inder(ja+j)
3346 nalc=max(nalc,ij) ! number of local parameters
3347 END DO
3348 END DO
3349101 FORMAT(' index measvalue corrvalue sigma')
3350102 FORMAT(i6,2x,2g12.4,' +-',g12.4)
3351
3352 IF(nalc <= 0) GO TO 90
3353
3354 ngg=(nalg*nalg+nalg)/2
3356 DO k=1,nalg*nalc
3358 END DO
3362 DO k=1,nalg
3363 iext=globalindexusage(ioffc+k)
3366 END DO
3367 DO k=1,ngg
3369 END DO
3370 END IF
3371 ! ----- iteration start and check ---------------------------------
3372
3373 nter=1 ! first loop without down-weighting
3375 localcorrections(1:neq) = 0._mpd
3376
3377 ! check matrix for bordered band structure (MBDR+MBND+1 <= NALC)
3378 mbnd=-1
3379 mbdr=nalc
3380 mside=-1 ! side (1: upper/left border, 2: lower/right border)
3381 DO i=1, 2*nalc
3382 ibandh(i)=0
3383 END DO
3384 irow=1
3385 idiag=1
3386 ndfsd=0
3387
3388 iter=0
3389 resmax=0.0
3390 DO WHILE(iter < nter) ! outlier suppresssion iteration loop
3391 iter=iter+1
3392 resmax=0.0
3393 IF(lprnt) THEN
3394 WRITE(1,*) ' '
3395 WRITE(1,*) 'Outlier-suppression iteration',iter,' of',nter
3396 WRITE(1,*) '=========================================='
3397 WRITE(1,*) ' '
3398 imeas=0
3399 END IF
3400
3401 ! ----- second loop ------------------------------------------------
3402 ! accumulate normal equations for local fit and determine solution
3403 DO i=1,nalc
3405 END DO
3406 DO i=1,(nalc*nalc+nalc)/2 ! GF: FIXME - not really, local parameter number...
3408 END DO
3409 neq=0
3410 ndown=0
3411 nweig=0
3412 ist=isfrst(ibuf)
3413 nst=islast(ibuf)
3414 DO ! loop over measurements
3416 IF(ja == 0) EXIT
3417 rmeas=real(glder(ja),mpd) ! data
3418 rerr =real(glder(jb),mpd) ! ... and the error
3419 wght =1.0_mpd/rerr**2 ! weight from error
3420 neq=neq+1 ! count equation
3421 nweig=nweig+1
3424 IF(iter <= 3) THEN
3425 IF(abs(resid) > chuber*rerr) THEN ! down-weighting
3426 wght=wght*chuber*rerr/abs(resid)
3427 ndown=ndown+1
3428 END IF
3429 ELSE ! Cauchy
3430 wght=wght/(1.0+(resid/rerr/cauchy)**2)
3431 END IF
3432 END IF
3433
3434 IF(lprnt.AND.iter /= 1.AND.nter /= 1) THEN
3435 chast=' '
3436 IF(abs(resid) > chuber*rerr) chast='* '
3437 IF(abs(resid) > 3.0*rerr) chast='** '
3438 IF(abs(resid) > 6.0*rerr) chast='***'
3439 IF(imeas == 0) WRITE(1,*) 'Second loop: accumulate'
3440 IF(imeas == 0) WRITE(1,103)
3441 imeas=imeas+1
3442 down=1.0/sqrt(wght)
3443 r1=resid/rerr
3444 r2=resid/down
3445 WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
3446 END IF
3447103 FORMAT(' index corrvalue residuum sigma', &
3448 ' nresid cnresid')
3449104 FORMAT(i6,2x,2g12.4,' +-',g12.4,f7.2,1x,a3,f8.2)
3450
3451 DO j=1,jb-ja-1 ! normal equations, local parameter loop
3452 ij=inder(ja+j) ! local parameter index J
3454 DO k=1,j
3455 ik=inder(ja+k) ! local parameter index K
3456 jk=ijsym(ij,ik) ! index in symmetric matrix
3458 +wght*real(glder(ja+j),mpd)*real(glder(ja+k),mpd)
3459 ! check for band matrix substructure
3460 IF (iter == 1) THEN
3461 id=iabs(ij-ik)+1
3462 im=min(ij,ik) ! upper/left border
3464 im=min(nalc+1-ij,nalc+1-ik) ! lower/rght border (mirrored)
3466 END IF
3467 END DO
3468 END DO
3469 END DO
3470 ! for non trivial fits check for bordered band matrix structure
3472 kx=-1
3473 kbdrx=0
3474 icmn=int(nalc,mpl)**3 ! cost (*6) should improve by at least factor 2
3475 ! upper/left border ?
3476 kbdr=0
3477 DO k=nalc,2,-1
3478 kbnd=k-2
3479 kbdr=max(kbdr,ibandh(k))
3480 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
3481 IF (icost < icmn) THEN
3482 icmn=icost
3483 kx=k
3484 kbdrx=kbdr
3485 mside=1
3486 END IF
3487 END DO
3488 IF (kx < 0) THEN
3489 ! lower/right border instead?
3490 kbdr=0
3491 DO k=nalc,2,-1
3492 kbnd=k-2
3493 kbdr=max(kbdr,ibandh(k+nalc))
3494 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
3495 IF (icost < icmn) THEN
3496 icmn=icost
3497 kx=k
3498 kbdrx=kbdr
3499 mside=2
3500 END IF
3501 END DO
3502 END IF
3503 IF (kx > 0) THEN
3504 mbnd=kx-2
3505 mbdr=kbdrx
3506 END IF
3507 END IF
3508
3509 IF (mbnd >= 0) THEN
3510 ! fast solution for border banded matrix (inverse for ICALCM>0)
3515 END IF
3516
3517 inv=0
3520 IF (mside == 1) THEN
3523 ELSE
3526 ENDIF
3527 ELSE
3528 ! full inversion and solution
3529 inv=2
3531 END IF
3532 ! check for NaNs
3533 nan=0
3534 DO k=1, nalc
3536 END DO
3537
3538 IF(lprnt) THEN
3539 WRITE(1,*) ' '
3540 WRITE(1,*) 'Parameter determination:',nalc,' parameters,', ' rank=',nrank
3541 WRITE(1,*) '-----------------------'
3542 IF(ndown /= 0) WRITE(1,*) ' ',ndown,' data down-weighted'
3543 WRITE(1,*) ' '
3544 END IF
3545
3546 ! ----- third loop -------------------------------------------------
3547 ! calculate single residuals remaining after local fit and chi^2
3548
3549 summ=0.0_mpd
3550 suwt=0.0
3551 neq=0
3552 imeas=0
3553 ist=isfrst(ibuf)
3554 nst=islast(ibuf)
3555 DO ! loop over measurements
3557 IF(ja == 0) EXIT
3558 rmeas=real(glder(ja),mpd) ! data (global contrib. subtracted)
3559 rerr =real(glder(jb),mpd) ! ... and the error
3560 wght =1.0_mpd/rerr**2 ! weight from error
3561 neq=neq+1 ! count equation
3562 rmloc=0.0 ! local fit result reset
3563 DO j=1,jb-ja-1 ! local parameter loop
3564 ij=inder(ja+j)
3566 END DO
3568 rmeas=rmeas-rmloc ! reduced to residual
3569
3570 ! calculate pulls? (needs covariance matrix)
3572 dvar=0.0_mpd
3573 DO j=1,jb-ja-1
3574 ij=inder(ja+j)
3575 DO k=1,jb-ja-1
3576 ik=inder(ja+k)
3577 jk=ijsym(ij,ik)
3578 dvar=dvar+clmat(jk)*real(glder(ja+j),mpd)*real(glder(ja+k),mpd)
3579 END DO
3580 END DO
3581 ! some variance left to define a pull?
3582 IF (0.999999_mpd/wght > dvar) THEN
3583 pull=rmeas/sqrt(1.0_mpd/wght-dvar)
3584 IF (lhist) THEN
3585 IF (jb < ist) THEN
3586 CALL hmpent(13,real(pull,mps)) ! histogram pull
3587 CALL gmpms(5,rec,real(pull,mps))
3588 ELSE
3589 CALL hmpent(14,real(pull,mps)) ! histogram pull
3590 END IF
3591 END IF
3592 ! monitoring
3594 IF (jb < ist) THEN
3595 ij=inder(jb+1) ! group by first global label
3598 ELSE
3600 ENDIF
3603 ENDIF
3604 ENDIF
3605 END IF
3606 END IF
3607
3608 IF(iter == 1.AND.jb < ist.AND.lhist) &
3609 CALL gmpms(4,rec,real(rmeas/rerr,mps)) ! residual (with global deriv.)
3610
3611 dchi2=wght*rmeas*rmeas
3612 ! DCHIT=DCHI2
3613 resid=rmeas
3615 IF(iter <= 3) THEN
3616 IF(abs(resid) > chuber*rerr) THEN ! down-weighting
3617 wght=wght*chuber*rerr/abs(resid)
3618 dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
3619 END IF
3620 ELSE
3621 wght=wght/(1.0_mpd+(resid/rerr/cauchy)**2)
3622 dchi2=log(1.0_mpd+(resid/rerr/cauchy)**2)*cauchy**2
3623 END IF
3624 END IF
3625
3626 down=1.0/sqrt(wght)
3627
3628 ! SUWT=SUWT+DCHI2/DCHIT
3629 suwt=suwt+rerr/down
3630 IF(lprnt) THEN
3631 chast=' '
3632 IF(abs(resid) > chuber*rerr) chast='* '
3633 IF(abs(resid) > 3.0*rerr) chast='** '
3634 IF(abs(resid) > 6.0*rerr) chast='***'
3635 IF(imeas == 0) WRITE(1,*) 'Third loop: single residuals'
3636 IF(imeas == 0) WRITE(1,105)
3637 imeas=imeas+1
3638 r1=resid/rerr
3639 r2=resid/down
3640 IF(resid < 0.0) r1=-r1
3641 IF(resid < 0.0) r2=-r2
3642 WRITE(1,106) imeas,glder(ja),rmeas,rerr,r1,chast,r2
3643 END IF
3644105 FORMAT(' index corrvalue residuum sigma', &
3645 ' nresid cnresid')
3646106 FORMAT(i6,2x,2g12.4,' +-',g12.4,f7.2,1x,a3,f8.2)
3647
3648 IF(iter == nter) THEN
3650 resmax=max(resmax,abs(rmeas)/rerr)
3651 END IF
3652
3653 IF(iter == 1.AND.lhist) THEN
3654 IF (jb < ist) THEN
3655 CALL hmpent( 3,real(rmeas/rerr,mps)) ! histogram norm residual
3656 ELSE
3657 CALL hmpent(12,real(rmeas/rerr,mps)) ! histogram norm residual
3658 END IF
3659 END IF
3660 summ=summ+dchi2 ! accumulate chi-square sum
3661 END DO
3662 ndf=neq-nrank
3663 resing=(real(nweig,mps)-real(suwt,mps))/real(nweig,mps)
3664 IF (lhist) THEN
3665 IF(iter == 1) CALL hmpent( 5,real(ndf,mps)) ! histogram Ndf
3666 IF(iter == 1) CALL hmpent(11,real(nalc,mps)) ! histogram Nlocal
3668 END IF
3669 IF(lprnt) THEN
3670 WRITE(1,*) ' '
3671 WRITE(1,*) 'Chi^2=',summ,' at',ndf,' degrees of freedom: ', &
3673 WRITE(1,*) suwt,' is sum of factors, compared to',nweig, &
3674 ' Downweight fraction:',resing
3675 END IF
3676 IF(nrank /= nalc.OR.nan > 0) THEN
3677 nrej(0)=nrej(0)+1 ! count cases
3679 IF(lprnt) THEN
3680 WRITE(1,*) ' rank deficit/NaN ', nalc, nrank, nan
3681 WRITE(1,*) ' ---> rejected!'
3682 END IF
3683 GO TO 90
3684 END IF
3685 IF(ndf <= 0) THEN
3686 nrej(1)=nrej(1)+1 ! count cases
3687 IF(lprnt) THEN
3688 WRITE(1,*) ' ---> rejected!'
3689 END IF
3690 GO TO 90
3691 END IF
3692
3693 chndf=real(summ/real(ndf,mpd),mps)
3694
3695 IF(iter == 1.AND.lhist) CALL hmpent(4,chndf) ! histogram chi^2/Ndf
3696 END DO ! outlier iteration loop
3697
3698 ndfs=ndfs+ndf ! (local) sum of Ndf
3699 sndf=sndf+real(ndf,mpd)*dw1 ! (local) weighted sum of Ndf
3700
3701 ! ----- reject eventually ------------------------------------------
3702
3705 writebufferdata(2,iproc+1)=chndf
3706 writebufferinfo(7,iproc+1)=nrc
3707 END IF
3708 END IF
3709
3710 chichi=chindl(3,ndf)*real(ndf,mps)
3711 ! GF IF(SUMM.GT.50.0*CHICHI) THEN ! huge
3712 ! CHK CHICUT<0: NO cut (1st iteration)
3715 nrej(2)=nrej(2)+1 ! count cases with huge chi^2
3716 IF(lprnt) THEN
3717 WRITE(1,*) ' ---> rejected!'
3718 END IF
3719 GO TO 90
3720 END IF
3721
3723 chlimt=chicut*chichi
3724 ! WRITE(*,*) 'chi^2 ',SUMM,CHLIMT,CHICUT,CHINDL(3,NDF),NDF
3725 IF(summ > chlimt) THEN
3726 IF(lprnt) THEN
3727 WRITE(1,*) ' ---> rejected!'
3728 END IF
3729 ! add to FVALUE
3730 dchi2=chlimt ! total contribution limit
3731 dchi2s=dchi2s+dchi2*dw1 ! add total contribution
3732 nrej(3)=nrej(3)+1 ! count cases with large chi^2
3733 GO TO 90
3734 END IF
3735 END IF
3736 END IF
3737
3739 ! add to FVALUE
3740 dchi2=summ ! total contribution
3741 dchi2s=dchi2s+dchi2*dw1 ! add total contribution
3742 nrej(3)=nrej(3)+1 ! count cases with large chi^2
3743 ! WRITE(*,*) 'Downweight fraction cut ',RESING,DWCUT,SUMM
3744 IF(lprnt) THEN
3745 WRITE(1,*) ' ---> rejected!'
3746 END IF
3747 GO TO 90
3748 END IF
3749
3752 writebufferdata(1,iproc+1)=real(resmax,mps)
3753 writebufferinfo(6,iproc+1)=nrc
3754 END IF
3755 END IF
3756 ! 'track quality' per binary file: accepted records
3757 naccf(kfl)=naccf(kfl)+1
3758 ndff(kfl) =ndff(kfl) +ndf
3759 chi2f(kfl)=chi2f(kfl)+chndf
3760
3761 ! ----- fourth loop ------------------------------------------------
3762 ! update of global matrix and vector according to the "Millepede"
3763 ! principle, from the global/local information
3764
3765 ist=isfrst(ibuf)
3766 nst=islast(ibuf)
3767 DO ! loop over measurements
3769 IF(ja <= 0) EXIT
3770
3771 rmeas=real(glder(ja),mpd) ! data residual
3772 rerr =real(glder(jb),mpd) ! ... and the error
3773 wght =1.0_mpd/rerr**2 ! weight from measurement error
3774 dchi2=wght*rmeas*rmeas ! least-square contribution
3775
3777 resid=abs(rmeas)
3778 IF(resid > chuber*rerr) THEN
3779 wght=wght*chuber*rerr/resid ! down-weighting
3780 dchi2=2.0*chuber*(resid/rerr-0.5*chuber) ! modified contribution
3781 END IF
3782 END IF
3783 dchi2s=dchi2s+dchi2*dw1 ! add to total objective function
3784
3785 ! global-global matrix contribution: add directly to gg-matrix
3786
3787 DO j=1,ist-jb
3789 IF(ivgbj > 0) THEN
3791 +dw1*wght*rmeas*real(glder(jb+j),mpd) ! vector !!! reverse
3795 DO k=1,j
3796 ivgbk=globalparlabelindex(2,inder(jb+k))
3797 IF(ivgbk > 0) THEN
3799 ia=max(ije,ike) ! larger
3800 ib=min(ije,ike) ! smaller
3801 ij=ib+(ia*ia-ia)/2
3803 -dw1*wght*real(glder(jb+j),mpd)*real(glder(jb+k),mpd)
3804 END IF
3805 END DO
3806 END IF
3807 END IF
3808 END DO
3809
3810 ! normal equations - rectangular matrix for global/local pars
3811 ! global-local matrix contribution: accumulate rectangular matrix
3813 DO j=1,ist-jb
3815 IF(ivgbj > 0) THEN
3817 DO k=1,jb-ja-1
3818 ik=inder(ja+k) ! local index
3819 jk=ik+(ije-1)*nalc ! matrix index
3820 localglobalmatrix(jk)=localglobalmatrix(jk)+dw2*wght*real(glder(jb+j),mpd)*real(glder(ja+k),mpd)
3821 END DO
3822 END IF
3823 END DO
3824 END DO
3825
3826
3827 ! ----- final matrix update ----------------------------------------
3828 ! update global matrices and vectors
3830 ! (inverse local matrix) * (rectang. matrix) -> CORM
3831 ! T
3832 ! resulting symmetrix matrix = G * Gamma^{-1} * G
3833
3835
3836 ! (rectang. matrix) * (local param vector) -> CORV
3837 ! resulting vector = G * q (q = local parameter)
3838 ! CALL DBGAX(DQ(IGLMA/2+1),BLVEC,DQ(ICORV/2+1),NALG,NALC) ! not done
3839 ! the vector update is not done, because after local fit it is zero!
3840
3841 ! update cache status
3845 ! check free space
3848 IF (nfred < 0.OR.nfrei < 0) THEN ! need to flush
3849 nb=writebufferinfo(1,iproc+1)
3856 !$OMP CRITICAL
3859
3860 DO ib=1,nb
3861 ijn=0
3863 i=writebufferindices(joffi+in)
3864 ! DQ(IGVEC/2+I)=DQ(IGVEC/2+I)+DQ(ICORV/2+IN) ! not done: = zero
3865 DO jn=1,in
3866 ijn=ijn+1
3867 j=writebufferindices(joffi+jn)
3869 END DO
3870 END DO
3871 joffd=joffd+ijn
3872 joffi=joffi+writebufferindices(joffi)+2
3873 END DO
3874 !$OMP END CRITICAL
3875 ! reset counter, pointers
3876 DO k=1,3
3877 writebufferinfo(k,iproc+1)=0
3878 END DO
3879 END IF
3880
388190 IF(lprnt) THEN
3882 WRITE(1,*) ' '
3883 WRITE(1,*) '------------------ End of printout for record',nrc
3884 WRITE(1,*) ' '
3885 END IF
3886
3887 DO i=1,nalg ! reset global index array
3888 iext=globalindexusage(ioffc+i)
3889 backindexusage(ioffe+iext)=0
3890 END DO
3891
3892 END DO
3893 !$OMP END PARALLEL DO
3894
3896 ! flush remaining matrices
3903 writebufferinfo(4,k)=0
3909 writebufferinfo(5,k)=0
3910 END DO
3911
3912 !$OMP PARALLEL &
3913 !$OMP DEFAULT(PRIVATE) &
3914 !$OMP SHARED(writeBufferHeader,writeBufferInfo,writeBufferIndices,writeBufferUpdates,MTHRD)
3915 iproc=0
3916 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
3918 nb=writebufferinfo(1,jproc+1)
3919 ! print *, ' flush end ', JPROC, NRC, NB
3922 DO ib=1,nb
3923 ! print *, ' buf end ', JPROC,IB,writeBufferIndices(JOFFI-1),writeBufferIndices(JOFFI)
3924 ijn=0
3926 i=writebufferindices(joffi+in)
3927 !$ IF (MOD(I,MTHRD).EQ.IPROC) THEN
3928 DO jn=1,in
3929 ijn=ijn+1
3930 j=writebufferindices(joffi+jn)
3932 END DO
3933 !$ ELSE
3934 !$ IJN=IJN+IN
3935 !$ ENDIF
3936 END DO
3937 joffd=joffd+ijn
3938 joffi=joffi+writebufferindices(joffi)+2
3939 END DO
3940 END DO
3941 !$OMP END PARALLEL
3942 END IF
3943
3950 END IF
3951 END DO
3952 END IF
3958 END IF
3959 END DO
3960 END IF
3961 END IF
3962
subroutine sqmibb2(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Band bordered matrix.
Definition mpnum.f90:2611
subroutine sqmibb(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Bordered band matrix.
Definition mpnum.f90:2355
integer(mpi) nrecer
record with error (rank deficit or Not-a-Number) for printout
Definition mpmod.f90:77
real(mpd), dimension(:), allocatable vbk
local fit 'matrix for border solution'
Definition mpmod.f90:251
real(mps) chhuge
cut in terms of 3-sigma for unreasonable data, all iterations
Definition mpmod.f90:49
real(mpd), dimension(:), allocatable localcorrections
local fit corrections (to residuals)
Definition mpmod.f90:255
integer(mpi), dimension(:), allocatable globalcounter
global counter (entries in 'x')
Definition mpmod.f90:186
integer(mpi), dimension(2) nbndr
number of records with bordered band matrix for local fit (upper/left, lower/right)
Definition mpmod.f90:74
integer(mpi) imonmd
monitoring mode: 0:residuals (normalized to average error), 1:pulls
Definition mpmod.f90:108
real(mpd), dimension(:), allocatable aux
local fit 'solutions for border rows'
Definition mpmod.f90:250
integer(mpi), dimension(:), allocatable ibandh
local fit 'band width histogram' (band size autodetection)
Definition mpmod.f90:246
real(mpd), dimension(:), allocatable blvec
local fit vector 'b' (in A*x=b), replaced by 'x'
Definition mpmod.f90:244
◆ loopn()
| subroutine loopn |
Loop with fits and sums.
Loop over all binary files. Perform local fits to calculate Chi2, ndf and r.h.s. 'b' of linear equation system A*x=b. In first iteration(s) fill matrix A.
Definition at line 2315 of file pede.f90.
2317
2318 IMPLICIT NONE
2319 REAL(mpd) :: dsum
2320 REAL(mps) :: elmt
2321 REAL(mpd) :: factrj
2322 REAL(mpd) :: factrk
2323 REAL(mpr8) :: glder
2324 REAL(mps) :: peakd
2325 REAL(mps) :: peaki
2326 REAL(mps) :: ratae
2327 REAL(mpd) :: rhs
2328 REAL(mps) :: rloop
2329 REAL(mpd) :: sgm
2330 REAL(mps) :: used
2331 REAL(mps) :: usei
2332 REAL(mpd) :: weight
2333 INTEGER(mpi) :: i
2334 INTEGER(mpi) :: ia
2335 INTEGER(mpi) :: ib
2336 INTEGER(mpi) :: ibuf
2337 INTEGER(mpi) :: inder
2338 INTEGER(mpi) :: ioffb
2339 INTEGER(mpi) :: ipr
2340 INTEGER(mpi) :: isfrst
2341 INTEGER(mpi) :: islast
2342 INTEGER(mpi) :: itgbi
2343 INTEGER(mpi) :: itgbij
2344 INTEGER(mpi) :: itgbik
2345 INTEGER(mpi) :: ivgb
2346 INTEGER(mpi) :: ivgbij
2347 INTEGER(mpi) :: ivgbik
2348 INTEGER(mpi) :: j
2349 INTEGER(mpi) :: k
2350 INTEGER(mpi) :: lastit
2351 INTEGER(mpi) :: lun
2352 INTEGER(mpi) :: ncrit
2353 INTEGER(mpi) :: ndfs
2354 INTEGER(mpi) :: ngras
2355 INTEGER(mpi) :: nparl
2356 INTEGER(mpi) :: nr
2357 INTEGER(mpi) :: nrej
2358 INTEGER(mpi) :: inone
2359 INTEGER(mpi) :: ilow
2360 INTEGER(mpi) :: nlow
2361 INTEGER(mpi) :: nzero
2362 LOGICAL :: btest
2363
2364 REAL(mpd):: adder
2365 REAL(mpd)::funref
2366 REAL(mpd)::dchi2s
2367 REAL(mpd)::sndf
2368 INTEGER(mpl):: ii
2369 SAVE
2370 ! ...
2371 isfrst(ibuf)=readbufferpointer(ibuf)+1
2373 inder(i)=readbufferdatai(i)
2374 glder(i)=readbufferdatad(i)
2375 ! ----- book and reset ---------------------------------------------
2377 lastit=-1
2378 iitera=0
2379 END IF
2380
2382 ndfsum=0
2383 sumndf=0.0_mpd
2384 funref=0.0_mpd
2385
2390 END IF
2397 END IF
2398
2399
2401 'Normalized residuals of single (global) measurement')
2403 'Normalized residuals of single (local) measurement')
2405 'Pulls of single (global) measurement')
2407 'Pulls of single (local) measurement')
2411
2412 ! WRITE(*,*) 'LOOPN ', NLOOPN, ' executing ICALCM=', ICALCM
2413
2414 ! reset
2415
2417 globalcounter=0
2419 globalmatd=0.0_mpd
2420 globalmatf=0.
2422 END IF
2423
2425
2426 newite=.false.
2428 newite=.true.
2429 funref=fvalue
2432 ! CALL MEND
2439 END IF
2440 END IF
2441 ! WRITE(*,111) ! header line
2442 END IF
2443
2444 DO i=0,3
2446 END DO
2447 DO k=3,6
2449 writebufferheader(-k)=0
2450 END DO
2451 ! statistics per binary file
2453 jfd(i)=0
2454 cfd(i)=0.0
2455 dfd(i)=0
2456 END DO
2457
2459
2460 ! ----- read next data ----------------------------------------------
2461 DO
2464 ndfs =0
2465 sndf =0.0_mpd
2466 dchi2s=0.0_mpd
2471 IF (nr <= 0) EXIT ! next block of events ?
2472 END DO
2473 ! sum up RHS (over threads) once (reduction in LOOPBF: summation for each block)
2474 ioffb=0
2476 ioffb=ioffb+lenglobalvec
2480 END DO
2481 END DO
2482
2484 ! PRINT *, ' cache/w ',(writeBufferHeader(-K),K=3,6),(writeBufferHeader(K),K=3,6)
2485 nparl=writebufferheader(3)
2486 ncrit=writebufferheader(4)
2489 peakd=real(writebufferheader(-6),mps)*0.1
2490 peaki=real(writebufferheader(6),mps)*0.1
2491 WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
2492111 FORMAT(' Write cache usage (#flush,#overrun,<levels>,', &
2493 'peak(levels))'/2i7,',',4(f6.1,'%'))
2494 ! fill second half (j>i) of global matric for extended storage
2496 END IF
2497
2498 ! check entries/counters
2499 nlow=0
2500 ilow=1
2501 nzero=0
2505 nlow=nlow+1
2507 END IF
2508 END DO
2509 IF(nlow > 0) THEN
2511 itgbi=globalparvartototal(ilow)
2512 print *
2513 print *, " ... warning ..."
2514 print *, " global parameters with too few (< MREQENA) accepted entries: ", nlow
2515 print *, " minimum entries: ", globalcounter(ilow), " for label ", globalparlabelindex(1,itgbi)
2516 print *
2517 END IF
2525 isubit=1
2526 WRITE(*,*) ' --> enforcing SUBITO mode'
2527 WRITE(lun,*) ' --> enforcing SUBITO mode'
2528 END IF
2529 END IF
2530
2531 ! ----- after end-of-data add contributions from pre-sigma ---------
2532
2534 ! plot diagonal elements
2535 elmt=0.0
2537 ii=0
2539 ii=i
2540 ii=(ii*ii+ii)/2
2541 END IF
2544 IF(ii /= 0) THEN
2545 elmt=real(globalmatd(ii),mps)
2546 IF(elmt > 0.0) CALL hmpent(23,1.0/sqrt(elmt))
2547 END IF
2548 END DO
2549 END IF
2550
2551
2552
2553 ! add pre-sigma contributions to matrix diagonal
2554
2555 ! WRITE(*,*) 'Adding to diagonal ICALCM IND6',ICALCM,IND6
2556
2559 ! WRITE(*,*) 'Index ',IVGB,IVGB,QM(IND6+IVGB)
2562 END IF
2563 END DO
2564 END IF
2565
2566 CALL hmpwrt(23)
2567 CALL hmpwrt(24)
2568 CALL hmpwrt(25)
2569 CALL hmpwrt(26)
2570
2571
2572 ! add regularization term to F and to rhs --------------------------
2573
2574 ! WRITE(*,*) 'NREGUL ',NREGUL,NLOOPN
2575
2582 END DO
2583 END IF
2584
2585
2586 ! ----- add contributions from "measurement" -----------------------
2587
2588
2589 i=1
2593 i=i+2
2594 weight=0.0
2595 IF(sgm > 0.0) weight=1.0/sgm**2
2596
2597 dsum=-rhs
2598
2599 ! loop over label/factor pairs
2600 ia=i
2601 DO
2602 i=i+1
2605 END DO
2606 ib=i-1
2607
2608 DO j=ia,ib
2609 factrj=listmeasurements(j)%value
2611 IF(itgbij /= 0) THEN
2613 END IF
2614 ! add to vector
2615 ivgbij=0
2617 IF(ivgbij > 0) THEN
2620 END IF
2621
2623 DO k=ia,j
2624 factrk=listmeasurements(k)%value
2626 ! add to matrix
2627 ivgbik=0
2629 IF(ivgbij > 0.AND.ivgbik > 0) THEN !
2631 END IF
2632 END DO
2633 END IF
2634 END DO
2635
2636 adder=weight*dsum**2
2638
2639 END DO
2640
2641 ! ----- printout ---------------------------------------------------
2642
2643
2649 ngras=nint(angras,mpi)
2650 ratae=0.0 !!!
2653 ratae=max(-99.9,ratae)
2654 END IF
2655
2656 ! rejects ...
2657
2660 IF(nrej /= 0) THEN
2661 WRITE(*,*) ' '
2662 WRITE(*,*) 'Data rejected in initial loop:'
2663 WRITE(*,*) ' ', &
2666 END IF
2667 END IF
2668 ! IF(NREJEC(1)+NREJEC(2)+NREJEC(3).NE.0) THEN
2669 ! WRITE(LUNLOG,*) 'Data rejected in initial loop:',NREJEC(1),
2670 ! + ' (Ndf=0) ',NREJEC(2),' (huge) ',NREJEC(3),' (large)'
2671 ! END IF
2672
2673
2678 END IF
2681 END IF
2682 END IF
2683
2686 ! CALL HMPRNT(3) ! scaled residual of single measurement
2687 ! CALL HMPRNT(12) ! scaled residual of single measurement
2688 ! CALL HMPRNT(4) ! chi^2/Ndf
2689 END IF
2690 CALL hmpwrt(3)
2691 CALL hmpwrt(12)
2692 CALL hmpwrt(4)
2693 CALL gmpwrt(4) ! location, dispersion (res.) as a function of record nr
2695 CALL hmpwrt(13)
2696 CALL hmpwrt(14)
2697 CALL gmpwrt(5) ! location, dispersion (pull) as a function of record nr
2698 END IF
2699 END IF
2700 ! IF(NLOOPN.EQ.2.AND.NHISTP.NE.0) CALL HMPRNT(6)
2703 ! IF(NHISTP.NE.0) CALL HMPRNT(5) ! number of degrees of freedom
2704 ! IF(NHISTP.NE.0) CALL HMPRNT(11) ! Nlocal
2705 CALL hmpwrt(5)
2706 CALL hmpwrt(11)
2707 END IF
2708
2709 ! local fit: band matrix structure !?
2711 DO lun=6,8,2
2712 WRITE(lun,*) ' '
2713 WRITE(lun,*) ' === local fits have bordered band matrix structure ==='
2718 END DO
2719 END IF
2720
2721 lastit=iterat
2722
2723 ! monitoring of residuals
2725
2726101 FORMAT(1x,a8,' =',i10,' = ',a)
2727! 101 FORMAT(' LOOPN',I6,' Function value',F22.8,10X,I6,' records')
2728! 102 FORMAT(' incl. constraint penalty',F22.8)
2729! 103 FORMAT(I13,3X,A,G12.4)
real(mps) prange
range (-PRANGE..PRANGE) for histograms of pulls, norm. residuals
Definition mpmod.f90:97
subroutine loopbf(nrej, ndfs, sndf, dchi2s, numfil, naccf, chi2f, ndff)
Loop over records in read buffer (block), fits and sums.
Definition pede.f90:3064
◆ mcsolv()
| subroutine mcsolv | ( | integer(mpi), intent(in) | n, |
| real(mpd), dimension(n), intent(in) | x, | ||
| real(mpd), dimension(n), intent(out) | y | ||
| ) |
Solution for zero band width preconditioner.
Used by MINRES.
- Parameters
-
[in] n size of vectors [in] x rhs vector [out] y result vector
Definition at line 6579 of file pede.f90.
6581
6582 IMPLICIT NONE
6583 INTEGER(mpi),INTENT(IN) :: n
6584 REAL(mpd), INTENT(IN) :: x(n)
6585 REAL(mpd), INTENT(OUT) :: y(n)
6586 SAVE
6587 ! ...
6588 CALL presol(nprecond(1),nprecond(2),matprecond,matprecond(1+nvgb),matprecond(1+nvgb+ncgb*nvgb),y,x)
subroutine presol(p, n, cu, a, s, x, y)
Constrained preconditioner, solution.
Definition mpnum.f90:2246
◆ mdiags()
| subroutine mdiags |
Solution by diagonalization.
Definition at line 6206 of file pede.f90.
6208
6209 IMPLICIT NONE
6210 REAL(mps) :: evalue
6211 INTEGER(mpi) :: i
6212 INTEGER(mpi) :: iast
6213 INTEGER(mpi) :: idia
6214 INTEGER(mpi) :: imin
6215 INTEGER(mpl) :: ioff1
6216 INTEGER(mpi) :: j
6217 INTEGER(mpi) :: lun
6218 INTEGER(mpi) :: nmax
6219 INTEGER(mpi) :: nmin
6220 INTEGER(mpi) :: ntop
6221 !
6222 INTEGER(mpl) :: ii
6224
6225 SAVE
6226 ! ...
6227
6230
6231 ! save diagonal (for global correlation)
6234 ii=i
6236 END DO
6237 ENDIF
6238
6239 !use elimination for constraints ?
6242 ! solve L^t*y=d by backward substitution
6244 ! transform, reduce rhs
6246 ! correction from eliminated part
6251 ioff1=ioff1+nfgb+j
6252 END DO
6253 END DO
6254 END IF
6255
6257 ! eigenvalues eigenvectors symm_input
6258 workspaceeigenvalues=0.0_mpd
6261
6262 ! histogram of positive eigenvalues
6263
6265 imin=1
6268 imin=i ! index of smallest pos. eigenvalue
6269 EXIT
6270 END IF
6271 END DO
6273 ntop=nmin+6
6274 DO WHILE(ntop < nmax)
6275 ntop=ntop+3
6276 END DO
6277
6281 evalue=log10(real(workspaceeigenvalues(idia),mps))
6282 CALL hmpent(7,evalue)
6283 END IF
6284 END DO
6286 CALL hmpwrt(7)
6287
6288 iast=max(1,imin-60)
6291 evalue=real(workspaceeigenvalues(i),mps)
6292 CALL gmpxy(3,real(i,mps),evalue)
6293 END DO
6295 CALL gmpwrt(3)
6296
6298 workspacediagonalization(i)=0.0_mpd
6302 END IF
6303 END DO
6304 WRITE(lun,*) ' '
6305 WRITE(lun,*) 'The first (largest) eigenvalues ...'
6307 WRITE(lun,*) ' '
6310 WRITE(lun,*) ' '
6314 WRITE(lun,*) ' '
6315 ENDIF
6317 WRITE(lun,*) '(for Eigenvalue < 0.001 the value 0.0 is shown)'
6320 'printed for negative eigenvalues'
6321 CALL devsig(nfgb,workspaceeigenvalues,workspaceeigenvectors,globalvector,workspacediagonalization)
6322 WRITE(lun,*) ' '
6324 'compatible with N(0,1)'
6326
6327
6328101 FORMAT(10f7.1)
6329102 FORMAT(5e14.6)
6330
6331 END IF
6332
6333 ! solution ---------------------------------------------------------
6335 ! eigenvalues eigenvectors
6336 CALL devsol(nfgb,workspaceeigenvalues,workspaceeigenvectors,workspaced,globalcorrections,workspacediagonalization)
6337
6338 !use elimination for constraints ?
6340 ! extend, transform back solution
6343 END IF
6344
real(mpd), dimension(:), allocatable workspaceeigenvectors
workspace eigen vectors
Definition mpmod.f90:199
real(mpd), dimension(:), allocatable globalcorrections
correction x_i (from A*x_i=b_i in iteration i)
Definition mpmod.f90:176
real(mpd), dimension(:), allocatable workspacediagonalization
workspace diag.
Definition mpmod.f90:197
real(mpd), dimension(:), allocatable workspaceeigenvalues
workspace eigen values
Definition mpmod.f90:198
real(mpd), dimension(:), allocatable workspacediag
diagonal of global matrix (for global corr.)
Definition mpmod.f90:195
◆ mend()
| subroutine mend |
◆ mhalf2()
| subroutine mhalf2 |
Fill 2nd half of matrix for extended storage.
Definition at line 4535 of file pede.f90.
4537
4538 IMPLICIT NONE
4539 INTEGER(mpi) :: i
4540 INTEGER(mpi) :: ichunk
4541 INTEGER(mpi) :: iencdb
4542 INTEGER(mpi) :: iencdm
4543 INTEGER(mpi) :: ir
4544 INTEGER(mpi) :: ispc
4545 INTEGER(mpi) :: j
4546 INTEGER(mpi) :: jtem
4547 INTEGER(mpi) :: jtemc
4548 INTEGER(mpi) :: jtemn
4549 INTEGER(mpi) :: nd
4550
4551 INTEGER(mpl) :: ij
4552 INTEGER(mpl) :: ijadd
4553 INTEGER(mpl) :: k
4554 INTEGER(mpl) :: kk
4555 INTEGER(mpl) :: kl
4556 INTEGER(mpl) :: ku
4557 INTEGER(mpl) :: indid
4558 INTEGER(mpl) :: ll
4559 INTEGER(mpl) :: k8
4560 ! ...
4561
4564
4566 iencdm=ishft(1,iencdb)-1
4568 ! parallelize row loop
4569 ! slot of 1024 'I' for next idle thread
4570 !$OMP PARALLEL DO &
4571 !$OMP PRIVATE(I,IR,K,KK,LL,KL,KU,K8,INDID,IJ,J,JTEMC,JTEM,JTEMN) &
4572 !$OMP SCHEDULE(DYNAMIC,ichunk)
4573 DO i=1,nd
4574 ir=i+(ispc-1)*(nd+1)
4578
4579 ku=sparsematrixoffsets(1,ir+1)-kk
4580 !IF (sparseMatrixColumns(kk) == 0) THEN ! compression ?
4581 ku=(ku*8)/9-1 ! number of regions (-1)
4582 indid=kk+ku/8+1 ! index of first region (after group offsets)
4583 k8=kl/8 ! region group (-1)
4585 DO k=k8*8,kl-1
4587 END DO
4588 DO k=kl,ku
4590 jtem =ishft(jtemc,-iencdb) ! first column of region
4591 jtemn=jtem+iand(jtemc,iencdm) ! first column after region
4592 DO j=jtem,jtemn-1
4593 ij=ijadd(i,j)
4594 IF (ispc==1) THEN
4596 ELSE
4598 END IF
4599 ll=ll+1
4600 END DO
4601 END DO
4602 END DO
4603 !$OMP END PARALLEL DO
4604 END DO
4605
◆ minver()
| subroutine minver |
Solution by matrix inversion.
Parallelized (SQMINL), solve A*x=b.
Definition at line 6129 of file pede.f90.
6131
6132 IMPLICIT NONE
6133 INTEGER(mpi) :: i
6134 INTEGER(mpl) :: ioff1
6135 INTEGER(mpi) :: j
6136 INTEGER(mpi) :: lun
6137 INTEGER(mpi) :: nrank
6138 INTEGER(mpl) :: ii
6140
6141 SAVE
6142 ! ...
6145
6146 ! save diagonal (for global correlation)
6149 ii=i
6151 END DO
6152 ENDIF
6153
6154 ! WRITE(*,*) 'MINVER ICALCM=',ICALCM
6155 !use elimination for constraints ?
6158 ! solve L^t*y=d by backward substitution
6160 ! transform, reduce rhs
6162 ! correction from eliminated part
6167 ioff1=ioff1+nfgb+j
6168 END DO
6169 END DO
6170 END IF
6171
6173 ! invert and solve
6182 isubit=1
6183 WRITE(*,*) ' --> enforcing SUBITO mode'
6184 WRITE(lun,*) ' --> enforcing SUBITO mode'
6185 END IF
6187 WRITE(lun,*) 'No rank defect of the symmetric matrix'
6188 END IF
6190
6191 ELSE ! multiply gradient by inverse matrix
6194 END IF
6195
6196 !use elimination for constraints ?
6198 ! extend, transform back solution
6201 END IF
6202
subroutine sqminl(v, b, n, nrank, diag, next)
Matrix inversion for LARGE matrices.
Definition mpnum.f90:225
◆ mminrs()
| subroutine mminrs |
Solution with MINRES.
Solve A*x=b by minimizing |A*x-b| iteratively. Parallelized (AVPROD).
Use preconditioner with zero (precon) or finite (equdec) band width.
Definition at line 6385 of file pede.f90.
6388
6389 IMPLICIT NONE
6390 INTEGER(mpi) :: istop
6391 INTEGER(mpi) :: itn
6392 INTEGER(mpi) :: itnlim
6393 INTEGER(mpi) :: lun
6394 INTEGER(mpi) :: nout
6395 INTEGER(mpi) :: nrkd
6396 INTEGER(mpi) :: nrkd2
6397
6398 REAL(mpd) :: shift
6399 REAL(mpd) :: rtol
6400 REAL(mpd) :: anorm
6401 REAL(mpd) :: acond
6402 REAL(mpd) :: arnorm
6403 REAL(mpd) :: rnorm
6404 REAL(mpd) :: ynorm
6405 LOGICAL :: checka
6407 SAVE
6408 ! ...
6411
6412 nout=lun
6413 itnlim=2000 ! iteration limit
6414 shift =0.0_mpd ! not used
6416 checka=.false.
6417
6419 !use elimination for constraints ?
6421 ! solve L^t*y=d by backward substitution
6423 ! input to AVPRD0
6427 ! calclulate vecBav=globalMat*vecXav
6429 ! correction from eliminated part
6431 ! transform, reduce rhs
6433 END IF
6434
6440 END IF
6442 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
6448 WRITE(lun,*) 'MMINRS: EQUDEC ended ', nrkd, nrkd2
6449 END IF
6451 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
6452 ELSE
6454 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
6455 END IF
6456
6457 !use elimination for constraints ?
6459 ! extend, transform back solution
6462 END IF
6463
6464 iitera=itn
6465 istopa=istop
6467
6469
subroutine precon(p, n, c, cu, a, s, nrkd)
Constrained preconditioner, decomposition.
Definition mpnum.f90:2172
subroutine equdec(n, m, ls, c, india, nrkd, nrkd2)
Decomposition of equilibrium systems.
Definition mpnum.f90:2035
subroutine qlpssq(aprod, b, m, t)
Partial similarity transformation by Q(t).
Definition mpqldec.f90:498
MINRES solves symmetric systems Ax = b or min ||Ax - b||_2, where the matrix A may be indefinite and/...
Definition minresModule.f90:39
subroutine, public minres(n, aprod, msolve, b, shift, checka, precon, x, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
Solution of linear equation system.
Definition minresModule.f90:297
integer(mpi) lprecm
additional flag for preconditioner (band) matrix (>0: preserve rank by skyline matrix)
Definition mpmod.f90:44
integer(mpi), dimension(:), allocatable indprecond
preconditioner pointer array
Definition mpmod.f90:192
◆ mminrsqlp()
| subroutine mminrsqlp |
Solution with MINRES-QLP.
Solve A*x=b by minimizing |A*x-b| iteratively. Parallelized (AVPROD).
Use preconditioner with zero (precon) or finite (equdec) band width.
Definition at line 6478 of file pede.f90.
6481
6482 IMPLICIT NONE
6483 INTEGER(mpi) :: istop
6484 INTEGER(mpi) :: itn
6485 INTEGER(mpi) :: itnlim
6486 INTEGER(mpi) :: lun
6487 INTEGER(mpi) :: nout
6488 INTEGER(mpi) :: nrkd
6489 INTEGER(mpi) :: nrkd2
6490
6491 REAL(mpd) :: rtol
6492 REAL(mpd) :: mxxnrm
6493 REAL(mpd) :: trcond
6494
6496 SAVE
6497 ! ...
6500
6501 nout=lun
6502 itnlim=2000 ! iteration limit
6504 mxxnrm = real(nagb,mpd)/sqrt(epsilon(mxxnrm))
6506 trcond = 1.0_mpd/epsilon(trcond) ! only QR
6508 trcond = 1.0_mpd ! only QLP
6509 ELSE
6511 END IF
6512
6514 !use elimination for constraints ?
6516 ! solve L^t*y=d by backward substitution
6518 ! input to AVPRD0
6522 ! calclulate vecBav=globalMat*vecXav
6524 ! correction from eliminated part
6526 ! transform, reduce rhs
6528 END IF
6529
6535 END IF
6537 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
6538 x=globalcorrections, istop=istop, itn=itn)
6544 WRITE(lun,*) 'MMINRS: EQUDEC ended ', nrkd, nrkd2
6545 END IF
6546
6548 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
6549 x=globalcorrections, istop=istop, itn=itn)
6550 ELSE
6552 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
6553 x=globalcorrections, istop=istop, itn=itn)
6554 END IF
6555
6556 !use elimination for constraints ?
6558 ! extend, transform back solution
6561 END IF
6562
6563 iitera=itn
6564 istopa=istop
6566
6568
MINRESQLP solves symmetric systems Ax = b or min ||Ax - b||_2, where the matrix A may be indefinite a...
Definition minresqlpModule.f90:42
subroutine, public minresqlp(n, aprod, b, shift, msolve, disable, nout, itnlim, rtol, maxxnorm, trancond, acondlim, x, istop, itn, rnorm, arnorm, xnorm, anorm, acond)
Solution of linear equation system or least squares problem.
Definition minresqlpModule.f90:414
integer(mpi) mrmode
MINRES-QLP mode (0: QR+QLP, 1: only QR, 2: only QLP factorization)
Definition mpmod.f90:62
◆ monres()
| subroutine monres |
Monitor input residuals.
Read all data files again to monitor input residuals
Definition at line 5940 of file pede.f90.
5943
5944 IMPLICIT NONE
5945 INTEGER(mpi) :: i
5946 INTEGER(mpi) :: ij
5947 INTEGER(mpi) :: imed
5948 INTEGER(mpi) :: j
5949 INTEGER(mpi) :: k
5950 INTEGER(mpi) :: nent
5951 INTEGER(mpi), DIMENSION(measBins) :: isuml ! location
5952 INTEGER(mpi), DIMENSION(measBins) :: isums ! scale
5953 REAL(mps) :: amed
5954 REAL(mps) :: amad
5955
5956 INTEGER(mpl) :: ioff
5957 LOGICAL :: lfirst
5958 SAVE
5959 DATA lfirst /.true./
5960
5961 ! combine data from threads
5962 ioff=0
5967 END DO
5968 END DO
5969
5970 IF (lfirst) THEN
5972 WRITE(lunmon,'(A)') '*** Normalized residuals grouped by first global label (per local fit cycle) ***'
5973 ELSE
5975 ENDIF
5977 lfirst=.false.
5978 END IF
5979
5980 ! analyze histograms
5981 ioff=0
5984 isuml=0
5985 ! sum up content
5986 isuml(1)=meashists(ioff+1)
5988 isuml(j)=isuml(j-1)+meashists(ioff+j)
5989 END DO
5990 nent=isuml(measbins)
5991 ! get median (for location)
5993 IF (2*isuml(j) > nent) EXIT
5994 END DO
5995 imed=j
5996 amed=real(j,mps)
5997 IF (isuml(j) > isuml(j-1)) amed=amed+real(nent-2*isuml(j-1),mps)/real(2*isuml(j)-2*isuml(j-1),mps)
5999 ! sum up differences
6000 isums = 0
6002 k=j-imed+1
6003 isums(k)=isums(k)+meashists(ioff+j)
6004 END DO
6005 DO j=imed-1,1,-1
6006 k=imed-j
6007 isums(k)=isums(k)+meashists(ioff+j)
6008 END DO
6010 isums(j)=isums(j)+isums(j-1)
6011 END DO
6012 ! get median (for scale)
6014 IF (2*isums(j) > nent) EXIT
6015 END DO
6016 amad=real(j-1,mps)
6017 IF (isums(j) > isums(j-1)) amad=amad+real(nent-2*isums(j-1),mps)/real(2*isums(j)-2*isums(j-1),mps)
6018 amad=real(measbinsize,mps)*amad
6019 ij=globalparlabelindex(1,i)
6021 !
6022 ioff=ioff+measbins
6023 END IF
6024 END DO
6025
6026110 FORMAT(i5,2i10,3g14.5)
◆ mptwo()
◆ mstart()
| subroutine mstart | ( | character (len=*), intent(in) | text | ) |
Start of 'module' printout.
Definition at line 8749 of file pede.f90.
8752
8753 IMPLICIT NONE
8754 INTEGER(mpi) :: i
8755 INTEGER(mpi) :: ka
8756 INTEGER(mpi) :: kb
8757 INTEGER(mpi) :: l
8758 CHARACTER (LEN=*), INTENT(IN) :: text
8759 CHARACTER (LEN=16) :: textc
8760 SAVE
8761 ! ...
8762 DO i=1,74
8764 END DO
8765 l=len(text)
8766 ka=(74-l)/2
8767 kb=ka+l-1
8768 textl(ka:kb)=text(1:l)
8769 WRITE(*,*) ' '
8771 WRITE(*,*) ' '
8772 textc=text(1:l)//'-end'
8773
8774 DO i=1,74
8776 END DO
8777 l=l+4
8778 ka=(74-l)/2
8779 kb=ka+l-1
8780 textl(ka:kb)=textc(1:l)
8781 RETURN
◆ mupdat()
| subroutine mupdat | ( | integer(mpi), intent(in) | i, |
| integer(mpi), intent(in) | j, | ||
| real(mpd), intent(in) | add | ||
| ) |
Update element of global matrix.
Add value ADD to matrix element (I,J).
- Parameters
-
[in] i first index [in] j second index [in] add summand
Definition at line 2977 of file pede.f90.
2979
2980 IMPLICIT NONE
2981
2982 INTEGER(mpi), INTENT(IN) :: i
2983 INTEGER(mpi), INTENT(IN) :: j
2984 REAL(mpd), INTENT(IN) :: add
2985
2986 INTEGER(mpl):: ijadd
2987 INTEGER(mpl):: ia
2988 INTEGER(mpl):: ja
2989 INTEGER(mpl):: ij
2990 ! ...
2991 IF(i <= 0.OR.j <= 0) RETURN
2992 ia=max(i,j) ! larger
2993 ja=min(i,j) ! smaller
2994 ij=0
2996 ij=ja+(ia*ia-ia)/2 ! ISYM index
3000 IF (ij == 0) RETURN ! pair is suppressed
3001 IF (ij > 0) THEN
3003 ELSE
3005 END IF
3006 END IF
3010 ij=indprecond(ia)-ia+ja
3015 stop 'mupdat: bad index'
3016 END IF
3017 ELSE ! Lagrange multiplier
3020 END IF
3024 ELSE ! Lagrange multiplier
3027 END IF
3028 END IF
3029 END IF
◆ mvopen()
| subroutine mvopen | ( | integer(mpi), intent(in) | lun, |
| character (len=*), intent(in) | fname | ||
| ) |
Open file.
Evtl. move existing file to ~ version.
- Parameters
-
[in] lun unit number [in] fname file name
Definition at line 8802 of file pede.f90.
8804
8805 IMPLICIT NONE
8806 INTEGER(mpi) :: l
8807 INTEGER(mpi), INTENT(IN) :: lun
8808 CHARACTER (LEN=*), INTENT(IN) :: fname
8809 CHARACTER (LEN=33) :: nafile
8810 CHARACTER (LEN=33) :: nbfile
8811 LOGICAL :: ex
8812 SAVE
8813 ! ...
8814 l=len(fname)
8815 IF(l > 32) THEN
8817 stop 'File name too long '
8818 END IF
8819 nafile=fname
8820 nafile(l+1:l+1)='~'
8821
8822 INQUIRE(file=nafile(1:l),exist=ex)
8823 IF(ex) THEN
8824 INQUIRE(file=nafile(1:l+1),exist=ex)
8825 IF(ex) THEN
8826 CALL system('rm '//nafile)
8827 END IF
8828 nbfile=nafile
8829 nafile(l+1:l+1)=' '
8830 CALL system('mv '//nafile//nbfile)
8831 END IF
8832 OPEN(unit=lun,file=fname)
◆ mvsolv()
| subroutine mvsolv | ( | integer(mpi), intent(in) | n, |
| real(mpd), dimension(n), intent(in) | x, | ||
| real(mpd), dimension(n), intent(out) | y | ||
| ) |
Solution for finite band width preconditioner.
Used by MINRES.
- Parameters
-
[in] n size of vectors [in] x rhs vector [out] y result vector
Definition at line 6599 of file pede.f90.
6601
6602 IMPLICIT NONE
6603
6604 INTEGER(mpi), INTENT(IN) :: n
6605 REAL(mpd), INTENT(IN) :: x(n)
6606 REAL(mpd), INTENT(OUT) :: y(n)
6607
6608 SAVE
6609 ! ...
6610 y=x ! copy to output vector
6611
subroutine equslv(n, m, c, india, x)
Solution of equilibrium systems (after decomposition).
Definition mpnum.f90:2101
◆ nufile()
| integer(mpi) function nufile | ( | character (len=*), intent(inout) | fname | ) |
Inquire on file.
Result = 1 for existing binary file, =2 for existing text file, else =0, < 0 open error.
Text file names are recognized by the filename extension, which should contain 'xt' or 'tx'.
- Parameters
-
[in,out] fname file name, optionaly strip prefix.
Definition at line 7965 of file pede.f90.
7967
7968 IMPLICIT NONE
7969 INTEGER(mpi) :: ios
7970 INTEGER(mpi) :: l1
7971 INTEGER(mpi) :: ll
7972 INTEGER(mpi) :: nm
7973 INTEGER(mpi) :: npat
7974 INTEGER(mpi) :: ntext
7975 INTEGER(mpi) :: nuprae
7976 INTEGER(mpi) :: matint
7977
7978 CHARACTER (LEN=*), INTENT(INOUT) :: fname
7979 LOGICAL :: ex
7980 SAVE
7981 ! ...
7982 nufile=0
7983 IF(fname(1:5) == 'rfio:') nuprae=1
7984 IF(fname(1:5) == 'dcap:') nuprae=2
7985 IF(fname(1:5) == 'root:') nuprae=3
7986 IF(nuprae == 0) THEN
7987 INQUIRE(file=fname,iostat=ios,exist=ex)
7989 IF(ios /= 0) RETURN
7990 ELSE IF(nuprae == 1) THEN ! rfio:
7991 ll=len(fname)
7992 fname=fname(6:ll)
7993 ex=.true.
7994 nufile=1
7995 RETURN
7996 ELSE
7997 ex=.true. ! assume file existence
7998 END IF
7999 IF(ex) THEN
8001 ll=len(fname)
8002 l1=max(1,ll-3)
8005 IF(nm < 2) THEN
8008 END IF
8009 END IF
◆ pechk()
| subroutine pechk | ( | integer(mpi), intent(in) | ibuf, |
| integer(mpi), intent(out) | nerr | ||
| ) |
Check Millepede record.
Check integer structure of labels and markers (zeros). Check floats for NaNs.
- Parameters
-
[in] ibuf buffer number [out] nerr error flags
Definition at line 2161 of file pede.f90.
2163
2164 IMPLICIT NONE
2165 REAL(mpr8) :: glder
2166 INTEGER(mpi) :: i
2167 INTEGER(mpi) :: is
2168 INTEGER(mpi) :: ist
2169 INTEGER(mpi) :: inder
2170 INTEGER(mpi) :: ioff
2171 INTEGER(mpi) :: isfrst
2172 INTEGER(mpi) :: islast
2173 INTEGER(mpi) :: ja
2174 INTEGER(mpi) :: jb
2175 INTEGER(mpi) :: jsp
2176 INTEGER(mpi) :: nan
2177 INTEGER(mpi) :: nst
2178
2179 INTEGER(mpi), INTENT(IN) :: ibuf
2180 INTEGER(mpi), INTENT(OUT) :: nerr
2181 SAVE
2182 ! ...
2183 inder(i)=readbufferdatai(i)
2184 glder(i)=readbufferdatad(i)
2185 isfrst(ibuf)=readbufferpointer(ibuf)+1
2187
2188 ist=isfrst(ibuf)
2189 nst=islast(ibuf)
2190 nerr=0
2191 is=ist
2192 jsp=0
2193 outer: DO WHILE(is < nst)
2194 ja=0
2195 jb=0
2196 inner1: DO
2197 is=is+1
2198 IF(is > nst) EXIT outer
2199 IF(inder(is) == 0) EXIT inner1 ! found 1. marker
2200 END DO inner1
2201 ja=is
2202 inner2: DO
2203 is=is+1
2204 IF(is > nst) EXIT outer
2205 IF(inder(is) == 0) EXIT inner2 ! found 2. marker
2206 END DO inner2
2207 jb=is
2208 IF(ja+1 == jb.AND.glder(jb) < 0.0_mpr8) THEN
2209 ! special data
2210 jsp=jb ! pointer to special data
2211 is=is+nint(-glder(jb),mpi) ! skip NSP words
2212 cycle outer
2213 END IF
2214 DO WHILE(inder(is+1) /= 0.AND.is < nst)
2215 is=is+1
2216 END DO
2217 END DO outer
2218 IF(is > nst) THEN
2219 ioff = readbufferpointer(ibuf)
2221100 FORMAT(' PEREAD: record ', i8,' in file ',i6, ' is broken !!!')
2222 nerr=nerr+1
2223 ENDIF
2224 nan=0
2225 DO i=ist, nst
2227 END DO
2228 IF(nan > 0) THEN
2229 ioff = readbufferpointer(ibuf)
2231101 FORMAT(' PEREAD: record ', i8,' in file ',i6, ' contains ', i6, ' NaNs !!!')
2232 nerr= nerr+2
2233 ENDIF
2234
◆ peend()
| subroutine peend | ( | integer(mpi), intent(in) | icode, |
| character (len=*), intent(in) | cmessage | ||
| ) |
Print exit code.
Print exit code and message.
- Parameters
-
[in] icode exit code [in] cmessage exit massage
◆ peprep()
| subroutine peprep | ( | integer(mpi), intent(in) | mode | ) |
Prepare records.
For global parameters replace label by index (INONE).
- Parameters
-
[in] mode <=0: build index table (INONE) for global variables;
>0: use index table, can be parallelized, optional scale errors
Definition at line 2066 of file pede.f90.
2068
2069 IMPLICIT NONE
2070
2071 INTEGER(mpi), INTENT(IN) :: mode
2072
2073 INTEGER(mpi) :: ibuf
2074 INTEGER(mpi) :: ichunk
2075 INTEGER(mpi) :: iproc
2076 INTEGER(mpi) :: isfrst
2077 INTEGER(mpi) :: islast
2078 INTEGER(mpi) :: ist
2079 INTEGER(mpi) :: j
2080 INTEGER(mpi) :: ja
2081 INTEGER(mpi) :: jb
2082 INTEGER(mpi) :: jsp
2083 INTEGER(mpi) :: nst
2084 INTEGER(mpi), PARAMETER :: maxbad = 100 ! max number of bad records with print out
2085 INTEGER(mpi) :: nbad
2086 INTEGER(mpi) :: nerr
2087 INTEGER(mpi) :: inone
2088 !$ INTEGER(mpi) :: OMP_GET_THREAD_NUM
2089
2090
2091 isfrst(ibuf)=readbufferpointer(ibuf)+1
2093
2094 IF (mode > 0) THEN
2096 ! parallelize record loop
2097 !$OMP PARALLEL DO &
2098 !$OMP DEFAULT(PRIVATE) &
2099 !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI,readBufferDataD,ICHUNK,iscerr,dscerr) &
2100 !$OMP SCHEDULE(DYNAMIC,ICHUNK)
2102 iproc=0
2103 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
2104 ist=isfrst(ibuf)
2105 nst=islast(ibuf)
2106 DO ! loop over measurements
2108 IF(jb == 0) EXIT
2109 DO j=1,ist-jb
2111 END DO
2112 ! scale error ?
2114 IF (jb < ist) THEN
2116 ELSE
2118 END IF
2119 END IF
2120 END DO
2121 END DO
2122 !$OMP END PARALLEL DO
2123 END IF
2124
2125 !$POMP INST BEGIN(peprep)
2126 IF (mode <= 0) THEN
2127 nbad=0
2130 IF(nerr > 0) THEN
2131 nbad=nbad+1
2132 IF(nbad >= maxbad) EXIT
2133 ELSE
2134 ist=isfrst(ibuf)
2135 nst=islast(ibuf)
2136 DO ! loop over measurements
2138 IF(jb == 0) EXIT
2139 DO j=1,ist-jb
2141 END DO
2142 END DO
2143 END IF
2144 END DO
2145 IF(nbad > 0) THEN
2147 stop 'PEREAD: stopping due to bad records'
2148 END IF
2149 END IF
2150 !$POMP INST END(peprep)
2151
real(mpd), dimension(2) dscerr
scaling factors for errors of 'global' and 'local' measurement
Definition mpmod.f90:110
◆ peread()
| subroutine peread | ( | integer(mpi), intent(out) | more | ) |
Read (block of) records from binary files.
Optionally using several threads (each file read by single thread). Records larger than the read buffer (ndimbuf) are skipped. In case of skipped events in the first loop over all binary files the buffer size is adapted to the maximum record size (and the initial loop (LOOP1) is repeated). C binary files are handled by readc.c and may be gzipped.
- Parameters
-
[out] more more records to come
The records consist of parallel integer and real arrays:
real array integer array
1 0.0 error count (this record)
2 RMEAS, measured value 0 JA
3 local derivative index of local derivative
4 local derivative index of local derivative
5 ...
6 SIGMA, error (>0) 0 JB
global derivative label of global derivative
global derivative label of global derivative IST
RMEAS, measured value 0
local derivative index of local derivative
local derivative index of local derivative
...
SIGMA, error 0
global derivative label of global derivative
global derivative label of global derivative
...
NR global derivative label of global derivative
Definition at line 1737 of file pede.f90.
1739
1740 IMPLICIT NONE
1741 INTEGER(mpi) :: i
1742 INTEGER(mpi) :: iact
1743 INTEGER(mpi) :: ierrc
1744 INTEGER(mpi) :: ierrf
1745 INTEGER(mpi) :: inder
1746 INTEGER(mpi) :: ioffp
1747 INTEGER(mpi) :: ios
1748 INTEGER(mpi) :: ithr
1749 INTEGER(mpi) :: jfile
1750 INTEGER(mpi) :: jrec
1751 INTEGER(mpi) :: k
1752 INTEGER(mpi) :: kfile
1753 INTEGER(mpi) :: l
1754 INTEGER(mpi) :: lun
1755 INTEGER(mpi) :: mpri
1756 INTEGER(mpi) :: n
1757 INTEGER(mpi) :: nact
1758 INTEGER(mpi) :: nbuf
1759 INTEGER(mpi) :: ndata
1760 INTEGER(mpi) :: noff
1761 INTEGER(mpi) :: noffs
1762 INTEGER(mpi) :: npointer
1763 INTEGER(mpi) :: npri
1764 INTEGER(mpi) :: nr
1765 INTEGER(mpi) :: nrc
1766 INTEGER(mpi) :: nrd
1767 INTEGER(mpi) :: nrpr
1768 INTEGER(mpi) :: nthr
1769 INTEGER(mpi) :: ntot
1770 INTEGER(mpi) :: maxRecordSize
1771 INTEGER(mpi) :: maxRecordFile
1772
1773 INTEGER(mpi), INTENT(OUT) :: more
1774
1775 LOGICAL :: lprint
1776 LOGICAL :: floop
1777 LOGICAL :: eof
1778 REAL(mpd) :: ds0
1779 REAL(mpd) :: ds1
1780 REAL(mpd) :: ds2
1781 REAL(mpd) :: dw
1782 !$ INTEGER(mpi) :: OMP_GET_THREAD_NUM
1783 CHARACTER (LEN=7) :: cfile
1784 SAVE
1785
1786 inder(i)=readbufferdatai(i)
1787
1788 DATA lprint/.true./
1789 DATA floop/.true./
1790 DATA npri / 0 /, mpri / 1000 /
1791 ! ...
1793 nrec=0
1794 nrecd=0
1795 ntot=0
1796 sumrecords=0
1797 skippedrecords=0
1798 numblocks=0
1800 maxrecordsinblock=0
1802 nrpr=1
1803 nthr=mthrdr
1804 nact=0 ! active threads (have something still to read)
1805 DO k=1,nthr
1809 readbufferinfo(2,k)=nact
1810 nact=nact+1
1811 END IF
1812 END DO
1813 END IF
1816 more=-1
1817 DO k=1,nthr
1818 iact=readbufferinfo(2,k)
1821 END DO
1823
1824 !$OMP PARALLEL &
1825 !$OMP DEFAULT(PRIVATE) &
1826 !$OMP SHARED(readBufferInfo,readBufferPointer,readBufferDataI,readBufferDataD, &
1827 !$OMP readBufferDataF,nPointer,nData,skippedRecords,ndimbuf,NTHR,NFILF,FLOOP, &
1828 !$OMP IFD,KFD,IFILE,NFILB,WFD,XFD,icheck,keepOpen) &
1829 !$OMP NUM_THREADS(NTHR)
1830
1831 ithr=1
1832 !$ ITHR=OMP_GET_THREAD_NUM()+1 ! thread number
1836 ioffp=iact*npointer
1838
1839 files: DO WHILE (jfile > 0)
1840 kfile=kfd(2,jfile)
1841 ! open again
1844 END IF
1845 records: DO
1846 nbuf=readbufferinfo(4,ithr)+1
1848 nr=ndimbuf
1850 lun=kfile+10
1852 (readbufferdatai(noff+i),i=1,min(n/2,nr))
1853 nr=n/2
1854 ! convert to double
1855 DO i=1,nr
1857 END DO
1858 ! IF (ierrf < 0) REWIND lun ! end-of-file ! CHK use binrwd()
1859 eof=(ierrf /= 0)
1860 ELSE ! C file
1861 lun=kfile-nfilf
1863#ifdef READ_C_FILES
1864 CALL readc(readbufferdatad(noff+1),readbufferdataf(noffs+1),readbufferdatai(noff+1),nr,lun,ierrc)
1865 n=nr+nr
1867#else
1868 ierrc=0
1869#endif
1870 eof=(ierrc <= 0.AND.ierrc /= -4) ! allow buffer overruns -> skip record
1871 IF(eof.AND.ierrc < 0) THEN
1872 WRITE(*,*) 'Read error for binary Cfile', kfile, 'record', jrec+1, ':', ierrc
1873 WRITE(8,*) 'Read error for binary Cfile', kfile, 'record', jrec+1, ':', ierrc
1875 WRITE(cfile,'(I7)') kfile
1877 stop 'PEREAD: stopping due to read errors'
1878 ENDIF
1879 END IF
1880 END IF
1881 IF(eof) EXIT records ! end-of-files or error
1882
1883 jrec=jrec+1
1884 readbufferinfo(3,ithr)=jrec
1885 IF(floop) THEN
1887 IF(ithr == 1) THEN
1888 CALL hmplnt(1,n)
1889 IF(inder(noff+1) /= 0) CALL hmpent(8,real(inder(noff+1),mps))
1890 END IF
1891 END IF
1892
1894 readbufferinfo(4,ithr)=nbuf
1895 readbufferinfo(5,ithr)=noff+nr
1896
1902
1904 ELSE
1905 !$OMP ATOMIC
1907 cycle records
1908 END IF
1909
1910 END DO records
1911
1915 ELSE ! rewind
1917 END IF
1919 print *, 'PEREAD: file ', kfile, 'read the first time, found',jrec,' records'
1920 kfd(1,jfile)=-jrec
1921 END IF
1922 ! take next file
1923 !$OMP CRITICAL
1926 jrec=0
1928 readbufferinfo(3,ithr)=jrec
1929 END IF
1930 !$OMP END CRITICAL
1931 jfile=readbufferinfo(1,ithr)
1932
1933 END DO files
1934 !$OMP END PARALLEL
1935 ! compress pointers
1937 DO k=2,nthr
1938 iact =readbufferinfo(2,k)
1939 ioffp=iact*npointer
1940 nbuf=readbufferinfo(4,k)
1941 DO l=1,nbuf
1943 END DO
1944 nrd=nrd+nbuf
1945 END DO
1946
1947 more=0
1948 DO k=1,nthr
1949 jfile=readbufferinfo(1,k)
1950 IF (jfile > 0) THEN ! no eof yet
1951 readbufferinfo(2,k)=more
1952 more=more+1
1953 ELSE
1954 ! no more files, thread retires
1955 readbufferinfo(1,k)=0
1956 readbufferinfo(2,k)=-1
1957 readbufferinfo(3,k)=0
1959 readbufferinfo(6,k)=0
1960 END IF
1961 END DO
1962 ! record limit ?
1964 nrd=mxrec-ntot
1965 more=-1
1966 DO k=1,nthr
1967 jfile=readbufferinfo(1,k)
1968 IF (jfile > 0) THEN ! rewind or close files
1969 nrc=readbufferinfo(3,k)
1971 kfile=kfd(2,jfile)
1974 ELSE ! rewind
1976 END IF
1977 END IF
1978 END DO
1979 END IF
1980
1981 ntot=ntot+nrd
1982 nrec=ntot
1983 numreadbuffer=nrd
1984
1988
1990 WRITE(*,*) ' Record ',nrpr
1991 IF (nrpr < 100000) THEN
1992 nrpr=nrpr*10
1993 ELSE
1994 nrpr=nrpr+100000
1995 END IF
1996 END DO
1997
1999 npri=npri+1
2000 IF (npri == 1) WRITE(*,100)
2002100 FORMAT(/' PeRead records active file' &
2003 /' total block threads number')
2004101 FORMAT(' PeRead',4i10)
2005 END IF
2006
2007 IF (more <= 0) THEN
2008 ifile=0
2009 IF (floop) THEN
2010 ! check for file weights
2011 ds0=0.0_mpd
2012 ds1=0.0_mpd
2013 ds2=0.0_mpd
2014 maxrecordsize=0
2015 maxrecordfile=0
2018 maxrecordsize=xfd(k)
2019 maxrecordfile=k
2020 END IF
2021 dw=real(-kfd(1,k),mpd)
2023 ds0=ds0+dw
2024 ds1=ds1+dw*real(wfd(k),mpd)
2025 ds2=ds2+dw*real(wfd(k)**2,mpd)
2026 END DO
2027 print *, 'PEREAD: file ', maxrecordfile, 'with max record size ', maxrecordsize
2029 ds1=ds1/ds0
2030 ds2=ds2/ds0-ds1*ds1
2033177 FORMAT(/' !!!!!',i4,' weighted binary files', &
2034 /' !!!!! mean, variance of weights =',2g12.4)
2035 END DO
2036 END IF
2037 ! integrate record numbers
2040 END DO
2041 ! sort
2046 END IF
2047 END IF
2048 lprint=.false.
2049 floop=.false.
2052179 FORMAT(/' Read cache usage (#blocks, #records, ', &
2053 'min,max records/block'/17x,4i10)
2054 END IF
2055 RETURN
2056
integer(mpi), dimension(:,:), allocatable readbufferinfo
buffer management (per thread)
Definition mpmod.f90:233
◆ petime()
| subroutine petime |
Print times.
Print the elapsed and total time.
Definition at line 8839 of file pede.f90.
8841
8842 IMPLICIT NONE
8843 REAL, DIMENSION(2) :: ta
8844 REAL :: rst
8845 REAL :: delta
8846 REAL :: rstp
8847 REAL :: secnd1
8848 REAL :: secnd2
8849 INTEGER :: ncount
8850 INTEGER :: nhour1
8851 INTEGER :: minut1
8852 INTEGER :: nsecd1
8853 INTEGER :: nhour2
8854 INTEGER :: minut2
8855 INTEGER :: nsecd2
8856
8857 SAVE
8858 DATA ncount/0/
8859 ! ...
8860 ncount=ncount+1
8861 CALL etime(ta,rst)
8862 IF(ncount > 1) THEN
8863 delta=rst
8865 nhour1=nsecd1/3600
8866 minut1=nsecd1/60-60*nhour1
8867 secnd1=delta-60*(minut1+60*nhour1)
8868 delta=rst-rstp
8870 nhour2=nsecd2/3600
8871 minut2=nsecd2/60-60*nhour2
8872 secnd2=delta-60*(minut2+60*nhour2)
8873 WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2
8874 END IF
8875
8876 rstp=rst
8877 RETURN
8878101 FORMAT(i4,' h',i3,' min',f5.1,' sec total',18x,'elapsed', &
8879 i4,' h',i3,' min',f5.1,' sec')
◆ ploopa()
| subroutine ploopa | ( | integer(mpi), intent(in) | lunp | ) |
Print title for iteration.
- Parameters
-
[in] lunp unit number
Definition at line 2736 of file pede.f90.
2738
2739 IMPLICIT NONE
2740
2741 INTEGER(mpi), INTENT(IN) :: lunp
2742 ! ..
2743 WRITE(lunp,*) ' '
2744 WRITE(lunp,101) ! header line
2745 WRITE(lunp,102) ! header line
2746101 FORMAT(' it fc',' fcn_value dfcn_exp slpr costh iit st', &
2747 ' ls step cutf',1x,'rejects hhmmss FMS')
2748102 FORMAT(' -- --',' ----------- -------- ---- ----- --- --', &
2749 ' -- ----- ----',1x,'------- ------ ---')
2750 RETURN
◆ ploopb()
| subroutine ploopb | ( | integer(mpi), intent(in) | lunp | ) |
Print iteration line.
- Parameters
-
[in] lunp unit number
Definition at line 2757 of file pede.f90.
2759
2760 IMPLICIT NONE
2761 INTEGER(mpi) :: ma
2762 INTEGER :: minut
2763 INTEGER(mpi) :: nfa
2764 INTEGER :: nhour
2765 INTEGER(mpi) :: nrej
2766 INTEGER(mpi) :: nsecnd
2767 REAL(mps) :: ratae
2768 REAL :: rstb
2769 REAL(mps) :: secnd
2770 REAL(mps) :: slopes(3)
2771 REAL(mps) :: steps(3)
2772 REAL, DIMENSION(2) :: ta
2773
2774 INTEGER(mpi), INTENT(IN) :: lunp
2775
2776 CHARACTER (LEN=4):: ccalcm(4)
2777 DATA ccalcm / ' end',' S', ' F ',' FMS' /
2778 SAVE
2779
2781 IF(nrej > 9999999) nrej=9999999
2782 CALL etime(ta,rstb)
2785 nsecnd=nint(secnd,mpi)
2789 ELSE
2793 ELSE
2795 ratae=max(-99.9,min(99.9,slopes(2)/slopes(1)))
2796 stepl=steps(2)
2799 ENDIF
2800 END IF
2801103 FORMAT(i3,i3,e12.5,38x,f5.1, 1x,i7, i3,i2.2,i2.2,a4)
2802104 FORMAT(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.3,f5.1, &
2803 1x,i7, i3,i2.2,i2.2,a4)
2804105 FORMAT(i3,i3,e12.5,1x,e8.2,12x,i5,i3,9x,f5.1, &
2805 1x,i7, i3,i2.2,i2.2,a4)
2806 RETURN
◆ ploopc()
| subroutine ploopc | ( | integer(mpi), intent(in) | lunp | ) |
Print sub-iteration line.
- Parameters
-
[in] lunp unit number
Definition at line 2813 of file pede.f90.
2815
2816 IMPLICIT NONE
2817 INTEGER(mpi) :: ma
2818 INTEGER(mpi) :: minut
2819 INTEGER(mpi) :: nfa
2820 INTEGER(mpi) :: nhour
2821 INTEGER(mpi) :: nrej
2822 INTEGER(mpi) :: nsecnd
2823 REAL(mps) :: ratae
2824 REAL :: rstb
2825 REAL(mps) :: secnd
2826 REAL(mps) :: slopes(3)
2827 REAL(mps) :: steps(3)
2828 REAL, DIMENSION(2) :: ta
2829
2830 INTEGER(mpi), INTENT(IN) :: lunp
2831 CHARACTER (LEN=4):: ccalcm(4)
2832 DATA ccalcm / ' end',' S', ' F ',' FMS' /
2833 SAVE
2834
2836 IF(nrej > 9999999) nrej=9999999
2837 CALL etime(ta,rstb)
2840 nsecnd=nint(secnd,mpi)
2843 ELSE
2845 ratae=abs(slopes(2)/slopes(1))
2846 stepl=steps(2)
2849 END IF
2850104 FORMAT(3x,i3,e12.5,9x, 35x, i7, i3,i2.2,i2.2,a4)
2851105 FORMAT(3x,i3,e12.5,9x, f6.3,14x,i3,f6.3,6x, i7, i3,i2.2,i2.2,a4)
2852 RETURN
2853
◆ ploopd()
| subroutine ploopd | ( | integer(mpi), intent(in) | lunp | ) |
Print solution line.
- Parameters
-
[in] lunp unit number
Definition at line 2860 of file pede.f90.
2862 IMPLICIT NONE
2863 INTEGER :: minut
2864 INTEGER :: nhour
2865 INTEGER(mpi) :: nsecnd
2866 REAL :: rstb
2867 REAL(mps) :: secnd
2868 REAL, DIMENSION(2) :: ta
2869
2870
2871 INTEGER(mpi), INTENT(IN) :: lunp
2872 CHARACTER (LEN=4):: ccalcm(4)
2873 DATA ccalcm / ' end',' S', ' F ',' FMS' /
2874 SAVE
2875 CALL etime(ta,rstb)
2878 nsecnd=nint(secnd,mpi)
2879
2881106 FORMAT(69x,i3,i2.2,i2.2,a4)
2882 RETURN
◆ prpcon()
| subroutine prpcon |
Prepare constraints.
Count, sort constraints and split into disjoint blocks.
Definition at line 1243 of file pede.f90.
1246
1247 IMPLICIT NONE
1248 INTEGER(mpi) :: i
1249 INTEGER(mpi) :: icgb
1250 INTEGER(mpi) :: isblck
1251 INTEGER(mpi) :: ilast
1252 INTEGER(mpi) :: itgbi
1253 INTEGER(mpi) :: ivgb
1254 INTEGER(mpi) :: jcgb
1255 INTEGER(mpi) :: label
1256 INTEGER(mpi) :: labelf
1257 INTEGER(mpi) :: labell
1258 INTEGER(mpi) :: ncon
1259 INTEGER(mpi) :: npar
1260 INTEGER(mpi) :: nconmx
1261 INTEGER(mpi) :: nparmx
1262 INTEGER(mpi) :: inone
1263 INTEGER(mpi) :: itype
1264 INTEGER(mpi) :: ncgbw
1265 INTEGER(mpi) :: newlen
1266 INTEGER(mpi) :: nvar
1267 INTEGER(mpi) :: last
1268 INTEGER(mpi) :: lastlen
1269
1270 INTEGER(mpl):: length
1271 INTEGER(mpl) :: rows
1272
1273 ncgb=0
1274 ncgbw=0
1275 ncgbe=0
1277
1278 newlen=0
1279 lastlen=0
1280 nvar=-1
1281 i=0
1282 last=-1
1283 itype=0
1284 ! find next constraint header and count nr of constraints
1286 i=i+1
1287 label=listconstraints(i)%label
1288 IF(last == 0.AND.label < 0) THEN
1292 ! overwrite
1293 newlen=lastlen
1294 ! copy previous value (for label 0)
1295 newlen=newlen+1
1297 ELSE
1298 lastlen=newlen-1 ! end of last accepted constraint
1299 END IF
1301 itype=-label
1302 IF(itype == 2) ncgbw=ncgbw+1
1303 nvar=0
1304 END IF
1305 last=label
1306 IF(label > 0) THEN
1309 IF (ivgb > 0) nvar=nvar+1
1310 END IF
1311 IF(label > 0.AND.itype == 2) THEN ! weighted constraints
1314 END IF
1315 newlen=newlen+1
1318 END DO
1322 lenconstraints=newlen
1323
1327 END IF
1328 IF (ncgbw == 0) THEN
1330 ELSE
1332 END IF
1333 WRITE(*,*)
1334
1336
1337 ! keys and index for sorting of constraints
1338 length=ncgb+1; rows=3
1341 ! start of constraint in list
1344 ! start and parameter range of constraint blocks
1346
1347 ! prepare
1348 i=1
1350 ! new constraint
1351 vecconsstart(icgb)=i
1355 i=i+2
1356 DO
1357 label=listconstraints(i)%label
1360 IF(ivgb > 0) THEN
1363 END IF
1364 i=i+1
1367 END DO
1368 END DO
1369
1370 ! sort constraints
1372
1373 ! loop over sorted constraints, try to split into blocks
1374 ncblck=0
1375 nconmx=0
1376 nparmx=0
1377 mszcon=0
1378 mszprd=0
1379 isblck=1
1380 ilast=0
1382 ! index in list
1383 icgb=matconssort(3,jcgb)
1384 ! split into disjoint blocks
1385 ilast=max(ilast, matconssort(2,jcgb))
1390 END IF
1396 ncon=jcgb+1-isblck
1397 npar=ilast+1-matconssort(1,isblck)
1398 nconmx=max(nconmx,ncon)
1399 nparmx=max(nparmx,npar)
1406 ENDIF
1407 ! reset for new block
1408 isblck=jcgb+1
1409 END IF
1410 END DO
1412
1414 WRITE(*,*)
1416 WRITE(*,*) ' max block size (cons., par.) ', nconmx, nparmx
1418 END IF
1419113 FORMAT(' constraint',i6,' : ',i9,' parameters,',i9,' variable')
1420
integer(mpi) iskpec
flag for skipping empty constraints (no variable parameters)
Definition mpmod.f90:105
◆ prtglo()
| subroutine prtglo |
Print final log file.
For each global parameter:
- label (I10)
- parameter value (G14.5)
- presigma (G14.5)
- difference of parameters values (G14.5)
- difference at last iteration (G14.5)
- error (standard deviation) (G14.5)
- global correlation (F8.3), on request only
- Entries from binary files or during iterations
Definition at line 3982 of file pede.f90.
3984
3985 IMPLICIT NONE
3986 REAL(mps):: dpa
3987 REAL(mps):: err
3988 REAL(mps):: gcor
3989 INTEGER(mpi) :: i
3990 INTEGER(mpi) :: icount
3991 INTEGER(mpi) :: ie
3992 INTEGER(mpi) :: iev
3993 INTEGER(mpi) :: ij
3994 INTEGER(mpi) :: imin
3995 INTEGER(mpi) :: iprlim
3996 INTEGER(mpi) :: isub
3997 INTEGER(mpi) :: itgbi
3998 INTEGER(mpi) :: itgbl
3999 INTEGER(mpi) :: ivgbi
4000 INTEGER(mpi) :: j
4001 INTEGER(mpi) :: label
4002 INTEGER(mpi) :: lup
4003 REAL(mps):: par
4004
4005 REAL(mpd):: diag
4006 REAL(mpd)::gmati
4007 REAL(mpd)::gcor2
4008 INTEGER(mpi) :: labele(3)
4009 INTEGER(mpl):: ii
4010 REAL(mps):: compnt(3)
4011 SAVE
4012 ! ...
4013
4014 lup=09
4016
4017 WRITE(*,*) ' '
4018 WRITE(*,*) ' Result of fit for global parameters'
4019 WRITE(*,*) ' ==================================='
4020 WRITE(*,*) ' '
4021
4022 WRITE(*,101)
4023
4024 WRITE(lup,*) 'Parameter ! first 3 elements per line are', &
4025 ' significant (if used as input)'
4026
4027
4028 iprlim=10
4030 itgbl=globalparlabelindex(1,itgbi)
4031 ivgbi=globalparlabelindex(2,itgbi)
4033 icount=0 ! counts
4034 IF(ivgbi > 0) THEN
4038 ii=ivgbi
4039 ii=(ii*ii+ii)/2
4040 gmati=globalmatd(ii)
4041 err=sqrt(abs(real(gmati,mps)))
4042 IF(gmati < 0.0_mpd) err=-err
4043 diag=workspacediag(ivgbi)
4044 gcor=-1.0
4045 IF(gmati*diag > 0.0_mpd) THEN ! global correlation
4046 gcor2=1.0_mpd-1.0_mpd/(gmati*diag)
4047 IF(gcor2 >= 0.0_mpd.AND.gcor2 <= 1.0_mpd) gcor=real(sqrt(gcor2),mps)
4048 END IF
4049 END IF
4050 END IF
4052 IF(itgbi <= iprlim) THEN
4053 IF(ivgbi <= 0) THEN
4055 ELSE
4059 ELSE
4061 END IF
4062 ELSE
4064 END IF
4065 END IF
4066 ELSE IF(itgbi == iprlim+1) THEN
4067 WRITE(* ,*) '... (further printout suppressed, but see log file)'
4068 END IF
4069
4070 ! file output
4071 IF(ivgbi <= 0) THEN
4074 ELSE
4076 END IF
4077 ELSE
4083 ELSE
4085 END IF
4086 ELSE
4089 ELSE
4091 END IF
4092 END IF
4093 END IF
4094 END DO
4095 rewind lup
4096 CLOSE(unit=lup)
4097
4100 imin=1
4103 imin=i ! index of smallest pos. eigenvalue
4104 EXIT
4105 ENDIF
4106 END DO
4107 iev=0
4108
4109 DO isub=0,min(15,imin-1)
4110 IF(isub < 10) THEN
4111 i=imin-isub
4112 ELSE
4113 i=isub-9
4114 END IF
4115
4116 ! DO I=IMIN,MAX(1,IMIN-9),-1 ! backward loop, up to 10 vectors
4122 itgbi=globalparvartototal(j)
4123 label=globalparlabelindex(1,itgbi)
4124 ELSE
4126 END IF
4127 iev=iev+1
4128 labele(iev)=label
4130 IF(iev == 3) THEN
4131 WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
4132 iev=0
4133 END IF
4134 END DO
4135 IF(iev /= 0) WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
4136 iev=0
4137 WRITE(lup,*) ' '
4138 END DO
4139
4140 END IF
4141
4142101 FORMAT(1x,' label parameter presigma differ', &
4143 ' error'/ 1x,'-----------',4x,4('-------------'))
4144102 FORMAT(i10,2x,4g14.5,f8.3)
4145103 FORMAT(3(i11,f11.7,2x))
4146110 FORMAT(i10,2x,2g14.5,28x,i12)
4147111 FORMAT(i10,2x,3g14.5,14x,i12)
4148112 FORMAT(i10,2x,4g14.5,i12)
integer(mpi) igcorr
flag for output of global correlations for inversion, =0: none
Definition mpmod.f90:95
integer(mpi) ipcntr
flag for output of global parameter counts (entries), =0: none, =1: local fits, >1: binary files
Definition mpmod.f90:100
◆ prtstat()
| subroutine prtstat |
Print input statistic.
For each global parameter:
- label (I10)
- parameter value (G14.5)
- presigma (G14.5)
- difference of parameters values (G14.5), = 0.
- Entries from binary files
Definition at line 4162 of file pede.f90.
4164
4165 IMPLICIT NONE
4166 REAL(mps):: par
4167 REAL(mps):: presig
4168 INTEGER(mpi) :: icount
4169 INTEGER(mpi) :: itgbi
4170 INTEGER(mpi) :: itgbl
4171 INTEGER(mpi) :: ivgbi
4172 INTEGER(mpi) :: lup
4173 INTEGER(mpi) :: ncon
4174 INTEGER(mpi) :: k
4175
4176 SAVE
4177 ! ...
4178
4179 lup=09
4181 WRITE(lup,*) '*** Results of checking input only, no solution performed ***'
4182 WRITE(lup,*) '! fixed-1: by pre-sigma, -2: by entries cut, -3: by iterated entries cut'
4183 WRITE(lup,*) '! Label Value Pre-sigma Entries Constraints Status '
4184 !iprlim=10
4186 itgbl=globalparlabelindex(1,itgbi)
4187 ivgbi=globalparlabelindex(2,itgbi)
4192
4193 IF (ivgbi <= 0) THEN
4194 WRITE(lup,110) itgbl,par,presig,icount,ncon,ivgbi
4195 ELSE
4196 WRITE(lup,111) itgbl,par,presig,icount,ncon
4197 END IF
4198 END DO
4199 ! appearance statistics
4201 WRITE(lup,*) '! '
4202 WRITE(lup,*) '! Appearance statistics '
4203 WRITE(lup,*) '! Label First file and record Last file and record #files #paired-par'
4205 WRITE(lup,112) globalparlabelindex(1,itgbi), (appearancecounter(itgbi*5+k), k=-4,0), paircounter(itgbi)
4206 END DO
4207 END IF
4208 rewind lup
4209 CLOSE(unit=lup)
4210
4211110 FORMAT(' ! ',i10,2x,2g14.5,2i12,' fixed',i2)
4212111 FORMAT(' ! ',i10,2x,2g14.5,2i12,' variable')
4213112 FORMAT(' !.',i10,6i11)
◆ sechms()
| subroutine sechms | ( | real(mps), intent(in) | deltat, |
| integer(mpi), intent(out) | nhour, | ||
| integer(mpi), intent(out) | minut, | ||
| real(mps), intent(out) | secnd | ||
| ) |
Time conversion.
Convert from seconds to hours, minues, seconds
- Parameters
-
[in] deltat time in seconds [out] nhour hours [out] minut minutes [out] secnd seconds
Definition at line 4617 of file pede.f90.
4619
4620 IMPLICIT NONE
4621 REAL(mps), INTENT(IN) :: deltat
4622 INTEGER(mpi), INTENT(OUT) :: minut
4623 INTEGER(mpi), INTENT(OUT):: nhour
4624 REAL(mps), INTENT(OUT):: secnd
4625 INTEGER(mpi) :: nsecd
4626 ! DELTAT = time in sec -> NHOUR,MINUT,SECND
4627 ! ...
4629 nhour=nsecd/3600
4630 minut=nsecd/60-60*nhour
4631 secnd=deltat-60*(minut+60*nhour)
◆ solglo()
| subroutine solglo | ( | integer(mpi), intent(in) | ivgbi | ) |
Error for single global parameter from MINRES.
Calculate single row 'x_i' from inverse matrix by solving A*x_i=b with b=0 except b_i=1.
- Parameters
-
[in] ivgbi index of variable parameter
Definition at line 978 of file pede.f90.
981
982 IMPLICIT NONE
983 REAL(mps) :: par
984 REAL(mps) :: dpa
985 REAL(mps) :: err
986 REAL(mps) :: gcor2
987 INTEGER(mpi) :: iph
988 INTEGER(mpi) :: istop
989 INTEGER(mpi) :: itgbi
990 INTEGER(mpi) :: itgbl
991 INTEGER(mpi) :: itn
992 INTEGER(mpi) :: itnlim
993 INTEGER(mpi) :: nout
994
995 INTEGER(mpi), INTENT(IN) :: ivgbi
996
997 REAL(mpd) :: shift
998 REAL(mpd) :: rtol
999 REAL(mpd) :: anorm
1000 REAL(mpd) :: acond
1001 REAL(mpd) :: arnorm
1002 REAL(mpd) :: rnorm
1003 REAL(mpd) :: ynorm
1004 REAL(mpd) :: gmati
1005 REAL(mpd) :: diag
1006 INTEGER(mpl) :: ijadd
1007 INTEGER(mpl) :: jk
1008 INTEGER(mpl) :: ii
1009 LOGICAL :: checka
1011 SAVE
1012 DATA iph/0/
1013 ! ...
1014 IF(iph == 0) THEN
1015 iph=1
1016 WRITE(*,101)
1017 END IF
1018 itgbi=globalparvartototal(ivgbi)
1019 itgbl=globalparlabelindex(1,itgbi)
1020
1022 globalvector(ivgbi)=1.0_mpd
1023
1024 ! NOUT =6
1025 nout =0
1026 itnlim=200
1027 shift =0.0_mpd
1029 checka=.false.
1030
1031
1034 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1035
1038 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1039 ELSE
1041 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1042 END IF
1043
1044 par=real(globalparameter(itgbi),mps)
1045 dpa=real(par-globalparstart(itgbi),mps)
1046 gmati=globalcorrections(ivgbi)
1047 err=sqrt(abs(real(gmati,mps)))
1048 IF(gmati < 0.0_mpd) err=-err
1050 ii=ivgbi
1051 jk=(ii*ii+ii)/2
1053 jk=ijadd(ivgbi,ivgbi)
1054 END IF
1055 IF (jk > 0) THEN
1056 diag=globalmatd(jk)
1057 ELSE
1058 diag=real(globalmatf(-jk),mpd)
1059 END IF
1060 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps) ! global correlation (squared)
1062101 FORMAT(1x,' label parameter presigma differ', &
1063 ' Error gcor^2 iit'/ 1x,'---------',2x,5('-----------'),2x,'----')
1064102 FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
◆ solgloqlp()
| subroutine solgloqlp | ( | integer(mpi), intent(in) | ivgbi | ) |
Error for single global parameter from MINRES-QLP.
Calculate single row 'x_i' from inverse matrix by solving A*x_i=b with b=0 except b_i=1.
- Parameters
-
[in] ivgbi index of variable parameter
Definition at line 1074 of file pede.f90.
1077
1078 IMPLICIT NONE
1079 REAL(mps) :: par
1080 REAL(mps) :: dpa
1081 REAL(mps) :: err
1082 REAL(mps) :: gcor2
1083 INTEGER(mpi) :: iph
1084 INTEGER(mpi) :: istop
1085 INTEGER(mpi) :: itgbi
1086 INTEGER(mpi) :: itgbl
1087 INTEGER(mpi) :: itn
1088 INTEGER(mpi) :: itnlim
1089 INTEGER(mpi) :: nout
1090
1091 INTEGER(mpi), INTENT(IN) :: ivgbi
1092
1093 REAL(mpd) :: shift
1094 REAL(mpd) :: rtol
1095 REAL(mpd) :: mxxnrm
1096 REAL(mpd) :: trcond
1097 REAL(mpd) :: gmati
1098 REAL(mpd) :: diag
1099 INTEGER(mpl) :: ijadd
1100 INTEGER(mpl) :: jk
1101 INTEGER(mpl) :: ii
1102
1104 SAVE
1105 DATA iph/0/
1106 ! ...
1107 IF(iph == 0) THEN
1108 iph=1
1109 WRITE(*,101)
1110 END IF
1111 itgbi=globalparvartototal(ivgbi)
1112 itgbl=globalparlabelindex(1,itgbi)
1113
1115 globalvector(ivgbi)=1.0_mpd
1116
1117 ! NOUT =6
1118 nout =0
1119 itnlim=200
1120 shift =0.0_mpd
1122 mxxnrm = real(nagb,mpd)/sqrt(epsilon(mxxnrm))
1124 trcond = 1.0_mpd/epsilon(trcond) ! only QR
1126 trcond = 1.0_mpd ! only QLP
1127 ELSE
1129 END IF
1130
1133 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1134 x=globalcorrections, istop=istop, itn=itn)
1137 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1138 x=globalcorrections, istop=istop, itn=itn)
1139 ELSE
1141 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1142 x=globalcorrections, istop=istop, itn=itn)
1143 END IF
1144
1145 par=real(globalparameter(itgbi),mps)
1146 dpa=real(par-globalparstart(itgbi),mps)
1147 gmati=globalcorrections(ivgbi)
1148 err=sqrt(abs(real(gmati,mps)))
1149 IF(gmati < 0.0_mpd) err=-err
1151 ii=ivgbi
1152 jk=(ii*ii+ii)/2
1154 jk=ijadd(ivgbi,ivgbi)
1155 END IF
1156 IF (jk > 0) THEN
1157 diag=globalmatd(jk)
1158 ELSE
1159 diag=real(globalmatf(-jk),mpd)
1160 END IF
1161 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps) ! global correlation (squared)
1163101 FORMAT(1x,' label parameter presigma differ', &
1164 ' Error gcor^2 iit'/ 1x,'---------',2x,5('-----------'),2x,'----')
1165102 FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
◆ upone()
| subroutine upone |
Update, redefine hash indices.
Definition at line 4725 of file pede.f90.
4728
4729 IMPLICIT NONE
4730 INTEGER(mpi) :: i
4731 INTEGER(mpi) :: j
4732 INTEGER(mpi) :: k
4733 INTEGER(mpi) :: iprime
4734 INTEGER(mpi) :: nused
4735 LOGICAL :: finalUpdate
4736 INTEGER(mpl) :: oldLength
4737 INTEGER(mpl) :: newLength
4738 INTEGER(mpl), PARAMETER :: two = 2
4739 INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: tempArr
4740 SAVE
4741 ! ...
4743 IF(finalupdate) THEN ! final (cleanup) call
4745 END IF
4746 ! save old LabelIndex
4747 nused = globalparheader(-1)
4748 oldlength = globalparheader(-0)
4750 temparr(:,1:nused)=globalparlabelindex(:,1:nused)
4753 ! create new LabelIndex
4754 newlength = globalparheader(-3)
4757 globalparhashtable = 0
4763 ! redefine hash
4766 inner: DO
4767 k=j
4768 j=globalparhashtable(k)
4769 IF(j == 0) EXIT inner ! unused hash code
4770 IF(j == i) cycle outer ! found
4771 ENDDO inner
4772 globalparhashtable(k)=i
4773 END DO outer
4774 IF(.NOT.finalupdate) RETURN
4775
4781 END IF
◆ vmprep()
| subroutine vmprep | ( | integer(mpl), dimension(2), intent(in) | msize | ) |
Prepare storage for vectors and matrices.
- Parameters
-
[in] msize number of words for storage of global matrix (double, single prec.)
Definition at line 6034 of file pede.f90.
6037
6038 IMPLICIT NONE
6039 INTEGER(mpi) :: i
6040 INTEGER(mpi) :: ncon
6041 !
6042 INTEGER(mpl), INTENT(IN) :: msize(2)
6043
6044 INTEGER(mpl) :: length
6045 SAVE
6046 ! ...
6047 ! Vector/matrix storage
6052 length=naeqn
6061 length=nalcn
6067
6070
6072 ! array space is:
6073 ! variable-width band matrix or diagonal matrix for parameters
6074 ! followed by rectangular matrix for constraints
6075 ! followed by symmetric matrix for constraints
6078 length=nagb
6082 END DO
6085 END DO
6087 indprecond(i)=0
6088 END DO
6091 ELSE ! default preconditioner
6094 END IF
6095 END IF
6096
6097
6098 length=nagb
6100
6104
6107 ! CALL MEGARR('t D',2*NAGB,'auxiliary array') ! double aux 8
6108 END IF
6109
6116 END IF
6117
6121 END IF
6122
real(mpd), dimension(:), allocatable workspacelinesearch
workspace line search
Definition mpmod.f90:196
◆ xloopn()
| subroutine xloopn |
Standard solution algorithm.
Iterative solution. In current iteration calculate:
ICALCM = +1 Matrix, gradient, Function value & solution ICALCM = 0 gradient, Function value ICALCM = -1 solution ICALCM = -2 end
Solution is obtained by selected method and improved by line search.
Definition at line 6632 of file pede.f90.
6634
6635 IMPLICIT NONE
6636 REAL(mps) :: catio
6637 REAL(mps) :: concu2
6638 REAL(mps) :: concut
6639 REAL, DIMENSION(2) :: ta
6640 INTEGER(mpi) :: i
6641 INTEGER(mpi) :: iact
6642 INTEGER(mpi) :: iagain
6643 INTEGER(mpi) :: idx
6644 INTEGER(mpi) :: info
6645 INTEGER(mpl) :: ioff
6646 INTEGER(mpi) :: itgbi
6647 INTEGER(mpi) :: ivgbi
6648 INTEGER(mpi) :: jcalcm
6649 INTEGER(mpi) :: k
6650 INTEGER(mpi) :: labelg
6651 INTEGER(mpi) :: litera
6652 INTEGER(mpi) :: lrej
6653 INTEGER(mpi) :: lun
6654 INTEGER(mpi) :: lunp
6655 INTEGER(mpi) :: minf
6656 INTEGER(mpi) :: mrati
6657 INTEGER(mpi) :: nan
6658 INTEGER(mpi) :: nfaci
6659 INTEGER(mpi) :: nloopsol
6660 INTEGER(mpi) :: nrati
6661 INTEGER(mpi) :: nrej
6662 INTEGER(mpi) :: nsol
6663 INTEGER(mpi) :: inone
6664
6665 REAL(mpd) :: stp
6666 REAL(mpd) :: dratio
6667 REAL(mpd) :: dwmean
6668 REAL(mpd) :: db
6669 REAL(mpd) :: db1
6670 REAL(mpd) :: db2
6671 REAL(mpd) :: dbdot
6672 LOGICAL :: btest
6673 LOGICAL :: warner
6674 LOGICAL :: warners
6675 LOGICAL :: warnerss
6676 LOGICAL :: lsflag
6677 CHARACTER (LEN=7) :: cratio
6678 CHARACTER (LEN=7) :: cfacin
6679 CHARACTER (LEN=7) :: crjrat
6681 SAVE
6682 ! ...
6683
6684 ! Printout of algorithm for solution and important parameters ------
6685
6688
6690 WRITE(lunp,*) ' '
6691 WRITE(lunp,*) 'Solution algorithm: '
6692 WRITE(lunp,121) '=================================================== '
6693
6695 WRITE(lunp,121) 'solution method:','matrix inversion'
6697 WRITE(lunp,121) 'solution method:','diagonalization'
6699 WRITE(lunp,121) 'solution method:', 'minres (Paige/Saunders)'
6701 WRITE(lunp,121) 'solution method:', 'minres-qlp (Choi/Paige/Saunders)'
6703 WRITE(lunp,121) ' ', ' using QR factorization' ! only QR
6705 WRITE(lunp,121) ' ', ' using QLP factorization' ! only QLP
6706 ELSE
6707 WRITE(lunp,121) ' ', ' using QR and QLP factorization' ! QR followed by QLP
6709 END IF
6711 WRITE(lunp,121) 'solution method:', &
6712 'gmres (generalized minimzation of residuals)'
6713 END IF
6719 END IF
6722 WRITE(lunp,121) 'matrix storage:','full'
6724 WRITE(lunp,121) 'matrix storage:','sparse'
6725 END IF
6728 WRITE(lunp,121) 'pre-conditioning:','default'
6730 WRITE(lunp,121) 'pre-conditioning:','none!'
6733 WRITE(lunp,121) 'pre-conditioning=','skyline-matrix (rank preserving)'
6734 ELSE
6735 WRITE(lunp,121) 'pre-conditioning=','band-matrix'
6736 ENDIF
6737 END IF
6738 END IF
6740 WRITE(lunp,121) 'using pre-sigmas:','no'
6741 ELSE
6742 ! FIXME: NPRESG contains parameters that failed the 'entries' cut...
6743 WRITE(lunp,124) 'pre-sigmas defined for', &
6744 REAL(100*npresg,mps)/REAL(nvgb,mps),' % of variable parameters'
6746 END IF
6748 WRITE(lunp,121) 'regularization:','no'
6749 ELSE
6750 WRITE(lunp,121) 'regularization:','yes'
6752 END IF
6753
6755 WRITE(lunp,121) 'Chi square cut equiv 3 st.dev applied'
6758 WRITE(lunp,121) ' (reduced by sqrt in next iterations)'
6759 END IF
6761 WRITE(lunp,121) 'Scaling of measurement errors applied'
6764 END IF
6768 END IF
6769
6770
6771121 FORMAT(1x,a40,3x,a)
6772122 FORMAT(1x,a40,3x,i0,a)
6773123 FORMAT(1x,a40,2x,e9.2)
6774124 FORMAT(1x,a40,3x,f5.1,a)
6775 END DO
6776
6777 ! initialization of iterations -------------------------------------
6778
6779 iitera=0
6780 nsol =0 ! counter for solutions
6781 info =0
6782 lsinfo=0
6783 stp =0.0_mpd
6784 stepl =real(stp,mps)
6785 concut=1.0e-12 ! initial constraint accuracy
6786 concu2=1.0e-06 ! constraint accuracy
6789 iterat=-1
6790 litera=-2
6791 nloopsol=0 ! (new) solution from this nloopn
6792 nrej=0 ! reset number of rejects
6795 minf=1
6798 minf=2
6801 minf=3
6804 minf=3
6807 minf=3
6808 END IF
6809
6810 ! check initial feasibility of constraint equations ----------------
6811
6812 WRITE(*,*) ' '
6815 WRITE(*,*) 'Checking feasibility of parameters:'
6817 IF(iact /= 0) THEN ! done ...
6818 WRITE(*,102) concut
6819 WRITE(*,*) ' parameters are made feasible'
6822 ELSE ! ... was OK
6823 WRITE(*,*) ' parameters are feasible (i.e. satisfy constraints)'
6825 END IF
6826 concut=concu2 ! cut for constraint check
6827 END IF
6828 iact=1 ! set flag for new data loop
6830
6831 WRITE(*,*) ' '
6832 WRITE(*,*)'Reading files and accumulating vectors/matrices ...'
6833 WRITE(*,*) ' '
6834
6836 iterat=-1
6837 litera= 0
6838 jcalcm=-1
6839 iagain= 0
6840
6841 icalcm=1
6842
6843 ! Block 1: data loop with vector (and matrix) calculation ----------
6844
6845 DO
6848 IF(jcalcm+1 /= 0) THEN
6854 iterat=1
6856 ELSE
6859 ! CALL PLOOPA(LUNLOG)
6861 litera=iterat
6865 END IF
6866 ELSE
6870 END IF
6871 END IF
6872 ELSE
6875 END IF
6877 ! CHK
6878 IF (iabs(jcalcm) <= 1) THEN
6879 idx=jcalcm+4
6882 END IF
6883 END IF
6884 jcalcm=icalcm
6885
6889 lrej=nrej
6892 WRITE(*,*) ' '
6893 WRITE(*,*) 'Data rejected in previous loop: '
6894 WRITE(*,*) ' ', &
6897 WRITE(*,*) 'Too many rejects (>33.3%) - stop'
6899 stop
6900 END IF
6901 END IF
6902 ! Block 2: new iteration with calculation of solution --------------
6903
6907 END DO
6909 itgbi=globalparvartototal(i)
6911 END DO
6922 WRITE(*,*) '... reserved for GMRES (not yet!)'
6924 END IF
6926
6927 ! check feasibility and evtl. make step vector feasible
6928
6930 itgbi=globalparvartototal(i)
6933 END DO
6935 concut=concu2 ! new cut for constraint check
6937 itgbi=globalparvartototal(i)
6940 END DO
6941
6945 delfun=real(db,mps)
6946 angras=real(db/sqrt(db1*db2),mps)
6947
6948 ! do line search for this iteration/solution ?
6949 ! lsearch >2: all, =2: all with (next) chicut =1., =1: last, <1: none
6952 IF (lsflag) THEN
6953 ! initialize line search based on slopes and prepare next
6957 minf=3
6960 minf=3
6963 minf=4
6966 minf=4
6969 minf=4
6970 END IF
6971 ENDIF
6972
6973 ! change significantly negative ?
6974 IF(db <= -16.0_mpd*sqrt(max(db1,db2))*epsilon(db)) THEN
6975 WRITE(*,*) 'Function not decreasing:',db
6976 IF(db > -1.0e-3_mpd) THEN ! 100311, VB/CK: allow some margin for numerics
6977 iagain=iagain+1
6978 IF (iagain <= 1) THEN
6979 WRITE(*,*) '... again matrix calculation'
6980 icalcm=1
6981 cycle
6982 ELSE
6983 WRITE(*,*) '... aborting iterations'
6984 GO TO 90
6985 END IF
6986 ELSE
6987 WRITE(*,*) '... stopping iterations'
6988 iagain=-1
6989 GO TO 90
6990 END IF
6991 ELSE
6992 iagain=0
6993 END IF
6995 ENDIF
6996 ! Block 3: line searching ------------------------------------------
6997
6999 IF (lsflag) THEN
7004 stp, & ! returned step factor
7005 info) ! returned information
7006 ! WRITE(*,*) 'PTLINE returns INFO, STP=',INFO, STP
7007 ELSE ! skip line search
7008 info=10
7009 stepl=1.0
7012 ENDIF
7013 ENDIF
7014 lsinfo=info
7015
7016 stepl=real(stp,mps)
7017 nan=0
7019 itgbi=globalparvartototal(i)
7021 (.NOT.(workspacelinesearch(i) > 0.0_mpd))) nan=nan+1
7023 END DO
7024
7025 IF (nan > 0) THEN
7026 WRITE(*,*) 'Result vector containes ', nan,' NaNs - stop'
7028 stop
7029 END IF
7030
7031 ! subito exit, if required -----------------------------------------
7032
7034 WRITE(*,*) 'Subito! Exit after first step.'
7035 GO TO 90
7036 END IF
7037
7038 IF(info == 0) THEN
7039 WRITE(*,*) 'INFO=0 should not happen (line search input err)'
7040 IF (iagain <= 0) THEN
7041 icalcm=1
7042 cycle
7043 ENDIF
7044 END IF
7046 ! Block 4: line search convergence ---------------------------------
7047
7051 icalcm=-1
7053 cycle ! iterate
7054 END IF
7057 icalcm=-1
7059 cycle ! next iteration
7060
7061 ! Block 5: iteration ending ----------------------------------------
7062
706390 icalcm=-2
7064 END DO
7066 WRITE(*,*) ' '
7067 WRITE(*,*) 'Data rejected in last loop: '
7068 WRITE(*,*) ' ', &
7071 END IF
7072
7073 ! monitoring of residuals
7076
7079 catio=real(dratio,mps)
7081 catio=catio/0.9326 ! correction Huber downweighting (in global chi2)
7082 END IF
7083 mrati=nint(100.0*catio,mpi)
7084
7086 WRITE(lunp,*) ' '
7090 WRITE(lunp,*) ' =',dratio
7091 ELSE
7094 WRITE(lunp,*) ' /',dwmean
7095 WRITE(lunp,*) ' =',dratio
7096 END IF
7097 WRITE(lunp,*) ' '
7099 ' with correction for down-weighting ',catio
7100 END DO
7102
7103 ! ... the end with exit code ???????????????????????????????????????
7104
7105 ! WRITE(*,199) ! write exit code
7106 ! + '-----------------------------------------------------------'
7107 ! IF(ITEXIT.EQ.0) WRITE(*,199)
7108 ! + 'Exit code = 0: Convergence reached'
7109 ! IF(ITEXIT.EQ.1) WRITE(*,199)
7110 ! + 'Exit code = 1: No improvement in last iteration'
7111 ! IF(ITEXIT.EQ.2) WRITE(*,199)
7112 ! + 'Exit code = 2: Maximum number of iterations reached'
7113 ! IF(ITEXIT.EQ.3) WRITE(*,199)
7114 ! + 'Exit code = 3: Failure'
7115 ! WRITE(*,199)
7116 ! + '-----------------------------------------------------------'
7117 ! WRITE(*,199) ' '
7118
7119
7120 nrati=nint(10000.0*real(nrej,mps)/real(nrecal,mps),mpi)
7121 WRITE(crjrat,197) 0.01_mpd*real(nrati,mpd)
7122 nfaci=nint(100.0*sqrt(catio),mpi)
7123
7124 WRITE(cratio,197) 0.01_mpd*real(mrati,mpd)
7125 WRITE(cfacin,197) 0.01_mpd*real(nfaci,mpd)
7126
7127 warner=.false. ! warnings
7128 IF(mrati < 90.OR.mrati > 110) warner=.true.
7129 IF(nrati > 100) warner=.true.
7131 warners = .false. ! severe warnings
7133 warnerss = .false. ! more severe warnings
7135 IF(iagain /= 0) warnerss=.true.
7137
7138 IF(warner.OR.warners.OR.warnerss) THEN
7139 WRITE(*,199) ' '
7140 WRITE(*,199) ' '
7141 WRITE(*,199) 'WarningWarningWarningWarningWarningWarningWarningWarningWar'
7142 WRITE(*,199) 'arningWarningWarningWarningWarningWarningWarningWarningWarn'
7143 WRITE(*,199) 'rningWarningWarningWarningWarningWarningWarningWarningWarni'
7144 WRITE(*,199) 'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
7145 WRITE(*,199) 'ingWarningWarningWarningWarningWarningWarningWarningWarning'
7146 WRITE(*,199) 'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
7147 WRITE(*,199) 'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
7148
7149 IF(mrati < 90.OR.mrati > 110) THEN
7150 WRITE(*,199) ' '
7151 WRITE(*,*) ' Chi^2/Ndf = ',cratio, ' (should be close to 1)'
7152 WRITE(*,*) ' => multiply all input standard ', &
7153 'deviations by factor',cfacin
7154 END IF
7155
7156 IF(nrati > 100) THEN
7157 WRITE(*,199) ' '
7158 WRITE(*,*) ' Fraction of rejects =',crjrat,' %', &
7159 ' (should be far below 1 %)'
7160 WRITE(*,*) ' => please provide correct mille data'
7161 END IF
7162
7163 IF(iagain /= 0) THEN
7164 WRITE(*,199) ' '
7165 WRITE(*,*) ' Matrix not positiv definite '// &
7166 '(function not decreasing)'
7167 WRITE(*,*) ' => please provide correct mille data'
7168 END IF
7169
7171 WRITE(*,199) ' '
7173 ' for global matrix, should be 0'
7174 WRITE(*,*) ' => please provide correct mille data'
7175 END IF
7176
7178 WRITE(*,199) ' '
7180 ' for constraint equations, should be 0'
7181 WRITE(*,*) ' => please correct constraint definition'
7182 END IF
7183
7185 WRITE(*,199) ' '
7187 WRITE(*,*) ' => please check constraint definition, mille data'
7188 END IF
7189
7191 WRITE(*,199) ' '
7193 ' for global vector (too few entries)'
7194 WRITE(*,*) ' => please check mille data and ENTRIES cut'
7195 END IF
7196
7197 WRITE(*,199) ' '
7198 WRITE(*,199) 'WarningWarningWarningWarningWarningWarningWarningWarningWar'
7199 WRITE(*,199) 'arningWarningWarningWarningWarningWarningWarningWarningWarn'
7200 WRITE(*,199) 'rningWarningWarningWarningWarningWarningWarningWarningWarni'
7201 WRITE(*,199) 'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
7202 WRITE(*,199) 'ingWarningWarningWarningWarningWarningWarningWarningWarning'
7203 WRITE(*,199) 'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
7204 WRITE(*,199) 'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
7205 WRITE(*,199) ' '
7206
7207 ENDIF
7208
7210
7211 ! ------------------------------------------------------------------
7212
7214
7218 ! errors and correlations from MINRES
7220 labelg=lbmnrs(k)
7221 IF(labelg == 0) cycle
7222 itgbi=inone(labelg)
7223 ivgbi=0
7225 IF(ivgbi < 0) ivgbi=0
7226 IF(ivgbi == 0) cycle
7227 ! determine error and global correlation for parameter IVGBI
7230 ELSE
7232 ENDIF
7233 END DO
7234
7236
7237 END IF
7238
7240 !use elimination for constraints ?
7242 ! extend, transform matrix
7244 ioff=((i-1)*i)/2
7245 globalmatd(ioff+1:ioff+i)=0.0_mpd
7246 END DO
7248 END IF
7249 END IF
7250
7252
7253 IF (warnerss) THEN
7255 ELSE IF (warners) THEN
7257 ELSE IF (warner) THEN
7259 ELSE
7261 END IF
7262
7263102 FORMAT(' Call FEASIB with cut=',g10.3)
7264 ! 103 FORMAT(1X,A,G12.4)
7265197 FORMAT(f7.2)
7266199 FORMAT(7x,a)
integer(mpi) mnrsel
number of MINRES error labels in LBMNRS (calc err, corr with SOLGLO)
Definition mpmod.f90:89
subroutine zdiags
Covariance matrix for diagonalization (,correction of eigenvectors).
Definition pede.f90:6349
subroutine solgloqlp(ivgbi)
Error for single global parameter from MINRES-QLP.
Definition pede.f90:1075
◆ zdiags()
| subroutine zdiags |
Covariance matrix for diagonalization (,correction of eigenvectors).
Definition at line 6348 of file pede.f90.
6350
6351 IMPLICIT NONE
6352 INTEGER(mpi) :: i
6353 INTEGER(mpl) :: ioff1
6354 INTEGER(mpl) :: ioff2
6355 INTEGER(mpi) :: j
6356
6357 ! eigenvalue eigenvectors cov.matrix
6359
6360 !use elimination for constraints ?
6362 ! extend, transform eigenvectors
6365 workspaceeigenvectors(ioff2+1:)=0.0_mpd
6367 ioff1=ioff1-nfgb
6368 ioff2=ioff2-nvgb
6371 END DO
6373 END DO
6375 END IF
6376
subroutine devinv(n, diag, u, v)
Inversion by diagonalization. Get inverse matrix V from DIAG and U.
Definition mpnum.f90:705
