LAPACK  3.10.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 *> \ingroup doubleOTHERcomputational
168 *
169 *> \par Contributors:
170 * ==================
171 *>
172 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
173 *> California at Berkeley, USA \n
174 *> Osni Marques, LBNL/NERSC, USA \n
175 *
176 * =====================================================================
177  SUBROUTINE dlalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
178  $ RANK, WORK, IWORK, INFO )
179 *
180 * -- LAPACK computational routine --
181 * -- LAPACK is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 *
184 * .. Scalar Arguments ..
185  CHARACTER UPLO
186  INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
187  DOUBLE PRECISION RCOND
188 * ..
189 * .. Array Arguments ..
190  INTEGER IWORK( * )
191  DOUBLE PRECISION B( LDB, * ), D( * ), E( * ), WORK( * )
192 * ..
193 *
194 * =====================================================================
195 *
196 * .. Parameters ..
197  DOUBLE PRECISION ZERO, ONE, TWO
198  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
199 * ..
200 * .. Local Scalars ..
201  INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
202  $ givptr, i, icmpq1, icmpq2, iwk, j, k, nlvl,
203  $ nm1, nsize, nsub, nwork, perm, poles, s, sizei,
204  $ smlszp, sqre, st, st1, u, vt, z
205  DOUBLE PRECISION CS, EPS, ORGNRM, R, RCND, SN, TOL
206 * ..
207 * .. External Functions ..
208  INTEGER IDAMAX
209  DOUBLE PRECISION DLAMCH, DLANST
210  EXTERNAL idamax, dlamch, dlanst
211 * ..
212 * .. External Subroutines ..
213  EXTERNAL dcopy, dgemm, dlacpy, dlalsa, dlartg, dlascl,
215 * ..
216 * .. Intrinsic Functions ..
217  INTRINSIC abs, dble, int, log, sign
218 * ..
219 * .. Executable Statements ..
220 *
221 * Test the input parameters.
222 *
223  info = 0
224 *
225  IF( n.LT.0 ) THEN
226  info = -3
227  ELSE IF( nrhs.LT.1 ) THEN
228  info = -4
229  ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) ) THEN
230  info = -8
231  END IF
232  IF( info.NE.0 ) THEN
233  CALL xerbla( 'DLALSD', -info )
234  RETURN
235  END IF
236 *
237  eps = dlamch( 'Epsilon' )
238 *
239 * Set up the tolerance.
240 *
241  IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) ) THEN
242  rcnd = eps
243  ELSE
244  rcnd = rcond
245  END IF
246 *
247  rank = 0
248 *
249 * Quick return if possible.
250 *
251  IF( n.EQ.0 ) THEN
252  RETURN
253  ELSE IF( n.EQ.1 ) THEN
254  IF( d( 1 ).EQ.zero ) THEN
255  CALL dlaset( 'A', 1, nrhs, zero, zero, b, ldb )
256  ELSE
257  rank = 1
258  CALL dlascl( 'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
259  d( 1 ) = abs( d( 1 ) )
260  END IF
261  RETURN
262  END IF
263 *
264 * Rotate the matrix if it is lower bidiagonal.
265 *
266  IF( uplo.EQ.'L' ) THEN
267  DO 10 i = 1, n - 1
268  CALL dlartg( d( i ), e( i ), cs, sn, r )
269  d( i ) = r
270  e( i ) = sn*d( i+1 )
271  d( i+1 ) = cs*d( i+1 )
272  IF( nrhs.EQ.1 ) THEN
273  CALL drot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
274  ELSE
275  work( i*2-1 ) = cs
276  work( i*2 ) = sn
277  END IF
278  10 CONTINUE
279  IF( nrhs.GT.1 ) THEN
280  DO 30 i = 1, nrhs
281  DO 20 j = 1, n - 1
282  cs = work( j*2-1 )
283  sn = work( j*2 )
284  CALL drot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
285  20 CONTINUE
286  30 CONTINUE
287  END IF
288  END IF
289 *
290 * Scale.
291 *
292  nm1 = n - 1
293  orgnrm = dlanst( 'M', n, d, e )
294  IF( orgnrm.EQ.zero ) THEN
295  CALL dlaset( 'A', n, nrhs, zero, zero, b, ldb )
296  RETURN
297  END IF
298 *
299  CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
300  CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
301 *
302 * If N is smaller than the minimum divide size SMLSIZ, then solve
303 * the problem with another solver.
304 *
305  IF( n.LE.smlsiz ) THEN
306  nwork = 1 + n*n
307  CALL dlaset( 'A', n, n, zero, one, work, n )
308  CALL dlasdq( 'U', 0, n, n, 0, nrhs, d, e, work, n, work, n, b,
309  $ ldb, work( nwork ), info )
310  IF( info.NE.0 ) THEN
311  RETURN
312  END IF
313  tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
314  DO 40 i = 1, n
315  IF( d( i ).LE.tol ) THEN
316  CALL dlaset( 'A', 1, nrhs, zero, zero, b( i, 1 ), ldb )
317  ELSE
318  CALL dlascl( 'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
319  $ ldb, info )
320  rank = rank + 1
321  END IF
322  40 CONTINUE
323  CALL dgemm( 'T', 'N', n, nrhs, n, one, work, n, b, ldb, zero,
324  $ work( nwork ), n )
325  CALL dlacpy( 'A', n, nrhs, work( nwork ), n, b, ldb )
326 *
327 * Unscale.
328 *
329  CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
330  CALL dlasrt( 'D', n, d, info )
331  CALL dlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
332 *
333  RETURN
334  END IF
335 *
336 * Book-keeping and setting up some constants.
337 *
338  nlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
339 *
340  smlszp = smlsiz + 1
341 *
342  u = 1
343  vt = 1 + smlsiz*n
344  difl = vt + smlszp*n
345  difr = difl + nlvl*n
346  z = difr + nlvl*n*2
347  c = z + nlvl*n
348  s = c + n
349  poles = s + n
350  givnum = poles + 2*nlvl*n
351  bx = givnum + 2*nlvl*n
352  nwork = bx + n*nrhs
353 *
354  sizei = 1 + n
355  k = sizei + n
356  givptr = k + n
357  perm = givptr + n
358  givcol = perm + nlvl*n
359  iwk = givcol + nlvl*n*2
360 *
361  st = 1
362  sqre = 0
363  icmpq1 = 1
364  icmpq2 = 0
365  nsub = 0
366 *
367  DO 50 i = 1, n
368  IF( abs( d( i ) ).LT.eps ) THEN
369  d( i ) = sign( eps, d( i ) )
370  END IF
371  50 CONTINUE
372 *
373  DO 60 i = 1, nm1
374  IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
375  nsub = nsub + 1
376  iwork( nsub ) = st
377 *
378 * Subproblem found. First determine its size and then
379 * apply divide and conquer on it.
380 *
381  IF( i.LT.nm1 ) THEN
382 *
383 * A subproblem with E(I) small for I < NM1.
384 *
385  nsize = i - st + 1
386  iwork( sizei+nsub-1 ) = nsize
387  ELSE IF( abs( e( i ) ).GE.eps ) THEN
388 *
389 * A subproblem with E(NM1) not too small but I = NM1.
390 *
391  nsize = n - st + 1
392  iwork( sizei+nsub-1 ) = nsize
393  ELSE
394 *
395 * A subproblem with E(NM1) small. This implies an
396 * 1-by-1 subproblem at D(N), which is not solved
397 * explicitly.
398 *
399  nsize = i - st + 1
400  iwork( sizei+nsub-1 ) = nsize
401  nsub = nsub + 1
402  iwork( nsub ) = n
403  iwork( sizei+nsub-1 ) = 1
404  CALL dcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
405  END IF
406  st1 = st - 1
407  IF( nsize.EQ.1 ) THEN
408 *
409 * This is a 1-by-1 subproblem and is not solved
410 * explicitly.
411 *
412  CALL dcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
413  ELSE IF( nsize.LE.smlsiz ) THEN
414 *
415 * This is a small subproblem and is solved by DLASDQ.
416 *
417  CALL dlaset( 'A', nsize, nsize, zero, one,
418  $ work( vt+st1 ), n )
419  CALL dlasdq( 'U', 0, nsize, nsize, 0, nrhs, d( st ),
420  $ e( st ), work( vt+st1 ), n, work( nwork ),
421  $ n, b( st, 1 ), ldb, work( nwork ), info )
422  IF( info.NE.0 ) THEN
423  RETURN
424  END IF
425  CALL dlacpy( 'A', nsize, nrhs, b( st, 1 ), ldb,
426  $ work( bx+st1 ), n )
427  ELSE
428 *
429 * A large problem. Solve it using divide and conquer.
430 *
431  CALL dlasda( icmpq1, smlsiz, nsize, sqre, d( st ),
432  $ e( st ), work( u+st1 ), n, work( vt+st1 ),
433  $ iwork( k+st1 ), work( difl+st1 ),
434  $ work( difr+st1 ), work( z+st1 ),
435  $ work( poles+st1 ), iwork( givptr+st1 ),
436  $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
437  $ work( givnum+st1 ), work( c+st1 ),
438  $ work( s+st1 ), work( nwork ), iwork( iwk ),
439  $ info )
440  IF( info.NE.0 ) THEN
441  RETURN
442  END IF
443  bxst = bx + st1
444  CALL dlalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
445  $ ldb, work( bxst ), n, work( u+st1 ), n,
446  $ work( vt+st1 ), iwork( k+st1 ),
447  $ work( difl+st1 ), work( difr+st1 ),
448  $ work( z+st1 ), work( poles+st1 ),
449  $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
450  $ iwork( perm+st1 ), work( givnum+st1 ),
451  $ work( c+st1 ), work( s+st1 ), work( nwork ),
452  $ iwork( iwk ), info )
453  IF( info.NE.0 ) THEN
454  RETURN
455  END IF
456  END IF
457  st = i + 1
458  END IF
459  60 CONTINUE
460 *
461 * Apply the singular values and treat the tiny ones as zero.
462 *
463  tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
464 *
465  DO 70 i = 1, n
466 *
467 * Some of the elements in D can be negative because 1-by-1
468 * subproblems were not solved explicitly.
469 *
470  IF( abs( d( i ) ).LE.tol ) THEN
471  CALL dlaset( 'A', 1, nrhs, zero, zero, work( bx+i-1 ), n )
472  ELSE
473  rank = rank + 1
474  CALL dlascl( 'G', 0, 0, d( i ), one, 1, nrhs,
475  $ work( bx+i-1 ), n, info )
476  END IF
477  d( i ) = abs( d( i ) )
478  70 CONTINUE
479 *
480 * Now apply back the right singular vectors.
481 *
482  icmpq2 = 1
483  DO 80 i = 1, nsub
484  st = iwork( i )
485  st1 = st - 1
486  nsize = iwork( sizei+i-1 )
487  bxst = bx + st1
488  IF( nsize.EQ.1 ) THEN
489  CALL dcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
490  ELSE IF( nsize.LE.smlsiz ) THEN
491  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
492  $ work( vt+st1 ), n, work( bxst ), n, zero,
493  $ b( st, 1 ), ldb )
494  ELSE
495  CALL dlalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
496  $ b( st, 1 ), ldb, work( u+st1 ), n,
497  $ work( vt+st1 ), iwork( k+st1 ),
498  $ work( difl+st1 ), work( difr+st1 ),
499  $ work( z+st1 ), work( poles+st1 ),
500  $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
501  $ iwork( perm+st1 ), work( givnum+st1 ),
502  $ work( c+st1 ), work( s+st1 ), work( nwork ),
503  $ iwork( iwk ), info )
504  IF( info.NE.0 ) THEN
505  RETURN
506  END IF
507  END IF
508  80 CONTINUE
509 *
510 * Unscale and sort the singular values.
511 *
512  CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
513  CALL dlasrt( 'D', n, d, info )
514  CALL dlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
515 *
516  RETURN
517 *
518 * End of DLALSD
519 *
520  END
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:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
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:110
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:273
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:211
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:88
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
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:267
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:179