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

some fixes

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@141 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent 760bee22
......@@ -1785,7 +1785,10 @@ SUBROUTINE lltdec(n,c,india,nrkd)
REAL(mpd), INTENT(IN OUT) :: c(*)
INTEGER(mpi), INTENT(IN) :: india(n)
INTEGER(mpi), INTENT(OUT) :: nrkd
REAL(mpd) eps
! ...
eps = 16.0_mpd * epsilon(eps) ! 16 * precision(mpd)
! ...
nrkd=0
diag=0.0_mpd
......@@ -1812,7 +1815,7 @@ SUBROUTINE lltdec(n,c,india,nrkd)
IF(j /= k) c(kj)=c(kj)*diag
END DO ! J
IF(diag+c(india(k)) > diag) THEN ! test for linear dependence
IF(c(india(k)) > eps*diag) THEN ! test for linear dependence
c(india(k))=1.0_mpd/SQRT(c(india(k))) ! square root
ELSE
DO j=mk,k ! reset row K
......@@ -1862,6 +1865,7 @@ SUBROUTINE lltfwd(n,c,india,x)
END DO ! J
x(k)=x(k)*c(india(k))
END DO ! K
RETURN
END SUBROUTINE lltfwd
......@@ -1882,13 +1886,13 @@ SUBROUTINE lltbwd(n,c,india,x)
IMPLICIT NONE
INTEGER(mpi) :: j
INTEGER(mpi) :: k
INTEGER(mpi) :: k
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(IN) :: c(*)
INTEGER(mpi), INTENT(IN) :: india(n)
REAL(mpd), INTENT(IN OUT) :: x(n)
DO k=n,2,-1 ! backward loop
x(k)=x(k)*c(india(k))
DO j=k-india(k)+india(k-1)+1,k-1
......@@ -1896,6 +1900,8 @@ SUBROUTINE lltbwd(n,c,india,x)
END DO ! J
END DO ! K
x(1)=x(1)*c(india(1))
RETURN
END SUBROUTINE lltbwd
!> Decomposition of equilibrium systems.
......@@ -1922,7 +1928,6 @@ SUBROUTINE equdec(n,m,c,india,nrkd,nrkd2)
INTEGER(mpi) :: j
INTEGER(mpi) :: jk
INTEGER(mpi) :: k
INTEGER(mpi) :: ntotal
INTEGER(mpi), INTENT(IN) :: n
INTEGER(mpi), INTENT(IN) :: m
......@@ -1931,33 +1936,35 @@ SUBROUTINE equdec(n,m,c,india,nrkd,nrkd2)
INTEGER(mpi), INTENT(OUT) :: nrkd
INTEGER(mpi), INTENT(OUT) :: nrkd2
! ...
ntotal=n+n*m+(m*m+m)/2
! ...
CALL lltdec(n,c,india,nrkd) ! decomposition G G^T
DO i=1,m
CALL lltfwd(n,c,india,c(india(n)+(i-1)*n+1)) ! forward solution K
END DO
IF (m>0) THEN
DO i=1,m
CALL lltfwd(n,c,india,c(india(n)+(i-1)*n+1)) ! forward solution K
END DO
jk=india(n)+n*m
DO j=1,m
DO k=1,j
jk=jk+1
c(jk)=0.0_mpd ! product K K^T
DO i=1,n
c(jk)=c(jk)+c(india(n)+(j-1)*n+i)*c(india(n)+(k-1)*n+i)
jk=india(n)+n*m
DO j=1,m
DO k=1,j
jk=jk+1
c(jk)=0.0_mpd ! product K K^T
DO i=1,n
c(jk)=c(jk)+c(india(n)+(j-1)*n+i)*c(india(n)+(k-1)*n+i)
END DO
END DO
END DO
END DO
india(n+1)=1
DO i=2,m
india(n+i)=india(n+i-1)+MIN(i,m) ! pointer for K K^T
END DO
CALL lltdec(m,c(india(n)+n*m+1),india(n+1),nrkd2) ! decomp. H H^T
india(n+1)=1
DO i=2,m
india(n+i)=india(n+i-1)+MIN(i,m) ! pointer for K K^T
END DO
ntotal=n+n*m+(m*m+m)/2
CALL lltdec(m,c(india(n)+n*m+1),india(n+1),nrkd2) ! decomp. H H^T
ELSE
nrkd2=0
ENDIF
RETURN
END SUBROUTINE equdec
......@@ -1989,27 +1996,33 @@ SUBROUTINE equslv(n,m,c,india,x) ! solution vector
REAL(mpd), INTENT(IN) :: c(*)
INTEGER(mpi), INTENT(IN) :: india(n+m)
REAL(mpd), INTENT(IN OUT) :: x(n+m)
CALL lltfwd(n,c,india,x) ! result is u
DO i=1,m
DO j=1,n
x(n+i)=x(n+i)-x(j)*c(india(n)+(i-1)*n+j) ! g - K u
IF (m>0) THEN
DO i=1,m
DO j=1,n
x(n+i)=x(n+i)-x(j)*c(india(n)+(i-1)*n+j) ! g - K u
END DO
END DO
END DO
CALL lltfwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is v
CALL lltfwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is v
CALL lltbwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is -y
DO i=1,m
x(n+i)=-x(n+i) ! result is +y
END DO
CALL lltbwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is -y
DO i=1,m
x(n+i)=-x(n+i) ! result is +y
END DO
DO i=1,n
DO j=1,m
x(i)=x(i)-x(n+j)*c(india(n)+(j-1)*n+i) ! u - K^T y
DO i=1,n
DO j=1,m
x(i)=x(i)-x(n+j)*c(india(n)+(j-1)*n+i) ! u - K^T y
END DO
END DO
END DO
ENDIF
CALL lltbwd(n,c,india,x) ! result is x
RETURN
END SUBROUTINE equslv
!> Constrained preconditioner, decomposition.
......@@ -2040,8 +2053,9 @@ END SUBROUTINE equslv
!! \param [out] cu 1/sqrt(c)
!! \param [in,out] a constraint matrix (size n*p), modified
!! \param [out] s Cholesky decomposed symmetric (P,P) matrix
!! \param [out] nrkd removed components
SUBROUTINE precon(p,n,c,cu,a,s)
SUBROUTINE precon(p,n,c,cu,a,s,nrkd)
USE mpdef
IMPLICIT NONE
......@@ -2059,10 +2073,12 @@ SUBROUTINE precon(p,n,c,cu,a,s)
REAL(mpd), INTENT(OUT) :: cu(n)
REAL(mpd), INTENT(IN OUT) :: a(n,p)
REAL(mpd), INTENT(OUT) :: s((p*p+p)/2)
INTEGER(mpi), INTENT(OUT) :: nrkd
REAL(mpd) :: div
REAL(mpd) :: ratio
nrkd=0
DO i=1,(p*p+p)/2
s(i)=0.0_mpd
END DO
......@@ -2073,6 +2089,7 @@ SUBROUTINE precon(p,n,c,cu,a,s)
cu(i)=1.0_mpd/SQRT(div)
ELSE
cu(i)=0.0_mpd
nrkd=nrkd+1
END IF
DO j=1,p
a(i,j)=a(i,j)*cu(i) ! K = A C^{-1/2}
......
......@@ -4935,7 +4935,7 @@ SUBROUTINE loop2
IF (nvar == 0 .AND. iskpec > 0) newlen=lastlen
END IF
lenConstraints=newlen
IF (ncgbe > 0) THEN
IF (ncgbe > 0 .AND. iskpec > 0) THEN
WRITE(*,*) 'LOOP2:',ncgbe,' empty constraints skipped'
ncgb=ncgb-ncgbe
END IF
......@@ -5884,18 +5884,18 @@ SUBROUTINE mminrs
IF(mbandw == 0) THEN ! default preconditioner
IF(icalcm == 1) THEN
IF(nfgb < nvgb) CALL qlpssq(avprd0,matPreCond,mbandw,.true.) ! transform preconditioner matrix
IF(nfgb < nvgb) CALL qlpssq(avprd0,matPreCond,1,.true.) ! transform preconditioner matrix
CALL precon(nprecond(1),nprecond(2),matPreCond,matPreCond, matPreCond(1+nvgb), &
matPreCond(1+nvgb+ncgb*nvgb))
matPreCond(1+nvgb+ncgb*nvgb),nrkd)
END IF
CALL minres(nfgb, avprod, mcsolv, workspaceD, shift, checka ,.TRUE. , &
globalCorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
ELSE IF(mbandw > 0) THEN ! band matrix preconditioner
IF(icalcm == 1) THEN
IF(nfgb < nvgb) CALL qlpssq(avprd0,matPreCond,mbandw,.true.) ! transform preconditioner matrix
WRITE(lun,*) 'MMINRS: EQUDEC started'
WRITE(lun,*) 'MMINRS: EQUDEC started', nprecond(2), nprecond(1)
CALL equdec(nprecond(2),nprecond(1),matPreCond,indPreCond,nrkd,nrkd2)
WRITE(lun,*) 'MMINRS: EQUDEC ended'
WRITE(lun,*) 'MMINRS: EQUDEC ended ', nrkd, nrkd2
END IF
CALL minres(nfgb, avprod, mvsolv, workspaceD, shift, checka ,.TRUE. , &
globalCorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
......@@ -5979,9 +5979,9 @@ SUBROUTINE mminrsqlp
IF(mbandw == 0) THEN ! default preconditioner
IF(icalcm == 1) THEN
IF(nfgb < nvgb) CALL qlpssq(avprd0,matPreCond,mbandw,.true.) ! transform preconditioner matrix
IF(nfgb < nvgb) CALL qlpssq(avprd0,matPreCond,1,.true.) ! transform preconditioner matrix
CALL precon(nprecond(1),nprecond(2),matPreCond,matPreCond, matPreCond(1+nvgb), &
matPreCond(1+nvgb+ncgb*nvgb))
matPreCond(1+nvgb+ncgb*nvgb),nrkd)
END IF
CALL minresqlp( n=nfgb, Aprod=avprod, b=workspaceD, Msolve=mcsolv, nout=nout, &
itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
......@@ -5989,9 +5989,9 @@ SUBROUTINE mminrsqlp
ELSE IF(mbandw > 0) THEN ! band matrix preconditioner
IF(icalcm == 1) THEN
IF(nfgb < nvgb) CALL qlpssq(avprd0,matPreCond,mbandw,.true.) ! transform preconditioner matrix
WRITE(lun,*) 'MMINRS: EQUDEC started'
WRITE(lun,*) 'MMINRS: EQUDEC started', nprecond(2), nprecond(1)
CALL equdec(nprecond(2),nprecond(1),matPreCond,indPreCond,nrkd,nrkd2)
WRITE(lun,*) 'MMINRS: EQUDEC ended'
WRITE(lun,*) 'MMINRS: EQUDEC ended ', nrkd, nrkd2
END IF
CALL minresqlp( n=nfgb, Aprod=avprod, b=workspaceD, Msolve=mvsolv, nout=nout, &
......
......@@ -3,7 +3,7 @@
## \file
# Read millepede binary file and print records
#
# \author Claus Kleinwort, DESY, 2009 (Claus.Kleinwort@desy.de)
# \author Claus Kleinwort, DESY, 2009 (Claus.Kleinwort@desy.de)
#
# \copyright
# Copyright (c) 2009 - 2015 Deutsches Elektronen-Synchroton,
......@@ -23,8 +23,9 @@
#
# Hardcoded defaults can be replaced by command line arguments for
# - Name of binary file
# - Number of records to print
# - Number of records to print (-1: all)
# - Number of records to skip (optional)
# - Mininum value to print derivative
#
# Description of the output from readMilleBinary.py
# - Records (tracks) start with \c '===' followed by record number and length
......@@ -60,6 +61,8 @@ fname = "milleBinaryISN.dat"
mrec = 10
## number of records (track) to skip before
skiprec = 0
## minimum value to print derivatives
minval = 0.
#
# ## C. Kleinwort - DESY ########################
......@@ -67,22 +70,25 @@ skiprec = 0
narg = len(sys.argv)
if narg > 1:
if narg < 3:
print " usage: readMilleBinary.py <file name> <number of records> [<number of records to skip>]"
print " usage: readMilleBinary.py <file name> <number of records> [<number of records to skip> <minimum value to print derivative>]"
sys.exit(2)
else:
fname = sys.argv[1]
mrec = int(sys.argv[2])
if narg > 3:
skiprec = int(sys.argv[3])
if narg > 4:
minval = float(sys.argv[4])
#print " input ", fname, mrec, skiprec
f = open(fname, "rb")
nrec = 0
try:
while (nrec < mrec + skiprec):
# read 1 record
while (nrec < mrec + skiprec) or (mrec < 0):
# read 1 record
nr = 0
if (Cfiles == 0):
lenf = array.array(intfmt)
lenf.fromfile(f, 2)
......@@ -141,15 +147,31 @@ try:
# measurement without global derivatives
print ' -l- meas. ', nh, inder[ja + 1], jb - ja - 1, i - jb, glder[ja], glder[jb]
if (ja + 1 < jb):
print " local ", inder[ja + 1:jb]
print " local ", glder[ja + 1:jb]
lab = []
val = []
for k in range(ja+1, jb):
if abs(glder[k]) >= minval:
lab.append(inder[k])
val.append(glder[k])
print " local ", lab
print " local ", val
if (jb + 1 < i + 1):
print " global ", inder[jb + 1:i + 1]
print " global ", glder[jb + 1:i + 1]
lab = []
val = []
for k in range(jb+1, i+1):
if abs(glder[k]) >= minval:
lab.append(inder[k])
val.append(glder[k])
print " global ", lab
print " global ", val
except EOFError:
pass
# print "end of file"
print
if (nr > 0):
print " >>> error: end of file before end of record", nrec
else:
print " end of file after", nrec, "records"
f.close()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment