Commit 22a9e135 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

inversion for large problems: little speedup, fix

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@195 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent e81131cc
......@@ -200,6 +200,7 @@ MODULE mpmod
! 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 :: workspaceRow !< (pivot) row 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
......
......@@ -222,8 +222,9 @@ END SUBROUTINE sqminv
!! \param [out] NRANK rank of matrix V
!! \param [out] DIAG double precision scratch array
!! \param [out] NEXT integer aux array
!! \param [out] VK double precision scratch array (pivot)
SUBROUTINE sqminl(v,b,n,nrank,diag,next) !
SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk) !
USE mpdef
IMPLICIT NONE
......@@ -236,10 +237,12 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next) !
REAL(mpd), INTENT(IN OUT) :: v(*)
REAL(mpd), INTENT(OUT) :: b(n)
INTEGER(mpi), INTENT(IN) :: n
INTEGER(mpi), INTENT(OUT) :: nrank
INTEGER(mpi), INTENT(IN) :: n
INTEGER(mpi), INTENT(OUT) :: nrank
REAL(mpd), INTENT(OUT) :: diag(n)
INTEGER(mpi), INTENT(OUT) :: next(n)
INTEGER(mpi), INTENT(OUT) :: next(n)
REAL(mpd), INTENT(OUT) :: vk(n)
INTEGER(mpl) :: i8
INTEGER(mpl) :: j8
INTEGER(mpl) :: jj
......@@ -248,8 +251,6 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next) !
INTEGER(mpl) :: kkmk
INTEGER(mpl) :: jk
INTEGER(mpl) :: jl
INTEGER(mpl) :: llk
INTEGER(mpl) :: ljl
REAL(mpd) :: vkk
REAL(mpd) :: vjk
......@@ -303,6 +304,7 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next) !
DO j=1,n
IF(j == k) THEN
jk=kk
vk(j)=0.
ELSE
IF(j < k) THEN
jk=jk+1
......@@ -310,41 +312,23 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next) !
jk=jk+int8(j)-1
END IF
v(jk)=v(jk)*vkk
vk(j)=v(jk)
END IF
END DO
! parallelize row loop
! slot of 128 'J' for next idle thread
!$OMP PARALLEL DO &
!$OMP PRIVATE(JL,JK,L,LJL,LLK,VJK,J8) &
!$OMP PRIVATE(JL,VJK,J8) &
!$OMP SCHEDULE(DYNAMIC,128)
DO j=n,1,-1
IF(j == k) CYCLE
j8=int8(j)
jl=j8*(j8-1)/2
IF(j /= k) THEN
IF(j < k) THEN
jk=kkmk+j8
ELSE
jk=k8+jl
END IF
vjk =v(jk)/vkk
b(j) =b(j)-b(k)*vjk
ljl=jl
llk=kkmk
DO l=1,MIN(j,k-1)
ljl=ljl+1
llk=llk+1
v(ljl)=v(ljl)-v(llk)*vjk
END DO
ljl=ljl+1
llk=kk
DO l=k+1,j
ljl=ljl+1
llk=llk+l-1
v(ljl)=v(ljl)-v(llk)*vjk
END DO
END IF
vjk =vk(j)/vkk
b(j) =b(j)-b(k)*vjk
v(jl+1:jl+j)=v(jl+1:jl+j)-vk(1:j)*vjk
END DO
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO
ELSE
DO k=1,n
k8=int8(k)
......
......@@ -510,6 +510,8 @@ SUBROUTINE qlssq(aprod,A,t)
REAL(mpd), INTENT(OUT) :: y(n)
END SUBROUTINE aprod
END INTERFACE
!$POMP INST BEGIN(qlssq)
length=npar
CALL mpalloc(Av,length,'qlssq: A*v')
......@@ -553,6 +555,7 @@ SUBROUTINE qlssq(aprod,A,t)
END DO
CALL mpdealloc(Av)
!$POMP INST END(qlssq)
END SUBROUTINE qlssq
......
......@@ -51,7 +51,7 @@
!! 1. Download the software package from the DESY \c svn server to
!! \a target directory, e.g.:
!!
!! svn checkout http://svnsrv.desy.de/public/MillepedeII/tags/V04-07-03 target
!! svn checkout http://svnsrv.desy.de/public/MillepedeII/tags/V04-07-04 target
!!
!! 2. Create **Pede** executable (in \a target directory):
!!
......@@ -7364,6 +7364,7 @@ SUBROUTINE vmprep(msize)
IF(metsol == 1) THEN
CALL mpalloc(workspaceDiag,length,'diagonal of global matrix)') ! double aux 1
CALL mpalloc(workspaceRow,length,'(pivot) row of global matrix)')
! CALL MEGARR('t D',2*NAGB,'auxiliary array') ! double aux 8
END IF
......@@ -7470,7 +7471,7 @@ SUBROUTINE minver
IF(icalcm == 1) THEN
! invert and solve
CALL sqminl(globalMatD(imoff+1:), globalCorrections(ipoff+1:),nfit,nrank, &
workspaceD,workspaceI)
workspaceD,workspaceI,workspaceRow)
IF(nfit /= nrank) THEN
WRITE(*,*) 'Warning: the rank defect of the symmetric',nfit, &
'-by-',nfit,' matrix is ',nfit-nrank,' (should be zero).'
......@@ -8387,6 +8388,27 @@ SUBROUTINE xloopn !
! monitoring of residuals
IF (imonit > 0 .AND. btest(imonit,1)) CALL monres
IF (lunmon > 0) CLOSE(UNIT=lunmon)
IF(metsol <= 2) THEN ! inversion or diagonalization ?
!use elimination for constraints ?
IF(nfgb < nvgb) THEN
! extend, transform matrix
! loop over blocks
DO ib=1,npblck
ipoff=matParBlockOffsets(1,ib)
imoff=vecParBlockOffsets(ib)
icboff=matParBlockOffsets(2,ib) ! constraint block offset
icblst=matParBlockOffsets(2,ib+1) ! constraint block offset
npar=matParBlockOffsets(1,ib+1)-ipoff ! size of block (number of parameters)
ncon=matConsBlocks(1,icblst+1)-matConsBlocks(1,icboff+1) ! number of constraints in (parameter) block
DO i=npar-ncon+1,npar
ioff=((INT8(i)-1)*INT8(i))/2+imoff
globalMatD(ioff+1:ioff+i)=0.0_mpd
END DO
END DO
CALL qlssq(avprd0,globalMatD,.false.) ! Q^t*A*Q
END IF
END IF
dwmean=sumndf/REAL(ndfsum,mpd)
dratio=fvalue/dwmean/REAL(ndfsum-nfgb,mpd)
......@@ -8558,27 +8580,6 @@ SUBROUTINE xloopn !
END IF
IF(metsol <= 2) THEN ! inversion or diagonalization ?
!use elimination for constraints ?
IF(nfgb < nvgb) THEN
! extend, transform matrix
! loop over blocks
DO ib=1,npblck
ipoff=matParBlockOffsets(1,ib)
imoff=vecParBlockOffsets(ib)
icboff=matParBlockOffsets(2,ib) ! constraint block offset
icblst=matParBlockOffsets(2,ib+1) ! constraint block offset
npar=matParBlockOffsets(1,ib+1)-ipoff ! size of block (number of parameters)
ncon=matConsBlocks(1,icblst+1)-matConsBlocks(1,icboff+1) ! number of constraints in (parameter) block
DO i=npar-ncon+1,npar
ioff=((i-1)*i)/2+imoff
globalMatD(ioff+1:ioff+i)=0.0_mpd
END DO
END DO
CALL qlssq(avprd0,globalMatD,.false.) ! Q^t*A*Q
END IF
END IF
CALL prtglo ! print result
IF (warners3) THEN
......
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