LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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 ZUNHR_COL using ZLATSQR, ZGEMQRT and ZUNGTSQR.
 Therefore, ZLATSQR (part of ZGEQR), ZGEMQRT (part ZGEMQR), ZUNGTSQR
 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.
          ( C is a M-by-N random matrix, D is a N-by-M random matrix )

          RESULT(1) = | A - Q * R | / (eps * m * |A|)
          RESULT(2) = | I - (Q**H) * Q | / (eps * m )
          RESULT(3) = | Q * C - Q * C | / (eps * m * |C|)
          RESULT(4) = | (Q**H) * C - (Q**H) * C | / (eps * m * |C|)
          RESULT(5) = | (D * Q) - D * Q | / (eps * m * |D|)
          RESULT(6) = | D * (Q**H) - D * (Q**H) | / (eps * m * |D|)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2019

Definition at line 89 of file zunhr_col01.f.

89  IMPLICIT NONE
90 *
91 * -- LAPACK test routine (version 3.9.0) --
92 * -- LAPACK is a software package provided by Univ. of Tennessee, --
93 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
94 * November 2019
95 *
96 * .. Scalar Arguments ..
97  INTEGER M, N, MB1, NB1, NB2
98 * .. Return values ..
99  DOUBLE PRECISION RESULT(6)
100 *
101 * =====================================================================
102 *
103 * ..
104 * .. Local allocatable arrays
105  COMPLEX*16, ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
106  $ WORK( : ), T1(:,:), T2(:,:), DIAG(:),
107  $ C(:,:), CF(:,:), D(:,:), DF(:,:)
108  DOUBLE PRECISION, ALLOCATABLE :: RWORK(:)
109 *
110 * .. Parameters ..
111  DOUBLE PRECISION ZERO
112  parameter( zero = 0.0d+0 )
113  COMPLEX*16 CONE, CZERO
114  parameter( cone = ( 1.0d+0, 0.0d+0 ),
115  $ czero = ( 0.0d+0, 0.0d+0 ) )
116 * ..
117 * .. Local Scalars ..
118  LOGICAL TESTZEROS
119  INTEGER INFO, I, J, K, L, LWORK, NB1_UB, NB2_UB, NRB
120  DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
121 * ..
122 * .. Local Arrays ..
123  INTEGER ISEED( 4 )
124  COMPLEX*16 WORKQUERY( 1 )
125 * ..
126 * .. External Functions ..
127  DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY
128  EXTERNAL dlamch, zlange, zlansy
129 * ..
130 * .. External Subroutines ..
131  EXTERNAL zlacpy, zlarnv, zlaset, zlatsqr, zunhr_col,
132  $ zungtsqr, zscal, zgemm, zgemqrt, zherk
133 * ..
134 * .. Intrinsic Functions ..
135  INTRINSIC ceiling, dble, max, min
136 * ..
137 * .. Scalars in Common ..
138  CHARACTER(LEN=32) SRNAMT
139 * ..
140 * .. Common blocks ..
141  COMMON / srmnamc / srnamt
142 * ..
143 * .. Data statements ..
144  DATA iseed / 1988, 1989, 1990, 1991 /
145 *
146 * TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
147 *
148  testzeros = .false.
149 *
150  eps = dlamch( 'Epsilon' )
151  k = min( m, n )
152  l = max( m, n, 1)
153 *
154 * Dynamically allocate local arrays
155 *
156  ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
157  $ c(m,n), cf(m,n),
158  $ d(n,m), df(n,m) )
159 *
160 * Put random numbers into A and copy to AF
161 *
162  DO j = 1, n
163  CALL zlarnv( 2, iseed, m, a( 1, j ) )
164  END DO
165  IF( testzeros ) THEN
166  IF( m.GE.4 ) THEN
167  DO j = 1, n
168  CALL zlarnv( 2, iseed, m/2, a( m/4, j ) )
169  END DO
170  END IF
171  END IF
172  CALL zlacpy( 'Full', m, n, a, m, af, m )
173 *
174 * Number of row blocks in ZLATSQR
175 *
176  nrb = max( 1, ceiling( dble( m - n ) / dble( mb1 - n ) ) )
177 *
178  ALLOCATE ( t1( nb1, n * nrb ) )
179  ALLOCATE ( t2( nb2, n ) )
180  ALLOCATE ( diag( n ) )
181 *
182 * Begin determine LWORK for the array WORK and allocate memory.
183 *
184 * ZLATSQR requires NB1 to be bounded by N.
185 *
186  nb1_ub = min( nb1, n)
187 *
188 * ZGEMQRT requires NB2 to be bounded by N.
189 *
190  nb2_ub = min( nb2, n)
191 *
192  CALL zlatsqr( m, n, mb1, nb1_ub, af, m, t1, nb1,
193  $ workquery, -1, info )
194  lwork = int( workquery( 1 ) )
195  CALL zungtsqr( m, n, mb1, nb1, af, m, t1, nb1, workquery, -1,
196  $ info )
197 
198  lwork = max( lwork, int( workquery( 1 ) ) )
199 *
200 * In ZGEMQRT, WORK is N*NB2_UB if SIDE = 'L',
201 * or M*NB2_UB if SIDE = 'R'.
202 *
203  lwork = max( lwork, nb2_ub * n, nb2_ub * m )
204 *
205  ALLOCATE ( work( lwork ) )
206 *
207 * End allocate memory for WORK.
208 *
209 *
210 * Begin Householder reconstruction routines
211 *
212 * Factor the matrix A in the array AF.
213 *
214  srnamt = 'ZLATSQR'
215  CALL zlatsqr( m, n, mb1, nb1_ub, af, m, t1, nb1, work, lwork,
216  $ info )
217 *
218 * Copy the factor R into the array R.
219 *
220  srnamt = 'ZLACPY'
221  CALL zlacpy( 'U', m, n, af, m, r, m )
222 *
223 * Reconstruct the orthogonal matrix Q.
224 *
225  srnamt = 'ZUNGTSQR'
226  CALL zungtsqr( m, n, mb1, nb1, af, m, t1, nb1, work, lwork,
227  $ info )
228 *
229 * Perform the Householder reconstruction, the result is stored
230 * the arrays AF and T2.
231 *
232  srnamt = 'ZUNHR_COL'
233  CALL zunhr_col( m, n, nb2, af, m, t2, nb2, diag, info )
234 *
235 * Compute the factor R_hr corresponding to the Householder
236 * reconstructed Q_hr and place it in the upper triangle of AF to
237 * match the Q storage format in ZGEQRT. R_hr = R_tsqr * S,
238 * this means changing the sign of I-th row of the matrix R_tsqr
239 * according to sign of of I-th diagonal element DIAG(I) of the
240 * matrix S.
241 *
242  srnamt = 'ZLACPY'
243  CALL zlacpy( 'U', m, n, r, m, af, m )
244 *
245  DO i = 1, n
246  IF( diag( i ).EQ.-cone ) THEN
247  CALL zscal( n+1-i, -cone, af( i, i ), m )
248  END IF
249  END DO
250 *
251 * End Householder reconstruction routines.
252 *
253 *
254 * Generate the m-by-m matrix Q
255 *
256  CALL zlaset( 'Full', m, m, czero, cone, q, m )
257 *
258  srnamt = 'ZGEMQRT'
259  CALL zgemqrt( 'L', 'N', m, m, k, nb2_ub, af, m, t2, nb2, q, m,
260  $ work, info )
261 *
262 * Copy R
263 *
264  CALL zlaset( 'Full', m, n, czero, czero, r, m )
265 *
266  CALL zlacpy( 'Upper', m, n, af, m, r, m )
267 *
268 * TEST 1
269 * Compute |R - (Q**H)*A| / ( eps * m * |A| ) and store in RESULT(1)
270 *
271  CALL zgemm( 'C', 'N', m, n, m, -cone, q, m, a, m, cone, r, m )
272 *
273  anorm = zlange( '1', m, n, a, m, rwork )
274  resid = zlange( '1', m, n, r, m, rwork )
275  IF( anorm.GT.zero ) THEN
276  result( 1 ) = resid / ( eps * max( 1, m ) * anorm )
277  ELSE
278  result( 1 ) = zero
279  END IF
280 *
281 * TEST 2
282 * Compute |I - (Q**H)*Q| / ( eps * m ) and store in RESULT(2)
283 *
284  CALL zlaset( 'Full', m, m, czero, cone, r, m )
285  CALL zherk( 'U', 'C', m, m, -cone, q, m, cone, r, m )
286  resid = zlansy( '1', 'Upper', m, r, m, rwork )
287  result( 2 ) = resid / ( eps * max( 1, m ) )
288 *
289 * Generate random m-by-n matrix C
290 *
291  DO j = 1, n
292  CALL zlarnv( 2, iseed, m, c( 1, j ) )
293  END DO
294  cnorm = zlange( '1', m, n, c, m, rwork )
295  CALL zlacpy( 'Full', m, n, c, m, cf, m )
296 *
297 * Apply Q to C as Q*C = CF
298 *
299  srnamt = 'ZGEMQRT'
300  CALL zgemqrt( 'L', 'N', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
301  $ work, info )
302 *
303 * TEST 3
304 * Compute |CF - Q*C| / ( eps * m * |C| )
305 *
306  CALL zgemm( 'N', 'N', m, n, m, -cone, q, m, c, m, cone, cf, m )
307  resid = zlange( '1', m, n, cf, m, rwork )
308  IF( cnorm.GT.zero ) THEN
309  result( 3 ) = resid / ( eps * max( 1, m ) * cnorm )
310  ELSE
311  result( 3 ) = zero
312  END IF
313 *
314 * Copy C into CF again
315 *
316  CALL zlacpy( 'Full', m, n, c, m, cf, m )
317 *
318 * Apply Q to C as (Q**H)*C = CF
319 *
320  srnamt = 'ZGEMQRT'
321  CALL zgemqrt( 'L', 'C', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
322  $ work, info )
323 *
324 * TEST 4
325 * Compute |CF - (Q**H)*C| / ( eps * m * |C|)
326 *
327  CALL zgemm( 'C', 'N', m, n, m, -cone, q, m, c, m, cone, cf, m )
328  resid = zlange( '1', m, n, cf, m, rwork )
329  IF( cnorm.GT.zero ) THEN
330  result( 4 ) = resid / ( eps * max( 1, m ) * cnorm )
331  ELSE
332  result( 4 ) = zero
333  END IF
334 *
335 * Generate random n-by-m matrix D and a copy DF
336 *
337  DO j = 1, m
338  CALL zlarnv( 2, iseed, n, d( 1, j ) )
339  END DO
340  dnorm = zlange( '1', n, m, d, n, rwork )
341  CALL zlacpy( 'Full', n, m, d, n, df, n )
342 *
343 * Apply Q to D as D*Q = DF
344 *
345  srnamt = 'ZGEMQRT'
346  CALL zgemqrt( 'R', 'N', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
347  $ work, info )
348 *
349 * TEST 5
350 * Compute |DF - D*Q| / ( eps * m * |D| )
351 *
352  CALL zgemm( 'N', 'N', n, m, m, -cone, d, n, q, m, cone, df, n )
353  resid = zlange( '1', n, m, df, n, rwork )
354  IF( dnorm.GT.zero ) THEN
355  result( 5 ) = resid / ( eps * max( 1, m ) * dnorm )
356  ELSE
357  result( 5 ) = zero
358  END IF
359 *
360 * Copy D into DF again
361 *
362  CALL zlacpy( 'Full', n, m, d, n, df, n )
363 *
364 * Apply Q to D as D*QT = DF
365 *
366  srnamt = 'ZGEMQRT'
367  CALL zgemqrt( 'R', 'C', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
368  $ work, info )
369 *
370 * TEST 6
371 * Compute |DF - D*(Q**H)| / ( eps * m * |D| )
372 *
373  CALL zgemm( 'N', 'C', n, m, m, -cone, d, n, q, m, cone, df, n )
374  resid = zlange( '1', n, m, df, n, rwork )
375  IF( dnorm.GT.zero ) THEN
376  result( 6 ) = resid / ( eps * max( 1, m ) * dnorm )
377  ELSE
378  result( 6 ) = zero
379  END IF
380 *
381 * Deallocate all arrays
382 *
383  DEALLOCATE ( a, af, q, r, rwork, work, t1, t2, diag,
384  $ c, d, cf, df )
385 *
386  RETURN
387 *
388 * End of ZUNHR_COL01
389 *
Here is the call graph for this function:
Here is the caller graph for this function:
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:125
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:117
zlatsqr
subroutine zlatsqr(M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK, INFO)
ZLATSQR
Definition: zlatsqr.f:166
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:260
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
zherk
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:175
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:108
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:105
zlarnv
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: zlarnv.f:101
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:70
zgemqrt
subroutine zgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
ZGEMQRT
Definition: zgemqrt.f:170