212 SUBROUTINE dorbdb4( 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 DOUBLE PRECISION phantom(*), taup1(*), taup2(*), tauq1(*),
227 $ work(*), x11(ldx11,*), x21(ldx21,*)
233 DOUBLE PRECISION negone, one, zero
234 parameter( negone = -1.0d0, one = 1.0d0, zero = 0.0d0 )
237 DOUBLE PRECISION c, s
238 INTEGER childinfo, i, ilarf, iorbdb5,
j, llarf,
239 $ lorbdb5, lworkmin, lworkopt
246 DOUBLE PRECISION dnrm2
250 INTRINSIC atan2, cos, max, sin, sqrt
257 lquery = lwork .EQ. -1
261 ELSE IF( p .LT. m-q .OR. m-p .LT. m-q )
THEN
263 ELSE IF( q .LT. m-q .OR. q .GT. m )
THEN
265 ELSE IF( ldx11 .LT. max( 1, p ) )
THEN
267 ELSE IF( ldx21 .LT. max( 1, m-p ) )
THEN
273 IF( info .EQ. 0 )
THEN
275 llarf = max( q-1, p-1, m-p-1 )
278 lworkopt = ilarf + llarf - 1
279 lworkopt = max( lworkopt, iorbdb5 + lorbdb5 - 1 )
282 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
286 IF( info .NE. 0 )
THEN
287 CALL
xerbla(
'DORBDB4', -info )
289 ELSE IF( lquery )
THEN
301 CALL
dorbdb5( p, m-p, q, phantom(1), 1, phantom(p+1), 1,
302 $ x11, ldx11, x21, ldx21, work(iorbdb5),
303 $ lorbdb5, childinfo )
304 CALL
dscal( p, negone, phantom(1), 1 )
305 CALL
dlarfgp( p, phantom(1), phantom(2), 1, taup1(1) )
306 CALL
dlarfgp( m-p, phantom(p+1), phantom(p+2), 1, taup2(1) )
307 theta(i) = atan2( phantom(1), phantom(p+1) )
312 CALL
dlarf(
'L', p, q, phantom(1), 1, taup1(1), x11, ldx11,
314 CALL
dlarf(
'L', m-p, q, phantom(p+1), 1, taup2(1), x21,
315 $ ldx21, work(ilarf) )
317 CALL
dorbdb5( p-i+1, m-p-i+1, q-i+1, x11(i,i-1), 1,
318 $ x21(i,i-1), 1, x11(i,i), ldx11, x21(i,i),
319 $ ldx21, work(iorbdb5), lorbdb5, childinfo )
320 CALL
dscal( p-i+1, negone, x11(i,i-1), 1 )
321 CALL
dlarfgp( p-i+1, x11(i,i-1), x11(i+1,i-1), 1, taup1(i) )
322 CALL
dlarfgp( m-p-i+1, x21(i,i-1), x21(i+1,i-1), 1,
324 theta(i) = atan2( x11(i,i-1), x21(i,i-1) )
329 CALL
dlarf(
'L', p-i+1, q-i+1, x11(i,i-1), 1, taup1(i),
330 $ x11(i,i), ldx11, work(ilarf) )
331 CALL
dlarf(
'L', m-p-i+1, q-i+1, x21(i,i-1), 1, taup2(i),
332 $ x21(i,i), ldx21, work(ilarf) )
335 CALL
drot( q-i+1, x11(i,i), ldx11, x21(i,i), ldx21, s, -c )
336 CALL
dlarfgp( q-i+1, x21(i,i), x21(i,i+1), ldx21, tauq1(i) )
339 CALL
dlarf(
'R', p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
340 $ x11(i+1,i), ldx11, work(ilarf) )
341 CALL
dlarf(
'R', m-p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
342 $ x21(i+1,i), ldx21, work(ilarf) )
343 IF( i .LT. m-q )
THEN
344 s = sqrt(
dnrm2( p-i, x11(i+1,i), 1, x11(i+1,i),
345 $ 1 )**2 +
dnrm2( m-p-i, x21(i+1,i), 1, x21(i+1,i),
347 phi(i) = atan2( s, c )
355 CALL
dlarfgp( q-i+1, x11(i,i), x11(i,i+1), ldx11, tauq1(i) )
357 CALL
dlarf(
'R', p-i, q-i+1, x11(i,i), ldx11, tauq1(i),
358 $ x11(i+1,i), ldx11, work(ilarf) )
359 CALL
dlarf(
'R', q-p, q-i+1, x11(i,i), ldx11, tauq1(i),
360 $ x21(m-q+1,i), ldx21, work(ilarf) )
366 CALL
dlarfgp( q-i+1, x21(m-q+i-p,i), x21(m-q+i-p,i+1), ldx21,
369 CALL
dlarf(
'R', q-i, q-i+1, x21(m-q+i-p,i), ldx21, tauq1(i),
370 $ x21(m-q+i-p+1,i), ldx21, work(ilarf) )