LAPACK 3.11.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 complexOTHERcomputational
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 rwork( i ) = poles( i, 2 )*z( i ) /
396 $ ( slamc3( poles( i, 2 ), dsigj )-
397 $ diflj ) / ( poles( i, 2 )+dj )
398 END IF
399 30 CONTINUE
400 DO 40 i = j + 1, k
401 IF( ( z( i ).EQ.zero ) .OR.
402 $ ( poles( i, 2 ).EQ.zero ) ) THEN
403 rwork( i ) = zero
404 ELSE
405 rwork( i ) = poles( i, 2 )*z( i ) /
406 $ ( slamc3( poles( i, 2 ), dsigjp )+
407 $ difrj ) / ( poles( i, 2 )+dj )
408 END IF
409 40 CONTINUE
410 rwork( 1 ) = negone
411 temp = snrm2( k, rwork, 1 )
412*
413* Since B and BX are complex, the following call to SGEMV
414* is performed in two steps (real and imaginary parts).
415*
416* CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
417* $ B( J, 1 ), LDB )
418*
419 i = k + nrhs*2
420 DO 60 jcol = 1, nrhs
421 DO 50 jrow = 1, k
422 i = i + 1
423 rwork( i ) = real( bx( jrow, jcol ) )
424 50 CONTINUE
425 60 CONTINUE
426 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
427 $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
428 i = k + nrhs*2
429 DO 80 jcol = 1, nrhs
430 DO 70 jrow = 1, k
431 i = i + 1
432 rwork( i ) = aimag( bx( jrow, jcol ) )
433 70 CONTINUE
434 80 CONTINUE
435 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
436 $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
437 DO 90 jcol = 1, nrhs
438 b( j, jcol ) = cmplx( rwork( jcol+k ),
439 $ rwork( jcol+k+nrhs ) )
440 90 CONTINUE
441 CALL clascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
442 $ ldb, info )
443 100 CONTINUE
444 END IF
445*
446* Move the deflated rows of BX to B also.
447*
448 IF( k.LT.max( m, n ) )
449 $ CALL clacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
450 $ b( k+1, 1 ), ldb )
451 ELSE
452*
453* Apply back the right orthogonal transformations.
454*
455* Step (1R): apply back the new right singular vector matrix
456* to B.
457*
458 IF( k.EQ.1 ) THEN
459 CALL ccopy( nrhs, b, ldb, bx, ldbx )
460 ELSE
461 DO 180 j = 1, k
462 dsigj = poles( j, 2 )
463 IF( z( j ).EQ.zero ) THEN
464 rwork( j ) = zero
465 ELSE
466 rwork( j ) = -z( j ) / difl( j ) /
467 $ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
468 END IF
469 DO 110 i = 1, j - 1
470 IF( z( j ).EQ.zero ) THEN
471 rwork( i ) = zero
472 ELSE
473 rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i+1,
474 $ 2 ) )-difr( i, 1 ) ) /
475 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
476 END IF
477 110 CONTINUE
478 DO 120 i = j + 1, k
479 IF( z( j ).EQ.zero ) THEN
480 rwork( i ) = zero
481 ELSE
482 rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i,
483 $ 2 ) )-difl( i ) ) /
484 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
485 END IF
486 120 CONTINUE
487*
488* Since B and BX are complex, the following call to SGEMV
489* is performed in two steps (real and imaginary parts).
490*
491* CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
492* $ BX( J, 1 ), LDBX )
493*
494 i = k + nrhs*2
495 DO 140 jcol = 1, nrhs
496 DO 130 jrow = 1, k
497 i = i + 1
498 rwork( i ) = real( b( jrow, jcol ) )
499 130 CONTINUE
500 140 CONTINUE
501 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
502 $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
503 i = k + nrhs*2
504 DO 160 jcol = 1, nrhs
505 DO 150 jrow = 1, k
506 i = i + 1
507 rwork( i ) = aimag( b( jrow, jcol ) )
508 150 CONTINUE
509 160 CONTINUE
510 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
511 $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
512 DO 170 jcol = 1, nrhs
513 bx( j, jcol ) = cmplx( rwork( jcol+k ),
514 $ rwork( jcol+k+nrhs ) )
515 170 CONTINUE
516 180 CONTINUE
517 END IF
518*
519* Step (2R): if SQRE = 1, apply back the rotation that is
520* related to the right null space of the subproblem.
521*
522 IF( sqre.EQ.1 ) THEN
523 CALL ccopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
524 CALL csrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
525 END IF
526 IF( k.LT.max( m, n ) )
527 $ CALL clacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb,
528 $ bx( k+1, 1 ), ldbx )
529*
530* Step (3R): permute rows of B.
531*
532 CALL ccopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
533 IF( sqre.EQ.1 ) THEN
534 CALL ccopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
535 END IF
536 DO 190 i = 2, n
537 CALL ccopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
538 190 CONTINUE
539*
540* Step (4R): apply back the Givens rotations performed.
541*
542 DO 200 i = givptr, 1, -1
543 CALL csrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
544 $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
545 $ -givnum( i, 1 ) )
546 200 CONTINUE
547 END IF
548*
549 RETURN
550*
551* End of CLALS0
552*
553 END
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
csrot
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
Definition: csrot.f:98
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
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
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156