Skip to content
GitLab
Commit 4c587a8c authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Solution by elimination for constraints implemented 

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@138 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent 1d88ddca
Expand all Hide whitespace changes
Inline Side-by-side
Makefile
View file @ 4c587a8c
......@@ -2,7 +2,7 @@
# Makefile for MillePede II (Fortran90) with possible input from C
#
# Author Volker Blobel, University Hamburg, 2005-2009
# Updates by Gero Flucke, CERN, Claus Kleinwort, DESY
# Updates by Gero Flucke, CERN, Claus Kleinwort, DESY
#
# Tested on - 64-bit SL5 with gcc version 4.4.4.
# - 64 bit Ubuntu 11.10 with gcc version 4.6.1.
......@@ -84,7 +84,7 @@ L_FLAGS = -Wall -O3
#
# objects for this project
#
USER_OBJ_PEDE = mpdef.o mpdalc.o mpmod.o mpbits.o mptest1.o mptest2.o mille.o mpnum.o mptext.o mphistab.o \
USER_OBJ_PEDE = mpdef.o mpdalc.o mpmod.o mpbits.o mpqldec.o mptest1.o mptest2.o mille.o mpnum.o mptext.o mphistab.o \
minresDataModule.o minresModule.o minresqlpDataModule.o minresqlpBlasModule.o minresqlpModule.o \
randoms.o vertpr.o linesrch.o Dbandmatrix.o pede.o
#
......@@ -138,6 +138,7 @@ mpbits.o: mpdef.o mpdalc.o
mpdalc.o: mpdef.o
mpmod.o: mpdef.o
mpnum.o: mpdef.o
mpqldec.o: mpdef.o mpdalc.o
pede.o: mpdef.o mpmod.o mpdalc.o mptest1.o mptest2.o mptext.o
# ##################################################################
# END
......
mpmod.f90
View file @ 4c587a8c
......@@ -98,13 +98,16 @@ MODULE mpmod
!! >2: all, =2: all with (next) Chi2 cut scaling factor =1., =1: last, <1: none
INTEGER(mpi) :: ipcntr=0 !< flag for output of global parameter counts (entries), =0: none, =1: local fits, >1: binary files
INTEGER(mpi) :: iwcons=0 !< flag for weighting of constraints (>0: weighting with \ref globalparcounts "globalParCounts", else: none)
INTEGER(mpi) :: icelim=1 !< flag for using elimination (instead of multipliers) for constraints
! variables
INTEGER(mpi) :: lunlog !< unit for logfile
INTEGER(mpi) :: lvllog !< log level
INTEGER(mpi) :: ntgb !< total number of global parameters
INTEGER(mpi) :: nvgb !< number of variable global parameters
INTEGER(mpi) :: nagb !< number of fit parameters (global par. + Lagrange mult.)
INTEGER(mpi) :: nagb !< number of all parameters (global par. + Lagrange mult.)
INTEGER(mpi) :: nfgb !< number of fit parameters
INTEGER(mpi) :: ncgb !< number of constraints
INTEGER(mpi), DIMENSION(2) :: nprecond !< number of constraints, matrix size for preconditioner
INTEGER(mpi) :: nagbn !< max number of global paramters per record
INTEGER(mpi) :: nalcn !< max number of local paramters per record
INTEGER(mpi) :: naeqn !< max number of equations (measurements) per record
......@@ -162,11 +165,15 @@ MODULE mpmod
REAL(mps), DIMENSION(:), ALLOCATABLE :: globalMatF !< global matrix 'A' (float part for compressed sparse)
REAL(mpd), DIMENSION(:), ALLOCATABLE :: globalVector !< global vector 'x' (in A*x=b)
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: globalCounter !< global counter (entries in 'x')
! AVPROD (A*x=b) by MINRES
REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecXav !< vector x for AVPROD (A*x=b)
REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecBav !< vector b for AVPROD (A*x=b)
! preconditioning
REAL(mpd), DIMENSION(:), ALLOCATABLE :: matPreCond !< preconditioner (band) matrix
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: indPreCond !< preconditioner pointer array
! auxiliary vectors
REAL(mpd), DIMENSION(:), ALLOCATABLE :: workspaceD !< (general) workspace (D)
REAL(mpd), DIMENSION(:), ALLOCATABLE :: workspaceDiag !< diagonal of global matrix (for global corr.)
REAL(mpd), DIMENSION(:), ALLOCATABLE :: workspaceLinesearch !< workspace line search
REAL(mpd), DIMENSION(:), ALLOCATABLE :: workspaceDiagonalization !< workspace diag.
REAL(mpd), DIMENSION(:), ALLOCATABLE :: workspaceEigenValues !< workspace eigen values
......@@ -175,6 +182,7 @@ MODULE mpmod
! constraint matrix, residuals
REAL(mpd), DIMENSION(:), ALLOCATABLE :: matConsProduct !< product matrix of constraints
REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecConsResiduals !< residuals of constraints
REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecConsSolution !< solution for constraint elimination
! global parameter mapping
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: globalParLabelIndex !< global parameters label, total -> var. index
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: globalParHashTable !< global parameters hash table
......
mpqldec.f90 0 → 100644
View file @ 4c587a8c
!> \file
!! QL decompostion.
!!
!! \author Claus Kleinwort, DESY, 2015 (Claus.Kleinwort@desy.de)
!!
!! \copyright
!! Copyright (c) 2015 Deutsches Elektronen-Synchroton,
!! Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
!! This library is free software; you can redistribute it and/or modify
!! it under the terms of the GNU Library General Public License as
!! published by the Free Software Foundation; either version 2 of the
!! License, or (at your option) any later version. \n\n
!! This library is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU Library General Public License for more details. \n\n
!! You should have received a copy of the GNU Library General Public
!! License along with this program (see the file COPYING.LIB for more
!! details); if not, write to the Free Software Foundation, Inc.,
!! 675 Mass Ave, Cambridge, MA 02139, USA.
!!
!! QL decomposition of constraints matrix by Householder transformations
!! for solution by elimination.
!!
!> QL data.
MODULE mpqldec
USE mpdef
IMPLICIT NONE
INTEGER(mpi) :: npar !< number of parameters
INTEGER(mpi) :: ncon !< number of constraints
REAL(mpd), DIMENSION(:), ALLOCATABLE :: matV !< unit normals (v_i) of Householder reflectors
REAL(mpd), DIMENSION(:), ALLOCATABLE :: matL !< lower diagonal matrix L
END MODULE mpqldec
!> Initialize QL decomposition.
!!
!! \param [in] n number of rows (parameters)
!! \param [in] m number of columns (constraints)
!!
SUBROUTINE qlini(n,m)
USE mpqldec
USE mpdalc
IMPLICIT NONE
INTEGER(mpl) :: length
INTEGER(mpi), INTENT(IN) :: n
INTEGER(mpi), INTENT(IN) :: m
npar=n
ncon=m
! allocate
length=npar*ncon
CALL mpalloc(matV,length,'QLDEC: V')
length=ncon*ncon
CALL mpalloc(matL,length,'QLDEC: L')
END SUBROUTINE qlini
! 141217 C. Kleinwort, DESY-FH1
!> QL decomposition.
!!
!! QL decomposition with Householder transformations.
!! Decompose N-By-M matrix A into orthogonal N-by-N matrix Q and a
!! N-by-M matrix containing zeros except for a lower triangular
!! M-by-M matrix L (at the bottom):
!!
!! | 0 |
!! A = Q * | |
!! | L |
!!
!! The decomposition is stored in a N-by-M matrix matV containing the unit
!! normal vectors v_i of the hyperplanes (Householder reflectors) defining Q.
!! The lower triangular matrix L is stored in the M-by-M matrix matL.
!!
!! \param [in] a Npar-by-Ncon matrix
!!
SUBROUTINE qldec(a)
USE mpqldec
USE mpdalc
! cost[dot ops] ~= Npar*Ncon*Ncon
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpl) :: ioff1
INTEGER(mpl) :: ioff2
INTEGER(mpl) :: ioff3
INTEGER(mpi) :: k
INTEGER(mpi) :: kn
INTEGER(mpl) :: length
REAL(mpd) :: nrm
REAL(mpd) :: sp
REAL(mpd), INTENT(IN) :: a(*)
REAL(mpd) :: v(npar)
! prepare
length=npar*ncon
matV=a(1:length)
matL=0.0_mpd
! Householder procedure
DO k=ncon,1,-1
kn=npar+k-ncon
! column offset
ioff1=(k-1)*npar
! get column
v(1:kn)=matV(ioff1+1:ioff1+kn)
nrm = SQRT(dot_product(v(1:kn),v(1:kn)))
IF (nrm == 0.0_mpd) CYCLE
!
IF (v(kn) >= 0.0_mpd) THEN
v(kn)=v(kn)+nrm
ELSE
v(kn)=v(kn)-nrm
END IF
! create normal vector
nrm = SQRT(dot_product(v(1:kn),v(1:kn)))
v(1:kn)=v(1:kn)/nrm
! transformation
ioff2=0
DO i=1,k
sp=dot_product(v(1:kn),matV(ioff2+1:ioff2+kn))
matV(ioff2+1:ioff2+kn)=matV(ioff2+1:ioff2+kn)-2.0_mpd*v(1:kn)*sp
ioff2=ioff2+npar
END DO
! store column of L
ioff3=(k-1)*ncon
matL(ioff3+k:ioff3+ncon)=matV(ioff1+kn:ioff1+npar)
! store normal vector
matV(ioff1+1:ioff1+kn)=v(1:kn)
matV(ioff1+kn+1:ioff1+npar)=0.0_mpd
END DO
END SUBROUTINE qldec
!> Multiply left by Q(t).
!!
!! Multiply left by Q(t) from QL decomposition.
!!
!! \param [in,out] x Npar-by-M matrix, overwritten with Q*X (t=false) or Q^t*X (t=true)
!! \param [in] m number of columns
!! \param [in] t use transposed of Q
!!
SUBROUTINE qlmlq(x,m,t)
USE mpqldec
! cost[dot ops] ~= N*M*Nhr
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpl) :: ioff1
INTEGER(mpl) :: ioff2
INTEGER(mpi) :: j
INTEGER(mpi) :: k
INTEGER(mpi) :: kn
REAL(mpd) :: sp
REAL(mpd), INTENT(IN OUT) :: x(*)
INTEGER(mpi), INTENT(IN) :: m
LOGICAL, INTENT(IN) :: t
DO j=1,ncon
k=j
IF (t) k=ncon+1-j
kn=npar+k-ncon
! column offset
ioff1=(k-1)*npar
! transformation
ioff2=0
DO i=1,m
sp=dot_product(matV(ioff1+1:ioff1+kn),x(ioff2+1:ioff2+kn))
x(ioff2+1:ioff2+kn)=x(ioff2+1:ioff2+kn)-2.0_mpd*matV(ioff1+1:ioff1+kn)*sp
ioff2=ioff2+npar
END DO
END DO
END SUBROUTINE qlmlq
!> Multiply right by Q(t).
!!
!! Multiply right by Q(t) from QL decomposition.
!!
!! \param [in,out] x M-by-Npar matrix, overwritten with X*Q (t=false) or X*Q^t (t=true)
!! \param [in] m number of rows
!! \param [in] t use transposed of Q
!!
SUBROUTINE qlmrq(x,m,t)
USE mpqldec
! cost[dot ops] ~= N*M*Nhr
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpl) :: ioff1
INTEGER(mpl) :: iend
INTEGER(mpi) :: j
INTEGER(mpi) :: k
INTEGER(mpi) :: kn
REAL(mpd) :: sp
REAL(mpd), INTENT(IN OUT) :: x(*)
INTEGER(mpi), INTENT(IN) :: m
LOGICAL, INTENT(IN) :: t
DO j=1,ncon
k=j
IF (.not.t) k=ncon+1-j
kn=npar+k-ncon
! column offset
ioff1=(k-1)*npar
! transformation
iend=m*kn
DO i=1,npar
sp=dot_product(matV(ioff1+1:ioff1+kn),x(i:iend:m))
x(i:iend:m)=x(i:iend:m)-2.0_mpd*matV(ioff1+1:ioff1+kn)*sp
END DO
END DO
END SUBROUTINE qlmrq
!> Similarity transformation by Q(t).
!!
!! Similarity transformation by Q from QL decomposition.
!!
!! \param [in,out] x Npar-by-Npar matrix, overwritten with Q*X*Q^t (t=false) or Q^t*X*Q (t=true)
!! \param [in] t use transposed of Q
!!
SUBROUTINE qlsmq(x,t)
USE mpqldec
! cost[dot ops] ~= N*N*Nhr
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpl) :: ioff1
INTEGER(mpl) :: ioff2
INTEGER(mpl) :: iend
INTEGER(mpi) :: j
INTEGER(mpi) :: k
INTEGER(mpi) :: kn
REAL(mpd) :: sp
REAL(mpd), INTENT(IN OUT) :: x(*)
LOGICAL, INTENT(IN) :: t
DO j=1,ncon
k=j
IF (t) k=ncon+1-j
kn=npar+k-ncon
! column offset
ioff1=(k-1)*npar
! transformation
iend=npar*kn
DO i=1,npar
sp=dot_product(matV(ioff1+1:ioff1+kn),x(i:iend:npar))
x(i:iend:npar)=x(i:iend:npar)-2.0_mpd*matV(ioff1+1:ioff1+kn)*sp
END DO
ioff2=0
DO i=1,npar
sp=dot_product(matV(ioff1+1:ioff1+kn),x(ioff2+1:ioff2+kn))
x(ioff2+1:ioff2+kn)=x(ioff2+1:ioff2+kn)-2.0_mpd*matV(ioff1+1:ioff1+kn)*sp
ioff2=ioff2+npar
END DO
END DO
END SUBROUTINE qlsmq
!> Similarity transformation by Q(t).
!!
!! Similarity transformation for symmetric matrix by Q from QL decomposition.
!!
!! \param [in] aprod external procedure to calculate A*v
!! \param [in,out] A symmetric Npar-by-Npar matrix A in symmetric storage mode
!! (V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, ...),
!! overwritten with Q*A*Q^t (t=false) or Q^t*A*Q (t=true)
!! \param [in] t use transposed of Q
!!
SUBROUTINE qlssq(aprod,A,t)
USE mpqldec
USE mpdalc
! cost[dot ops] ~= N*N*Nhr
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpl) :: ioff1
INTEGER(mpl) :: ioff2
INTEGER(mpi) :: j
INTEGER(mpi) :: k
INTEGER(mpi) :: kn
INTEGER(mpi) :: l
INTEGER(mpl) :: length
REAL(mpd) :: vtAv
REAL(mpd), DIMENSION(:), ALLOCATABLE :: Av
REAL(mpd), INTENT(IN OUT) :: A(*)
LOGICAL, INTENT(IN) :: t
INTERFACE
SUBROUTINE aprod(n,x,y) ! y=A*x
USE mpdef
INTEGER(mpi), INTENT(in) :: n
REAL(mpd), INTENT(IN) :: x(n)
REAL(mpd), INTENT(OUT) :: y(n)
END SUBROUTINE aprod
END INTERFACE
length=npar
CALL mpalloc(Av,length,'qlssq: A*v')
DO j=1,ncon
k=j
IF (t) k=ncon+1-j
kn=npar+k-ncon
! column offset
ioff1=(k-1)*npar
! A*v
CALL aprod(npar,matV(ioff1+1:ioff1+npar),Av(1:npar))
! transformation
! diagonal block
! v^t*A*v
vtAv=dot_product(matV(ioff1+1:ioff1+kn),Av(1:kn))
! update
ioff2=0
DO i=1,kn
! correct with 2*(2v*vtAv*v^t - Av*v^t - (Av*v^t)^t)
DO l=1,i
ioff2=ioff2+1
A(ioff2)=A(ioff2)+2.0_mpd*((2.0_mpd*matV(ioff1+i)*vtAv-Av(i))*matV(ioff1+l)-Av(l)*matV(ioff1+i))
END DO
END DO
! off diagonal block
DO i=kn+1,npar
! correct with -2Av*v^t
A(ioff2+1:ioff2+kn)=A(ioff2+1:ioff2+kn)-2.0_mpd*matV(ioff1+1:ioff1+kn)*Av(i)
ioff2=ioff2+i
END DO
END DO
CALL mpdealloc(Av)
END SUBROUTINE qlssq
!> Partial similarity transformation by Q(t).
!!
!! Partial similarity transformation for symmetric matrix by Q from QL decomposition.
!! Calculate corrections to band part of matrix.
!!
!! \param [in] aprod external procedure to calculate A*v
!! \param [in,out] B band part of symmetric Npar-by-Npar matrix A in symmetric storage mode,
!! overwritten with band part of Q^t*A*Q (t=false) or Q^t*A*Q (t=true)
!! \param [in] m band width (including diagonal)
!! \param [in] t use transposed of Q
!!
SUBROUTINE qlpssq(aprod,B,m,t)
USE mpqldec
USE mpdalc
! cost[dot ops] ~= N*N*Nhr
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpl) :: ioff1
INTEGER(mpl) :: ioff2
INTEGER(mpi) :: j
INTEGER(mpi) :: j2
INTEGER(mpi) :: k
INTEGER(mpi) :: k2
INTEGER(mpi) :: kn
INTEGER(mpi) :: l
INTEGER(mpl) :: length
INTEGER(mpi) :: mbnd
REAL(mpd) :: vtAv
REAL(mpd) :: vtAvp
REAL(mpd) :: vtvp
REAL(mpd), DIMENSION(:), ALLOCATABLE :: Av ! A*v
REAL(mpd), INTENT(IN OUT) :: B(*)
INTEGER(mpi), INTENT(IN) :: m
LOGICAL, INTENT(IN) :: t
INTERFACE
SUBROUTINE aprod(n,x,y) ! y=A*x
USE mpdef
INTEGER(mpi), INTENT(in) :: n
REAL(mpd), INTENT(IN) :: x(n)
REAL(mpd), INTENT(OUT) :: y(n)
END SUBROUTINE aprod
END INTERFACE
length=npar
length=npar*ncon
CALL mpalloc(Av,length,'qlpssq: Av')
mbnd=max(0,m-1) ! band width without diagonal
! A*V
ioff1=0
DO i=1,ncon
CALL aprod(npar,matV(ioff1+1:ioff1+npar),Av(ioff1+1:ioff1+npar))
ioff1=ioff1+npar
END DO
DO j=1,ncon
k=j
IF (t) k=ncon+1-j
kn=npar+k-ncon
! column offset
ioff1=(k-1)*npar
! transformation (diagonal block)
! diagonal block
! v^t*A*v
vtAv=dot_product(matV(ioff1+1:ioff1+kn),Av(ioff1+1:ioff1+kn))
! update
ioff2=0
DO i=1,kn
! correct with 2*(2v*vtAv*v^t - Av*v^t - (Av*v^t)^t)
DO l=max(1,i-mbnd),i
ioff2=ioff2+1
B(ioff2)=B(ioff2)+2.0_mpd*((2.0_mpd*matV(ioff1+i)*vtAv-Av(ioff1+i))*matV(ioff1+l)-Av(ioff1+l)*matV(ioff1+i))
END DO
END DO
! off diagonal block
DO i=kn+1,npar
! correct with -2Av*v^t
DO l=max(1,i-mbnd),i
ioff2=ioff2+1
B(ioff2)=B(ioff2)-2.0_mpd*Av(ioff1+i)*matV(ioff1+l)
END DO
END DO
! correct A*v for the remainung v
DO j2=j+1,ncon
k2=j2
IF (t) k2=ncon+1-j2
ioff2=(k2-1)*npar
vtvp=dot_product(matV(ioff1+1:ioff1+npar),matV(ioff2+1:ioff2+npar)) ! v^t*v'
vtAvp=dot_product(matV(ioff1+1:ioff1+npar),Av(ioff2+1:ioff2+npar)) ! v^t*(A*v')
DO i=1,kn
Av(ioff2+i)=Av(ioff2+i)+2.0_mpd*((2.0_mpd*matV(ioff1+i)*vtAv-Av(ioff1+i))*vtvp-matV(ioff1+i)*vtAvp)
END DO
DO i=kn+1,npar
Av(ioff2+i)=Av(ioff2+i)-2.0_mpd*Av(ioff1+i)*vtvp
END DO
END DO
END DO
CALL mpdealloc(Av)
END SUBROUTINE qlpssq
!> Get eigenvalues.
!!
!! Get smallest and largest eingenvalue of L.
!!
!! \param [out] emin smallest eigenvalue
!! \param [out] emax largest eigenvalue
!!
SUBROUTINE qlgete(emin,emax)
USE mpqldec
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpl) :: idiag
REAL(mpd), INTENT(OUT) :: emin
REAL(mpd), INTENT(OUT) :: emax
idiag=1
emax=matL(1)
emin=emax
DO i=2,ncon
idiag=idiag+ncon+1
IF (ABS(emax) < ABS(matL(idiag))) emax=matL(idiag)
IF (ABS(emin) > ABS(matL(idiag))) emin=matL(idiag)
END DO
END SUBROUTINE qlgete
!> Backward substitution.
!!
!! Get y from L^t*y=d.
!!
!! \param [in] d Ncon vector, resdiduals
!! \param [out] y Ncon vector, solution
!!
SUBROUTINE qlbsub(d,y)
USE mpqldec
IMPLICIT NONE
INTEGER(mpl) :: idiag
INTEGER(mpi) :: k
REAL(mpd), INTENT(IN) :: d(ncon)
REAL(mpd), INTENT(OUT) :: y(ncon)
! solve L*y=d by forward substitution
idiag=ncon*ncon
DO k=ncon,1,-1
y(k)=(d(k)-dot_product(matL(idiag+1:idiag+ncon-k),y(k+1:ncon)))/matL(idiag)
idiag=idiag-ncon-1
END DO
END SUBROUTINE qlbsub
pede.f90
View file @ 4c587a8c
This diff is collapsed.
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

imprint | privacy | Conditions of Use