Skip to content
GitLab
Commit 24ce62c3 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

MINRES updated to f90 version (2007) 

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@116 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent 863de05a
Expand all Hide whitespace changes
Inline Side-by-side
Makefile
View file @ 24ce62c3
......@@ -82,7 +82,7 @@ L_FLAGS = -Wall -O3
# objects for this project
#
USER_OBJ_PEDE = mpdef.o mpdalc.o mpmod.o mpbits.o mptest1.o mptest2.o mille.o mpnum.o mptext.o mphistab.o \
minresblas.o minres.o randoms.o vertpr.o linesrch.o Dbandmatrix.o pede.o
minresDataModule.o minresModule.o randoms.o vertpr.o linesrch.o Dbandmatrix.o pede.o
#
# Chose flags/object files for C-binary support:
#
......
minres.f90 deleted 100644 → 0
View file @ 863de05a
This diff is collapsed.
minresDataModule.f90 0 → 100644
View file @ 24ce62c3
!> \file
!! MINRES (data) definitions.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! File minresDataModule.f90
!
! Defines real(kind=dp) and a few constants for use in other modules.
!
! 14 Oct 2007: First version implemented after realizing -r8 is not
! a standard compiler option.
! 15 Oct 2007: Temporarily used real(8) everywhere.
! 16 Oct 2007: Found that we need
! use minresDataModule
! at the beginning of modules AND inside interfaces.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Defines real(kind=dp) and a few constants for use in other modules.
module minresDataModule
use mpdef, only: mpd
implicit none
intrinsic :: selected_real_kind
integer, parameter, public :: dp = mpd !selected_real_kind(15)
real(kind=dp), parameter, public :: zero = 0.0_dp, one = 1.0_dp
end module minresDataModule
minresModule.f90 0 → 100644
View file @ 24ce62c3
This diff is collapsed.
minresblas.f90 deleted 100644 → 0
View file @ 863de05a
! Code converted using TO_F90 by Alan Miller
! Date: 2012-03-16 Time: 11:07:10
!> \file
!! BLAS routines for MINRES.
!!\verbatim
!! minresblas.f
!!
!! This file contains Level 1 BLAS from netlib, Thu May 16 1991
!! (with declarations of the form dx(1) changed to dx(*)):
!! daxpy dcopy ddot
!! Also
!! dnrm2 (from NAG,I think).
!!
!! Also a few utilities to avoid some of the
!! loops in MINRES (so the debugger can step past them quickly):
!! daxpy2 dload2 dscal2
!!
!! 15 Jul 2003: dnrm2 is now the NAG version.
!!\endverbatim
!> Constant times a vector plus a vector.
!!
!! Uses unrolled loops for increments equal to one.
!! (jack dongarra, linpack, 3/11/78)
!!
!! \param[in] n size of vectors
!! \param[in] da scalar constant
!! \param[in] dx input vector
!! \param[in] incx increment for dx
!! \param[out] dy output vector
!! \param[in] incy increment for dy
SUBROUTINE daxpy(n,da,dx,incx,dy,incy)
USE mpdef
REAL(mpd) :: dx(*),dy(*),da
INTEGER(mpi) :: i,incx,incy,ix,iy,m,mp1,n
IF(n <= 0)RETURN
IF (da == 0.0_mpd) RETURN
IF(incx == 1.AND.incy == 1)GO TO 20
! code for unequal increments or equal increments
! not equal to 1
ix = 1
iy = 1
IF(incx < 0)ix = (-n+1)*incx + 1
IF(incy < 0)iy = (-n+1)*incy + 1
DO i = 1,n
dy(iy) = dy(iy) + da*dx(ix)
ix = ix + incx
iy = iy + incy
END DO
RETURN
! code for both increments equal to 1
! clean-up loop
20 m = MOD(n,4)
IF( m == 0 ) GO TO 40
DO i = 1,m
dy(i) = dy(i) + da*dx(i)
END DO
IF( n < 4 ) RETURN
40 mp1 = m + 1
DO i = mp1,n,4
dy(i) = dy(i) + da*dx(i)
dy(i + 1) = dy(i + 1) + da*dx(i + 1)
dy(i + 2) = dy(i + 2) + da*dx(i + 2)
dy(i + 3) = dy(i + 3) + da*dx(i + 3)
END DO
END SUBROUTINE daxpy
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Copies a vector, x, to a vector, y.
!!
!! Uses unrolled loops for increments equal to one.
!! (jack dongarra, linpack, 3/11/78)
!!
!! \param[in] n size of vectors
!! \param[in] dx input vector
!! \param[in] incx increment for dx
!! \param[out] dy output vector
!! \param[in] incy increment for dy
SUBROUTINE dcopy(n,dx,incx,dy,incy)
USE mpdef
REAL(mpd) :: dx(*),dy(*)
INTEGER(mpi) :: i,incx,incy,ix,iy,m,mp1,n
IF(n <= 0)RETURN
IF(incx == 1.AND.incy == 1)GO TO 20
! code for unequal increments or equal increments
! not equal to 1
ix = 1
iy = 1
IF(incx < 0)ix = (-n+1)*incx + 1
IF(incy < 0)iy = (-n+1)*incy + 1
DO i = 1,n
dy(iy) = dx(ix)
ix = ix + incx
iy = iy + incy
END DO
RETURN
! code for both increments equal to 1
! clean-up loop
20 m = MOD(n,7)
IF( m == 0 ) GO TO 40
DO i = 1,m
dy(i) = dx(i)
END DO
IF( n < 7 ) RETURN
40 mp1 = m + 1
DO i = mp1,n,7
dy(i) = dx(i)
dy(i + 1) = dx(i + 1)
dy(i + 2) = dx(i + 2)
dy(i + 3) = dx(i + 3)
dy(i + 4) = dx(i + 4)
dy(i + 5) = dx(i + 5)
dy(i + 6) = dx(i + 6)
END DO
END SUBROUTINE dcopy
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Forms the dot product of two vectors.
!!
!! Uses unrolled loops for increments equal to one.
!! (jack dongarra, linpack, 3/11/78)
!!
!! \param[in] n size of vectors
!! \param[in] dx input vector
!! \param[in] incx increment for dx
!! \param[in] dy input vector
!! \param[in] incy increment for dy
!! \return dot porduct dx*dy
REAL(mpd) FUNCTION ddot(n,dx,incx,dy,incy)
USE mpdef
REAL(mpd) :: dx(*),dy(*),dtemp
INTEGER(mpi) :: i,incx,incy,ix,iy,m,mp1,n
ddot = 0.0_mpd
dtemp = 0.0_mpd
IF(n <= 0)RETURN
IF(incx == 1.AND.incy == 1)GO TO 20
! code for unequal increments or equal increments
! not equal to 1
ix = 1
iy = 1
IF(incx < 0)ix = (-n+1)*incx + 1
IF(incy < 0)iy = (-n+1)*incy + 1
DO i = 1,n
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
END DO
ddot = dtemp
RETURN
! code for both increments equal to 1
! clean-up loop
20 m = MOD(n,5)
IF( m == 0 ) GO TO 40
DO i = 1,m
dtemp = dtemp + dx(i)*dy(i)
END DO
IF( n < 5 ) GO TO 60
40 mp1 = m + 1
DO i = mp1,n,5
dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + &
dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
END DO
60 ddot = dtemp
END FUNCTION ddot
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Euclidean vector norm.
!!
!! dnrm2 returns the Euclidean norm of a vector via the function
!! name, so that dnrm2 := sqrt( x'*x ).
!!
!! 15 Jul 2003: dnrm2 obtained from SNOPT src (probably from NAG).
!! s1flmx replaced by safe large number.
!!
!! \param[in] n size of vectors
!! \param[in] x input vector
!! \param[in] incx increment for x
!! \return vector norm
REAL(mpd) FUNCTION dnrm2 ( n, x, incx )
USE mpdef
IMPLICIT REAL(mpd) (A-H,O-Z)
INTEGER(mpi) :: incx, n
REAL(mpd) :: x(*)
!!! REAL(mpd) s1flmx
PARAMETER (one = 1.0_mpd, zero = 0.0_mpd )
REAL(mpd) :: norm
INTRINSIC ABS
! ------------------------------------------------------------------
! flmax = s1flmx( )
flmax = 1.0E+50_mpd
IF ( n < 1) THEN
norm = zero
ELSE IF (n == 1) THEN
norm = ABS( x(1) )
ELSE
scale = zero
ssq = one
DO ix = 1, 1+(n-1)*incx, incx
IF (x(ix) /= zero) THEN
absxi = ABS( x(ix) )
IF (scale < absxi) THEN
ssq = one + ssq*(scale/absxi)**2
scale = absxi
ELSE
ssq = ssq + (absxi/scale)**2
END IF
END IF
END DO
sqt = SQRT( ssq )
IF (scale < flmax/sqt) THEN
norm = scale*sqt
ELSE
norm = flmax
END IF
END IF
dnrm2 = norm
END FUNCTION dnrm2
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Set z = a*x + y.
!!
!! 31 May 1999: First version written for MINRES.
!!
!! \param[in] n size of vectors
!! \param[in] a scalar constant
!! \param[in] x input vector
!! \param[in] y input vector
!! \param[out] z output vector
SUBROUTINE daxpy2( n, a, x, y, z )
USE mpdef
IMPLICIT NONE
INTEGER(mpi) :: n
REAL(mpd) :: a, x(n), y(n), z(n)
INTEGER(mpi) :: i
DO i = 1, n
z(i) = a*x(i) + y(i)
END DO
END SUBROUTINE daxpy2
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Set x = constant.
!!
!! Load all elements of x with const.
!!
!! \param[in] n size of vectors
!! \param[in] const scalar constant
!! \param[out] x output vector
SUBROUTINE dload2( n, const, x )
USE mpdef
IMPLICIT NONE
INTEGER(mpi) :: n
REAL(mpd) :: const, x(n)
! ------------------------------------------------------------------
! dload2
! ------------------------------------------------------------------
INTEGER(mpi) :: i
DO i = 1, n
x(i) = const
END DO
END SUBROUTINE dload2
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Set y = a*x.
!!
!! \param[in] n size of vectors
!! \param[in] a scalar constant
!! \param[in] x input vector
!! \param[out] y output vector
SUBROUTINE dscal2( n, a, x, y )
USE mpdef
IMPLICIT NONE
INTEGER(mpi) :: n
REAL(mpd) :: a, x(n), y(n)
INTEGER(mpi) :: i
DO i = 1, n
y(i) = a*x(i)
END DO
END SUBROUTINE dscal2
mpmod.f90
View file @ 24ce62c3
......@@ -12,7 +12,7 @@ MODULE mpmod
SAVE
! steering parameters
INTEGER(mpi) :: ictest=0 !< test mode '-t'
INTEGER(mpi) :: metsol=0 !< solution method (1: inversion, 2: diagonalization, 3: \ref minres "MINRES")
INTEGER(mpi) :: metsol=0 !< solution method (1: inversion, 2: diagonalization, 3: \ref minresmodule::minres "MINRES")
INTEGER(mpi) :: matsto=2 !< (global) matrix storage mode (1: full, 2: sparse)
INTEGER(mpi) :: mprint=1 !< print flag (0: minimal, 1: normal, >1: more)
INTEGER(mpi) :: mdebug=0 !< debug flag (number of records to print)
......@@ -142,7 +142,6 @@ MODULE mpmod
REAL(mpd), DIMENSION(:), ALLOCATABLE :: workspaceDiagonalization !< workspace diag.
REAL(mpd), DIMENSION(:), ALLOCATABLE :: workspaceEigenValues !< workspace eigen values
REAL(mpd), DIMENSION(:), ALLOCATABLE :: workspaceEigenVectors !< workspace eigen vectors
REAL(mpd), DIMENSION(:), ALLOCATABLE :: workspaceMinres !< workspace MINRES
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: workspaceI !< (general) workspace (I)
! constraint matrix, residuals
REAL(mpd), DIMENSION(:), ALLOCATABLE :: matConsProduct !< product matrix of constraints
......
pede.f90
View file @ 24ce62c3
......@@ -108,7 +108,7 @@
!! eigenvalue (and eigenvector) for all global parameters.
!! \subsection ch-minres Minimal Residual Method (MINRES)
!! The solution is obtained by minimizing \f$\Vert\Vek{A}\cdot\Vek{x}-\Vek{b}\Vert_2\f$
!! iteratively. \ref minres "MINRES" is a special case of the
!! iteratively. \ref minresmodule::minres "MINRES" is a special case of the
!! generalized minimal residual method (\ref an-gmres "GMRES") for symmetric matrices.
!! Preconditioning with a band matrix of zero or finite
!! \ref mpmod::mbandw "bandwidth" is possible.
......@@ -236,7 +236,7 @@
!!
!! \subsection cmd-bandwidth bandwidth
!! Set band width \ref mpmod::mbandw "mbandw" for
!! \ref minres "MINRES" preconditioner to \a number1.
!! \ref minresmodule::minres "MINRES" preconditioner to \a number1.
!! \subsection cmd-cache cache
!! Set (read+write) cache size \ref mpmod::ncache "ncache" to \a number1.
!! Define cache size and average fill level.
......@@ -308,7 +308,7 @@
!! (minimum) number of iterations \ref mpmod::mitera "mitera" to \a number1,
!! convergence limit \ref mpmod::dflim "dflim" to \a number2.
!! \subsection cmd-mrestol mrestol
!! Set tolerance criterion \ref mpmod::mrestl "mrestl" for \ref minres "MINRES"
!! Set tolerance criterion \ref mpmod::mrestl "mrestl" for \ref minresmodule::minres "MINRES"
!! to \a number1 (\f$10^{-10}\f$ .. \f$10^{-4}\f$).
!! \subsection cmd-nofeasiblestart nofeasiblestart
!! Set flag \ref mpmod::nofeas "nofeas" for \ref an-nofeas "skipping"
......@@ -758,7 +758,7 @@ PROGRAM mptwo
END PROGRAM mptwo ! Mille
!> Error for single global parameter from \ref minres "MINRES".
!> Error for single global parameter from \ref minresmodule::minres "MINRES".
!!
!! Calculate single row 'x_i' from inverse matrix by solving A*x_i=b
!! with b=0 except b_i=1.
......@@ -767,6 +767,7 @@ END PROGRAM mptwo ! Mille
SUBROUTINE solglo(ivgbi)
USE mpmod
USE minresModule, ONLY: minres
IMPLICIT NONE
REAL(mps) :: dpa
......@@ -787,6 +788,7 @@ SUBROUTINE solglo(ivgbi)
REAL(mpd) :: rtol
REAL(mpd) :: anorm
REAL(mpd) :: acond
REAL(mpd) :: arnorm
REAL(mpd) :: rnorm
REAL(mpd) :: ynorm
REAL(mpd) :: gmati
......@@ -816,31 +818,19 @@ SUBROUTINE solglo(ivgbi)
rtol = mrestl ! from steering
checka=.FALSE.
IF(mbandw == 0) THEN ! default preconditioner
CALL minres(nagb,globalVector, &
workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), &
workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), &
globalCorrections,workspaceMinres, avprod, mcsolv, checka ,.TRUE. , shift, &
nout , itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
ELSE IF(mbandw > 0) THEN ! band matrix preconditioner
CALL minres(nagb,globalVector, &
workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), &
workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), &
globalCorrections,workspaceMinres, avprod, mvsolv, checka ,.TRUE. , shift, &
nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
CALL minres(nagb, avprod, mcsolv, globalVector, shift, checka ,.TRUE. , &
globalCorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
ELSE IF(mbandw > 0) THEN ! band matrix preconditioner
CALL minres(nagb, avprod, mvsolv, globalVector, shift, checka ,.TRUE. , &
globalCorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
ELSE
CALL minres(nagb,globalVector, &
workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), &
workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), &
globalCorrections,workspaceMinres, avprod, mvsolv, checka ,.FALSE., shift, &
nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
CALL minres(nagb, avprod, mvsolv, globalVector, shift, checka ,.FALSE. , &
globalCorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
END IF
! subroutine MINRES( n, b, r1, r2, v, w, w1, w2, x, y,
! $ AVPROD, Msolve, .TRUE. , .FALSE. , SHIFT,
! $ NOUT , ITNLIM, rtol,
! $ ISTOP, ITN, Anorm, Acond, rnorm, ynorm )
par=REAL(globalParameter(itgbi),mps)
dpa=par-globalParStart(itgbi)
gmati=globalCorrections(ivgbi)
......@@ -3441,7 +3431,7 @@ END SUBROUTINE prtglo ! print final log file
!> Product symmetric matrix times vector.
!!
!! A(sym) * X => B. Used by \ref minres "MINRES" method (Is most CPU intensive part).
!! A(sym) * X => B. Used by \ref minresmodule::minres "MINRES" method (Is most CPU intensive part).
!! The matrix A is the global matrix in full symmetric or (compressed) sparse storage.
!!
!! \param [in] n size of matrix
......@@ -4815,7 +4805,6 @@ SUBROUTINE vmprep(msize)
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpi) :: nmats
!
INTEGER(mpl), INTENT(IN) :: msize(2)
......@@ -4876,7 +4865,6 @@ SUBROUTINE vmprep(msize)
CALL mpalloc(workspaceLinesearch,length,'auxiliary array (D2)') ! double aux 2
CALL mpalloc(workspaceI, length,'auxiliary array (I)') ! int aux 1
CALL mpalloc(workspaceMinres,length,'auxiliary array (D7)') ! double aux 7
IF(metsol == 1) THEN
CALL mpalloc(workspaceD,length,'auxiliary array (D1)') ! double aux 1
......@@ -4891,12 +4879,6 @@ SUBROUTINE vmprep(msize)
CALL mpalloc(workspaceEigenVectors,length,'(rotation) matrix U') ! rotation matrix
END IF
IF(metsol >= 3) THEN
nmats=6
length=nmats*nagb
CALL mpalloc(workspaceD,length,'auxiliary array (D1)') ! double aux 1
END IF
END SUBROUTINE vmprep
!> Solution by matrix inversion.
......@@ -5060,7 +5042,7 @@ SUBROUTINE zdiags
END SUBROUTINE zdiags
!> Solution with \ref minres "MINRES".
!> Solution with \ref minresmodule::minres "MINRES".
!!
!! Solve A*x=b by minimizing |A*x-b| iteratively. Parallelized (AVPROD).
!!
......@@ -5068,6 +5050,7 @@ END SUBROUTINE zdiags
SUBROUTINE mminrs
USE mpmod
USE minresModule, ONLY: minres
IMPLICIT NONE
INTEGER(mpi) :: istop
......@@ -5082,6 +5065,7 @@ SUBROUTINE mminrs
REAL(mpd) :: rtol
REAL(mpd) :: anorm
REAL(mpd) :: acond
REAL(mpd) :: arnorm
REAL(mpd) :: rnorm
REAL(mpd) :: ynorm
LOGICAL :: checka
......@@ -5102,29 +5086,20 @@ SUBROUTINE mminrs
CALL precon(ncgb,nvgb,matPreCond,matPreCond, matPreCond(1+nvgb), &
matPreCond(1+nvgb+ncgb*nvgb))
END IF
CALL minres(nagb,globalVector, &
workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), &
workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), &
globalCorrections,workspaceMinres, avprod, mcsolv, checka ,.TRUE. , shift, &
nout , itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
CALL minres(nagb, avprod, mcsolv, globalVector, shift, checka ,.TRUE. , &
globalCorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
ELSE IF(mbandw > 0) THEN ! band matrix preconditioner
IF(icalcm == 1) THEN
WRITE(lun,*) 'MMINRS: EQUDEC started'
CALL equdec(nvgb,ncgb,matPreCond,indPreCond,nrkd,nrkd2)
WRITE(lun,*) 'MMINRS: EQUDEC ended'
END IF
CALL minres(nagb,globalVector, &
workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), &
workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), &
globalCorrections,workspaceMinres, avprod, mvsolv, checka ,.TRUE. , shift, &
nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
CALL minres(nagb, avprod, mvsolv, globalVector, shift, checka ,.TRUE. , &
globalCorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
ELSE
CALL minres(nagb,globalVector, &
workspaceD(1), workspaceD(1+nagb), workspaceD(1+2*nagb), &
workspaceD(1+3*nagb),workspaceD(1+4*nagb),workspaceD(1+5*nagb), &
globalCorrections,workspaceMinres, avprod, mvsolv, checka ,.FALSE., shift, &
nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
CALL minres(nagb, avprod, mvsolv, globalVector, shift, checka ,.FALSE. , &
globalCorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
END IF
iitera=itn
istopa=istop
......@@ -5136,7 +5111,7 @@ END SUBROUTINE mminrs
!> Solution for zero band width preconditioner.
!!
!! Used by \ref minres "MINRES".
!! Used by \ref minresmodule::minres "MINRES".
!!
!! \param[in] n size of vectors
!! \param [in] x rhs vector
......@@ -5156,7 +5131,7 @@ END SUBROUTINE mcsolv
!> Solution for finite band width preconditioner.
!!
!! Used by \ref minres "MINRES".
!! Used by \ref minresmodule::minres "MINRES".
!!
!! \param[in] n size of vectors