LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
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

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.```
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*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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
subroutine slarfgp(n, alpha, x, incx, tau)
SLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition slarfgp.f:104
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
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: