LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \ingroup complex_lin
130*
131* =====================================================================
132 SUBROUTINE clahilb( 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 REAL WORK(N)
143 COMPLEX 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 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 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 claset, lsamen
176 INTRINSIC real
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('CLAHILB', -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, take
221* 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) * (real(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) * (real(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 = real(m)
241 CALL claset('Full', n, nrhs, (0.0,0.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) =
258 $ 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) =
267 $ invd2(mod(j,size_d)+1) *
268 $ ((work(i)*work(j)) / (i + j - 1))
269 $ * invd1(mod(i,size_d)+1)
270 END DO
271 END DO
272 END IF
273 END
274
subroutine clahilb(n, nrhs, a, lda, x, ldx, b, ldb, work, info, path)
CLAHILB
Definition clahilb.f:134
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
logical function lsamen(n, ca, cb)
LSAMEN
Definition lsamen.f:74