LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zlahilb.f
Go to the documentation of this file.
00001 *> \brief \b ZLAHILB
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE ZLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
00012 *            INFO, PATH)
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER N, NRHS, LDA, LDX, LDB, INFO
00016 *       .. Array Arguments ..
00017 *       DOUBLE PRECISION WORK(N)
00018 *       COMPLEX*16 A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
00019 *       CHARACTER*3 PATH
00020 *       ..
00021 *  
00022 *
00023 *> \par Purpose:
00024 *  =============
00025 *>
00026 *> \verbatim
00027 *>
00028 *> ZLAHILB generates an N by N scaled Hilbert matrix in A along with
00029 *> NRHS right-hand sides in B and solutions in X such that A*X=B.
00030 *>
00031 *> The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
00032 *> entries are integers.  The right-hand sides are the first NRHS
00033 *> columns of M * the identity matrix, and the solutions are the
00034 *> first NRHS columns of the inverse Hilbert matrix.
00035 *>
00036 *> The condition number of the Hilbert matrix grows exponentially with
00037 *> its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
00038 *> Hilbert matrices beyond a relatively small dimension cannot be
00039 *> generated exactly without extra precision.  Precision is exhausted
00040 *> when the largest entry in the inverse Hilbert matrix is greater than
00041 *> 2 to the power of the number of bits in the fraction of the data type
00042 *> used plus one, which is 24 for single precision.
00043 *>
00044 *> In single, the generated solution is exact for N <= 6 and has
00045 *> small componentwise error for 7 <= N <= 11.
00046 *> \endverbatim
00047 *
00048 *  Arguments:
00049 *  ==========
00050 *
00051 *> \param[in] N
00052 *> \verbatim
00053 *>          N is INTEGER
00054 *>          The dimension of the matrix A.
00055 *> \endverbatim
00056 *>
00057 *> \param[in] NRHS
00058 *> \verbatim
00059 *>          NRHS is INTEGER
00060 *>          The requested number of right-hand sides.
00061 *> \endverbatim
00062 *>
00063 *> \param[out] A
00064 *> \verbatim
00065 *>          A is COMPLEX array, dimension (LDA, N)
00066 *>          The generated scaled Hilbert matrix.
00067 *> \endverbatim
00068 *>
00069 *> \param[in] LDA
00070 *> \verbatim
00071 *>          LDA is INTEGER
00072 *>          The leading dimension of the array A.  LDA >= N.
00073 *> \endverbatim
00074 *>
00075 *> \param[out] X
00076 *> \verbatim
00077 *>          X is COMPLEX array, dimension (LDX, NRHS)
00078 *>          The generated exact solutions.  Currently, the first NRHS
00079 *>          columns of the inverse Hilbert matrix.
00080 *> \endverbatim
00081 *>
00082 *> \param[in] LDX
00083 *> \verbatim
00084 *>          LDX is INTEGER
00085 *>          The leading dimension of the array X.  LDX >= N.
00086 *> \endverbatim
00087 *>
00088 *> \param[out] B
00089 *> \verbatim
00090 *>          B is REAL array, dimension (LDB, NRHS)
00091 *>          The generated right-hand sides.  Currently, the first NRHS
00092 *>          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] LDB
00096 *> \verbatim
00097 *>          LDB is INTEGER
00098 *>          The leading dimension of the array B.  LDB >= N.
00099 *> \endverbatim
00100 *>
00101 *> \param[out] WORK
00102 *> \verbatim
00103 *>          WORK is REAL array, dimension (N)
00104 *> \endverbatim
00105 *>
00106 *> \param[out] INFO
00107 *> \verbatim
00108 *>          INFO is INTEGER
00109 *>          = 0: successful exit
00110 *>          = 1: N is too large; the data is still generated but may not
00111 *>               be not exact.
00112 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00113 *> \endverbatim
00114 *>
00115 *> \param[in] PATH
00116 *> \verbatim
00117 *>          PATH is CHARACTER*3
00118 *>          The LAPACK path name.
00119 *> \endverbatim
00120 *
00121 *  Authors:
00122 *  ========
00123 *
00124 *> \author Univ. of Tennessee 
00125 *> \author Univ. of California Berkeley 
00126 *> \author Univ. of Colorado Denver 
00127 *> \author NAG Ltd. 
00128 *
00129 *> \date November 2011
00130 *
00131 *> \ingroup complex16_matgen
00132 *
00133 *  =====================================================================
00134       SUBROUTINE ZLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
00135      $     INFO, PATH)
00136 *
00137 *  -- LAPACK auxiliary routine (version 3.4.0) --
00138 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00139 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00140 *     November 2011
00141 *
00142 *     .. Scalar Arguments ..
00143       INTEGER N, NRHS, LDA, LDX, LDB, INFO
00144 *     .. Array Arguments ..
00145       DOUBLE PRECISION WORK(N)
00146       COMPLEX*16 A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
00147       CHARACTER*3 PATH
00148 *     ..
00149 *
00150 *  =====================================================================
00151 *     .. Local Scalars ..
00152       INTEGER TM, TI, R
00153       INTEGER M
00154       INTEGER I, J
00155       COMPLEX*16 TMP
00156       CHARACTER*2 C2
00157 
00158 *     .. Parameters ..
00159 *     NMAX_EXACT   the largest dimension where the generated data is
00160 *                  exact.
00161 *     NMAX_APPROX  the largest dimension where the generated data has
00162 *                  a small componentwise relative error.
00163 *     ??? complex uses how many bits ???
00164       INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D
00165       PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11, SIZE_D = 8)
00166 
00167 *     d's are generated from random permuation of those eight elements.
00168       COMPLEX*16 d1(8), d2(8), invd1(8), invd2(8)
00169       DATA D1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/
00170       DATA D2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/
00171       
00172       DATA INVD1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0),
00173      $     (-.5,-.5),(.5,-.5),(.5,.5)/
00174       DATA INVD2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0),
00175      $     (-.5,.5),(.5,.5),(.5,-.5)/
00176 *     ..
00177 *     .. External Functions
00178       EXTERNAL ZLASET, LSAMEN
00179       INTRINSIC DBLE
00180       LOGICAL LSAMEN
00181 *     ..
00182 *     .. Executable Statements ..
00183       C2 = PATH( 2: 3 )
00184 *
00185 *     Test the input arguments
00186 *
00187       INFO = 0
00188       IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
00189          INFO = -1
00190       ELSE IF (NRHS .LT. 0) THEN
00191          INFO = -2
00192       ELSE IF (LDA .LT. N) THEN
00193          INFO = -4
00194       ELSE IF (LDX .LT. N) THEN
00195          INFO = -6
00196       ELSE IF (LDB .LT. N) THEN
00197          INFO = -8
00198       END IF
00199       IF (INFO .LT. 0) THEN
00200          CALL XERBLA('ZLAHILB', -INFO)
00201          RETURN
00202       END IF
00203       IF (N .GT. NMAX_EXACT) THEN
00204          INFO = 1
00205       END IF
00206 
00207 *     Compute M = the LCM of the integers [1, 2*N-1].  The largest
00208 *     reasonable N is small enough that integers suffice (up to N = 11).
00209       M = 1
00210       DO I = 2, (2*N-1)
00211          TM = M
00212          TI = I
00213          R = MOD(TM, TI)
00214          DO WHILE (R .NE. 0)
00215             TM = TI
00216             TI = R
00217             R = MOD(TM, TI)
00218          END DO
00219          M = (M / TI) * I
00220       END DO
00221 
00222 *     Generate the scaled Hilbert matrix in A
00223 *     If we are testing SY routines,
00224 *        take D1_i = D2_i, else, D1_i = D2_i*
00225       IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
00226          DO J = 1, N
00227             DO I = 1, N
00228                A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1))
00229      $              * D1(MOD(I,SIZE_D)+1)
00230             END DO
00231          END DO
00232       ELSE
00233          DO J = 1, N
00234             DO I = 1, N
00235                A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1))
00236      $              * D2(MOD(I,SIZE_D)+1)
00237             END DO
00238          END DO
00239       END IF
00240       
00241 *     Generate matrix B as simply the first NRHS columns of M * the
00242 *     identity.
00243       TMP = DBLE(M)
00244       CALL ZLASET('Full', N, NRHS, (0.0D+0,0.0D+0), TMP, B, LDB)
00245 
00246 *     Generate the true solutions in X.  Because B = the first NRHS
00247 *     columns of M*I, the true solutions are just the first NRHS columns
00248 *     of the inverse Hilbert matrix.
00249       WORK(1) = N
00250       DO J = 2, N
00251          WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
00252      $        * (N +J -1)
00253       END DO
00254 
00255 *     If we are testing SY routines, 
00256 *           take D1_i = D2_i, else, D1_i = D2_i*
00257       IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
00258          DO J = 1, NRHS
00259             DO I = 1, N
00260                X(I, J) = INVD1(MOD(J,SIZE_D)+1) *
00261      $              ((WORK(I)*WORK(J)) / (I + J - 1))
00262      $              * INVD1(MOD(I,SIZE_D)+1)
00263             END DO
00264          END DO
00265       ELSE
00266          DO J = 1, NRHS
00267             DO I = 1, N
00268                X(I, J) = INVD2(MOD(J,SIZE_D)+1) *
00269      $              ((WORK(I)*WORK(J)) / (I + J - 1))
00270      $              * INVD1(MOD(I,SIZE_D)+1)
00271             END DO
00272          END DO
00273       END IF
00274       END
 All Files Functions