LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
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 December 2016
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.7.0) --
274 * -- LAPACK is a software package provided by Univ. of Tennessee, --
275 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276 * December 2016
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  n = nl + nr + 1
318 *
319  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
320  info = -1
321  ELSE IF( nl.LT.1 ) THEN
322  info = -2
323  ELSE IF( nr.LT.1 ) THEN
324  info = -3
325  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
326  info = -4
327  ELSE 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( 'CLALS0', -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 csrot( 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 ccopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
365  DO 20 i = 2, n
366  CALL ccopy( 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 ccopy( nrhs, bx, ldbx, b, ldb )
374  IF( z( 1 ).LT.zero ) THEN
375  CALL csscal( nrhs, negone, b, ldb )
376  END IF
377  ELSE
378  DO 100 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  rwork( j ) = zero
389  ELSE
390  rwork( 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  rwork( i ) = zero
397  ELSE
398  rwork( i ) = poles( i, 2 )*z( i ) /
399  $ ( slamc3( 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  rwork( i ) = zero
407  ELSE
408  rwork( i ) = poles( i, 2 )*z( i ) /
409  $ ( slamc3( poles( i, 2 ), dsigjp )+
410  $ difrj ) / ( poles( i, 2 )+dj )
411  END IF
412  40 CONTINUE
413  rwork( 1 ) = negone
414  temp = snrm2( k, rwork, 1 )
415 *
416 * Since B and BX are complex, the following call to SGEMV
417 * is performed in two steps (real and imaginary parts).
418 *
419 * CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
420 * $ B( J, 1 ), LDB )
421 *
422  i = k + nrhs*2
423  DO 60 jcol = 1, nrhs
424  DO 50 jrow = 1, k
425  i = i + 1
426  rwork( i ) = real( bx( jrow, jcol ) )
427  50 CONTINUE
428  60 CONTINUE
429  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
430  $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
431  i = k + nrhs*2
432  DO 80 jcol = 1, nrhs
433  DO 70 jrow = 1, k
434  i = i + 1
435  rwork( i ) = aimag( bx( jrow, jcol ) )
436  70 CONTINUE
437  80 CONTINUE
438  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
439  $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
440  DO 90 jcol = 1, nrhs
441  b( j, jcol ) = cmplx( rwork( jcol+k ),
442  $ rwork( jcol+k+nrhs ) )
443  90 CONTINUE
444  CALL clascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
445  $ ldb, info )
446  100 CONTINUE
447  END IF
448 *
449 * Move the deflated rows of BX to B also.
450 *
451  IF( k.LT.max( m, n ) )
452  $ CALL clacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
453  $ b( k+1, 1 ), ldb )
454  ELSE
455 *
456 * Apply back the right orthogonal transformations.
457 *
458 * Step (1R): apply back the new right singular vector matrix
459 * to B.
460 *
461  IF( k.EQ.1 ) THEN
462  CALL ccopy( nrhs, b, ldb, bx, ldbx )
463  ELSE
464  DO 180 j = 1, k
465  dsigj = poles( j, 2 )
466  IF( z( j ).EQ.zero ) THEN
467  rwork( j ) = zero
468  ELSE
469  rwork( j ) = -z( j ) / difl( j ) /
470  $ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
471  END IF
472  DO 110 i = 1, j - 1
473  IF( z( j ).EQ.zero ) THEN
474  rwork( i ) = zero
475  ELSE
476  rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i+1,
477  $ 2 ) )-difr( i, 1 ) ) /
478  $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
479  END IF
480  110 CONTINUE
481  DO 120 i = j + 1, k
482  IF( z( j ).EQ.zero ) THEN
483  rwork( i ) = zero
484  ELSE
485  rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i,
486  $ 2 ) )-difl( i ) ) /
487  $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
488  END IF
489  120 CONTINUE
490 *
491 * Since B and BX are complex, the following call to SGEMV
492 * is performed in two steps (real and imaginary parts).
493 *
494 * CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
495 * $ BX( J, 1 ), LDBX )
496 *
497  i = k + nrhs*2
498  DO 140 jcol = 1, nrhs
499  DO 130 jrow = 1, k
500  i = i + 1
501  rwork( i ) = real( b( jrow, jcol ) )
502  130 CONTINUE
503  140 CONTINUE
504  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
505  $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
506  i = k + nrhs*2
507  DO 160 jcol = 1, nrhs
508  DO 150 jrow = 1, k
509  i = i + 1
510  rwork( i ) = aimag( b( jrow, jcol ) )
511  150 CONTINUE
512  160 CONTINUE
513  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
514  $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
515  DO 170 jcol = 1, nrhs
516  bx( j, jcol ) = cmplx( rwork( jcol+k ),
517  $ rwork( jcol+k+nrhs ) )
518  170 CONTINUE
519  180 CONTINUE
520  END IF
521 *
522 * Step (2R): if SQRE = 1, apply back the rotation that is
523 * related to the right null space of the subproblem.
524 *
525  IF( sqre.EQ.1 ) THEN
526  CALL ccopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
527  CALL csrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
528  END IF
529  IF( k.LT.max( m, n ) )
530  $ CALL clacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb,
531  $ bx( k+1, 1 ), ldbx )
532 *
533 * Step (3R): permute rows of B.
534 *
535  CALL ccopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
536  IF( sqre.EQ.1 ) THEN
537  CALL ccopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
538  END IF
539  DO 190 i = 2, n
540  CALL ccopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
541  190 CONTINUE
542 *
543 * Step (4R): apply back the Givens rotations performed.
544 *
545  DO 200 i = givptr, 1, -1
546  CALL csrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
547  $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
548  $ -givnum( i, 1 ) )
549  200 CONTINUE
550  END IF
551 *
552  RETURN
553 *
554 * End of CLALS0
555 *
556  END
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:80
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
csrot
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
Definition: csrot.f:100
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
clals0
subroutine clals0(ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO)
CLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer...
Definition: clals0.f:272
clascl
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:145