LAPACK 3.3.1 Linear Algebra PACKage

# MATGEN/clahilb.f

Go to the documentation of this file.
```00001       SUBROUTINE CLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
00002      \$     INFO, PATH)
00003 !
00004 !  -- LAPACK auxiliary test routine (version 3.2.2) --
00005 !     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
00006 !     Courant Institute, Argonne National Lab, and Rice University
00007 *     June 2010
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       REAL WORK(N)
00019       COMPLEX A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
00020       CHARACTER*3        PATH
00021 !     ..
00022 !
00023 !  Purpose
00024 !  =======
00025 !
00026 !  CLAHILB 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) INTEGER
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 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 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 !     ..
00112 !     .. External Functions
00113       EXTERNAL CLASET, LSAMEN
00114       INTRINSIC REAL
00115       LOGICAL LSAMEN
00116 !     ..
00117 !     .. Executable Statements ..
00118       C2 = PATH( 2: 3 )
00119 !
00120 !     Test the input arguments
00121 !
00122       INFO = 0
00123       IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
00124          INFO = -1
00125       ELSE IF (NRHS .LT. 0) THEN
00126          INFO = -2
00127       ELSE IF (LDA .LT. N) THEN
00128          INFO = -4
00129       ELSE IF (LDX .LT. N) THEN
00130          INFO = -6
00131       ELSE IF (LDB .LT. N) THEN
00132          INFO = -8
00133       END IF
00134       IF (INFO .LT. 0) THEN
00135          CALL XERBLA('CLAHILB', -INFO)
00136          RETURN
00137       END IF
00138       IF (N .GT. NMAX_EXACT) THEN
00139          INFO = 1
00140       END IF
00141
00142 !     Compute M = the LCM of the integers [1, 2*N-1].  The largest
00143 !     reasonable N is small enough that integers suffice (up to N = 11).
00144       M = 1
00145       DO I = 2, (2*N-1)
00146          TM = M
00147          TI = I
00148          R = MOD(TM, TI)
00149          DO WHILE (R .NE. 0)
00150             TM = TI
00151             TI = R
00152             R = MOD(TM, TI)
00153          END DO
00154          M = (M / TI) * I
00155       END DO
00156
00157 !     Generate the scaled Hilbert matrix in A
00158 !     If we are testing SY routines, take
00159 !          D1_i = D2_i, else, D1_i = D2_i*
00160       IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
00161          DO J = 1, N
00162             DO I = 1, N
00163                A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1))
00164      \$              * D1(MOD(I,SIZE_D)+1)
00165             END DO
00166          END DO
00167       ELSE
00168          DO J = 1, N
00169             DO I = 1, N
00170                A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1))
00171      \$              * D2(MOD(I,SIZE_D)+1)
00172             END DO
00173          END DO
00174       END IF
00175
00176 !     Generate matrix B as simply the first NRHS columns of M * the
00177 !     identity.
00178       TMP = REAL(M)
00179       CALL CLASET('Full', N, NRHS, (0.0,0.0), TMP, B, LDB)
00180
00181 !     Generate the true solutions in X.  Because B = the first NRHS
00182 !     columns of M*I, the true solutions are just the first NRHS columns
00183 !     of the inverse Hilbert matrix.
00184       WORK(1) = N
00185       DO J = 2, N
00186          WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
00187      \$        * (N +J -1)
00188       END DO
00189
00190 !     If we are testing SY routines,
00191 !            take D1_i = D2_i, else, D1_i = D2_i*
00192       IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
00193          DO J = 1, NRHS
00194             DO I = 1, N
00195                X(I, J) =
00196      \$              INVD1(MOD(J,SIZE_D)+1) *
00197      \$              ((WORK(I)*WORK(J)) / (I + J - 1))
00198      \$              * INVD1(MOD(I,SIZE_D)+1)
00199             END DO
00200          END DO
00201       ELSE
00202          DO J = 1, NRHS
00203             DO I = 1, N
00204                X(I, J) =
00205      \$              INVD2(MOD(J,SIZE_D)+1) *
00206      \$              ((WORK(I)*WORK(J)) / (I + J - 1))
00207      \$              * INVD1(MOD(I,SIZE_D)+1)
00208             END DO
00209          END DO
00210       END IF
00211       END
00212
```