Commit 3f65b7d6 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

preconditioning with skyline matrix added

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@142 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent a327d641
......@@ -41,6 +41,7 @@ MODULE mpmod
INTEGER(mpi) :: mitera=1 !< number of iterations
INTEGER(mpi) :: nloopn=0 !< number of data reading, fitting loops
INTEGER(mpi) :: mbandw=0 !< band width of preconditioner matrix
INTEGER(mpi) :: lprecm=0 !< additional flag for preconditioner (band) matrix (>0: preserve rank by skyline matrix)
INTEGER(mpi) :: lunkno=0 !< flag for unkown keywords
INTEGER(mpi) :: lhuber=0 !< Huber down-weighting flag
REAL(mps) :: chicut=0.0 !< cut in terms of 3-sigma cut, first iteration
......
......@@ -1750,6 +1750,7 @@ END FUNCTION chindl
!! decomposition. No fill-in is created ahead in any row or ahead of the
!! first entry in any column, but existing zero-values will become
!! non-zero. The decomposition is done "in-place".
!! (The diagonal will contain the inverse of the diaginal of L).
!!
!! - NRKD = 0 no component removed
!!
......@@ -1763,13 +1764,16 @@ END FUNCTION chindl
!! by about a word length (see line "test for linear dependence"),
!! then the pivot is assumed as zero and the entire row/column is
!! reset to zero, removing the corresponding element from the solution.
!! Optionally use only diagonal element in this case to preserve rank
!! (changing band to skyline matrix).
!!
!! \param [in] n size of matrix
!! \param [in,out] c variable-band matrix, replaced by L
!! \param [in] india pointer array
!! \param [out] nrkd removed components
!! \param [in] iopt >0: use diagonal to preserve rank ('skyline')
SUBROUTINE lltdec(n,c,india,nrkd)
SUBROUTINE lltdec(n,c,india,nrkd,iopt)
USE mpdef
IMPLICIT NONE
......@@ -1785,11 +1789,12 @@ SUBROUTINE lltdec(n,c,india,nrkd)
REAL(mpd), INTENT(IN OUT) :: c(*)
INTEGER(mpi), INTENT(IN) :: india(n)
INTEGER(mpi), INTENT(OUT) :: nrkd
INTEGER(mpi), INTENT(IN) :: iopt
REAL(mpd) eps
! ...
eps = 16.0_mpd * epsilon(eps) ! 16 * precision(mpd)
! ...
! ..
nrkd=0
diag=0.0_mpd
IF(c(india(1)) > 0.0) THEN
......@@ -1821,11 +1826,15 @@ SUBROUTINE lltdec(n,c,india,nrkd)
DO j=mk,k ! reset row K
c(india(k)-k+j)=0.0_mpd
END DO ! J
IF(nrkd == 0) THEN
nrkd=-k
IF (iopt > 0 .and. diag > 0.0) THEN ! skyline
c(india(k))=1.0_mpd/SQRT(diag) ! square root
ELSE
IF(nrkd < 0) nrkd=1
nrkd=nrkd+1
IF(nrkd == 0) THEN
nrkd=-k
ELSE
IF(nrkd < 0) nrkd=1
nrkd=nrkd+1
END IF
END IF
END IF
......@@ -1915,12 +1924,13 @@ END SUBROUTINE lltbwd
!!
!! \param [in] n size of symmetric matrix
!! \param [in] m number of constrains
!! \param [in] ls flag for skyline decomposition
!! \param [in,out] c combined variable-band + constraints matrix, replaced by decomposition
!! \param [in,out] india pointer array
!! \param [out] nrkd removed components
!! \param [out] nrkd2 removed components
!!
SUBROUTINE equdec(n,m,c,india,nrkd,nrkd2)
SUBROUTINE equdec(n,m,ls,c,india,nrkd,nrkd2)
USE mpdef
IMPLICIT NONE
......@@ -1931,6 +1941,7 @@ SUBROUTINE equdec(n,m,c,india,nrkd,nrkd2)
INTEGER(mpi), INTENT(IN) :: n
INTEGER(mpi), INTENT(IN) :: m
INTEGER(mpi), INTENT(IN) :: ls
REAL(mpd), INTENT(IN OUT) :: c(*)
INTEGER(mpi), INTENT(IN OUT) :: india(n+m)
INTEGER(mpi), INTENT(OUT) :: nrkd
......@@ -1938,7 +1949,10 @@ SUBROUTINE equdec(n,m,c,india,nrkd,nrkd2)
! ...
CALL lltdec(n,c,india,nrkd) ! decomposition G G^T
nrkd=0
nrkd2=0
CALL lltdec(n,c,india,nrkd,ls) ! decomposition G G^T
IF (m>0) THEN
DO i=1,m
......@@ -1961,9 +1975,7 @@ SUBROUTINE equdec(n,m,c,india,nrkd,nrkd2)
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
ELSE
nrkd2=0
CALL lltdec(m,c(india(n)+n*m+1),india(n+1),nrkd2,0) ! decomp. H H^T
ENDIF
RETURN
......
......@@ -318,7 +318,8 @@
!!
!! \subsection cmd-bandwidth bandwidth
!! Set band width \ref mpmod::mbandw "mbandw" for
!! \ref minresmodule::minres "MINRES" preconditioner to \a number1.
!! \ref minresmodule::minres "MINRES" preconditioner to \a number1 [0]
!! and additional flag \ref mpmod::lprecm "lprecm" to \number2 [0].
!! \subsection cmd-cache cache
!! Set (read+write) cache size \ref mpmod::ncache "ncache" to \a number1.
!! Define cache size and average fill level.
......@@ -5489,6 +5490,7 @@ SUBROUTINE vmprep(msize)
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpi) :: ncon
!
INTEGER(mpl), INTENT(IN) :: msize(2)
......@@ -5524,6 +5526,7 @@ SUBROUTINE vmprep(msize)
! variable-width band matrix or diagonal matrix for parameters
! followed by rectangular matrix for constraints
! followed by symmetric matrix for constraints
ncon=nagb-nvgb
IF(mbandw > 0) THEN ! variable-width band matrix
length=nagb
CALL mpalloc(indPreCond,length,'pointer-array variable-band matrix')
......@@ -5536,10 +5539,10 @@ SUBROUTINE vmprep(msize)
DO i=nvgb+1,nagb ! reset
indPreCond(i)=0
END DO
length=indPreCond(nvgb)+ncgb*nvgb+(ncgb*ncgb+ncgb)/2
length=indPreCond(nvgb)+ncon*nvgb+(ncon*ncon+ncon)/2
CALL mpalloc(matPreCond,length,'variable-band matrix')
ELSE ! default preconditioner
length=nvgb+ncgb*nvgb+(ncgb*ncgb+ncgb)/2
length=nvgb+ncon*nvgb+(ncon*ncon+ncon)/2
CALL mpalloc(matPreCond,length,'default preconditioner matrix')
END IF
END IF
......@@ -5894,7 +5897,7 @@ SUBROUTINE mminrs
IF(icalcm == 1) THEN
IF(nfgb < nvgb) CALL qlpssq(avprd0,matPreCond,mbandw,.true.) ! transform preconditioner matrix
WRITE(lun,*) 'MMINRS: EQUDEC started', nprecond(2), nprecond(1)
CALL equdec(nprecond(2),nprecond(1),matPreCond,indPreCond,nrkd,nrkd2)
CALL equdec(nprecond(2),nprecond(1),lprecm,matPreCond,indPreCond,nrkd,nrkd2)
WRITE(lun,*) 'MMINRS: EQUDEC ended ', nrkd, nrkd2
END IF
CALL minres(nfgb, avprod, mvsolv, workspaceD, shift, checka ,.TRUE. , &
......@@ -5990,7 +5993,7 @@ SUBROUTINE mminrsqlp
IF(icalcm == 1) THEN
IF(nfgb < nvgb) CALL qlpssq(avprd0,matPreCond,mbandw,.true.) ! transform preconditioner matrix
WRITE(lun,*) 'MMINRS: EQUDEC started', nprecond(2), nprecond(1)
CALL equdec(nprecond(2),nprecond(1),matPreCond,indPreCond,nrkd,nrkd2)
CALL equdec(nprecond(2),nprecond(1),lprecm,matPreCond,indPreCond,nrkd,nrkd2)
WRITE(lun,*) 'MMINRS: EQUDEC ended ', nrkd, nrkd2
END IF
......@@ -6054,6 +6057,7 @@ SUBROUTINE mvsolv(n,x,y) ! solve M*y = x
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(IN) :: x(n)
REAL(mpd), INTENT(OUT) :: y(n)
SAVE
! ...
y=x ! copy to output vector
......@@ -6177,7 +6181,11 @@ SUBROUTINE xloopn !
ELSE IF(mbandw < 0) THEN
WRITE(lunp,121) 'pre-conditioning:','none!'
ELSE IF(mbandw > 0) THEN
WRITE(lunp,121) 'pre-conditioning=','band-matrix'
IF(lprecm > 0) THEN
WRITE(lunp,121) 'pre-conditioning=','skyline-matrix (rank preserving)'
ELSE
WRITE(lunp,121) 'pre-conditioning=','band-matrix'
ENDIF
END IF
END IF
IF(regpre == 0.0_mpd.AND.npresg == 0) THEN
......@@ -6209,7 +6217,7 @@ SUBROUTINE xloopn !
121 FORMAT(1X,a40,3X,a)
122 FORMAT(1X,a40,2X,i2,a)
122 FORMAT(1X,a40,3X,i0,a)
123 FORMAT(1X,a40,2X,e9.2)
124 FORMAT(1X,a40,3X,f5.1,a)
END DO
......@@ -7658,6 +7666,7 @@ SUBROUTINE intext(text,nline)
IF(mat >= (npat-npat/5)) THEN
IF(nums > 0) mbandw=NINT(dnum(1),mpi)
IF(mbandw < 0) mbandw=-1
IF(nums > 1) lprecm=NINT(dnum(2),mpi)
RETURN
END IF
......
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