LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clals0.f
Go to the documentation of this file.
1*> \brief \b CLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer SVD approach. Used by sgelsd.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CLALS0 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clals0.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clals0.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clals0.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
22* PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
23* POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO )
24*
25* .. Scalar Arguments ..
26* INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
27* $ LDGNUM, NL, NR, NRHS, SQRE
28* REAL C, S
29* ..
30* .. Array Arguments ..
31* INTEGER GIVCOL( LDGCOL, * ), PERM( * )
32* REAL DIFL( * ), DIFR( LDGNUM, * ),
33* $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
34* $ RWORK( * ), Z( * )
35* COMPLEX B( LDB, * ), BX( LDBX, * )
36* ..
37*
38*
39*> \par Purpose:
40* =============
41*>
42*> \verbatim
43*>
44*> CLALS0 applies back the multiplying factors of either the left or the
45*> right singular vector matrix of a diagonal matrix appended by a row
46*> to the right hand side matrix B in solving the least squares problem
47*> using the divide-and-conquer SVD approach.
48*>
49*> For the left singular vector matrix, three types of orthogonal
50*> matrices are involved:
51*>
52*> (1L) Givens rotations: the number of such rotations is GIVPTR; the
53*> pairs of columns/rows they were applied to are stored in GIVCOL;
54*> and the C- and S-values of these rotations are stored in GIVNUM.
55*>
56*> (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
57*> row, and for J=2:N, PERM(J)-th row of B is to be moved to the
58*> J-th row.
59*>
60*> (3L) The left singular vector matrix of the remaining matrix.
61*>
62*> For the right singular vector matrix, four types of orthogonal
63*> matrices are involved:
64*>
65*> (1R) The right singular vector matrix of the remaining matrix.
66*>
67*> (2R) If SQRE = 1, one extra Givens rotation to generate the right
68*> null space.
69*>
70*> (3R) The inverse transformation of (2L).
71*>
72*> (4R) The inverse transformation of (1L).
73*> \endverbatim
74*
75* Arguments:
76* ==========
77*
78*> \param[in] ICOMPQ
79*> \verbatim
80*> ICOMPQ is INTEGER
81*> Specifies whether singular vectors are to be computed in
82*> factored form:
83*> = 0: Left singular vector matrix.
84*> = 1: Right singular vector matrix.
85*> \endverbatim
86*>
87*> \param[in] NL
88*> \verbatim
89*> NL is INTEGER
90*> The row dimension of the upper block. NL >= 1.
91*> \endverbatim
92*>
93*> \param[in] NR
94*> \verbatim
95*> NR is INTEGER
96*> The row dimension of the lower block. NR >= 1.
97*> \endverbatim
98*>
99*> \param[in] SQRE
100*> \verbatim
101*> SQRE is INTEGER
102*> = 0: the lower block is an NR-by-NR square matrix.
103*> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
104*>
105*> The bidiagonal matrix has row dimension N = NL + NR + 1,
106*> and column dimension M = N + SQRE.
107*> \endverbatim
108*>
109*> \param[in] NRHS
110*> \verbatim
111*> NRHS is INTEGER
112*> The number of columns of B and BX. NRHS must be at least 1.
113*> \endverbatim
114*>
115*> \param[in,out] B
116*> \verbatim
117*> B is COMPLEX array, dimension ( LDB, NRHS )
118*> On input, B contains the right hand sides of the least
119*> squares problem in rows 1 through M. On output, B contains
120*> the solution X in rows 1 through N.
121*> \endverbatim
122*>
123*> \param[in] LDB
124*> \verbatim
125*> LDB is INTEGER
126*> The leading dimension of B. LDB must be at least
127*> max(1,MAX( M, N ) ).
128*> \endverbatim
129*>
130*> \param[out] BX
131*> \verbatim
132*> BX is COMPLEX array, dimension ( LDBX, NRHS )
133*> \endverbatim
134*>
135*> \param[in] LDBX
136*> \verbatim
137*> LDBX is INTEGER
138*> The leading dimension of BX.
139*> \endverbatim
140*>
141*> \param[in] PERM
142*> \verbatim
143*> PERM is INTEGER array, dimension ( N )
144*> The permutations (from deflation and sorting) applied
145*> to the two blocks.
146*> \endverbatim
147*>
148*> \param[in] GIVPTR
149*> \verbatim
150*> GIVPTR is INTEGER
151*> The number of Givens rotations which took place in this
152*> subproblem.
153*> \endverbatim
154*>
155*> \param[in] GIVCOL
156*> \verbatim
157*> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
158*> Each pair of numbers indicates a pair of rows/columns
159*> involved in a Givens rotation.
160*> \endverbatim
161*>
162*> \param[in] LDGCOL
163*> \verbatim
164*> LDGCOL is INTEGER
165*> The leading dimension of GIVCOL, must be at least N.
166*> \endverbatim
167*>
168*> \param[in] GIVNUM
169*> \verbatim
170*> GIVNUM is REAL array, dimension ( LDGNUM, 2 )
171*> Each number indicates the C or S value used in the
172*> corresponding Givens rotation.
173*> \endverbatim
174*>
175*> \param[in] LDGNUM
176*> \verbatim
177*> LDGNUM is INTEGER
178*> The leading dimension of arrays DIFR, POLES and
179*> GIVNUM, must be at least K.
180*> \endverbatim
181*>
182*> \param[in] POLES
183*> \verbatim
184*> POLES is REAL array, dimension ( LDGNUM, 2 )
185*> On entry, POLES(1:K, 1) contains the new singular
186*> values obtained from solving the secular equation, and
187*> POLES(1:K, 2) is an array containing the poles in the secular
188*> equation.
189*> \endverbatim
190*>
191*> \param[in] DIFL
192*> \verbatim
193*> DIFL is REAL array, dimension ( K ).
194*> On entry, DIFL(I) is the distance between I-th updated
195*> (undeflated) singular value and the I-th (undeflated) old
196*> singular value.
197*> \endverbatim
198*>
199*> \param[in] DIFR
200*> \verbatim
201*> DIFR is REAL array, dimension ( LDGNUM, 2 ).
202*> On entry, DIFR(I, 1) contains the distances between I-th
203*> updated (undeflated) singular value and the I+1-th
204*> (undeflated) old singular value. And DIFR(I, 2) is the
205*> normalizing factor for the I-th right singular vector.
206*> \endverbatim
207*>
208*> \param[in] Z
209*> \verbatim
210*> Z is REAL array, dimension ( K )
211*> Contain the components of the deflation-adjusted updating row
212*> vector.
213*> \endverbatim
214*>
215*> \param[in] K
216*> \verbatim
217*> K is INTEGER
218*> Contains the dimension of the non-deflated matrix,
219*> This is the order of the related secular equation. 1 <= K <=N.
220*> \endverbatim
221*>
222*> \param[in] C
223*> \verbatim
224*> C is REAL
225*> C contains garbage if SQRE =0 and the C-value of a Givens
226*> rotation related to the right null space if SQRE = 1.
227*> \endverbatim
228*>
229*> \param[in] S
230*> \verbatim
231*> S is REAL
232*> S contains garbage if SQRE =0 and the S-value of a Givens
233*> rotation related to the right null space if SQRE = 1.
234*> \endverbatim
235*>
236*> \param[out] RWORK
237*> \verbatim
238*> RWORK is REAL array, dimension
239*> ( K*(1+NRHS) + 2*NRHS )
240*> \endverbatim
241*>
242*> \param[out] INFO
243*> \verbatim
244*> INFO is INTEGER
245*> = 0: successful exit.
246*> < 0: if INFO = -i, the i-th argument had an illegal value.
247*> \endverbatim
248*
249* Authors:
250* ========
251*
252*> \author Univ. of Tennessee
253*> \author Univ. of California Berkeley
254*> \author Univ. of Colorado Denver
255*> \author NAG Ltd.
256*
257*> \ingroup lals0
258*
259*> \par Contributors:
260* ==================
261*>
262*> Ming Gu and Ren-Cang Li, Computer Science Division, University of
263*> California at Berkeley, USA \n
264*> Osni Marques, LBNL/NERSC, USA \n
265*
266* =====================================================================
267 SUBROUTINE clals0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
268 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
269 $ POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO )
270*
271* -- LAPACK computational routine --
272* -- LAPACK is a software package provided by Univ. of Tennessee, --
273* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
274*
275* .. Scalar Arguments ..
276 INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
277 $ LDGNUM, NL, NR, NRHS, SQRE
278 REAL C, S
279* ..
280* .. Array Arguments ..
281 INTEGER GIVCOL( LDGCOL, * ), PERM( * )
282 REAL DIFL( * ), DIFR( LDGNUM, * ),
283 $ givnum( ldgnum, * ), poles( ldgnum, * ),
284 $ rwork( * ), z( * )
285 COMPLEX B( LDB, * ), BX( LDBX, * )
286* ..
287*
288* =====================================================================
289*
290* .. Parameters ..
291 REAL ONE, ZERO, NEGONE
292 PARAMETER ( ONE = 1.0e0, zero = 0.0e0, negone = -1.0e0 )
293* ..
294* .. Local Scalars ..
295 INTEGER I, J, JCOL, JROW, M, N, NLP1
296 REAL DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
297* ..
298* .. External Subroutines ..
299 EXTERNAL ccopy, clacpy, clascl, csrot, csscal, sgemv,
300 $ xerbla
301* ..
302* .. External Functions ..
303 REAL SLAMC3, SNRM2
304 EXTERNAL SLAMC3, SNRM2
305* ..
306* .. Intrinsic Functions ..
307 INTRINSIC aimag, cmplx, max, real
308* ..
309* .. Executable Statements ..
310*
311* Test the input parameters.
312*
313 info = 0
314 n = nl + nr + 1
315*
316 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
317 info = -1
318 ELSE IF( nl.LT.1 ) THEN
319 info = -2
320 ELSE IF( nr.LT.1 ) THEN
321 info = -3
322 ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
323 info = -4
324 ELSE IF( nrhs.LT.1 ) THEN
325 info = -5
326 ELSE IF( ldb.LT.n ) THEN
327 info = -7
328 ELSE IF( ldbx.LT.n ) THEN
329 info = -9
330 ELSE IF( givptr.LT.0 ) THEN
331 info = -11
332 ELSE IF( ldgcol.LT.n ) THEN
333 info = -13
334 ELSE IF( ldgnum.LT.n ) THEN
335 info = -15
336 ELSE IF( k.LT.1 ) THEN
337 info = -20
338 END IF
339 IF( info.NE.0 ) THEN
340 CALL xerbla( 'CLALS0', -info )
341 RETURN
342 END IF
343*
344 m = n + sqre
345 nlp1 = nl + 1
346*
347 IF( icompq.EQ.0 ) THEN
348*
349* Apply back orthogonal transformations from the left.
350*
351* Step (1L): apply back the Givens rotations performed.
352*
353 DO 10 i = 1, givptr
354 CALL csrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
355 $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
356 $ givnum( i, 1 ) )
357 10 CONTINUE
358*
359* Step (2L): permute rows of B.
360*
361 CALL ccopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
362 DO 20 i = 2, n
363 CALL ccopy( nrhs, b( perm( i ), 1 ), ldb, bx( i, 1 ), ldbx )
364 20 CONTINUE
365*
366* Step (3L): apply the inverse of the left singular vector
367* matrix to BX.
368*
369 IF( k.EQ.1 ) THEN
370 CALL ccopy( nrhs, bx, ldbx, b, ldb )
371 IF( z( 1 ).LT.zero ) THEN
372 CALL csscal( nrhs, negone, b, ldb )
373 END IF
374 ELSE
375 DO 100 j = 1, k
376 diflj = difl( j )
377 dj = poles( j, 1 )
378 dsigj = -poles( j, 2 )
379 IF( j.LT.k ) THEN
380 difrj = -difr( j, 1 )
381 dsigjp = -poles( j+1, 2 )
382 END IF
383 IF( ( z( j ).EQ.zero ) .OR. ( poles( j, 2 ).EQ.zero ) )
384 $ THEN
385 rwork( j ) = zero
386 ELSE
387 rwork( j ) = -poles( j, 2 )*z( j ) / diflj /
388 $ ( poles( j, 2 )+dj )
389 END IF
390 DO 30 i = 1, j - 1
391 IF( ( z( i ).EQ.zero ) .OR.
392 $ ( poles( i, 2 ).EQ.zero ) ) THEN
393 rwork( i ) = zero
394 ELSE
395*
396* Use calls to the subroutine SLAMC3 to enforce the
397* parentheses (x+y)+z. The goal is to prevent
398* optimizing compilers from doing x+(y+z).
399*
400 rwork( i ) = poles( i, 2 )*z( i ) /
401 $ ( slamc3( poles( i, 2 ), dsigj )-
402 $ diflj ) / ( poles( i, 2 )+dj )
403 END IF
404 30 CONTINUE
405 DO 40 i = j + 1, k
406 IF( ( z( i ).EQ.zero ) .OR.
407 $ ( poles( i, 2 ).EQ.zero ) ) THEN
408 rwork( i ) = zero
409 ELSE
410 rwork( i ) = poles( i, 2 )*z( i ) /
411 $ ( slamc3( poles( i, 2 ), dsigjp )+
412 $ difrj ) / ( poles( i, 2 )+dj )
413 END IF
414 40 CONTINUE
415 rwork( 1 ) = negone
416 temp = snrm2( k, rwork, 1 )
417*
418* Since B and BX are complex, the following call to SGEMV
419* is performed in two steps (real and imaginary parts).
420*
421* CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
422* $ B( J, 1 ), LDB )
423*
424 i = k + nrhs*2
425 DO 60 jcol = 1, nrhs
426 DO 50 jrow = 1, k
427 i = i + 1
428 rwork( i ) = real( bx( jrow, jcol ) )
429 50 CONTINUE
430 60 CONTINUE
431 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
432 $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
433 i = k + nrhs*2
434 DO 80 jcol = 1, nrhs
435 DO 70 jrow = 1, k
436 i = i + 1
437 rwork( i ) = aimag( bx( jrow, jcol ) )
438 70 CONTINUE
439 80 CONTINUE
440 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
441 $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
442 DO 90 jcol = 1, nrhs
443 b( j, jcol ) = cmplx( rwork( jcol+k ),
444 $ rwork( jcol+k+nrhs ) )
445 90 CONTINUE
446 CALL clascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
447 $ ldb, info )
448 100 CONTINUE
449 END IF
450*
451* Move the deflated rows of BX to B also.
452*
453 IF( k.LT.max( m, n ) )
454 $ CALL clacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
455 $ b( k+1, 1 ), ldb )
456 ELSE
457*
458* Apply back the right orthogonal transformations.
459*
460* Step (1R): apply back the new right singular vector matrix
461* to B.
462*
463 IF( k.EQ.1 ) THEN
464 CALL ccopy( nrhs, b, ldb, bx, ldbx )
465 ELSE
466 DO 180 j = 1, k
467 dsigj = poles( j, 2 )
468 IF( z( j ).EQ.zero ) THEN
469 rwork( j ) = zero
470 ELSE
471 rwork( j ) = -z( j ) / difl( j ) /
472 $ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
473 END IF
474 DO 110 i = 1, j - 1
475 IF( z( j ).EQ.zero ) THEN
476 rwork( i ) = zero
477 ELSE
478*
479* Use calls to the subroutine SLAMC3 to enforce the
480* parentheses (x+y)+z. The goal is to prevent optimizing
481* compilers from doing x+(y+z).
482*
483 rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i+1,
484 $ 2 ) )-difr( i, 1 ) ) /
485 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
486 END IF
487 110 CONTINUE
488 DO 120 i = j + 1, k
489 IF( z( j ).EQ.zero ) THEN
490 rwork( i ) = zero
491 ELSE
492 rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i,
493 $ 2 ) )-difl( i ) ) /
494 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
495 END IF
496 120 CONTINUE
497*
498* Since B and BX are complex, the following call to SGEMV
499* is performed in two steps (real and imaginary parts).
500*
501* CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
502* $ BX( J, 1 ), LDBX )
503*
504 i = k + nrhs*2
505 DO 140 jcol = 1, nrhs
506 DO 130 jrow = 1, k
507 i = i + 1
508 rwork( i ) = real( b( jrow, jcol ) )
509 130 CONTINUE
510 140 CONTINUE
511 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
512 $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
513 i = k + nrhs*2
514 DO 160 jcol = 1, nrhs
515 DO 150 jrow = 1, k
516 i = i + 1
517 rwork( i ) = aimag( b( jrow, jcol ) )
518 150 CONTINUE
519 160 CONTINUE
520 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
521 $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
522 DO 170 jcol = 1, nrhs
523 bx( j, jcol ) = cmplx( rwork( jcol+k ),
524 $ rwork( jcol+k+nrhs ) )
525 170 CONTINUE
526 180 CONTINUE
527 END IF
528*
529* Step (2R): if SQRE = 1, apply back the rotation that is
530* related to the right null space of the subproblem.
531*
532 IF( sqre.EQ.1 ) THEN
533 CALL ccopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
534 CALL csrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
535 END IF
536 IF( k.LT.max( m, n ) )
537 $ CALL clacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb,
538 $ bx( k+1, 1 ), ldbx )
539*
540* Step (3R): permute rows of B.
541*
542 CALL ccopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
543 IF( sqre.EQ.1 ) THEN
544 CALL ccopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
545 END IF
546 DO 190 i = 2, n
547 CALL ccopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
548 190 CONTINUE
549*
550* Step (4R): apply back the Givens rotations performed.
551*
552 DO 200 i = givptr, 1, -1
553 CALL csrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
554 $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
555 $ -givnum( i, 1 ) )
556 200 CONTINUE
557 END IF
558*
559 RETURN
560*
561* End of CLALS0
562*
563 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
clals0
subroutine clals0(icompq, nl, nr, sqre, nrhs, b, ldb, bx, ldbx, perm, givptr, givcol, ldgcol, givnum, ldgnum, poles, difl, difr, z, k, c, s, rwork, info)
CLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer...
Definition clals0.f:270
clascl
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:143
csrot
subroutine csrot(n, cx, incx, cy, incy, c, s)
CSROT
Definition csrot.f:98
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78