LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slahilb.f
Go to the documentation of this file.
00001 *> \brief \b SLAHILB
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 SLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       INTEGER N, NRHS, LDA, LDX, LDB, INFO
00015 *       .. Array Arguments ..
00016 *       REAL A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
00017 *       ..
00018 *  
00019 *
00020 *> \par Purpose:
00021 *  =============
00022 *>
00023 *> \verbatim
00024 *>
00025 *> SLAHILB generates an N by N scaled Hilbert matrix in A along with
00026 *> NRHS right-hand sides in B and solutions in X such that A*X=B.
00027 *>
00028 *> The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
00029 *> entries are integers.  The right-hand sides are the first NRHS 
00030 *> columns of M * the identity matrix, and the solutions are the 
00031 *> first NRHS columns of the inverse Hilbert matrix.
00032 *>
00033 *> The condition number of the Hilbert matrix grows exponentially with
00034 *> its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
00035 *> Hilbert matrices beyond a relatively small dimension cannot be
00036 *> generated exactly without extra precision.  Precision is exhausted
00037 *> when the largest entry in the inverse Hilbert matrix is greater than
00038 *> 2 to the power of the number of bits in the fraction of the data type
00039 *> used plus one, which is 24 for single precision.  
00040 *>
00041 *> In single, the generated solution is exact for N <= 6 and has
00042 *> small componentwise error for 7 <= N <= 11.
00043 *> \endverbatim
00044 *
00045 *  Arguments:
00046 *  ==========
00047 *
00048 *> \param[in] N
00049 *> \verbatim
00050 *>          N is INTEGER
00051 *>          The dimension of the matrix A.
00052 *> \endverbatim
00053 *>      
00054 *> \param[in] NRHS
00055 *> \verbatim
00056 *>          NRHS is INTEGER
00057 *>          The requested number of right-hand sides.
00058 *> \endverbatim
00059 *>
00060 *> \param[out] A
00061 *> \verbatim
00062 *>          A is REAL array, dimension (LDA, N)
00063 *>          The generated scaled Hilbert matrix.
00064 *> \endverbatim
00065 *>
00066 *> \param[in] LDA
00067 *> \verbatim
00068 *>          LDA is INTEGER
00069 *>          The leading dimension of the array A.  LDA >= N.
00070 *> \endverbatim
00071 *>
00072 *> \param[out] X
00073 *> \verbatim
00074 *>          X is REAL array, dimension (LDX, NRHS)
00075 *>          The generated exact solutions.  Currently, the first NRHS
00076 *>          columns of the inverse Hilbert matrix.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] LDX
00080 *> \verbatim
00081 *>          LDX is INTEGER
00082 *>          The leading dimension of the array X.  LDX >= N.
00083 *> \endverbatim
00084 *>
00085 *> \param[out] B
00086 *> \verbatim
00087 *>          B is REAL array, dimension (LDB, NRHS)
00088 *>          The generated right-hand sides.  Currently, the first NRHS
00089 *>          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] LDB
00093 *> \verbatim
00094 *>          LDB is INTEGER
00095 *>          The leading dimension of the array B.  LDB >= N.
00096 *> \endverbatim
00097 *>
00098 *> \param[out] WORK
00099 *> \verbatim
00100 *>          WORK is REAL array, dimension (N)
00101 *> \endverbatim
00102 *>
00103 *> \param[out] INFO
00104 *> \verbatim
00105 *>          INFO is INTEGER
00106 *>          = 0: successful exit
00107 *>          = 1: N is too large; the data is still generated but may not
00108 *>               be not exact.
00109 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00110 *> \endverbatim
00111 *
00112 *  Authors:
00113 *  ========
00114 *
00115 *> \author Univ. of Tennessee 
00116 *> \author Univ. of California Berkeley 
00117 *> \author Univ. of Colorado Denver 
00118 *> \author NAG Ltd. 
00119 *
00120 *> \date November 2011
00121 *
00122 *> \ingroup real_matgen
00123 *
00124 *  =====================================================================
00125       SUBROUTINE SLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
00126 *
00127 *  -- LAPACK auxiliary routine (version 3.4.0) --
00128 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00129 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00130 *     November 2011
00131 *
00132 *     .. Scalar Arguments ..
00133       INTEGER N, NRHS, LDA, LDX, LDB, INFO
00134 *     .. Array Arguments ..
00135       REAL A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
00136 *     ..
00137 *
00138 *  =====================================================================
00139 *     .. Local Scalars ..
00140       INTEGER TM, TI, R
00141       INTEGER M
00142       INTEGER I, J
00143 
00144 *     .. Parameters ..
00145 *     NMAX_EXACT   the largest dimension where the generated data is
00146 *                  exact.
00147 *     NMAX_APPROX  the largest dimension where the generated data has
00148 *                  a small componentwise relative error.
00149       INTEGER NMAX_EXACT, NMAX_APPROX
00150       PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11)
00151 
00152 *     ..
00153 *     .. External Functions
00154       EXTERNAL SLASET
00155       INTRINSIC REAL
00156 *     ..
00157 *     .. Executable Statements ..
00158 *
00159 *     Test the input arguments
00160 *
00161       INFO = 0
00162       IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
00163          INFO = -1
00164       ELSE IF (NRHS .LT. 0) THEN
00165          INFO = -2
00166       ELSE IF (LDA .LT. N) THEN
00167          INFO = -4
00168       ELSE IF (LDX .LT. N) THEN
00169          INFO = -6
00170       ELSE IF (LDB .LT. N) THEN
00171          INFO = -8
00172       END IF
00173       IF (INFO .LT. 0) THEN
00174          CALL XERBLA('SLAHILB', -INFO)
00175          RETURN
00176       END IF
00177       IF (N .GT. NMAX_EXACT) THEN
00178          INFO = 1
00179       END IF
00180 
00181 *     Compute M = the LCM of the integers [1, 2*N-1].  The largest
00182 *     reasonable N is small enough that integers suffice (up to N = 11).
00183       M = 1
00184       DO I = 2, (2*N-1)
00185          TM = M
00186          TI = I
00187          R = MOD(TM, TI)
00188          DO WHILE (R .NE. 0)
00189             TM = TI
00190             TI = R
00191             R = MOD(TM, TI)
00192          END DO
00193          M = (M / TI) * I
00194       END DO
00195 
00196 *     Generate the scaled Hilbert matrix in A
00197       DO J = 1, N
00198          DO I = 1, N
00199             A(I, J) = REAL(M) / (I + J - 1)
00200          END DO
00201       END DO
00202 
00203 *     Generate matrix B as simply the first NRHS columns of M * the
00204 *     identity.
00205       CALL SLASET('Full', N, NRHS, 0.0, REAL(M), B, LDB)
00206 
00207 *     Generate the true solutions in X.  Because B = the first NRHS
00208 *     columns of M*I, the true solutions are just the first NRHS columns
00209 *     of the inverse Hilbert matrix.
00210       WORK(1) = N
00211       DO J = 2, N
00212          WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
00213      $        * (N +J -1)
00214       END DO
00215       
00216       DO J = 1, NRHS
00217          DO I = 1, N
00218             X(I, J) = (WORK(I)*WORK(J)) / (I + J - 1)
00219          END DO
00220       END DO
00221 
00222       END
00223 
 All Files Functions