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 December 2016
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.7.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*     December 2016
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      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( 'DLALS0', -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 DROT( 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 DCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
362         DO 20 I = 2, N
363            CALL DCOPY( 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 DCOPY( NRHS, BX, LDBX, B, LDB )
371            IF( Z( 1 ).LT.ZERO ) THEN
372               CALL DSCAL( 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     $                           ( 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                     WORK( I ) = ZERO
404                  ELSE
405                     WORK( 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               WORK( 1 ) = NEGONE
411               TEMP = DNRM2( K, WORK, 1 )
412               CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
413     $                     B( J, 1 ), LDB )
414               CALL DLASCL( '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 DLACPY( '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 DCOPY( 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 ) / ( DLAMC3( 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 ) / ( DLAMC3( DSIGJ, -POLES( I,
456     $                           2 ) )-DIFL( I ) ) /
457     $                           ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
458                  END IF
459   70          CONTINUE
460               CALL DGEMV( '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 DCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
470            CALL DROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
471         END IF
472         IF( K.LT.MAX( M, N ) )
473     $      CALL DLACPY( '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 DCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
479         IF( SQRE.EQ.1 ) THEN
480            CALL DCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
481         END IF
482         DO 90 I = 2, N
483            CALL DCOPY( 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 DROT( 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 DLALS0
498*
499      END
500