![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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 NRHS 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_lin 00132 * 00133 * ===================================================================== 00134 SUBROUTINE ZLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, 00135 $ INFO, PATH) 00136 * 00137 * -- LAPACK test 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, take D1_i = D2_i, else, D1_i = D2_i* 00224 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 00225 DO J = 1, N 00226 DO I = 1, N 00227 A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1)) 00228 $ * D1(MOD(I,SIZE_D)+1) 00229 END DO 00230 END DO 00231 ELSE 00232 DO J = 1, N 00233 DO I = 1, N 00234 A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1)) 00235 $ * D2(MOD(I,SIZE_D)+1) 00236 END DO 00237 END DO 00238 END IF 00239 * 00240 * Generate matrix B as simply the first NRHS columns of M * the 00241 * identity. 00242 TMP = DBLE(M) 00243 CALL ZLASET('Full', N, NRHS, (0.0D+0,0.0D+0), TMP, B, LDB) 00244 * 00245 * Generate the true solutions in X. Because B = the first NRHS 00246 * columns of M*I, the true solutions are just the first NRHS columns 00247 * of the inverse Hilbert matrix. 00248 WORK(1) = N 00249 DO J = 2, N 00250 WORK(J) = ( ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1) ) 00251 $ * (N +J -1) 00252 END DO 00253 * 00254 * If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i* 00255 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 00256 DO J = 1, NRHS 00257 DO I = 1, N 00258 X(I, J) = INVD1(MOD(J,SIZE_D)+1) * 00259 $ ((WORK(I)*WORK(J)) / (I + J - 1)) 00260 $ * INVD1(MOD(I,SIZE_D)+1) 00261 END DO 00262 END DO 00263 ELSE 00264 DO J = 1, NRHS 00265 DO I = 1, N 00266 X(I, J) = INVD2(MOD(J,SIZE_D)+1) * 00267 $ ((WORK(I)*WORK(J)) / (I + J - 1)) 00268 $ * INVD1(MOD(I,SIZE_D)+1) 00269 END DO 00270 END DO 00271 END IF 00272 END