LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
dlahilb.f
Go to the documentation of this file.
1 *> \brief \b DLAHILB
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 DLAHILB(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 * DOUBLE PRECISION A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
17 * ..
18 *
19 *
20 *> \par Purpose:
21 * =============
22 *>
23 *> \verbatim
24 *>
25 *> DLAHILB 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 double_lin
123 *
124 * =====================================================================
125  SUBROUTINE dlahilb(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  DOUBLE PRECISION 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  COMPLEX*16 TMP
144 * ..
145 * .. Parameters ..
146 * NMAX_EXACT the largest dimension where the generated data is
147 * exact.
148 * NMAX_APPROX the largest dimension where the generated data has
149 * a small componentwise relative error.
150  INTEGER NMAX_EXACT, NMAX_APPROX
151  parameter(nmax_exact = 6, nmax_approx = 11)
152
153 * ..
154 * .. External Functions
155  EXTERNAL dlaset
156  INTRINSIC dble
157 * ..
158 * .. Executable Statements ..
159 *
160 * Test the input arguments
161 *
162  info = 0
163  IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
164  info = -1
165  ELSE IF (nrhs .LT. 0) THEN
166  info = -2
167  ELSE IF (lda .LT. n) THEN
168  info = -4
169  ELSE IF (ldx .LT. n) THEN
170  info = -6
171  ELSE IF (ldb .LT. n) THEN
172  info = -8
173  END IF
174  IF (info .LT. 0) THEN
175  CALL xerbla('DLAHILB', -info)
176  RETURN
177  END IF
178  IF (n .GT. nmax_exact) THEN
179  info = 1
180  END IF
181 *
182 * Compute M = the LCM of the integers [1, 2*N-1]. The largest
183 * reasonable N is small enough that integers suffice (up to N = 11).
184  m = 1
185  DO i = 2, (2*n-1)
186  tm = m
187  ti = i
188  r = mod(tm, ti)
189  DO WHILE (r .NE. 0)
190  tm = ti
191  ti = r
192  r = mod(tm, ti)
193  END DO
194  m = (m / ti) * i
195  END DO
196 *
197 * Generate the scaled Hilbert matrix in A
198  DO j = 1, n
199  DO i = 1, n
200  a(i, j) = dble(m) / (i + j - 1)
201  END DO
202  END DO
203 *
204 * Generate matrix B as simply the first NRHS columns of M * the
205 * identity.
206  tmp = dble(m)
207  CALL dlaset('Full', n, nrhs, 0.0d+0, tmp, b, ldb)
208 *
209 * Generate the true solutions in X. Because B = the first NRHS
210 * columns of M*I, the true solutions are just the first NRHS columns
211 * of the inverse Hilbert matrix.
212  work(1) = n
213  DO j = 2, n
214  work(j) = ( ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1) )
215  \$ * (n +j -1)
216  END DO
217 *
218  DO j = 1, nrhs
219  DO i = 1, n
220  x(i, j) = (work(i)*work(j)) / (i + j - 1)
221  END DO
222  END DO
223 *
224  END
225
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
DLAHILB
Definition: dlahilb.f:126