Skip to content

GitLab

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

Implementation of parameter groups, sparse similarity operations 

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@187 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent 17e0d7fc
Expand all Hide whitespace changes
Inline Side-by-side
mpbits.f90
View file @ f263758b
This diff is collapsed.
mpmod.f90
View file @ f263758b
......@@ -4,7 +4,7 @@
!! \author Claus Kleinwort, DESY, 2012 (Claus.Kleinwort@desy.de)
!!
!! \copyright
!! Copyright (c) 2012 - 2019 Deutsches Elektronen-Synchroton,
!! Copyright (c) 2012 - 2020 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
......@@ -79,7 +79,6 @@ MODULE mpmod
INTEGER(mpi) :: mreqpe=1 !< min number of pair entries
INTEGER(mpi) :: mhispe=0 !< upper bound for pair entry histogrammimg
INTEGER(mpi) :: msngpe=-1 !< upper bound for pair entry single precision storage
INTEGER(mpi) :: mcmprs=0 !< compression flag for sparsity (column indices)
INTEGER(mpi) :: mextnd=0 !< flag for extended storage (both 'halves' of sym. mat. for improved access patterns)
INTEGER(mpi) :: mthrd =1 !< number of (OpenMP) threads
INTEGER(mpi) :: mxrec =0 !< max number of records
......@@ -117,10 +116,13 @@ MODULE mpmod
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 all parameters (global par. + Lagrange mult.)
INTEGER(mpi) :: nagb !< number of all parameters (var. global par. + Lagrange mult.)
INTEGER(mpi) :: nfgb !< number of fit parameters
INTEGER(mpi) :: ncgb !< number of constraints
INTEGER(mpi) :: ncgbe !< number of empty constraints (no variable parameters)
INTEGER(mpi) :: ntpgrp !< number of parameter groups
INTEGER(mpi) :: nvpgrp !< number of variable parameter groups
INTEGER(mpi) :: napgrp !< number of all parameter groups (variable + Lagrange mult.)
INTEGER(mpi) :: npblck !< number of (disjoint) parameter blocks
INTEGER(mpi) :: ncblck !< number of (disjoint) constraint blocks
INTEGER(mpi) :: mszcon !< (integrated block) matrix size for constraint matrix
......@@ -218,7 +220,8 @@ MODULE mpmod
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: globalParLabelIndex !< global parameters label, total -> var. index
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: globalParHashTable !< global parameters hash table
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: globalParVarToTotal !< global parameters variable -> total index
INTEGER(mpi), DIMENSION(-7:0) :: globalParHeader = 0 !< global parameters (mapping) header
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: globalAllParToGroup !< all parameters variable -> group index
INTEGER(mpi), DIMENSION(-8:0) :: globalParHeader = 0 !< global parameters (mapping) header
!!
!! 0: length of labels/indices; \n
!! -1: number of stored items; \n
......@@ -228,9 +231,11 @@ MODULE mpmod
!! -5: number of overflows; \n
!! -6: nr of variable parameters; \n
!! -7: call counter for build-up;
!! -8: number of sorted items;
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: globalTotIndexGroups !< 1. (total) index, size per group
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: globalAllIndexGroups !< 1. (all variable) index per group
! row information for sparse matrix
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: sparseMatrixCompression !< compression info (per row)
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: sparseMatrixColumns !< (compressed) list of columns for sparse matrix
INTEGER(mpl), DIMENSION(:,:), ALLOCATABLE :: sparseMatrixOffsets !< row offsets for column list, sparse matrix elements
! read buffer
......@@ -260,8 +265,11 @@ MODULE mpmod
REAL(mpd), DIMENSION(:), ALLOCATABLE::vzru !< local fit 'border solution'
REAL(mpd), DIMENSION(:), ALLOCATABLE::scdiag !< local fit workspace (D)
INTEGER(mpi), DIMENSION(:), ALLOCATABLE:: scflag !< local fit workspace (I)
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: localEquations !< indices (ISJAJB) for local equations (measurements)
REAL(mpd), DIMENSION(:), ALLOCATABLE :: localCorrections !< local fit corrections (to residuals)
REAL(mpd), DIMENSION(:), ALLOCATABLE :: localGlobalMatrix !< matrix correlating local and global par
REAL(mpd), DIMENSION(:), ALLOCATABLE :: localGlobalMatrix !< matrix correlating local and global par, content
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: localGlobalMap !< matrix correlating local and global par, map (counts)
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: localGlobalStructure !< matrix correlating local and global par, (sparsity) structure
! update of global matrix
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: writeBufferInfo !< write buffer management (per thread)
REAL(mps), DIMENSION(:,:), ALLOCATABLE :: writeBufferData !< write buffer data (largest residual, Chi2/ndf, per thread)
......
mpnum.f90
View file @ f263758b
......@@ -9,7 +9,7 @@
!! \author Claus Kleinwort, DESY (maintenance and developement)
!!
!! \copyright
!! Copyright (c) 2009 - 2019 Deutsches Elektronen-Synchroton,
!! Copyright (c) 2009 - 2020 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
......@@ -54,6 +54,7 @@
!! DBSVX LARGE symmetric matrix vector (CHK)
!! DBGAX general matrix vector
!! DBAVAT AVAT product
!! DBAVATS AVAT product for sparse A (CHK)
!! DBMPRV print parameter and matrix
!! DBPRV print matrix (CHK)
!!
......@@ -68,6 +69,7 @@
!! SORT1K sort 1-dim key-array (CHK)
!! SORT2K sort 2-dim key-array
!! SORT2I sort 2-dim key-array with index (CHK)
!! SORT22 sort 2-dim key-array with two additional values (CHK)
!!
!----------------------------------------------------------------------
......@@ -1335,6 +1337,89 @@ SUBROUTINE dbavat(v,a,w,n,ms)
END DO ! end do I
END SUBROUTINE dbavat
!> A V AT product (similarity, sparse).
!!
!! Multiply symmetric N-by-N matrix from the left with sparse M-by-N
!! matrix and from the right with the transposed of the same general
!! matrix to form symmetric M-by-M matrix (used for error propagation).
!!
!! \param [in] V symmetric N-by-N matrix
!! \param [in] A sparse M-by-N matrix, content
!! \param [in] IS sparse M-by-N matrix, structure
!! \param [in,out] W symmetric M-by-M matrix
!! \param [in] MS rows of A (-rows: don't reset W)
!! \param [in] N columns of A
!! \param [in] SC scratch array
!!
!! Sparsity structure:
!! - IS(1..M): row offsets
!! - IS(M+1..N+M+1): column offsets
!! - IS(IS(1)+1..IS(M+1)): non-zero columns (column number, index for A)
!! - IS(IS(M+1)+1..IS(M+N+1)): non-zero rows (row number, index for A)
!!
SUBROUTINE dbavats(v,a,is,w,n,ms,sc)
USE mpdef
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpi) :: ic
INTEGER(mpi) :: ij
INTEGER(mpi) :: ijs
INTEGER(mpi) :: in
INTEGER(mpi) :: ir
INTEGER(mpi) :: j
INTEGER(mpi) :: k
INTEGER(mpi) :: l
INTEGER(mpi) :: lk
INTEGER(mpi) :: m
REAL(mpd), INTENT(IN) :: v(*)
REAL(mpd), INTENT(IN) :: a(*)
INTEGER(mpi), INTENT(IN) :: is(*)
REAL(mpd), INTENT(INOUT) :: w(*)
INTEGER(mpi), INTENT(IN) :: n
INTEGER(mpi), INTENT(IN) :: ms
INTEGER(mpi), INTENT(OUT) :: sc(*)
REAL(mpd) :: cik
! ...
m=ms
IF (m > 0) THEN
DO i=1,(m*m+m)/2
w(i)=0.0_mpd ! reset output matrix
END DO
ELSE
m=-m
END IF
! offsets in V
sc(1)=0
DO k=2,n
sc(k)=sc(k-1)+k-1
END DO
ijs=0
DO i=1,m
ijs=ijs+i-1
DO k=1,n
cik=0.0_mpd
DO l=is(i)+1,is(i+1),2
ic=is(l)
in=is(l+1)
lk=sc(max(k,ic))+min(k,ic)
cik=cik+a(in)*v(lk)
END DO
DO j=is(m+k)+1,is(m+k+1),2
ir=is(j)
in=is(j+1)
IF (ir > i) EXIT
ij=ijs+ir
w(ij)=w(ij)+cik*a(in)
END DO
END DO
END DO
END SUBROUTINE dbavats
!> Print symmetric matrix, vector.
!!
!! Prints the n-vector X and the symmetric N-by-N covariance matrix
......@@ -1713,7 +1798,7 @@ SUBROUTINE sort2i(a,n)
INTEGER(mpi) ::maxlev
INTEGER(mpi) ::a1 ! pivot key
INTEGER(mpi) ::a2 ! pivot key
INTEGER(mpi) ::at ! pivot key
INTEGER(mpi) ::at(3)
INTEGER(mpi), INTENT(IN OUT) :: a(3,*)
INTEGER(mpi), INTENT(IN) :: n
......@@ -1724,15 +1809,9 @@ SUBROUTINE sort2i(a,n)
r=n
10 IF(r-l == 1) THEN ! sort two elements L and R
IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
at=a(1,l) ! exchange L <-> R
a(1,l)=a(1,r)
a(1,r)=at
at=a(2,l)
a(2,l)=a(2,r)
a(2,r)=at
at=a(3,l)
a(3,l)=a(3,r)
a(3,r)=at
at=a(:,l) ! exchange L <-> R
a(:,l)=a(:,r)
a(:,r)=at
END IF
r=l
END IF
......@@ -1759,15 +1838,9 @@ SUBROUTINE sort2i(a,n)
IF(a(1,j) > a1) GO TO 30
IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
IF(i <= j) THEN
at=a(1,i) ! exchange I <-> J
a(1,i)=a(1,j)
a(1,j)=at
at=a(2,i)
a(2,i)=a(2,j)
a(2,j)=at
at=a(3,i)
a(3,i)=a(3,j)
a(3,j)=at
at=a(:,i) ! exchange I <-> J
a(:,i)=a(:,j)
a(:,j)=at
GO TO 20
END IF
IF(lev+2 > nlev) THEN
......@@ -1789,6 +1862,94 @@ SUBROUTINE sort2i(a,n)
GO TO 10
END SUBROUTINE sort2i
!> Quick sort 2 with index.
!!
!! Quick sort of A(4,N) integer.
!!
!! \param[in,out] a vector (pair) of integers, sorted at return and an index
!! \param[in] n size of vector
SUBROUTINE sort22(a,n)
USE mpdef
IMPLICIT NONE
INTEGER(mpi) :: nlev ! stack size
PARAMETER (nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
INTEGER(mpi) :: i
INTEGER(mpi) ::j
INTEGER(mpi) ::l
INTEGER(mpi) ::r
INTEGER(mpi) ::lev
INTEGER(mpi) ::lr(nlev)
INTEGER(mpi) ::lrh
INTEGER(mpi) ::maxlev
INTEGER(mpi) ::a1 ! pivot key
INTEGER(mpi) ::a2 ! pivot key
INTEGER(mpi) ::at(4)
INTEGER(mpi), INTENT(IN OUT) :: a(4,*)
INTEGER(mpi), INTENT(IN) :: n
! ...
maxlev=0
lev=0
l=1
r=n
IF (n<=1) RETURN
10 IF(r-l == 1) THEN ! sort two elements L and R
IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
at=a(:,l) ! exchange L <-> R
a(:,l)=a(:,r)
a(:,r)=at
END IF
r=l
END IF
IF(r == l) THEN
IF(lev <= 0) THEN
WRITE(*,*) 'SORT22 (quicksort): maxlevel used/available =', maxlev,'/64'
RETURN
END IF
lev=lev-2
l=lr(lev+1)
r=lr(lev+2)
ELSE
! LRH=(L+R)/2
lrh=(l/2)+(r/2) ! avoid bit overflow
IF(MOD(l,2) == 1.AND.MOD(r,2) == 1) lrh=lrh+1
a1=a(1,lrh) ! middle
a2=a(2,lrh)
i=l-1 ! find limits [J,I] with [L,R]
j=r+1
20 i=i+1
IF(a(1,i) < a1) GO TO 20
IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20
30 j=j-1
IF(a(1,j) > a1) GO TO 30
IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
IF(i <= j) THEN
at=a(:,i) ! exchange I <-> J
a(:,i)=a(:,j)
a(:,j)=at
GO TO 20
END IF
IF(lev+2 > nlev) THEN
CALL peend(33,'Aborted, stack overflow in quicksort')
STOP 'SORT22 (quicksort): stack overflow'
END IF
IF(r-i < j-l) THEN
lr(lev+1)=l
lr(lev+2)=j
l=i
ELSE
lr(lev+1)=i
lr(lev+2)=r
r=j
END IF
lev=lev+2
maxlev=MAX(maxlev,lev)
END IF
GO TO 10
END SUBROUTINE sort22
!> Chi2/ndf cuts.
!!
!! Return limit in Chi^2/ndf for N sigmas (N=1, 2 or 3).
......
mpqldec.f90
View file @ f263758b
......@@ -4,7 +4,7 @@
!! \author Claus Kleinwort, DESY, 2015 (Claus.Kleinwort@desy.de)
!!
!! \copyright
!! Copyright (c) 2015-2019 Deutsches Elektronen-Synchroton,
!! Copyright (c) 2015-2020 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
......@@ -32,6 +32,7 @@ MODULE mpqldec
INTEGER(mpi) :: ncon !< number of constraints
INTEGER(mpi) :: nblock !< number of blocks
INTEGER(mpi) :: iblock !< active block
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: sparseV !< sparsity structure matV
REAL(mpd), DIMENSION(:), ALLOCATABLE :: matV !< unit normals (v_i) of Householder reflectors
REAL(mpd), DIMENSION(:), ALLOCATABLE :: matL !< lower diagonal matrix L
REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecN !< normal vector
......@@ -62,6 +63,8 @@ SUBROUTINE qlini(n,m,l)
nblock=l
iblock=1
! allocate
length=5*ncon
CALL mpalloc(sparseV,length,'QLDEC: sparsity structure of V')
length=npar*ncon
CALL mpalloc(matV,length,'QLDEC: V')
length=ncon*ncon
......@@ -153,6 +156,11 @@ SUBROUTINE qldec(a)
! store normal vector
matV(ioff1+1:ioff1+kn)=vecN(1:kn)
matV(ioff1+kn+1:ioff1+npar)=0.0_mpd
! sparsity structure
ioff3=(k-1)*5
sparseV(ioff3+1)=1 ! number of non-zero regions
sparseV(ioff3+2)=1 ! start
sparseV(ioff3+3)=kn ! end
END DO
END SUBROUTINE qldec
......@@ -300,6 +308,13 @@ SUBROUTINE qldecb(a,bpar,bcon)
matV(ioff1+1:ioff1+npar)=0.0_mpd
matV(ioff1+ifirst-ipoff:ioff1+ilast-ipoff)=vecN(ifirst:ilast)
matV(ioff1+kn-ipoff)=vecN(kn)
! sparsity structure
ioff3=(k-1)*5
sparseV(ioff3+1)=2 ! number of non-zero regions
sparseV(ioff3+2)=ifirst-ipoff ! start 1
sparseV(ioff3+3)=ilast-ipoff ! end 1
sparseV(ioff3+4)=kn-ipoff ! start 2
sparseV(ioff3+5)=kn-ipoff ! end 2
END DO
END DO
......@@ -495,7 +510,6 @@ SUBROUTINE qlssq(aprod,A,t)
REAL(mpd), INTENT(OUT) :: y(n)
END SUBROUTINE aprod
END INTERFACE
length=npar
CALL mpalloc(Av,length,'qlssq: A*v')
......@@ -582,14 +596,16 @@ SUBROUTINE qlpssq(aprod,B,m,t)
LOGICAL, INTENT(IN) :: t
INTERFACE
SUBROUTINE aprod(n,l,x,y) ! y=A*x
SUBROUTINE aprod(n,l,x,y,is) ! y=A*x
USE mpdef
INTEGER(mpi), INTENT(in) :: n
INTEGER(mpl), INTENT(in) :: l
REAL(mpd), INTENT(IN) :: x(n)
REAL(mpd), INTENT(OUT) :: y(n)
INTEGER(mpi), INTENT(in) :: is(*)
END SUBROUTINE aprod
END INTERFACE
!$POMP INST BEGIN(qlpssq)
length=npar
length=npar*ncon
......@@ -599,7 +615,7 @@ SUBROUTINE qlpssq(aprod,B,m,t)
! A*V
ioff1=0
DO i=1,ncon
CALL aprod(npar,0_mpl,matV(ioff1+1:ioff1+npar),Av(ioff1+1:ioff1+npar))
CALL aprod(npar,0_mpl,matV(ioff1+1:ioff1+npar),Av(ioff1+1:ioff1+npar),sparseV(i*5-4))
ioff1=ioff1+npar
END DO
......@@ -648,6 +664,7 @@ SUBROUTINE qlpssq(aprod,B,m,t)
END DO
CALL mpdealloc(Av)
!$POMP INST END(qlpssq)
END SUBROUTINE qlpssq
......
pede.f90
View file @ f263758b
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