TESTING/LIN/clahilb.f

   1 *> \brief \b CLAHILB
   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 CLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
  12 *            INFO, PATH)
  13 *
  14 *       .. Scalar Arguments ..
  15 *       INTEGER T, N, NRHS, LDA, LDX, LDB, INFO
  16 *       .. Array Arguments ..
  17 *       REAL WORK(N)
  18 *       COMPLEX A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
  19 *       CHARACTER*3        PATH
  20 *       ..
  21 *
  22 *
  23 *> \par Purpose:
  24 *  =============
  25 *>
  26 *> \verbatim
  27 *>
  28 *> CLAHILB generates an N by N scaled Hilbert matrix in A along with
  29 *> NRHS right-hand sides in B and solutions in X such that A*X=B.
  30 *>
  31 *> The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
  32 *> entries are integers.  The right-hand sides are the first NRHS
  33 *> columns of M * the identity matrix, and the solutions are the
  34 *> first NRHS columns of the inverse Hilbert matrix.
  35 *>
  36 *> The condition number of the Hilbert matrix grows exponentially with
  37 *> its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
  38 *> Hilbert matrices beyond a relatively small dimension cannot be
  39 *> generated exactly without extra precision.  Precision is exhausted
  40 *> when the largest entry in the inverse Hilbert matrix is greater than
  41 *> 2 to the power of the number of bits in the fraction of the data type
  42 *> used plus one, which is 24 for single precision.
  43 *>
  44 *> In single, the generated solution is exact for N <= 6 and has
  45 *> small componentwise error for 7 <= N <= 11.
  46 *> \endverbatim
  47 *
  48 *  Arguments:
  49 *  ==========
  50 *
  51 *> \param[in] N
  52 *> \verbatim
  53 *>          N is INTEGER
  54 *>          The dimension of the matrix A.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] NRHS
  58 *> \verbatim
  59 *>          NRHS is NRHS
  60 *>          The requested number of right-hand sides.
  61 *> \endverbatim
  62 *>
  63 *> \param[out] A
  64 *> \verbatim
  65 *>          A is COMPLEX array, dimension (LDA, N)
  66 *>          The generated scaled Hilbert matrix.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] LDA
  70 *> \verbatim
  71 *>          LDA is INTEGER
  72 *>          The leading dimension of the array A.  LDA >= N.
  73 *> \endverbatim
  74 *>
  75 *> \param[out] X
  76 *> \verbatim
  77 *>          X is COMPLEX array, dimension (LDX, NRHS)
  78 *>          The generated exact solutions.  Currently, the first NRHS
  79 *>          columns of the inverse Hilbert matrix.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] LDX
  83 *> \verbatim
  84 *>          LDX is INTEGER
  85 *>          The leading dimension of the array X.  LDX >= N.
  86 *> \endverbatim
  87 *>
  88 *> \param[out] B
  89 *> \verbatim
  90 *>          B is REAL array, dimension (LDB, NRHS)
  91 *>          The generated right-hand sides.  Currently, the first NRHS
  92 *>          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] LDB
  96 *> \verbatim
  97 *>          LDB is INTEGER
  98 *>          The leading dimension of the array B.  LDB >= N.
  99 *> \endverbatim
 100 *>
 101 *> \param[out] WORK
 102 *> \verbatim
 103 *>          WORK is REAL array, dimension (N)
 104 *> \endverbatim
 105 *>
 106 *> \param[out] INFO
 107 *> \verbatim
 108 *>          INFO is INTEGER
 109 *>          = 0: successful exit
 110 *>          = 1: N is too large; the data is still generated but may not
 111 *>               be not exact.
 112 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 113 *> \endverbatim
 114 *>
 115 *> \param[in] PATH
 116 *> \verbatim
 117 *>          PATH is CHARACTER*3
 118 *>          The LAPACK path name.
 119 *> \endverbatim
 120 *
 121 *  Authors:
 122 *  ========
 123 *
 124 *> \author Univ. of Tennessee
 125 *> \author Univ. of California Berkeley
 126 *> \author Univ. of Colorado Denver
 127 *> \author NAG Ltd.
 128 *
 129 *> \date December 2016
 130 *
 131 *> \ingroup complex_lin
 132 *
 133 *  =====================================================================
 134       SUBROUTINE CLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
 135      $     INFO, PATH)
 136 *
 137 *  -- LAPACK test routine (version 3.7.0) --
 138 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 139 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 140 *     December 2016
 141 *
 142 *     .. Scalar Arguments ..
 143       INTEGER T, N, NRHS, LDA, LDX, LDB, INFO
 144 *     .. Array Arguments ..
 145       REAL WORK(N)
 146       COMPLEX A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
 147       CHARACTER*3        PATH
 148 *     ..
 149 *
 150 *  =====================================================================
 151 *     .. Local Scalars ..
 152       INTEGER TM, TI, R
 153       INTEGER M
 154       INTEGER I, J
 155       COMPLEX TMP
 156       CHARACTER*2 C2
 157 *     ..
 158 *     .. Parameters ..
 159 *     NMAX_EXACT   the largest dimension where the generated data is
 160 *                  exact.
 161 *     NMAX_APPROX  the largest dimension where the generated data has
 162 *                  a small componentwise relative error.
 163 *     ??? complex uses how many bits ???
 164       INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D
 165       PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11, SIZE_D = 8)
 166 *
 167 *     d's are generated from random permuation of those eight elements.
 168       COMPLEX D1(8), D2(8), INVD1(8), INVD2(8)
 169       DATA D1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/
 170       DATA D2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/
 171
 172       DATA INVD1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0),
 173      $     (-.5,-.5),(.5,-.5),(.5,.5)/
 174       DATA INVD2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0),
 175      $     (-.5,.5),(.5,.5),(.5,-.5)/
 176 *     ..
 177 *     .. External Functions
 178       EXTERNAL CLASET, LSAMEN
 179       INTRINSIC REAL
 180       LOGICAL LSAMEN
 181 *     ..
 182 *     .. Executable Statements ..
 183       C2 = PATH( 2: 3 )
 184 *
 185 *     Test the input arguments
 186 *
 187       INFO = 0
 188       IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
 189          INFO = -1
 190       ELSE IF (NRHS .LT. 0) THEN
 191          INFO = -2
 192       ELSE IF (LDA .LT. N) THEN
 193          INFO = -4
 194       ELSE IF (LDX .LT. N) THEN
 195          INFO = -6
 196       ELSE IF (LDB .LT. N) THEN
 197          INFO = -8
 198       END IF
 199       IF (INFO .LT. 0) THEN
 200          CALL XERBLA('CLAHILB', -INFO)
 201          RETURN
 202       END IF
 203       IF (N .GT. NMAX_EXACT) THEN
 204          INFO = 1
 205       END IF
 206 *
 207 *     Compute M = the LCM of the integers [1, 2*N-1].  The largest
 208 *     reasonable N is small enough that integers suffice (up to N = 11).
 209       M = 1
 210       DO I = 2, (2*N-1)
 211          TM = M
 212          TI = I
 213          R = MOD(TM, TI)
 214          DO WHILE (R .NE. 0)
 215             TM = TI
 216             TI = R
 217             R = MOD(TM, TI)
 218          END DO
 219          M = (M / TI) * I
 220       END DO
 221 *
 222 *     Generate the scaled Hilbert matrix in A
 223 *     If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i*
 224       IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
 225          DO J = 1, N
 226             DO I = 1, N
 227                A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1))
 228      $              * D1(MOD(I,SIZE_D)+1)
 229             END DO
 230          END DO
 231       ELSE
 232          DO J = 1, N
 233             DO I = 1, N
 234                A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1))
 235      $              * D2(MOD(I,SIZE_D)+1)
 236             END DO
 237          END DO
 238       END IF
 239 *
 240 *     Generate matrix B as simply the first NRHS columns of M * the
 241 *     identity.
 242       TMP = REAL(M)
 243       CALL CLASET('Full', N, NRHS, (0.0,0.0), TMP, B, LDB)
 244 *
 245 *     Generate the true solutions in X.  Because B = the first NRHS
 246 *     columns of M*I, the true solutions are just the first NRHS columns
 247 *     of the inverse Hilbert matrix.
 248       WORK(1) = N
 249       DO J = 2, N
 250          WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
 251      $        * (N +J -1)
 252       END DO
 253 *
 254 *     If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i*
 255       IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
 256          DO J = 1, NRHS
 257             DO I = 1, N
 258                X(I, J) =
 259      $              INVD1(MOD(J,SIZE_D)+1) *
 260      $              ((WORK(I)*WORK(J)) / (I + J - 1))
 261      $              * INVD1(MOD(I,SIZE_D)+1)
 262             END DO
 263          END DO
 264       ELSE
 265          DO J = 1, NRHS
 266             DO I = 1, N
 267                X(I, J) =
 268      $              INVD2(MOD(J,SIZE_D)+1) *
 269      $              ((WORK(I)*WORK(J)) / (I + J - 1))
 270      $              * INVD1(MOD(I,SIZE_D)+1)
 271             END DO
 272          END DO
 273       END IF
 274       END
 275