001:       SUBROUTINE CLALS0( 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: *     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               DIFL( * ), DIFR( LDGNUM, * ),
017:      $                   GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
018:      $                   RWORK( * ), Z( * )
019:       COMPLEX            B( LDB, * ), BX( LDBX, * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  CLALS0 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) COMPLEX 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) COMPLEX 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) REAL 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) REAL 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) REAL 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) REAL 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) REAL 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) REAL
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) REAL
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: *  RWORK  (workspace) REAL array, dimension
151: *         ( K*(1+NRHS) + 2*NRHS )
152: *
153: *  INFO   (output) INTEGER
154: *          = 0:  successful exit.
155: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
156: *
157: *  Further Details
158: *  ===============
159: *
160: *  Based on contributions by
161: *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
162: *       California at Berkeley, USA
163: *     Osni Marques, LBNL/NERSC, USA
164: *
165: *  =====================================================================
166: *
167: *     .. Parameters ..
168:       REAL               ONE, ZERO, NEGONE
169:       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0, NEGONE = -1.0E0 )
170: *     ..
171: *     .. Local Scalars ..
172:       INTEGER            I, J, JCOL, JROW, M, N, NLP1
173:       REAL               DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
174: *     ..
175: *     .. External Subroutines ..
176:       EXTERNAL           CCOPY, CLACPY, CLASCL, CSROT, CSSCAL, SGEMV,
177:      $                   XERBLA
178: *     ..
179: *     .. External Functions ..
180:       REAL               SLAMC3, SNRM2
181:       EXTERNAL           SLAMC3, SNRM2
182: *     ..
183: *     .. Intrinsic Functions ..
184:       INTRINSIC          AIMAG, CMPLX, MAX, REAL
185: *     ..
186: *     .. Executable Statements ..
187: *
188: *     Test the input parameters.
189: *
190:       INFO = 0
191: *
192:       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
193:          INFO = -1
194:       ELSE IF( NL.LT.1 ) THEN
195:          INFO = -2
196:       ELSE IF( NR.LT.1 ) THEN
197:          INFO = -3
198:       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
199:          INFO = -4
200:       END IF
201: *
202:       N = NL + NR + 1
203: *
204:       IF( NRHS.LT.1 ) THEN
205:          INFO = -5
206:       ELSE IF( LDB.LT.N ) THEN
207:          INFO = -7
208:       ELSE IF( LDBX.LT.N ) THEN
209:          INFO = -9
210:       ELSE IF( GIVPTR.LT.0 ) THEN
211:          INFO = -11
212:       ELSE IF( LDGCOL.LT.N ) THEN
213:          INFO = -13
214:       ELSE IF( LDGNUM.LT.N ) THEN
215:          INFO = -15
216:       ELSE IF( K.LT.1 ) THEN
217:          INFO = -20
218:       END IF
219:       IF( INFO.NE.0 ) THEN
220:          CALL XERBLA( 'CLALS0', -INFO )
221:          RETURN
222:       END IF
223: *
224:       M = N + SQRE
225:       NLP1 = NL + 1
226: *
227:       IF( ICOMPQ.EQ.0 ) THEN
228: *
229: *        Apply back orthogonal transformations from the left.
230: *
231: *        Step (1L): apply back the Givens rotations performed.
232: *
233:          DO 10 I = 1, GIVPTR
234:             CALL CSROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
235:      $                  B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
236:      $                  GIVNUM( I, 1 ) )
237:    10    CONTINUE
238: *
239: *        Step (2L): permute rows of B.
240: *
241:          CALL CCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
242:          DO 20 I = 2, N
243:             CALL CCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
244:    20    CONTINUE
245: *
246: *        Step (3L): apply the inverse of the left singular vector
247: *        matrix to BX.
248: *
249:          IF( K.EQ.1 ) THEN
250:             CALL CCOPY( NRHS, BX, LDBX, B, LDB )
251:             IF( Z( 1 ).LT.ZERO ) THEN
252:                CALL CSSCAL( NRHS, NEGONE, B, LDB )
253:             END IF
254:          ELSE
255:             DO 100 J = 1, K
256:                DIFLJ = DIFL( J )
257:                DJ = POLES( J, 1 )
258:                DSIGJ = -POLES( J, 2 )
259:                IF( J.LT.K ) THEN
260:                   DIFRJ = -DIFR( J, 1 )
261:                   DSIGJP = -POLES( J+1, 2 )
262:                END IF
263:                IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
264:      $              THEN
265:                   RWORK( J ) = ZERO
266:                ELSE
267:                   RWORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
268:      $                         ( POLES( J, 2 )+DJ )
269:                END IF
270:                DO 30 I = 1, J - 1
271:                   IF( ( Z( I ).EQ.ZERO ) .OR.
272:      $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
273:                      RWORK( I ) = ZERO
274:                   ELSE
275:                      RWORK( I ) = POLES( I, 2 )*Z( I ) /
276:      $                            ( SLAMC3( POLES( I, 2 ), DSIGJ )-
277:      $                            DIFLJ ) / ( POLES( I, 2 )+DJ )
278:                   END IF
279:    30          CONTINUE
280:                DO 40 I = J + 1, K
281:                   IF( ( Z( I ).EQ.ZERO ) .OR.
282:      $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
283:                      RWORK( I ) = ZERO
284:                   ELSE
285:                      RWORK( I ) = POLES( I, 2 )*Z( I ) /
286:      $                            ( SLAMC3( POLES( I, 2 ), DSIGJP )+
287:      $                            DIFRJ ) / ( POLES( I, 2 )+DJ )
288:                   END IF
289:    40          CONTINUE
290:                RWORK( 1 ) = NEGONE
291:                TEMP = SNRM2( K, RWORK, 1 )
292: *
293: *              Since B and BX are complex, the following call to SGEMV
294: *              is performed in two steps (real and imaginary parts).
295: *
296: *              CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
297: *    $                     B( J, 1 ), LDB )
298: *
299:                I = K + NRHS*2
300:                DO 60 JCOL = 1, NRHS
301:                   DO 50 JROW = 1, K
302:                      I = I + 1
303:                      RWORK( I ) = REAL( BX( JROW, JCOL ) )
304:    50             CONTINUE
305:    60          CONTINUE
306:                CALL SGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
307:      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
308:                I = K + NRHS*2
309:                DO 80 JCOL = 1, NRHS
310:                   DO 70 JROW = 1, K
311:                      I = I + 1
312:                      RWORK( I ) = AIMAG( BX( JROW, JCOL ) )
313:    70             CONTINUE
314:    80          CONTINUE
315:                CALL SGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
316:      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
317:                DO 90 JCOL = 1, NRHS
318:                   B( J, JCOL ) = CMPLX( RWORK( JCOL+K ),
319:      $                           RWORK( JCOL+K+NRHS ) )
320:    90          CONTINUE
321:                CALL CLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
322:      $                      LDB, INFO )
323:   100       CONTINUE
324:          END IF
325: *
326: *        Move the deflated rows of BX to B also.
327: *
328:          IF( K.LT.MAX( M, N ) )
329:      $      CALL CLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
330:      $                   B( K+1, 1 ), LDB )
331:       ELSE
332: *
333: *        Apply back the right orthogonal transformations.
334: *
335: *        Step (1R): apply back the new right singular vector matrix
336: *        to B.
337: *
338:          IF( K.EQ.1 ) THEN
339:             CALL CCOPY( NRHS, B, LDB, BX, LDBX )
340:          ELSE
341:             DO 180 J = 1, K
342:                DSIGJ = POLES( J, 2 )
343:                IF( Z( J ).EQ.ZERO ) THEN
344:                   RWORK( J ) = ZERO
345:                ELSE
346:                   RWORK( J ) = -Z( J ) / DIFL( J ) /
347:      $                         ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
348:                END IF
349:                DO 110 I = 1, J - 1
350:                   IF( Z( J ).EQ.ZERO ) THEN
351:                      RWORK( I ) = ZERO
352:                   ELSE
353:                      RWORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I+1,
354:      $                            2 ) )-DIFR( I, 1 ) ) /
355:      $                            ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
356:                   END IF
357:   110          CONTINUE
358:                DO 120 I = J + 1, K
359:                   IF( Z( J ).EQ.ZERO ) THEN
360:                      RWORK( I ) = ZERO
361:                   ELSE
362:                      RWORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I,
363:      $                            2 ) )-DIFL( I ) ) /
364:      $                            ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
365:                   END IF
366:   120          CONTINUE
367: *
368: *              Since B and BX are complex, the following call to SGEMV
369: *              is performed in two steps (real and imaginary parts).
370: *
371: *              CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
372: *    $                     BX( J, 1 ), LDBX )
373: *
374:                I = K + NRHS*2
375:                DO 140 JCOL = 1, NRHS
376:                   DO 130 JROW = 1, K
377:                      I = I + 1
378:                      RWORK( I ) = REAL( B( JROW, JCOL ) )
379:   130             CONTINUE
380:   140          CONTINUE
381:                CALL SGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
382:      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
383:                I = K + NRHS*2
384:                DO 160 JCOL = 1, NRHS
385:                   DO 150 JROW = 1, K
386:                      I = I + 1
387:                      RWORK( I ) = AIMAG( B( JROW, JCOL ) )
388:   150             CONTINUE
389:   160          CONTINUE
390:                CALL SGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
391:      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
392:                DO 170 JCOL = 1, NRHS
393:                   BX( J, JCOL ) = CMPLX( RWORK( JCOL+K ),
394:      $                            RWORK( JCOL+K+NRHS ) )
395:   170          CONTINUE
396:   180       CONTINUE
397:          END IF
398: *
399: *        Step (2R): if SQRE = 1, apply back the rotation that is
400: *        related to the right null space of the subproblem.
401: *
402:          IF( SQRE.EQ.1 ) THEN
403:             CALL CCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
404:             CALL CSROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
405:          END IF
406:          IF( K.LT.MAX( M, N ) )
407:      $      CALL CLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB,
408:      $                   BX( K+1, 1 ), LDBX )
409: *
410: *        Step (3R): permute rows of B.
411: *
412:          CALL CCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
413:          IF( SQRE.EQ.1 ) THEN
414:             CALL CCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
415:          END IF
416:          DO 190 I = 2, N
417:             CALL CCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
418:   190    CONTINUE
419: *
420: *        Step (4R): apply back the Givens rotations performed.
421: *
422:          DO 200 I = GIVPTR, 1, -1
423:             CALL CSROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
424:      $                  B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
425:      $                  -GIVNUM( I, 1 ) )
426:   200    CONTINUE
427:       END IF
428: *
429:       RETURN
430: *
431: *     End of CLALS0
432: *
433:       END
434: