TESTING/MATGEN/dlahilb.f

   1 C> \brief \b DLAHILB
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE DLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
  12 *
  13 *       .. Scalar Arguments ..
  14 *       INTEGER N, NRHS, LDA, LDX, LDB, INFO
  15 *       .. Array Arguments ..
  16 *       DOUBLE PRECISION A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
  17 *       ..
  18 *
  19 *
  20 *> \par Purpose:
  21 *  =============
  22 *>
  23 *> \verbatim
  24 *>
  25 *> DLAHILB generates an N by N scaled Hilbert matrix in A along with
  26 *> NRHS right-hand sides in B and solutions in X such that A*X=B.
  27 *>
  28 *> The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
  29 *> entries are integers.  The right-hand sides are the first NRHS
  30 *> columns of M * the identity matrix, and the solutions are the
  31 *> first NRHS columns of the inverse Hilbert matrix.
  32 *>
  33 *> The condition number of the Hilbert matrix grows exponentially with
  34 *> its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
  35 *> Hilbert matrices beyond a relatively small dimension cannot be
  36 *> generated exactly without extra precision.  Precision is exhausted
  37 *> when the largest entry in the inverse Hilbert matrix is greater than
  38 *> 2 to the power of the number of bits in the fraction of the data type
  39 *> used plus one, which is 24 for single precision.
  40 *>
  41 *> In single, the generated solution is exact for N <= 6 and has
  42 *> small componentwise error for 7 <= N <= 11.
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in] N
  49 *> \verbatim
  50 *>          N is INTEGER
  51 *>          The dimension of the matrix A.
  52 *> \endverbatim
  53 *>
  54 *> \param[in] NRHS
  55 *> \verbatim
  56 *>          NRHS is INTEGER
  57 *>          The requested number of right-hand sides.
  58 *> \endverbatim
  59 *>
  60 *> \param[out] A
  61 *> \verbatim
  62 *>          A is DOUBLE PRECISION array, dimension (LDA, N)
  63 *>          The generated scaled Hilbert matrix.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] LDA
  67 *> \verbatim
  68 *>          LDA is INTEGER
  69 *>          The leading dimension of the array A.  LDA >= N.
  70 *> \endverbatim
  71 *>
  72 *> \param[out] X
  73 *> \verbatim
  74 *>          X is DOUBLE PRECISION array, dimension (LDX, NRHS)
  75 *>          The generated exact solutions.  Currently, the first NRHS
  76 *>          columns of the inverse Hilbert matrix.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] LDX
  80 *> \verbatim
  81 *>          LDX is INTEGER
  82 *>          The leading dimension of the array X.  LDX >= N.
  83 *> \endverbatim
  84 *>
  85 *> \param[out] B
  86 *> \verbatim
  87 *>          B is DOUBLE PRECISION array, dimension (LDB, NRHS)
  88 *>          The generated right-hand sides.  Currently, the first NRHS
  89 *>          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] LDB
  93 *> \verbatim
  94 *>          LDB is INTEGER
  95 *>          The leading dimension of the array B.  LDB >= N.
  96 *> \endverbatim
  97 *>
  98 *> \param[out] WORK
  99 *> \verbatim
 100 *>          WORK is DOUBLE PRECISION array, dimension (N)
 101 *> \endverbatim
 102 *>
 103 *> \param[out] INFO
 104 *> \verbatim
 105 *>          INFO is INTEGER
 106 *>          = 0: successful exit
 107 *>          = 1: N is too large; the data is still generated but may not
 108 *>               be not exact.
 109 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 110 *> \endverbatim
 111 *
 112 *  Authors:
 113 *  ========
 114 *
 115 *> \author Univ. of Tennessee
 116 *> \author Univ. of California Berkeley
 117 *> \author Univ. of Colorado Denver
 118 *> \author NAG Ltd.
 119 *
 120 *> \date November 2015
 121 *
 122 *> \ingroup double_matgen
 123 *
 124 *  =====================================================================
 125       SUBROUTINE DLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
 126 *
 127 *  -- LAPACK test routine (version 3.6.0) --
 128 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 129 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 130 *     November 2015
 131 *
 132 *     .. Scalar Arguments ..
 133       INTEGER N, NRHS, LDA, LDX, LDB, INFO
 134 *     .. Array Arguments ..
 135       DOUBLE PRECISION A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
 136 *     ..
 137 *
 138 *  =====================================================================
 139 *     .. Local Scalars ..
 140       INTEGER TM, TI, R
 141       INTEGER M
 142       INTEGER I, J
 143
 144 *     .. Parameters ..
 145 *     NMAX_EXACT   the largest dimension where the generated data is
 146 *                  exact.
 147 *     NMAX_APPROX  the largest dimension where the generated data has
 148 *                  a small componentwise relative error.
 149       INTEGER NMAX_EXACT, NMAX_APPROX
 150       PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11)
 151
 152 *     ..
 153 *     .. External Functions
 154       EXTERNAL DLASET
 155       INTRINSIC DBLE
 156 *     ..
 157 *     .. Executable Statements ..
 158 *
 159 *     Test the input arguments
 160 *
 161       INFO = 0
 162       IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
 163          INFO = -1
 164       ELSE IF (NRHS .LT. 0) THEN
 165          INFO = -2
 166       ELSE IF (LDA .LT. N) THEN
 167          INFO = -4
 168       ELSE IF (LDX .LT. N) THEN
 169          INFO = -6
 170       ELSE IF (LDB .LT. N) THEN
 171          INFO = -8
 172       END IF
 173       IF (INFO .LT. 0) THEN
 174          CALL XERBLA('DLAHILB', -INFO)
 175          RETURN
 176       END IF
 177       IF (N .GT. NMAX_EXACT) THEN
 178          INFO = 1
 179       END IF
 180
 181 *     Compute M = the LCM of the integers [1, 2*N-1].  The largest
 182 *     reasonable N is small enough that integers suffice (up to N = 11).
 183       M = 1
 184       DO I = 2, (2*N-1)
 185          TM = M
 186          TI = I
 187          R = MOD(TM, TI)
 188          DO WHILE (R .NE. 0)
 189             TM = TI
 190             TI = R
 191             R = MOD(TM, TI)
 192          END DO
 193          M = (M / TI) * I
 194       END DO
 195
 196 *     Generate the scaled Hilbert matrix in A
 197       DO J = 1, N
 198          DO I = 1, N
 199             A(I, J) = DBLE(M) / (I + J - 1)
 200          END DO
 201       END DO
 202
 203 *     Generate matrix B as simply the first NRHS columns of M * the
 204 *     identity.
 205       CALL DLASET('Full', N, NRHS, 0.0D+0, DBLE(M), B, LDB)
 206
 207 *     Generate the true solutions in X.  Because B = the first NRHS
 208 *     columns of M*I, the true solutions are just the first NRHS columns
 209 *     of the inverse Hilbert matrix.
 210       WORK(1) = N
 211       DO J = 2, N
 212          WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
 213      $        * (N +J -1)
 214       END DO
 215
 216       DO J = 1, NRHS
 217          DO I = 1, N
 218             X(I, J) = (WORK(I)*WORK(J)) / (I + J - 1)
 219          END DO
 220       END DO
 221
 222       END
 223