LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
dlalsd.f
Go to the documentation of this file.
1 *> \brief \b DLALSD 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 DLALSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlalsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlalsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlalsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLALSD( 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 * DOUBLE PRECISION RCOND
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * )
31 * DOUBLE PRECISION B( LDB, * ), D( * ), E( * ), WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> DLALSD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 doubleOTHERcomputational
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 dlalsd( 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  DOUBLE PRECISION RCOND
191 * ..
192 * .. Array Arguments ..
193  INTEGER IWORK( * )
194  DOUBLE PRECISION B( ldb, * ), D( * ), E( * ), WORK( * )
195 * ..
196 *
197 * =====================================================================
198 *
199 * .. Parameters ..
200  DOUBLE PRECISION ZERO, ONE, TWO
201  parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
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  DOUBLE PRECISION CS, EPS, ORGNRM, R, RCND, SN, TOL
209 * ..
210 * .. External Functions ..
211  INTEGER IDAMAX
212  DOUBLE PRECISION DLAMCH, DLANST
213  EXTERNAL idamax, dlamch, dlanst
214 * ..
215 * .. External Subroutines ..
216  EXTERNAL dcopy, dgemm, dlacpy, dlalsa, dlartg, dlascl,
218 * ..
219 * .. Intrinsic Functions ..
220  INTRINSIC abs, dble, int, log, 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( 'DLALSD', -info )
237  RETURN
238  END IF
239 *
240  eps = dlamch( '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 dlaset( 'A', 1, nrhs, zero, zero, b, ldb )
259  ELSE
260  rank = 1
261  CALL dlascl( '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 dlartg( 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 drot( 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 drot( 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 = dlanst( 'M', n, d, e )
297  IF( orgnrm.EQ.zero ) THEN
298  CALL dlaset( 'A', n, nrhs, zero, zero, b, ldb )
299  RETURN
300  END IF
301 *
302  CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
303  CALL dlascl( '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 dlaset( 'A', n, n, zero, one, work, n )
311  CALL dlasdq( '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( idamax( n, d, 1 ) ) )
317  DO 40 i = 1, n
318  IF( d( i ).LE.tol ) THEN
319  CALL dlaset( 'A', 1, nrhs, zero, zero, b( i, 1 ), ldb )
320  ELSE
321  CALL dlascl( '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 dgemm( 'T', 'N', n, nrhs, n, one, work, n, b, ldb, zero,
327  $ work( nwork ), n )
328  CALL dlacpy( 'A', n, nrhs, work( nwork ), n, b, ldb )
329 *
330 * Unscale.
331 *
332  CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
333  CALL dlasrt( 'D', n, d, info )
334  CALL dlascl( '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( dble( n ) / dble( 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 dcopy( 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 dcopy( 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 DLASDQ.
419 *
420  CALL dlaset( 'A', nsize, nsize, zero, one,
421  $ work( vt+st1 ), n )
422  CALL dlasdq( '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 dlacpy( '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 dlasda( 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 dlalsa( 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( idamax( 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 dlaset( 'A', 1, nrhs, zero, zero, work( bx+i-1 ), n )
475  ELSE
476  rank = rank + 1
477  CALL dlascl( '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 dcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
493  ELSE IF( nsize.LE.smlsiz ) THEN
494  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
495  $ work( vt+st1 ), n, work( bxst ), n, zero,
496  $ b( st, 1 ), ldb )
497  ELSE
498  CALL dlalsa( 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 dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
516  CALL dlasrt( 'D', n, d, info )
517  CALL dlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
518 *
519  RETURN
520 *
521 * End of DLALSD
522 *
523  END
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:90
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
subroutine dlasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
Definition: dlasda.f:276
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:145
subroutine dlalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, IWORK, INFO)
DLALSD uses the singular value decomposition of A to solve the least squares problem.
Definition: dlalsd.f:181
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlalsa(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)
DLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
Definition: dlalsa.f:271
subroutine dlasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e...
Definition: dlasdq.f:213
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99