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 ! Yozo Hida ! Jason Riedy ! D. Halligan ! 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