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