1*> \brief \b SLAHILB
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 SLAHILB( 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*       REAL A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
17*       ..
18*
19*
20*> \par Purpose:
21*  =============
22*>
23*> \verbatim
24*>
25*> SLAHILB 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 REAL 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 REAL 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 REAL 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 REAL 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*> \ingroup single_lin
121*
122*  =====================================================================
123      SUBROUTINE SLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
124*
125*  -- LAPACK test routine --
126*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
127*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
128*
129*     .. Scalar Arguments ..
130      INTEGER N, NRHS, LDA, LDX, LDB, INFO
131*     .. Array Arguments ..
132      REAL A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
133*     ..
134*
135*  =====================================================================
136*     .. Local Scalars ..
137      INTEGER TM, TI, R
138      INTEGER M
139      INTEGER I, J
140*     ..
141*     .. Parameters ..
142*     NMAX_EXACT   the largest dimension where the generated data is
143*                  exact.
144*     NMAX_APPROX  the largest dimension where the generated data has
145*                  a small componentwise relative error.
146      INTEGER NMAX_EXACT, NMAX_APPROX
147      PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11)
148*     ..
149*     .. External Functions
150      EXTERNAL SLASET
151      INTRINSIC REAL
152*     ..
153*     .. Executable Statements ..
154*
155*     Test the input arguments
156*
157      INFO = 0
158      IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
159         INFO = -1
160      ELSE IF (NRHS .LT. 0) THEN
161         INFO = -2
162      ELSE IF (LDA .LT. N) THEN
163         INFO = -4
164      ELSE IF (LDX .LT. N) THEN
165         INFO = -6
166      ELSE IF (LDB .LT. N) THEN
167         INFO = -8
168      END IF
169      IF (INFO .LT. 0) THEN
170         CALL XERBLA('SLAHILB', -INFO)
171         RETURN
172      END IF
173      IF (N .GT. NMAX_EXACT) THEN
174         INFO = 1
175      END IF
176*
177*     Compute M = the LCM of the integers [1, 2*N-1].  The largest
178*     reasonable N is small enough that integers suffice (up to N = 11).
179      M = 1
180      DO I = 2, (2*N-1)
181         TM = M
182         TI = I
183         R = MOD(TM, TI)
184         DO WHILE (R .NE. 0)
185            TM = TI
186            TI = R
187            R = MOD(TM, TI)
188         END DO
189         M = (M / TI) * I
190      END DO
191*
192*     Generate the scaled Hilbert matrix in A
193      DO J = 1, N
194         DO I = 1, N
195            A(I, J) = REAL(M) / (I + J - 1)
196         END DO
197      END DO
198*
199*     Generate matrix B as simply the first NRHS columns of M * the
200*     identity.
201      CALL SLASET('Full', N, NRHS, 0.0, REAL(M), B, LDB)
202*
203*     Generate the true solutions in X.  Because B = the first NRHS
204*     columns of M*I, the true solutions are just the first NRHS columns
205*     of the inverse Hilbert matrix.
206      WORK(1) = N
207      DO J = 2, N
208         WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
209     $        * (N +J -1)
210      END DO
211*
212      DO J = 1, NRHS
213         DO I = 1, N
214            X(I, J) = (WORK(I)*WORK(J)) / (I + J - 1)
215         END DO
216      END DO
217*
218      END
219
220