LAPACK 3.3.0

slahilb.f

Go to the documentation of this file.
00001       SUBROUTINE SLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
00002 !
00003 !  -- LAPACK auxiliary test routine (version 3.0) --
00004 !     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
00005 !     Courant Institute, Argonne National Lab, and Rice University
00006 !     28 August, 2006
00007 !
00008 !     David Vu <dtv@cs.berkeley.edu>      
00009 !     Yozo Hida <yozo@cs.berkeley.edu>      
00010 !     Jason Riedy <ejr@cs.berkeley.edu>
00011 !     D. Halligan <dhalligan@berkeley.edu>
00012 !
00013       IMPLICIT NONE
00014 !     .. Scalar Arguments ..
00015       INTEGER N, NRHS, LDA, LDX, LDB, INFO
00016 !     .. Array Arguments ..
00017       REAL A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
00018 !     ..
00019 !
00020 !  Purpose
00021 !  =======
00022 !
00023 !  SLAHILB generates an N by N scaled Hilbert matrix in A along with
00024 !  NRHS right-hand sides in B and solutions in X such that A*X=B.
00025 !
00026 !  The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
00027 !  entries are integers.  The right-hand sides are the first NRHS 
00028 !  columns of M * the identity matrix, and the solutions are the 
00029 !  first NRHS columns of the inverse Hilbert matrix.
00030 !
00031 !  The condition number of the Hilbert matrix grows exponentially with
00032 !  its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
00033 !  Hilbert matrices beyond a relatively small dimension cannot be
00034 !  generated exactly without extra precision.  Precision is exhausted
00035 !  when the largest entry in the inverse Hilbert matrix is greater than
00036 !  2 to the power of the number of bits in the fraction of the data type
00037 !  used plus one, which is 24 for single precision.  
00038 !
00039 !  In single, the generated solution is exact for N <= 6 and has
00040 !  small componentwise error for 7 <= N <= 11.
00041 !
00042 !  Arguments
00043 !  =========
00044 !
00045 !  N       (input) INTEGER
00046 !          The dimension of the matrix A.
00047 !      
00048 !  NRHS    (input) NRHS
00049 !          The requested number of right-hand sides.
00050 !
00051 !  A       (output) REAL array, dimension (LDA, N)
00052 !          The generated scaled Hilbert matrix.
00053 !
00054 !  LDA     (input) INTEGER
00055 !          The leading dimension of the array A.  LDA >= N.
00056 !
00057 !  X       (output) REAL array, dimension (LDX, NRHS)
00058 !          The generated exact solutions.  Currently, the first NRHS
00059 !          columns of the inverse Hilbert matrix.
00060 !
00061 !  LDX     (input) INTEGER
00062 !          The leading dimension of the array X.  LDX >= N.
00063 !
00064 !  B       (output) REAL array, dimension (LDB, NRHS)
00065 !          The generated right-hand sides.  Currently, the first NRHS
00066 !          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
00067 !
00068 !  LDB     (input) INTEGER
00069 !          The leading dimension of the array B.  LDB >= N.
00070 !
00071 !  WORK    (workspace) REAL array, dimension (N)
00072 !
00073 !
00074 !  INFO    (output) INTEGER
00075 !          = 0: successful exit
00076 !          = 1: N is too large; the data is still generated but may not
00077 !               be not exact.
00078 !          < 0: if INFO = -i, the i-th argument had an illegal value
00079 !
00080 !  =====================================================================
00081 
00082 !     .. Local Scalars ..
00083       INTEGER TM, TI, R
00084       INTEGER M
00085       INTEGER I, J
00086 
00087 !     .. Parameters ..
00088 !     NMAX_EXACT   the largest dimension where the generated data is
00089 !                  exact.
00090 !     NMAX_APPROX  the largest dimension where the generated data has
00091 !                  a small componentwise relative error.
00092       INTEGER NMAX_EXACT, NMAX_APPROX
00093       PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11)
00094 
00095 !     ..
00096 !     .. External Functions
00097       EXTERNAL SLASET
00098       INTRINSIC REAL
00099 !     ..
00100 !     .. Executable Statements ..
00101 !
00102 !     Test the input arguments
00103 !
00104       INFO = 0
00105       IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
00106          INFO = -1
00107       ELSE IF (NRHS .LT. 0) THEN
00108          INFO = -2
00109       ELSE IF (LDA .LT. N) THEN
00110          INFO = -4
00111       ELSE IF (LDX .LT. N) THEN
00112          INFO = -6
00113       ELSE IF (LDB .LT. N) THEN
00114          INFO = -8
00115       END IF
00116       IF (INFO .LT. 0) THEN
00117          CALL XERBLA('SLAHILB', -INFO)
00118          RETURN
00119       END IF
00120       IF (N .GT. NMAX_EXACT) THEN
00121          INFO = 1
00122       END IF
00123 
00124 !     Compute M = the LCM of the integers [1, 2*N-1].  The largest
00125 !     reasonable N is small enough that integers suffice (up to N = 11).
00126       M = 1
00127       DO I = 2, (2*N-1)
00128          TM = M
00129          TI = I
00130          R = MOD(TM, TI)
00131          DO WHILE (R .NE. 0)
00132             TM = TI
00133             TI = R
00134             R = MOD(TM, TI)
00135          END DO
00136          M = (M / TI) * I
00137       END DO
00138 
00139 !     Generate the scaled Hilbert matrix in A
00140       DO J = 1, N
00141          DO I = 1, N
00142             A(I, J) = REAL(M) / (I + J - 1)
00143          END DO
00144       END DO
00145 
00146 !     Generate matrix B as simply the first NRHS columns of M * the
00147 !     identity.
00148       CALL SLASET('Full', N, NRHS, 0.0, REAL(M), B, LDB)
00149 
00150 !     Generate the true solutions in X.  Because B = the first NRHS
00151 !     columns of M*I, the true solutions are just the first NRHS columns
00152 !     of the inverse Hilbert matrix.
00153       WORK(1) = N
00154       DO J = 2, N
00155          WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
00156      $        * (N +J -1)
00157       END DO
00158       
00159       DO J = 1, NRHS
00160          DO I = 1, N
00161             X(I, J) = (WORK(I)*WORK(J)) / (I + J - 1)
00162          END DO
00163       END DO
00164 
00165       END
00166 
 All Files Functions