LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zlahilb.f
Go to the documentation of this file.
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 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 complex16_lin
132 *
133 * =====================================================================
134  SUBROUTINE zlahilb(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 n, nrhs, lda, ldx, ldb, info
144 * .. Array Arguments ..
145  DOUBLE PRECISION work(n)
146  COMPLEX*16 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*16 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*16 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 zlaset, lsamen
179  INTRINSIC dble
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('ZLAHILB', -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) * (dble(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) * (dble(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 = dble(m)
243  CALL zlaset('Full', n, nrhs, (0.0d+0,0.0d+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) = invd1(mod(j,size_d)+1) *
259  $ ((work(i)*work(j)) / (i + j - 1))
260  $ * invd1(mod(i,size_d)+1)
261  END DO
262  END DO
263  ELSE
264  DO j = 1, nrhs
265  DO i = 1, n
266  x(i, j) = invd2(mod(j,size_d)+1) *
267  $ ((work(i)*work(j)) / (i + j - 1))
268  $ * invd1(mod(i,size_d)+1)
269  END DO
270  END DO
271  END IF
272  END