LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
dglmts.f
Go to the documentation of this file.
1 *> \brief \b DGLMTS
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 DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U,
12 * WORK, LWORK, RWORK, RESULT )
13 *
14 * .. Scalar Arguments ..
15 * INTEGER LDA, LDB, LWORK, M, N, P
16 * DOUBLE PRECISION RESULT
17 * ..
18 * .. Array Arguments ..
19 *
20 *
21 *> \par Purpose:
22 * =============
23 *>
24 *> \verbatim
25 *>
26 *> DGLMTS tests DGGGLM - a subroutine for solving the generalized
27 *> linear model problem.
28 *> \endverbatim
29 *
30 * Arguments:
31 * ==========
32 *
33 *> \param[in] N
34 *> \verbatim
35 *> N is INTEGER
36 *> The number of rows of the matrices A and B. N >= 0.
37 *> \endverbatim
38 *>
39 *> \param[in] M
40 *> \verbatim
41 *> M is INTEGER
42 *> The number of columns of the matrix A. M >= 0.
43 *> \endverbatim
44 *>
45 *> \param[in] P
46 *> \verbatim
47 *> P is INTEGER
48 *> The number of columns of the matrix B. P >= 0.
49 *> \endverbatim
50 *>
51 *> \param[in] A
52 *> \verbatim
53 *> A is DOUBLE PRECISION array, dimension (LDA,M)
54 *> The N-by-M matrix A.
55 *> \endverbatim
56 *>
57 *> \param[out] AF
58 *> \verbatim
59 *> AF is DOUBLE PRECISION array, dimension (LDA,M)
60 *> \endverbatim
61 *>
62 *> \param[in] LDA
63 *> \verbatim
64 *> LDA is INTEGER
65 *> The leading dimension of the arrays A, AF. LDA >= max(M,N).
66 *> \endverbatim
67 *>
68 *> \param[in] B
69 *> \verbatim
70 *> B is DOUBLE PRECISION array, dimension (LDB,P)
71 *> The N-by-P matrix A.
72 *> \endverbatim
73 *>
74 *> \param[out] BF
75 *> \verbatim
76 *> BF is DOUBLE PRECISION array, dimension (LDB,P)
77 *> \endverbatim
78 *>
79 *> \param[in] LDB
80 *> \verbatim
81 *> LDB is INTEGER
82 *> The leading dimension of the arrays B, BF. LDB >= max(P,N).
83 *> \endverbatim
84 *>
85 *> \param[in] D
86 *> \verbatim
87 *> D is DOUBLE PRECISION array, dimension( N )
88 *> On input, the left hand side of the GLM.
89 *> \endverbatim
90 *>
91 *> \param[out] DF
92 *> \verbatim
93 *> DF is DOUBLE PRECISION array, dimension( N )
94 *> \endverbatim
95 *>
96 *> \param[out] X
97 *> \verbatim
98 *> X is DOUBLE PRECISION array, dimension( M )
99 *> solution vector X in the GLM problem.
100 *> \endverbatim
101 *>
102 *> \param[out] U
103 *> \verbatim
104 *> U is DOUBLE PRECISION array, dimension( P )
105 *> solution vector U in the GLM problem.
106 *> \endverbatim
107 *>
108 *> \param[out] WORK
109 *> \verbatim
110 *> WORK is DOUBLE PRECISION array, dimension (LWORK)
111 *> \endverbatim
112 *>
113 *> \param[in] LWORK
114 *> \verbatim
115 *> LWORK is INTEGER
116 *> The dimension of the array WORK.
117 *> \endverbatim
118 *>
119 *> \param[out] RWORK
120 *> \verbatim
121 *> RWORK is DOUBLE PRECISION array, dimension (M)
122 *> \endverbatim
123 *>
124 *> \param[out] RESULT
125 *> \verbatim
126 *> RESULT is DOUBLE PRECISION
127 *> The test ratio:
128 *> norm( d - A*x - B*u )
129 *> RESULT = -----------------------------------------
130 *> (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
131 *> \endverbatim
132 *
133 * Authors:
134 * ========
135 *
136 *> \author Univ. of Tennessee
137 *> \author Univ. of California Berkeley
138 *> \author Univ. of Colorado Denver
139 *> \author NAG Ltd.
140 *
141 *> \date November 2011
142 *
143 *> \ingroup double_eig
144 *
145 * =====================================================================
146  SUBROUTINE dglmts( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U,
147  \$ work, lwork, rwork, result )
148 *
149 * -- LAPACK test routine (version 3.4.0) --
150 * -- LAPACK is a software package provided by Univ. of Tennessee, --
151 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
152 * November 2011
153 *
154 * .. Scalar Arguments ..
155  INTEGER lda, ldb, lwork, m, n, p
156  DOUBLE PRECISION result
157 * ..
158 * .. Array Arguments ..
159 *
160 * ====================================================================
161 *
162  DOUBLE PRECISION a( lda, * ), af( lda, * ), b( ldb, * ),
163  \$ bf( ldb, * ), d( * ), df( * ), rwork( * ),
164  \$ u( * ), work( lwork ), x( * )
165 * ..
166 * .. Parameters ..
167  DOUBLE PRECISION zero, one
168  parameter( zero = 0.0d+0, one = 1.0d+0 )
169 * ..
170 * .. Local Scalars ..
171  INTEGER info
172  DOUBLE PRECISION anorm, bnorm, dnorm, eps, unfl, xnorm, ynorm
173 * ..
174 * .. External Functions ..
175  DOUBLE PRECISION dasum, dlamch, dlange
176  EXTERNAL dasum, dlamch, dlange
177 * ..
178 * .. External Subroutines ..
179 *
180  EXTERNAL dcopy, dgemv, dggglm, dlacpy
181 * ..
182 * .. Intrinsic Functions ..
183  INTRINSIC max
184 * ..
185 * .. Executable Statements ..
186 *
187  eps = dlamch( 'Epsilon' )
188  unfl = dlamch( 'Safe minimum' )
189  anorm = max( dlange( '1', n, m, a, lda, rwork ), unfl )
190  bnorm = max( dlange( '1', n, p, b, ldb, rwork ), unfl )
191 *
192 * Copy the matrices A and B to the arrays AF and BF,
193 * and the vector D the array DF.
194 *
195  CALL dlacpy( 'Full', n, m, a, lda, af, lda )
196  CALL dlacpy( 'Full', n, p, b, ldb, bf, ldb )
197  CALL dcopy( n, d, 1, df, 1 )
198 *
199 * Solve GLM problem
200 *
201  CALL dggglm( n, m, p, af, lda, bf, ldb, df, x, u, work, lwork,
202  \$ info )
203 *
204 * Test the residual for the solution of LSE
205 *
206 * norm( d - A*x - B*u )
207 * RESULT = -----------------------------------------
208 * (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
209 *
210  CALL dcopy( n, d, 1, df, 1 )
211  CALL dgemv( 'No transpose', n, m, -one, a, lda, x, 1, one, df, 1 )
212 *
213  CALL dgemv( 'No transpose', n, p, -one, b, ldb, u, 1, one, df, 1 )
214 *
215  dnorm = dasum( n, df, 1 )
216  xnorm = dasum( m, x, 1 ) + dasum( p, u, 1 )
217  ynorm = anorm + bnorm
218 *
219  IF( xnorm.LE.zero ) THEN
220  result = zero
221  ELSE
222  result = ( ( dnorm / ynorm ) / xnorm ) / eps
223  END IF
224 *
225  return
226 *
227 * End of DGLMTS
228 *
229  END