00001 SUBROUTINE ZLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, 00002 $ INFO, PATH) 00003 ! 00004 ! -- LAPACK auxiliary test routine (version 3.0) -- 00005 ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00006 ! Courant Institute, Argonne National Lab, and Rice University 00007 ! 28 August, 2006 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 DOUBLE PRECISION WORK(N) 00019 COMPLEX*16 A(LDA,N), X(LDX, NRHS), B(LDB, NRHS) 00020 CHARACTER*3 PATH 00021 ! .. 00022 ! 00023 ! Purpose 00024 ! ======= 00025 ! 00026 ! ZLAHILB 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) NRHS 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*16 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*16 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 ! .. External Functions 00112 EXTERNAL ZLASET, LSAMEN 00113 INTRINSIC DBLE 00114 LOGICAL LSAMEN 00115 ! .. 00116 ! .. Executable Statements .. 00117 C2 = PATH( 2: 3 ) 00118 ! 00119 ! Test the input arguments 00120 ! 00121 INFO = 0 00122 IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN 00123 INFO = -1 00124 ELSE IF (NRHS .LT. 0) THEN 00125 INFO = -2 00126 ELSE IF (LDA .LT. N) THEN 00127 INFO = -4 00128 ELSE IF (LDX .LT. N) THEN 00129 INFO = -6 00130 ELSE IF (LDB .LT. N) THEN 00131 INFO = -8 00132 END IF 00133 IF (INFO .LT. 0) THEN 00134 CALL XERBLA('ZLAHILB', -INFO) 00135 RETURN 00136 END IF 00137 IF (N .GT. NMAX_EXACT) THEN 00138 INFO = 1 00139 END IF 00140 00141 ! Compute M = the LCM of the integers [1, 2*N-1]. The largest 00142 ! reasonable N is small enough that integers suffice (up to N = 11). 00143 M = 1 00144 DO I = 2, (2*N-1) 00145 TM = M 00146 TI = I 00147 R = MOD(TM, TI) 00148 DO WHILE (R .NE. 0) 00149 TM = TI 00150 TI = R 00151 R = MOD(TM, TI) 00152 END DO 00153 M = (M / TI) * I 00154 END DO 00155 00156 ! Generate the scaled Hilbert matrix in A 00157 ! If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i* 00158 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 00159 DO J = 1, N 00160 DO I = 1, N 00161 A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1)) 00162 $ * D1(MOD(I,SIZE_D)+1) 00163 END DO 00164 END DO 00165 ELSE 00166 DO J = 1, N 00167 DO I = 1, N 00168 A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1)) 00169 $ * D2(MOD(I,SIZE_D)+1) 00170 END DO 00171 END DO 00172 END IF 00173 00174 ! Generate matrix B as simply the first NRHS columns of M * the 00175 ! identity. 00176 TMP = DBLE(M) 00177 CALL ZLASET('Full', N, NRHS, (0.0D+0,0.0D+0), TMP, B, LDB) 00178 00179 ! Generate the true solutions in X. Because B = the first NRHS 00180 ! columns of M*I, the true solutions are just the first NRHS columns 00181 ! of the inverse Hilbert matrix. 00182 WORK(1) = N 00183 DO J = 2, N 00184 WORK(J) = ( ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1) ) 00185 $ * (N +J -1) 00186 END DO 00187 00188 ! If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i* 00189 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 00190 DO J = 1, NRHS 00191 DO I = 1, N 00192 X(I, J) = INVD1(MOD(J,SIZE_D)+1) * 00193 $ ((WORK(I)*WORK(J)) / (I + J - 1)) 00194 $ * INVD1(MOD(I,SIZE_D)+1) 00195 END DO 00196 END DO 00197 ELSE 00198 DO J = 1, NRHS 00199 DO I = 1, N 00200 X(I, J) = INVD2(MOD(J,SIZE_D)+1) * 00201 $ ((WORK(I)*WORK(J)) / (I + J - 1)) 00202 $ * INVD1(MOD(I,SIZE_D)+1) 00203 END DO 00204 END DO 00205 END IF 00206 END
1.7.2