83 SUBROUTINE ztsqr01(TSSW, M, N, MB, NB, RESULT)
95 DOUBLE PRECISION RESULT(6)
101 COMPLEX*16,
ALLOCATABLE :: AF(:,:), Q(:,:),
102 $ r(:,:), rwork(:), work( : ), t(:),
103 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:), lq(:,:)
106 DOUBLE PRECISION ZERO
107 COMPLEX*16 ONE, CZERO
108 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
111 LOGICAL TESTZEROS, TS
112 INTEGER INFO, J, K, L, LWORK, TSIZE, MNB
113 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
117 COMPLEX*16 TQUERY( 5 ), WORKQUERY
120 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY
123 EXTERNAL dlamch, zlange, zlansy, lsame, ilaenv
131 COMMON / srnamc / srnamt
134 DATA iseed / 1988, 1989, 1990, 1991 /
138 ts = lsame(tssw,
'TS')
144 eps = dlamch(
'Epsilon' )
152 ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
154 $ d(n,m), df(n,m), lq(l,n) )
159 CALL zlarnv( 2, iseed, m, a( 1, j ) )
164 CALL zlarnv( 2, iseed, m/2, a( m/4, j ) )
168 CALL zlacpy(
'Full', m, n, a, m, af, m )
174 CALL zgeqr( m, n, af, m, tquery, -1, workquery, -1, info )
175 tsize = int( tquery( 1 ) )
176 lwork = int( workquery )
177 CALL zgemqr(
'L',
'N', m, m, k, af, m, tquery, tsize, cf, m,
178 $ workquery, -1, info)
179 lwork = max( lwork, int( workquery ) )
180 CALL zgemqr(
'L',
'N', m, n, k, af, m, tquery, tsize, cf, m,
181 $ workquery, -1, info)
182 lwork = max( lwork, int( workquery ) )
183 CALL zgemqr(
'L',
'C', m, n, k, af, m, tquery, tsize, cf, m,
184 $ workquery, -1, info)
185 lwork = max( lwork, int( workquery ) )
186 CALL zgemqr(
'R',
'N', n, m, k, af, m, tquery, tsize, df, n,
187 $ workquery, -1, info)
188 lwork = max( lwork, int( workquery ) )
189 CALL zgemqr(
'R',
'C', n, m, k, af, m, tquery, tsize, df, n,
190 $ workquery, -1, info)
191 lwork = max( lwork, int( workquery ) )
192 ALLOCATE ( t( tsize ) )
193 ALLOCATE ( work( lwork ) )
195 CALL zgeqr( m, n, af, m, t, tsize, work, lwork, info )
199 CALL zlaset(
'Full', m, m, czero, one, q, m )
201 CALL zgemqr(
'L',
'N', m, m, k, af, m, t, tsize, q, m,
202 $ work, lwork, info )
206 CALL zlaset(
'Full', m, n, czero, czero, r, m )
207 CALL zlacpy(
'Upper', m, n, af, m, r, m )
211 CALL zgemm(
'C',
'N', m, n, m, -one, q, m, a, m, one, r, m )
212 anorm = zlange(
'1', m, n, a, m, rwork )
213 resid = zlange(
'1', m, n, r, m, rwork )
214 IF( anorm.GT.zero )
THEN
215 result( 1 ) = resid / (eps*max(1,m)*anorm)
222 CALL zlaset(
'Full', m, m, czero, one, r, m )
223 CALL zherk(
'U',
'C', m, m, dreal(-one), q, m, dreal(one), r, m )
224 resid = zlansy(
'1',
'Upper', m, r, m, rwork )
225 result( 2 ) = resid / (eps*max(1,m))
230 CALL zlarnv( 2, iseed, m, c( 1, j ) )
232 cnorm = zlange(
'1', m, n, c, m, rwork)
233 CALL zlacpy(
'Full', m, n, c, m, cf, m )
238 CALL zgemqr(
'L',
'N', m, n, k, af, m, t, tsize, cf, m,
243 CALL zgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
244 resid = zlange(
'1', m, n, cf, m, rwork )
245 IF( cnorm.GT.zero )
THEN
246 result( 3 ) = resid / (eps*max(1,m)*cnorm)
253 CALL zlacpy(
'Full', m, n, c, m, cf, m )
258 CALL zgemqr(
'L',
'C', m, n, k, af, m, t, tsize, cf, m,
263 CALL zgemm(
'C',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
264 resid = zlange(
'1', m, n, cf, m, rwork )
265 IF( cnorm.GT.zero )
THEN
266 result( 4 ) = resid / (eps*max(1,m)*cnorm)
274 CALL zlarnv( 2, iseed, n, d( 1, j ) )
276 dnorm = zlange(
'1', n, m, d, n, rwork)
277 CALL zlacpy(
'Full', n, m, d, n, df, n )
282 CALL zgemqr(
'R',
'N', n, m, k, af, m, t, tsize, df, n,
287 CALL zgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
288 resid = zlange(
'1', n, m, df, n, rwork )
289 IF( dnorm.GT.zero )
THEN
290 result( 5 ) = resid / (eps*max(1,m)*dnorm)
297 CALL zlacpy(
'Full', n, m, d, n, df, n )
301 CALL zgemqr(
'R',
'C', n, m, k, af, m, t, tsize, df, n,
306 CALL zgemm(
'N',
'C', n, m, m, -one, d, n, q, m, one, df, n )
307 resid = zlange(
'1', n, m, df, n, rwork )
308 IF( cnorm.GT.zero )
THEN
309 result( 6 ) = resid / (eps*max(1,m)*dnorm)
317 CALL zgelq( m, n, af, m, tquery, -1, workquery, -1, info )
318 tsize = int( tquery( 1 ) )
319 lwork = int( workquery )
320 CALL zgemlq(
'R',
'N', n, n, k, af, m, tquery, tsize, q, n,
321 $ workquery, -1, info )
322 lwork = max( lwork, int( workquery ) )
323 CALL zgemlq(
'L',
'N', n, m, k, af, m, tquery, tsize, df, n,
324 $ workquery, -1, info)
325 lwork = max( lwork, int( workquery ) )
326 CALL zgemlq(
'L',
'C', n, m, k, af, m, tquery, tsize, df, n,
327 $ workquery, -1, info)
328 lwork = max( lwork, int( workquery ) )
329 CALL zgemlq(
'R',
'N', m, n, k, af, m, tquery, tsize, cf, m,
330 $ workquery, -1, info)
331 lwork = max( lwork, int( workquery ) )
332 CALL zgemlq(
'R',
'C', m, n, k, af, m, tquery, tsize, cf, m,
333 $ workquery, -1, info)
334 lwork = max( lwork, int( workquery ) )
335 ALLOCATE ( t( tsize ) )
336 ALLOCATE ( work( lwork ) )
338 CALL zgelq( m, n, af, m, t, tsize, work, lwork, info )
343 CALL zlaset(
'Full', n, n, czero, one, q, n )
345 CALL zgemlq(
'R',
'N', n, n, k, af, m, t, tsize, q, n,
346 $ work, lwork, info )
350 CALL zlaset(
'Full', m, n, czero, czero, lq, l )
351 CALL zlacpy(
'Lower', m, n, af, m, lq, l )
355 CALL zgemm(
'N',
'C', m, n, n, -one, a, m, q, n, one, lq, l )
356 anorm = zlange(
'1', m, n, a, m, rwork )
357 resid = zlange(
'1', m, n, lq, l, rwork )
358 IF( anorm.GT.zero )
THEN
359 result( 1 ) = resid / (eps*max(1,n)*anorm)
366 CALL zlaset(
'Full', n, n, czero, one, lq, l )
367 CALL zherk(
'U',
'C', n, n, dreal(-one), q, n, dreal(one), lq, l)
368 resid = zlansy(
'1',
'Upper', n, lq, l, rwork )
369 result( 2 ) = resid / (eps*max(1,n))
374 CALL zlarnv( 2, iseed, n, d( 1, j ) )
376 dnorm = zlange(
'1', n, m, d, n, rwork)
377 CALL zlacpy(
'Full', n, m, d, n, df, n )
381 CALL zgemlq(
'L',
'N', n, m, k, af, m, t, tsize, df, n,
386 CALL zgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
387 resid = zlange(
'1', n, m, df, n, rwork )
388 IF( dnorm.GT.zero )
THEN
389 result( 3 ) = resid / (eps*max(1,n)*dnorm)
396 CALL zlacpy(
'Full', n, m, d, n, df, n )
400 CALL zgemlq(
'L',
'C', n, m, k, af, m, t, tsize, df, n,
405 CALL zgemm(
'C',
'N', n, m, n, -one, q, n, d, n, one, df, n )
406 resid = zlange(
'1', n, m, df, n, rwork )
407 IF( dnorm.GT.zero )
THEN
408 result( 4 ) = resid / (eps*max(1,n)*dnorm)
416 CALL zlarnv( 2, iseed, m, c( 1, j ) )
418 cnorm = zlange(
'1', m, n, c, m, rwork)
419 CALL zlacpy(
'Full', m, n, c, m, cf, m )
423 CALL zgemlq(
'R',
'N', m, n, k, af, m, t, tsize, cf, m,
428 CALL zgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
429 resid = zlange(
'1', n, m, df, n, rwork )
430 IF( cnorm.GT.zero )
THEN
431 result( 5 ) = resid / (eps*max(1,n)*cnorm)
438 CALL zlacpy(
'Full', m, n, c, m, cf, m )
442 CALL zgemlq(
'R',
'C', m, n, k, af, m, t, tsize, cf, m,
447 CALL zgemm(
'N',
'C', m, n, n, -one, c, m, q, n, one, cf, m )
448 resid = zlange(
'1', m, n, cf, m, rwork )
449 IF( cnorm.GT.zero )
THEN
450 result( 6 ) = resid / (eps*max(1,n)*cnorm)
459 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
subroutine zgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine ztsqr01(TSSW, M, N, MB, NB, RESULT)
ZTSQR01
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...
subroutine zgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
subroutine zgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
subroutine zgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.