LAPACK 3.3.0

dlahilb.f

Go to the documentation of this file.
00001       SUBROUTINE DLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
00002 !
00003 !  -- LAPACK auxiliary test routine (version 3.2.2) --
00004 !     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
00005 !     Courant Institute, Argonne National Lab, and Rice University
00006 *     June 2010
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       DOUBLE PRECISION A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
00018 !     ..
00019 !
00020 !  Purpose
00021 !  =======
00022 !
00023 !  DLAHILB 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) INTEGER
00049 !          The requested number of right-hand sides.
00050 !
00051 !  A       (output) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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       COMPLEX*16 TMP
00087 
00088 !     .. Parameters ..
00089 !     NMAX_EXACT   the largest dimension where the generated data is
00090 !                  exact.
00091 !     NMAX_APPROX  the largest dimension where the generated data has
00092 !                  a small componentwise relative error.
00093       INTEGER NMAX_EXACT, NMAX_APPROX
00094       PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11)
00095 
00096 !     ..
00097 !     .. External Functions
00098       EXTERNAL DLASET
00099       INTRINSIC DBLE
00100 !     ..
00101 !     .. Executable Statements ..
00102 !
00103 !     Test the input arguments
00104 !
00105       INFO = 0
00106       IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
00107          INFO = -1
00108       ELSE IF (NRHS .LT. 0) THEN
00109          INFO = -2
00110       ELSE IF (LDA .LT. N) THEN
00111          INFO = -4
00112       ELSE IF (LDX .LT. N) THEN
00113          INFO = -6
00114       ELSE IF (LDB .LT. N) THEN
00115          INFO = -8
00116       END IF
00117       IF (INFO .LT. 0) THEN
00118          CALL XERBLA('DLAHILB', -INFO)
00119          RETURN
00120       END IF
00121       IF (N .GT. NMAX_EXACT) THEN
00122          INFO = 1
00123       END IF
00124 
00125 !     Compute M = the LCM of the integers [1, 2*N-1].  The largest
00126 !     reasonable N is small enough that integers suffice (up to N = 11).
00127       M = 1
00128       DO I = 2, (2*N-1)
00129          TM = M
00130          TI = I
00131          R = MOD(TM, TI)
00132          DO WHILE (R .NE. 0)
00133             TM = TI
00134             TI = R
00135             R = MOD(TM, TI)
00136          END DO
00137          M = (M / TI) * I
00138       END DO
00139 
00140 !     Generate the scaled Hilbert matrix in A
00141       DO J = 1, N
00142          DO I = 1, N
00143             A(I, J) = DBLE(M) / (I + J - 1)
00144          END DO
00145       END DO
00146 
00147 !     Generate matrix B as simply the first NRHS columns of M * the
00148 !     identity.
00149       TMP = DBLE(M)
00150       CALL DLASET('Full', N, NRHS, 0.0D+0, TMP, B, LDB)
00151 
00152 !     Generate the true solutions in X.  Because B = the first NRHS
00153 !     columns of M*I, the true solutions are just the first NRHS columns
00154 !     of the inverse Hilbert matrix.
00155       WORK(1) = N
00156       DO J = 2, N
00157          WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
00158      $        * (N +J -1)
00159       END DO
00160       
00161       DO J = 1, NRHS
00162          DO I = 1, N
00163             X(I, J) = (WORK(I)*WORK(J)) / (I + J - 1)
00164          END DO
00165       END DO
00166 
00167       END
00168 
 All Files Functions