LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ sorbdb4()

subroutine sorbdb4 ( integer  m,
integer  p,
integer  q,
real, dimension(ldx11,*)  x11,
integer  ldx11,
real, dimension(ldx21,*)  x21,
integer  ldx21,
real, dimension(*)  theta,
real, dimension(*)  phi,
real, dimension(*)  taup1,
real, dimension(*)  taup2,
real, dimension(*)  tauq1,
real, dimension(*)  phantom,
real, dimension(*)  work,
integer  lwork,
integer  info 
)

SORBDB4

Download SORBDB4 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SORBDB4 simultaneously bidiagonalizes the blocks of a tall and skinny
 matrix X with orthonormal columns:

                            [ B11 ]
      [ X11 ]   [ P1 |    ] [  0  ]
      [-----] = [---------] [-----] Q1**T .
      [ X21 ]   [    | P2 ] [ B21 ]
                            [  0  ]

 X11 is P-by-Q, and X21 is (M-P)-by-Q. M-Q must be no larger than P,
 M-P, or Q. Routines SORBDB1, SORBDB2, and SORBDB3 handle cases in
 which M-Q is not the minimum dimension.

 The orthogonal matrices P1, P2, and Q1 are P-by-P, (M-P)-by-(M-P),
 and (M-Q)-by-(M-Q), respectively. They are represented implicitly by
 Householder vectors.

 B11 and B12 are (M-Q)-by-(M-Q) bidiagonal matrices represented
 implicitly by angles THETA, PHI.
Parameters
[in]M
          M is INTEGER
           The number of rows X11 plus the number of rows in X21.
[in]P
          P is INTEGER
           The number of rows in X11. 0 <= P <= M.
[in]Q
          Q is INTEGER
           The number of columns in X11 and X21. 0 <= Q <= M and
           M-Q <= min(P,M-P,Q).
[in,out]X11
          X11 is REAL array, dimension (LDX11,Q)
           On entry, the top block of the matrix X to be reduced. On
           exit, the columns of tril(X11) specify reflectors for P1 and
           the rows of triu(X11,1) specify reflectors for Q1.
[in]LDX11
          LDX11 is INTEGER
           The leading dimension of X11. LDX11 >= P.
[in,out]X21
          X21 is REAL array, dimension (LDX21,Q)
           On entry, the bottom block of the matrix X to be reduced. On
           exit, the columns of tril(X21) specify reflectors for P2.
[in]LDX21
          LDX21 is INTEGER
           The leading dimension of X21. LDX21 >= M-P.
[out]THETA
          THETA is REAL array, dimension (Q)
           The entries of the bidiagonal blocks B11, B21 are defined by
           THETA and PHI. See Further Details.
[out]PHI
          PHI is REAL array, dimension (Q-1)
           The entries of the bidiagonal blocks B11, B21 are defined by
           THETA and PHI. See Further Details.
[out]TAUP1
          TAUP1 is REAL array, dimension (M-Q)
           The scalar factors of the elementary reflectors that define
           P1.
[out]TAUP2
          TAUP2 is REAL array, dimension (M-Q)
           The scalar factors of the elementary reflectors that define
           P2.
[out]TAUQ1
          TAUQ1 is REAL array, dimension (Q)
           The scalar factors of the elementary reflectors that define
           Q1.
[out]PHANTOM
          PHANTOM is REAL array, dimension (M)
           The routine computes an M-by-1 column vector Y that is
           orthogonal to the columns of [ X11; X21 ]. PHANTOM(1:P) and
           PHANTOM(P+1:M) contain Householder vectors for Y(1:P) and
           Y(P+1:M), respectively.
[out]WORK
          WORK is REAL array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
           The dimension of the array WORK. LWORK >= M-Q.

           If LWORK = -1, then a workspace query is assumed; the routine
           only calculates the optimal size of the WORK array, returns
           this value as the first entry of the WORK array, and no error
           message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
           = 0:  successful exit.
           < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  The upper-bidiagonal blocks B11, B21 are represented implicitly by
  angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
  in each bidiagonal band is a product of a sine or cosine of a THETA
  with a sine or cosine of a PHI. See [1] or SORCSD for details.

  P1, P2, and Q1 are represented as products of elementary reflectors.
  See SORCSD2BY1 for details on generating P1, P2, and Q1 using SORGQR
  and SORGLQ.
References:
[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.

Definition at line 211 of file sorbdb4.f.

214*
215* -- LAPACK computational routine --
216* -- LAPACK is a software package provided by Univ. of Tennessee, --
217* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
218*
219* .. Scalar Arguments ..
220 INTEGER INFO, LWORK, M, P, Q, LDX11, LDX21
221* ..
222* .. Array Arguments ..
223 REAL PHI(*), THETA(*)
224 REAL PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
225 $ WORK(*), X11(LDX11,*), X21(LDX21,*)
226* ..
227*
228* ====================================================================
229*
230* .. Parameters ..
231 REAL NEGONE, ONE, ZERO
232 parameter( negone = -1.0e0, one = 1.0e0, zero = 0.0e0 )
233* ..
234* .. Local Scalars ..
235 REAL C, S
236 INTEGER CHILDINFO, I, ILARF, IORBDB5, J, LLARF,
237 $ LORBDB5, LWORKMIN, LWORKOPT
238 LOGICAL LQUERY
239* ..
240* .. External Subroutines ..
241 EXTERNAL slarf, slarfgp, sorbdb5, srot, sscal, xerbla
242* ..
243* .. External Functions ..
244 REAL SNRM2
245 EXTERNAL snrm2
246* ..
247* .. Intrinsic Function ..
248 INTRINSIC atan2, cos, max, sin, sqrt
249* ..
250* .. Executable Statements ..
251*
252* Test input arguments
253*
254 info = 0
255 lquery = lwork .EQ. -1
256*
257 IF( m .LT. 0 ) THEN
258 info = -1
259 ELSE IF( p .LT. m-q .OR. m-p .LT. m-q ) THEN
260 info = -2
261 ELSE IF( q .LT. m-q .OR. q .GT. m ) THEN
262 info = -3
263 ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
264 info = -5
265 ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
266 info = -7
267 END IF
268*
269* Compute workspace
270*
271 IF( info .EQ. 0 ) THEN
272 ilarf = 2
273 llarf = max( q-1, p-1, m-p-1 )
274 iorbdb5 = 2
275 lorbdb5 = q
276 lworkopt = ilarf + llarf - 1
277 lworkopt = max( lworkopt, iorbdb5 + lorbdb5 - 1 )
278 lworkmin = lworkopt
279 work(1) = lworkopt
280 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
281 info = -14
282 END IF
283 END IF
284 IF( info .NE. 0 ) THEN
285 CALL xerbla( 'SORBDB4', -info )
286 RETURN
287 ELSE IF( lquery ) THEN
288 RETURN
289 END IF
290*
291* Reduce columns 1, ..., M-Q of X11 and X21
292*
293 DO i = 1, m-q
294*
295 IF( i .EQ. 1 ) THEN
296 DO j = 1, m
297 phantom(j) = zero
298 END DO
299 CALL sorbdb5( p, m-p, q, phantom(1), 1, phantom(p+1), 1,
300 $ x11, ldx11, x21, ldx21, work(iorbdb5),
301 $ lorbdb5, childinfo )
302 CALL sscal( p, negone, phantom(1), 1 )
303 CALL slarfgp( p, phantom(1), phantom(2), 1, taup1(1) )
304 CALL slarfgp( m-p, phantom(p+1), phantom(p+2), 1, taup2(1) )
305 theta(i) = atan2( phantom(1), phantom(p+1) )
306 c = cos( theta(i) )
307 s = sin( theta(i) )
308 phantom(1) = one
309 phantom(p+1) = one
310 CALL slarf( 'L', p, q, phantom(1), 1, taup1(1), x11, ldx11,
311 $ work(ilarf) )
312 CALL slarf( 'L', m-p, q, phantom(p+1), 1, taup2(1), x21,
313 $ ldx21, work(ilarf) )
314 ELSE
315 CALL sorbdb5( p-i+1, m-p-i+1, q-i+1, x11(i,i-1), 1,
316 $ x21(i,i-1), 1, x11(i,i), ldx11, x21(i,i),
317 $ ldx21, work(iorbdb5), lorbdb5, childinfo )
318 CALL sscal( p-i+1, negone, x11(i,i-1), 1 )
319 CALL slarfgp( p-i+1, x11(i,i-1), x11(i+1,i-1), 1, taup1(i) )
320 CALL slarfgp( m-p-i+1, x21(i,i-1), x21(i+1,i-1), 1,
321 $ taup2(i) )
322 theta(i) = atan2( x11(i,i-1), x21(i,i-1) )
323 c = cos( theta(i) )
324 s = sin( theta(i) )
325 x11(i,i-1) = one
326 x21(i,i-1) = one
327 CALL slarf( 'L', p-i+1, q-i+1, x11(i,i-1), 1, taup1(i),
328 $ x11(i,i), ldx11, work(ilarf) )
329 CALL slarf( 'L', m-p-i+1, q-i+1, x21(i,i-1), 1, taup2(i),
330 $ x21(i,i), ldx21, work(ilarf) )
331 END IF
332*
333 CALL srot( q-i+1, x11(i,i), ldx11, x21(i,i), ldx21, s, -c )
334 CALL slarfgp( q-i+1, x21(i,i), x21(i,i+1), ldx21, tauq1(i) )
335 c = x21(i,i)
336 x21(i,i) = one
337 CALL slarf( 'R', p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
338 $ x11(i+1,i), ldx11, work(ilarf) )
339 CALL slarf( 'R', m-p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
340 $ x21(i+1,i), ldx21, work(ilarf) )
341 IF( i .LT. m-q ) THEN
342 s = sqrt( snrm2( p-i, x11(i+1,i), 1 )**2
343 $ + snrm2( m-p-i, x21(i+1,i), 1 )**2 )
344 phi(i) = atan2( s, c )
345 END IF
346*
347 END DO
348*
349* Reduce the bottom-right portion of X11 to [ I 0 ]
350*
351 DO i = m - q + 1, p
352 CALL slarfgp( q-i+1, x11(i,i), x11(i,i+1), ldx11, tauq1(i) )
353 x11(i,i) = one
354 CALL slarf( 'R', p-i, q-i+1, x11(i,i), ldx11, tauq1(i),
355 $ x11(i+1,i), ldx11, work(ilarf) )
356 CALL slarf( 'R', q-p, q-i+1, x11(i,i), ldx11, tauq1(i),
357 $ x21(m-q+1,i), ldx21, work(ilarf) )
358 END DO
359*
360* Reduce the bottom-right portion of X21 to [ 0 I ]
361*
362 DO i = p + 1, q
363 CALL slarfgp( q-i+1, x21(m-q+i-p,i), x21(m-q+i-p,i+1), ldx21,
364 $ tauq1(i) )
365 x21(m-q+i-p,i) = one
366 CALL slarf( 'R', q-i, q-i+1, x21(m-q+i-p,i), ldx21, tauq1(i),
367 $ x21(m-q+i-p+1,i), ldx21, work(ilarf) )
368 END DO
369*
370 RETURN
371*
372* End of SORBDB4
373*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
slarf
subroutine slarf(side, m, n, v, incv, tau, c, ldc, work)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition slarf.f:124
slarfgp
subroutine slarfgp(n, alpha, x, incx, tau)
SLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition slarfgp.f:104
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
srot
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
sorbdb5
subroutine sorbdb5(m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2, ldq2, work, lwork, info)
SORBDB5
Definition sorbdb5.f:156
Here is the call graph for this function:
Here is the caller graph for this function: