remove the fortran code from jrl-walkgen

src/Fortran/CMakeLists.txt

-enable_language (Fortran)
-
-SET(SOURCES_FORTRAN
-  dgesvd.f
-  dgesvx.f
-  dgetrf.f
-  gdgges.f
-  dlamch.f
-  dlapy2.f
-  dgetrf2.f
-  gdtgsen.f
-  gdtgexc.f
-  gdtgex2.f
-  gdgemm.f
-  gsgemm.f
-)
-
-ADD_LIBRARY(sublapack SHARED ${SOURCES_FORTRAN})
-
-INSTALL(TARGETS sublapack DESTINATION ${CMAKE_INSTALL_LIBDIR})

src/Fortran/dgetrf.f

-*> \brief \b DGETRF
-*
-*  =========== DOCUMENTATION ===========
-*
-* Online html documentation available at 
-*            http://www.netlib.org/lapack/explore-html/ 
-*
-*> \htmlonly
-*> Download DGETRF + dependencies 
-*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgetrf.f"> 
-*> [TGZ]</a> 
-*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgetrf.f"> 
-*> [ZIP]</a> 
-*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgetrf.f"> 
-*> [TXT]</a>
-*> \endhtmlonly 
-*
-*  Definition:
-*  ===========
-*
-*       SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
-* 
-*       .. Scalar Arguments ..
-*       INTEGER            INFO, LDA, M, N
-*       ..
-*       .. Array Arguments ..
-*       INTEGER            IPIV( * )
-*       DOUBLE PRECISION   A( LDA, * )
-*       ..
-*  
-*
-*> \par Purpose:
-*  =============
-*>
-*> \verbatim
-*>
-*> DGETRF computes an LU factorization of a general M-by-N matrix A
-*> using partial pivoting with row interchanges.
-*>
-*> The factorization has the form
-*>    A = P * L * U
-*> where P is a permutation matrix, L is lower triangular with unit
-*> diagonal elements (lower trapezoidal if m > n), and U is upper
-*> triangular (upper trapezoidal if m < n).
-*>
-*> This is the right-looking Level 3 BLAS version of the algorithm.
-*> \endverbatim
-*
-*  Arguments:
-*  ==========
-*
-*> \param[in] M
-*> \verbatim
-*>          M is INTEGER
-*>          The number of rows of the matrix A.  M >= 0.
-*> \endverbatim
-*>
-*> \param[in] N
-*> \verbatim
-*>          N is INTEGER
-*>          The number of columns of the matrix A.  N >= 0.
-*> \endverbatim
-*>
-*> \param[in,out] A
-*> \verbatim
-*>          A is DOUBLE PRECISION array, dimension (LDA,N)
-*>          On entry, the M-by-N matrix to be factored.
-*>          On exit, the factors L and U from the factorization
-*>          A = P*L*U; the unit diagonal elements of L are not stored.
-*> \endverbatim
-*>
-*> \param[in] LDA
-*> \verbatim
-*>          LDA is INTEGER
-*>          The leading dimension of the array A.  LDA >= max(1,M).
-*> \endverbatim
-*>
-*> \param[out] IPIV
-*> \verbatim
-*>          IPIV is INTEGER array, dimension (min(M,N))
-*>          The pivot indices; for 1 <= i <= min(M,N), row i of the
-*>          matrix was interchanged with row IPIV(i).
-*> \endverbatim
-*>
-*> \param[out] INFO
-*> \verbatim
-*>          INFO is INTEGER
-*>          = 0:  successful exit
-*>          < 0:  if INFO = -i, the i-th argument had an illegal value
-*>          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
-*>                has been completed, but the factor U is exactly
-*>                singular, and division by zero will occur if it is used
-*>                to solve a system of equations.
-*> \endverbatim
-*
-*  Authors:
-*  ========
-*
-*> \author Univ. of Tennessee 
-*> \author Univ. of California Berkeley 
-*> \author Univ. of Colorado Denver 
-*> \author NAG Ltd. 
-*
-*> \date November 2015
-*
-*> \ingroup doubleGEcomputational
-*
-*  =====================================================================
-      SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
-*
-*  -- LAPACK computational routine (version 3.6.0) --
-*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
-*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2015
-*
-*     .. Scalar Arguments ..
-      INTEGER            INFO, LDA, M, N
-*     ..
-*     .. Array Arguments ..
-      INTEGER            IPIV( * )
-      DOUBLE PRECISION   A( LDA, * )
-*     ..
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION   ONE
-      PARAMETER          ( ONE = 1.0D+0 )
-*     ..
-*     .. Local Scalars ..
-      INTEGER            I, IINFO, J, JB, NB
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL           GDGEMM, DGETRF2, DLASWP, DTRSM, XERBLA
-*     ..
-*     .. External Functions ..
-      INTEGER            ILAENV
-      EXTERNAL           ILAENV
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          MAX, MIN
-*     ..
-*     .. Executable Statements ..
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF( M.LT.0 ) THEN
-         INFO = -1
-      ELSE IF( N.LT.0 ) THEN
-         INFO = -2
-      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
-         INFO = -4
-      END IF
-      IF( INFO.NE.0 ) THEN
-         CALL XERBLA( 'DGETRF', -INFO )
-         RETURN
-      END IF
-*
-*     Quick return if possible
-*
-      IF( M.EQ.0 .OR. N.EQ.0 )
-     $   RETURN
-*
-*     Determine the block size for this environment.
-*
-      NB = ILAENV( 1, 'DGETRF', ' ', M, N, -1, -1 )
-      IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
-*
-*        Use unblocked code.
-*
-         CALL DGETRF2( M, N, A, LDA, IPIV, INFO )
-      ELSE
-*
-*        Use blocked code.
-*
-         DO 20 J = 1, MIN( M, N ), NB
-            JB = MIN( MIN( M, N )-J+1, NB )
-*
-*           Factor diagonal and subdiagonal blocks and test for exact
-*           singularity.
-*
-            CALL DGETRF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO )
-*
-*           Adjust INFO and the pivot indices.
-*
-            IF( INFO.EQ.0 .AND. IINFO.GT.0 )
-     $         INFO = IINFO + J - 1
-            DO 10 I = J, MIN( M, J+JB-1 )
-               IPIV( I ) = J - 1 + IPIV( I )
-   10       CONTINUE
-*
-*           Apply interchanges to columns 1:J-1.
-*
-            CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
-*
-            IF( J+JB.LE.N ) THEN
-*
-*              Apply interchanges to columns J+JB:N.
-*
-               CALL DLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1,
-     $                      IPIV, 1 )
-*
-*              Compute block row of U.
-*
-               CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
-     $                     N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
-     $                     LDA )
-               IF( J+JB.LE.M ) THEN
-*
-*                 Update trailing submatrix.
-*
-                  CALL GDGEMM( 'No transpose', 'No transpose', M-J-JB+1,
-     $                        N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA,
-     $                        A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ),
-     $                        LDA )
-               END IF
-            END IF
-   20    CONTINUE
-      END IF
-      RETURN
-*
-*     End of DGETRF
-*
-      END

