LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
November 2017

Definition at line 136 of file zlahilb.f.

136 *
137 * -- LAPACK test routine (version 3.8.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 2017
141 *
142 * .. Scalar Arguments ..
143  INTEGER n, nrhs, lda, ldx, ldb, info
144 * .. Array Arguments ..
145  DOUBLE PRECISION work(n)
146  COMPLEX*16 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*16 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*16 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 * .. External Subroutines ..
178  EXTERNAL xerbla
179 * ..
180 * .. External Functions
181  EXTERNAL zlaset, lsamen
182  INTRINSIC dble
183  LOGICAL lsamen
184 * ..
185 * .. Executable Statements ..
186  c2 = path( 2: 3 )
187 *
188 * Test the input arguments
189 *
190  info = 0
191  IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
192  info = -1
193  ELSE IF (nrhs .LT. 0) THEN
194  info = -2
195  ELSE IF (lda .LT. n) THEN
196  info = -4
197  ELSE IF (ldx .LT. n) THEN
198  info = -6
199  ELSE IF (ldb .LT. n) THEN
200  info = -8
201  END IF
202  IF (info .LT. 0) THEN
203  CALL xerbla('ZLAHILB', -info)
204  RETURN
205  END IF
206  IF (n .GT. nmax_exact) THEN
207  info = 1
208  END IF
209 *
210 * Compute M = the LCM of the integers [1, 2*N-1]. The largest
211 * reasonable N is small enough that integers suffice (up to N = 11).
212  m = 1
213  DO i = 2, (2*n-1)
214  tm = m
215  ti = i
216  r = mod(tm, ti)
217  DO WHILE (r .NE. 0)
218  tm = ti
219  ti = r
220  r = mod(tm, ti)
221  END DO
222  m = (m / ti) * i
223  END DO
224 *
225 * Generate the scaled Hilbert matrix in A
226 * If we are testing SY routines,
227 * take D1_i = D2_i, else, D1_i = D2_i*
228  IF ( lsamen( 2, c2, 'SY' ) ) THEN
229  DO j = 1, n
230  DO i = 1, n
231  a(i, j) = d1(mod(j,size_d)+1) * (dble(m) / (i + j - 1))
232  $ * d1(mod(i,size_d)+1)
233  END DO
234  END DO
235  ELSE
236  DO j = 1, n
237  DO i = 1, n
238  a(i, j) = d1(mod(j,size_d)+1) * (dble(m) / (i + j - 1))
239  $ * d2(mod(i,size_d)+1)
240  END DO
241  END DO
242  END IF
243 *
244 * Generate matrix B as simply the first NRHS columns of M * the
245 * identity.
246  tmp = dble(m)
247  CALL zlaset('Full', n, nrhs, (0.0d+0,0.0d+0), tmp, b, ldb)
248 *
249 * Generate the true solutions in X. Because B = the first NRHS
250 * columns of M*I, the true solutions are just the first NRHS columns
251 * of the inverse Hilbert matrix.
252  work(1) = n
253  DO j = 2, n
254  work(j) = ( ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1) )
255  $ * (n +j -1)
256  END DO
257 
258 * If we are testing SY routines,
259 * take D1_i = D2_i, else, D1_i = D2_i*
260  IF ( lsamen( 2, c2, 'SY' ) ) THEN
261  DO j = 1, nrhs
262  DO i = 1, n
263  x(i, j) = invd1(mod(j,size_d)+1) *
264  $ ((work(i)*work(j)) / (i + j - 1))
265  $ * invd1(mod(i,size_d)+1)
266  END DO
267  END DO
268  ELSE
269  DO j = 1, nrhs
270  DO i = 1, n
271  x(i, j) = 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
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:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsamen(N, CA, CB)
LSAMEN
Definition: lsamen.f:76
Here is the call graph for this function: