1      SUBROUTINE SLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
2     $                   LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
3     $                   GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
4     $                   IWORK, INFO )
5*
6*  -- LAPACK routine (version 3.0) --
7*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
8*     Courant Institute, Argonne National Lab, and Rice University
9*     June 30, 1999
10*
11*     .. Scalar Arguments ..
12      INTEGER            ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
13     $                   SMLSIZ
14*     ..
15*     .. Array Arguments ..
16      INTEGER            GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
17     $                   K( * ), PERM( LDGCOL, * )
18      REAL               B( LDB, * ), BX( LDBX, * ), C( * ),
19     $                   DIFL( LDU, * ), DIFR( LDU, * ),
20     $                   GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
21     $                   U( LDU, * ), VT( LDU, * ), WORK( * ),
22     $                   Z( LDU, * )
23*     ..
24*
25*  Purpose
26*  =======
27*
28*  SLALSA is an itermediate step in solving the least squares problem
29*  by computing the SVD of the coefficient matrix in compact form (The
30*  singular vectors are computed as products of simple orthorgonal
31*  matrices.).
32*
33*  If ICOMPQ = 0, SLALSA applies the inverse of the left singular vector
34*  matrix of an upper bidiagonal matrix to the right hand side; and if
35*  ICOMPQ = 1, SLALSA applies the right singular vector matrix to the
36*  right hand side. The singular vector matrices were generated in
37*  compact form by SLALSA.
38*
39*  Arguments
40*  =========
41*
42*
43*  ICOMPQ (input) INTEGER
44*         Specifies whether the left or the right singular vector
45*         matrix is involved.
46*         = 0: Left singular vector matrix
47*         = 1: Right singular vector matrix
48*
49*  SMLSIZ (input) INTEGER
50*         The maximum size of the subproblems at the bottom of the
51*         computation tree.
52*
53*  N      (input) INTEGER
54*         The row and column dimensions of the upper bidiagonal matrix.
55*
56*  NRHS   (input) INTEGER
57*         The number of columns of B and BX. NRHS must be at least 1.
58*
59*  B      (input) REAL array, dimension ( LDB, NRHS )
60*         On input, B contains the right hand sides of the least
61*         squares problem in rows 1 through M. On output, B contains
62*         the solution X in rows 1 through N.
63*
64*  LDB    (input) INTEGER
65*         The leading dimension of B in the calling subprogram.
66*         LDB must be at least max(1,MAX( M, N ) ).
67*
68*  BX     (output) REAL array, dimension ( LDBX, NRHS )
69*         On exit, the result of applying the left or right singular
70*         vector matrix to B.
71*
72*  LDBX   (input) INTEGER
73*         The leading dimension of BX.
74*
75*  U      (input) REAL array, dimension ( LDU, SMLSIZ ).
76*         On entry, U contains the left singular vector matrices of all
77*         subproblems at the bottom level.
78*
79*  LDU    (input) INTEGER, LDU = > N.
80*         The leading dimension of arrays U, VT, DIFL, DIFR,
81*         POLES, GIVNUM, and Z.
82*
83*  VT     (input) REAL array, dimension ( LDU, SMLSIZ+1 ).
84*         On entry, VT' contains the right singular vector matrices of
85*         all subproblems at the bottom level.
86*
87*  K      (input) INTEGER array, dimension ( N ).
88*
89*  DIFL   (input) REAL array, dimension ( LDU, NLVL ).
90*         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
91*
92*  DIFR   (input) REAL array, dimension ( LDU, 2 * NLVL ).
93*         On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
94*         distances between singular values on the I-th level and
95*         singular values on the (I -1)-th level, and DIFR(*, 2 * I)
96*         record the normalizing factors of the right singular vectors
97*         matrices of subproblems on I-th level.
98*
99*  Z      (input) REAL array, dimension ( LDU, NLVL ).
100*         On entry, Z(1, I) contains the components of the deflation-
101*         adjusted updating row vector for subproblems on the I-th
102*         level.
103*
104*  POLES  (input) REAL array, dimension ( LDU, 2 * NLVL ).
105*         On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
106*         singular values involved in the secular equations on the I-th
107*         level.
108*
109*  GIVPTR (input) INTEGER array, dimension ( N ).
110*         On entry, GIVPTR( I ) records the number of Givens
111*         rotations performed on the I-th problem on the computation
112*         tree.
113*
114*  GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
115*         On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
116*         locations of Givens rotations performed on the I-th level on
117*         the computation tree.
118*
119*  LDGCOL (input) INTEGER, LDGCOL = > N.
120*         The leading dimension of arrays GIVCOL and PERM.
121*
122*  PERM   (input) INTEGER array, dimension ( LDGCOL, NLVL ).
123*         On entry, PERM(*, I) records permutations done on the I-th
124*         level of the computation tree.
125*
126*  GIVNUM (input) REAL array, dimension ( LDU, 2 * NLVL ).
127*         On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
128*         values of Givens rotations performed on the I-th level on the
129*         computation tree.
130*
131*  C      (input) REAL array, dimension ( N ).
132*         On entry, if the I-th subproblem is not square,
133*         C( I ) contains the C-value of a Givens rotation related to
134*         the right null space of the I-th subproblem.
135*
136*  S      (input) REAL array, dimension ( N ).
137*         On entry, if the I-th subproblem is not square,
138*         S( I ) contains the S-value of a Givens rotation related to
139*         the right null space of the I-th subproblem.
140*
141*  WORK   (workspace) REAL array.
142*         The dimension must be at least N.
143*
144*  IWORK  (workspace) INTEGER array.
145*         The dimension must be at least 3 * N
146*
147*  INFO   (output) INTEGER
148*          = 0:  successful exit.
149*          < 0:  if INFO = -i, the i-th argument had an illegal value.
150*
151*  Further Details
152*  ===============
153*
154*  Based on contributions by
155*     Ming Gu and Ren-Cang Li, Computer Science Division, University of
156*       California at Berkeley, USA
157*     Osni Marques, LBNL/NERSC, USA
158*
159*  =====================================================================
160*
161*     .. Parameters ..
162      REAL               ZERO, ONE
163      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
164*     ..
165*     .. Local Scalars ..
166      INTEGER            I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
167     $                   ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
168     $                   NR, NRF, NRP1, SQRE
169*     ..
170*     .. External Subroutines ..
171      EXTERNAL           SCOPY, SGEMM, SLALS0, SLASDT, XERBLA
172*     ..
173*     .. Executable Statements ..
174*
175*     Test the input parameters.
176*
177      INFO = 0
178*
179      IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
180         INFO = -1
181      ELSE IF( SMLSIZ.LT.3 ) THEN
182         INFO = -2
183      ELSE IF( N.LT.SMLSIZ ) THEN
184         INFO = -3
185      ELSE IF( NRHS.LT.1 ) THEN
186         INFO = -4
187      ELSE IF( LDB.LT.N ) THEN
188         INFO = -6
189      ELSE IF( LDBX.LT.N ) THEN
190         INFO = -8
191      ELSE IF( LDU.LT.N ) THEN
192         INFO = -10
193      ELSE IF( LDGCOL.LT.N ) THEN
194         INFO = -19
195      END IF
196      IF( INFO.NE.0 ) THEN
197         CALL XERBLA( 'SLALSA', -INFO )
198         RETURN
199      END IF
200*
201*     Book-keeping and  setting up the computation tree.
202*
203      INODE = 1
204      NDIML = INODE + N
205      NDIMR = NDIML + N
206*
207      CALL SLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
208     $             IWORK( NDIMR ), SMLSIZ )
209*
210*     The following code applies back the left singular vector factors.
211*     For applying back the right singular vector factors, go to 50.
212*
213      IF( ICOMPQ.EQ.1 ) THEN
214         GO TO 50
215      END IF
216*
217*     The nodes on the bottom level of the tree were solved
218*     by SLASDQ. The corresponding left and right singular vector
219*     matrices are in explicit form. First apply back the left
220*     singular vector matrices.
221*
222      NDB1 = ( ND+1 ) / 2
223      DO 10 I = NDB1, ND
224*
225*        IC : center row of each node
226*        NL : number of rows of left  subproblem
227*        NR : number of rows of right subproblem
228*        NLF: starting row of the left   subproblem
229*        NRF: starting row of the right  subproblem
230*
231         I1 = I - 1
232         IC = IWORK( INODE+I1 )
233         NL = IWORK( NDIML+I1 )
234         NR = IWORK( NDIMR+I1 )
235         NLF = IC - NL
236         NRF = IC + 1
237         CALL SGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
238     $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
239         CALL SGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
240     $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
241   10 CONTINUE
242*
243*     Next copy the rows of B that correspond to unchanged rows
244*     in the bidiagonal matrix to BX.
245*
246      DO 20 I = 1, ND
247         IC = IWORK( INODE+I-1 )
248         CALL SCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
249   20 CONTINUE
250*
251*     Finally go through the left singular vector matrices of all
252*     the other subproblems bottom-up on the tree.
253*
254      J = 2**NLVL
255      SQRE = 0
256*
257      DO 40 LVL = NLVL, 1, -1
258         LVL2 = 2*LVL - 1
259*
260*        find the first node LF and last node LL on
261*        the current level LVL
262*
263         IF( LVL.EQ.1 ) THEN
264            LF = 1
265            LL = 1
266         ELSE
267            LF = 2**( LVL-1 )
268            LL = 2*LF - 1
269         END IF
270         DO 30 I = LF, LL
271            IM1 = I - 1
272            IC = IWORK( INODE+IM1 )
273            NL = IWORK( NDIML+IM1 )
274            NR = IWORK( NDIMR+IM1 )
275            NLF = IC - NL
276            NRF = IC + 1
277            J = J - 1
278            CALL SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
279     $                   B( NLF, 1 ), LDB, PERM( NLF, LVL ),
280     $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
281     $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
282     $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
283     $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
284     $                   INFO )
285   30    CONTINUE
286   40 CONTINUE
287      GO TO 90
288*
289*     ICOMPQ = 1: applying back the right singular vector factors.
290*
291   50 CONTINUE
292*
293*     First now go through the right singular vector matrices of all
294*     the tree nodes top-down.
295*
296      J = 0
297      DO 70 LVL = 1, NLVL
298         LVL2 = 2*LVL - 1
299*
300*        Find the first node LF and last node LL on
301*        the current level LVL.
302*
303         IF( LVL.EQ.1 ) THEN
304            LF = 1
305            LL = 1
306         ELSE
307            LF = 2**( LVL-1 )
308            LL = 2*LF - 1
309         END IF
310         DO 60 I = LL, LF, -1
311            IM1 = I - 1
312            IC = IWORK( INODE+IM1 )
313            NL = IWORK( NDIML+IM1 )
314            NR = IWORK( NDIMR+IM1 )
315            NLF = IC - NL
316            NRF = IC + 1
317            IF( I.EQ.LL ) THEN
318               SQRE = 0
319            ELSE
320               SQRE = 1
321            END IF
322            J = J + 1
323            CALL SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
324     $                   BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
325     $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
326     $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
327     $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
328     $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
329     $                   INFO )
330   60    CONTINUE
331   70 CONTINUE
332*
333*     The nodes on the bottom level of the tree were solved
334*     by SLASDQ. The corresponding right singular vector
335*     matrices are in explicit form. Apply them back.
336*
337      NDB1 = ( ND+1 ) / 2
338      DO 80 I = NDB1, ND
339         I1 = I - 1
340         IC = IWORK( INODE+I1 )
341         NL = IWORK( NDIML+I1 )
342         NR = IWORK( NDIMR+I1 )
343         NLP1 = NL + 1
344         IF( I.EQ.ND ) THEN
345            NRP1 = NR
346         ELSE
347            NRP1 = NR + 1
348         END IF
349         NLF = IC - NL
350         NRF = IC + 1
351         CALL SGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
352     $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
353         CALL SGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
354     $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
355   80 CONTINUE
356*
357   90 CONTINUE
358*
359      RETURN
360*
361*     End of SLALSA
362*
363      END
364