001:       SUBROUTINE SLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
002:      $                   LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
003:      $                   GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
004:      $                   IWORK, INFO )
005: *
006: *  -- LAPACK routine (version 3.2) --
007: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
008: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
009: *     November 2006
010: *
011: *     .. Scalar Arguments ..
012:       INTEGER            ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
013:      $                   SMLSIZ
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
017:      $                   K( * ), PERM( LDGCOL, * )
018:       REAL               B( LDB, * ), BX( LDBX, * ), C( * ),
019:      $                   DIFL( LDU, * ), DIFR( LDU, * ),
020:      $                   GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
021:      $                   U( LDU, * ), VT( LDU, * ), WORK( * ),
022:      $                   Z( LDU, * )
023: *     ..
024: *
025: *  Purpose
026: *  =======
027: *
028: *  SLALSA is an itermediate step in solving the least squares problem
029: *  by computing the SVD of the coefficient matrix in compact form (The
030: *  singular vectors are computed as products of simple orthorgonal
031: *  matrices.).
032: *
033: *  If ICOMPQ = 0, SLALSA applies the inverse of the left singular vector
034: *  matrix of an upper bidiagonal matrix to the right hand side; and if
035: *  ICOMPQ = 1, SLALSA applies the right singular vector matrix to the
036: *  right hand side. The singular vector matrices were generated in
037: *  compact form by SLALSA.
038: *
039: *  Arguments
040: *  =========
041: *
042: *
043: *  ICOMPQ (input) INTEGER
044: *         Specifies whether the left or the right singular vector
045: *         matrix is involved.
046: *         = 0: Left singular vector matrix
047: *         = 1: Right singular vector matrix
048: *
049: *  SMLSIZ (input) INTEGER
050: *         The maximum size of the subproblems at the bottom of the
051: *         computation tree.
052: *
053: *  N      (input) INTEGER
054: *         The row and column dimensions of the upper bidiagonal matrix.
055: *
056: *  NRHS   (input) INTEGER
057: *         The number of columns of B and BX. NRHS must be at least 1.
058: *
059: *  B      (input/output) REAL array, dimension ( LDB, NRHS )
060: *         On input, B contains the right hand sides of the least
061: *         squares problem in rows 1 through M.
062: *         On output, B contains the solution X in rows 1 through N.
063: *
064: *  LDB    (input) INTEGER
065: *         The leading dimension of B in the calling subprogram.
066: *         LDB must be at least max(1,MAX( M, N ) ).
067: *
068: *  BX     (output) REAL array, dimension ( LDBX, NRHS )
069: *         On exit, the result of applying the left or right singular
070: *         vector matrix to B.
071: *
072: *  LDBX   (input) INTEGER
073: *         The leading dimension of BX.
074: *
075: *  U      (input) REAL array, dimension ( LDU, SMLSIZ ).
076: *         On entry, U contains the left singular vector matrices of all
077: *         subproblems at the bottom level.
078: *
079: *  LDU    (input) INTEGER, LDU = > N.
080: *         The leading dimension of arrays U, VT, DIFL, DIFR,
081: *         POLES, GIVNUM, and Z.
082: *
083: *  VT     (input) REAL array, dimension ( LDU, SMLSIZ+1 ).
084: *         On entry, VT' contains the right singular vector matrices of
085: *         all subproblems at the bottom level.
086: *
087: *  K      (input) INTEGER array, dimension ( N ).
088: *
089: *  DIFL   (input) REAL array, dimension ( LDU, NLVL ).
090: *         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
091: *
092: *  DIFR   (input) REAL array, dimension ( LDU, 2 * NLVL ).
093: *         On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
094: *         distances between singular values on the I-th level and
095: *         singular values on the (I -1)-th level, and DIFR(*, 2 * I)
096: *         record the normalizing factors of the right singular vectors
097: *         matrices of subproblems on I-th level.
098: *
099: *  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: