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

Checking for lower/right border for local fit too

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@163 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent 513410e0
......@@ -72,7 +72,7 @@ MODULE mpmod
INTEGER(mpi) :: matrit=0 !< matrix calculation up to iteration MATRIT
INTEGER(mpi) :: icalcm=0 !< calculation mode (for \ref xloopn "XLOOPN") , >0: calculate matrix
INTEGER(mpi) :: numbit=1 !< number of bits for pair counters
INTEGER(mpi) :: nbndr =0 !< number of records with bordered band matrix for local fit
INTEGER(mpi), DIMENSION(2) :: nbndr =0 !< number of records with bordered band matrix for local fit (upper/left, lower/right)
INTEGER(mpi) :: nbdrx =0 !< max border size for local fit
INTEGER(mpi) :: nbndx =0 !< max band width for local fit
INTEGER(mpi) :: nrecer=0 !< record with error (rank deficit or Not-a-Number) for printout
......
......@@ -44,7 +44,8 @@
!! VABDEC
!!
!! Solution by Cholesky decomposition of bordered band matrix
!! SQMIBB (CHK)
!! SQMIBB upper/left border (CHK)
!! SQMIBB2 lower/right border (CHK)
!!
!! Matrix/vector products
!! DBDOT dot vector product
......@@ -2473,3 +2474,262 @@ SUBROUTINE sqmibb(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag
END IF
END SUBROUTINE sqmibb
! 181105 C. Kleinwort, DESY-BELLE
!> Band bordered matrix.
!!
!! Obtain solution of a system of linear equations with symmetric
!! band bordered matrix (V * X = B), on request inverse is calculated.
!! For band part root-free Cholesky decomposition and forward/backward
!! substitution is used.
!!
!! Use decomposition in band and border part for block matrix algebra:
!!
!! | A Ct | | x1 | | b1 | , A is the band part
!! | | * | | = | | , Ct is the mixed part
!! | C D | | x2 | | b2 | , D is the border part
!!
!! \param [in,out] v symmetric N-by-N matrix in symmetric storage mode
!! (V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, ...),
!! replaced by inverse matrix
!! \param [in,out] b N-vector, replaced by solution vector
!! \param [in] n size of V, B
!! \param [in] nbdr border size
!! \param [in] nbnd band width
!! \param [in] inv =1 calculate band part of inverse (for pulls),
!! >1 calculate complete inverse
!! \param [out] nrank rank of matrix V
!! \param [out] vbnd band part of V
!! \param [out] vbdr border part of V
!! \param [out] aux solutions for border rows
!! \param [out] vbk matrix for border solution
!! \param [out] vzru border solution
!! \param [out] scdiag workspace (D)
!! \param [out] scflag workspace (I)
!!
SUBROUTINE sqmibb2(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
USE mpdef
! REAL(mpd) scratch arrays:
! VBND(N*(NBND+1)) = storage of band part
! VBDR(N* NBDR) = storage of border part
! AUX (N* NBDR) = intermediate results
! cost[dot ops] ~= (N-NBDR)*(NBDR+NBND+1)**2 + NBDR**3/3 (leading term, solution only)
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpi) :: ib
INTEGER(mpi) :: ij
INTEGER(mpi) :: ioff
INTEGER(mpi) :: ip
INTEGER(mpi) :: ip1
INTEGER(mpi) :: is
INTEGER(mpi) :: j
INTEGER(mpi) :: j0
INTEGER(mpi) :: jb
INTEGER(mpi) :: joff
INTEGER(mpi) :: koff
INTEGER(mpi) :: mp1
INTEGER(mpi) :: nb1
INTEGER(mpi) :: nmb
INTEGER(mpi) :: npri
INTEGER(mpi) :: nrankb
REAL(mpd), INTENT(IN OUT) :: v(*)
REAL(mpd), INTENT(OUT) :: b(n)
INTEGER(mpi), INTENT(IN) :: n
INTEGER(mpi), INTENT(IN) :: nbdr
INTEGER(mpi), INTENT(IN) :: nbnd
INTEGER(mpi), INTENT(IN) :: inv
INTEGER(mpi), INTENT(OUT) :: nrank
REAL(mpd), INTENT(OUT) :: vbnd(n*(nbnd+1))
REAL(mpd), INTENT(OUT) :: vbdr(n*nbdr)
REAL(mpd), INTENT(OUT) :: aux(n*nbdr)
REAL(mpd), INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2)
REAL(mpd), INTENT(OUT) :: vzru(nbdr)
REAL(mpd), INTENT(OUT) :: scdiag(nbdr)
INTEGER(mpi), INTENT(OUT) :: scflag(nbdr)
SAVE npri
DATA npri / 100 /
! ...
nrank=0
mp1=nbnd+1
nmb=n-nbdr
nb1=nmb+1
! copy band part
DO i=1,nmb
ip=(i*(i+1))/2
is=0
DO j=i,MIN(nmb,i+nbnd)
ip=ip+is
is=j
ib=j-i+1
vbnd(ib+(i-1)*mp1)=v(ip)
END DO
END DO
! copy border part
IF (nbdr > 0) THEN
ioff=0
DO i=nb1,n
ip=(i*(i-1))/2
DO j=1,i
vbdr(ioff+j)=v(ip+j)
END DO
ioff=ioff+n
END DO
END IF
CALL dbcdec(vbnd,mp1,nmb,aux)
! use? CALL DBFDEC(VBND,MP1,NMB) ! modified decomp., numerically more stable
! CALL DBCPRB(VBND,MP1,NMB)
ip=1
DO i=1, nmb
IF (vbnd(ip) <= 0.0_mpd) THEN
npri=npri-1
IF (npri >= 0) THEN
IF (vbnd(ip) == 0.0_mpd) THEN
PRINT *, ' SQMIBB2 matrix singular', n, nbdr, nbnd
ELSE
PRINT *, ' SQMIBB2 matrix not positive definite', n, nbdr, nbnd
END IF
END IF
! return zeros
DO ip=1,n
b(ip)=0.0_mpd
END DO
DO ip=1,(n*n+n)/2
v(ip)=0.0_mpd
END DO
RETURN
END IF
ip=ip+mp1
END DO
nrank=nmb
IF (nbdr == 0) THEN ! special case NBDR=0
CALL dbcslv(vbnd,mp1,nmb,b,b)
IF (inv > 0) THEN
IF (inv > 1) THEN
CALL dbcinv(vbnd,mp1,nmb,v)
ELSE
CALL dbcinb(vbnd,mp1,nmb,v)
END IF
END IF
ELSE ! general case NBDR>0
ioff=0
DO ib=1,nbdr
! solve for aux. vectors
CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff+1),aux(ioff+1))
! zT ru
vzru(ib)=b(nmb+ib)
DO i=1,nmb
vzru(ib)=vzru(ib)-b(i)*aux(ioff+i)
END DO
ioff=ioff+n
END DO
! solve for band part only
CALL dbcslv(vbnd,mp1,nmb,b,b)
! Ck - cT z
ip=0
ioff=0
koff=nmb
DO ib=1,nbdr
joff=0
DO jb=1,ib
ip=ip+1
vbk(ip)=vbdr(koff+jb)
DO i=1,nmb
vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i)
END DO
joff=joff+n
END DO
ioff=ioff+n
koff=koff+n
END DO
! solve border part
CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag)
IF (nrankb == nbdr) THEN
nrank=nrank+nbdr
ELSE
npri=npri-1
IF (npri >= 0) PRINT *, ' SQMIBB2 undef border ', n, nbdr, nbnd, nrankb
DO ib=1,nbdr
vzru(ib)=0.0_mpd
END DO
DO ip=(nbdr*nbdr+nbdr)/2,1,-1
vbk(ip)=0.0_mpd
END DO
END IF
! smoothed data points
ioff=0
DO ib=1, nbdr
DO i=1,nmb
b(i)=b(i)-vzru(ib)*aux(ioff+i)
END DO
ioff=ioff+n
b(nmb+ib)=vzru(ib)
END DO
! inverse requested ?
IF (inv > 0) THEN
IF (inv > 1) THEN
CALL dbcinv(vbnd,mp1,nmb,v)
ELSE
CALL dbcinb(vbnd,mp1,nmb,v)
END IF
! assemble band and border
IF (nbdr > 0) THEN
! band part
ip1=(nmb*nmb+nmb)/2
DO i=nmb-1,0,-1
j0=0
IF (inv == 1) j0=MAX(0,i-nbnd)
DO j=i,j0,-1
ioff=1
DO ib=1,nbdr
joff=1
DO jb=1,nbdr
ij=MAX(ib,jb)
ij=(ij*ij-ij)/2+MIN(ib,jb)
v(ip1)=v(ip1)+vbk(ij)*aux(ioff+i)*aux(joff+j)
joff=joff+n
END DO
ioff=ioff+n
END DO
ip1=ip1-1
END DO
ip1=ip1-j0
END DO
! band part
ip1=(nmb*nmb+nmb)/2
ip=0
DO ib=1,nbdr
DO i=1,nmb
ip1=ip1+1
v(ip1)=0.0_mpd
joff=0
DO jb=1,nbdr
ij=MAX(ib,jb)
ij=(ij*ij-ij)/2+MIN(ib,jb)
v(ip1)=v(ip1)-vbk(ij)*aux(i+joff)
joff=joff+n
END DO
END DO
DO jb=1,ib
ip1=ip1+1
ip=ip+1
v(ip1)=vbk(ip)
END DO
END DO
END IF
END IF
END IF
END SUBROUTINE sqmibb2
......@@ -2467,11 +2467,12 @@ SUBROUTINE loopn
END IF
! local fit: band matrix structure !?
IF (nloopn == 1.AND.nbndr > 0) THEN
IF (nloopn == 1.AND.nbndr(1)+nbndr(2) > 0) THEN
DO lun=6,8,2
WRITE(lun,*) ' '
WRITE(lun,*) ' === local fits have bordered band matrix structure ==='
WRITE(lun,101) ' NBNDR',nbndr,'number of records'
IF (nbndr(1) > 0 ) WRITE(lun,101) ' NBNDR',nbndr(1),'number of records (upper/left border)'
IF (nbndr(2) > 0 ) WRITE(lun,101) ' NBNDR',nbndr(2),'number of records (lower/right border)'
WRITE(lun,101) ' NBDRX',nbdrx,'max border size'
WRITE(lun,101) ' NBNDX',nbndx,'max band width'
END DO
......@@ -2893,6 +2894,7 @@ SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff)
INTEGER(mpi) :: kx
INTEGER(mpi) :: mbdr
INTEGER(mpi) :: mbnd
INTEGER(mpi) :: mside
INTEGER(mpi) :: nalc
INTEGER(mpi) :: nalg
INTEGER(mpi) :: nan
......@@ -3130,7 +3132,8 @@ SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff)
! check matrix for bordered band structure (MBDR+MBND+1 <= NALC)
mbnd=-1
mbdr=nalc
DO i=1, nalc
mside=-1 ! side (1: upper/left border, 2: lower/right border)
DO i=1, 2*nalc
ibandh(i)=0
END DO
irow=1
......@@ -3211,8 +3214,10 @@ SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff)
! check for band matrix substructure
IF (iter == 1) THEN
id=IABS(ij-ik)+1
im=MIN(ij,ik)
im=MIN(ij,ik) ! upper/left border
ibandh(id)=MAX(ibandh(id),im)
im=MIN(nalc+1-ij,nalc+1-ik) ! lower/rght border (mirrored)
ibandh(nalc+id)=MAX(ibandh(nalc+id),im)
END IF
END DO
END DO
......@@ -3220,9 +3225,10 @@ SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff)
! for non trivial fits check for bordered band matrix structure
IF (iter == 1.AND.nalc > 5.AND.lfitbb > 0) THEN
kx=-1
kbdr=0
kbdrx=0
icmn=INT(nalc,mpl)**3 ! cost (*6) should improve by at least factor 2
! upper/left border ?
kbdr=0
DO k=nalc,2,-1
kbnd=k-2
kbdr=MAX(kbdr,ibandh(k))
......@@ -3231,8 +3237,24 @@ SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff)
icmn=icost
kx=k
kbdrx=kbdr
mside=1
END IF
END DO
END DO
IF (kx < 0) THEN
! lower/right border instead?
kbdr=0
DO k=nalc,2,-1
kbnd=k-2
kbdr=MAX(kbdr,ibandh(k+nalc))
icost=6*INT(nalc-kbdr,mpl)*INT(kbnd+kbdr+1,mpl)**2+2*INT(kbdr,mpl)**3
IF (icost < icmn) THEN
icmn=icost
kx=k
kbdrx=kbdr
mside=2
END IF
END DO
END IF
IF (kx > 0) THEN
mbnd=kx-2
mbdr=kbdrx
......@@ -3242,7 +3264,7 @@ SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff)
IF (mbnd >= 0) THEN
! fast solution for border banded matrix (inverse for ICALCM>0)
IF (nloopn == 1) THEN
nbndr=nbndr+1
nbndr(mside)=nbndr(mside)+1
nbdrx=MAX(nbdrx,mbdr)
nbndx=MAX(nbndx,mbnd)
END IF
......@@ -3250,8 +3272,13 @@ SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff)
inv=0
IF (nloopn <= lfitnp.AND.iter == 1) inv=1 ! band part of inverse (for pulls)
IF (icalcm == 1.OR.lprnt) inv=2 ! complete inverse
CALL sqmibb(clmat,blvec,nalc,mbdr,mbnd,inv,nrank, &
vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
IF (mside == 1) THEN
CALL sqmibb(clmat,blvec,nalc,mbdr,mbnd,inv,nrank, &
vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
ELSE
CALL sqmibb2(clmat,blvec,nalc,mbdr,mbnd,inv,nrank, &
vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
ENDIF
ELSE
! full inversion and solution
inv=2
......@@ -5753,7 +5780,7 @@ SUBROUTINE vmprep(msize)
CALL mpalloc(vzru,length,' local fit scratch array: vzru')
CALL mpalloc(scdiag,length,' local fit scratch array: scdiag')
CALL mpalloc(scflag,length,' local fit scratch array: scflag')
CALL mpalloc(ibandh,length,' local fit band width hist.: ibandh')
CALL mpalloc(ibandh,2*length,' local fit band width hist.: ibandh')
CALL mpalloc(globalMatD,msize(1),'global matrix (D)' )
CALL mpalloc(globalMatF,msize(2),'global matrix (F)')
......
Markdown is supported
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