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