1*> \brief \b CLAHILB 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 CLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, 12* INFO, PATH) 13* 14* .. Scalar Arguments .. 15* INTEGER T, N, NRHS, LDA, LDX, LDB, INFO 16* .. Array Arguments .. 17* REAL WORK(N) 18* COMPLEX 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*> CLAHILB 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 NRHS 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*> \date November 2011 130* 131*> \ingroup complex_lin 132* 133* ===================================================================== 134 SUBROUTINE CLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, 135 $ INFO, PATH) 136* 137* -- LAPACK test routine (version 3.4.0) -- 138* -- LAPACK is a software package provided by Univ. of Tennessee, -- 139* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 140* November 2011 141* 142* .. Scalar Arguments .. 143 INTEGER T, N, NRHS, LDA, LDX, LDB, INFO 144* .. Array Arguments .. 145 REAL WORK(N) 146 COMPLEX A(LDA,N), X(LDX, NRHS), B(LDB, NRHS) 147 CHARACTER*3 PATH 148* .. 149* 150* ===================================================================== 151* .. Local Scalars .. 152 INTEGER TM, TI, R 153 INTEGER M 154 INTEGER I, J 155 COMPLEX TMP 156 CHARACTER*2 C2 157* .. 158* .. Parameters .. 159* NMAX_EXACT the largest dimension where the generated data is 160* exact. 161* NMAX_APPROX the largest dimension where the generated data has 162* a small componentwise relative error. 163* ??? complex uses how many bits ??? 164 INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D 165 PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11, SIZE_D = 8) 166* 167* d's are generated from random permuation of those eight elements. 168 COMPLEX D1(8), D2(8), INVD1(8), INVD2(8) 169 DATA D1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/ 170 DATA D2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/ 171 172 DATA INVD1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0), 173 $ (-.5,-.5),(.5,-.5),(.5,.5)/ 174 DATA INVD2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0), 175 $ (-.5,.5),(.5,.5),(.5,-.5)/ 176* .. 177* .. External Functions 178 EXTERNAL CLASET, LSAMEN 179 INTRINSIC REAL 180 LOGICAL LSAMEN 181* .. 182* .. Executable Statements .. 183 C2 = PATH( 2: 3 ) 184* 185* Test the input arguments 186* 187 INFO = 0 188 IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN 189 INFO = -1 190 ELSE IF (NRHS .LT. 0) THEN 191 INFO = -2 192 ELSE IF (LDA .LT. N) THEN 193 INFO = -4 194 ELSE IF (LDX .LT. N) THEN 195 INFO = -6 196 ELSE IF (LDB .LT. N) THEN 197 INFO = -8 198 END IF 199 IF (INFO .LT. 0) THEN 200 CALL XERBLA('CLAHILB', -INFO) 201 RETURN 202 END IF 203 IF (N .GT. NMAX_EXACT) THEN 204 INFO = 1 205 END IF 206* 207* Compute M = the LCM of the integers [1, 2*N-1]. The largest 208* reasonable N is small enough that integers suffice (up to N = 11). 209 M = 1 210 DO I = 2, (2*N-1) 211 TM = M 212 TI = I 213 R = MOD(TM, TI) 214 DO WHILE (R .NE. 0) 215 TM = TI 216 TI = R 217 R = MOD(TM, TI) 218 END DO 219 M = (M / TI) * I 220 END DO 221* 222* Generate the scaled Hilbert matrix in A 223* If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i* 224 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 225 DO J = 1, N 226 DO I = 1, N 227 A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1)) 228 $ * D1(MOD(I,SIZE_D)+1) 229 END DO 230 END DO 231 ELSE 232 DO J = 1, N 233 DO I = 1, N 234 A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1)) 235 $ * D2(MOD(I,SIZE_D)+1) 236 END DO 237 END DO 238 END IF 239* 240* Generate matrix B as simply the first NRHS columns of M * the 241* identity. 242 TMP = REAL(M) 243 CALL CLASET('Full', N, NRHS, (0.0,0.0), TMP, B, LDB) 244* 245* Generate the true solutions in X. Because B = the first NRHS 246* columns of M*I, the true solutions are just the first NRHS columns 247* of the inverse Hilbert matrix. 248 WORK(1) = N 249 DO J = 2, N 250 WORK(J) = ( ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1) ) 251 $ * (N +J -1) 252 END DO 253* 254* If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i* 255 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 256 DO J = 1, NRHS 257 DO I = 1, N 258 X(I, J) = 259 $ INVD1(MOD(J,SIZE_D)+1) * 260 $ ((WORK(I)*WORK(J)) / (I + J - 1)) 261 $ * INVD1(MOD(I,SIZE_D)+1) 262 END DO 263 END DO 264 ELSE 265 DO J = 1, NRHS 266 DO I = 1, N 267 X(I, J) = 268 $ INVD2(MOD(J,SIZE_D)+1) * 269 $ ((WORK(I)*WORK(J)) / (I + J - 1)) 270 $ * INVD1(MOD(I,SIZE_D)+1) 271 END DO 272 END DO 273 END IF 274 END 275 276