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 November 2011
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.4.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*     November 2011
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
276