LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slalsd()

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.
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 slalsd.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  REAL RCOND
188 * ..
189 * .. Array Arguments ..
190  INTEGER IWORK( * )
191  REAL B( LDB, * ), D( * ), E( * ), WORK( * )
192 * ..
193 *
194 * =====================================================================
195 *
196 * .. Parameters ..
197  REAL ZERO, ONE, TWO
198  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
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  REAL CS, EPS, ORGNRM, R, RCND, SN, TOL
206 * ..
207 * .. External Functions ..
208  INTEGER ISAMAX
209  REAL SLAMCH, SLANST
210  EXTERNAL isamax, slamch, slanst
211 * ..
212 * .. External Subroutines ..
213  EXTERNAL scopy, sgemm, slacpy, slalsa, slartg, slascl,
215 * ..
216 * .. Intrinsic Functions ..
217  INTRINSIC abs, int, log, real, 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( 'SLALSD', -info )
234  RETURN
235  END IF
236 *
237  eps = slamch( '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 slaset( 'A', 1, nrhs, zero, zero, b, ldb )
256  ELSE
257  rank = 1
258  CALL slascl( '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 slartg( 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 srot( 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 srot( 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 = slanst( 'M', n, d, e )
294  IF( orgnrm.EQ.zero ) THEN
295  CALL slaset( 'A', n, nrhs, zero, zero, b, ldb )
296  RETURN
297  END IF
298 *
299  CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
300  CALL slascl( '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 slaset( 'A', n, n, zero, one, work, n )
308  CALL slasdq( '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( isamax( n, d, 1 ) ) )
314  DO 40 i = 1, n
315  IF( d( i ).LE.tol ) THEN
316  CALL slaset( 'A', 1, nrhs, zero, zero, b( i, 1 ), ldb )
317  ELSE
318  CALL slascl( '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 sgemm( 'T', 'N', n, nrhs, n, one, work, n, b, ldb, zero,
324  $ work( nwork ), n )
325  CALL slacpy( 'A', n, nrhs, work( nwork ), n, b, ldb )
326 *
327 * Unscale.
328 *
329  CALL slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
330  CALL slasrt( 'D', n, d, info )
331  CALL slascl( '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( real( n ) / real( 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 scopy( 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 scopy( 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 SLASDQ.
416 *
417  CALL slaset( 'A', nsize, nsize, zero, one,
418  $ work( vt+st1 ), n )
419  CALL slasdq( '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 slacpy( '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 slasda( 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 slalsa( 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( isamax( 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 slaset( 'A', 1, nrhs, zero, zero, work( bx+i-1 ), n )
472  ELSE
473  rank = rank + 1
474  CALL slascl( '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 scopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
490  ELSE IF( nsize.LE.smlsiz ) THEN
491  CALL sgemm( 'T', 'N', nsize, nrhs, nsize, one,
492  $ work( vt+st1 ), n, work( bxst ), n, zero,
493  $ b( st, 1 ), ldb )
494  ELSE
495  CALL slalsa( 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 slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
513  CALL slasrt( 'D', n, d, info )
514  CALL slascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
515 *
516  RETURN
517 *
518 * End of SLALSD
519 *
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:143
real function slanst(NORM, N, D, E)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: slanst.f:100
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:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f90:113
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:211
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:273
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:88
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:267
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: