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

◆ zunhr_col01()

subroutine zunhr_col01 ( integer  m,
integer  n,
integer  mb1,
integer  nb1,
integer  nb2,
double precision, dimension(6)  result 
)

ZUNHR_COL01

Purpose:
 ZUNHR_COL01 tests ZUNGTSQR and ZUNHR_COL using ZLATSQR, ZGEMQRT.
 Therefore, ZLATSQR (part of ZGEQR), ZGEMQRT (part of ZGEMQR)
 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 unitary 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 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 ZGEMQRT,

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

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