LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dorhr_col02()

 subroutine dorhr_col02 ( integer M, integer N, integer MB1, integer NB1, integer NB2, double precision, dimension(6) RESULT )

DORHR_COL02

Purpose:
``` DORHR_COL02 tests DORGTSQR_ROW and DORHR_COL inside DGETSQRHRT
(which calls DLATSQR, DORGTSQR_ROW and DORHR_COL) using DGEMQRT.
Therefore, DLATSQR (part of DGEQR), DGEMQRT (part of DGEMQR)
have to be tested before this test.```
Parameters
 [in] M ``` M is INTEGER Number of rows in test matrix.``` [in] N ``` N is INTEGER Number of columns in test matrix.``` [in] MB1 ``` MB1 is INTEGER Number of row in row block in an input test matrix.``` [in] NB1 ``` NB1 is INTEGER Number of columns in column block an input test matrix.``` [in] NB2 ``` NB2 is INTEGER Number of columns in column block in an output test matrix.``` [out] RESULT ``` RESULT is DOUBLE PRECISION array, dimension (6) Results of each of the six tests below. A is a m-by-n test input matrix to be factored. so that A = Q_gr * ( R ) ( 0 ), Q_qr is an implicit m-by-m orthogonal Q matrix, the result of factorization in blocked WY-representation, stored in ZGEQRT output format. R is a n-by-n upper-triangular matrix, 0 is a (m-n)-by-n zero matrix, Q is an explicit m-by-m orthogonal matrix Q = Q_gr * I C is an m-by-n random matrix, D is an n-by-m random matrix. The six tests are: RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| ) is equivalent to test for | A - Q * R | / (eps * m * |A|), RESULT(2) = |I - (Q**H) * Q| / ( eps * m ), RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|), RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|) RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|) RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|), where: Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are computed using DGEMQRT, Q * C, (Q**H) * C, D * Q, D * (Q**H) are computed using DGEMM.```

Definition at line 119 of file dorhr_col02.f.

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,
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*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
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
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:114
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
double precision function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: dlansy.f:122
Here is the call graph for this function:
Here is the caller graph for this function: