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

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

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

Purpose: 
 SLALSD 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 REAL 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 REAL array, dimension (N-1)
         Contains the super-diagonal entries of the bidiagonal matrix.
         On exit, E has been destroyed.
[in,out]B
          B is REAL 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 REAL
         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 REAL 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 slalsd.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  REAL rcond
191 * ..
192 * .. Array Arguments ..
193  INTEGER iwork( * )
194  REAL b( ldb, * ), d( * ), e( * ), work( * )
195 * ..
196 *
197 * =====================================================================
198 *
199 * .. Parameters ..
200  REAL zero, one, two
201  parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
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  REAL cs, eps, orgnrm, r, rcnd, sn, tol
209 * ..
210 * .. External Functions ..
211  INTEGER isamax
212  REAL slamch, slanst
213  EXTERNAL isamax, slamch, slanst
214 * ..
215 * .. External Subroutines ..
216  EXTERNAL scopy, sgemm, slacpy, slalsa, slartg, slascl,
217  $ slasda, slasdq, slaset, slasrt, srot, xerbla
218 * ..
219 * .. Intrinsic Functions ..
220  INTRINSIC abs, int, log, REAL, 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( 'SLALSD', -info )
237  RETURN
238  END IF
239 *
240  eps = slamch( '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 slaset( 'A', 1, nrhs, zero, zero, b, ldb )
259  ELSE
260  rank = 1
261  CALL slascl( '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 slartg( 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 srot( 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 srot( 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 = slanst( 'M', n, d, e )
297  IF( orgnrm.EQ.zero ) THEN
298  CALL slaset( 'A', n, nrhs, zero, zero, b, ldb )
299  RETURN
300  END IF
301 *
302  CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
303  CALL slascl( '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 slaset( 'A', n, n, zero, one, work, n )
311  CALL slasdq( '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( isamax( n, d, 1 ) ) )
317  DO 40 i = 1, n
318  IF( d( i ).LE.tol ) THEN
319  CALL slaset( 'A', 1, nrhs, zero, zero, b( i, 1 ), ldb )
320  ELSE
321  CALL slascl( '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 sgemm( 'T', 'N', n, nrhs, n, one, work, n, b, ldb, zero,
327  $ work( nwork ), n )
328  CALL slacpy( 'A', n, nrhs, work( nwork ), n, b, ldb )
329 *
330 * Unscale.
331 *
332  CALL slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
333  CALL slasrt( 'D', n, d, info )
334  CALL slascl( '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( REAL( N ) / REAL( 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 scopy( 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 scopy( 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 SLASDQ.
419 *
420  CALL slaset( 'A', nsize, nsize, zero, one,
421  $ work( vt+st1 ), n )
422  CALL slasdq( '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 slacpy( '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 slasda( 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 slalsa( 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( isamax( 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 slaset( 'A', 1, nrhs, zero, zero, work( bx+i-1 ), n )
475  ELSE
476  rank = rank + 1
477  CALL slascl( '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 scopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
493  ELSE IF( nsize.LE.smlsiz ) THEN
494  CALL sgemm( 'T', 'N', nsize, nrhs, nsize, one,
495  $ work( vt+st1 ), n, work( bxst ), n, zero,
496  $ b( st, 1 ), ldb )
497  ELSE
498  CALL slalsa( 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 slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
516  CALL slasrt( 'D', n, d, info )
517  CALL slascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
518 *
519  RETURN
520 *
521 * End of SLALSD
522 *
slanst
real function slanst(NORM, N, D, E)
SLANST 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: slanst.f:102
isamax
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
slascl
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
slalsa
subroutine slalsa(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)
SLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
Definition: slalsa.f:271
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
slacpy
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
srot
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:53
slartg
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
slaset
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
slasrt
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:90
slasdq
subroutine slasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e...
Definition: slasdq.f:213
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
slasda
subroutine slasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
Definition: slasda.f:275

Here is the call graph for this function:

Here is the caller graph for this function: