dlahilb.f

Go to the documentation of this file.

*> \brief \b DLAHILB
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
* 
*       .. Scalar Arguments ..
*       INTEGER N, NRHS, LDA, LDX, LDB, INFO
*       .. Array Arguments ..
*       DOUBLE PRECISION A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLAHILB generates an N by N scaled Hilbert matrix in A along with
*> NRHS right-hand sides in B and solutions in X such that A*X=B.
*>
*> The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
*> entries are integers.  The right-hand sides are the first NRHS 
*> columns of M * the identity matrix, and the solutions are the 
*> first NRHS columns of the inverse Hilbert matrix.
*>
*> The condition number of the Hilbert matrix grows exponentially with
*> its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
*> Hilbert matrices beyond a relatively small dimension cannot be
*> generated exactly without extra precision.  Precision is exhausted
*> when the largest entry in the inverse Hilbert matrix is greater than
*> 2 to the power of the number of bits in the fraction of the data type
*> used plus one, which is 24 for single precision.  
*>
*> In single, the generated solution is exact for N <= 6 and has
*> small componentwise error for 7 <= N <= 11.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The dimension of the matrix A.
*> \endverbatim
*>      
*> \param[in] NRHS
*> \verbatim
*>          NRHS is NRHS
*>          The requested number of right-hand sides.
*> \endverbatim
*>
*> \param[out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA, N)
*>          The generated scaled Hilbert matrix.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= N.
*> \endverbatim
*>
*> \param[out] X
*> \verbatim
*>          X is DOUBLE PRECISION array, dimension (LDX, NRHS)
*>          The generated exact solutions.  Currently, the first NRHS
*>          columns of the inverse Hilbert matrix.
*> \endverbatim
*>
*> \param[in] LDX
*> \verbatim
*>          LDX is INTEGER
*>          The leading dimension of the array X.  LDX >= N.
*> \endverbatim
*>
*> \param[out] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB, NRHS)
*>          The generated right-hand sides.  Currently, the first NRHS
*>          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B.  LDB >= N.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0: successful exit
*>          = 1: N is too large; the data is still generated but may not
*>               be not exact.
*>          < 0: if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup double_lin
*
*  =====================================================================
      SUBROUTINE dlahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER N, NRHS, LDA, LDX, LDB, INFO
*     .. Array Arguments ..
      DOUBLE PRECISION A(lda, n), X(ldx, nrhs), B(ldb, nrhs), WORK(n)
*     ..
*
*  =====================================================================
*     .. Local Scalars ..
      INTEGER TM, TI, R
      INTEGER M
      INTEGER I, J
      COMPLEX*16 TMP
*     ..
*     .. Parameters ..
*     NMAX_EXACT   the largest dimension where the generated data is
*                  exact.
*     NMAX_APPROX  the largest dimension where the generated data has
*                  a small componentwise relative error.
      INTEGER NMAX_EXACT, NMAX_APPROX
      parameter(nmax_exact = 6, nmax_approx = 11)

*     ..
*     .. External Functions
      EXTERNAL dlaset
      INTRINSIC dble
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments
*
      info = 0
      IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
         info = -1
      ELSE IF (nrhs .LT. 0) THEN
         info = -2
      ELSE IF (lda .LT. n) THEN
         info = -4
      ELSE IF (ldx .LT. n) THEN
         info = -6
      ELSE IF (ldb .LT. n) THEN
         info = -8
      END IF
      IF (info .LT. 0) THEN
         CALL xerbla('DLAHILB', -info)
         RETURN
      END IF
      IF (n .GT. nmax_exact) THEN
         info = 1
      END IF
*
*     Compute M = the LCM of the integers [1, 2*N-1].  The largest
*     reasonable N is small enough that integers suffice (up to N = 11).
      m = 1
      DO i = 2, (2*n-1)
         tm = m
         ti = i
         r = mod(tm, ti)
         DO WHILE (r .NE. 0)
            tm = ti
            ti = r
            r = mod(tm, ti)
         END DO
         m = (m / ti) * i
      END DO
*
*     Generate the scaled Hilbert matrix in A
      DO j = 1, n
         DO i = 1, n
            a(i, j) = dble(m) / (i + j - 1)
         END DO
      END DO
*
*     Generate matrix B as simply the first NRHS columns of M * the
*     identity.
      tmp = dble(m)
      CALL dlaset('Full', n, nrhs, 0.0d+0, tmp, b, ldb)
*
*     Generate the true solutions in X.  Because B = the first NRHS
*     columns of M*I, the true solutions are just the first NRHS columns
*     of the inverse Hilbert matrix.
      work(1) = n
      DO j = 2, n
         work(j) = (  ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1)  )
     $        * (n +j -1)
      END DO
*
      DO j = 1, nrhs
         DO i = 1, n
            x(i, j) = (work(i)*work(j)) / (i + j - 1)
         END DO
      END DO
*
      END