OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ztsqr01.f
Go to the documentation of this file.
1*> \brief \b ZTSQR01
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 ZTSQR01(TSSW, M,N, MB, NB, RESULT)
12*
13* .. Scalar Arguments ..
14* INTEGER M, N, MB
15* .. Return values ..
16* DOUBLE PRECISION RESULT(6)
17*
18*
19*> \par Purpose:
20* =============
21*>
22*> \verbatim
23*>
24*> ZTSQR01 tests ZGEQR , ZGELQ, ZGEMLQ and ZGEMQR.
25*> \endverbatim
26*
27* Arguments:
28* ==========
29*
30*> \param[in] TSSW
31*> \verbatim
32*> TSSW is CHARACTER
33*> 'TS' for testing tall skinny QR
34*> and anything else for testing short wide LQ
35*> \endverbatim
36*> \param[in] M
37*> \verbatim
38*> M is INTEGER
39*> Number of rows in test matrix.
40*> \endverbatim
41*>
42*> \param[in] N
43*> \verbatim
44*> N is INTEGER
45*> Number of columns in test matrix.
46*> \endverbatim
47*> \param[in] MB
48*> \verbatim
49*> MB is INTEGER
50*> Number of row in row block in test matrix.
51*> \endverbatim
52*>
53*> \param[in] NB
54*> \verbatim
55*> NB is INTEGER
56*> Number of columns in column block test matrix.
57*> \endverbatim
58*>
59*> \param[out] RESULT
60*> \verbatim
61*> RESULT is DOUBLE PRECISION array, dimension (6)
62*> Results of each of the six tests below.
63*>
64*> RESULT(1) = | A - Q R | or | A - L Q |
65*> RESULT(2) = | I - Q^H Q | or | I - Q Q^H |
66*> RESULT(3) = | Q C - Q C |
67*> RESULT(4) = | Q^H C - Q^H C |
68*> RESULT(5) = | C Q - C Q |
69*> RESULT(6) = | C Q^H - C Q^H |
70*> \endverbatim
71*
72* Authors:
73* ========
74*
75*> \author Univ. of Tennessee
76*> \author Univ. of California Berkeley
77*> \author Univ. of Colorado Denver
78*> \author NAG Ltd.
79*
80* =====================================================================
81 SUBROUTINE ztsqr01(TSSW, M, N, MB, NB, RESULT)
82 IMPLICIT NONE
83*
84* -- LAPACK test routine --
85* -- LAPACK is a software package provided by Univ. of Tennessee, --
86* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
87*
88* .. Scalar Arguments ..
89 CHARACTER TSSW
90 INTEGER M, N, MB, NB
91* .. Return values ..
92 DOUBLE PRECISION RESULT(6)
93*
94* =====================================================================
95*
96* ..
97* .. Local allocatable arrays
98 COMPLEX*16, ALLOCATABLE :: AF(:,:), Q(:,:),
99 $ R(:,:), RWORK(:), WORK( : ), T(:),
100 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:), LQ(:,:)
101*
102* .. Parameters ..
103 DOUBLE PRECISION ZERO
104 COMPLEX*16 ONE, CZERO
105 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
106* ..
107* .. Local Scalars ..
108 LOGICAL TESTZEROS, TS
109 INTEGER INFO, J, K, L, LWORK, TSIZE, MNB
110 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
111* ..
112* .. Local Arrays ..
113 INTEGER ISEED( 4 )
114 COMPLEX*16 TQUERY( 5 ), WORKQUERY( 1 )
115* ..
116* .. External Functions ..
117 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY
118 LOGICAL LSAME
119 INTEGER ILAENV
120 EXTERNAL dlamch, zlange, zlansy, lsame, ilaenv
121* ..
122* .. Intrinsic Functions ..
123 INTRINSIC max, min
124* .. Scalars in Common ..
125 CHARACTER*32 srnamt
126* ..
127* .. Common blocks ..
128 COMMON / srnamc / srnamt
129* ..
130* .. Data statements ..
131 DATA iseed / 1988, 1989, 1990, 1991 /
132*
133* TEST TALL SKINNY OR SHORT WIDE
134*
135 ts = lsame(tssw, 'TS')
136*
137* TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
138*
139 testzeros = .false.
140*
141 eps = dlamch( 'Epsilon' )
142 k = min(m,n)
143 l = max(m,n,1)
144 mnb = max( mb, nb)
145 lwork = max(3,l)*mnb
146*
147* Dynamically allocate local arrays
148*
149 ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
150 $ c(m,n), cf(m,n),
151 $ d(n,m), df(n,m), lq(l,n) )
152*
153* Put random numbers into A and copy to AF
154*
155 DO j=1,n
156 CALL zlarnv( 2, iseed, m, a( 1, j ) )
157 END DO
158 IF (testzeros) THEN
159 IF (m.GE.4) THEN
160 DO j=1,n
161 CALL zlarnv( 2, iseed, m/2, a( m/4, j ) )
162 END DO
163 END IF
164 END IF
165 CALL zlacpy( 'Full', m, n, a, m, af, m )
166*
167 IF (ts) THEN
168*
169* Factor the matrix A in the array AF.
170*
171 CALL zgeqr( m, n, af, m, tquery, -1, workquery, -1, info )
172 tsize = int( tquery( 1 ) )
173 lwork = int( workquery( 1 ) )
174 CALL zgemqr( 'L', 'N', m, m, k, af, m, tquery, tsize, cf, m,
175 $ workquery, -1, info)
176 lwork = max( lwork, int( workquery( 1 ) ) )
177 CALL zgemqr( 'L', 'N', m, n, k, af, m, tquery, tsize, cf, m,
178 $ workquery, -1, info)
179 lwork = max( lwork, int( workquery( 1 ) ) )
180 CALL zgemqr( 'L', 'C', m, n, k, af, m, tquery, tsize, cf, m,
181 $ workquery, -1, info)
182 lwork = max( lwork, int( workquery( 1 ) ) )
183 CALL zgemqr( 'R', 'N', n, m, k, af, m, tquery, tsize, df, n,
184 $ workquery, -1, info)
185 lwork = max( lwork, int( workquery( 1 ) ) )
186 CALL zgemqr( 'R', 'C', n, m, k, af, m, tquery, tsize, df, n,
187 $ workquery, -1, info)
188 lwork = max( lwork, int( workquery( 1 ) ) )
189 ALLOCATE ( t( tsize ) )
190 ALLOCATE ( work( lwork ) )
191 srnamt = 'ZGEQR'
192 CALL zgeqr( m, n, af, m, t, tsize, work, lwork, info )
193*
194* Generate the m-by-m matrix Q
195*
196 CALL zlaset( 'Full', m, m, czero, one, q, m )
197 srnamt = 'zgemqr'
198 CALL ZGEMQR( 'l', 'n', M, M, K, AF, M, T, TSIZE, Q, M,
199 $ WORK, LWORK, INFO )
200*
201* Copy R
202*
203 CALL ZLASET( 'full', M, N, CZERO, CZERO, R, M )
204 CALL ZLACPY( 'upper', M, N, AF, M, R, M )
205*
206* Compute |R - Q'*A| / |A| and store in RESULT(1)
207*
208 CALL ZGEMM( 'c', 'n', M, N, M, -ONE, Q, M, A, M, ONE, R, M )
209 ANORM = ZLANGE( '1', M, N, A, M, RWORK )
210 RESID = ZLANGE( '1', M, N, R, M, RWORK )
211.GT. IF( ANORMZERO ) THEN
212 RESULT( 1 ) = RESID / (EPS*MAX(1,M)*ANORM)
213 ELSE
214 RESULT( 1 ) = ZERO
215 END IF
216*
217* Compute |I - Q'*Q| and store in RESULT(2)
218*
219 CALL ZLASET( 'full', M, M, CZERO, ONE, R, M )
220 CALL ZHERK( 'u', 'c', M, M, DREAL(-ONE), Q, M, DREAL(ONE), R, M )
221 RESID = ZLANSY( '1', 'upper', M, R, M, RWORK )
222 RESULT( 2 ) = RESID / (EPS*MAX(1,M))
223*
224* Generate random m-by-n matrix C and a copy CF
225*
226 DO J=1,N
227 CALL ZLARNV( 2, ISEED, M, C( 1, J ) )
228 END DO
229 CNORM = ZLANGE( '1', M, N, C, M, RWORK)
230 CALL ZLACPY( 'full', M, N, C, M, CF, M )
231*
232* Apply Q to C as Q*C
233*
234 srnamt = 'zgemqr'
235 CALL ZGEMQR( 'l', 'n', M, N, K, AF, M, T, TSIZE, CF, M,
236 $ WORK, LWORK, INFO)
237*
238* Compute |Q*C - Q*C| / |C|
239*
240 CALL ZGEMM( 'n', 'n', M, N, M, -ONE, Q, M, C, M, ONE, CF, M )
241 RESID = ZLANGE( '1', M, N, CF, M, RWORK )
242.GT. IF( CNORMZERO ) THEN
243 RESULT( 3 ) = RESID / (EPS*MAX(1,M)*CNORM)
244 ELSE
245 RESULT( 3 ) = ZERO
246 END IF
247*
248* Copy C into CF again
249*
250 CALL ZLACPY( 'full', M, N, C, M, CF, M )
251*
252* Apply Q to C as QT*C
253*
254 srnamt = 'zgemqr'
255 CALL ZGEMQR( 'l', 'c', M, N, K, AF, M, T, TSIZE, CF, M,
256 $ WORK, LWORK, INFO)
257*
258* Compute |QT*C - QT*C| / |C|
259*
260 CALL ZGEMM( 'c', 'n', M, N, M, -ONE, Q, M, C, M, ONE, CF, M )
261 RESID = ZLANGE( '1', M, N, CF, M, RWORK )
262.GT. IF( CNORMZERO ) THEN
263 RESULT( 4 ) = RESID / (EPS*MAX(1,M)*CNORM)
264 ELSE
265 RESULT( 4 ) = ZERO
266 END IF
267*
268* Generate random n-by-m matrix D and a copy DF
269*
270 DO J=1,M
271 CALL ZLARNV( 2, ISEED, N, D( 1, J ) )
272 END DO
273 DNORM = ZLANGE( '1', N, M, D, N, RWORK)
274 CALL ZLACPY( 'full', N, M, D, N, DF, N )
275*
276* Apply Q to D as D*Q
277*
278 srnamt = 'zgemqr'
279 CALL ZGEMQR( 'r', 'n', N, M, K, AF, M, T, TSIZE, DF, N,
280 $ WORK, LWORK, INFO)
281*
282* Compute |D*Q - D*Q| / |D|
283*
284 CALL ZGEMM( 'n', 'n', N, M, M, -ONE, D, N, Q, M, ONE, DF, N )
285 RESID = ZLANGE( '1', N, M, DF, N, RWORK )
286.GT. IF( DNORMZERO ) THEN
287 RESULT( 5 ) = RESID / (EPS*MAX(1,M)*DNORM)
288 ELSE
289 RESULT( 5 ) = ZERO
290 END IF
291*
292* Copy D into DF again
293*
294 CALL ZLACPY( 'full', N, M, D, N, DF, N )
295*
296* Apply Q to D as D*QT
297*
298 CALL ZGEMQR( 'r', 'c', N, M, K, AF, M, T, TSIZE, DF, N,
299 $ WORK, LWORK, INFO)
300*
301* Compute |D*QT - D*QT| / |D|
302*
303 CALL ZGEMM( 'n', 'c', N, M, M, -ONE, D, N, Q, M, ONE, DF, N )
304 RESID = ZLANGE( '1', N, M, DF, N, RWORK )
305.GT. IF( CNORMZERO ) THEN
306 RESULT( 6 ) = RESID / (EPS*MAX(1,M)*DNORM)
307 ELSE
308 RESULT( 6 ) = ZERO
309 END IF
310*
311* Short and wide
312*
313 ELSE
314 CALL ZGELQ( M, N, AF, M, TQUERY, -1, WORKQUERY, -1, INFO )
315 TSIZE = INT( TQUERY( 1 ) )
316 LWORK = INT( WORKQUERY( 1 ) )
317 CALL ZGEMLQ( 'r', 'n', N, N, K, AF, M, TQUERY, TSIZE, Q, N,
318 $ WORKQUERY, -1, INFO )
319 LWORK = MAX( LWORK, INT( WORKQUERY( 1 ) ) )
320 CALL ZGEMLQ( 'l', 'n', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
321 $ WORKQUERY, -1, INFO)
322 LWORK = MAX( LWORK, INT( WORKQUERY( 1 ) ) )
323 CALL ZGEMLQ( 'l', 'c', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
324 $ WORKQUERY, -1, INFO)
325 LWORK = MAX( LWORK, INT( WORKQUERY( 1 ) ) )
326 CALL ZGEMLQ( 'r', 'n', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
327 $ WORKQUERY, -1, INFO)
328 LWORK = MAX( LWORK, INT( WORKQUERY( 1 ) ) )
329 CALL ZGEMLQ( 'r', 'c', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
330 $ WORKQUERY, -1, INFO)
331 LWORK = MAX( LWORK, INT( WORKQUERY( 1 ) ) )
332 ALLOCATE ( T( TSIZE ) )
333 ALLOCATE ( WORK( LWORK ) )
334 srnamt = 'zgelq'
335 CALL ZGELQ( M, N, AF, M, T, TSIZE, WORK, LWORK, INFO )
336*
337*
338* Generate the n-by-n matrix Q
339*
340 CALL ZLASET( 'full', N, N, CZERO, ONE, Q, N )
341 srnamt = 'zgemlq'
342 CALL ZGEMLQ( 'r', 'n', N, N, K, AF, M, T, TSIZE, Q, N,
343 $ WORK, LWORK, INFO )
344*
345* Copy R
346*
347 CALL ZLASET( 'full', M, N, CZERO, CZERO, LQ, L )
348 CALL ZLACPY( 'lower', M, N, AF, M, LQ, L )
349*
350* Compute |L - A*Q'| / |A| and store in RESULT(1)
351*
352 CALL ZGEMM( 'n', 'c', M, N, N, -ONE, A, M, Q, N, ONE, LQ, L )
353 ANORM = ZLANGE( '1', M, N, A, M, RWORK )
354 RESID = ZLANGE( '1', M, N, LQ, L, RWORK )
355.GT. IF( ANORMZERO ) THEN
356 RESULT( 1 ) = RESID / (EPS*MAX(1,N)*ANORM)
357 ELSE
358 RESULT( 1 ) = ZERO
359 END IF
360*
361* Compute |I - Q'*Q| and store in RESULT(2)
362*
363 CALL ZLASET( 'full', N, N, CZERO, ONE, LQ, L )
364 CALL ZHERK( 'u', 'c', N, N, DREAL(-ONE), Q, N, DREAL(ONE), LQ, L)
365 RESID = ZLANSY( '1', 'upper', N, LQ, L, RWORK )
366 RESULT( 2 ) = RESID / (EPS*MAX(1,N))
367*
368* Generate random m-by-n matrix C and a copy CF
369*
370 DO J=1,M
371 CALL ZLARNV( 2, ISEED, N, D( 1, J ) )
372 END DO
373 DNORM = ZLANGE( '1', N, M, D, N, RWORK)
374 CALL ZLACPY( 'full', N, M, D, N, DF, N )
375*
376* Apply Q to C as Q*C
377*
378 CALL ZGEMLQ( 'l', 'n', N, M, K, AF, M, T, TSIZE, DF, N,
379 $ WORK, LWORK, INFO)
380*
381* Compute |Q*D - Q*D| / |D|
382*
383 CALL ZGEMM( 'n', 'n', N, M, N, -ONE, Q, N, D, N, ONE, DF, N )
384 RESID = ZLANGE( '1', N, M, DF, N, RWORK )
385.GT. IF( DNORMZERO ) THEN
386 RESULT( 3 ) = RESID / (EPS*MAX(1,N)*DNORM)
387 ELSE
388 RESULT( 3 ) = ZERO
389 END IF
390*
391* Copy D into DF again
392*
393 CALL ZLACPY( 'full', N, M, D, N, DF, N )
394*
395* Apply Q to D as QT*D
396*
397 CALL ZGEMLQ( 'l', 'c', N, M, K, AF, M, T, TSIZE, DF, N,
398 $ WORK, LWORK, INFO)
399*
400* Compute |QT*D - QT*D| / |D|
401*
402 CALL ZGEMM( 'c', 'n', N, M, N, -ONE, Q, N, D, N, ONE, DF, N )
403 RESID = ZLANGE( '1', N, M, DF, N, RWORK )
404.GT. IF( DNORMZERO ) THEN
405 RESULT( 4 ) = RESID / (EPS*MAX(1,N)*DNORM)
406 ELSE
407 RESULT( 4 ) = ZERO
408 END IF
409*
410* Generate random n-by-m matrix D and a copy DF
411*
412 DO J=1,N
413 CALL ZLARNV( 2, ISEED, M, C( 1, J ) )
414 END DO
415 CNORM = ZLANGE( '1', M, N, C, M, RWORK)
416 CALL ZLACPY( 'full', M, N, C, M, CF, M )
417*
418* Apply Q to C as C*Q
419*
420 CALL ZGEMLQ( 'r', 'n', M, N, K, AF, M, T, TSIZE, CF, M,
421 $ WORK, LWORK, INFO)
422*
423* Compute |C*Q - C*Q| / |C|
424*
425 CALL ZGEMM( 'n', 'n', M, N, N, -ONE, C, M, Q, N, ONE, CF, M )
426 RESID = ZLANGE( '1', N, M, DF, N, RWORK )
427.GT. IF( CNORMZERO ) THEN
428 RESULT( 5 ) = RESID / (EPS*MAX(1,N)*CNORM)
429 ELSE
430 RESULT( 5 ) = ZERO
431 END IF
432*
433* Copy C into CF again
434*
435 CALL ZLACPY( 'full', M, N, C, M, CF, M )
436*
437* Apply Q to D as D*QT
438*
439 CALL ZGEMLQ( 'r', 'c', M, N, K, AF, M, T, TSIZE, CF, M,
440 $ WORK, LWORK, INFO)
441*
442* Compute |C*QT - C*QT| / |C|
443*
444 CALL ZGEMM( 'n', 'c', M, N, N, -ONE, C, M, Q, N, ONE, CF, M )
445 RESID = ZLANGE( '1', M, N, CF, M, RWORK )
446.GT. IF( CNORMZERO ) THEN
447 RESULT( 6 ) = RESID / (EPS*MAX(1,N)*CNORM)
448 ELSE
449 RESULT( 6 ) = ZERO
450 END IF
451*
452 END IF
453*
454* Deallocate all arrays
455*
456 DEALLOCATE ( A, AF, Q, R, RWORK, WORK, T, C, D, CF, DF)
457*
458 RETURN
459 END
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
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
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
zgelq
subroutine zgelq(m, n, a, lda, t, tsize, work, lwork, info)
ZGELQ
Definition zgelq.f:172
zgemlq
subroutine zgemlq(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
ZGEMLQ
Definition zgemlq.f:169
zgemqr
subroutine zgemqr(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
ZGEMQR
Definition zgemqr.f:172
zgeqr
subroutine zgeqr(m, n, a, lda, t, tsize, work, lwork, info)
ZGEQR
Definition zgeqr.f:174
ztsqr01
subroutine ztsqr01(tssw, m, n, mb, nb, result)
ZTSQR01
Definition ztsqr01.f:82