LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
slahilb.f
Go to the documentation of this file.
1 *> \brief \b SLAHILB
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 SLAHILB( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
12 *
13 * .. Scalar Arguments ..
14 * INTEGER N, NRHS, LDA, LDX, LDB, INFO
15 * .. Array Arguments ..
16 * REAL A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
17 * ..
18 *
19 *
20 *> \par Purpose:
21 * =============
22 *>
23 *> \verbatim
24 *>
25 *> SLAHILB generates an N by N scaled Hilbert matrix in A along with
26 *> NRHS right-hand sides in B and solutions in X such that A*X=B.
27 *>
28 *> The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
29 *> entries are integers. The right-hand sides are the first NRHS
30 *> columns of M * the identity matrix, and the solutions are the
31 *> first NRHS columns of the inverse Hilbert matrix.
32 *>
33 *> The condition number of the Hilbert matrix grows exponentially with
34 *> its size, roughly as O(e ** (3.5*N)). Additionally, the inverse
35 *> Hilbert matrices beyond a relatively small dimension cannot be
36 *> generated exactly without extra precision. Precision is exhausted
37 *> when the largest entry in the inverse Hilbert matrix is greater than
38 *> 2 to the power of the number of bits in the fraction of the data type
39 *> used plus one, which is 24 for single precision.
40 *>
41 *> In single, the generated solution is exact for N <= 6 and has
42 *> small componentwise error for 7 <= N <= 11.
43 *> \endverbatim
44 *
45 * Arguments:
46 * ==========
47 *
48 *> \param[in] N
49 *> \verbatim
50 *> N is INTEGER
51 *> The dimension of the matrix A.
52 *> \endverbatim
53 *>
54 *> \param[in] NRHS
55 *> \verbatim
56 *> NRHS is INTEGER
57 *> The requested number of right-hand sides.
58 *> \endverbatim
59 *>
60 *> \param[out] A
61 *> \verbatim
62 *> A is REAL array, dimension (LDA, N)
63 *> The generated scaled Hilbert matrix.
64 *> \endverbatim
65 *>
66 *> \param[in] LDA
67 *> \verbatim
68 *> LDA is INTEGER
69 *> The leading dimension of the array A. LDA >= N.
70 *> \endverbatim
71 *>
72 *> \param[out] X
73 *> \verbatim
74 *> X is REAL array, dimension (LDX, NRHS)
75 *> The generated exact solutions. Currently, the first NRHS
76 *> columns of the inverse Hilbert matrix.
77 *> \endverbatim
78 *>
79 *> \param[in] LDX
80 *> \verbatim
81 *> LDX is INTEGER
82 *> The leading dimension of the array X. LDX >= N.
83 *> \endverbatim
84 *>
85 *> \param[out] B
86 *> \verbatim
87 *> B is REAL array, dimension (LDB, NRHS)
88 *> The generated right-hand sides. Currently, the first NRHS
89 *> columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
90 *> \endverbatim
91 *>
92 *> \param[in] LDB
93 *> \verbatim
94 *> LDB is INTEGER
95 *> The leading dimension of the array B. LDB >= N.
96 *> \endverbatim
97 *>
98 *> \param[out] WORK
99 *> \verbatim
100 *> WORK is REAL array, dimension (N)
101 *> \endverbatim
102 *>
103 *> \param[out] INFO
104 *> \verbatim
105 *> INFO is INTEGER
106 *> = 0: successful exit
107 *> = 1: N is too large; the data is still generated but may not
108 *> be not exact.
109 *> < 0: if INFO = -i, the i-th argument had an illegal value
110 *> \endverbatim
111 *
112 * Authors:
113 * ========
114 *
115 *> \author Univ. of Tennessee
116 *> \author Univ. of California Berkeley
117 *> \author Univ. of Colorado Denver
118 *> \author NAG Ltd.
119 *
120 *> \date November 2011
121 *
122 *> \ingroup real_matgen
123 *
124 * =====================================================================
125  SUBROUTINE slahilb( N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
126 *
127 * -- LAPACK auxiliary routine (version 3.4.0) --
128 * -- LAPACK is a software package provided by Univ. of Tennessee, --
129 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
130 * November 2011
131 *
132 * .. Scalar Arguments ..
133  INTEGER n, nrhs, lda, ldx, ldb, info
134 * .. Array Arguments ..
135  REAL a(lda, n), x(ldx, nrhs), b(ldb, nrhs), work(n)
136 * ..
137 *
138 * =====================================================================
139 * .. Local Scalars ..
140  INTEGER tm, ti, r
141  INTEGER m
142  INTEGER i, j
143 
144 * .. Parameters ..
145 * NMAX_EXACT the largest dimension where the generated data is
146 * exact.
147 * NMAX_APPROX the largest dimension where the generated data has
148 * a small componentwise relative error.
149  INTEGER nmax_exact, nmax_approx
150  parameter(nmax_exact = 6, nmax_approx = 11)
151 
152 * ..
153 * .. External Functions
154  EXTERNAL slaset
155  INTRINSIC real
156 * ..
157 * .. Executable Statements ..
158 *
159 * Test the input arguments
160 *
161  info = 0
162  IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
163  info = -1
164  ELSE IF (nrhs .LT. 0) THEN
165  info = -2
166  ELSE IF (lda .LT. n) THEN
167  info = -4
168  ELSE IF (ldx .LT. n) THEN
169  info = -6
170  ELSE IF (ldb .LT. n) THEN
171  info = -8
172  END IF
173  IF (info .LT. 0) THEN
174  CALL xerbla('SLAHILB', -info)
175  RETURN
176  END IF
177  IF (n .GT. nmax_exact) THEN
178  info = 1
179  END IF
180 
181 * Compute M = the LCM of the integers [1, 2*N-1]. The largest
182 * reasonable N is small enough that integers suffice (up to N = 11).
183  m = 1
184  DO i = 2, (2*n-1)
185  tm = m
186  ti = i
187  r = mod(tm, ti)
188  DO WHILE (r .NE. 0)
189  tm = ti
190  ti = r
191  r = mod(tm, ti)
192  END DO
193  m = (m / ti) * i
194  END DO
195 
196 * Generate the scaled Hilbert matrix in A
197  DO j = 1, n
198  DO i = 1, n
199  a(i, j) = REAL(M) / (i + j - 1)
200  END DO
201  END DO
202 
203 * Generate matrix B as simply the first NRHS columns of M * the
204 * identity.
205  CALL slaset('Full', n, nrhs, 0.0, REAL(M), b, ldb)
206 
207 * Generate the true solutions in X. Because B = the first NRHS
208 * columns of M*I, the true solutions are just the first NRHS columns
209 * of the inverse Hilbert matrix.
210  work(1) = n
211  DO j = 2, n
212  work(j) = ( ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1) )
213  $ * (n +j -1)
214  END DO
215 
216  DO j = 1, nrhs
217  DO i = 1, n
218  x(i, j) = (work(i)*work(j)) / (i + j - 1)
219  END DO
220  END DO
221 
222  END
223