001:       SUBROUTINE SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
002:      $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
003:      $                   POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO )
004: *
005: *  -- LAPACK routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
011:      $                   LDGNUM, NL, NR, NRHS, SQRE
012:       REAL               C, S
013: *     ..
014: *     .. Array Arguments ..
015:       INTEGER            GIVCOL( LDGCOL, * ), PERM( * )
016:       REAL               B( LDB, * ), BX( LDBX, * ), DIFL( * ),
017:      $                   DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
018:      $                   POLES( LDGNUM, * ), WORK( * ), Z( * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  SLALS0 applies back the multiplying factors of either the left or the
025: *  right singular vector matrix of a diagonal matrix appended by a row
026: *  to the right hand side matrix B in solving the least squares problem
027: *  using the divide-and-conquer SVD approach.
028: *
029: *  For the left singular vector matrix, three types of orthogonal
030: *  matrices are involved:
031: *
032: *  (1L) Givens rotations: the number of such rotations is GIVPTR; the
033: *       pairs of columns/rows they were applied to are stored in GIVCOL;
034: *       and the C- and S-values of these rotations are stored in GIVNUM.
035: *
036: *  (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
037: *       row, and for J=2:N, PERM(J)-th row of B is to be moved to the
038: *       J-th row.
039: *
040: *  (3L) The left singular vector matrix of the remaining matrix.
041: *
042: *  For the right singular vector matrix, four types of orthogonal
043: *  matrices are involved:
044: *
045: *  (1R) The right singular vector matrix of the remaining matrix.
046: *
047: *  (2R) If SQRE = 1, one extra Givens rotation to generate the right
048: *       null space.
049: *
050: *  (3R) The inverse transformation of (2L).
051: *
052: *  (4R) The inverse transformation of (1L).
053: *
054: *  Arguments
055: *  =========
056: *
057: *  ICOMPQ (input) INTEGER
058: *         Specifies whether singular vectors are to be computed in
059: *         factored form:
060: *         = 0: Left singular vector matrix.
061: *         = 1: Right singular vector matrix.
062: *
063: *  NL     (input) INTEGER
064: *         The row dimension of the upper block. NL >= 1.
065: *
066: *  NR     (input) INTEGER
067: *         The row dimension of the lower block. NR >= 1.
068: *
069: *  SQRE   (input) INTEGER
070: *         = 0: the lower block is an NR-by-NR square matrix.
071: *         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
072: *
073: *         The bidiagonal matrix has row dimension N = NL + NR + 1,
074: *         and column dimension M = N + SQRE.
075: *
076: *  NRHS   (input) INTEGER
077: *         The number of columns of B and BX. NRHS must be at least 1.
078: *
079: *  B      (input/output) REAL array, dimension ( LDB, NRHS )
080: *         On input, B contains the right hand sides of the least
081: *         squares problem in rows 1 through M. On output, B contains
082: *         the solution X in rows 1 through N.
083: *
084: *  LDB    (input) INTEGER
085: *         The leading dimension of B. LDB must be at least
086: *         max(1,MAX( M, N ) ).
087: *
088: *  BX     (workspace) REAL array, dimension ( LDBX, NRHS )
089: *
090: *  LDBX   (input) INTEGER
091: *         The leading dimension of BX.
092: *
093: *  PERM   (input) INTEGER array, dimension ( N )
094: *         The permutations (from deflation and sorting) applied
095: *         to the two blocks.
096: *
097: *  GIVPTR (input) INTEGER
098: *         The number of Givens rotations which took place in this
099: *         subproblem.
100: *
101: *  GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 )
102: *         Each pair of numbers indicates a pair of rows/columns
103: *         involved in a Givens rotation.
104: *
105: *  LDGCOL (input) INTEGER
106: *         The leading dimension of GIVCOL, must be at least N.
107: *
108: *  GIVNUM (input) REAL array, dimension ( LDGNUM, 2 )
109: *         Each number indicates the C or S value used in the
110: *         corresponding Givens rotation.
111: *
112: *  LDGNUM (input) INTEGER
113: *         The leading dimension of arrays DIFR, POLES and
114: *         GIVNUM, must be at least K.
115: *
116: *  POLES  (input) REAL array, dimension ( LDGNUM, 2 )
117: *         On entry, POLES(1:K, 1) contains the new singular
118: *         values obtained from solving the secular equation, and
119: *         POLES(1:K, 2) is an array containing the poles in the secular
120: *         equation.
121: *
122: *  DIFL   (input) REAL array, dimension ( K ).
123: *         On entry, DIFL(I) is the distance between I-th updated
124: *         (undeflated) singular value and the I-th (undeflated) old
125: *         singular value.
126: *
127: *  DIFR   (input) REAL array, dimension ( LDGNUM, 2 ).
128: *         On entry, DIFR(I, 1) contains the distances between I-th
129: *         updated (undeflated) singular value and the I+1-th
130: *         (undeflated) old singular value. And DIFR(I, 2) is the
131: *         normalizing factor for the I-th right singular vector.
132: *
133: *  Z      (input) REAL array, dimension ( K )
134: *         Contain the components of the deflation-adjusted updating row
135: *         vector.
136: *
137: *  K      (input) INTEGER
138: *         Contains the dimension of the non-deflated matrix,
139: *         This is the order of the related secular equation. 1 <= K <=N.
140: *
141: *  C      (input) REAL
142: *         C contains garbage if SQRE =0 and the C-value of a Givens
143: *         rotation related to the right null space if SQRE = 1.
144: *
145: *  S      (input) REAL
146: *         S contains garbage if SQRE =0 and the S-value of a Givens
147: *         rotation related to the right null space if SQRE = 1.
148: *
149: *  WORK   (workspace) REAL array, dimension ( K )
150: *
151: *  INFO   (output) INTEGER
152: *          = 0:  successful exit.
153: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
154: *
155: *  Further Details
156: *  ===============
157: *
158: *  Based on contributions by
159: *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
160: *       California at Berkeley, USA
161: *     Osni Marques, LBNL/NERSC, USA
162: *
163: *  =====================================================================
164: *
165: *     .. Parameters ..
166:       REAL               ONE, ZERO, NEGONE
167:       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0, NEGONE = -1.0E0 )
168: *     ..
169: *     .. Local Scalars ..
170:       INTEGER            I, J, M, N, NLP1
171:       REAL               DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
172: *     ..
173: *     .. External Subroutines ..
174:       EXTERNAL           SCOPY, SGEMV, SLACPY, SLASCL, SROT, SSCAL,
175:      $                   XERBLA
176: *     ..
177: *     .. External Functions ..
178:       REAL               SLAMC3, SNRM2
179:       EXTERNAL           SLAMC3, SNRM2
180: *     ..
181: *     .. Intrinsic Functions ..
182:       INTRINSIC          MAX
183: *     ..
184: *     .. Executable Statements ..
185: *
186: *     Test the input parameters.
187: *
188:       INFO = 0
189: *
190:       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
191:          INFO = -1
192:       ELSE IF( NL.LT.1 ) THEN
193:          INFO = -2
194:       ELSE IF( NR.LT.1 ) THEN
195:          INFO = -3
196:       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
197:          INFO = -4
198:       END IF
199: *
200:       N = NL + NR + 1
201: *
202:       IF( NRHS.LT.1 ) THEN
203:          INFO = -5
204:       ELSE IF( LDB.LT.N ) THEN
205:          INFO = -7
206:       ELSE IF( LDBX.LT.N ) THEN
207:          INFO = -9
208:       ELSE IF( GIVPTR.LT.0 ) THEN
209:          INFO = -11
210:       ELSE IF( LDGCOL.LT.N ) THEN
211:          INFO = -13
212:       ELSE IF( LDGNUM.LT.N ) THEN
213:          INFO = -15
214:       ELSE IF( K.LT.1 ) THEN
215:          INFO = -20
216:       END IF
217:       IF( INFO.NE.0 ) THEN
218:          CALL XERBLA( 'SLALS0', -INFO )
219:          RETURN
220:       END IF
221: *
222:       M = N + SQRE
223:       NLP1 = NL + 1
224: *
225:       IF( ICOMPQ.EQ.0 ) THEN
226: *
227: *        Apply back orthogonal transformations from the left.
228: *
229: *        Step (1L): apply back the Givens rotations performed.
230: *
231:          DO 10 I = 1, GIVPTR
232:             CALL SROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
233:      $                 B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
234:      $                 GIVNUM( I, 1 ) )
235:    10    CONTINUE
236: *
237: *        Step (2L): permute rows of B.
238: *
239:          CALL SCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
240:          DO 20 I = 2, N
241:             CALL SCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
242:    20    CONTINUE
243: *
244: *        Step (3L): apply the inverse of the left singular vector
245: *        matrix to BX.
246: *
247:          IF( K.EQ.1 ) THEN
248:             CALL SCOPY( NRHS, BX, LDBX, B, LDB )
249:             IF( Z( 1 ).LT.ZERO ) THEN
250:                CALL SSCAL( NRHS, NEGONE, B, LDB )
251:             END IF
252:          ELSE
253:             DO 50 J = 1, K
254:                DIFLJ = DIFL( J )
255:                DJ = POLES( J, 1 )
256:                DSIGJ = -POLES( J, 2 )
257:                IF( J.LT.K ) THEN
258:                   DIFRJ = -DIFR( J, 1 )
259:                   DSIGJP = -POLES( J+1, 2 )
260:                END IF
261:                IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
262:      $              THEN
263:                   WORK( J ) = ZERO
264:                ELSE
265:                   WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
266:      $                        ( POLES( J, 2 )+DJ )
267:                END IF
268:                DO 30 I = 1, J - 1
269:                   IF( ( Z( I ).EQ.ZERO ) .OR.
270:      $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
271:                      WORK( I ) = ZERO
272:                   ELSE
273:                      WORK( I ) = POLES( I, 2 )*Z( I ) /
274:      $                           ( SLAMC3( POLES( I, 2 ), DSIGJ )-
275:      $                           DIFLJ ) / ( POLES( I, 2 )+DJ )
276:                   END IF
277:    30          CONTINUE
278:                DO 40 I = J + 1, K
279:                   IF( ( Z( I ).EQ.ZERO ) .OR.
280:      $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
281:                      WORK( I ) = ZERO
282:                   ELSE
283:                      WORK( I ) = POLES( I, 2 )*Z( I ) /
284:      $                           ( SLAMC3( POLES( I, 2 ), DSIGJP )+
285:      $                           DIFRJ ) / ( POLES( I, 2 )+DJ )
286:                   END IF
287:    40          CONTINUE
288:                WORK( 1 ) = NEGONE
289:                TEMP = SNRM2( K, WORK, 1 )
290:                CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
291:      $                     B( J, 1 ), LDB )
292:                CALL SLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
293:      $                      LDB, INFO )
294:    50       CONTINUE
295:          END IF
296: *
297: *        Move the deflated rows of BX to B also.
298: *
299:          IF( K.LT.MAX( M, N ) )
300:      $      CALL SLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
301:      $                   B( K+1, 1 ), LDB )
302:       ELSE
303: *
304: *        Apply back the right orthogonal transformations.
305: *
306: *        Step (1R): apply back the new right singular vector matrix
307: *        to B.
308: *
309:          IF( K.EQ.1 ) THEN
310:             CALL SCOPY( NRHS, B, LDB, BX, LDBX )
311:          ELSE
312:             DO 80 J = 1, K
313:                DSIGJ = POLES( J, 2 )
314:                IF( Z( J ).EQ.ZERO ) THEN
315:                   WORK( J ) = ZERO
316:                ELSE
317:                   WORK( J ) = -Z( J ) / DIFL( J ) /
318:      $                        ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
319:                END IF
320:                DO 60 I = 1, J - 1
321:                   IF( Z( J ).EQ.ZERO ) THEN
322:                      WORK( I ) = ZERO
323:                   ELSE
324:                      WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I+1,
325:      $                           2 ) )-DIFR( I, 1 ) ) /
326:      $                           ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
327:                   END IF
328:    60          CONTINUE
329:                DO 70 I = J + 1, K
330:                   IF( Z( J ).EQ.ZERO ) THEN
331:                      WORK( I ) = ZERO
332:                   ELSE
333:                      WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I,
334:      $                           2 ) )-DIFL( I ) ) /
335:      $                           ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
336:                   END IF
337:    70          CONTINUE
338:                CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
339:      $                     BX( J, 1 ), LDBX )
340:    80       CONTINUE
341:          END IF
342: *
343: *        Step (2R): if SQRE = 1, apply back the rotation that is
344: *        related to the right null space of the subproblem.
345: *
346:          IF( SQRE.EQ.1 ) THEN
347:             CALL SCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
348:             CALL SROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
349:          END IF
350:          IF( K.LT.MAX( M, N ) )
351:      $      CALL SLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
352:      $                   LDBX )
353: *
354: *        Step (3R): permute rows of B.
355: *
356:          CALL SCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
357:          IF( SQRE.EQ.1 ) THEN
358:             CALL SCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
359:          END IF
360:          DO 90 I = 2, N
361:             CALL SCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
362:    90    CONTINUE
363: *
364: *        Step (4R): apply back the Givens rotations performed.
365: *
366:          DO 100 I = GIVPTR, 1, -1
367:             CALL SROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
368:      $                 B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
369:      $                 -GIVNUM( I, 1 ) )
370:   100    CONTINUE
371:       END IF
372: *
373:       RETURN
374: *
375: *     End of SLALS0
376: *
377:       END
378: