1*> \brief \b ZLAHILB 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE ZLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, 12* INFO, PATH) 13* 14* .. Scalar Arguments .. 15* INTEGER N, NRHS, LDA, LDX, LDB, INFO 16* .. Array Arguments .. 17* DOUBLE PRECISION WORK(N) 18* COMPLEX*16 A(LDA,N), X(LDX, NRHS), B(LDB, NRHS) 19* CHARACTER*3 PATH 20* .. 21* 22* 23*> \par Purpose: 24* ============= 25*> 26*> \verbatim 27*> 28*> ZLAHILB generates an N by N scaled Hilbert matrix in A along with 29*> NRHS right-hand sides in B and solutions in X such that A*X=B. 30*> 31*> The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all 32*> entries are integers. The right-hand sides are the first NRHS 33*> columns of M * the identity matrix, and the solutions are the 34*> first NRHS columns of the inverse Hilbert matrix. 35*> 36*> The condition number of the Hilbert matrix grows exponentially with 37*> its size, roughly as O(e ** (3.5*N)). Additionally, the inverse 38*> Hilbert matrices beyond a relatively small dimension cannot be 39*> generated exactly without extra precision. Precision is exhausted 40*> when the largest entry in the inverse Hilbert matrix is greater than 41*> 2 to the power of the number of bits in the fraction of the data type 42*> used plus one, which is 24 for single precision. 43*> 44*> In single, the generated solution is exact for N <= 6 and has 45*> small componentwise error for 7 <= N <= 11. 46*> \endverbatim 47* 48* Arguments: 49* ========== 50* 51*> \param[in] N 52*> \verbatim 53*> N is INTEGER 54*> The dimension of the matrix A. 55*> \endverbatim 56*> 57*> \param[in] NRHS 58*> \verbatim 59*> NRHS is INTEGER 60*> The requested number of right-hand sides. 61*> \endverbatim 62*> 63*> \param[out] A 64*> \verbatim 65*> A is COMPLEX array, dimension (LDA, N) 66*> The generated scaled Hilbert matrix. 67*> \endverbatim 68*> 69*> \param[in] LDA 70*> \verbatim 71*> LDA is INTEGER 72*> The leading dimension of the array A. LDA >= N. 73*> \endverbatim 74*> 75*> \param[out] X 76*> \verbatim 77*> X is COMPLEX array, dimension (LDX, NRHS) 78*> The generated exact solutions. Currently, the first NRHS 79*> columns of the inverse Hilbert matrix. 80*> \endverbatim 81*> 82*> \param[in] LDX 83*> \verbatim 84*> LDX is INTEGER 85*> The leading dimension of the array X. LDX >= N. 86*> \endverbatim 87*> 88*> \param[out] B 89*> \verbatim 90*> B is REAL array, dimension (LDB, NRHS) 91*> The generated right-hand sides. Currently, the first NRHS 92*> columns of LCM(1, 2, ..., 2*N-1) * the identity matrix. 93*> \endverbatim 94*> 95*> \param[in] LDB 96*> \verbatim 97*> LDB is INTEGER 98*> The leading dimension of the array B. LDB >= N. 99*> \endverbatim 100*> 101*> \param[out] WORK 102*> \verbatim 103*> WORK is REAL array, dimension (N) 104*> \endverbatim 105*> 106*> \param[out] INFO 107*> \verbatim 108*> INFO is INTEGER 109*> = 0: successful exit 110*> = 1: N is too large; the data is still generated but may not 111*> be not exact. 112*> < 0: if INFO = -i, the i-th argument had an illegal value 113*> \endverbatim 114*> 115*> \param[in] PATH 116*> \verbatim 117*> PATH is CHARACTER*3 118*> The LAPACK path name. 119*> \endverbatim 120* 121* Authors: 122* ======== 123* 124*> \author Univ. of Tennessee 125*> \author Univ. of California Berkeley 126*> \author Univ. of Colorado Denver 127*> \author NAG Ltd. 128* 129*> \ingroup complex16_lin 130* 131* ===================================================================== 132 SUBROUTINE ZLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, 133 $ INFO, PATH) 134* 135* -- LAPACK test routine -- 136* -- LAPACK is a software package provided by Univ. of Tennessee, -- 137* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 138* 139* .. Scalar Arguments .. 140 INTEGER N, NRHS, LDA, LDX, LDB, INFO 141* .. Array Arguments .. 142 DOUBLE PRECISION WORK(N) 143 COMPLEX*16 A(LDA,N), X(LDX, NRHS), B(LDB, NRHS) 144 CHARACTER*3 PATH 145* .. 146* 147* ===================================================================== 148* .. Local Scalars .. 149 INTEGER TM, TI, R 150 INTEGER M 151 INTEGER I, J 152 COMPLEX*16 TMP 153 CHARACTER*2 C2 154* .. 155* .. Parameters .. 156* NMAX_EXACT the largest dimension where the generated data is 157* exact. 158* NMAX_APPROX the largest dimension where the generated data has 159* a small componentwise relative error. 160* ??? complex uses how many bits ??? 161 INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D 162 PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11, SIZE_D = 8) 163* 164* d's are generated from random permutation of those eight elements. 165 COMPLEX*16 d1(8), d2(8), invd1(8), invd2(8) 166 DATA D1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/ 167 DATA D2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/ 168 169 DATA INVD1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0), 170 $ (-.5,-.5),(.5,-.5),(.5,.5)/ 171 DATA INVD2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0), 172 $ (-.5,.5),(.5,.5),(.5,-.5)/ 173* .. 174* .. External Functions 175 EXTERNAL ZLASET, LSAMEN 176 INTRINSIC DBLE 177 LOGICAL LSAMEN 178* .. 179* .. Executable Statements .. 180 C2 = PATH( 2: 3 ) 181* 182* Test the input arguments 183* 184 INFO = 0 185 IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN 186 INFO = -1 187 ELSE IF (NRHS .LT. 0) THEN 188 INFO = -2 189 ELSE IF (LDA .LT. N) THEN 190 INFO = -4 191 ELSE IF (LDX .LT. N) THEN 192 INFO = -6 193 ELSE IF (LDB .LT. N) THEN 194 INFO = -8 195 END IF 196 IF (INFO .LT. 0) THEN 197 CALL XERBLA('ZLAHILB', -INFO) 198 RETURN 199 END IF 200 IF (N .GT. NMAX_EXACT) THEN 201 INFO = 1 202 END IF 203* 204* Compute M = the LCM of the integers [1, 2*N-1]. The largest 205* reasonable N is small enough that integers suffice (up to N = 11). 206 M = 1 207 DO I = 2, (2*N-1) 208 TM = M 209 TI = I 210 R = MOD(TM, TI) 211 DO WHILE (R .NE. 0) 212 TM = TI 213 TI = R 214 R = MOD(TM, TI) 215 END DO 216 M = (M / TI) * I 217 END DO 218* 219* Generate the scaled Hilbert matrix in A 220* If we are testing SY routines, 221* take D1_i = D2_i, else, D1_i = D2_i* 222 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 223 DO J = 1, N 224 DO I = 1, N 225 A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1)) 226 $ * D1(MOD(I,SIZE_D)+1) 227 END DO 228 END DO 229 ELSE 230 DO J = 1, N 231 DO I = 1, N 232 A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1)) 233 $ * D2(MOD(I,SIZE_D)+1) 234 END DO 235 END DO 236 END IF 237* 238* Generate matrix B as simply the first NRHS columns of M * the 239* identity. 240 TMP = DBLE(M) 241 CALL ZLASET('Full', N, NRHS, (0.0D+0,0.0D+0), TMP, B, LDB) 242* 243* Generate the true solutions in X. Because B = the first NRHS 244* columns of M*I, the true solutions are just the first NRHS columns 245* of the inverse Hilbert matrix. 246 WORK(1) = N 247 DO J = 2, N 248 WORK(J) = ( ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1) ) 249 $ * (N +J -1) 250 END DO 251 252* If we are testing SY routines, 253* take D1_i = D2_i, else, D1_i = D2_i* 254 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 255 DO J = 1, NRHS 256 DO I = 1, N 257 X(I, J) = INVD1(MOD(J,SIZE_D)+1) * 258 $ ((WORK(I)*WORK(J)) / (I + J - 1)) 259 $ * INVD1(MOD(I,SIZE_D)+1) 260 END DO 261 END DO 262 ELSE 263 DO J = 1, NRHS 264 DO I = 1, N 265 X(I, J) = INVD2(MOD(J,SIZE_D)+1) * 266 $ ((WORK(I)*WORK(J)) / (I + J - 1)) 267 $ * INVD1(MOD(I,SIZE_D)+1) 268 END DO 269 END DO 270 END IF 271 END 272