src/Fortran/dgetrf2.f

-*> \brief \b DGETRF2
-*
-*  =========== DOCUMENTATION ===========
-*
-* Online html documentation available at 
-*            http://www.netlib.org/lapack/explore-html/ 
-*
-*  Definition:
-*  ===========
-*
-*       RECURSIVE SUBROUTINE DGETRF2( M, N, A, LDA, IPIV, INFO )
-* 
-*       .. Scalar Arguments ..
-*       INTEGER            INFO, LDA, M, N
-*       ..
-*       .. Array Arguments ..
-*       INTEGER            IPIV( * )
-*       DOUBLE PRECISION   A( LDA, * )
-*       ..
-*  
-*
-*> \par Purpose:
-*  =============
-*>
-*> \verbatim
-*>
-*> DGETRF2 computes an LU factorization of a general M-by-N matrix A
-*> using partial pivoting with row interchanges.
-*>
-*> The factorization has the form
-*>    A = P * L * U
-*> where P is a permutation matrix, L is lower triangular with unit
-*> diagonal elements (lower trapezoidal if m > n), and U is upper
-*> triangular (upper trapezoidal if m < n).
-*>
-*> This is the recursive version of the algorithm. It divides
-*> the matrix into four submatrices:
-*>            
-*>        [  A11 | A12  ]  where A11 is n1 by n1 and A22 is n2 by n2
-*>    A = [ -----|----- ]  with n1 = min(m,n)/2
-*>        [  A21 | A22  ]       n2 = n-n1
-*>            
-*>                                       [ A11 ]
-*> The subroutine calls itself to factor [ --- ],
-*>                                       [ A12 ]
-*>                 [ A12 ]
-*> do the swaps on [ --- ], solve A12, update A22,
-*>                 [ A22 ]
-*>
-*> then calls itself to factor A22 and do the swaps on A21.
-*>
-*> \endverbatim
-*
-*  Arguments:
-*  ==========
-*
-*> \param[in] M
-*> \verbatim
-*>          M is INTEGER
-*>          The number of rows of the matrix A.  M >= 0.
-*> \endverbatim
-*>
-*> \param[in] N
-*> \verbatim
-*>          N is INTEGER
-*>          The number of columns of the matrix A.  N >= 0.
-*> \endverbatim
-*>
-*> \param[in,out] A
-*> \verbatim
-*>          A is DOUBLE PRECISION array, dimension (LDA,N)
-*>          On entry, the M-by-N matrix to be factored.
-*>          On exit, the factors L and U from the factorization
-*>          A = P*L*U; the unit diagonal elements of L are not stored.
-*> \endverbatim
-*>
-*> \param[in] LDA
-*> \verbatim
-*>          LDA is INTEGER
-*>          The leading dimension of the array A.  LDA >= max(1,M).
-*> \endverbatim
-*>
-*> \param[out] IPIV
-*> \verbatim
-*>          IPIV is INTEGER array, dimension (min(M,N))
-*>          The pivot indices; for 1 <= i <= min(M,N), row i of the
-*>          matrix was interchanged with row IPIV(i).
-*> \endverbatim
-*>
-*> \param[out] INFO
-*> \verbatim
-*>          INFO is INTEGER
-*>          = 0:  successful exit
-*>          < 0:  if INFO = -i, the i-th argument had an illegal value
-*>          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
-*>                has been completed, but the factor U is exactly
-*>                singular, and division by zero will occur if it is used
-*>                to solve a system of equations.
-*> \endverbatim
-*
-*  Authors:
-*  ========
-*
-*> \author Univ. of Tennessee 
-*> \author Univ. of California Berkeley 
-*> \author Univ. of Colorado Denver 
-*> \author NAG Ltd. 
-*
-*> \date June 2016
-*
-*> \ingroup doubleGEcomputational
-*
-*  =====================================================================
-      RECURSIVE SUBROUTINE DGETRF2( M, N, A, LDA, IPIV, INFO )
-*
-*  -- LAPACK computational routine (version 3.6.1) --
-*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
-*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     June 2016
-*
-*     .. Scalar Arguments ..
-      INTEGER            INFO, LDA, M, N
-*     ..
-*     .. Array Arguments ..
-      INTEGER            IPIV( * )
-      DOUBLE PRECISION   A( LDA, * )
-*     ..
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION   ONE, ZERO
-      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
-*     ..
-*     .. Local Scalars ..
-      DOUBLE PRECISION   SFMIN, TEMP
-      INTEGER            I, IINFO, N1, N2
-*     ..
-*     .. External Functions ..
-      DOUBLE PRECISION   DLAMCH
-      INTEGER            IDAMAX
-      EXTERNAL           DLAMCH, IDAMAX
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL           GDGEMM, DSCAL, DLASWP, DTRSM, XERBLA
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          MAX, MIN
-*     ..
-*     .. Executable Statements ..
-*
-*     Test the input parameters
-*
-      INFO = 0
-      IF( M.LT.0 ) THEN
-         INFO = -1
-      ELSE IF( N.LT.0 ) THEN
-         INFO = -2
-      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
-         INFO = -4
-      END IF
-      IF( INFO.NE.0 ) THEN
-         CALL XERBLA( 'DGETRF2', -INFO )
-         RETURN
-      END IF
-*
-*     Quick return if possible
-*
-      IF( M.EQ.0 .OR. N.EQ.0 )
-     $   RETURN
-
-      IF ( M.EQ.1 ) THEN
-*
-*        Use unblocked code for one row case
-*        Just need to handle IPIV and INFO
-*
-         IPIV( 1 ) = 1
-         IF ( A(1,1).EQ.ZERO )
-     $      INFO = 1
-*
-      ELSE IF( N.EQ.1 ) THEN
-*
-*        Use unblocked code for one column case
-*
-*
-*        Compute machine safe minimum
-*
-         SFMIN = DLAMCH('S')
-*
-*        Find pivot and test for singularity
-*
-         I = IDAMAX( M, A( 1, 1 ), 1 )
-         IPIV( 1 ) = I
-         IF( A( I, 1 ).NE.ZERO ) THEN
-*
-*           Apply the interchange
-*
-            IF( I.NE.1 ) THEN
-               TEMP = A( 1, 1 )
-               A( 1, 1 ) = A( I, 1 )
-               A( I, 1 ) = TEMP
-            END IF
-*
-*           Compute elements 2:M of the column
-*
-            IF( ABS(A( 1, 1 )) .GE. SFMIN ) THEN
-               CALL DSCAL( M-1, ONE / A( 1, 1 ), A( 2, 1 ), 1 )
-            ELSE
-               DO 10 I = 1, M-1
-                  A( 1+I, 1 ) = A( 1+I, 1 ) / A( 1, 1 )
-   10          CONTINUE
-            END IF
-*
-         ELSE
-            INFO = 1
-         END IF
-*
-      ELSE
-*
-*        Use recursive code
-*
-         N1 = MIN( M, N ) / 2
-         N2 = N-N1
-*
-*               [ A11 ]
-*        Factor [ --- ]
-*               [ A21 ]
-*
-         CALL DGETRF2( M, N1, A, LDA, IPIV, IINFO )
-
-         IF ( INFO.EQ.0 .AND. IINFO.GT.0 )
-     $      INFO = IINFO
-*
-*                              [ A12 ]
-*        Apply interchanges to [ --- ]
-*                              [ A22 ]
-*
-         CALL DLASWP( N2, A( 1, N1+1 ), LDA, 1, N1, IPIV, 1 )
-*
-*        Solve A12
-*
-         CALL DTRSM( 'L', 'L', 'N', 'U', N1, N2, ONE, A, LDA, 
-     $               A( 1, N1+1 ), LDA )
-*
-*        Update A22
-*
-         CALL GDGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( N1+1, 1 ), LDA, 
-     $               A( 1, N1+1 ), LDA, ONE, A( N1+1, N1+1 ), LDA )
-*
-*        Factor A22
-*
-         CALL DGETRF2( M-N1, N2, A( N1+1, N1+1 ), LDA, IPIV( N1+1 ),
-     $                 IINFO )
-*
-*        Adjust INFO and the pivot indices
-*
-         IF ( INFO.EQ.0 .AND. IINFO.GT.0 )
-     $      INFO = IINFO + N1
-         DO 20 I = N1+1, MIN( M, N )
-            IPIV( I ) = IPIV( I ) + N1
-   20    CONTINUE
-*
-*        Apply interchanges to A21
-*
-         CALL DLASWP( N1, A( 1, 1 ), LDA, N1+1, MIN( M, N), IPIV, 1 )
-*
-      END IF
-      RETURN
-*
-*     End of DGETRF2
-*
-      END

