Commit 85c88543 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Cleanup (diagonalization with elimination)

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@204 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent 6cc4fbaf
......@@ -321,7 +321,7 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk,mon) !
END IF
END DO
! parallelize row loop
! slot of 128 'J' for next idle thread
! slot of 128 'J' for next idle thread (optimized on Intel Xeon)
!$OMP PARALLEL DO &
!$OMP PRIVATE(JL,VJK,J8) &
!$OMP SCHEDULE(DYNAMIC,128)
......@@ -924,18 +924,19 @@ SUBROUTINE chdec2(g,n,nrank,evmax,evmin,mon)
END IF
ii=ii-i
! parallelize row loop
! slot of 128 'J' for next idle thread
! slot of 32 'J' for next idle thread (optimized on Intel Xeon)
!$OMP PARALLEL DO &
!$OMP PRIVATE(RATIO,JJ) &
!$OMP SCHEDULE(DYNAMIC,128)
!$OMP SCHEDULE(DYNAMIC,32)
DO j=1,i-1
jj=(INT(j-1,mpl)*INT(j,mpl))/2
ratio=g(ii+j)*g(ii+i) ! (I,J) (I,I)
IF (ratio == 0.0_mpd) CYCLE
jj=(INT(j-1,mpl)*INT(j,mpl))/2
g(jj+1:jj+j)=g(jj+1:jj+j)-g(ii+1:ii+j)*ratio ! (K,J) (K,I)
END DO ! J
!$OMP END PARALLEL DO
g(ii+1:ii+i-1)=g(ii+1:ii+i-1)*g(ii+i) ! (I,J)
END DO ! I
END DO ! I
END SUBROUTINE chdec2
......
......@@ -549,10 +549,9 @@ SUBROUTINE qlssq(aprod,A,t)
ioff2=ioffb
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
A(ioff2+1:ioff2+i)=A(ioff2+1:ioff2+i)+2.0_mpd* &
((2.0_mpd*matV(ioff1+i)*vtAv-Av(i))*matV(ioff1+1:ioff1+i)-Av(1:i)*matV(ioff1+i))
ioff2=ioff2+i
END DO
! off diagonal block
DO i=kn+1,nparb
......
......@@ -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-08-02 target
!! svn checkout http://svnsrv.desy.de/public/MillepedeII/tags/V04-08-03 target
!!
!! 2. Create **Pede** executable (in \a target directory):
!!
......@@ -7710,6 +7710,7 @@ SUBROUTINE mdiags
INTEGER(mpi) :: imin
INTEGER(mpl) :: ioff1
INTEGER(mpi) :: j
INTEGER(mpi) :: last
INTEGER(mpi) :: lun
INTEGER(mpi) :: nmax
INTEGER(mpi) :: nmin
......@@ -7758,7 +7759,7 @@ SUBROUTINE mdiags
nmax=INT(1.0+LOG10(REAL(workspaceEigenValues(1),mps)),mpi) ! > log of largest eigenvalue
imin=1
DO i=nagb,1,-1
DO i=nfgb,1,-1
IF(workspaceEigenValues(i) > 0.0_mpd) THEN
imin=i ! index of smallest pos. eigenvalue
EXIT
......@@ -7771,7 +7772,7 @@ SUBROUTINE mdiags
END DO
CALL hmpdef(7,REAL(nmin,mps),REAL(ntop,mps), 'log10 of positive eigenvalues')
DO idia=1,nagb
DO idia=1,nfgb
IF(workspaceEigenValues(idia) > 0.0_mpd) THEN ! positive
evalue=LOG10(REAL(workspaceEigenValues(idia),mps))
CALL hmpent(7,evalue)
......@@ -7782,7 +7783,7 @@ SUBROUTINE mdiags
iast=MAX(1,imin-60)
CALL gmpdef(3,2,'low-value end of eigenvalues')
DO i=iast,nagb
DO i=iast,nfgb
evalue=REAL(workspaceEigenValues(i),mps)
CALL gmpxy(3,REAL(i,mps),evalue)
END DO
......@@ -7796,28 +7797,29 @@ SUBROUTINE mdiags
IF(workspaceEigenValues(i) < 0.0_mpd) workspaceDiagonalization(i)=-workspaceDiagonalization(i)
END IF
END DO
last=min(nfgb,nvgb)
WRITE(lun,*) ' '
WRITE(lun,*) 'The first (largest) eigenvalues ...'
WRITE(lun,102) (workspaceEigenValues(i),i=1,MIN(20,nagb))
WRITE(lun,*) ' '
WRITE(lun,*) 'The last eigenvalues ... up to',nvgb
WRITE(lun,102) (workspaceEigenValues(i),i=MAX(1,nvgb-19),nvgb)
WRITE(lun,*) 'The last eigenvalues ... up to',last
WRITE(lun,102) (workspaceEigenValues(i),i=MAX(1,last-19),last)
WRITE(lun,*) ' '
IF(nagb > nvgb) THEN
WRITE(lun,*) 'The eigenvalues from',nvgb+1,' to',nagb
WRITE(lun,102) (workspaceEigenValues(i),i=nvgb+1,nagb)
WRITE(lun,*) ' '
ENDIF
WRITE(lun,*) 'Log10 + 3 of ',nagb,' eigenvalues in decreasing', ' order'
WRITE(lun,*) 'Log10 + 3 of ',nfgb,' eigenvalues in decreasing', ' order'
WRITE(lun,*) '(for Eigenvalue < 0.001 the value 0.0 is shown)'
WRITE(lun,101) (workspaceDiagonalization(i),i=1,nagb)
WRITE(lun,101) (workspaceDiagonalization(i),i=1,nfgb)
IF(workspaceDiagonalization(nfgb) < 0) WRITE(lun,*) 'Negative values are ', &
'printed for negative eigenvalues'
CALL devsig(nfgb,workspaceEigenValues,workspaceEigenVectors,globalVector,workspaceDiagonalization)
WRITE(lun,*) ' '
WRITE(lun,*) nvgb,' significances: insignificant if ', &
WRITE(lun,*) last,' significances: insignificant if ', &
'compatible with N(0,1)'
WRITE(lun,101) (workspaceDiagonalization(i),i=1,nvgb)
WRITE(lun,101) (workspaceDiagonalization(i),i=1,last)
101 FORMAT(10F7.1)
......
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