LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date September 2012
258 *
259 *> \ingroup complexOTHERcomputational
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 clals0( 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  REAL c, s
282 * ..
283 * .. Array Arguments ..
284  INTEGER givcol( ldgcol, * ), perm( * )
285  REAL difl( * ), difr( ldgnum, * ),
286  $ givnum( ldgnum, * ), poles( ldgnum, * ),
287  $ rwork( * ), z( * )
288  COMPLEX b( ldb, * ), bx( ldbx, * )
289 * ..
290 *
291 * =====================================================================
292 *
293 * .. Parameters ..
294  REAL one, zero, negone
295  parameter( one = 1.0e0, zero = 0.0e0, negone = -1.0e0 )
296 * ..
297 * .. Local Scalars ..
298  INTEGER i, j, jcol, jrow, m, n, nlp1
299  REAL diflj, difrj, dj, dsigj, dsigjp, temp
300 * ..
301 * .. External Subroutines ..
302  EXTERNAL ccopy, clacpy, clascl, csrot, csscal, sgemv,
303  $ xerbla
304 * ..
305 * .. External Functions ..
306  REAL slamc3, snrm2
307  EXTERNAL slamc3, snrm2
308 * ..
309 * .. Intrinsic Functions ..
310  INTRINSIC aimag, cmplx, max, real
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( 'CLALS0', -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 csrot( 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 ccopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
368  DO 20 i = 2, n
369  CALL ccopy( 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 ccopy( nrhs, bx, ldbx, b, ldb )
377  IF( z( 1 ).LT.zero ) THEN
378  CALL csscal( 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  $ ( slamc3( 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  $ ( slamc3( poles( i, 2 ), dsigjp )+
413  $ difrj ) / ( poles( i, 2 )+dj )
414  END IF
415  40 continue
416  rwork( 1 ) = negone
417  temp = snrm2( k, rwork, 1 )
418 *
419 * Since B and BX are complex, the following call to SGEMV
420 * is performed in two steps (real and imaginary parts).
421 *
422 * CALL SGEMV( '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 ) = REAL( BX( JROW, JCOL ) )
430  50 continue
431  60 continue
432  CALL sgemv( '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 ) = aimag( bx( jrow, jcol ) )
439  70 continue
440  80 continue
441  CALL sgemv( '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 ) = cmplx( rwork( jcol+k ),
445  $ rwork( jcol+k+nrhs ) )
446  90 continue
447  CALL clascl( '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 clacpy( '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 ccopy( 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 ) / ( slamc3( 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 ) / ( slamc3( 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 SGEMV
495 * is performed in two steps (real and imaginary parts).
496 *
497 * CALL SGEMV( '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 ) = REAL( B( JROW, JCOL ) )
505  130 continue
506  140 continue
507  CALL sgemv( '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 ) = aimag( b( jrow, jcol ) )
514  150 continue
515  160 continue
516  CALL sgemv( '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 ) = cmplx( 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 ccopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
530  CALL csrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
531  END IF
532  IF( k.LT.max( m, n ) )
533  $ CALL clacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb,
534  $ bx( k+1, 1 ), ldbx )
535 *
536 * Step (3R): permute rows of B.
537 *
538  CALL ccopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
539  IF( sqre.EQ.1 ) THEN
540  CALL ccopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
541  END IF
542  DO 190 i = 2, n
543  CALL ccopy( 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 csrot( 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 CLALS0
558 *
559  END