src/Fortran/dlamch.f

-*> \brief \b DLAMCH
-*
-*  =========== DOCUMENTATION ===========
-*
-* Online html documentation available at 
-*            http://www.netlib.org/lapack/explore-html/ 
-*
-*  Definition:
-*  ===========
-*
-*      DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
-*  
-*
-*> \par Purpose:
-*  =============
-*>
-*> \verbatim
-*>
-*> DLAMCH determines double precision machine parameters.
-*> \endverbatim
-*
-*  Arguments:
-*  ==========
-*
-*> \param[in] CMACH
-*> \verbatim
-*>          Specifies the value to be returned by DLAMCH:
-*>          = 'E' or 'e',   DLAMCH := eps
-*>          = 'S' or 's ,   DLAMCH := sfmin
-*>          = 'B' or 'b',   DLAMCH := base
-*>          = 'P' or 'p',   DLAMCH := eps*base
-*>          = 'N' or 'n',   DLAMCH := t
-*>          = 'R' or 'r',   DLAMCH := rnd
-*>          = 'M' or 'm',   DLAMCH := emin
-*>          = 'U' or 'u',   DLAMCH := rmin
-*>          = 'L' or 'l',   DLAMCH := emax
-*>          = 'O' or 'o',   DLAMCH := rmax
-*>          where
-*>          eps   = relative machine precision
-*>          sfmin = safe minimum, such that 1/sfmin does not overflow
-*>          base  = base of the machine
-*>          prec  = eps*base
-*>          t     = number of (base) digits in the mantissa
-*>          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
-*>          emin  = minimum exponent before (gradual) underflow
-*>          rmin  = underflow threshold - base**(emin-1)
-*>          emax  = largest exponent before overflow
-*>          rmax  = overflow threshold  - (base**emax)*(1-eps)
-*> \endverbatim
-*
-*  Authors:
-*  ========
-*
-*> \author Univ. of Tennessee 
-*> \author Univ. of California Berkeley 
-*> \author Univ. of Colorado Denver 
-*> \author NAG Ltd. 
-*
-*> \date November 2015
-*
-*> \ingroup auxOTHERauxiliary
-*
-*  =====================================================================
-      DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
-*
-*  -- LAPACK auxiliary routine (version 3.6.0) --
-*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
-*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2015
-*
-*     .. Scalar Arguments ..
-      CHARACTER          CMACH
-*     ..
-*
-* =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION   ONE, ZERO
-      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
-*     ..
-*     .. Local Scalars ..
-      DOUBLE PRECISION   RND, EPS, SFMIN, SMALL, RMACH
-*     ..
-*     .. External Functions ..
-      LOGICAL            LSAME
-      EXTERNAL           LSAME
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          DIGITS, EPSILON, HUGE, MAXEXPONENT,
-     $                   MINEXPONENT, RADIX, TINY
-*     ..
-*     .. Executable Statements ..
-*
-*
-*     Assume rounding, not chopping. Always.
-*
-      RND = ONE
-*
-      IF( ONE.EQ.RND ) THEN
-         EPS = EPSILON(ZERO) * 0.5
-      ELSE
-         EPS = EPSILON(ZERO)
-      END IF
-*
-      IF( LSAME( CMACH, 'E' ) ) THEN
-         RMACH = EPS
-      ELSE IF( LSAME( CMACH, 'S' ) ) THEN
-         SFMIN = TINY(ZERO)
-         SMALL = ONE / HUGE(ZERO)
-         IF( SMALL.GE.SFMIN ) THEN
-*
-*           Use SMALL plus a bit, to avoid the possibility of rounding
-*           causing overflow when computing  1/sfmin.
-*
-            SFMIN = SMALL*( ONE+EPS )
-         END IF
-         RMACH = SFMIN
-      ELSE IF( LSAME( CMACH, 'B' ) ) THEN
-         RMACH = RADIX(ZERO)
-      ELSE IF( LSAME( CMACH, 'P' ) ) THEN
-         RMACH = EPS * RADIX(ZERO)
-      ELSE IF( LSAME( CMACH, 'N' ) ) THEN
-         RMACH = DIGITS(ZERO)
-      ELSE IF( LSAME( CMACH, 'R' ) ) THEN
-         RMACH = RND
-      ELSE IF( LSAME( CMACH, 'M' ) ) THEN
-         RMACH = MINEXPONENT(ZERO)
-      ELSE IF( LSAME( CMACH, 'U' ) ) THEN
-         RMACH = tiny(zero)
-      ELSE IF( LSAME( CMACH, 'L' ) ) THEN
-         RMACH = MAXEXPONENT(ZERO)
-      ELSE IF( LSAME( CMACH, 'O' ) ) THEN
-         RMACH = HUGE(ZERO)
-      ELSE
-         RMACH = ZERO
-      END IF
-*
-      DLAMCH = RMACH
-      RETURN
-*
-*     End of DLAMCH
-*
-      END
-************************************************************************
-*> \brief \b DLAMC3
-*> \details
-*> \b Purpose:
-*> \verbatim
-*> DLAMC3  is intended to force  A  and  B  to be stored prior to doing
-*> the addition of  A  and  B ,  for use in situations where optimizers
-*> might hold one of these in a register.
-*> \endverbatim
-*> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
-*> \date November 2015
-*> \ingroup auxOTHERauxiliary
-*>
-*> \param[in] A
-*> \verbatim
-*>          A is a DOUBLE PRECISION
-*> \endverbatim
-*>
-*> \param[in] B
-*> \verbatim
-*>          B is a DOUBLE PRECISION
-*>          The values A and B.
-*> \endverbatim
-*>
-      DOUBLE PRECISION FUNCTION DLAMC3( A, B )
-*
-*  -- LAPACK auxiliary routine (version 3.6.0) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2010
-*
-*     .. Scalar Arguments ..
-      DOUBLE PRECISION   A, B
-*     ..
-* =====================================================================
-*
-*     .. Executable Statements ..
-*
-      DLAMC3 = A + B
-*
-      RETURN
-*
-*     End of DLAMC3
-*
-      END
-*
-************************************************************************

