LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dlalsd()

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.
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 177 of file dlalsd.f.

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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
double precision function dlanst(NORM, N, D, E)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: dlanst.f:100
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
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
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
Here is the call graph for this function:
Here is the caller graph for this function: