1*> \brief \b ZLAHILB
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 ZLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
12*            INFO, PATH)
13*
14*       .. Scalar Arguments ..
15*       INTEGER N, NRHS, LDA, LDX, LDB, INFO
16*       .. Array Arguments ..
17*       DOUBLE PRECISION WORK(N)
18*       COMPLEX*16 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*> ZLAHILB 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 INTEGER
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*> \ingroup complex16_lin
130*
131*  =====================================================================
132      SUBROUTINE ZLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
133     $     INFO, PATH)
134*
135*  -- LAPACK test routine --
136*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
137*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138*
139*     .. Scalar Arguments ..
140      INTEGER N, NRHS, LDA, LDX, LDB, INFO
141*     .. Array Arguments ..
142      DOUBLE PRECISION WORK(N)
143      COMPLEX*16 A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
144      CHARACTER*3 PATH
145*     ..
146*
147*  =====================================================================
148*     .. Local Scalars ..
149      INTEGER TM, TI, R
150      INTEGER M
151      INTEGER I, J
152      COMPLEX*16 TMP
153      CHARACTER*2 C2
154*     ..
155*     .. Parameters ..
156*     NMAX_EXACT   the largest dimension where the generated data is
157*                  exact.
158*     NMAX_APPROX  the largest dimension where the generated data has
159*                  a small componentwise relative error.
160*     ??? complex uses how many bits ???
161      INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D
162      PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11, SIZE_D = 8)
163*
164*     d's are generated from random permutation of those eight elements.
165      COMPLEX*16 d1(8), d2(8), invd1(8), invd2(8)
166      DATA D1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/
167      DATA D2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/
168
169      DATA INVD1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0),
170     $     (-.5,-.5),(.5,-.5),(.5,.5)/
171      DATA INVD2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0),
172     $     (-.5,.5),(.5,.5),(.5,-.5)/
173*     ..
174*     .. External Functions
175      EXTERNAL ZLASET, LSAMEN
176      INTRINSIC DBLE
177      LOGICAL LSAMEN
178*     ..
179*     .. Executable Statements ..
180      C2 = PATH( 2: 3 )
181*
182*     Test the input arguments
183*
184      INFO = 0
185      IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
186         INFO = -1
187      ELSE IF (NRHS .LT. 0) THEN
188         INFO = -2
189      ELSE IF (LDA .LT. N) THEN
190         INFO = -4
191      ELSE IF (LDX .LT. N) THEN
192         INFO = -6
193      ELSE IF (LDB .LT. N) THEN
194         INFO = -8
195      END IF
196      IF (INFO .LT. 0) THEN
197         CALL XERBLA('ZLAHILB', -INFO)
198         RETURN
199      END IF
200      IF (N .GT. NMAX_EXACT) THEN
201         INFO = 1
202      END IF
203*
204*     Compute M = the LCM of the integers [1, 2*N-1].  The largest
205*     reasonable N is small enough that integers suffice (up to N = 11).
206      M = 1
207      DO I = 2, (2*N-1)
208         TM = M
209         TI = I
210         R = MOD(TM, TI)
211         DO WHILE (R .NE. 0)
212            TM = TI
213            TI = R
214            R = MOD(TM, TI)
215         END DO
216         M = (M / TI) * I
217      END DO
218*
219*     Generate the scaled Hilbert matrix in A
220*     If we are testing SY routines,
221*        take D1_i = D2_i, else, D1_i = D2_i*
222      IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
223         DO J = 1, N
224            DO I = 1, N
225               A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1))
226     $              * D1(MOD(I,SIZE_D)+1)
227            END DO
228         END DO
229      ELSE
230         DO J = 1, N
231            DO I = 1, N
232               A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1))
233     $              * D2(MOD(I,SIZE_D)+1)
234            END DO
235         END DO
236      END IF
237*
238*     Generate matrix B as simply the first NRHS columns of M * the
239*     identity.
240      TMP = DBLE(M)
241      CALL ZLASET('Full', N, NRHS, (0.0D+0,0.0D+0), TMP, B, LDB)
242*
243*     Generate the true solutions in X.  Because B = the first NRHS
244*     columns of M*I, the true solutions are just the first NRHS columns
245*     of the inverse Hilbert matrix.
246      WORK(1) = N
247      DO J = 2, N
248         WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
249     $        * (N +J -1)
250      END DO
251
252*     If we are testing SY routines,
253*           take D1_i = D2_i, else, D1_i = D2_i*
254      IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
255         DO J = 1, NRHS
256            DO I = 1, N
257               X(I, J) = INVD1(MOD(J,SIZE_D)+1) *
258     $              ((WORK(I)*WORK(J)) / (I + J - 1))
259     $              * INVD1(MOD(I,SIZE_D)+1)
260            END DO
261         END DO
262      ELSE
263         DO J = 1, NRHS
264            DO I = 1, N
265               X(I, J) = INVD2(MOD(J,SIZE_D)+1) *
266     $              ((WORK(I)*WORK(J)) / (I + J - 1))
267     $              * INVD1(MOD(I,SIZE_D)+1)
268            END DO
269         END DO
270      END IF
271      END
272