LAPACK  3.8.0 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 November 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.8.0) --
138 * -- LAPACK is a software package provided by Univ. of Tennessee, --
139 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140 * November 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 Subroutines ..
178  EXTERNAL xerbla
179 * ..
180 * .. External Functions
181  EXTERNAL claset, lsamen
182  INTRINSIC real
183  LOGICAL LSAMEN
184 * ..
185 * .. Executable Statements ..
186  c2 = path( 2: 3 )
187 *
188 * Test the input arguments
189 *
190  info = 0
191  IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
192  info = -1
193  ELSE IF (nrhs .LT. 0) THEN
194  info = -2
195  ELSE IF (lda .LT. n) THEN
196  info = -4
197  ELSE IF (ldx .LT. n) THEN
198  info = -6
199  ELSE IF (ldb .LT. n) THEN
200  info = -8
201  END IF
202  IF (info .LT. 0) THEN
203  CALL xerbla('CLAHILB', -info)
204  RETURN
205  END IF
206  IF (n .GT. nmax_exact) THEN
207  info = 1
208  END IF
209 *
210 * Compute M = the LCM of the integers [1, 2*N-1]. The largest
211 * reasonable N is small enough that integers suffice (up to N = 11).
212  m = 1
213  DO i = 2, (2*n-1)
214  tm = m
215  ti = i
216  r = mod(tm, ti)
217  DO WHILE (r .NE. 0)
218  tm = ti
219  ti = r
220  r = mod(tm, ti)
221  END DO
222  m = (m / ti) * i
223  END DO
224 *
225 * Generate the scaled Hilbert matrix in A
226 * If we are testing SY routines, take
227 * D1_i = D2_i, else, D1_i = D2_i*
228  IF ( lsamen( 2, c2, 'SY' ) ) THEN
229  DO j = 1, n
230  DO i = 1, n
231  a(i, j) = d1(mod(j,size_d)+1) * (REAL(M) / (i + j - 1))
232  \$ * d1(mod(i,size_d)+1)
233  END DO
234  END DO
235  ELSE
236  DO j = 1, n
237  DO i = 1, n
238  a(i, j) = d1(mod(j,size_d)+1) * (REAL(M) / (i + j - 1))
239  \$ * d2(mod(i,size_d)+1)
240  END DO
241  END DO
242  END IF
243 *
244 * Generate matrix B as simply the first NRHS columns of M * the
245 * identity.
246  tmp = REAL(m)
247  CALL claset('Full', n, nrhs, (0.0,0.0), tmp, b, ldb)
248 *
249 * Generate the true solutions in X. Because B = the first NRHS
250 * columns of M*I, the true solutions are just the first NRHS columns
251 * of the inverse Hilbert matrix.
252  work(1) = n
253  DO j = 2, n
254  work(j) = ( ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1) )
255  \$ * (n +j -1)
256  END DO
257
258 * If we are testing SY routines,
259 * take D1_i = D2_i, else, D1_i = D2_i*
260  IF ( lsamen( 2, c2, 'SY' ) ) THEN
261  DO j = 1, nrhs
262  DO i = 1, n
263  x(i, j) =
264  \$ invd1(mod(j,size_d)+1) *
265  \$ ((work(i)*work(j)) / (i + j - 1))
266  \$ * invd1(mod(i,size_d)+1)
267  END DO
268  END DO
269  ELSE
270  DO j = 1, nrhs
271  DO i = 1, n
272  x(i, j) =
273  \$ invd2(mod(j,size_d)+1) *
274  \$ ((work(i)*work(j)) / (i + j - 1))
275  \$ * invd1(mod(i,size_d)+1)
276  END DO
277  END DO
278  END IF
279  END
280
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
subroutine clahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO, PATH)
CLAHILB
Definition: clahilb.f:136
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsamen(N, CA, CB)
LSAMEN
Definition: lsamen.f:76