Commit 0590f1f1 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

solution method MINRES-QLP added

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@117 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent b67b3d88
......@@ -82,7 +82,8 @@ 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 \
minresDataModule.o minresModule.o randoms.o vertpr.o linesrch.o Dbandmatrix.o pede.o
minresDataModule.o minresModule.o minresqlpDataModule.o minresqlpBlasModule.o minresqlpModule.o \
randoms.o vertpr.o linesrch.o Dbandmatrix.o pede.o
#
# Chose flags/object files for C-binary support:
#
......
......@@ -576,7 +576,6 @@ contains
Arnorml = rnorml*rootl ! ||A r_{k-1} ||
relArnorml = rootl / Anorm; ! ||Ar|| / (||A|| ||r||)
!relArnorml = Arnorml / Anorm; ! ||Ar|| / ||A||
print *, ' iter ', itn, anorm, ynorm, qrnorm, relArnorml, rtol
! Estimate cond(A).
! In this version we look at the diagonals of R in the
......
!> \file
!! MINRES-QLP BLAS subroutines.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! File minresqlpBlasModule.f90
!
! This file contains the following BLAS subroutines
! ddot, dnrm2
! required by subroutine MINRESQLP.
!
! Contributors:
! Sou-Cheng Choi <sctchoi@uchicago.edu>
! Computation Institute (CI)
! University of Chicago
! Chicago, IL 60637, USA
!
! Michael Saunders <saunders@stanford.edu>
! Systems Optimization Laboratory (SOL)
! Stanford University
! Stanford, CA 94305-4026, USA
!
! History:
! 24 Sep 2007: All parameters declared with correct intent
! to avoid compiler warnings.
! 24 Oct 2007: Use real(8) instead of double precision or -r8.
! 24 May 2011: Use a module to package the BLAS subroutines. Use real(dp)
! instead of real(8), where dp is a constant defined in
! minresqlpDataModule and used in other program units.
! 12 Jul 2011: Created complex version zminresqlpBlasModule.f90
! from real version minresqlpBlasModule.f90.
! 03 Aug 2013: dp constants 0.d0 and 1.d0 defined with _dp.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
module minresqlpBlasModule
use minresqlpDataModule, only : dp, ip, zero, one
implicit none
public :: ddot, dnrm2
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!> DDOT forms the dot product of two vectors.
!
! Discussion:
! This routine uses double precision real arithmetic.
! This routine uses unrolled loops for increments equal to one.
!
! Modified:
! 16 May 2005
!
! Author:
! Jack Dongarra
! Fortran90 translation by John Burkardt.
!
! Reference:
! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
! LINPACK User's Guide,
! SIAM, 1979,
! ISBN13: 978-0-898711-72-1,
! LC: QA214.L56.
!
! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
! Algorithm 539,
! Basic Linear Algebra Subprograms for Fortran Usage,
! ACM Transactions on Mathematical Software,
! Volume 5, Number 3, September 1979, pages 308-323.
!
! Parameters:
!
! Input, integer N, the number of entries in the vectors.
!
! Input, real ( kind = dp ) DX(*), the first vector.
!
! Input, integer INCX, the increment between successive entries in DX.
!
! Input, real ( kind = dp ) DY(*), the second vector.
!
! Input, integer INCY, the increment between successive entries in DY.
!
! Output, real ( kind = dp ) DDOT, the sum of the product of the
! corresponding entries of DX and DY.
real(dp) function ddot(n,dx,incx,dy,incy)
implicit none
integer(ip), intent(in) :: n,incx,incy
real(dp), intent(in) :: dx(*),dy(*)
real(dp) :: dtemp
integer(ip) :: i,ix,iy,m
ddot = zero
dtemp = zero
if ( n <= 0 ) then
return
end if
! Code for unequal increments or equal increments
! not equal to 1.
if ( incx /= 1 .or. incy /= 1 ) then
if ( 0 <= incx ) then
ix = 1
else
ix = ( - n + 1 ) * incx + 1
end if
if ( 0 <= incy ) then
iy = 1
else
iy = ( - n + 1 ) * incy + 1
end if
do i = 1, n
dtemp = dtemp + dx(ix) * dy(iy)
ix = ix + incx
iy = iy + incy
end do
! Code for both increments equal to 1.
else
m = mod ( n, 5 )
do i = 1, m
dtemp = dtemp + dx(i) * dy(i)
end do
do i = m+1, 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
end if
ddot = dtemp
return
end function ddot
!*****************************************************************************
!
!> DNRM2 returns the euclidean norm of a vector.
!
! Discussion:
! This routine uses real(dp) real arithmetic.
! DNRM2 ( X ) = sqrt ( X' * X )
!
! Modified:
! 16 May 2005
!
! Author:
! Sven Hammarling
! Fortran90 translation by John Burkardt.
!
! Reference:
! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
! LINPACK User's Guide,
! SIAM, 1979,
! ISBN13: 978-0-898711-72-1,
! LC: QA214.L56.
!
! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
! Algorithm 539,
! Basic Linear Algebra Subprograms for Fortran Usage,
! ACM Transactions on Mathematical Software,
! Volume 5, Number 3, September 1979, pages 308-323.
!
! Parameters:
!
! Input, integer N, the number of entries in the vector.
!
! Input, real ( kind = dp ) X(*), the vector whose norm is to be computed.
!
! Input, integer INCX, the increment between successive entries of X.
!
! Output, real ( kind = dp ) DNRM2, the Euclidean norm of X.
real(dp) function dnrm2 ( n, x, incx )
implicit none
integer(ip), intent(in) :: n,incx
real(dp), intent(in) :: x(*)
integer(ip) :: ix
real(dp) :: ssq,absxi,norm,scale
if ( n < 1 .or. incx < 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 = 1_dp + ssq * ( scale / absxi )**2
scale = absxi
else
ssq = ssq + ( absxi / scale )**2
end if
end if
end do
norm = scale * sqrt ( ssq )
end if
dnrm2 = norm
return
end function dnrm2
end module minresqlpBlasModule
!> \file
!! MINRES-QLP (data) definitions.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! File minresqlpDataModule.f90
!
!> Defines precision and range in real(kind=dp) and integer(kind=ip) for
!! portability and a few constants for use in other modules.
!
!
! Authors:
! Sou-Cheng Choi <sctchoi@uchicago.edu>
! Computation Institute (CI)
! University of Chicago
! Chicago, IL 60637, USA
!
! Michael Saunders <saunders@stanford.edu>
! Systems Optimization Laboratory (SOL)
! Stanford University
! Stanford, CA 94305-4026, USA
!
! History:
! 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 minresqlpDataModule
! at the beginning of modules AND inside interfaces.
! 20 Aug 2012: (1) Added single real kind 'sp' and integer kind 'ip'.
! (2) Added smallest and largest real positive 'realmin'
! and 'realmax'.
! (3) Added single precision kind 'sp'.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
module minresqlpDataModule
use mpdef, only: mpi, mps, mpd
implicit none
intrinsic :: selected_real_kind, selected_int_kind, tiny, huge
! The following reals are provided for portability. Do not use 'DOUBLE PRECISION'.
integer, parameter, public :: dp = mpd !selected_real_kind(15,307) ! 64-bit real, default
integer, parameter, public :: sp = mps !selected_real_kind(6,37) ! 32-bit real
!integer, parameter, public :: qp = selected_real_kind(33,4931) !128-bit real
integer, parameter, public :: ip = mpi !selected_int_kind(9) ! R: (-10^R, 10^R)
real(dp), parameter, public :: zero = 0.0_dp, one = 1.0_dp, eps = epsilon(zero)
real(dp), parameter, public :: realmin = tiny(one), realmax = huge(one)
integer, parameter, public :: prcsn = precision(zero) ! first argument of selected_real_kind
! WARN: turning on debug could significantly slow down the program due to file output
logical, public :: debug = .false.
logical, public :: testSymortho = .true., testMtx = .true.
end module minresqlpDataModule
This diff is collapsed.
......@@ -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 minresmodule::minres "MINRES")
INTEGER(mpi) :: metsol=0 !< solution method (1: inversion, 2: diagonalization, 3: \ref minresqlpmodule::minresqlp "MINRES-QLP")
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)
......@@ -36,7 +36,9 @@ MODULE mpmod
INTEGER(mpi) :: isubit=0 !< subito flag '-s'
REAL(mps) :: wolfc1=0.0!< C_1 of strong Wolfe condition
REAL(mps) :: wolfc2=0.0!< C_2 of strong Wolfe condition
REAL(mpd) :: mrestl=1.0D-06 !< tolerance criterion for MINRES
REAL(mpd) :: mrestl=1.0E-06 !< tolerance criterion for MINRES-QLP
REAL(mpd) :: mrtcnd=1.0E+07 !< transition (QR -> QLP) (matrix) condition for MINRES-QLP
INTEGER(mpi) :: mrmode=0 !< MINRES-QLP mode (0: QR+QLP, 1: only QR, 2: only QLP factorization)
INTEGER(mpi) :: nofeas=0 !< flag for skipping making parameters feasible
INTEGER(mpi) :: nhistp=0 !< flag for histogram printout
REAL(mps) :: delfun=0.0!< expected function change
......
This diff is collapsed.
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