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