001:       SUBROUTINE SLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
002:      $                   RANK, WORK, IWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          UPLO
010:       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
011:       REAL               RCOND
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       REAL               B( LDB, * ), D( * ), E( * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SLALSD uses the singular value decomposition of A to solve the least
022: *  squares problem of finding X to minimize the Euclidean norm of each
023: *  column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
024: *  are N-by-NRHS. The solution X overwrites B.
025: *
026: *  The singular values of A smaller than RCOND times the largest
027: *  singular value are treated as zero in solving the least squares
028: *  problem; in this case a minimum norm solution is returned.
029: *  The actual singular values are returned in D in ascending order.
030: *
031: *  This code makes very mild assumptions about floating point
032: *  arithmetic. It will work on machines with a guard digit in
033: *  add/subtract, or on those binary machines without guard digits
034: *  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
035: *  It could conceivably fail on hexadecimal or decimal machines
036: *  without guard digits, but we know of none.
037: *
038: *  Arguments
039: *  =========
040: *
041: *  UPLO   (input) CHARACTER*1
042: *         = 'U': D and E define an upper bidiagonal matrix.
043: *         = 'L': D and E define a  lower bidiagonal matrix.
044: *
045: *  SMLSIZ (input) INTEGER
046: *         The maximum size of the subproblems at the bottom of the
047: *         computation tree.
048: *
049: *  N      (input) INTEGER
050: *         The dimension of the  bidiagonal matrix.  N >= 0.
051: *
052: *  NRHS   (input) INTEGER
053: *         The number of columns of B. NRHS must be at least 1.
054: *
055: *  D      (input/output) REAL array, dimension (N)
056: *         On entry D contains the main diagonal of the bidiagonal
057: *         matrix. On exit, if INFO = 0, D contains its singular values.
058: *
059: *  E      (input/output) REAL array, dimension (N-1)
060: *         Contains the super-diagonal entries of the bidiagonal matrix.
061: *         On exit, E has been destroyed.
062: *
063: *  B      (input/output) REAL array, dimension (LDB,NRHS)
064: *         On input, B contains the right hand sides of the least
065: *         squares problem. On output, B contains the solution X.
066: *
067: *  LDB    (input) INTEGER
068: *         The leading dimension of B in the calling subprogram.
069: *         LDB must be at least max(1,N).
070: *
071: *  RCOND  (input) REAL
072: *         The singular values of A less than or equal to RCOND times
073: *         the largest singular value are treated as zero in solving
074: *         the least squares problem. If RCOND is negative,
075: *         machine precision is used instead.
076: *         For example, if diag(S)*X=B were the least squares problem,
077: *         where diag(S) is a diagonal matrix of singular values, the
078: *         solution would be X(i) = B(i) / S(i) if S(i) is greater than
079: *         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
080: *         RCOND*max(S).
081: *
082: *  RANK   (output) INTEGER
083: *         The number of singular values of A greater than RCOND times
084: *         the largest singular value.
085: *
086: *  WORK   (workspace) REAL array, dimension at least
087: *         (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
088: *         where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).
089: *
090: *  IWORK  (workspace) INTEGER array, dimension at least
091: *         (3*N*NLVL + 11*N)
092: *
093: *  INFO   (output) INTEGER
094: *         = 0:  successful exit.
095: *         < 0:  if INFO = -i, the i-th argument had an illegal value.
096: *         > 0:  The algorithm failed to compute an singular value while
097: *               working on the submatrix lying in rows and columns
098: *               INFO/(N+1) through MOD(INFO,N+1).
099: *
100: *  Further Details
101: *  ===============
102: *
103: *  Based on contributions by
104: *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
105: *       California at Berkeley, USA
106: *     Osni Marques, LBNL/NERSC, USA
107: *
108: *  =====================================================================
109: *
110: *     .. Parameters ..
111:       REAL               ZERO, ONE, TWO
112:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
113: *     ..
114: *     .. Local Scalars ..
115:       INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
116:      $                   GIVPTR, I, ICMPQ1, ICMPQ2, IWK, J, K, NLVL,
117:      $                   NM1, NSIZE, NSUB, NWORK, PERM, POLES, S, SIZEI,
118:      $                   SMLSZP, SQRE, ST, ST1, U, VT, Z
119:       REAL               CS, EPS, ORGNRM, R, RCND, SN, TOL
120: *     ..
121: *     .. External Functions ..
122:       INTEGER            ISAMAX
123:       REAL               SLAMCH, SLANST
124:       EXTERNAL           ISAMAX, SLAMCH, SLANST
125: *     ..
126: *     .. External Subroutines ..
127:       EXTERNAL           SCOPY, SGEMM, SLACPY, SLALSA, SLARTG, SLASCL,
128:      $                   SLASDA, SLASDQ, SLASET, SLASRT, SROT, XERBLA
129: *     ..
130: *     .. Intrinsic Functions ..
131:       INTRINSIC          ABS, INT, LOG, REAL, SIGN
132: *     ..
133: *     .. Executable Statements ..
134: *
135: *     Test the input parameters.
136: *
137:       INFO = 0
138: *
139:       IF( N.LT.0 ) THEN
140:          INFO = -3
141:       ELSE IF( NRHS.LT.1 ) THEN
142:          INFO = -4
143:       ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
144:          INFO = -8
145:       END IF
146:       IF( INFO.NE.0 ) THEN
147:          CALL XERBLA( 'SLALSD', -INFO )
148:          RETURN
149:       END IF
150: *
151:       EPS = SLAMCH( 'Epsilon' )
152: *
153: *     Set up the tolerance.
154: *
155:       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
156:          RCND = EPS
157:       ELSE
158:          RCND = RCOND
159:       END IF
160: *
161:       RANK = 0
162: *
163: *     Quick return if possible.
164: *
165:       IF( N.EQ.0 ) THEN
166:          RETURN
167:       ELSE IF( N.EQ.1 ) THEN
168:          IF( D( 1 ).EQ.ZERO ) THEN
169:             CALL SLASET( 'A', 1, NRHS, ZERO, ZERO, B, LDB )
170:          ELSE
171:             RANK = 1
172:             CALL SLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
173:             D( 1 ) = ABS( D( 1 ) )
174:          END IF
175:          RETURN
176:       END IF
177: *
178: *     Rotate the matrix if it is lower bidiagonal.
179: *
180:       IF( UPLO.EQ.'L' ) THEN
181:          DO 10 I = 1, N - 1
182:             CALL SLARTG( D( I ), E( I ), CS, SN, R )
183:             D( I ) = R
184:             E( I ) = SN*D( I+1 )
185:             D( I+1 ) = CS*D( I+1 )
186:             IF( NRHS.EQ.1 ) THEN
187:                CALL SROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
188:             ELSE
189:                WORK( I*2-1 ) = CS
190:                WORK( I*2 ) = SN
191:             END IF
192:    10    CONTINUE
193:          IF( NRHS.GT.1 ) THEN
194:             DO 30 I = 1, NRHS
195:                DO 20 J = 1, N - 1
196:                   CS = WORK( J*2-1 )
197:                   SN = WORK( J*2 )
198:                   CALL SROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
199:    20          CONTINUE
200:    30       CONTINUE
201:          END IF
202:       END IF
203: *
204: *     Scale.
205: *
206:       NM1 = N - 1
207:       ORGNRM = SLANST( 'M', N, D, E )
208:       IF( ORGNRM.EQ.ZERO ) THEN
209:          CALL SLASET( 'A', N, NRHS, ZERO, ZERO, B, LDB )
210:          RETURN
211:       END IF
212: *
213:       CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
214:       CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
215: *
216: *     If N is smaller than the minimum divide size SMLSIZ, then solve
217: *     the problem with another solver.
218: *
219:       IF( N.LE.SMLSIZ ) THEN
220:          NWORK = 1 + N*N
221:          CALL SLASET( 'A', N, N, ZERO, ONE, WORK, N )
222:          CALL SLASDQ( 'U', 0, N, N, 0, NRHS, D, E, WORK, N, WORK, N, B,
223:      $                LDB, WORK( NWORK ), INFO )
224:          IF( INFO.NE.0 ) THEN
225:             RETURN
226:          END IF
227:          TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
228:          DO 40 I = 1, N
229:             IF( D( I ).LE.TOL ) THEN
230:                CALL SLASET( 'A', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
231:             ELSE
232:                CALL SLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
233:      $                      LDB, INFO )
234:                RANK = RANK + 1
235:             END IF
236:    40    CONTINUE
237:          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
238:      $               WORK( NWORK ), N )
239:          CALL SLACPY( 'A', N, NRHS, WORK( NWORK ), N, B, LDB )
240: *
241: *        Unscale.
242: *
243:          CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
244:          CALL SLASRT( 'D', N, D, INFO )
245:          CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
246: *
247:          RETURN
248:       END IF
249: *
250: *     Book-keeping and setting up some constants.
251: *
252:       NLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
253: *
254:       SMLSZP = SMLSIZ + 1
255: *
256:       U = 1
257:       VT = 1 + SMLSIZ*N
258:       DIFL = VT + SMLSZP*N
259:       DIFR = DIFL + NLVL*N
260:       Z = DIFR + NLVL*N*2
261:       C = Z + NLVL*N
262:       S = C + N
263:       POLES = S + N
264:       GIVNUM = POLES + 2*NLVL*N
265:       BX = GIVNUM + 2*NLVL*N
266:       NWORK = BX + N*NRHS
267: *
268:       SIZEI = 1 + N
269:       K = SIZEI + N
270:       GIVPTR = K + N
271:       PERM = GIVPTR + N
272:       GIVCOL = PERM + NLVL*N
273:       IWK = GIVCOL + NLVL*N*2
274: *
275:       ST = 1
276:       SQRE = 0
277:       ICMPQ1 = 1
278:       ICMPQ2 = 0
279:       NSUB = 0
280: *
281:       DO 50 I = 1, N
282:          IF( ABS( D( I ) ).LT.EPS ) THEN
283:             D( I ) = SIGN( EPS, D( I ) )
284:          END IF
285:    50 CONTINUE
286: *
287:       DO 60 I = 1, NM1
288:          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
289:             NSUB = NSUB + 1
290:             IWORK( NSUB ) = ST
291: *
292: *           Subproblem found. First determine its size and then
293: *           apply divide and conquer on it.
294: *
295:             IF( I.LT.NM1 ) THEN
296: *
297: *              A subproblem with E(I) small for I < NM1.
298: *
299:                NSIZE = I - ST + 1
300:                IWORK( SIZEI+NSUB-1 ) = NSIZE
301:             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
302: *
303: *              A subproblem with E(NM1) not too small but I = NM1.
304: *
305:                NSIZE = N - ST + 1
306:                IWORK( SIZEI+NSUB-1 ) = NSIZE
307:             ELSE
308: *
309: *              A subproblem with E(NM1) small. This implies an
310: *              1-by-1 subproblem at D(N), which is not solved
311: *              explicitly.
312: *
313:                NSIZE = I - ST + 1
314:                IWORK( SIZEI+NSUB-1 ) = NSIZE
315:                NSUB = NSUB + 1
316:                IWORK( NSUB ) = N
317:                IWORK( SIZEI+NSUB-1 ) = 1
318:                CALL SCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
319:             END IF
320:             ST1 = ST - 1
321:             IF( NSIZE.EQ.1 ) THEN
322: *
323: *              This is a 1-by-1 subproblem and is not solved
324: *              explicitly.
325: *
326:                CALL SCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
327:             ELSE IF( NSIZE.LE.SMLSIZ ) THEN
328: *
329: *              This is a small subproblem and is solved by SLASDQ.
330: *
331:                CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
332:      $                      WORK( VT+ST1 ), N )
333:                CALL SLASDQ( 'U', 0, NSIZE, NSIZE, 0, NRHS, D( ST ),
334:      $                      E( ST ), WORK( VT+ST1 ), N, WORK( NWORK ),
335:      $                      N, B( ST, 1 ), LDB, WORK( NWORK ), INFO )
336:                IF( INFO.NE.0 ) THEN
337:                   RETURN
338:                END IF
339:                CALL SLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
340:      $                      WORK( BX+ST1 ), N )
341:             ELSE
342: *
343: *              A large problem. Solve it using divide and conquer.
344: *
345:                CALL SLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
346:      $                      E( ST ), WORK( U+ST1 ), N, WORK( VT+ST1 ),
347:      $                      IWORK( K+ST1 ), WORK( DIFL+ST1 ),
348:      $                      WORK( DIFR+ST1 ), WORK( Z+ST1 ),
349:      $                      WORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
350:      $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
351:      $                      WORK( GIVNUM+ST1 ), WORK( C+ST1 ),
352:      $                      WORK( S+ST1 ), WORK( NWORK ), IWORK( IWK ),
353:      $                      INFO )
354:                IF( INFO.NE.0 ) THEN
355:                   RETURN
356:                END IF
357:                BXST = BX + ST1
358:                CALL SLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
359:      $                      LDB, WORK( BXST ), N, WORK( U+ST1 ), N,
360:      $                      WORK( VT+ST1 ), IWORK( K+ST1 ),
361:      $                      WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
362:      $                      WORK( Z+ST1 ), WORK( POLES+ST1 ),
363:      $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
364:      $                      IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
365:      $                      WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
366:      $                      IWORK( IWK ), INFO )
367:                IF( INFO.NE.0 ) THEN
368:                   RETURN
369:                END IF
370:             END IF
371:             ST = I + 1
372:          END IF
373:    60 CONTINUE
374: *
375: *     Apply the singular values and treat the tiny ones as zero.
376: *
377:       TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
378: *
379:       DO 70 I = 1, N
380: *
381: *        Some of the elements in D can be negative because 1-by-1
382: *        subproblems were not solved explicitly.
383: *
384:          IF( ABS( D( I ) ).LE.TOL ) THEN
385:             CALL SLASET( 'A', 1, NRHS, ZERO, ZERO, WORK( BX+I-1 ), N )
386:          ELSE
387:             RANK = RANK + 1
388:             CALL SLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
389:      $                   WORK( BX+I-1 ), N, INFO )
390:          END IF
391:          D( I ) = ABS( D( I ) )
392:    70 CONTINUE
393: *
394: *     Now apply back the right singular vectors.
395: *
396:       ICMPQ2 = 1
397:       DO 80 I = 1, NSUB
398:          ST = IWORK( I )
399:          ST1 = ST - 1
400:          NSIZE = IWORK( SIZEI+I-1 )
401:          BXST = BX + ST1
402:          IF( NSIZE.EQ.1 ) THEN
403:             CALL SCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
404:          ELSE IF( NSIZE.LE.SMLSIZ ) THEN
405:             CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
406:      $                  WORK( VT+ST1 ), N, WORK( BXST ), N, ZERO,
407:      $                  B( ST, 1 ), LDB )
408:          ELSE
409:             CALL SLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
410:      $                   B( ST, 1 ), LDB, WORK( U+ST1 ), N,
411:      $                   WORK( VT+ST1 ), IWORK( K+ST1 ),
412:      $                   WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
413:      $                   WORK( Z+ST1 ), WORK( POLES+ST1 ),
414:      $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
415:      $                   IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
416:      $                   WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
417:      $                   IWORK( IWK ), INFO )
418:             IF( INFO.NE.0 ) THEN
419:                RETURN
420:             END IF
421:          END IF
422:    80 CONTINUE
423: *
424: *     Unscale and sort the singular values.
425: *
426:       CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
427:       CALL SLASRT( 'D', N, D, INFO )
428:       CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
429: *
430:       RETURN
431: *
432: *     End of SLALSD
433: *
434:       END
435: