subroutine slahilb	(	integer	N,
		integer	NRHS,
		real, dimension(lda, n)	A,
		integer	LDA,
		real, dimension(ldx, nrhs)	X,
		integer	LDX,
		real, dimension(ldb, nrhs)	B,
		integer	LDB,
		real, dimension(n)	WORK,
		integer	INFO
	)

SLAHILB

Purpose:

 SLAHILB 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 REAL 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 REAL 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

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Date

November 2015

Definition at line 126 of file slahilb.f.

*
*  -- LAPACK test routine (version 3.6.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2015
*
*     .. Scalar Arguments ..
      INTEGER n, nrhs, lda, ldx, ldb, info
*     .. Array Arguments ..
      REAL a(lda, n), x(ldx, nrhs), b(ldb, nrhs), work(n)
*     ..
*
*  =====================================================================
*     .. Local Scalars ..
      INTEGER tm, ti, r
      INTEGER m
      INTEGER i, j

*     .. 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.
      INTEGER nmax_exact, nmax_approx
      parameter(nmax_exact = 6, nmax_approx = 11)

*     ..
*     .. External Functions
      EXTERNAL slaset
      INTRINSIC real
*     ..
*     .. Executable Statements ..
*
*     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('SLAHILB', -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
      DO j = 1, n
         DO i = 1, n
            a(i, j) = REAL(M) / (i + j - 1)
         END DO
      END DO

*     Generate matrix B as simply the first NRHS columns of M * the
*     identity.
      CALL slaset('Full', n, nrhs, 0.0, REAL(M), 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
      
      DO j = 1, nrhs
         DO i = 1, n
            x(i, j) = (work(i)*work(j)) / (i + j - 1)
         END DO
      END DO

Here is the call graph for this function: