LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zlals0.f
Go to the documentation of this file.
1*> \brief \b ZLALS0 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 ZLALS0 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlals0.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlals0.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlals0.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZLALS0( 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* DOUBLE PRECISION C, S
27* ..
28* .. Array Arguments ..
29* INTEGER GIVCOL( LDGCOL, * ), PERM( * )
30* DOUBLE PRECISION DIFL( * ), DIFR( LDGNUM, * ),
31* $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
32* $ RWORK( * ), Z( * )
33* COMPLEX*16 B( LDB, * ), BX( LDBX, * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> ZLALS0 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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 zlals0( 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 DOUBLE PRECISION C, S
278* ..
279* .. Array Arguments ..
280 INTEGER GIVCOL( LDGCOL, * ), PERM( * )
281 DOUBLE PRECISION DIFL( * ), DIFR( LDGNUM, * ),
282 $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
283 $ rwork( * ), z( * )
284 COMPLEX*16 B( LDB, * ), BX( LDBX, * )
285* ..
286*
287* =====================================================================
288*
289* .. Parameters ..
290 DOUBLE PRECISION ONE, ZERO, NEGONE
291 PARAMETER ( ONE = 1.0d0, zero = 0.0d0, negone = -1.0d0 )
292* ..
293* .. Local Scalars ..
294 INTEGER I, J, JCOL, JROW, M, N, NLP1
295 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
296* ..
297* .. External Subroutines ..
298 EXTERNAL dgemv, xerbla, zcopy, zdrot, zdscal,
299 $ zlacpy,
300 $ zlascl
301* ..
302* .. External Functions ..
303 DOUBLE PRECISION DLAMC3, DNRM2
304 EXTERNAL DLAMC3, DNRM2
305* ..
306* .. Intrinsic Functions ..
307 INTRINSIC dble, dcmplx, dimag, max
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( 'ZLALS0', -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 zdrot( 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 zcopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
362 DO 20 i = 2, n
363 CALL zcopy( 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 zcopy( nrhs, bx, ldbx, b, ldb )
372 IF( z( 1 ).LT.zero ) THEN
373 CALL zdscal( 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 DLAMC3 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 $ ( dlamc3( 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 $ ( dlamc3( poles( i, 2 ), dsigjp )+
413 $ difrj ) / ( poles( i, 2 )+dj )
414 END IF
415 40 CONTINUE
416 rwork( 1 ) = negone
417 temp = dnrm2( k, rwork, 1 )
418*
419* Since B and BX are complex, the following call to DGEMV
420* is performed in two steps (real and imaginary parts).
421*
422* CALL DGEMV( '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 ) = dble( bx( jrow, jcol ) )
430 50 CONTINUE
431 60 CONTINUE
432 CALL dgemv( '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 ) = dimag( bx( jrow, jcol ) )
439 70 CONTINUE
440 80 CONTINUE
441 CALL dgemv( '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 ) = dcmplx( rwork( jcol+k ),
445 $ rwork( jcol+k+nrhs ) )
446 90 CONTINUE
447 CALL zlascl( '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 zlacpy( '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 zcopy( 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 DLAMC3 to enforce the
481* parentheses (x+y)+z. The goal is to prevent
482* optimizing compilers from doing x+(y+z).
483*
484 rwork( i ) = z( j ) / ( dlamc3( 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 ) / ( dlamc3( 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 DGEMV
502* is performed in two steps (real and imaginary parts).
503*
504* CALL DGEMV( '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 ) = dble( b( jrow, jcol ) )
512 130 CONTINUE
513 140 CONTINUE
514 CALL dgemv( '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 ) = dimag( b( jrow, jcol ) )
521 150 CONTINUE
522 160 CONTINUE
523 CALL dgemv( '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 ) = dcmplx( 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 zcopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
537 CALL zdrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c,
538 $ s )
539 END IF
540 IF( k.LT.max( m, n ) )
541 $ CALL zlacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb, bx( k+1,
542 $ 1 ),
543 $ ldbx )
544*
545* Step (3R): permute rows of B.
546*
547 CALL zcopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
548 IF( sqre.EQ.1 ) THEN
549 CALL zcopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
550 END IF
551 DO 190 i = 2, n
552 CALL zcopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ),
553 $ ldb )
554 190 CONTINUE
555*
556* Step (4R): apply back the Givens rotations performed.
557*
558 DO 200 i = givptr, 1, -1
559 CALL zdrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
560 $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
561 $ -givnum( i, 1 ) )
562 200 CONTINUE
563 END IF
564*
565 RETURN
566*
567* End of ZLALS0
568*
569 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
zlals0
subroutine zlals0(icompq, nl, nr, sqre, nrhs, b, ldb, bx, ldbx, perm, givptr, givcol, ldgcol, givnum, ldgnum, poles, difl, difr, z, k, c, s, rwork, info)
ZLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer...
Definition zlals0.f:269
zlascl
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition zlascl.f:142
zdrot
subroutine zdrot(n, zx, incx, zy, incy, c, s)
ZDROT
Definition zdrot.f:98
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78