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