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

◆ cunhr_col02()

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

CUNHR_COL02

Purpose:
 CUNHR_COL02 tests CUNGTSQR_ROW and CUNHR_COL inside CGETSQRHRT
 (which calls CLATSQR, CUNGTSQR_ROW and CUNHR_COL) using CGEMQRT.
 Therefore, CLATSQR (part of CGEQR), CGEMQRT (part of CGEMQR)
 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 unitary Q matrix, the result
            of factorization in blocked WY-representation,
            stored in CGEQRT 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 unitary 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 CGEMQRT,

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

Definition at line 119 of file cunhr_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 REAL RESULT(6)
130*
131* =====================================================================
132*
133* ..
134* .. Local allocatable arrays
135 COMPLEX , ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
136 $ WORK( : ), T1(:,:), T2(:,:), DIAG(:),
137 $ C(:,:), CF(:,:), D(:,:), DF(:,:)
138 REAL , ALLOCATABLE :: RWORK(:)
139*
140* .. Parameters ..
141 REAL ZERO
142 parameter( zero = 0.0e+0 )
143 COMPLEX CONE, CZERO
144 parameter( cone = ( 1.0e+0, 0.0e+0 ),
145 $ czero = ( 0.0e+0, 0.0e+0 ) )
146* ..
147* .. Local Scalars ..
148 LOGICAL TESTZEROS
149 INTEGER INFO, J, K, L, LWORK, NB2_UB, NRB
150 REAL ANORM, EPS, RESID, CNORM, DNORM
151* ..
152* .. Local Arrays ..
153 INTEGER ISEED( 4 )
154 COMPLEX WORKQUERY( 1 )
155* ..
156* .. External Functions ..
157 REAL SLAMCH, CLANGE, CLANSY
158 EXTERNAL slamch, clange, clansy
159* ..
160* .. External Subroutines ..
161 EXTERNAL clacpy, clarnv, claset, cgetsqrhrt,
162 $ cscal, cgemm, cgemqrt, cherk
163* ..
164* .. Intrinsic Functions ..
165 INTRINSIC ceiling, real, max, min
166* ..
167* .. Scalars in Common ..
168 CHARACTER(LEN=32) SRNAMT
169* ..
170* .. Common blocks ..
171 COMMON / srmnamc / srnamt
172* ..
173* .. Data statements ..
174 DATA iseed / 1988, 1989, 1990, 1991 /
175*
176* TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
177*
178 testzeros = .false.
179*
180 eps = slamch( 'Epsilon' )
181 k = min( m, n )
182 l = max( m, n, 1)
183*
184* Dynamically allocate local arrays
185*
186 ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
187 $ c(m,n), cf(m,n),
188 $ d(n,m), df(n,m) )
189*
190* Put random numbers into A and copy to AF
191*
192 DO j = 1, n
193 CALL clarnv( 2, iseed, m, a( 1, j ) )
194 END DO
195 IF( testzeros ) THEN
196 IF( m.GE.4 ) THEN
197 DO j = 1, n
198 CALL clarnv( 2, iseed, m/2, a( m/4, j ) )
199 END DO
200 END IF
201 END IF
202 CALL clacpy( 'Full', m, n, a, m, af, m )
203*
204* Number of row blocks in CLATSQR
205*
206 nrb = max( 1, ceiling( real( m - n ) / real( mb1 - n ) ) )
207*
208 ALLOCATE ( t1( nb1, n * nrb ) )
209 ALLOCATE ( t2( nb2, n ) )
210 ALLOCATE ( diag( n ) )
211*
212* Begin determine LWORK for the array WORK and allocate memory.
213*
214* CGEMQRT requires NB2 to be bounded by N.
215*
216 nb2_ub = min( nb2, n)
217*
218*
219 CALL cgetsqrhrt( m, n, mb1, nb1, nb2, af, m, t2, nb2,
220 $ workquery, -1, info )
221*
222 lwork = int( workquery( 1 ) )
223*
224* In CGEMQRT, WORK is N*NB2_UB if SIDE = 'L',
225* or M*NB2_UB if SIDE = 'R'.
226*
227 lwork = max( lwork, nb2_ub * n, nb2_ub * m )
228*
229 ALLOCATE ( work( lwork ) )
230*
231* End allocate memory for WORK.
232*
233*
234* Begin Householder reconstruction routines
235*
236* Factor the matrix A in the array AF.
237*
238 srnamt = 'CGETSQRHRT'
239 CALL cgetsqrhrt( m, n, mb1, nb1, nb2, af, m, t2, nb2,
240 $ work, lwork, info )
241*
242* End Householder reconstruction routines.
243*
244*
245* Generate the m-by-m matrix Q
246*
247 CALL claset( 'Full', m, m, czero, cone, q, m )
248*
249 srnamt = 'CGEMQRT'
250 CALL cgemqrt( 'L', 'N', m, m, k, nb2_ub, af, m, t2, nb2, q, m,
251 $ work, info )
252*
253* Copy R
254*
255 CALL claset( 'Full', m, n, czero, czero, r, m )
256*
257 CALL clacpy( 'Upper', m, n, af, m, r, m )
258*
259* TEST 1
260* Compute |R - (Q**T)*A| / ( eps * m * |A| ) and store in RESULT(1)
261*
262 CALL cgemm( 'C', 'N', m, n, m, -cone, q, m, a, m, cone, r, m )
263*
264 anorm = clange( '1', m, n, a, m, rwork )
265 resid = clange( '1', m, n, r, m, rwork )
266 IF( anorm.GT.zero ) THEN
267 result( 1 ) = resid / ( eps * max( 1, m ) * anorm )
268 ELSE
269 result( 1 ) = zero
270 END IF
271*
272* TEST 2
273* Compute |I - (Q**T)*Q| / ( eps * m ) and store in RESULT(2)
274*
275 CALL claset( 'Full', m, m, czero, cone, r, m )
276 CALL cherk( 'U', 'C', m, m, real(-cone), q, m, real(cone), r, m )
277 resid = clansy( '1', 'Upper', m, r, m, rwork )
278 result( 2 ) = resid / ( eps * max( 1, m ) )
279*
280* Generate random m-by-n matrix C
281*
282 DO j = 1, n
283 CALL clarnv( 2, iseed, m, c( 1, j ) )
284 END DO
285 cnorm = clange( '1', m, n, c, m, rwork )
286 CALL clacpy( 'Full', m, n, c, m, cf, m )
287*
288* Apply Q to C as Q*C = CF
289*
290 srnamt = 'CGEMQRT'
291 CALL cgemqrt( 'L', 'N', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
292 $ work, info )
293*
294* TEST 3
295* Compute |CF - Q*C| / ( eps * m * |C| )
296*
297 CALL cgemm( 'N', 'N', m, n, m, -cone, q, m, c, m, cone, cf, m )
298 resid = clange( '1', m, n, cf, m, rwork )
299 IF( cnorm.GT.zero ) THEN
300 result( 3 ) = resid / ( eps * max( 1, m ) * cnorm )
301 ELSE
302 result( 3 ) = zero
303 END IF
304*
305* Copy C into CF again
306*
307 CALL clacpy( 'Full', m, n, c, m, cf, m )
308*
309* Apply Q to C as (Q**T)*C = CF
310*
311 srnamt = 'CGEMQRT'
312 CALL cgemqrt( 'L', 'C', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
313 $ work, info )
314*
315* TEST 4
316* Compute |CF - (Q**T)*C| / ( eps * m * |C|)
317*
318 CALL cgemm( 'C', 'N', m, n, m, -cone, q, m, c, m, cone, cf, m )
319 resid = clange( '1', m, n, cf, m, rwork )
320 IF( cnorm.GT.zero ) THEN
321 result( 4 ) = resid / ( eps * max( 1, m ) * cnorm )
322 ELSE
323 result( 4 ) = zero
324 END IF
325*
326* Generate random n-by-m matrix D and a copy DF
327*
328 DO j = 1, m
329 CALL clarnv( 2, iseed, n, d( 1, j ) )
330 END DO
331 dnorm = clange( '1', n, m, d, n, rwork )
332 CALL clacpy( 'Full', n, m, d, n, df, n )
333*
334* Apply Q to D as D*Q = DF
335*
336 srnamt = 'CGEMQRT'
337 CALL cgemqrt( 'R', 'N', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
338 $ work, info )
339*
340* TEST 5
341* Compute |DF - D*Q| / ( eps * m * |D| )
342*
343 CALL cgemm( 'N', 'N', n, m, m, -cone, d, n, q, m, cone, df, n )
344 resid = clange( '1', n, m, df, n, rwork )
345 IF( dnorm.GT.zero ) THEN
346 result( 5 ) = resid / ( eps * max( 1, m ) * dnorm )
347 ELSE
348 result( 5 ) = zero
349 END IF
350*
351* Copy D into DF again
352*
353 CALL clacpy( 'Full', n, m, d, n, df, n )
354*
355* Apply Q to D as D*QT = DF
356*
357 srnamt = 'CGEMQRT'
358 CALL cgemqrt( 'R', 'C', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
359 $ work, info )
360*
361* TEST 6
362* Compute |DF - D*(Q**T)| / ( eps * m * |D| )
363*
364 CALL cgemm( 'N', 'C', n, m, m, -cone, d, n, q, m, cone, df, n )
365 resid = clange( '1', n, m, df, n, rwork )
366 IF( dnorm.GT.zero ) THEN
367 result( 6 ) = resid / ( eps * max( 1, m ) * dnorm )
368 ELSE
369 result( 6 ) = zero
370 END IF
371*
372* Deallocate all arrays
373*
374 DEALLOCATE ( a, af, q, r, rwork, work, t1, t2, diag,
375 $ c, d, cf, df )
376*
377 RETURN
378*
379* End of CUNHR_COL02
380*
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
cgemqrt
subroutine cgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
CGEMQRT
Definition cgemqrt.f:168
cgetsqrhrt
subroutine cgetsqrhrt(m, n, mb1, nb1, nb2, a, lda, t, ldt, work, lwork, info)
CGETSQRHRT
Definition cgetsqrhrt.f:179
cherk
subroutine cherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
CHERK
Definition cherk.f:173
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clange
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:115
clansy
real function clansy(norm, uplo, n, a, lda, work)
CLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clansy.f:123
clarnv
subroutine clarnv(idist, iseed, n, x)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition clarnv.f:99
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
Here is the call graph for this function:
Here is the caller graph for this function: