LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ sorhr_col01()

subroutine sorhr_col01 ( integer  m,
integer  n,
integer  mb1,
integer  nb1,
integer  nb2,
real, dimension(6)  result 
)

SORHR_COL01

Purpose:
 SORHR_COL01 tests SORGTSQR and SORHR_COL using SLATSQR, SGEMQRT.
 Therefore, SLATSQR (part of SGEQR), SGEMQRT (part of SGEMQR)
 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 REAL 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 SGEQRT 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 SGEMQRT,

            Q * C, (Q**H) * C, D * Q, D * (Q**H)  are
            computed using SGEMM.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 118 of file sorhr_col01.f.

119 IMPLICIT NONE
120*
121* -- LAPACK test routine --
122* -- LAPACK is a software package provided by Univ. of Tennessee, --
123* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
124*
125* .. Scalar Arguments ..
126 INTEGER M, N, MB1, NB1, NB2
127* .. Return values ..
128 REAL RESULT(6)
129*
130* =====================================================================
131*
132* ..
133* .. Local allocatable arrays
134 REAL , ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
135 $ RWORK(:), WORK( : ), T1(:,:), T2(:,:), DIAG(:),
136 $ C(:,:), CF(:,:), D(:,:), DF(:,:)
137*
138* .. Parameters ..
139 REAL ONE, ZERO
140 parameter( zero = 0.0e+0, one = 1.0e+0 )
141* ..
142* .. Local Scalars ..
143 LOGICAL TESTZEROS
144 INTEGER INFO, I, J, K, L, LWORK, NB1_UB, NB2_UB, NRB
145 REAL ANORM, EPS, RESID, CNORM, DNORM
146* ..
147* .. Local Arrays ..
148 INTEGER ISEED( 4 )
149 REAL WORKQUERY( 1 )
150* ..
151* .. External Functions ..
152 REAL SLAMCH, SLANGE, SLANSY
153 EXTERNAL slamch, slange, slansy
154* ..
155* .. External Subroutines ..
156 EXTERNAL slacpy, slarnv, slaset, slatsqr, sorhr_col,
157 $ sorgtsqr, sscal, sgemm, sgemqrt, ssyrk
158* ..
159* .. Intrinsic Functions ..
160 INTRINSIC ceiling, real, max, min
161* ..
162* .. Scalars in Common ..
163 CHARACTER(LEN=32) SRNAMT
164* ..
165* .. Common blocks ..
166 COMMON / srmnamc / srnamt
167* ..
168* .. Data statements ..
169 DATA iseed / 1988, 1989, 1990, 1991 /
170*
171* TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
172*
173 testzeros = .false.
174*
175 eps = slamch( 'Epsilon' )
176 k = min( m, n )
177 l = max( m, n, 1)
178*
179* Dynamically allocate local arrays
180*
181 ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
182 $ c(m,n), cf(m,n),
183 $ d(n,m), df(n,m) )
184*
185* Put random numbers into A and copy to AF
186*
187 DO j = 1, n
188 CALL slarnv( 2, iseed, m, a( 1, j ) )
189 END DO
190 IF( testzeros ) THEN
191 IF( m.GE.4 ) THEN
192 DO j = 1, n
193 CALL slarnv( 2, iseed, m/2, a( m/4, j ) )
194 END DO
195 END IF
196 END IF
197 CALL slacpy( 'Full', m, n, a, m, af, m )
198*
199* Number of row blocks in SLATSQR
200*
201 nrb = max( 1, ceiling( real( m - n ) / real( mb1 - n ) ) )
202*
203 ALLOCATE ( t1( nb1, n * nrb ) )
204 ALLOCATE ( t2( nb2, n ) )
205 ALLOCATE ( diag( n ) )
206*
207* Begin determine LWORK for the array WORK and allocate memory.
208*
209* SLATSQR requires NB1 to be bounded by N.
210*
211 nb1_ub = min( nb1, n)
212*
213* SGEMQRT requires NB2 to be bounded by N.
214*
215 nb2_ub = min( nb2, n)
216*
217 CALL slatsqr( m, n, mb1, nb1_ub, af, m, t1, nb1,
218 $ workquery, -1, info )
219 lwork = int( workquery( 1 ) )
220 CALL sorgtsqr( m, n, mb1, nb1, af, m, t1, nb1, workquery, -1,
221 $ info )
222
223 lwork = max( lwork, int( workquery( 1 ) ) )
224*
225* In SGEMQRT, WORK is N*NB2_UB if SIDE = 'L',
226* or M*NB2_UB if SIDE = 'R'.
227*
228 lwork = max( lwork, nb2_ub * n, nb2_ub * m )
229*
230 ALLOCATE ( work( lwork ) )
231*
232* End allocate memory for WORK.
233*
234*
235* Begin Householder reconstruction routines
236*
237* Factor the matrix A in the array AF.
238*
239 srnamt = 'SLATSQR'
240 CALL slatsqr( m, n, mb1, nb1_ub, af, m, t1, nb1, work, lwork,
241 $ info )
242*
243* Copy the factor R into the array R.
244*
245 srnamt = 'SLACPY'
246 CALL slacpy( 'U', n, n, af, m, r, m )
247*
248* Reconstruct the orthogonal matrix Q.
249*
250 srnamt = 'SORGTSQR'
251 CALL sorgtsqr( m, n, mb1, nb1, af, m, t1, nb1, work, lwork,
252 $ info )
253*
254* Perform the Householder reconstruction, the result is stored
255* the arrays AF and T2.
256*
257 srnamt = 'SORHR_COL'
258 CALL sorhr_col( m, n, nb2, af, m, t2, nb2, diag, info )
259*
260* Compute the factor R_hr corresponding to the Householder
261* reconstructed Q_hr and place it in the upper triangle of AF to
262* match the Q storage format in SGEQRT. R_hr = R_tsqr * S,
263* this means changing the sign of I-th row of the matrix R_tsqr
264* according to sign of of I-th diagonal element DIAG(I) of the
265* matrix S.
266*
267 srnamt = 'SLACPY'
268 CALL slacpy( 'U', n, n, r, m, af, m )
269*
270 DO i = 1, n
271 IF( diag( i ).EQ.-one ) THEN
272 CALL sscal( n+1-i, -one, af( i, i ), m )
273 END IF
274 END DO
275*
276* End Householder reconstruction routines.
277*
278*
279* Generate the m-by-m matrix Q
280*
281 CALL slaset( 'Full', m, m, zero, one, q, m )
282*
283 srnamt = 'SGEMQRT'
284 CALL sgemqrt( 'L', 'N', m, m, k, nb2_ub, af, m, t2, nb2, q, m,
285 $ work, info )
286*
287* Copy R
288*
289 CALL slaset( 'Full', m, n, zero, zero, r, m )
290*
291 CALL slacpy( 'Upper', m, n, af, m, r, m )
292*
293* TEST 1
294* Compute |R - (Q**T)*A| / ( eps * m * |A| ) and store in RESULT(1)
295*
296 CALL sgemm( 'T', 'N', m, n, m, -one, q, m, a, m, one, r, m )
297*
298 anorm = slange( '1', m, n, a, m, rwork )
299 resid = slange( '1', m, n, r, m, rwork )
300 IF( anorm.GT.zero ) THEN
301 result( 1 ) = resid / ( eps * max( 1, m ) * anorm )
302 ELSE
303 result( 1 ) = zero
304 END IF
305*
306* TEST 2
307* Compute |I - (Q**T)*Q| / ( eps * m ) and store in RESULT(2)
308*
309 CALL slaset( 'Full', m, m, zero, one, r, m )
310 CALL ssyrk( 'U', 'T', m, m, -one, q, m, one, r, m )
311 resid = slansy( '1', 'Upper', m, r, m, rwork )
312 result( 2 ) = resid / ( eps * max( 1, m ) )
313*
314* Generate random m-by-n matrix C
315*
316 DO j = 1, n
317 CALL slarnv( 2, iseed, m, c( 1, j ) )
318 END DO
319 cnorm = slange( '1', m, n, c, m, rwork )
320 CALL slacpy( 'Full', m, n, c, m, cf, m )
321*
322* Apply Q to C as Q*C = CF
323*
324 srnamt = 'SGEMQRT'
325 CALL sgemqrt( 'L', 'N', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
326 $ work, info )
327*
328* TEST 3
329* Compute |CF - Q*C| / ( eps * m * |C| )
330*
331 CALL sgemm( 'N', 'N', m, n, m, -one, q, m, c, m, one, cf, m )
332 resid = slange( '1', m, n, cf, m, rwork )
333 IF( cnorm.GT.zero ) THEN
334 result( 3 ) = resid / ( eps * max( 1, m ) * cnorm )
335 ELSE
336 result( 3 ) = zero
337 END IF
338*
339* Copy C into CF again
340*
341 CALL slacpy( 'Full', m, n, c, m, cf, m )
342*
343* Apply Q to C as (Q**T)*C = CF
344*
345 srnamt = 'SGEMQRT'
346 CALL sgemqrt( 'L', 'T', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
347 $ work, info )
348*
349* TEST 4
350* Compute |CF - (Q**T)*C| / ( eps * m * |C|)
351*
352 CALL sgemm( 'T', 'N', m, n, m, -one, q, m, c, m, one, cf, m )
353 resid = slange( '1', m, n, cf, m, rwork )
354 IF( cnorm.GT.zero ) THEN
355 result( 4 ) = resid / ( eps * max( 1, m ) * cnorm )
356 ELSE
357 result( 4 ) = zero
358 END IF
359*
360* Generate random n-by-m matrix D and a copy DF
361*
362 DO j = 1, m
363 CALL slarnv( 2, iseed, n, d( 1, j ) )
364 END DO
365 dnorm = slange( '1', n, m, d, n, rwork )
366 CALL slacpy( 'Full', n, m, d, n, df, n )
367*
368* Apply Q to D as D*Q = DF
369*
370 srnamt = 'SGEMQRT'
371 CALL sgemqrt( 'R', 'N', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
372 $ work, info )
373*
374* TEST 5
375* Compute |DF - D*Q| / ( eps * m * |D| )
376*
377 CALL sgemm( 'N', 'N', n, m, m, -one, d, n, q, m, one, df, n )
378 resid = slange( '1', n, m, df, n, rwork )
379 IF( dnorm.GT.zero ) THEN
380 result( 5 ) = resid / ( eps * max( 1, m ) * dnorm )
381 ELSE
382 result( 5 ) = zero
383 END IF
384*
385* Copy D into DF again
386*
387 CALL slacpy( 'Full', n, m, d, n, df, n )
388*
389* Apply Q to D as D*QT = DF
390*
391 srnamt = 'SGEMQRT'
392 CALL sgemqrt( 'R', 'T', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
393 $ work, info )
394*
395* TEST 6
396* Compute |DF - D*(Q**T)| / ( eps * m * |D| )
397*
398 CALL sgemm( 'N', 'T', n, m, m, -one, d, n, q, m, one, df, n )
399 resid = slange( '1', n, m, df, n, rwork )
400 IF( dnorm.GT.zero ) THEN
401 result( 6 ) = resid / ( eps * max( 1, m ) * dnorm )
402 ELSE
403 result( 6 ) = zero
404 END IF
405*
406* Deallocate all arrays
407*
408 DEALLOCATE ( a, af, q, r, rwork, work, t1, t2, diag,
409 $ c, d, cf, df )
410*
411 RETURN
412*
413* End of SORHR_COL01
414*
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
sgemqrt
subroutine sgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
SGEMQRT
Definition sgemqrt.f:168
ssyrk
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
slacpy
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
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slange
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:114
slansy
real function slansy(norm, uplo, n, a, lda, work)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansy.f:122
slarnv
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition slarnv.f:97
slaset
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
slatsqr
subroutine slatsqr(m, n, mb, nb, a, lda, t, ldt, work, lwork, info)
SLATSQR
Definition slatsqr.f:169
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
sorgtsqr
subroutine sorgtsqr(m, n, mb, nb, a, lda, t, ldt, work, lwork, info)
SORGTSQR
Definition sorgtsqr.f:176
sorhr_col
subroutine sorhr_col(m, n, nb, a, lda, t, ldt, d, info)
SORHR_COL
Definition sorhr_col.f:259
Here is the call graph for this function:
Here is the caller graph for this function: