212 SUBROUTINE zunbdb4( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI,
213 $ taup1, taup2, tauq1, phantom, work, lwork,
222 INTEGER info, lwork, m, p, q, ldx11, ldx21
225 DOUBLE PRECISION phi(*), theta(*)
226 COMPLEX*16 phantom(*), taup1(*), taup2(*), tauq1(*),
227 $ work(*), x11(ldx11,*), x21(ldx21,*)
233 COMPLEX*16 negone, one, zero
234 parameter( negone = (-1.0d0,0.0d0), one = (1.0d0,0.0d0),
235 $ zero = (0.0d0,0.0d0) )
238 DOUBLE PRECISION c, s
239 INTEGER childinfo, i, ilarf, iorbdb5,
j, llarf,
240 $ lorbdb5, lworkmin, lworkopt
251 INTRINSIC atan2, cos, max, sin, sqrt
258 lquery = lwork .EQ. -1
262 ELSE IF( p .LT. m-q .OR. m-p .LT. m-q )
THEN
264 ELSE IF( q .LT. m-q .OR. q .GT. m )
THEN
266 ELSE IF( ldx11 .LT. max( 1, p ) )
THEN
268 ELSE IF( ldx21 .LT. max( 1, m-p ) )
THEN
274 IF( info .EQ. 0 )
THEN
276 llarf = max( q-1, p-1, m-p-1 )
279 lworkopt = ilarf + llarf - 1
280 lworkopt = max( lworkopt, iorbdb5 + lorbdb5 - 1 )
283 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
287 IF( info .NE. 0 )
THEN
288 CALL
xerbla(
'ZUNBDB4', -info )
290 ELSE IF( lquery )
THEN
302 CALL
zunbdb5( p, m-p, q, phantom(1), 1, phantom(p+1), 1,
303 $ x11, ldx11, x21, ldx21, work(iorbdb5),
304 $ lorbdb5, childinfo )
305 CALL
zscal( p, negone, phantom(1), 1 )
306 CALL
zlarfgp( p, phantom(1), phantom(2), 1, taup1(1) )
307 CALL
zlarfgp( m-p, phantom(p+1), phantom(p+2), 1, taup2(1) )
308 theta(i) = atan2( dble( phantom(1) ), dble( phantom(p+1) ) )
313 CALL
zlarf(
'L', p, q, phantom(1), 1, dconjg(taup1(1)), x11,
314 $ ldx11, work(ilarf) )
315 CALL
zlarf(
'L', m-p, q, phantom(p+1), 1, dconjg(taup2(1)),
316 $ x21, ldx21, work(ilarf) )
318 CALL
zunbdb5( p-i+1, m-p-i+1, q-i+1, x11(i,i-1), 1,
319 $ x21(i,i-1), 1, x11(i,i), ldx11, x21(i,i),
320 $ ldx21, work(iorbdb5), lorbdb5, childinfo )
321 CALL
zscal( p-i+1, negone, x11(i,i-1), 1 )
322 CALL
zlarfgp( p-i+1, x11(i,i-1), x11(i+1,i-1), 1, taup1(i) )
323 CALL
zlarfgp( m-p-i+1, x21(i,i-1), x21(i+1,i-1), 1,
325 theta(i) = atan2( dble( x11(i,i-1) ), dble( x21(i,i-1) ) )
330 CALL
zlarf(
'L', p-i+1, q-i+1, x11(i,i-1), 1,
331 $ dconjg(taup1(i)), x11(i,i), ldx11, work(ilarf) )
332 CALL
zlarf(
'L', m-p-i+1, q-i+1, x21(i,i-1), 1,
333 $ dconjg(taup2(i)), x21(i,i), ldx21, work(ilarf) )
336 CALL
zdrot( q-i+1, x11(i,i), ldx11, x21(i,i), ldx21, s, -c )
337 CALL
zlacgv( q-i+1, x21(i,i), ldx21 )
338 CALL
zlarfgp( q-i+1, x21(i,i), x21(i,i+1), ldx21, tauq1(i) )
341 CALL
zlarf(
'R', p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
342 $ x11(i+1,i), ldx11, work(ilarf) )
343 CALL
zlarf(
'R', m-p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
344 $ x21(i+1,i), ldx21, work(ilarf) )
345 CALL
zlacgv( q-i+1, x21(i,i), ldx21 )
346 IF( i .LT. m-q )
THEN
347 s = sqrt(
dznrm2( p-i, x11(i+1,i), 1, x11(i+1,i),
348 $ 1 )**2 +
dznrm2( m-p-i, x21(i+1,i), 1, x21(i+1,i),
350 phi(i) = atan2( s, c )
358 CALL
zlacgv( q-i+1, x11(i,i), ldx11 )
359 CALL
zlarfgp( q-i+1, x11(i,i), x11(i,i+1), ldx11, tauq1(i) )
361 CALL
zlarf(
'R', p-i, q-i+1, x11(i,i), ldx11, tauq1(i),
362 $ x11(i+1,i), ldx11, work(ilarf) )
363 CALL
zlarf(
'R', q-p, q-i+1, x11(i,i), ldx11, tauq1(i),
364 $ x21(m-q+1,i), ldx21, work(ilarf) )
365 CALL
zlacgv( q-i+1, x11(i,i), ldx11 )
371 CALL
zlacgv( q-i+1, x21(m-q+i-p,i), ldx21 )
372 CALL
zlarfgp( q-i+1, x21(m-q+i-p,i), x21(m-q+i-p,i+1), ldx21,
375 CALL
zlarf(
'R', q-i, q-i+1, x21(m-q+i-p,i), ldx21, tauq1(i),
376 $ x21(m-q+i-p+1,i), ldx21, work(ilarf) )
377 CALL
zlacgv( q-i+1, x21(m-q+i-p,i), ldx21 )