LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlals0.f
Go to the documentation of this file.
1 *> \brief \b DLALS0 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 DLALS0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlals0.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlals0.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlals0.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
22 * PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
23 * POLES, DIFL, DIFR, Z, K, C, S, WORK, 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 B( LDB, * ), BX( LDBX, * ), DIFL( * ),
33 * $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
34 * $ POLES( LDGNUM, * ), WORK( * ), Z( * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> DLALS0 applies back the multiplying factors of either the left or the
44 *> right singular vector matrix of a diagonal matrix appended by a row
45 *> to the right hand side matrix B in solving the least squares problem
46 *> using the divide-and-conquer SVD approach.
47 *>
48 *> For the left singular vector matrix, three types of orthogonal
49 *> matrices are involved:
50 *>
51 *> (1L) Givens rotations: the number of such rotations is GIVPTR; the
52 *> pairs of columns/rows they were applied to are stored in GIVCOL;
53 *> and the C- and S-values of these rotations are stored in GIVNUM.
54 *>
55 *> (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
56 *> row, and for J=2:N, PERM(J)-th row of B is to be moved to the
57 *> J-th row.
58 *>
59 *> (3L) The left singular vector matrix of the remaining matrix.
60 *>
61 *> For the right singular vector matrix, four types of orthogonal
62 *> matrices are involved:
63 *>
64 *> (1R) The right singular vector matrix of the remaining matrix.
65 *>
66 *> (2R) If SQRE = 1, one extra Givens rotation to generate the right
67 *> null space.
68 *>
69 *> (3R) The inverse transformation of (2L).
70 *>
71 *> (4R) The inverse transformation of (1L).
72 *> \endverbatim
73 *
74 * Arguments:
75 * ==========
76 *
77 *> \param[in] ICOMPQ
78 *> \verbatim
79 *> ICOMPQ is INTEGER
80 *> Specifies whether singular vectors are to be computed in
81 *> factored form:
82 *> = 0: Left singular vector matrix.
83 *> = 1: Right singular vector matrix.
84 *> \endverbatim
85 *>
86 *> \param[in] NL
87 *> \verbatim
88 *> NL is INTEGER
89 *> The row dimension of the upper block. NL >= 1.
90 *> \endverbatim
91 *>
92 *> \param[in] NR
93 *> \verbatim
94 *> NR is INTEGER
95 *> The row dimension of the lower block. NR >= 1.
96 *> \endverbatim
97 *>
98 *> \param[in] SQRE
99 *> \verbatim
100 *> SQRE is INTEGER
101 *> = 0: the lower block is an NR-by-NR square matrix.
102 *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
103 *>
104 *> The bidiagonal matrix has row dimension N = NL + NR + 1,
105 *> and column dimension M = N + SQRE.
106 *> \endverbatim
107 *>
108 *> \param[in] NRHS
109 *> \verbatim
110 *> NRHS is INTEGER
111 *> The number of columns of B and BX. NRHS must be at least 1.
112 *> \endverbatim
113 *>
114 *> \param[in,out] B
115 *> \verbatim
116 *> B is DOUBLE PRECISION array, dimension ( LDB, NRHS )
117 *> On input, B contains the right hand sides of the least
118 *> squares problem in rows 1 through M. On output, B contains
119 *> the solution X in rows 1 through N.
120 *> \endverbatim
121 *>
122 *> \param[in] LDB
123 *> \verbatim
124 *> LDB is INTEGER
125 *> The leading dimension of B. LDB must be at least
126 *> max(1,MAX( M, N ) ).
127 *> \endverbatim
128 *>
129 *> \param[out] BX
130 *> \verbatim
131 *> BX is DOUBLE PRECISION array, dimension ( LDBX, NRHS )
132 *> \endverbatim
133 *>
134 *> \param[in] LDBX
135 *> \verbatim
136 *> LDBX is INTEGER
137 *> The leading dimension of BX.
138 *> \endverbatim
139 *>
140 *> \param[in] PERM
141 *> \verbatim
142 *> PERM is INTEGER array, dimension ( N )
143 *> The permutations (from deflation and sorting) applied
144 *> to the two blocks.
145 *> \endverbatim
146 *>
147 *> \param[in] GIVPTR
148 *> \verbatim
149 *> GIVPTR is INTEGER
150 *> The number of Givens rotations which took place in this
151 *> subproblem.
152 *> \endverbatim
153 *>
154 *> \param[in] GIVCOL
155 *> \verbatim
156 *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
157 *> Each pair of numbers indicates a pair of rows/columns
158 *> involved in a Givens rotation.
159 *> \endverbatim
160 *>
161 *> \param[in] LDGCOL
162 *> \verbatim
163 *> LDGCOL is INTEGER
164 *> The leading dimension of GIVCOL, must be at least N.
165 *> \endverbatim
166 *>
167 *> \param[in] GIVNUM
168 *> \verbatim
169 *> GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
170 *> Each number indicates the C or S value used in the
171 *> corresponding Givens rotation.
172 *> \endverbatim
173 *>
174 *> \param[in] LDGNUM
175 *> \verbatim
176 *> LDGNUM is INTEGER
177 *> The leading dimension of arrays DIFR, POLES and
178 *> GIVNUM, must be at least K.
179 *> \endverbatim
180 *>
181 *> \param[in] POLES
182 *> \verbatim
183 *> POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
184 *> On entry, POLES(1:K, 1) contains the new singular
185 *> values obtained from solving the secular equation, and
186 *> POLES(1:K, 2) is an array containing the poles in the secular
187 *> equation.
188 *> \endverbatim
189 *>
190 *> \param[in] DIFL
191 *> \verbatim
192 *> DIFL is DOUBLE PRECISION array, dimension ( K ).
193 *> On entry, DIFL(I) is the distance between I-th updated
194 *> (undeflated) singular value and the I-th (undeflated) old
195 *> singular value.
196 *> \endverbatim
197 *>
198 *> \param[in] DIFR
199 *> \verbatim
200 *> DIFR is DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).
201 *> On entry, DIFR(I, 1) contains the distances between I-th
202 *> updated (undeflated) singular value and the I+1-th
203 *> (undeflated) old singular value. And DIFR(I, 2) is the
204 *> normalizing factor for the I-th right singular vector.
205 *> \endverbatim
206 *>
207 *> \param[in] Z
208 *> \verbatim
209 *> Z is DOUBLE PRECISION array, dimension ( K )
210 *> Contain the components of the deflation-adjusted updating row
211 *> vector.
212 *> \endverbatim
213 *>
214 *> \param[in] K
215 *> \verbatim
216 *> K is INTEGER
217 *> Contains the dimension of the non-deflated matrix,
218 *> This is the order of the related secular equation. 1 <= K <=N.
219 *> \endverbatim
220 *>
221 *> \param[in] C
222 *> \verbatim
223 *> C is DOUBLE PRECISION
224 *> C contains garbage if SQRE =0 and the C-value of a Givens
225 *> rotation related to the right null space if SQRE = 1.
226 *> \endverbatim
227 *>
228 *> \param[in] S
229 *> \verbatim
230 *> S is DOUBLE PRECISION
231 *> S contains garbage if SQRE =0 and the S-value of a Givens
232 *> rotation related to the right null space if SQRE = 1.
233 *> \endverbatim
234 *>
235 *> \param[out] WORK
236 *> \verbatim
237 *> WORK is DOUBLE PRECISION array, dimension ( K )
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 *> \date September 2012
256 *
257 *> \ingroup doubleOTHERcomputational
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 dlals0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
268  $ perm, givptr, givcol, ldgcol, givnum, ldgnum,
269  $ poles, difl, difr, z, k, c, s, work, info )
270 *
271 * -- LAPACK computational routine (version 3.4.2) --
272 * -- LAPACK is a software package provided by Univ. of Tennessee, --
273 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
274 * September 2012
275 *
276 * .. Scalar Arguments ..
277  INTEGER givptr, icompq, info, k, ldb, ldbx, ldgcol,
278  $ ldgnum, nl, nr, nrhs, sqre
279  DOUBLE PRECISION c, s
280 * ..
281 * .. Array Arguments ..
282  INTEGER givcol( ldgcol, * ), perm( * )
283  DOUBLE PRECISION b( ldb, * ), bx( ldbx, * ), difl( * ),
284  $ difr( ldgnum, * ), givnum( ldgnum, * ),
285  $ poles( ldgnum, * ), work( * ), z( * )
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, m, n, nlp1
296  DOUBLE PRECISION diflj, difrj, dj, dsigj, dsigjp, temp
297 * ..
298 * .. External Subroutines ..
299  EXTERNAL dcopy, dgemv, dlacpy, dlascl, drot, dscal,
300  $ xerbla
301 * ..
302 * .. External Functions ..
303  DOUBLE PRECISION dlamc3, dnrm2
304  EXTERNAL dlamc3, dnrm2
305 * ..
306 * .. Intrinsic Functions ..
307  INTRINSIC max
308 * ..
309 * .. Executable Statements ..
310 *
311 * Test the input parameters.
312 *
313  info = 0
314 *
315  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
316  info = -1
317  ELSE IF( nl.LT.1 ) THEN
318  info = -2
319  ELSE IF( nr.LT.1 ) THEN
320  info = -3
321  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
322  info = -4
323  END IF
324 *
325  n = nl + nr + 1
326 *
327  IF( nrhs.LT.1 ) THEN
328  info = -5
329  ELSE IF( ldb.LT.n ) THEN
330  info = -7
331  ELSE IF( ldbx.LT.n ) THEN
332  info = -9
333  ELSE IF( givptr.LT.0 ) THEN
334  info = -11
335  ELSE IF( ldgcol.LT.n ) THEN
336  info = -13
337  ELSE IF( ldgnum.LT.n ) THEN
338  info = -15
339  ELSE IF( k.LT.1 ) THEN
340  info = -20
341  END IF
342  IF( info.NE.0 ) THEN
343  CALL xerbla( 'DLALS0', -info )
344  return
345  END IF
346 *
347  m = n + sqre
348  nlp1 = nl + 1
349 *
350  IF( icompq.EQ.0 ) THEN
351 *
352 * Apply back orthogonal transformations from the left.
353 *
354 * Step (1L): apply back the Givens rotations performed.
355 *
356  DO 10 i = 1, givptr
357  CALL drot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
358  $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
359  $ givnum( i, 1 ) )
360  10 continue
361 *
362 * Step (2L): permute rows of B.
363 *
364  CALL dcopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
365  DO 20 i = 2, n
366  CALL dcopy( nrhs, b( perm( i ), 1 ), ldb, bx( i, 1 ), ldbx )
367  20 continue
368 *
369 * Step (3L): apply the inverse of the left singular vector
370 * matrix to BX.
371 *
372  IF( k.EQ.1 ) THEN
373  CALL dcopy( nrhs, bx, ldbx, b, ldb )
374  IF( z( 1 ).LT.zero ) THEN
375  CALL dscal( nrhs, negone, b, ldb )
376  END IF
377  ELSE
378  DO 50 j = 1, k
379  diflj = difl( j )
380  dj = poles( j, 1 )
381  dsigj = -poles( j, 2 )
382  IF( j.LT.k ) THEN
383  difrj = -difr( j, 1 )
384  dsigjp = -poles( j+1, 2 )
385  END IF
386  IF( ( z( j ).EQ.zero ) .OR. ( poles( j, 2 ).EQ.zero ) )
387  $ THEN
388  work( j ) = zero
389  ELSE
390  work( j ) = -poles( j, 2 )*z( j ) / diflj /
391  $ ( poles( j, 2 )+dj )
392  END IF
393  DO 30 i = 1, j - 1
394  IF( ( z( i ).EQ.zero ) .OR.
395  $ ( poles( i, 2 ).EQ.zero ) ) THEN
396  work( i ) = zero
397  ELSE
398  work( i ) = poles( i, 2 )*z( i ) /
399  $ ( dlamc3( poles( i, 2 ), dsigj )-
400  $ diflj ) / ( poles( i, 2 )+dj )
401  END IF
402  30 continue
403  DO 40 i = j + 1, k
404  IF( ( z( i ).EQ.zero ) .OR.
405  $ ( poles( i, 2 ).EQ.zero ) ) THEN
406  work( i ) = zero
407  ELSE
408  work( i ) = poles( i, 2 )*z( i ) /
409  $ ( dlamc3( poles( i, 2 ), dsigjp )+
410  $ difrj ) / ( poles( i, 2 )+dj )
411  END IF
412  40 continue
413  work( 1 ) = negone
414  temp = dnrm2( k, work, 1 )
415  CALL dgemv( 'T', k, nrhs, one, bx, ldbx, work, 1, zero,
416  $ b( j, 1 ), ldb )
417  CALL dlascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
418  $ ldb, info )
419  50 continue
420  END IF
421 *
422 * Move the deflated rows of BX to B also.
423 *
424  IF( k.LT.max( m, n ) )
425  $ CALL dlacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
426  $ b( k+1, 1 ), ldb )
427  ELSE
428 *
429 * Apply back the right orthogonal transformations.
430 *
431 * Step (1R): apply back the new right singular vector matrix
432 * to B.
433 *
434  IF( k.EQ.1 ) THEN
435  CALL dcopy( nrhs, b, ldb, bx, ldbx )
436  ELSE
437  DO 80 j = 1, k
438  dsigj = poles( j, 2 )
439  IF( z( j ).EQ.zero ) THEN
440  work( j ) = zero
441  ELSE
442  work( j ) = -z( j ) / difl( j ) /
443  $ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
444  END IF
445  DO 60 i = 1, j - 1
446  IF( z( j ).EQ.zero ) THEN
447  work( i ) = zero
448  ELSE
449  work( i ) = z( j ) / ( dlamc3( dsigj, -poles( i+1,
450  $ 2 ) )-difr( i, 1 ) ) /
451  $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
452  END IF
453  60 continue
454  DO 70 i = j + 1, k
455  IF( z( j ).EQ.zero ) THEN
456  work( i ) = zero
457  ELSE
458  work( i ) = z( j ) / ( dlamc3( dsigj, -poles( i,
459  $ 2 ) )-difl( i ) ) /
460  $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
461  END IF
462  70 continue
463  CALL dgemv( 'T', k, nrhs, one, b, ldb, work, 1, zero,
464  $ bx( j, 1 ), ldbx )
465  80 continue
466  END IF
467 *
468 * Step (2R): if SQRE = 1, apply back the rotation that is
469 * related to the right null space of the subproblem.
470 *
471  IF( sqre.EQ.1 ) THEN
472  CALL dcopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
473  CALL drot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
474  END IF
475  IF( k.LT.max( m, n ) )
476  $ CALL dlacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb, bx( k+1, 1 ),
477  $ ldbx )
478 *
479 * Step (3R): permute rows of B.
480 *
481  CALL dcopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
482  IF( sqre.EQ.1 ) THEN
483  CALL dcopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
484  END IF
485  DO 90 i = 2, n
486  CALL dcopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
487  90 continue
488 *
489 * Step (4R): apply back the Givens rotations performed.
490 *
491  DO 100 i = givptr, 1, -1
492  CALL drot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
493  $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
494  $ -givnum( i, 1 ) )
495  100 continue
496  END IF
497 *
498  return
499 *
500 * End of DLALS0
501 *
502  END