src/Fortran/dlapy2.f

-*> \brief \b DLAPY2 returns sqrt(x2+y2).
-*
-*  =========== DOCUMENTATION ===========
-*
-* Online html documentation available at 
-*            http://www.netlib.org/lapack/explore-html/ 
-*
-*> \htmlonly
-*> Download DLAPY2 + dependencies 
-*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlapy2.f"> 
-*> [TGZ]</a> 
-*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlapy2.f"> 
-*> [ZIP]</a> 
-*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlapy2.f"> 
-*> [TXT]</a>
-*> \endhtmlonly 
-*
-*  Definition:
-*  ===========
-*
-*       DOUBLE PRECISION FUNCTION DLAPY2( X, Y )
-* 
-*       .. Scalar Arguments ..
-*       DOUBLE PRECISION   X, Y
-*       ..
-*  
-*
-*> \par Purpose:
-*  =============
-*>
-*> \verbatim
-*>
-*> DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary
-*> overflow.
-*> \endverbatim
-*
-*  Arguments:
-*  ==========
-*
-*> \param[in] X
-*> \verbatim
-*>          X is DOUBLE PRECISION
-*> \endverbatim
-*>
-*> \param[in] Y
-*> \verbatim
-*>          Y is DOUBLE PRECISION
-*>          X and Y specify the values x and y.
-*> \endverbatim
-*
-*  Authors:
-*  ========
-*
-*> \author Univ. of Tennessee 
-*> \author Univ. of California Berkeley 
-*> \author Univ. of Colorado Denver 
-*> \author NAG Ltd. 
-*
-*> \date September 2012
-*
-*> \ingroup auxOTHERauxiliary
-*
-*  =====================================================================
-      DOUBLE PRECISION FUNCTION DLAPY2( X, Y )
-*
-*  -- LAPACK auxiliary routine (version 3.4.2) --
-*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
-*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     September 2012
-*
-*     .. Scalar Arguments ..
-      DOUBLE PRECISION   X, Y
-*     ..
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION   ZERO
-      PARAMETER          ( ZERO = 0.0D0 )
-      DOUBLE PRECISION   ONE
-      PARAMETER          ( ONE = 1.0D0 )
-*     ..
-*     .. Local Scalars ..
-      DOUBLE PRECISION   W, XABS, YABS, Z
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          ABS, MAX, MIN, SQRT
-*     ..
-*     .. Executable Statements ..
-*
-      XABS = ABS( X )
-      YABS = ABS( Y )
-      W = MAX( XABS, YABS )
-      Z = MIN( XABS, YABS )
-      IF( Z.EQ.ZERO ) THEN
-         DLAPY2 = W
-      ELSE
-         DLAPY2 = W*SQRT( ONE+( Z / W )**2 )
-      END IF
-      RETURN
-*
-*     End of DLAPY2
-*
-      END