LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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