LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sorhr_col02.f
Go to the documentation of this file.
1*> \brief \b SORHR_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 SORHR_COL02( M, N, MB1, NB1, NB2, RESULT )
12*
13* .. Scalar Arguments ..
14* INTEGER M, N, MB1, NB1, NB2
15* .. Return values ..
16* REAL RESULT(6)
17*
18*
19*> \par Purpose:
20* =============
21*>
22*> \verbatim
23*>
24*> SORHR_COL02 tests SORGTSQR_ROW and SORHR_COL inside SGETSQRHRT
25*> (which calls SLATSQR, SORGTSQR_ROW and SORHR_COL) using SGEMQRT.
26*> Therefore, SLATSQR (part of SGEQR), SGEMQRT (part of SGEMQR)
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 REAL 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 SGEQRT 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 SGEMQRT,
103*>
104*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
105*> computed using SGEMM.
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 single_lin
117*
118* =====================================================================
119 SUBROUTINE sorhr_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 REAL RESULT(6)
130*
131* =====================================================================
132*
133* ..
134* .. Local allocatable arrays
135 REAL , ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
136 $ RWORK(:), WORK( : ), T1(:,:), T2(:,:), DIAG(:),
137 $ C(:,:), CF(:,:), D(:,:), DF(:,:)
138*
139* .. Parameters ..
140 REAL ONE, ZERO
141 parameter( zero = 0.0e+0, one = 1.0e+0 )
142* ..
143* .. Local Scalars ..
144 LOGICAL TESTZEROS
145 INTEGER INFO, J, K, L, LWORK, NB2_UB, NRB
146 REAL ANORM, EPS, RESID, CNORM, DNORM
147* ..
148* .. Local Arrays ..
149 INTEGER ISEED( 4 )
150 REAL WORKQUERY( 1 )
151* ..
152* .. External Functions ..
153 REAL SLAMCH, SLANGE, SLANSY
154 EXTERNAL slamch, slange, slansy
155* ..
156* .. External Subroutines ..
157 EXTERNAL slacpy, slarnv, slaset, sgetsqrhrt,
159* ..
160* .. Intrinsic Functions ..
161 INTRINSIC ceiling, real, 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 = slamch( '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 slarnv( 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 slarnv( 2, iseed, m/2, a( m/4, j ) )
195 END DO
196 END IF
197 END IF
198 CALL slacpy( 'Full', m, n, a, m, af, m )
199*
200* Number of row blocks in SLATSQR
201*
202 nrb = max( 1, ceiling( real( m - n ) / real( 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* SGEMQRT requires NB2 to be bounded by N.
211*
212 nb2_ub = min( nb2, n)
213*
214 CALL sgetsqrhrt( m, n, mb1, nb1, nb2, af, m, t2, nb2,
215 $ workquery, -1, info )
216*
217 lwork = int( workquery( 1 ) )
218*
219* In SGEMQRT, WORK is N*NB2_UB if SIDE = 'L',
220* or M*NB2_UB if SIDE = 'R'.
221*
222 lwork = max( lwork, nb2_ub * n, nb2_ub * m )
223*
224 ALLOCATE ( work( lwork ) )
225*
226* End allocate memory for WORK.
227*
228*
229* Begin Householder reconstruction routines
230*
231* Factor the matrix A in the array AF.
232*
233 srnamt = 'SGETSQRHRT'
234 CALL sgetsqrhrt( m, n, mb1, nb1, nb2, af, m, t2, nb2,
235 $ work, lwork, info )
236*
237* End Householder reconstruction routines.
238*
239*
240* Generate the m-by-m matrix Q
241*
242 CALL slaset( 'Full', m, m, zero, one, q, m )
243*
244 srnamt = 'SGEMQRT'
245 CALL sgemqrt( 'L', 'N', m, m, k, nb2_ub, af, m, t2, nb2, q, m,
246 $ work, info )
247*
248* Copy R
249*
250 CALL slaset( 'Full', m, n, zero, zero, r, m )
251*
252 CALL slacpy( 'Upper', m, n, af, m, r, m )
253*
254* TEST 1
255* Compute |R - (Q**T)*A| / ( eps * m * |A| ) and store in RESULT(1)
256*
257 CALL sgemm( 'T', 'N', m, n, m, -one, q, m, a, m, one, r, m )
258*
259 anorm = slange( '1', m, n, a, m, rwork )
260 resid = slange( '1', m, n, r, m, rwork )
261 IF( anorm.GT.zero ) THEN
262 result( 1 ) = resid / ( eps * max( 1, m ) * anorm )
263 ELSE
264 result( 1 ) = zero
265 END IF
266*
267* TEST 2
268* Compute |I - (Q**T)*Q| / ( eps * m ) and store in RESULT(2)
269*
270 CALL slaset( 'Full', m, m, zero, one, r, m )
271 CALL ssyrk( 'U', 'T', m, m, -one, q, m, one, r, m )
272 resid = slansy( '1', 'Upper', m, r, m, rwork )
273 result( 2 ) = resid / ( eps * max( 1, m ) )
274*
275* Generate random m-by-n matrix C
276*
277 DO j = 1, n
278 CALL slarnv( 2, iseed, m, c( 1, j ) )
279 END DO
280 cnorm = slange( '1', m, n, c, m, rwork )
281 CALL slacpy( 'Full', m, n, c, m, cf, m )
282*
283* Apply Q to C as Q*C = CF
284*
285 srnamt = 'SGEMQRT'
286 CALL sgemqrt( 'L', 'N', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
287 $ work, info )
288*
289* TEST 3
290* Compute |CF - Q*C| / ( eps * m * |C| )
291*
292 CALL sgemm( 'N', 'N', m, n, m, -one, q, m, c, m, one, cf, m )
293 resid = slange( '1', m, n, cf, m, rwork )
294 IF( cnorm.GT.zero ) THEN
295 result( 3 ) = resid / ( eps * max( 1, m ) * cnorm )
296 ELSE
297 result( 3 ) = zero
298 END IF
299*
300* Copy C into CF again
301*
302 CALL slacpy( 'Full', m, n, c, m, cf, m )
303*
304* Apply Q to C as (Q**T)*C = CF
305*
306 srnamt = 'SGEMQRT'
307 CALL sgemqrt( 'L', 'T', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
308 $ work, info )
309*
310* TEST 4
311* Compute |CF - (Q**T)*C| / ( eps * m * |C|)
312*
313 CALL sgemm( 'T', 'N', m, n, m, -one, q, m, c, m, one, cf, m )
314 resid = slange( '1', m, n, cf, m, rwork )
315 IF( cnorm.GT.zero ) THEN
316 result( 4 ) = resid / ( eps * max( 1, m ) * cnorm )
317 ELSE
318 result( 4 ) = zero
319 END IF
320*
321* Generate random n-by-m matrix D and a copy DF
322*
323 DO j = 1, m
324 CALL slarnv( 2, iseed, n, d( 1, j ) )
325 END DO
326 dnorm = slange( '1', n, m, d, n, rwork )
327 CALL slacpy( 'Full', n, m, d, n, df, n )
328*
329* Apply Q to D as D*Q = DF
330*
331 srnamt = 'SGEMQRT'
332 CALL sgemqrt( 'R', 'N', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
333 $ work, info )
334*
335* TEST 5
336* Compute |DF - D*Q| / ( eps * m * |D| )
337*
338 CALL sgemm( 'N', 'N', n, m, m, -one, d, n, q, m, one, df, n )
339 resid = slange( '1', n, m, df, n, rwork )
340 IF( dnorm.GT.zero ) THEN
341 result( 5 ) = resid / ( eps * max( 1, m ) * dnorm )
342 ELSE
343 result( 5 ) = zero
344 END IF
345*
346* Copy D into DF again
347*
348 CALL slacpy( 'Full', n, m, d, n, df, n )
349*
350* Apply Q to D as D*QT = DF
351*
352 srnamt = 'SGEMQRT'
353 CALL sgemqrt( 'R', 'T', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
354 $ work, info )
355*
356* TEST 6
357* Compute |DF - D*(Q**T)| / ( eps * m * |D| )
358*
359 CALL sgemm( 'N', 'T', n, m, m, -one, d, n, q, m, one, df, n )
360 resid = slange( '1', n, m, df, n, rwork )
361 IF( dnorm.GT.zero ) THEN
362 result( 6 ) = resid / ( eps * max( 1, m ) * dnorm )
363 ELSE
364 result( 6 ) = zero
365 END IF
366*
367* Deallocate all arrays
368*
369 DEALLOCATE ( a, af, q, r, rwork, work, t1, t2, diag,
370 $ c, d, cf, df )
371*
372 RETURN
373*
374* End of SORHR_COL02
375*
376 END
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine sgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
SGEMQRT
Definition sgemqrt.f:168
subroutine sgetsqrhrt(m, n, mb1, nb1, nb2, a, lda, t, ldt, work, lwork, info)
SGETSQRHRT
Definition sgetsqrhrt.f:179
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition slarnv.f:97
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sorhr_col02(m, n, mb1, nb1, nb2, result)
SORHR_COL02