LAPACK 3.3.1
Linear Algebra PACKage

LIN/zlahilb.f

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