176 SUBROUTINE dlatm6( TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA,
177 $ beta, wx, wy, s, dif )
185 INTEGER LDA, LDX, LDY, N, TYPE
186 DOUBLE PRECISION ALPHA, BETA, WX, WY
189 DOUBLE PRECISION A( lda, * ), B( lda, * ), DIF( * ), S( * ),
190 $ x( ldx, * ), y( ldy, * )
196 DOUBLE PRECISION ZERO, ONE, TWO, THREE
197 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
204 DOUBLE PRECISION WORK( 100 ), Z( 12, 12 )
221 a( i, i ) = dble( i ) + alpha
233 CALL
dlacpy(
'F', n, n, b, lda, y, ldy )
241 CALL
dlacpy(
'F', n, n, b, lda, x, ldx )
258 a( 1, 3 ) = wx*a( 1, 1 ) + wy*a( 3, 3 )
259 a( 2, 3 ) = -wx*a( 2, 2 ) + wy*a( 3, 3 )
260 a( 1, 4 ) = wx*a( 1, 1 ) - wy*a( 4, 4 )
261 a( 2, 4 ) = wx*a( 2, 2 ) - wy*a( 4, 4 )
262 a( 1, 5 ) = -wx*a( 1, 1 ) + wy*a( 5, 5 )
263 a( 2, 5 ) = wx*a( 2, 2 ) + wy*a( 5, 5 )
264 ELSE IF( type.EQ.2 )
THEN
265 a( 1, 3 ) = two*wx + wy
267 a( 1, 4 ) = -wy*( two+alpha+beta )
268 a( 2, 4 ) = two*wx - wy*( two+alpha+beta )
269 a( 1, 5 ) = -two*wx + wy*( alpha-beta )
270 a( 2, 5 ) = wy*( alpha-beta )
274 a( 2, 2 ) = a( 1, 1 )
276 a( 4, 4 ) = one + alpha
277 a( 4, 5 ) = one + beta
278 a( 5, 4 ) = -a( 4, 5 )
279 a( 5, 5 ) = a( 4, 4 )
286 s( 1 ) = one / sqrt( ( one+three*wy*wy ) /
287 $ ( one+a( 1, 1 )*a( 1, 1 ) ) )
288 s( 2 ) = one / sqrt( ( one+three*wy*wy ) /
289 $ ( one+a( 2, 2 )*a( 2, 2 ) ) )
290 s( 3 ) = one / sqrt( ( one+two*wx*wx ) /
291 $ ( one+a( 3, 3 )*a( 3, 3 ) ) )
292 s( 4 ) = one / sqrt( ( one+two*wx*wx ) /
293 $ ( one+a( 4, 4 )*a( 4, 4 ) ) )
294 s( 5 ) = one / sqrt( ( one+two*wx*wx ) /
295 $ ( one+a( 5, 5 )*a( 5, 5 ) ) )
297 CALL
dlakf2( 1, 4, a, lda, a( 2, 2 ), b, b( 2, 2 ), z, 12 )
298 CALL
dgesvd(
'N',
'N', 8, 8, z, 12, work, work( 9 ), 1,
299 $ work( 10 ), 1, work( 11 ), 40, info )
302 CALL
dlakf2( 4, 1, a, lda, a( 5, 5 ), b, b( 5, 5 ), z, 12 )
303 CALL
dgesvd(
'N',
'N', 8, 8, z, 12, work, work( 9 ), 1,
304 $ work( 10 ), 1, work( 11 ), 40, info )
307 ELSE IF( type.EQ.2 )
THEN
309 s( 1 ) = one / sqrt( one / three+wy*wy )
311 s( 3 ) = one / sqrt( one / two+wx*wx )
312 s( 4 ) = one / sqrt( ( one+two*wx*wx ) /
313 $ ( one+( one+alpha )*( one+alpha )+( one+beta )*( one+
317 CALL
dlakf2( 2, 3, a, lda, a( 3, 3 ), b, b( 3, 3 ), z, 12 )
318 CALL
dgesvd(
'N',
'N', 12, 12, z, 12, work, work( 13 ), 1,
319 $ work( 14 ), 1, work( 15 ), 60, info )
320 dif( 1 ) = work( 12 )
322 CALL
dlakf2( 3, 2, a, lda, a( 4, 4 ), b, b( 4, 4 ), z, 12 )
323 CALL
dgesvd(
'N',
'N', 12, 12, z, 12, work, work( 13 ), 1,
324 $ work( 14 ), 1, work( 15 ), 60, info )
325 dif( 5 ) = work( 12 )
subroutine dlakf2(M, N, A, LDA, B, D, E, Z, LDZ)
DLAKF2
subroutine dlatm6(TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, BETA, WX, WY, S, DIF)
DLATM6
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
DGESVD computes the singular value decomposition (SVD) for GE matrices