LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
slalsd.f
Go to the documentation of this file.
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 September 2012
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.4.2) --
183 * -- LAPACK is a software package provided by Univ. of Tennessee, --
184 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
185 * September 2012
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,
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
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
subroutine slalsa(ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
SLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
Definition: slalsa.f:271
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine slalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, IWORK, INFO)
SLALSD uses the singular value decomposition of A to solve the least squares problem.
Definition: slalsd.f:181
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:53
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:90
subroutine slasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e...
Definition: slasdq.f:213
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine slasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
Definition: slasda.f:275