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