solution method MINRES-QLP added

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@117 3547b9b0-65b8-46d3-b95d-921b3f43af62

Makefile

@@ -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:
 #

minresModule.f90

@@ -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

minresqlpBlasModule.f90

+!> \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

minresqlpDataModule.f90

+!> \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

mpmod.f90

@@ -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