LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
clahilb.f
Go to the documentation of this file.
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 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 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 *> \date November 2011
130 *
131 *> \ingroup complex_matgen
132 *
133 * =====================================================================
134  SUBROUTINE clahilb( N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
135  $ info, path)
136 *
137 * -- LAPACK auxiliary 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  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 * ..
178 * .. External Functions
179  EXTERNAL claset, lsamen
180  INTRINSIC real
181  LOGICAL lsamen
182 * ..
183 * .. Executable Statements ..
184  c2 = path( 2: 3 )
185 *
186 * Test the input arguments
187 *
188  info = 0
189  IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
190  info = -1
191  ELSE IF (nrhs .LT. 0) THEN
192  info = -2
193  ELSE IF (lda .LT. n) THEN
194  info = -4
195  ELSE IF (ldx .LT. n) THEN
196  info = -6
197  ELSE IF (ldb .LT. n) THEN
198  info = -8
199  END IF
200  IF (info .LT. 0) THEN
201  CALL xerbla('CLAHILB', -info)
202  return
203  END IF
204  IF (n .GT. nmax_exact) THEN
205  info = 1
206  END IF
207 
208 * Compute M = the LCM of the integers [1, 2*N-1]. The largest
209 * reasonable N is small enough that integers suffice (up to N = 11).
210  m = 1
211  DO i = 2, (2*n-1)
212  tm = m
213  ti = i
214  r = mod(tm, ti)
215  DO WHILE (r .NE. 0)
216  tm = ti
217  ti = r
218  r = mod(tm, ti)
219  END DO
220  m = (m / ti) * i
221  END DO
222 
223 * Generate the scaled Hilbert matrix in A
224 * If we are testing SY routines, take
225 * D1_i = D2_i, else, D1_i = D2_i*
226  IF ( lsamen( 2, c2, 'SY' ) ) THEN
227  DO j = 1, n
228  DO i = 1, n
229  a(i, j) = d1(mod(j,size_d)+1) * (REAL(M) / (i + j - 1))
230  $ * d1(mod(i,size_d)+1)
231  END DO
232  END DO
233  ELSE
234  DO j = 1, n
235  DO i = 1, n
236  a(i, j) = d1(mod(j,size_d)+1) * (REAL(M) / (i + j - 1))
237  $ * d2(mod(i,size_d)+1)
238  END DO
239  END DO
240  END IF
241 
242 * Generate matrix B as simply the first NRHS columns of M * the
243 * identity.
244  tmp = REAL(m)
245  CALL claset('Full', n, nrhs, (0.0,0.0), tmp, b, ldb)
246 
247 * Generate the true solutions in X. Because B = the first NRHS
248 * columns of M*I, the true solutions are just the first NRHS columns
249 * of the inverse Hilbert matrix.
250  work(1) = n
251  DO j = 2, n
252  work(j) = ( ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1) )
253  $ * (n +j -1)
254  END DO
255 
256 * If we are testing SY routines,
257 * take D1_i = D2_i, else, D1_i = D2_i*
258  IF ( lsamen( 2, c2, 'SY' ) ) THEN
259  DO j = 1, nrhs
260  DO i = 1, n
261  x(i, j) =
262  $ invd1(mod(j,size_d)+1) *
263  $ ((work(i)*work(j)) / (i + j - 1))
264  $ * invd1(mod(i,size_d)+1)
265  END DO
266  END DO
267  ELSE
268  DO j = 1, nrhs
269  DO i = 1, n
270  x(i, j) =
271  $ invd2(mod(j,size_d)+1) *
272  $ ((work(i)*work(j)) / (i + j - 1))
273  $ * invd1(mod(i,size_d)+1)
274  END DO
275  END DO
276  END IF
277  END
278