clahilb()

subroutine clahilb	(	integer	n,
		integer	nrhs,
		complex, dimension(lda,n)	a,
		integer	lda,
		complex, dimension(ldx, nrhs)	x,
		integer	ldx,
		complex, dimension(ldb, nrhs)	b,
		integer	ldb,
		real, dimension(n)	work,
		integer	info,
		character*3	path
	)

CLAHILB

Purpose:

 CLAHILB 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 clahilb.f.

*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER N, NRHS, LDA, LDX, LDB, INFO
*     .. Array Arguments ..
      REAL WORK(N)
      COMPLEX A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
      CHARACTER*3        PATH
*     ..
*
*  =====================================================================
*     .. Local Scalars ..
      INTEGER TM, TI, R
      INTEGER M
      INTEGER I, J
      COMPLEX TMP
      CHARACTER*2 C2
*     ..
*     .. Parameters ..
*     NMAX_EXACT   the largest dimension where the generated data is
*                  exact.
*     NMAX_APPROX  the largest dimension where the generated data has
*                  a small componentwise relative error.
*     ??? complex uses how many bits ???
      INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D
      parameter(nmax_exact = 6, nmax_approx = 11, size_d = 8)
*
*     d's are generated from random permutation of those eight elements.
      COMPLEX D1(8), D2(8), INVD1(8), INVD2(8)
      DATA d1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/
      DATA d2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/
 
      DATA invd1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0),
     $     (-.5,-.5),(.5,-.5),(.5,.5)/
      DATA invd2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0),
     $     (-.5,.5),(.5,.5),(.5,-.5)/
*     ..
*     .. External Functions
      EXTERNAL claset, lsamen
      INTRINSIC real
      LOGICAL LSAMEN
*     ..
*     .. Executable Statements ..
      c2 = path( 2: 3 )
*
*     Test the input arguments
*
      info = 0
      IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
         info = -1
      ELSE IF (nrhs .LT. 0) THEN
         info = -2
      ELSE IF (lda .LT. n) THEN
         info = -4
      ELSE IF (ldx .LT. n) THEN
         info = -6
      ELSE IF (ldb .LT. n) THEN
         info = -8
      END IF
      IF (info .LT. 0) THEN
         CALL xerbla('CLAHILB', -info)
         RETURN
      END IF
      IF (n .GT. nmax_exact) THEN
         info = 1
      END IF
*
*     Compute M = the LCM of the integers [1, 2*N-1].  The largest
*     reasonable N is small enough that integers suffice (up to N = 11).
      m = 1
      DO i = 2, (2*n-1)
         tm = m
         ti = i
         r = mod(tm, ti)
         DO WHILE (r .NE. 0)
            tm = ti
            ti = r
            r = mod(tm, ti)
         END DO
         m = (m / ti) * i
      END DO
*
*     Generate the scaled Hilbert matrix in A
*     If we are testing SY routines, take
*          D1_i = D2_i, else, D1_i = D2_i*
      IF ( lsamen( 2, c2, 'SY' ) ) THEN
         DO j = 1, n
            DO i = 1, n
               a(i, j) = d1(mod(j,size_d)+1) * (real(m) / (i + j - 1))
     $              * d1(mod(i,size_d)+1)
            END DO
         END DO
      ELSE
         DO j = 1, n
            DO i = 1, n
               a(i, j) = d1(mod(j,size_d)+1) * (real(m) / (i + j - 1))
     $              * d2(mod(i,size_d)+1)
            END DO
         END DO
      END IF
*
*     Generate matrix B as simply the first NRHS columns of M * the
*     identity.
      tmp = real(m)
      CALL claset('Full', n, nrhs, (0.0,0.0), tmp, b, ldb)
*
*     Generate the true solutions in X.  Because B = the first NRHS
*     columns of M*I, the true solutions are just the first NRHS columns
*     of the inverse Hilbert matrix.
      work(1) = n
      DO j = 2, n
         work(j) = (  ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1)  )
     $        * (n +j -1)
      END DO
 
*     If we are testing SY routines,
*            take D1_i = D2_i, else, D1_i = D2_i*
      IF ( lsamen( 2, c2, 'SY' ) ) THEN
         DO j = 1, nrhs
            DO i = 1, n
               x(i, j) =
     $              invd1(mod(j,size_d)+1) *
     $              ((work(i)*work(j)) / (i + j - 1))
     $              * invd1(mod(i,size_d)+1)
            END DO
         END DO
      ELSE
         DO j = 1, nrhs
            DO i = 1, n
               x(i, j) =
     $              invd2(mod(j,size_d)+1) *
     $              ((work(i)*work(j)) / (i + j - 1))
     $              * invd1(mod(i,size_d)+1)
            END DO
         END DO
      END IF

Here is the call graph for this function:

Here is the caller graph for this function: