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: