LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
zlals0.f
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 *> \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 zlals0( 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  DOUBLE PRECISION C, S
279 * ..
280 * .. Array Arguments ..
281  INTEGER GIVCOL( ldgcol, * ), PERM( * )
282  DOUBLE PRECISION DIFL( * ), DIFR( ldgnum, * ),
283  $ givnum( ldgnum, * ), poles( ldgnum, * ),
284  $ rwork( * ), z( * )
285  COMPLEX*16 B( ldb, * ), BX( ldbx, * )
286 * ..
287 *
288 * =====================================================================
289 *
290 * .. Parameters ..
291  DOUBLE PRECISION ONE, ZERO, NEGONE
292  parameter( one = 1.0d0, zero = 0.0d0, negone = -1.0d0 )
293 * ..
294 * .. Local Scalars ..
295  INTEGER I, J, JCOL, JROW, M, N, NLP1
296  DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
297 * ..
298 * .. External Subroutines ..
299  EXTERNAL dgemv, xerbla, zcopy, zdrot, zdscal, 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 ), 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 zcopy( nrhs, bx, ldbx, b, ldb )
371  IF( z( 1 ).LT.zero ) THEN
372  CALL zdscal( 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  $ ( dlamc3( 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  $ ( dlamc3( poles( i, 2 ), dsigjp )+
407  $ difrj ) / ( poles( i, 2 )+dj )
408  END IF
409  40 CONTINUE
410  rwork( 1 ) = negone
411  temp = dnrm2( k, rwork, 1 )
412 *
413 * Since B and BX are complex, the following call to DGEMV
414 * is performed in two steps (real and imaginary parts).
415 *
416 * CALL DGEMV( '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 ) = dble( bx( jrow, jcol ) )
424  50 CONTINUE
425  60 CONTINUE
426  CALL dgemv( '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 ) = dimag( bx( jrow, jcol ) )
433  70 CONTINUE
434  80 CONTINUE
435  CALL dgemv( '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 ) = dcmplx( rwork( jcol+k ),
439  $ rwork( jcol+k+nrhs ) )
440  90 CONTINUE
441  CALL zlascl( '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 zlacpy( '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 zcopy( 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 ) / ( dlamc3( 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 ) / ( dlamc3( 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 DGEMV
489 * is performed in two steps (real and imaginary parts).
490 *
491 * CALL DGEMV( '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 ) = dble( b( jrow, jcol ) )
499  130 CONTINUE
500  140 CONTINUE
501  CALL dgemv( '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 ) = dimag( b( jrow, jcol ) )
508  150 CONTINUE
509  160 CONTINUE
510  CALL dgemv( '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 ) = dcmplx( 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 zcopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
524  CALL zdrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
525  END IF
526  IF( k.LT.max( m, n ) )
527  $ CALL zlacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb, bx( k+1, 1 ),
528  $ ldbx )
529 *
530 * Step (3R): permute rows of B.
531 *
532  CALL zcopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
533  IF( sqre.EQ.1 ) THEN
534  CALL zcopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
535  END IF
536  DO 190 i = 2, n
537  CALL zcopy( 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 zdrot( 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 ZLALS0
552 *
553  END
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:270
zdscal
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
zdrot
subroutine zdrot(N, ZX, INCX, ZY, INCY, C, S)
ZDROT
Definition: zdrot.f:98
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
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:143
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:103