LAPACK  3.7.1
LAPACK: Linear Algebra PACKage
clahilb.f
Go to the documentation of this file.
1 *> \brief \b CLAHILB
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE CLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
12 * INFO, PATH)
13 *
14 * .. Scalar Arguments ..
15 * INTEGER N, NRHS, LDA, LDX, LDB, INFO
16 * .. Array Arguments ..
17 * REAL WORK(N)
18 * COMPLEX A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
19 * CHARACTER*3 PATH
20 * ..
21 *
22 *
23 *> \par Purpose:
24 * =============
25 *>
26 *> \verbatim
27 *>
28 *> CLAHILB generates an N by N scaled Hilbert matrix in A along with
29 *> NRHS right-hand sides in B and solutions in X such that A*X=B.
30 *>
31 *> The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
32 *> entries are integers. The right-hand sides are the first NRHS
33 *> columns of M * the identity matrix, and the solutions are the
34 *> first NRHS columns of the inverse Hilbert matrix.
35 *>
36 *> The condition number of the Hilbert matrix grows exponentially with
37 *> its size, roughly as O(e ** (3.5*N)). Additionally, the inverse
38 *> Hilbert matrices beyond a relatively small dimension cannot be
39 *> generated exactly without extra precision. Precision is exhausted
40 *> when the largest entry in the inverse Hilbert matrix is greater than
41 *> 2 to the power of the number of bits in the fraction of the data type
42 *> used plus one, which is 24 for single precision.
43 *>
44 *> In single, the generated solution is exact for N <= 6 and has
45 *> small componentwise error for 7 <= N <= 11.
46 *> \endverbatim
47 *
48 * Arguments:
49 * ==========
50 *
51 *> \param[in] N
52 *> \verbatim
53 *> N is INTEGER
54 *> The dimension of the matrix A.
55 *> \endverbatim
56 *>
57 *> \param[in] NRHS
58 *> \verbatim
59 *> NRHS is INTEGER
60 *> The requested number of right-hand sides.
61 *> \endverbatim
62 *>
63 *> \param[out] A
64 *> \verbatim
65 *> A is COMPLEX array, dimension (LDA, N)
66 *> The generated scaled Hilbert matrix.
67 *> \endverbatim
68 *>
69 *> \param[in] LDA
70 *> \verbatim
71 *> LDA is INTEGER
72 *> The leading dimension of the array A. LDA >= N.
73 *> \endverbatim
74 *>
75 *> \param[out] X
76 *> \verbatim
77 *> X is COMPLEX array, dimension (LDX, NRHS)
78 *> The generated exact solutions. Currently, the first NRHS
79 *> columns of the inverse Hilbert matrix.
80 *> \endverbatim
81 *>
82 *> \param[in] LDX
83 *> \verbatim
84 *> LDX is INTEGER
85 *> The leading dimension of the array X. LDX >= N.
86 *> \endverbatim
87 *>
88 *> \param[out] B
89 *> \verbatim
90 *> B is REAL array, dimension (LDB, NRHS)
91 *> The generated right-hand sides. Currently, the first NRHS
92 *> columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
93 *> \endverbatim
94 *>
95 *> \param[in] LDB
96 *> \verbatim
97 *> LDB is INTEGER
98 *> The leading dimension of the array B. LDB >= N.
99 *> \endverbatim
100 *>
101 *> \param[out] WORK
102 *> \verbatim
103 *> WORK is REAL array, dimension (N)
104 *> \endverbatim
105 *>
106 *> \param[out] INFO
107 *> \verbatim
108 *> INFO is INTEGER
109 *> = 0: successful exit
110 *> = 1: N is too large; the data is still generated but may not
111 *> be not exact.
112 *> < 0: if INFO = -i, the i-th argument had an illegal value
113 *> \endverbatim
114 *>
115 *> \param[in] PATH
116 *> \verbatim
117 *> PATH is CHARACTER*3
118 *> The LAPACK path name.
119 *> \endverbatim
120 *
121 * Authors:
122 * ========
123 *
124 *> \author Univ. of Tennessee
125 *> \author Univ. of California Berkeley
126 *> \author Univ. of Colorado Denver
127 *> \author NAG Ltd.
128 *
129 *> \date June 2017
130 *
131 *> \ingroup complex_matgen
132 *
133 * =====================================================================
134  SUBROUTINE clahilb( N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
135  $ INFO, PATH)
136 *
137 * -- LAPACK test routine (version 3.7.1) --
138 * -- LAPACK is a software package provided by Univ. of Tennessee, --
139 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140 * June 2017
141 *
142 * .. Scalar Arguments ..
143  INTEGER N, NRHS, LDA, LDX, LDB, INFO
144 * .. Array Arguments ..
145  REAL WORK(n)
146  COMPLEX A(lda,n), X(ldx, nrhs), B(ldb, nrhs)
147  CHARACTER*3 PATH
148 * ..
149 *
150 * =====================================================================
151 * .. Local Scalars ..
152  INTEGER TM, TI, R
153  INTEGER M
154  INTEGER I, J
155  COMPLEX TMP
156  CHARACTER*2 C2
157 * ..
158 * .. Parameters ..
159 * NMAX_EXACT the largest dimension where the generated data is
160 * exact.
161 * NMAX_APPROX the largest dimension where the generated data has
162 * a small componentwise relative error.
163 * ??? complex uses how many bits ???
164  INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D
165  parameter(nmax_exact = 6, nmax_approx = 11, size_d = 8)
166 *
167 * d's are generated from random permuation of those eight elements.
168  COMPLEX D1(8), D2(8), INVD1(8), INVD2(8)
169  DATA d1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/
170  DATA d2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/
171 
172  DATA invd1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0),
173  $ (-.5,-.5),(.5,-.5),(.5,.5)/
174  DATA invd2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0),
175  $ (-.5,.5),(.5,.5),(.5,-.5)/
176 * ..
177 * .. External Functions
178  EXTERNAL claset, lsamen
179  INTRINSIC real
180  LOGICAL LSAMEN
181 * ..
182 * .. Executable Statements ..
183  c2 = path( 2: 3 )
184 *
185 * Test the input arguments
186 *
187  info = 0
188  IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
189  info = -1
190  ELSE IF (nrhs .LT. 0) THEN
191  info = -2
192  ELSE IF (lda .LT. n) THEN
193  info = -4
194  ELSE IF (ldx .LT. n) THEN
195  info = -6
196  ELSE IF (ldb .LT. n) THEN
197  info = -8
198  END IF
199  IF (info .LT. 0) THEN
200  CALL xerbla('CLAHILB', -info)
201  RETURN
202  END IF
203  IF (n .GT. nmax_exact) THEN
204  info = 1
205  END IF
206 *
207 * Compute M = the LCM of the integers [1, 2*N-1]. The largest
208 * reasonable N is small enough that integers suffice (up to N = 11).
209  m = 1
210  DO i = 2, (2*n-1)
211  tm = m
212  ti = i
213  r = mod(tm, ti)
214  DO WHILE (r .NE. 0)
215  tm = ti
216  ti = r
217  r = mod(tm, ti)
218  END DO
219  m = (m / ti) * i
220  END DO
221 *
222 * Generate the scaled Hilbert matrix in A
223 * If we are testing SY routines, take
224 * D1_i = D2_i, else, D1_i = D2_i*
225  IF ( lsamen( 2, c2, 'SY' ) ) THEN
226  DO j = 1, n
227  DO i = 1, n
228  a(i, j) = d1(mod(j,size_d)+1) * (REAL(M) / (i + j - 1))
229  $ * d1(mod(i,size_d)+1)
230  END DO
231  END DO
232  ELSE
233  DO j = 1, n
234  DO i = 1, n
235  a(i, j) = d1(mod(j,size_d)+1) * (REAL(M) / (i + j - 1))
236  $ * d2(mod(i,size_d)+1)
237  END DO
238  END DO
239  END IF
240 *
241 * Generate matrix B as simply the first NRHS columns of M * the
242 * identity.
243  tmp = REAL(m)
244  CALL claset('Full', n, nrhs, (0.0,0.0), tmp, b, ldb)
245 *
246 * Generate the true solutions in X. Because B = the first NRHS
247 * columns of M*I, the true solutions are just the first NRHS columns
248 * of the inverse Hilbert matrix.
249  work(1) = n
250  DO j = 2, n
251  work(j) = ( ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1) )
252  $ * (n +j -1)
253  END DO
254 
255 * If we are testing SY routines,
256 * take D1_i = D2_i, else, D1_i = D2_i*
257  IF ( lsamen( 2, c2, 'SY' ) ) THEN
258  DO j = 1, nrhs
259  DO i = 1, n
260  x(i, j) =
261  $ invd1(mod(j,size_d)+1) *
262  $ ((work(i)*work(j)) / (i + j - 1))
263  $ * invd1(mod(i,size_d)+1)
264  END DO
265  END DO
266  ELSE
267  DO j = 1, nrhs
268  DO i = 1, n
269  x(i, j) =
270  $ invd2(mod(j,size_d)+1) *
271  $ ((work(i)*work(j)) / (i + j - 1))
272  $ * invd1(mod(i,size_d)+1)
273  END DO
274  END DO
275  END IF
276  END
277 
subroutine clahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO, PATH)
CLAHILB
Definition: clahilb.f:136
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
logical function lsamen(N, CA, CB)
LSAMEN
Definition: lsamen.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62