LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
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 NRHS
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 single_lin
123 *
124 * =====================================================================
125  SUBROUTINE slahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
126 *
127 * -- LAPACK test 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 * .. External Functions
153  EXTERNAL slaset
154  INTRINSIC real
155 * ..
156 * .. Executable Statements ..
157 *
158 * Test the input arguments
159 *
160  info = 0
161  IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
162  info = -1
163  ELSE IF (nrhs .LT. 0) THEN
164  info = -2
165  ELSE IF (lda .LT. n) THEN
166  info = -4
167  ELSE IF (ldx .LT. n) THEN
168  info = -6
169  ELSE IF (ldb .LT. n) THEN
170  info = -8
171  END IF
172  IF (info .LT. 0) THEN
173  CALL xerbla('SLAHILB', -info)
174  return
175  END IF
176  IF (n .GT. nmax_exact) THEN
177  info = 1
178  END IF
179 *
180 * Compute M = the LCM of the integers [1, 2*N-1]. The largest
181 * reasonable N is small enough that integers suffice (up to N = 11).
182  m = 1
183  DO i = 2, (2*n-1)
184  tm = m
185  ti = i
186  r = mod(tm, ti)
187  DO WHILE (r .NE. 0)
188  tm = ti
189  ti = r
190  r = mod(tm, ti)
191  END DO
192  m = (m / ti) * i
193  END DO
194 *
195 * Generate the scaled Hilbert matrix in A
196  DO j = 1, n
197  DO i = 1, n
198  a(i, j) = REAL(M) / (i + j - 1)
199  END DO
200  END DO
201 *
202 * Generate matrix B as simply the first NRHS columns of M * the
203 * identity.
204  CALL slaset('Full', n, nrhs, 0.0, REAL(M), b, ldb)
205 *
206 * Generate the true solutions in X. Because B = the first NRHS
207 * columns of M*I, the true solutions are just the first NRHS columns
208 * of the inverse Hilbert matrix.
209  work(1) = n
210  DO j = 2, n
211  work(j) = ( ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1) )
212  \$ * (n +j -1)
213  END DO
214 *
215  DO j = 1, nrhs
216  DO i = 1, n
217  x(i, j) = (work(i)*work(j)) / (i + j - 1)
218  END DO
219  END DO
220 *
221  END
222