LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Files Functions Typedefs Macros
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 *> \htmlonly
9 *> Download ZLALS0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlals0.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlals0.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlals0.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLALS0( 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 * DOUBLE PRECISION C, S
29 * ..
30 * .. Array Arguments ..
31 * INTEGER GIVCOL( LDGCOL, * ), PERM( * )
32 * DOUBLE PRECISION DIFL( * ), DIFR( LDGNUM, * ),
33 * $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
34 * $ RWORK( * ), Z( * )
35 * COMPLEX*16 B( LDB, * ), BX( LDBX, * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> ZLALS0 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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 *> \date September 2012
258 *
259 *> \ingroup complex16OTHERcomputational
260 *
261 *> \par Contributors:
262 * ==================
263 *>
264 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
265 *> California at Berkeley, USA \n
266 *> Osni Marques, LBNL/NERSC, USA \n
267 *
268 * =====================================================================
269  SUBROUTINE zlals0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
270  $ perm, givptr, givcol, ldgcol, givnum, ldgnum,
271  $ poles, difl, difr, z, k, c, s, rwork, info )
272 *
273 * -- LAPACK computational routine (version 3.4.2) --
274 * -- LAPACK is a software package provided by Univ. of Tennessee, --
275 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276 * September 2012
277 *
278 * .. Scalar Arguments ..
279  INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
280  $ ldgnum, nl, nr, nrhs, sqre
281  DOUBLE PRECISION C, S
282 * ..
283 * .. Array Arguments ..
284  INTEGER GIVCOL( ldgcol, * ), PERM( * )
285  DOUBLE PRECISION DIFL( * ), DIFR( ldgnum, * ),
286  $ givnum( ldgnum, * ), poles( ldgnum, * ),
287  $ rwork( * ), z( * )
288  COMPLEX*16 B( ldb, * ), BX( ldbx, * )
289 * ..
290 *
291 * =====================================================================
292 *
293 * .. Parameters ..
294  DOUBLE PRECISION ONE, ZERO, NEGONE
295  parameter( one = 1.0d0, zero = 0.0d0, negone = -1.0d0 )
296 * ..
297 * .. Local Scalars ..
298  INTEGER I, J, JCOL, JROW, M, N, NLP1
299  DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
300 * ..
301 * .. External Subroutines ..
302  EXTERNAL dgemv, xerbla, zcopy, zdrot, zdscal, zlacpy,
303  $ zlascl
304 * ..
305 * .. External Functions ..
306  DOUBLE PRECISION DLAMC3, DNRM2
307  EXTERNAL dlamc3, dnrm2
308 * ..
309 * .. Intrinsic Functions ..
310  INTRINSIC dble, dcmplx, dimag, max
311 * ..
312 * .. Executable Statements ..
313 *
314 * Test the input parameters.
315 *
316  info = 0
317 *
318  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
319  info = -1
320  ELSE IF( nl.LT.1 ) THEN
321  info = -2
322  ELSE IF( nr.LT.1 ) THEN
323  info = -3
324  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
325  info = -4
326  END IF
327 *
328  n = nl + nr + 1
329 *
330  IF( nrhs.LT.1 ) THEN
331  info = -5
332  ELSE IF( ldb.LT.n ) THEN
333  info = -7
334  ELSE IF( ldbx.LT.n ) THEN
335  info = -9
336  ELSE IF( givptr.LT.0 ) THEN
337  info = -11
338  ELSE IF( ldgcol.LT.n ) THEN
339  info = -13
340  ELSE IF( ldgnum.LT.n ) THEN
341  info = -15
342  ELSE IF( k.LT.1 ) THEN
343  info = -20
344  END IF
345  IF( info.NE.0 ) THEN
346  CALL xerbla( 'ZLALS0', -info )
347  RETURN
348  END IF
349 *
350  m = n + sqre
351  nlp1 = nl + 1
352 *
353  IF( icompq.EQ.0 ) THEN
354 *
355 * Apply back orthogonal transformations from the left.
356 *
357 * Step (1L): apply back the Givens rotations performed.
358 *
359  DO 10 i = 1, givptr
360  CALL zdrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
361  $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
362  $ givnum( i, 1 ) )
363  10 CONTINUE
364 *
365 * Step (2L): permute rows of B.
366 *
367  CALL zcopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
368  DO 20 i = 2, n
369  CALL zcopy( nrhs, b( perm( i ), 1 ), ldb, bx( i, 1 ), ldbx )
370  20 CONTINUE
371 *
372 * Step (3L): apply the inverse of the left singular vector
373 * matrix to BX.
374 *
375  IF( k.EQ.1 ) THEN
376  CALL zcopy( nrhs, bx, ldbx, b, ldb )
377  IF( z( 1 ).LT.zero ) THEN
378  CALL zdscal( nrhs, negone, b, ldb )
379  END IF
380  ELSE
381  DO 100 j = 1, k
382  diflj = difl( j )
383  dj = poles( j, 1 )
384  dsigj = -poles( j, 2 )
385  IF( j.LT.k ) THEN
386  difrj = -difr( j, 1 )
387  dsigjp = -poles( j+1, 2 )
388  END IF
389  IF( ( z( j ).EQ.zero ) .OR. ( poles( j, 2 ).EQ.zero ) )
390  $ THEN
391  rwork( j ) = zero
392  ELSE
393  rwork( j ) = -poles( j, 2 )*z( j ) / diflj /
394  $ ( poles( j, 2 )+dj )
395  END IF
396  DO 30 i = 1, j - 1
397  IF( ( z( i ).EQ.zero ) .OR.
398  $ ( poles( i, 2 ).EQ.zero ) ) THEN
399  rwork( i ) = zero
400  ELSE
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  rwork( i ) = z( j ) / ( dlamc3( dsigj, -poles( i+1,
480  $ 2 ) )-difr( i, 1 ) ) /
481  $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
482  END IF
483  110 CONTINUE
484  DO 120 i = j + 1, k
485  IF( z( j ).EQ.zero ) THEN
486  rwork( i ) = zero
487  ELSE
488  rwork( i ) = z( j ) / ( dlamc3( dsigj, -poles( i,
489  $ 2 ) )-difl( i ) ) /
490  $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
491  END IF
492  120 CONTINUE
493 *
494 * Since B and BX are complex, the following call to DGEMV
495 * is performed in two steps (real and imaginary parts).
496 *
497 * CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
498 * $ BX( J, 1 ), LDBX )
499 *
500  i = k + nrhs*2
501  DO 140 jcol = 1, nrhs
502  DO 130 jrow = 1, k
503  i = i + 1
504  rwork( i ) = dble( b( jrow, jcol ) )
505  130 CONTINUE
506  140 CONTINUE
507  CALL dgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
508  $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
509  i = k + nrhs*2
510  DO 160 jcol = 1, nrhs
511  DO 150 jrow = 1, k
512  i = i + 1
513  rwork( i ) = dimag( b( jrow, jcol ) )
514  150 CONTINUE
515  160 CONTINUE
516  CALL dgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
517  $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
518  DO 170 jcol = 1, nrhs
519  bx( j, jcol ) = dcmplx( rwork( jcol+k ),
520  $ rwork( jcol+k+nrhs ) )
521  170 CONTINUE
522  180 CONTINUE
523  END IF
524 *
525 * Step (2R): if SQRE = 1, apply back the rotation that is
526 * related to the right null space of the subproblem.
527 *
528  IF( sqre.EQ.1 ) THEN
529  CALL zcopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
530  CALL zdrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
531  END IF
532  IF( k.LT.max( m, n ) )
533  $ CALL zlacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb, bx( k+1, 1 ),
534  $ ldbx )
535 *
536 * Step (3R): permute rows of B.
537 *
538  CALL zcopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
539  IF( sqre.EQ.1 ) THEN
540  CALL zcopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
541  END IF
542  DO 190 i = 2, n
543  CALL zcopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
544  190 CONTINUE
545 *
546 * Step (4R): apply back the Givens rotations performed.
547 *
548  DO 200 i = givptr, 1, -1
549  CALL zdrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
550  $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
551  $ -givnum( i, 1 ) )
552  200 CONTINUE
553  END IF
554 *
555  RETURN
556 *
557 * End of ZLALS0
558 *
559  END
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:104
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:53
subroutine zdrot(N, CX, INCX, CY, INCY, C, S)
ZDROT
Definition: zdrot.f:99
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
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:140
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:51
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
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:157