1*> \brief \b SLALSD uses the singular value decomposition of A to solve the least squares problem.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLALSD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slalsd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slalsd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slalsd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
22*                          RANK, WORK, IWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          UPLO
26*       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
27*       REAL               RCOND
28*       ..
29*       .. Array Arguments ..
30*       INTEGER            IWORK( * )
31*       REAL               B( LDB, * ), D( * ), E( * ), WORK( * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> SLALSD uses the singular value decomposition of A to solve the least
41*> squares problem of finding X to minimize the Euclidean norm of each
42*> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
43*> are N-by-NRHS. The solution X overwrites B.
44*>
45*> The singular values of A smaller than RCOND times the largest
46*> singular value are treated as zero in solving the least squares
47*> problem; in this case a minimum norm solution is returned.
48*> The actual singular values are returned in D in ascending order.
49*>
50*> This code makes very mild assumptions about floating point
51*> arithmetic. It will work on machines with a guard digit in
52*> add/subtract, or on those binary machines without guard digits
53*> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
54*> It could conceivably fail on hexadecimal or decimal machines
55*> without guard digits, but we know of none.
56*> \endverbatim
57*
58*  Arguments:
59*  ==========
60*
61*> \param[in] UPLO
62*> \verbatim
63*>          UPLO is CHARACTER*1
64*>         = 'U': D and E define an upper bidiagonal matrix.
65*>         = 'L': D and E define a  lower bidiagonal matrix.
66*> \endverbatim
67*>
68*> \param[in] SMLSIZ
69*> \verbatim
70*>          SMLSIZ is INTEGER
71*>         The maximum size of the subproblems at the bottom of the
72*>         computation tree.
73*> \endverbatim
74*>
75*> \param[in] N
76*> \verbatim
77*>          N is INTEGER
78*>         The dimension of the  bidiagonal matrix.  N >= 0.
79*> \endverbatim
80*>
81*> \param[in] NRHS
82*> \verbatim
83*>          NRHS is INTEGER
84*>         The number of columns of B. NRHS must be at least 1.
85*> \endverbatim
86*>
87*> \param[in,out] D
88*> \verbatim
89*>          D is REAL array, dimension (N)
90*>         On entry D contains the main diagonal of the bidiagonal
91*>         matrix. On exit, if INFO = 0, D contains its singular values.
92*> \endverbatim
93*>
94*> \param[in,out] E
95*> \verbatim
96*>          E is REAL array, dimension (N-1)
97*>         Contains the super-diagonal entries of the bidiagonal matrix.
98*>         On exit, E has been destroyed.
99*> \endverbatim
100*>
101*> \param[in,out] B
102*> \verbatim
103*>          B is REAL array, dimension (LDB,NRHS)
104*>         On input, B contains the right hand sides of the least
105*>         squares problem. On output, B contains the solution X.
106*> \endverbatim
107*>
108*> \param[in] LDB
109*> \verbatim
110*>          LDB is INTEGER
111*>         The leading dimension of B in the calling subprogram.
112*>         LDB must be at least max(1,N).
113*> \endverbatim
114*>
115*> \param[in] RCOND
116*> \verbatim
117*>          RCOND is REAL
118*>         The singular values of A less than or equal to RCOND times
119*>         the largest singular value are treated as zero in solving
120*>         the least squares problem. If RCOND is negative,
121*>         machine precision is used instead.
122*>         For example, if diag(S)*X=B were the least squares problem,
123*>         where diag(S) is a diagonal matrix of singular values, the
124*>         solution would be X(i) = B(i) / S(i) if S(i) is greater than
125*>         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
126*>         RCOND*max(S).
127*> \endverbatim
128*>
129*> \param[out] RANK
130*> \verbatim
131*>          RANK is INTEGER
132*>         The number of singular values of A greater than RCOND times
133*>         the largest singular value.
134*> \endverbatim
135*>
136*> \param[out] WORK
137*> \verbatim
138*>          WORK is REAL array, dimension at least
139*>         (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
140*>         where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).
141*> \endverbatim
142*>
143*> \param[out] IWORK
144*> \verbatim
145*>          IWORK is INTEGER array, dimension at least
146*>         (3*N*NLVL + 11*N)
147*> \endverbatim
148*>
149*> \param[out] INFO
150*> \verbatim
151*>          INFO is INTEGER
152*>         = 0:  successful exit.
153*>         < 0:  if INFO = -i, the i-th argument had an illegal value.
154*>         > 0:  The algorithm failed to compute a singular value while
155*>               working on the submatrix lying in rows and columns
156*>               INFO/(N+1) through MOD(INFO,N+1).
157*> \endverbatim
158*
159*  Authors:
160*  ========
161*
162*> \author Univ. of Tennessee
163*> \author Univ. of California Berkeley
164*> \author Univ. of Colorado Denver
165*> \author NAG Ltd.
166*
167*> \date December 2016
168*
169*> \ingroup realOTHERcomputational
170*
171*> \par Contributors:
172*  ==================
173*>
174*>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
175*>       California at Berkeley, USA \n
176*>     Osni Marques, LBNL/NERSC, USA \n
177*
178*  =====================================================================
179      SUBROUTINE SLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
180     $                   RANK, WORK, IWORK, INFO )
181*
182*  -- LAPACK computational routine (version 3.7.0) --
183*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
184*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
185*     December 2016
186*
187*     .. Scalar Arguments ..
188      CHARACTER          UPLO
189      INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
190      REAL               RCOND
191*     ..
192*     .. Array Arguments ..
193      INTEGER            IWORK( * )
194      REAL               B( LDB, * ), D( * ), E( * ), WORK( * )
195*     ..
196*
197*  =====================================================================
198*
199*     .. Parameters ..
200      REAL               ZERO, ONE, TWO
201      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
202*     ..
203*     .. Local Scalars ..
204      INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
205     $                   GIVPTR, I, ICMPQ1, ICMPQ2, IWK, J, K, NLVL,
206     $                   NM1, NSIZE, NSUB, NWORK, PERM, POLES, S, SIZEI,
207     $                   SMLSZP, SQRE, ST, ST1, U, VT, Z
208      REAL               CS, EPS, ORGNRM, R, RCND, SN, TOL
209*     ..
210*     .. External Functions ..
211      INTEGER            ISAMAX
212      REAL               SLAMCH, SLANST
213      EXTERNAL           ISAMAX, SLAMCH, SLANST
214*     ..
215*     .. External Subroutines ..
216      EXTERNAL           SCOPY, SGEMM, SLACPY, SLALSA, SLARTG, SLASCL,
217     $                   SLASDA, SLASDQ, SLASET, SLASRT, SROT, XERBLA
218*     ..
219*     .. Intrinsic Functions ..
220      INTRINSIC          ABS, INT, LOG, REAL, SIGN
221*     ..
222*     .. Executable Statements ..
223*
224*     Test the input parameters.
225*
226      INFO = 0
227*
228      IF( N.LT.0 ) THEN
229         INFO = -3
230      ELSE IF( NRHS.LT.1 ) THEN
231         INFO = -4
232      ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
233         INFO = -8
234      END IF
235      IF( INFO.NE.0 ) THEN
236         CALL XERBLA( 'SLALSD', -INFO )
237         RETURN
238      END IF
239*
240      EPS = SLAMCH( 'Epsilon' )
241*
242*     Set up the tolerance.
243*
244      IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
245         RCND = EPS
246      ELSE
247         RCND = RCOND
248      END IF
249*
250      RANK = 0
251*
252*     Quick return if possible.
253*
254      IF( N.EQ.0 ) THEN
255         RETURN
256      ELSE IF( N.EQ.1 ) THEN
257         IF( D( 1 ).EQ.ZERO ) THEN
258            CALL SLASET( 'A', 1, NRHS, ZERO, ZERO, B, LDB )
259         ELSE
260            RANK = 1
261            CALL SLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
262            D( 1 ) = ABS( D( 1 ) )
263         END IF
264         RETURN
265      END IF
266*
267*     Rotate the matrix if it is lower bidiagonal.
268*
269      IF( UPLO.EQ.'L' ) THEN
270         DO 10 I = 1, N - 1
271            CALL SLARTG( D( I ), E( I ), CS, SN, R )
272            D( I ) = R
273            E( I ) = SN*D( I+1 )
274            D( I+1 ) = CS*D( I+1 )
275            IF( NRHS.EQ.1 ) THEN
276               CALL SROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
277            ELSE
278               WORK( I*2-1 ) = CS
279               WORK( I*2 ) = SN
280            END IF
281   10    CONTINUE
282         IF( NRHS.GT.1 ) THEN
283            DO 30 I = 1, NRHS
284               DO 20 J = 1, N - 1
285                  CS = WORK( J*2-1 )
286                  SN = WORK( J*2 )
287                  CALL SROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
288   20          CONTINUE
289   30       CONTINUE
290         END IF
291      END IF
292*
293*     Scale.
294*
295      NM1 = N - 1
296      ORGNRM = SLANST( 'M', N, D, E )
297      IF( ORGNRM.EQ.ZERO ) THEN
298         CALL SLASET( 'A', N, NRHS, ZERO, ZERO, B, LDB )
299         RETURN
300      END IF
301*
302      CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
303      CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
304*
305*     If N is smaller than the minimum divide size SMLSIZ, then solve
306*     the problem with another solver.
307*
308      IF( N.LE.SMLSIZ ) THEN
309         NWORK = 1 + N*N
310         CALL SLASET( 'A', N, N, ZERO, ONE, WORK, N )
311         CALL SLASDQ( 'U', 0, N, N, 0, NRHS, D, E, WORK, N, WORK, N, B,
312     $                LDB, WORK( NWORK ), INFO )
313         IF( INFO.NE.0 ) THEN
314            RETURN
315         END IF
316         TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
317         DO 40 I = 1, N
318            IF( D( I ).LE.TOL ) THEN
319               CALL SLASET( 'A', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
320            ELSE
321               CALL SLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
322     $                      LDB, INFO )
323               RANK = RANK + 1
324            END IF
325   40    CONTINUE
326         CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
327     $               WORK( NWORK ), N )
328         CALL SLACPY( 'A', N, NRHS, WORK( NWORK ), N, B, LDB )
329*
330*        Unscale.
331*
332         CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
333         CALL SLASRT( 'D', N, D, INFO )
334         CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
335*
336         RETURN
337      END IF
338*
339*     Book-keeping and setting up some constants.
340*
341      NLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
342*
343      SMLSZP = SMLSIZ + 1
344*
345      U = 1
346      VT = 1 + SMLSIZ*N
347      DIFL = VT + SMLSZP*N
348      DIFR = DIFL + NLVL*N
349      Z = DIFR + NLVL*N*2
350      C = Z + NLVL*N
351      S = C + N
352      POLES = S + N
353      GIVNUM = POLES + 2*NLVL*N
354      BX = GIVNUM + 2*NLVL*N
355      NWORK = BX + N*NRHS
356*
357      SIZEI = 1 + N
358      K = SIZEI + N
359      GIVPTR = K + N
360      PERM = GIVPTR + N
361      GIVCOL = PERM + NLVL*N
362      IWK = GIVCOL + NLVL*N*2
363*
364      ST = 1
365      SQRE = 0
366      ICMPQ1 = 1
367      ICMPQ2 = 0
368      NSUB = 0
369*
370      DO 50 I = 1, N
371         IF( ABS( D( I ) ).LT.EPS ) THEN
372            D( I ) = SIGN( EPS, D( I ) )
373         END IF
374   50 CONTINUE
375*
376      DO 60 I = 1, NM1
377         IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
378            NSUB = NSUB + 1
379            IWORK( NSUB ) = ST
380*
381*           Subproblem found. First determine its size and then
382*           apply divide and conquer on it.
383*
384            IF( I.LT.NM1 ) THEN
385*
386*              A subproblem with E(I) small for I < NM1.
387*
388               NSIZE = I - ST + 1
389               IWORK( SIZEI+NSUB-1 ) = NSIZE
390            ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
391*
392*              A subproblem with E(NM1) not too small but I = NM1.
393*
394               NSIZE = N - ST + 1
395               IWORK( SIZEI+NSUB-1 ) = NSIZE
396            ELSE
397*
398*              A subproblem with E(NM1) small. This implies an
399*              1-by-1 subproblem at D(N), which is not solved
400*              explicitly.
401*
402               NSIZE = I - ST + 1
403               IWORK( SIZEI+NSUB-1 ) = NSIZE
404               NSUB = NSUB + 1
405               IWORK( NSUB ) = N
406               IWORK( SIZEI+NSUB-1 ) = 1
407               CALL SCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
408            END IF
409            ST1 = ST - 1
410            IF( NSIZE.EQ.1 ) THEN
411*
412*              This is a 1-by-1 subproblem and is not solved
413*              explicitly.
414*
415               CALL SCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
416            ELSE IF( NSIZE.LE.SMLSIZ ) THEN
417*
418*              This is a small subproblem and is solved by SLASDQ.
419*
420               CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
421     $                      WORK( VT+ST1 ), N )
422               CALL SLASDQ( 'U', 0, NSIZE, NSIZE, 0, NRHS, D( ST ),
423     $                      E( ST ), WORK( VT+ST1 ), N, WORK( NWORK ),
424     $                      N, B( ST, 1 ), LDB, WORK( NWORK ), INFO )
425               IF( INFO.NE.0 ) THEN
426                  RETURN
427               END IF
428               CALL SLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
429     $                      WORK( BX+ST1 ), N )
430            ELSE
431*
432*              A large problem. Solve it using divide and conquer.
433*
434               CALL SLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
435     $                      E( ST ), WORK( U+ST1 ), N, WORK( VT+ST1 ),
436     $                      IWORK( K+ST1 ), WORK( DIFL+ST1 ),
437     $                      WORK( DIFR+ST1 ), WORK( Z+ST1 ),
438     $                      WORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
439     $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
440     $                      WORK( GIVNUM+ST1 ), WORK( C+ST1 ),
441     $                      WORK( S+ST1 ), WORK( NWORK ), IWORK( IWK ),
442     $                      INFO )
443               IF( INFO.NE.0 ) THEN
444                  RETURN
445               END IF
446               BXST = BX + ST1
447               CALL SLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
448     $                      LDB, WORK( BXST ), N, WORK( U+ST1 ), N,
449     $                      WORK( VT+ST1 ), IWORK( K+ST1 ),
450     $                      WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
451     $                      WORK( Z+ST1 ), WORK( POLES+ST1 ),
452     $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
453     $                      IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
454     $                      WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
455     $                      IWORK( IWK ), INFO )
456               IF( INFO.NE.0 ) THEN
457                  RETURN
458               END IF
459            END IF
460            ST = I + 1
461         END IF
462   60 CONTINUE
463*
464*     Apply the singular values and treat the tiny ones as zero.
465*
466      TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
467*
468      DO 70 I = 1, N
469*
470*        Some of the elements in D can be negative because 1-by-1
471*        subproblems were not solved explicitly.
472*
473         IF( ABS( D( I ) ).LE.TOL ) THEN
474            CALL SLASET( 'A', 1, NRHS, ZERO, ZERO, WORK( BX+I-1 ), N )
475         ELSE
476            RANK = RANK + 1
477            CALL SLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
478     $                   WORK( BX+I-1 ), N, INFO )
479         END IF
480         D( I ) = ABS( D( I ) )
481   70 CONTINUE
482*
483*     Now apply back the right singular vectors.
484*
485      ICMPQ2 = 1
486      DO 80 I = 1, NSUB
487         ST = IWORK( I )
488         ST1 = ST - 1
489         NSIZE = IWORK( SIZEI+I-1 )
490         BXST = BX + ST1
491         IF( NSIZE.EQ.1 ) THEN
492            CALL SCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
493         ELSE IF( NSIZE.LE.SMLSIZ ) THEN
494            CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
495     $                  WORK( VT+ST1 ), N, WORK( BXST ), N, ZERO,
496     $                  B( ST, 1 ), LDB )
497         ELSE
498            CALL SLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
499     $                   B( ST, 1 ), LDB, WORK( U+ST1 ), N,
500     $                   WORK( VT+ST1 ), IWORK( K+ST1 ),
501     $                   WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
502     $                   WORK( Z+ST1 ), WORK( POLES+ST1 ),
503     $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
504     $                   IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
505     $                   WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
506     $                   IWORK( IWK ), INFO )
507            IF( INFO.NE.0 ) THEN
508               RETURN
509            END IF
510         END IF
511   80 CONTINUE
512*
513*     Unscale and sort the singular values.
514*
515      CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
516      CALL SLASRT( 'D', N, D, INFO )
517      CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
518*
519      RETURN
520*
521*     End of SLALSD
522*
523      END
524