1 SUBROUTINE DLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
2 !
3 ! -- LAPACK auxiliary test routine (version 3.0) --
4 ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 ! Courant Institute, Argonne National Lab, and Rice University
6 ! 28 August, 2006
7 !
8 ! David Vu <dtv@cs.berkeley.edu>
9 ! Yozo Hida <yozo@cs.berkeley.edu>
10 ! Jason Riedy <ejr@cs.berkeley.edu>
11 ! D. Halligan <dhalligan@berkeley.edu>
12 !
13 IMPLICIT NONE
14 ! .. Scalar Arguments ..
15 INTEGER N, NRHS, LDA, LDX, LDB, INFO
16 ! .. Array Arguments ..
17 DOUBLE PRECISION A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
18 ! ..
19 !
20 ! Purpose
21 ! =======
22 !
23 ! DLAHILB generates an N by N scaled Hilbert matrix in A along with
24 ! NRHS right-hand sides in B and solutions in X such that A*X=B.
25 !
26 ! The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
27 ! entries are integers. The right-hand sides are the first NRHS
28 ! columns of M * the identity matrix, and the solutions are the
29 ! first NRHS columns of the inverse Hilbert matrix.
30 !
31 ! The condition number of the Hilbert matrix grows exponentially with
32 ! its size, roughly as O(e ** (3.5*N)). Additionally, the inverse
33 ! Hilbert matrices beyond a relatively small dimension cannot be
34 ! generated exactly without extra precision. Precision is exhausted
35 ! when the largest entry in the inverse Hilbert matrix is greater than
36 ! 2 to the power of the number of bits in the fraction of the data type
37 ! used plus one, which is 24 for single precision.
38 !
39 ! In single, the generated solution is exact for N <= 6 and has
40 ! small componentwise error for 7 <= N <= 11.
41 !
42 ! Arguments
43 ! =========
44 !
45 ! N (input) INTEGER
46 ! The dimension of the matrix A.
47 !
48 ! NRHS (input) NRHS
49 ! The requested number of right-hand sides.
50 !
51 ! A (output) DOUBLE PRECISION array, dimension (LDA, N)
52 ! The generated scaled Hilbert matrix.
53 !
54 ! LDA (input) INTEGER
55 ! The leading dimension of the array A. LDA >= N.
56 !
57 ! X (output) DOUBLE PRECISION array, dimension (LDX, NRHS)
58 ! The generated exact solutions. Currently, the first NRHS
59 ! columns of the inverse Hilbert matrix.
60 !
61 ! LDX (input) INTEGER
62 ! The leading dimension of the array X. LDX >= N.
63 !
64 ! B (output) DOUBLE PRECISION array, dimension (LDB, NRHS)
65 ! The generated right-hand sides. Currently, the first NRHS
66 ! columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
67 !
68 ! LDB (input) INTEGER
69 ! The leading dimension of the array B. LDB >= N.
70 !
71 ! WORK (workspace) DOUBLE PRECISION array, dimension (N)
72 !
73 !
74 ! INFO (output) INTEGER
75 ! = 0: successful exit
76 ! = 1: N is too large; the data is still generated but may not
77 ! be not exact.
78 ! < 0: if INFO = -i, the i-th argument had an illegal value
79 !
80 ! =====================================================================
81
82 ! .. Local Scalars ..
83 INTEGER TM, TI, R
84 INTEGER M
85 INTEGER I, J
86 COMPLEX*16 TMP
87
88 ! .. Parameters ..
89 ! NMAX_EXACT the largest dimension where the generated data is
90 ! exact.
91 ! NMAX_APPROX the largest dimension where the generated data has
92 ! a small componentwise relative error.
93 INTEGER NMAX_EXACT, NMAX_APPROX
94 PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11)
95
96 ! ..
97 ! .. External Functions
98 EXTERNAL DLASET
99 INTRINSIC DBLE
100 ! ..
101 ! .. Executable Statements ..
102 !
103 ! Test the input arguments
104 !
105 INFO = 0
106 IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
107 INFO = -1
108 ELSE IF (NRHS .LT. 0) THEN
109 INFO = -2
110 ELSE IF (LDA .LT. N) THEN
111 INFO = -4
112 ELSE IF (LDX .LT. N) THEN
113 INFO = -6
114 ELSE IF (LDB .LT. N) THEN
115 INFO = -8
116 END IF
117 IF (INFO .LT. 0) THEN
118 CALL XERBLA('DLAHILB', -INFO)
119 RETURN
120 END IF
121 IF (N .GT. NMAX_EXACT) THEN
122 INFO = 1
123 END IF
124
125 ! Compute M = the LCM of the integers [1, 2*N-1]. The largest
126 ! reasonable N is small enough that integers suffice (up to N = 11).
127 M = 1
128 DO I = 2, (2*N-1)
129 TM = M
130 TI = I
131 R = MOD(TM, TI)
132 DO WHILE (R .NE. 0)
133 TM = TI
134 TI = R
135 R = MOD(TM, TI)
136 END DO
137 M = (M / TI) * I
138 END DO
139
140 ! Generate the scaled Hilbert matrix in A
141 DO J = 1, N
142 DO I = 1, N
143 A(I, J) = DBLE(M) / (I + J - 1)
144 END DO
145 END DO
146
147 ! Generate matrix B as simply the first NRHS columns of M * the
148 ! identity.
149 TMP = DBLE(M)
150 CALL DLASET('Full', N, NRHS, 0.0D+0, TMP, B, LDB)
151
152 ! Generate the true solutions in X. Because B = the first NRHS
153 ! columns of M*I, the true solutions are just the first NRHS columns
154 ! of the inverse Hilbert matrix.
155 WORK(1) = N
156 DO J = 2, N
157 WORK(J) = ( ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1) )
158 $ * (N +J -1)
159 END DO
160
161 DO J = 1, NRHS
162 DO I = 1, N
163 X(I, J) = (WORK(I)*WORK(J)) / (I + J - 1)
164 END DO
165 END DO
166
167 END
168
2 !
3 ! -- LAPACK auxiliary test routine (version 3.0) --
4 ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 ! Courant Institute, Argonne National Lab, and Rice University
6 ! 28 August, 2006
7 !
8 ! David Vu <dtv@cs.berkeley.edu>
9 ! Yozo Hida <yozo@cs.berkeley.edu>
10 ! Jason Riedy <ejr@cs.berkeley.edu>
11 ! D. Halligan <dhalligan@berkeley.edu>
12 !
13 IMPLICIT NONE
14 ! .. Scalar Arguments ..
15 INTEGER N, NRHS, LDA, LDX, LDB, INFO
16 ! .. Array Arguments ..
17 DOUBLE PRECISION A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
18 ! ..
19 !
20 ! Purpose
21 ! =======
22 !
23 ! DLAHILB generates an N by N scaled Hilbert matrix in A along with
24 ! NRHS right-hand sides in B and solutions in X such that A*X=B.
25 !
26 ! The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
27 ! entries are integers. The right-hand sides are the first NRHS
28 ! columns of M * the identity matrix, and the solutions are the
29 ! first NRHS columns of the inverse Hilbert matrix.
30 !
31 ! The condition number of the Hilbert matrix grows exponentially with
32 ! its size, roughly as O(e ** (3.5*N)). Additionally, the inverse
33 ! Hilbert matrices beyond a relatively small dimension cannot be
34 ! generated exactly without extra precision. Precision is exhausted
35 ! when the largest entry in the inverse Hilbert matrix is greater than
36 ! 2 to the power of the number of bits in the fraction of the data type
37 ! used plus one, which is 24 for single precision.
38 !
39 ! In single, the generated solution is exact for N <= 6 and has
40 ! small componentwise error for 7 <= N <= 11.
41 !
42 ! Arguments
43 ! =========
44 !
45 ! N (input) INTEGER
46 ! The dimension of the matrix A.
47 !
48 ! NRHS (input) NRHS
49 ! The requested number of right-hand sides.
50 !
51 ! A (output) DOUBLE PRECISION array, dimension (LDA, N)
52 ! The generated scaled Hilbert matrix.
53 !
54 ! LDA (input) INTEGER
55 ! The leading dimension of the array A. LDA >= N.
56 !
57 ! X (output) DOUBLE PRECISION array, dimension (LDX, NRHS)
58 ! The generated exact solutions. Currently, the first NRHS
59 ! columns of the inverse Hilbert matrix.
60 !
61 ! LDX (input) INTEGER
62 ! The leading dimension of the array X. LDX >= N.
63 !
64 ! B (output) DOUBLE PRECISION array, dimension (LDB, NRHS)
65 ! The generated right-hand sides. Currently, the first NRHS
66 ! columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
67 !
68 ! LDB (input) INTEGER
69 ! The leading dimension of the array B. LDB >= N.
70 !
71 ! WORK (workspace) DOUBLE PRECISION array, dimension (N)
72 !
73 !
74 ! INFO (output) INTEGER
75 ! = 0: successful exit
76 ! = 1: N is too large; the data is still generated but may not
77 ! be not exact.
78 ! < 0: if INFO = -i, the i-th argument had an illegal value
79 !
80 ! =====================================================================
81
82 ! .. Local Scalars ..
83 INTEGER TM, TI, R
84 INTEGER M
85 INTEGER I, J
86 COMPLEX*16 TMP
87
88 ! .. Parameters ..
89 ! NMAX_EXACT the largest dimension where the generated data is
90 ! exact.
91 ! NMAX_APPROX the largest dimension where the generated data has
92 ! a small componentwise relative error.
93 INTEGER NMAX_EXACT, NMAX_APPROX
94 PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11)
95
96 ! ..
97 ! .. External Functions
98 EXTERNAL DLASET
99 INTRINSIC DBLE
100 ! ..
101 ! .. Executable Statements ..
102 !
103 ! Test the input arguments
104 !
105 INFO = 0
106 IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
107 INFO = -1
108 ELSE IF (NRHS .LT. 0) THEN
109 INFO = -2
110 ELSE IF (LDA .LT. N) THEN
111 INFO = -4
112 ELSE IF (LDX .LT. N) THEN
113 INFO = -6
114 ELSE IF (LDB .LT. N) THEN
115 INFO = -8
116 END IF
117 IF (INFO .LT. 0) THEN
118 CALL XERBLA('DLAHILB', -INFO)
119 RETURN
120 END IF
121 IF (N .GT. NMAX_EXACT) THEN
122 INFO = 1
123 END IF
124
125 ! Compute M = the LCM of the integers [1, 2*N-1]. The largest
126 ! reasonable N is small enough that integers suffice (up to N = 11).
127 M = 1
128 DO I = 2, (2*N-1)
129 TM = M
130 TI = I
131 R = MOD(TM, TI)
132 DO WHILE (R .NE. 0)
133 TM = TI
134 TI = R
135 R = MOD(TM, TI)
136 END DO
137 M = (M / TI) * I
138 END DO
139
140 ! Generate the scaled Hilbert matrix in A
141 DO J = 1, N
142 DO I = 1, N
143 A(I, J) = DBLE(M) / (I + J - 1)
144 END DO
145 END DO
146
147 ! Generate matrix B as simply the first NRHS columns of M * the
148 ! identity.
149 TMP = DBLE(M)
150 CALL DLASET('Full', N, NRHS, 0.0D+0, TMP, B, LDB)
151
152 ! Generate the true solutions in X. Because B = the first NRHS
153 ! columns of M*I, the true solutions are just the first NRHS columns
154 ! of the inverse Hilbert matrix.
155 WORK(1) = N
156 DO J = 2, N
157 WORK(J) = ( ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1) )
158 $ * (N +J -1)
159 END DO
160
161 DO J = 1, NRHS
162 DO I = 1, N
163 X(I, J) = (WORK(I)*WORK(J)) / (I + J - 1)
164 END DO
165 END DO
166
167 END
168