LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 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
logical function lsamen(N, CA, CB)
LSAMEN
Definition: lsamen.f:74
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zlahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO, PATH)
ZLAHILB
Definition: zlahilb.f:134
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106