213 SUBROUTINE sorbdb4( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI,
214 $ taup1, taup2, tauq1, phantom, work, lwork,
223 INTEGER info, lwork, m, p, q, ldx11, ldx21
226 REAL phi(*), theta(*)
227 REAL phantom(*), taup1(*), taup2(*), tauq1(*),
228 $ work(*), x11(ldx11,*), x21(ldx21,*)
234 REAL negone, one, zero
235 parameter( negone = -1.0e0, one = 1.0e0, zero = 0.0e0 )
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(
'SORBDB4', -info )
290 ELSE IF( lquery )
THEN
302 CALL
sorbdb5( p, m-p, q, phantom(1), 1, phantom(p+1), 1,
303 $ x11, ldx11, x21, ldx21, work(iorbdb5),
304 $ lorbdb5, childinfo )
305 CALL
sscal( p, negone, phantom(1), 1 )
306 CALL
slarfgp( p, phantom(1), phantom(2), 1, taup1(1) )
307 CALL
slarfgp( m-p, phantom(p+1), phantom(p+2), 1, taup2(1) )
308 theta(i) = atan2( phantom(1), phantom(p+1) )
313 CALL
slarf(
'L', p, q, phantom(1), 1, taup1(1), x11, ldx11,
315 CALL
slarf(
'L', m-p, q, phantom(p+1), 1, taup2(1), x21,
316 $ ldx21, work(ilarf) )
318 CALL
sorbdb5( 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
sscal( p-i+1, negone, x11(i,i-1), 1 )
322 CALL
slarfgp( p-i+1, x11(i,i-1), x11(i+1,i-1), 1, taup1(i) )
323 CALL
slarfgp( m-p-i+1, x21(i,i-1), x21(i+1,i-1), 1,
325 theta(i) = atan2( x11(i,i-1), x21(i,i-1) )
330 CALL
slarf(
'L', p-i+1, q-i+1, x11(i,i-1), 1, taup1(i),
331 $ x11(i,i), ldx11, work(ilarf) )
332 CALL
slarf(
'L', m-p-i+1, q-i+1, x21(i,i-1), 1, taup2(i),
333 $ x21(i,i), ldx21, work(ilarf) )
336 CALL
srot( q-i+1, x11(i,i), ldx11, x21(i,i), ldx21, s, -c )
337 CALL
slarfgp( q-i+1, x21(i,i), x21(i,i+1), ldx21, tauq1(i) )
340 CALL
slarf(
'R', p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
341 $ x11(i+1,i), ldx11, work(ilarf) )
342 CALL
slarf(
'R', m-p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
343 $ x21(i+1,i), ldx21, work(ilarf) )
344 IF( i .LT. m-q )
THEN
345 s = sqrt(
snrm2( p-i, x11(i+1,i), 1, x11(i+1,i),
346 $ 1 )**2 +
snrm2( m-p-i, x21(i+1,i), 1, x21(i+1,i),
348 phi(i) = atan2( s, c )
356 CALL
slarfgp( q-i+1, x11(i,i), x11(i,i+1), ldx11, tauq1(i) )
358 CALL
slarf(
'R', p-i, q-i+1, x11(i,i), ldx11, tauq1(i),
359 $ x11(i+1,i), ldx11, work(ilarf) )
360 CALL
slarf(
'R', q-p, q-i+1, x11(i,i), ldx11, tauq1(i),
361 $ x21(m-q+1,i), ldx21, work(ilarf) )
367 CALL
slarfgp( q-i+1, x21(m-q+i-p,i), x21(m-q+i-p,i+1), ldx21,
370 CALL
slarf(
'R', q-i, q-i+1, x21(m-q+i-p,i), ldx21, tauq1(i),
371 $ x21(m-q+i-p+1,i), ldx21, work(ilarf) )