00001 SUBROUTINE CLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, 00002 $ INFO, PATH) 00003 ! 00004 ! -- LAPACK auxiliary test routine (version 3.2.2) -- 00005 ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00006 ! Courant Institute, Argonne National Lab, and Rice University 00007 * June 2010 00008 ! 00009 ! David Vu <dtv@cs.berkeley.edu> 00010 ! Yozo Hida <yozo@cs.berkeley.edu> 00011 ! Jason Riedy <ejr@cs.berkeley.edu> 00012 ! D. Halligan <dhalligan@berkeley.edu> 00013 ! 00014 IMPLICIT NONE 00015 ! .. Scalar Arguments .. 00016 INTEGER N, NRHS, LDA, LDX, LDB, INFO 00017 ! .. Array Arguments .. 00018 REAL WORK(N) 00019 COMPLEX A(LDA,N), X(LDX, NRHS), B(LDB, NRHS) 00020 CHARACTER*3 PATH 00021 ! .. 00022 ! 00023 ! Purpose 00024 ! ======= 00025 ! 00026 ! CLAHILB generates an N by N scaled Hilbert matrix in A along with 00027 ! NRHS right-hand sides in B and solutions in X such that A*X=B. 00028 ! 00029 ! The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all 00030 ! entries are integers. The right-hand sides are the first NRHS 00031 ! columns of M * the identity matrix, and the solutions are the 00032 ! first NRHS columns of the inverse Hilbert matrix. 00033 ! 00034 ! The condition number of the Hilbert matrix grows exponentially with 00035 ! its size, roughly as O(e ** (3.5*N)). Additionally, the inverse 00036 ! Hilbert matrices beyond a relatively small dimension cannot be 00037 ! generated exactly without extra precision. Precision is exhausted 00038 ! when the largest entry in the inverse Hilbert matrix is greater than 00039 ! 2 to the power of the number of bits in the fraction of the data type 00040 ! used plus one, which is 24 for single precision. 00041 ! 00042 ! In single, the generated solution is exact for N <= 6 and has 00043 ! small componentwise error for 7 <= N <= 11. 00044 ! 00045 ! Arguments 00046 ! ========= 00047 ! 00048 ! N (input) INTEGER 00049 ! The dimension of the matrix A. 00050 ! 00051 ! NRHS (input) INTEGER 00052 ! The requested number of right-hand sides. 00053 ! 00054 ! A (output) COMPLEX array, dimension (LDA, N) 00055 ! The generated scaled Hilbert matrix. 00056 ! 00057 ! LDA (input) INTEGER 00058 ! The leading dimension of the array A. LDA >= N. 00059 ! 00060 ! X (output) COMPLEX array, dimension (LDX, NRHS) 00061 ! The generated exact solutions. Currently, the first NRHS 00062 ! columns of the inverse Hilbert matrix. 00063 ! 00064 ! LDX (input) INTEGER 00065 ! The leading dimension of the array X. LDX >= N. 00066 ! 00067 ! B (output) REAL array, dimension (LDB, NRHS) 00068 ! The generated right-hand sides. Currently, the first NRHS 00069 ! columns of LCM(1, 2, ..., 2*N-1) * the identity matrix. 00070 ! 00071 ! LDB (input) INTEGER 00072 ! The leading dimension of the array B. LDB >= N. 00073 ! 00074 ! WORK (workspace) REAL array, dimension (N) 00075 ! 00076 ! 00077 ! INFO (output) INTEGER 00078 ! = 0: successful exit 00079 ! = 1: N is too large; the data is still generated but may not 00080 ! be not exact. 00081 ! < 0: if INFO = -i, the i-th argument had an illegal value 00082 ! 00083 ! ===================================================================== 00084 00085 ! .. Local Scalars .. 00086 INTEGER TM, TI, R 00087 INTEGER M 00088 INTEGER I, J 00089 COMPLEX TMP 00090 CHARACTER*2 C2 00091 00092 ! .. Parameters .. 00093 ! NMAX_EXACT the largest dimension where the generated data is 00094 ! exact. 00095 ! NMAX_APPROX the largest dimension where the generated data has 00096 ! a small componentwise relative error. 00097 ! ??? complex uses how many bits ??? 00098 INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D 00099 PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11, SIZE_D = 8) 00100 00101 ! d's are generated from random permuation of those eight elements. 00102 COMPLEX D1(8), D2(8), INVD1(8), INVD2(8) 00103 DATA D1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/ 00104 DATA D2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/ 00105 00106 DATA INVD1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0), 00107 $ (-.5,-.5),(.5,-.5),(.5,.5)/ 00108 DATA INVD2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0), 00109 $ (-.5,.5),(.5,.5),(.5,-.5)/ 00110 00111 ! .. 00112 ! .. External Functions 00113 EXTERNAL CLASET, LSAMEN 00114 INTRINSIC REAL 00115 LOGICAL LSAMEN 00116 ! .. 00117 ! .. Executable Statements .. 00118 C2 = PATH( 2: 3 ) 00119 ! 00120 ! Test the input arguments 00121 ! 00122 INFO = 0 00123 IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN 00124 INFO = -1 00125 ELSE IF (NRHS .LT. 0) THEN 00126 INFO = -2 00127 ELSE IF (LDA .LT. N) THEN 00128 INFO = -4 00129 ELSE IF (LDX .LT. N) THEN 00130 INFO = -6 00131 ELSE IF (LDB .LT. N) THEN 00132 INFO = -8 00133 END IF 00134 IF (INFO .LT. 0) THEN 00135 CALL XERBLA('CLAHILB', -INFO) 00136 RETURN 00137 END IF 00138 IF (N .GT. NMAX_EXACT) THEN 00139 INFO = 1 00140 END IF 00141 00142 ! Compute M = the LCM of the integers [1, 2*N-1]. The largest 00143 ! reasonable N is small enough that integers suffice (up to N = 11). 00144 M = 1 00145 DO I = 2, (2*N-1) 00146 TM = M 00147 TI = I 00148 R = MOD(TM, TI) 00149 DO WHILE (R .NE. 0) 00150 TM = TI 00151 TI = R 00152 R = MOD(TM, TI) 00153 END DO 00154 M = (M / TI) * I 00155 END DO 00156 00157 ! Generate the scaled Hilbert matrix in A 00158 ! If we are testing SY routines, take 00159 ! D1_i = D2_i, else, D1_i = D2_i* 00160 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 00161 DO J = 1, N 00162 DO I = 1, N 00163 A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1)) 00164 $ * D1(MOD(I,SIZE_D)+1) 00165 END DO 00166 END DO 00167 ELSE 00168 DO J = 1, N 00169 DO I = 1, N 00170 A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1)) 00171 $ * D2(MOD(I,SIZE_D)+1) 00172 END DO 00173 END DO 00174 END IF 00175 00176 ! Generate matrix B as simply the first NRHS columns of M * the 00177 ! identity. 00178 TMP = REAL(M) 00179 CALL CLASET('Full', N, NRHS, (0.0,0.0), TMP, B, LDB) 00180 00181 ! Generate the true solutions in X. Because B = the first NRHS 00182 ! columns of M*I, the true solutions are just the first NRHS columns 00183 ! of the inverse Hilbert matrix. 00184 WORK(1) = N 00185 DO J = 2, N 00186 WORK(J) = ( ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1) ) 00187 $ * (N +J -1) 00188 END DO 00189 00190 ! If we are testing SY routines, 00191 ! take D1_i = D2_i, else, D1_i = D2_i* 00192 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 00193 DO J = 1, NRHS 00194 DO I = 1, N 00195 X(I, J) = 00196 $ INVD1(MOD(J,SIZE_D)+1) * 00197 $ ((WORK(I)*WORK(J)) / (I + J - 1)) 00198 $ * INVD1(MOD(I,SIZE_D)+1) 00199 END DO 00200 END DO 00201 ELSE 00202 DO J = 1, NRHS 00203 DO I = 1, N 00204 X(I, J) = 00205 $ INVD2(MOD(J,SIZE_D)+1) * 00206 $ ((WORK(I)*WORK(J)) / (I + J - 1)) 00207 $ * INVD1(MOD(I,SIZE_D)+1) 00208 END DO 00209 END DO 00210 END IF 00211 END 00212
1.7.2