Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / slals0.f
1 *> \brief \b SLALS0 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 SLALS0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slals0.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slals0.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slals0.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE SLALS0( 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 *       REAL               C, S
29 *       ..
30 *       .. Array Arguments ..
31 *       INTEGER            GIVCOL( LDGCOL, * ), PERM( * )
32 *       REAL               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 *> SLALS0 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL
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 REAL
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 REAL 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 November 2015
256 *
257 *> \ingroup realOTHERcomputational
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 SLALS0( 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.6.0) --
272 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
273 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
274 *     November 2015
275 *
276 *     .. Scalar Arguments ..
277       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
278      $                   LDGNUM, NL, NR, NRHS, SQRE
279       REAL               C, S
280 *     ..
281 *     .. Array Arguments ..
282       INTEGER            GIVCOL( LDGCOL, * ), PERM( * )
283       REAL               B( LDB, * ), BX( LDBX, * ), DIFL( * ),
284      $                   DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
285      $                   POLES( LDGNUM, * ), WORK( * ), Z( * )
286 *     ..
287 *
288 *  =====================================================================
289 *
290 *     .. Parameters ..
291       REAL               ONE, ZERO, NEGONE
292       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0, NEGONE = -1.0E0 )
293 *     ..
294 *     .. Local Scalars ..
295       INTEGER            I, J, M, N, NLP1
296       REAL               DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
297 *     ..
298 *     .. External Subroutines ..
299       EXTERNAL           SCOPY, SGEMV, SLACPY, SLASCL, SROT, SSCAL,
300      $                   XERBLA
301 *     ..
302 *     .. External Functions ..
303       REAL               SLAMC3, SNRM2
304       EXTERNAL           SLAMC3, SNRM2
305 *     ..
306 *     .. Intrinsic Functions ..
307       INTRINSIC          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( 'SLALS0', -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 SROT( 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 SCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
362          DO 20 I = 2, N
363             CALL SCOPY( 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 SCOPY( NRHS, BX, LDBX, B, LDB )
371             IF( Z( 1 ).LT.ZERO ) THEN
372                CALL SSCAL( NRHS, NEGONE, B, LDB )
373             END IF
374          ELSE
375             DO 50 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                   WORK( J ) = ZERO
386                ELSE
387                   WORK( 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                      WORK( I ) = ZERO
394                   ELSE
395                      WORK( I ) = POLES( I, 2 )*Z( I ) /
396      $                           ( SLAMC3( 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                      WORK( I ) = ZERO
404                   ELSE
405                      WORK( I ) = POLES( I, 2 )*Z( I ) /
406      $                           ( SLAMC3( POLES( I, 2 ), DSIGJP )+
407      $                           DIFRJ ) / ( POLES( I, 2 )+DJ )
408                   END IF
409    40          CONTINUE
410                WORK( 1 ) = NEGONE
411                TEMP = SNRM2( K, WORK, 1 )
412                CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
413      $                     B( J, 1 ), LDB )
414                CALL SLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
415      $                      LDB, INFO )
416    50       CONTINUE
417          END IF
418 *
419 *        Move the deflated rows of BX to B also.
420 *
421          IF( K.LT.MAX( M, N ) )
422      $      CALL SLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
423      $                   B( K+1, 1 ), LDB )
424       ELSE
425 *
426 *        Apply back the right orthogonal transformations.
427 *
428 *        Step (1R): apply back the new right singular vector matrix
429 *        to B.
430 *
431          IF( K.EQ.1 ) THEN
432             CALL SCOPY( NRHS, B, LDB, BX, LDBX )
433          ELSE
434             DO 80 J = 1, K
435                DSIGJ = POLES( J, 2 )
436                IF( Z( J ).EQ.ZERO ) THEN
437                   WORK( J ) = ZERO
438                ELSE
439                   WORK( J ) = -Z( J ) / DIFL( J ) /
440      $                        ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
441                END IF
442                DO 60 I = 1, J - 1
443                   IF( Z( J ).EQ.ZERO ) THEN
444                      WORK( I ) = ZERO
445                   ELSE
446                      WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I+1,
447      $                           2 ) )-DIFR( I, 1 ) ) /
448      $                           ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
449                   END IF
450    60          CONTINUE
451                DO 70 I = J + 1, K
452                   IF( Z( J ).EQ.ZERO ) THEN
453                      WORK( I ) = ZERO
454                   ELSE
455                      WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I,
456      $                           2 ) )-DIFL( I ) ) /
457      $                           ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
458                   END IF
459    70          CONTINUE
460                CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
461      $                     BX( J, 1 ), LDBX )
462    80       CONTINUE
463          END IF
464 *
465 *        Step (2R): if SQRE = 1, apply back the rotation that is
466 *        related to the right null space of the subproblem.
467 *
468          IF( SQRE.EQ.1 ) THEN
469             CALL SCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
470             CALL SROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
471          END IF
472          IF( K.LT.MAX( M, N ) )
473      $      CALL SLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
474      $                   LDBX )
475 *
476 *        Step (3R): permute rows of B.
477 *
478          CALL SCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
479          IF( SQRE.EQ.1 ) THEN
480             CALL SCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
481          END IF
482          DO 90 I = 2, N
483             CALL SCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
484    90    CONTINUE
485 *
486 *        Step (4R): apply back the Givens rotations performed.
487 *
488          DO 100 I = GIVPTR, 1, -1
489             CALL SROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
490      $                 B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
491      $                 -GIVNUM( I, 1 ) )
492   100    CONTINUE
493       END IF
494 *
495       RETURN
496 *
497 *     End of SLALS0
498 *
499       END