LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zlahilb()

subroutine zlahilb ( integer  n,
integer  nrhs,
complex*16, dimension(lda,n)  a,
integer  lda,
complex*16, dimension(ldx, nrhs)  x,
integer  ldx,
complex*16, dimension(ldb, nrhs)  b,
integer  ldb,
double precision, dimension(n)  work,
integer  info,
character*3  path 
)

ZLAHILB

Purpose:
 ZLAHILB generates an N by N scaled Hilbert matrix in A along with
 NRHS right-hand sides in B and solutions in X such that A*X=B.

 The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
 entries are integers.  The right-hand sides are the first NRHS
 columns of M * the identity matrix, and the solutions are the
 first NRHS columns of the inverse Hilbert matrix.

 The condition number of the Hilbert matrix grows exponentially with
 its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
 Hilbert matrices beyond a relatively small dimension cannot be
 generated exactly without extra precision.  Precision is exhausted
 when the largest entry in the inverse Hilbert matrix is greater than
 2 to the power of the number of bits in the fraction of the data type
 used plus one, which is 24 for single precision.

 In single, the generated solution is exact for N <= 6 and has
 small componentwise error for 7 <= N <= 11.
Parameters
[in]N
          N is INTEGER
          The dimension of the matrix A.
[in]NRHS
          NRHS is INTEGER
          The requested number of right-hand sides.
[out]A
          A is COMPLEX array, dimension (LDA, N)
          The generated scaled Hilbert matrix.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= N.
[out]X
          X is COMPLEX array, dimension (LDX, NRHS)
          The generated exact solutions.  Currently, the first NRHS
          columns of the inverse Hilbert matrix.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= N.
[out]B
          B is REAL array, dimension (LDB, NRHS)
          The generated right-hand sides.  Currently, the first NRHS
          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= N.
[out]WORK
          WORK is REAL array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          = 1: N is too large; the data is still generated but may not
               be not exact.
          < 0: if INFO = -i, the i-th argument had an illegal value
[in]PATH
          PATH is CHARACTER*3
          The LAPACK path name.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 132 of file zlahilb.f.

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
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zlaset
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
lsamen
logical function lsamen(n, ca, cb)
LSAMEN
Definition lsamen.f:74
Here is the call graph for this function: