![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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 NRHS 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 single_lin 00123 * 00124 * ===================================================================== 00125 SUBROUTINE SLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO) 00126 * 00127 * -- LAPACK test 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 * .. External Functions 00153 EXTERNAL SLASET 00154 INTRINSIC REAL 00155 * .. 00156 * .. Executable Statements .. 00157 * 00158 * Test the input arguments 00159 * 00160 INFO = 0 00161 IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN 00162 INFO = -1 00163 ELSE IF (NRHS .LT. 0) THEN 00164 INFO = -2 00165 ELSE IF (LDA .LT. N) THEN 00166 INFO = -4 00167 ELSE IF (LDX .LT. N) THEN 00168 INFO = -6 00169 ELSE IF (LDB .LT. N) THEN 00170 INFO = -8 00171 END IF 00172 IF (INFO .LT. 0) THEN 00173 CALL XERBLA('SLAHILB', -INFO) 00174 RETURN 00175 END IF 00176 IF (N .GT. NMAX_EXACT) THEN 00177 INFO = 1 00178 END IF 00179 * 00180 * Compute M = the LCM of the integers [1, 2*N-1]. The largest 00181 * reasonable N is small enough that integers suffice (up to N = 11). 00182 M = 1 00183 DO I = 2, (2*N-1) 00184 TM = M 00185 TI = I 00186 R = MOD(TM, TI) 00187 DO WHILE (R .NE. 0) 00188 TM = TI 00189 TI = R 00190 R = MOD(TM, TI) 00191 END DO 00192 M = (M / TI) * I 00193 END DO 00194 * 00195 * Generate the scaled Hilbert matrix in A 00196 DO J = 1, N 00197 DO I = 1, N 00198 A(I, J) = REAL(M) / (I + J - 1) 00199 END DO 00200 END DO 00201 * 00202 * Generate matrix B as simply the first NRHS columns of M * the 00203 * identity. 00204 CALL SLASET('Full', N, NRHS, 0.0, REAL(M), B, LDB) 00205 * 00206 * Generate the true solutions in X. Because B = the first NRHS 00207 * columns of M*I, the true solutions are just the first NRHS columns 00208 * of the inverse Hilbert matrix. 00209 WORK(1) = N 00210 DO J = 2, N 00211 WORK(J) = ( ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1) ) 00212 $ * (N +J -1) 00213 END DO 00214 * 00215 DO J = 1, NRHS 00216 DO I = 1, N 00217 X(I, J) = (WORK(I)*WORK(J)) / (I + J - 1) 00218 END DO 00219 END DO 00220 * 00221 END 00222