LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
dorhr_col02.f
Go to the documentation of this file.
1 *> \brief \b DORHR_COL02
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 DORHR_COL02( M, N, MB1, NB1, NB2, RESULT )
12 *
13 * .. Scalar Arguments ..
14 * INTEGER M, N, MB1, NB1, NB2
15 * .. Return values ..
16 * DOUBLE PRECISION RESULT(6)
17 *
18 *
19 *> \par Purpose:
20 * =============
21 *>
22 *> \verbatim
23 *>
24 *> DORHR_COL02 tests DORGTSQR_ROW and DORHR_COL inside DGETSQRHRT
25 *> (which calls DLATSQR, DORGTSQR_ROW and DORHR_COL) using DGEMQRT.
26 *> Therefore, DLATSQR (part of DGEQR), DGEMQRT (part of DGEMQR)
27 *> have to be tested before this test.
28 *>
29 *> \endverbatim
30 *
31 * Arguments:
32 * ==========
33 *
34 *> \param[in] M
35 *> \verbatim
36 *> M is INTEGER
37 *> Number of rows in test matrix.
38 *> \endverbatim
39 *> \param[in] N
40 *> \verbatim
41 *> N is INTEGER
42 *> Number of columns in test matrix.
43 *> \endverbatim
44 *> \param[in] MB1
45 *> \verbatim
46 *> MB1 is INTEGER
47 *> Number of row in row block in an input test matrix.
48 *> \endverbatim
49 *>
50 *> \param[in] NB1
51 *> \verbatim
52 *> NB1 is INTEGER
53 *> Number of columns in column block an input test matrix.
54 *> \endverbatim
55 *>
56 *> \param[in] NB2
57 *> \verbatim
58 *> NB2 is INTEGER
59 *> Number of columns in column block in an output test matrix.
60 *> \endverbatim
61 *>
62 *> \param[out] RESULT
63 *> \verbatim
64 *> RESULT is DOUBLE PRECISION array, dimension (6)
65 *> Results of each of the six tests below.
66 *>
67 *> A is a m-by-n test input matrix to be factored.
68 *> so that A = Q_gr * ( R )
69 *> ( 0 ),
70 *>
71 *> Q_qr is an implicit m-by-m orthogonal Q matrix, the result
72 *> of factorization in blocked WY-representation,
73 *> stored in ZGEQRT output format.
74 *>
75 *> R is a n-by-n upper-triangular matrix,
76 *>
77 *> 0 is a (m-n)-by-n zero matrix,
78 *>
79 *> Q is an explicit m-by-m orthogonal matrix Q = Q_gr * I
80 *>
81 *> C is an m-by-n random matrix,
82 *>
83 *> D is an n-by-m random matrix.
84 *>
85 *> The six tests are:
86 *>
87 *> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| )
88 *> is equivalent to test for | A - Q * R | / (eps * m * |A|),
89 *>
90 *> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ),
91 *>
92 *> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|),
93 *>
94 *> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|)
95 *>
96 *> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|)
97 *>
98 *> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|),
99 *>
100 *> where:
101 *> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are
102 *> computed using DGEMQRT,
103 *>
104 *> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
105 *> computed using DGEMM.
106 *> \endverbatim
107 *
108 * Authors:
109 * ========
110 *
111 *> \author Univ. of Tennessee
112 *> \author Univ. of California Berkeley
113 *> \author Univ. of Colorado Denver
114 *> \author NAG Ltd.
115 *
116 *> \ingroup double_lin
117 *
118 * =====================================================================
119  SUBROUTINE dorhr_col02( M, N, MB1, NB1, NB2, RESULT )
120  IMPLICIT NONE
121 *
122 * -- LAPACK test routine --
123 * -- LAPACK is a software package provided by Univ. of Tennessee, --
124 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
125 *
126 * .. Scalar Arguments ..
127  INTEGER M, N, MB1, NB1, NB2
128 * .. Return values ..
129  DOUBLE PRECISION RESULT(6)
130 *
131 * =====================================================================
132 *
133 * ..
134 * .. Local allocatable arrays
135  DOUBLE PRECISION, ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
136  $ RWORK(:), WORK( : ), T1(:,:), T2(:,:), DIAG(:),
137  $ C(:,:), CF(:,:), D(:,:), DF(:,:)
138 *
139 * .. Parameters ..
140  DOUBLE PRECISION ONE, ZERO
141  parameter( zero = 0.0d+0, one = 1.0d+0 )
142 * ..
143 * .. Local Scalars ..
144  LOGICAL TESTZEROS
145  INTEGER INFO, J, K, L, LWORK, NB2_UB, NRB
146  DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
147 * ..
148 * .. Local Arrays ..
149  INTEGER ISEED( 4 )
150  DOUBLE PRECISION WORKQUERY( 1 )
151 * ..
152 * .. External Functions ..
153  DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
154  EXTERNAL dlamch, dlange, dlansy
155 * ..
156 * .. External Subroutines ..
157  EXTERNAL dlacpy, dlarnv, dlaset, dgetsqrhrt,
158  $ dscal, dgemm, dgemqrt, dsyrk
159 * ..
160 * .. Intrinsic Functions ..
161  INTRINSIC ceiling, dble, max, min
162 * ..
163 * .. Scalars in Common ..
164  CHARACTER(LEN=32) SRNAMT
165 * ..
166 * .. Common blocks ..
167  COMMON / srmnamc / srnamt
168 * ..
169 * .. Data statements ..
170  DATA iseed / 1988, 1989, 1990, 1991 /
171 *
172 * TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
173 *
174  testzeros = .false.
175 *
176  eps = dlamch( 'Epsilon' )
177  k = min( m, n )
178  l = max( m, n, 1)
179 *
180 * Dynamically allocate local arrays
181 *
182  ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
183  $ c(m,n), cf(m,n),
184  $ d(n,m), df(n,m) )
185 *
186 * Put random numbers into A and copy to AF
187 *
188  DO j = 1, n
189  CALL dlarnv( 2, iseed, m, a( 1, j ) )
190  END DO
191  IF( testzeros ) THEN
192  IF( m.GE.4 ) THEN
193  DO j = 1, n
194  CALL dlarnv( 2, iseed, m/2, a( m/4, j ) )
195  END DO
196  END IF
197  END IF
198  CALL dlacpy( 'Full', m, n, a, m, af, m )
199 *
200 * Number of row blocks in DLATSQR
201 *
202  nrb = max( 1, ceiling( dble( m - n ) / dble( mb1 - n ) ) )
203 *
204  ALLOCATE ( t1( nb1, n * nrb ) )
205  ALLOCATE ( t2( nb2, n ) )
206  ALLOCATE ( diag( n ) )
207 *
208 * Begin determine LWORK for the array WORK and allocate memory.
209 *
210 * DGEMQRT requires NB2 to be bounded by N.
211 *
212  nb2_ub = min( nb2, n)
213 *
214 *
215  CALL dgetsqrhrt( m, n, mb1, nb1, nb2, af, m, t2, nb2,
216  $ workquery, -1, info )
217 *
218  lwork = int( workquery( 1 ) )
219 *
220 * In DGEMQRT, WORK is N*NB2_UB if SIDE = 'L',
221 * or M*NB2_UB if SIDE = 'R'.
222 *
223  lwork = max( lwork, nb2_ub * n, nb2_ub * m )
224 *
225  ALLOCATE ( work( lwork ) )
226 *
227 * End allocate memory for WORK.
228 *
229 *
230 * Begin Householder reconstruction routines
231 *
232 * Factor the matrix A in the array AF.
233 *
234  srnamt = 'DGETSQRHRT'
235  CALL dgetsqrhrt( m, n, mb1, nb1, nb2, af, m, t2, nb2,
236  $ work, lwork, info )
237 *
238 * End Householder reconstruction routines.
239 *
240 *
241 * Generate the m-by-m matrix Q
242 *
243  CALL dlaset( 'Full', m, m, zero, one, q, m )
244 *
245  srnamt = 'DGEMQRT'
246  CALL dgemqrt( 'L', 'N', m, m, k, nb2_ub, af, m, t2, nb2, q, m,
247  $ work, info )
248 *
249 * Copy R
250 *
251  CALL dlaset( 'Full', m, n, zero, zero, r, m )
252 *
253  CALL dlacpy( 'Upper', m, n, af, m, r, m )
254 *
255 * TEST 1
256 * Compute |R - (Q**T)*A| / ( eps * m * |A| ) and store in RESULT(1)
257 *
258  CALL dgemm( 'T', 'N', m, n, m, -one, q, m, a, m, one, r, m )
259 *
260  anorm = dlange( '1', m, n, a, m, rwork )
261  resid = dlange( '1', m, n, r, m, rwork )
262  IF( anorm.GT.zero ) THEN
263  result( 1 ) = resid / ( eps * max( 1, m ) * anorm )
264  ELSE
265  result( 1 ) = zero
266  END IF
267 *
268 * TEST 2
269 * Compute |I - (Q**T)*Q| / ( eps * m ) and store in RESULT(2)
270 *
271  CALL dlaset( 'Full', m, m, zero, one, r, m )
272  CALL dsyrk( 'U', 'T', m, m, -one, q, m, one, r, m )
273  resid = dlansy( '1', 'Upper', m, r, m, rwork )
274  result( 2 ) = resid / ( eps * max( 1, m ) )
275 *
276 * Generate random m-by-n matrix C
277 *
278  DO j = 1, n
279  CALL dlarnv( 2, iseed, m, c( 1, j ) )
280  END DO
281  cnorm = dlange( '1', m, n, c, m, rwork )
282  CALL dlacpy( 'Full', m, n, c, m, cf, m )
283 *
284 * Apply Q to C as Q*C = CF
285 *
286  srnamt = 'DGEMQRT'
287  CALL dgemqrt( 'L', 'N', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
288  $ work, info )
289 *
290 * TEST 3
291 * Compute |CF - Q*C| / ( eps * m * |C| )
292 *
293  CALL dgemm( 'N', 'N', m, n, m, -one, q, m, c, m, one, cf, m )
294  resid = dlange( '1', m, n, cf, m, rwork )
295  IF( cnorm.GT.zero ) THEN
296  result( 3 ) = resid / ( eps * max( 1, m ) * cnorm )
297  ELSE
298  result( 3 ) = zero
299  END IF
300 *
301 * Copy C into CF again
302 *
303  CALL dlacpy( 'Full', m, n, c, m, cf, m )
304 *
305 * Apply Q to C as (Q**T)*C = CF
306 *
307  srnamt = 'DGEMQRT'
308  CALL dgemqrt( 'L', 'T', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
309  $ work, info )
310 *
311 * TEST 4
312 * Compute |CF - (Q**T)*C| / ( eps * m * |C|)
313 *
314  CALL dgemm( 'T', 'N', m, n, m, -one, q, m, c, m, one, cf, m )
315  resid = dlange( '1', m, n, cf, m, rwork )
316  IF( cnorm.GT.zero ) THEN
317  result( 4 ) = resid / ( eps * max( 1, m ) * cnorm )
318  ELSE
319  result( 4 ) = zero
320  END IF
321 *
322 * Generate random n-by-m matrix D and a copy DF
323 *
324  DO j = 1, m
325  CALL dlarnv( 2, iseed, n, d( 1, j ) )
326  END DO
327  dnorm = dlange( '1', n, m, d, n, rwork )
328  CALL dlacpy( 'Full', n, m, d, n, df, n )
329 *
330 * Apply Q to D as D*Q = DF
331 *
332  srnamt = 'DGEMQRT'
333  CALL dgemqrt( 'R', 'N', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
334  $ work, info )
335 *
336 * TEST 5
337 * Compute |DF - D*Q| / ( eps * m * |D| )
338 *
339  CALL dgemm( 'N', 'N', n, m, m, -one, d, n, q, m, one, df, n )
340  resid = dlange( '1', n, m, df, n, rwork )
341  IF( dnorm.GT.zero ) THEN
342  result( 5 ) = resid / ( eps * max( 1, m ) * dnorm )
343  ELSE
344  result( 5 ) = zero
345  END IF
346 *
347 * Copy D into DF again
348 *
349  CALL dlacpy( 'Full', n, m, d, n, df, n )
350 *
351 * Apply Q to D as D*QT = DF
352 *
353  srnamt = 'DGEMQRT'
354  CALL dgemqrt( 'R', 'T', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
355  $ work, info )
356 *
357 * TEST 6
358 * Compute |DF - D*(Q**T)| / ( eps * m * |D| )
359 *
360  CALL dgemm( 'N', 'T', n, m, m, -one, d, n, q, m, one, df, n )
361  resid = dlange( '1', n, m, df, n, rwork )
362  IF( dnorm.GT.zero ) THEN
363  result( 6 ) = resid / ( eps * max( 1, m ) * dnorm )
364  ELSE
365  result( 6 ) = zero
366  END IF
367 *
368 * Deallocate all arrays
369 *
370  DEALLOCATE ( a, af, q, r, rwork, work, t1, t2, diag,
371  $ c, d, cf, df )
372 *
373  RETURN
374 *
375 * End of DORHR_COL02
376 *
377  END
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:97
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
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:110
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
Definition: dsyrk.f:169
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dorhr_col02(M, N, MB1, NB1, NB2, RESULT)
DORHR_COL02
Definition: dorhr_col02.f:120
subroutine dgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
DGEMQRT
Definition: dgemqrt.f:168
subroutine dgetsqrhrt(M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK, LWORK, INFO)
DGETSQRHRT
Definition: dgetsqrhrt.f:179