LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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_matgen
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.0d0,0.0d0),(0.0d0,1.0d0),(-1.0d0,-1.0d0),
167 $ (0.0d0,-1.0d0),(1.0d0,0.0d0),(-1.0d0,1.0d0),(1.0d0,1.0d0),
168 $ (1.0d0,-1.0d0)/
169 DATA d2 /(-1.0d0,0.0d0),(0.0d0,-1.0d0),(-1.0d0,1.0d0),
170 $ (0.0d0,1.0d0),(1.0d0,0.0d0),(-1.0d0,-1.0d0),(1.0d0,-1.0d0),
171 $ (1.0d0,1.0d0)/
172
173 DATA invd1 /(-1.0d0,0.0d0),(0.0d0,-1.0d0),(-0.5d0,0.5d0),
174 $ (0.0d0,1.0d0),(1.0d0,0.0d0),(-0.5d0,-0.5d0),(0.5d0,-0.5d0),
175 $ (0.5d0,0.5d0)/
176 DATA invd2 /(-1.0d0,0.0d0),(0.0d0,1.0d0),(-0.5d0,-0.5d0),
177 $ (0.0d0,-1.0d0),(1.0d0,0.0d0),(-0.5d0,0.5d0),(0.5d0,0.5d0),
178 $ (0.5d0,-0.5d0)/
179* ..
180* .. External Subroutines ..
181 EXTERNAL xerbla
182* ..
183* .. External Functions
184 EXTERNAL zlaset, lsamen
185 INTRINSIC dble
186 LOGICAL LSAMEN
187* ..
188* .. Executable Statements ..
189 c2 = path( 2: 3 )
190*
191* Test the input arguments
192*
193 info = 0
194 IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
195 info = -1
196 ELSE IF (nrhs .LT. 0) THEN
197 info = -2
198 ELSE IF (lda .LT. n) THEN
199 info = -4
200 ELSE IF (ldx .LT. n) THEN
201 info = -6
202 ELSE IF (ldb .LT. n) THEN
203 info = -8
204 END IF
205 IF (info .LT. 0) THEN
206 CALL xerbla('ZLAHILB', -info)
207 RETURN
208 END IF
209 IF (n .GT. nmax_exact) THEN
210 info = 1
211 END IF
212*
213* Compute M = the LCM of the integers [1, 2*N-1]. The largest
214* reasonable N is small enough that integers suffice (up to N = 11).
215 m = 1
216 DO i = 2, (2*n-1)
217 tm = m
218 ti = i
219 r = mod(tm, ti)
220 DO WHILE (r .NE. 0)
221 tm = ti
222 ti = r
223 r = mod(tm, ti)
224 END DO
225 m = (m / ti) * i
226 END DO
227*
228* Generate the scaled Hilbert matrix in A
229* If we are testing SY routines,
230* take D1_i = D2_i, else, D1_i = D2_i*
231 IF ( lsamen( 2, c2, 'SY' ) ) THEN
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 $ * d1(mod(i,size_d)+1)
236 END DO
237 END DO
238 ELSE
239 DO j = 1, n
240 DO i = 1, n
241 a(i, j) = d1(mod(j,size_d)+1) * (dble(m) / (i + j - 1))
242 $ * d2(mod(i,size_d)+1)
243 END DO
244 END DO
245 END IF
246*
247* Generate matrix B as simply the first NRHS columns of M * the
248* identity.
249 tmp = dble(m)
250 CALL zlaset('Full', n, nrhs, (0.0d+0,0.0d+0), tmp, b, ldb)
251*
252* Generate the true solutions in X. Because B = the first NRHS
253* columns of M*I, the true solutions are just the first NRHS columns
254* of the inverse Hilbert matrix.
255 work(1) = n
256 DO j = 2, n
257 work(j) = ( ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1) )
258 $ * (n +j -1)
259 END DO
260
261* If we are testing SY routines,
262* take D1_i = D2_i, else, D1_i = D2_i*
263 IF ( lsamen( 2, c2, 'SY' ) ) THEN
264 DO j = 1, nrhs
265 DO i = 1, n
266 x(i, j) = invd1(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 ELSE
272 DO j = 1, nrhs
273 DO i = 1, n
274 x(i, j) = invd2(mod(j,size_d)+1) *
275 $ ((work(i)*work(j)) / (i + j - 1))
276 $ * invd1(mod(i,size_d)+1)
277 END DO
278 END DO
279 END IF
280 END
subroutine zlahilb(n, nrhs, a, lda, x, ldx, b, ldb, work, info, path)
ZLAHILB
Definition zlahilb.f:134
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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
logical function lsamen(n, ca, cb)
LSAMEN
Definition lsamen.f:74