LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlalsd ( character  UPLO,
integer  SMLSIZ,
integer  N,
integer  NRHS,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision  RCOND,
integer  RANK,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DLALSD uses the singular value decomposition of A to solve the least squares problem.

Download DLALSD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DLALSD uses the singular value decomposition of A to solve the least
 squares problem of finding X to minimize the Euclidean norm of each
 column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
 are N-by-NRHS. The solution X overwrites B.

 The singular values of A smaller than RCOND times the largest
 singular value are treated as zero in solving the least squares
 problem; in this case a minimum norm solution is returned.
 The actual singular values are returned in D in ascending order.

 This code makes very mild assumptions about floating point
 arithmetic. It will work on machines with a guard digit in
 add/subtract, or on those binary machines without guard digits
 which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
 It could conceivably fail on hexadecimal or decimal machines
 without guard digits, but we know of none.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
         = 'U': D and E define an upper bidiagonal matrix.
         = 'L': D and E define a  lower bidiagonal matrix.
[in]SMLSIZ
          SMLSIZ is INTEGER
         The maximum size of the subproblems at the bottom of the
         computation tree.
[in]N
          N is INTEGER
         The dimension of the  bidiagonal matrix.  N >= 0.
[in]NRHS
          NRHS is INTEGER
         The number of columns of B. NRHS must be at least 1.
[in,out]D
          D is DOUBLE PRECISION array, dimension (N)
         On entry D contains the main diagonal of the bidiagonal
         matrix. On exit, if INFO = 0, D contains its singular values.
[in,out]E
          E is DOUBLE PRECISION array, dimension (N-1)
         Contains the super-diagonal entries of the bidiagonal matrix.
         On exit, E has been destroyed.
[in,out]B
          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
         On input, B contains the right hand sides of the least
         squares problem. On output, B contains the solution X.
[in]LDB
          LDB is INTEGER
         The leading dimension of B in the calling subprogram.
         LDB must be at least max(1,N).
[in]RCOND
          RCOND is DOUBLE PRECISION
         The singular values of A less than or equal to RCOND times
         the largest singular value are treated as zero in solving
         the least squares problem. If RCOND is negative,
         machine precision is used instead.
         For example, if diag(S)*X=B were the least squares problem,
         where diag(S) is a diagonal matrix of singular values, the
         solution would be X(i) = B(i) / S(i) if S(i) is greater than
         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
         RCOND*max(S).
[out]RANK
          RANK is INTEGER
         The number of singular values of A greater than RCOND times
         the largest singular value.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension at least
         (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
         where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).
[out]IWORK
          IWORK is INTEGER array, dimension at least
         (3*N*NLVL + 11*N)
[out]INFO
          INFO is INTEGER
         = 0:  successful exit.
         < 0:  if INFO = -i, the i-th argument had an illegal value.
         > 0:  The algorithm failed to compute a singular value while
               working on the submatrix lying in rows and columns
               INFO/(N+1) through MOD(INFO,N+1).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 181 of file dlalsd.f.

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 *
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
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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 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
double precision function dlanst(NORM, N, D, E)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.
Definition: dlanst.f:102
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99

Here is the call graph for this function:

Here is the caller graph for this function: