LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.
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 171 of file dlalsd.f.

173*
174* -- LAPACK computational routine --
175* -- LAPACK is a software package provided by Univ. of Tennessee, --
176* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
177*
178* .. Scalar Arguments ..
179 CHARACTER UPLO
180 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
181 DOUBLE PRECISION RCOND
182* ..
183* .. Array Arguments ..
184 INTEGER IWORK( * )
185 DOUBLE PRECISION B( LDB, * ), D( * ), E( * ), WORK( * )
186* ..
187*
188* =====================================================================
189*
190* .. Parameters ..
191 DOUBLE PRECISION ZERO, ONE, TWO
192 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
193* ..
194* .. Local Scalars ..
195 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
196 $ GIVPTR, I, ICMPQ1, ICMPQ2, IWK, J, K, NLVL,
197 $ NM1, NSIZE, NSUB, NWORK, PERM, POLES, S, SIZEI,
198 $ SMLSZP, SQRE, ST, ST1, U, VT, Z
199 DOUBLE PRECISION CS, EPS, ORGNRM, R, RCND, SN, TOL
200* ..
201* .. External Functions ..
202 INTEGER IDAMAX
203 DOUBLE PRECISION DLAMCH, DLANST
204 EXTERNAL idamax, dlamch, dlanst
205* ..
206* .. External Subroutines ..
207 EXTERNAL dcopy, dgemm, dlacpy, dlalsa, dlartg, dlascl,
208 $ dlasda, dlasdq, dlaset, dlasrt, drot, xerbla
209* ..
210* .. Intrinsic Functions ..
211 INTRINSIC abs, dble, int, log, sign
212* ..
213* .. Executable Statements ..
214*
215* Test the input parameters.
216*
217 info = 0
218*
219 IF( n.LT.0 ) THEN
220 info = -3
221 ELSE IF( nrhs.LT.1 ) THEN
222 info = -4
223 ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) ) THEN
224 info = -8
225 END IF
226 IF( info.NE.0 ) THEN
227 CALL xerbla( 'DLALSD', -info )
228 RETURN
229 END IF
230*
231 eps = dlamch( 'Epsilon' )
232*
233* Set up the tolerance.
234*
235 IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) ) THEN
236 rcnd = eps
237 ELSE
238 rcnd = rcond
239 END IF
240*
241 rank = 0
242*
243* Quick return if possible.
244*
245 IF( n.EQ.0 ) THEN
246 RETURN
247 ELSE IF( n.EQ.1 ) THEN
248 IF( d( 1 ).EQ.zero ) THEN
249 CALL dlaset( 'A', 1, nrhs, zero, zero, b, ldb )
250 ELSE
251 rank = 1
252 CALL dlascl( 'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
253 d( 1 ) = abs( d( 1 ) )
254 END IF
255 RETURN
256 END IF
257*
258* Rotate the matrix if it is lower bidiagonal.
259*
260 IF( uplo.EQ.'L' ) THEN
261 DO 10 i = 1, n - 1
262 CALL dlartg( d( i ), e( i ), cs, sn, r )
263 d( i ) = r
264 e( i ) = sn*d( i+1 )
265 d( i+1 ) = cs*d( i+1 )
266 IF( nrhs.EQ.1 ) THEN
267 CALL drot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
268 ELSE
269 work( i*2-1 ) = cs
270 work( i*2 ) = sn
271 END IF
272 10 CONTINUE
273 IF( nrhs.GT.1 ) THEN
274 DO 30 i = 1, nrhs
275 DO 20 j = 1, n - 1
276 cs = work( j*2-1 )
277 sn = work( j*2 )
278 CALL drot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
279 20 CONTINUE
280 30 CONTINUE
281 END IF
282 END IF
283*
284* Scale.
285*
286 nm1 = n - 1
287 orgnrm = dlanst( 'M', n, d, e )
288 IF( orgnrm.EQ.zero ) THEN
289 CALL dlaset( 'A', n, nrhs, zero, zero, b, ldb )
290 RETURN
291 END IF
292*
293 CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
294 CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
295*
296* If N is smaller than the minimum divide size SMLSIZ, then solve
297* the problem with another solver.
298*
299 IF( n.LE.smlsiz ) THEN
300 nwork = 1 + n*n
301 CALL dlaset( 'A', n, n, zero, one, work, n )
302 CALL dlasdq( 'U', 0, n, n, 0, nrhs, d, e, work, n, work, n, b,
303 $ ldb, work( nwork ), info )
304 IF( info.NE.0 ) THEN
305 RETURN
306 END IF
307 tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
308 DO 40 i = 1, n
309 IF( d( i ).LE.tol ) THEN
310 CALL dlaset( 'A', 1, nrhs, zero, zero, b( i, 1 ), ldb )
311 ELSE
312 CALL dlascl( 'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
313 $ ldb, info )
314 rank = rank + 1
315 END IF
316 40 CONTINUE
317 CALL dgemm( 'T', 'N', n, nrhs, n, one, work, n, b, ldb, zero,
318 $ work( nwork ), n )
319 CALL dlacpy( 'A', n, nrhs, work( nwork ), n, b, ldb )
320*
321* Unscale.
322*
323 CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
324 CALL dlasrt( 'D', n, d, info )
325 CALL dlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
326*
327 RETURN
328 END IF
329*
330* Book-keeping and setting up some constants.
331*
332 nlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
333*
334 smlszp = smlsiz + 1
335*
336 u = 1
337 vt = 1 + smlsiz*n
338 difl = vt + smlszp*n
339 difr = difl + nlvl*n
340 z = difr + nlvl*n*2
341 c = z + nlvl*n
342 s = c + n
343 poles = s + n
344 givnum = poles + 2*nlvl*n
345 bx = givnum + 2*nlvl*n
346 nwork = bx + n*nrhs
347*
348 sizei = 1 + n
349 k = sizei + n
350 givptr = k + n
351 perm = givptr + n
352 givcol = perm + nlvl*n
353 iwk = givcol + nlvl*n*2
354*
355 st = 1
356 sqre = 0
357 icmpq1 = 1
358 icmpq2 = 0
359 nsub = 0
360*
361 DO 50 i = 1, n
362 IF( abs( d( i ) ).LT.eps ) THEN
363 d( i ) = sign( eps, d( i ) )
364 END IF
365 50 CONTINUE
366*
367 DO 60 i = 1, nm1
368 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
369 nsub = nsub + 1
370 iwork( nsub ) = st
371*
372* Subproblem found. First determine its size and then
373* apply divide and conquer on it.
374*
375 IF( i.LT.nm1 ) THEN
376*
377* A subproblem with E(I) small for I < NM1.
378*
379 nsize = i - st + 1
380 iwork( sizei+nsub-1 ) = nsize
381 ELSE IF( abs( e( i ) ).GE.eps ) THEN
382*
383* A subproblem with E(NM1) not too small but I = NM1.
384*
385 nsize = n - st + 1
386 iwork( sizei+nsub-1 ) = nsize
387 ELSE
388*
389* A subproblem with E(NM1) small. This implies an
390* 1-by-1 subproblem at D(N), which is not solved
391* explicitly.
392*
393 nsize = i - st + 1
394 iwork( sizei+nsub-1 ) = nsize
395 nsub = nsub + 1
396 iwork( nsub ) = n
397 iwork( sizei+nsub-1 ) = 1
398 CALL dcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
399 END IF
400 st1 = st - 1
401 IF( nsize.EQ.1 ) THEN
402*
403* This is a 1-by-1 subproblem and is not solved
404* explicitly.
405*
406 CALL dcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
407 ELSE IF( nsize.LE.smlsiz ) THEN
408*
409* This is a small subproblem and is solved by DLASDQ.
410*
411 CALL dlaset( 'A', nsize, nsize, zero, one,
412 $ work( vt+st1 ), n )
413 CALL dlasdq( 'U', 0, nsize, nsize, 0, nrhs, d( st ),
414 $ e( st ), work( vt+st1 ), n, work( nwork ),
415 $ n, b( st, 1 ), ldb, work( nwork ), info )
416 IF( info.NE.0 ) THEN
417 RETURN
418 END IF
419 CALL dlacpy( 'A', nsize, nrhs, b( st, 1 ), ldb,
420 $ work( bx+st1 ), n )
421 ELSE
422*
423* A large problem. Solve it using divide and conquer.
424*
425 CALL dlasda( icmpq1, smlsiz, nsize, sqre, d( st ),
426 $ e( st ), work( u+st1 ), n, work( vt+st1 ),
427 $ iwork( k+st1 ), work( difl+st1 ),
428 $ work( difr+st1 ), work( z+st1 ),
429 $ work( poles+st1 ), iwork( givptr+st1 ),
430 $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
431 $ work( givnum+st1 ), work( c+st1 ),
432 $ work( s+st1 ), work( nwork ), iwork( iwk ),
433 $ info )
434 IF( info.NE.0 ) THEN
435 RETURN
436 END IF
437 bxst = bx + st1
438 CALL dlalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
439 $ ldb, work( bxst ), n, work( u+st1 ), n,
440 $ work( vt+st1 ), iwork( k+st1 ),
441 $ work( difl+st1 ), work( difr+st1 ),
442 $ work( z+st1 ), work( poles+st1 ),
443 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
444 $ iwork( perm+st1 ), work( givnum+st1 ),
445 $ work( c+st1 ), work( s+st1 ), work( nwork ),
446 $ iwork( iwk ), info )
447 IF( info.NE.0 ) THEN
448 RETURN
449 END IF
450 END IF
451 st = i + 1
452 END IF
453 60 CONTINUE
454*
455* Apply the singular values and treat the tiny ones as zero.
456*
457 tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
458*
459 DO 70 i = 1, n
460*
461* Some of the elements in D can be negative because 1-by-1
462* subproblems were not solved explicitly.
463*
464 IF( abs( d( i ) ).LE.tol ) THEN
465 CALL dlaset( 'A', 1, nrhs, zero, zero, work( bx+i-1 ), n )
466 ELSE
467 rank = rank + 1
468 CALL dlascl( 'G', 0, 0, d( i ), one, 1, nrhs,
469 $ work( bx+i-1 ), n, info )
470 END IF
471 d( i ) = abs( d( i ) )
472 70 CONTINUE
473*
474* Now apply back the right singular vectors.
475*
476 icmpq2 = 1
477 DO 80 i = 1, nsub
478 st = iwork( i )
479 st1 = st - 1
480 nsize = iwork( sizei+i-1 )
481 bxst = bx + st1
482 IF( nsize.EQ.1 ) THEN
483 CALL dcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
484 ELSE IF( nsize.LE.smlsiz ) THEN
485 CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
486 $ work( vt+st1 ), n, work( bxst ), n, zero,
487 $ b( st, 1 ), ldb )
488 ELSE
489 CALL dlalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
490 $ b( st, 1 ), ldb, work( u+st1 ), n,
491 $ work( vt+st1 ), iwork( k+st1 ),
492 $ work( difl+st1 ), work( difr+st1 ),
493 $ work( z+st1 ), work( poles+st1 ),
494 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
495 $ iwork( perm+st1 ), work( givnum+st1 ),
496 $ work( c+st1 ), work( s+st1 ), work( nwork ),
497 $ iwork( iwk ), info )
498 IF( info.NE.0 ) THEN
499 RETURN
500 END IF
501 END IF
502 80 CONTINUE
503*
504* Unscale and sort the singular values.
505*
506 CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
507 CALL dlasrt( 'D', n, d, info )
508 CALL dlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
509*
510 RETURN
511*
512* End of DLALSD
513*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dlacpy
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
dlalsa
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
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanst
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
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlascl
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
dlasda
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
dlasdq
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
dlaset
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
dlasrt
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:88
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: