SUBROUTINE CLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
$ INFO, PATH)
!
! -- LAPACK auxiliary test routine (version 3.0) --
! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
! Courant Institute, Argonne National Lab, and Rice University
! 28 August, 2006
!
! David Vu <dtv@cs.berkeley.edu>
! Yozo Hida <yozo@cs.berkeley.edu>
! Jason Riedy <ejr@cs.berkeley.edu>
! D. Halligan <dhalligan@berkeley.edu>
!
IMPLICIT NONE
! .. 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
! ..
!
! 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.
!
! Arguments
! =========
!
! N (input) INTEGER
! The dimension of the matrix A.
!
! NRHS (input) NRHS
! The requested number of right-hand sides.
!
! A (output) COMPLEX array, dimension (LDA, N)
! The generated scaled Hilbert matrix.
!
! LDA (input) INTEGER
! The leading dimension of the array A. LDA >= N.
!
! X (output) COMPLEX array, dimension (LDX, NRHS)
! The generated exact solutions. Currently, the first NRHS
! columns of the inverse Hilbert matrix.
!
! LDX (input) INTEGER
! The leading dimension of the array X. LDX >= N.
!
! B (output) 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.
!
! LDB (input) INTEGER
! The leading dimension of the array B. LDB >= N.
!
! WORK (workspace) REAL array, dimension (N)
!
!
! INFO (output) 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
!
! =====================================================================
! .. 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 permuation 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
END