Annotation of rpl/lapack/lapack/dgelsd.f, revision 1.9
1.9 ! bertrand 1: *> \brief <b> DGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices</b>
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DGELSD + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsd.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsd.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsd.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
! 22: * WORK, LWORK, IWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
! 26: * DOUBLE PRECISION RCOND
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * INTEGER IWORK( * )
! 30: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DGELSD computes the minimum-norm solution to a real linear least
! 40: *> squares problem:
! 41: *> minimize 2-norm(| b - A*x |)
! 42: *> using the singular value decomposition (SVD) of A. A is an M-by-N
! 43: *> matrix which may be rank-deficient.
! 44: *>
! 45: *> Several right hand side vectors b and solution vectors x can be
! 46: *> handled in a single call; they are stored as the columns of the
! 47: *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
! 48: *> matrix X.
! 49: *>
! 50: *> The problem is solved in three steps:
! 51: *> (1) Reduce the coefficient matrix A to bidiagonal form with
! 52: *> Householder transformations, reducing the original problem
! 53: *> into a "bidiagonal least squares problem" (BLS)
! 54: *> (2) Solve the BLS using a divide and conquer approach.
! 55: *> (3) Apply back all the Householder tranformations to solve
! 56: *> the original least squares problem.
! 57: *>
! 58: *> The effective rank of A is determined by treating as zero those
! 59: *> singular values which are less than RCOND times the largest singular
! 60: *> value.
! 61: *>
! 62: *> The divide and conquer algorithm makes very mild assumptions about
! 63: *> floating point arithmetic. It will work on machines with a guard
! 64: *> digit in add/subtract, or on those binary machines without guard
! 65: *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
! 66: *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
! 67: *> without guard digits, but we know of none.
! 68: *> \endverbatim
! 69: *
! 70: * Arguments:
! 71: * ==========
! 72: *
! 73: *> \param[in] M
! 74: *> \verbatim
! 75: *> M is INTEGER
! 76: *> The number of rows of A. M >= 0.
! 77: *> \endverbatim
! 78: *>
! 79: *> \param[in] N
! 80: *> \verbatim
! 81: *> N is INTEGER
! 82: *> The number of columns of A. N >= 0.
! 83: *> \endverbatim
! 84: *>
! 85: *> \param[in] NRHS
! 86: *> \verbatim
! 87: *> NRHS is INTEGER
! 88: *> The number of right hand sides, i.e., the number of columns
! 89: *> of the matrices B and X. NRHS >= 0.
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] A
! 93: *> \verbatim
! 94: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 95: *> On entry, the M-by-N matrix A.
! 96: *> On exit, A has been destroyed.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] LDA
! 100: *> \verbatim
! 101: *> LDA is INTEGER
! 102: *> The leading dimension of the array A. LDA >= max(1,M).
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[in,out] B
! 106: *> \verbatim
! 107: *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
! 108: *> On entry, the M-by-NRHS right hand side matrix B.
! 109: *> On exit, B is overwritten by the N-by-NRHS solution
! 110: *> matrix X. If m >= n and RANK = n, the residual
! 111: *> sum-of-squares for the solution in the i-th column is given
! 112: *> by the sum of squares of elements n+1:m in that column.
! 113: *> \endverbatim
! 114: *>
! 115: *> \param[in] LDB
! 116: *> \verbatim
! 117: *> LDB is INTEGER
! 118: *> The leading dimension of the array B. LDB >= max(1,max(M,N)).
! 119: *> \endverbatim
! 120: *>
! 121: *> \param[out] S
! 122: *> \verbatim
! 123: *> S is DOUBLE PRECISION array, dimension (min(M,N))
! 124: *> The singular values of A in decreasing order.
! 125: *> The condition number of A in the 2-norm = S(1)/S(min(m,n)).
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[in] RCOND
! 129: *> \verbatim
! 130: *> RCOND is DOUBLE PRECISION
! 131: *> RCOND is used to determine the effective rank of A.
! 132: *> Singular values S(i) <= RCOND*S(1) are treated as zero.
! 133: *> If RCOND < 0, machine precision is used instead.
! 134: *> \endverbatim
! 135: *>
! 136: *> \param[out] RANK
! 137: *> \verbatim
! 138: *> RANK is INTEGER
! 139: *> The effective rank of A, i.e., the number of singular values
! 140: *> which are greater than RCOND*S(1).
! 141: *> \endverbatim
! 142: *>
! 143: *> \param[out] WORK
! 144: *> \verbatim
! 145: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 146: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 147: *> \endverbatim
! 148: *>
! 149: *> \param[in] LWORK
! 150: *> \verbatim
! 151: *> LWORK is INTEGER
! 152: *> The dimension of the array WORK. LWORK must be at least 1.
! 153: *> The exact minimum amount of workspace needed depends on M,
! 154: *> N and NRHS. As long as LWORK is at least
! 155: *> 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,
! 156: *> if M is greater than or equal to N or
! 157: *> 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,
! 158: *> if M is less than N, the code will execute correctly.
! 159: *> SMLSIZ is returned by ILAENV and is equal to the maximum
! 160: *> size of the subproblems at the bottom of the computation
! 161: *> tree (usually about 25), and
! 162: *> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
! 163: *> For good performance, LWORK should generally be larger.
! 164: *>
! 165: *> If LWORK = -1, then a workspace query is assumed; the routine
! 166: *> only calculates the optimal size of the WORK array, returns
! 167: *> this value as the first entry of the WORK array, and no error
! 168: *> message related to LWORK is issued by XERBLA.
! 169: *> \endverbatim
! 170: *>
! 171: *> \param[out] IWORK
! 172: *> \verbatim
! 173: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
! 174: *> LIWORK >= max(1, 3 * MINMN * NLVL + 11 * MINMN),
! 175: *> where MINMN = MIN( M,N ).
! 176: *> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
! 177: *> \endverbatim
! 178: *>
! 179: *> \param[out] INFO
! 180: *> \verbatim
! 181: *> INFO is INTEGER
! 182: *> = 0: successful exit
! 183: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 184: *> > 0: the algorithm for computing the SVD failed to converge;
! 185: *> if INFO = i, i off-diagonal elements of an intermediate
! 186: *> bidiagonal form did not converge to zero.
! 187: *> \endverbatim
! 188: *
! 189: * Authors:
! 190: * ========
! 191: *
! 192: *> \author Univ. of Tennessee
! 193: *> \author Univ. of California Berkeley
! 194: *> \author Univ. of Colorado Denver
! 195: *> \author NAG Ltd.
! 196: *
! 197: *> \date November 2011
! 198: *
! 199: *> \ingroup doubleGEsolve
! 200: *
! 201: *> \par Contributors:
! 202: * ==================
! 203: *>
! 204: *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
! 205: *> California at Berkeley, USA \n
! 206: *> Osni Marques, LBNL/NERSC, USA \n
! 207: *
! 208: * =====================================================================
1.1 bertrand 209: SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
210: $ WORK, LWORK, IWORK, INFO )
211: *
1.9 ! bertrand 212: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 213: * -- LAPACK is a software package provided by Univ. of Tennessee, --
214: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 215: * November 2011
1.1 bertrand 216: *
217: * .. Scalar Arguments ..
218: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
219: DOUBLE PRECISION RCOND
220: * ..
221: * .. Array Arguments ..
222: INTEGER IWORK( * )
223: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
224: * ..
225: *
226: * =====================================================================
227: *
228: * .. Parameters ..
229: DOUBLE PRECISION ZERO, ONE, TWO
230: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
231: * ..
232: * .. Local Scalars ..
233: LOGICAL LQUERY
234: INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
1.5 bertrand 235: $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK,
236: $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD
1.1 bertrand 237: DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
238: * ..
239: * .. External Subroutines ..
240: EXTERNAL DGEBRD, DGELQF, DGEQRF, DLABAD, DLACPY, DLALSD,
241: $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR, XERBLA
242: * ..
243: * .. External Functions ..
244: INTEGER ILAENV
245: DOUBLE PRECISION DLAMCH, DLANGE
246: EXTERNAL ILAENV, DLAMCH, DLANGE
247: * ..
248: * .. Intrinsic Functions ..
249: INTRINSIC DBLE, INT, LOG, MAX, MIN
250: * ..
251: * .. Executable Statements ..
252: *
253: * Test the input arguments.
254: *
255: INFO = 0
256: MINMN = MIN( M, N )
257: MAXMN = MAX( M, N )
258: MNTHR = ILAENV( 6, 'DGELSD', ' ', M, N, NRHS, -1 )
259: LQUERY = ( LWORK.EQ.-1 )
260: IF( M.LT.0 ) THEN
261: INFO = -1
262: ELSE IF( N.LT.0 ) THEN
263: INFO = -2
264: ELSE IF( NRHS.LT.0 ) THEN
265: INFO = -3
266: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
267: INFO = -5
268: ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
269: INFO = -7
270: END IF
271: *
272: SMLSIZ = ILAENV( 9, 'DGELSD', ' ', 0, 0, 0, 0 )
273: *
274: * Compute workspace.
275: * (Note: Comments in the code beginning "Workspace:" describe the
276: * minimal amount of workspace needed at that point in the code,
277: * as well as the preferred amount for good performance.
278: * NB refers to the optimal block size for the immediately
279: * following subroutine, as returned by ILAENV.)
280: *
281: MINWRK = 1
1.5 bertrand 282: LIWORK = 1
1.1 bertrand 283: MINMN = MAX( 1, MINMN )
284: NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ+1 ) ) /
285: $ LOG( TWO ) ) + 1, 0 )
286: *
287: IF( INFO.EQ.0 ) THEN
288: MAXWRK = 0
1.5 bertrand 289: LIWORK = 3*MINMN*NLVL + 11*MINMN
1.1 bertrand 290: MM = M
291: IF( M.GE.N .AND. M.GE.MNTHR ) THEN
292: *
293: * Path 1a - overdetermined, with many more rows than columns.
294: *
295: MM = N
296: MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, 'DGEQRF', ' ', M, N,
297: $ -1, -1 ) )
298: MAXWRK = MAX( MAXWRK, N+NRHS*
299: $ ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N, -1 ) )
300: END IF
301: IF( M.GE.N ) THEN
302: *
303: * Path 1 - overdetermined or exactly determined.
304: *
305: MAXWRK = MAX( MAXWRK, 3*N+( MM+N )*
306: $ ILAENV( 1, 'DGEBRD', ' ', MM, N, -1, -1 ) )
307: MAXWRK = MAX( MAXWRK, 3*N+NRHS*
308: $ ILAENV( 1, 'DORMBR', 'QLT', MM, NRHS, N, -1 ) )
309: MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
310: $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, N, -1 ) )
311: WLALSD = 9*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS+(SMLSIZ+1)**2
312: MAXWRK = MAX( MAXWRK, 3*N+WLALSD )
313: MINWRK = MAX( 3*N+MM, 3*N+NRHS, 3*N+WLALSD )
314: END IF
315: IF( N.GT.M ) THEN
316: WLALSD = 9*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS+(SMLSIZ+1)**2
317: IF( N.GE.MNTHR ) THEN
318: *
319: * Path 2a - underdetermined, with many more columns
320: * than rows.
321: *
322: MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
323: MAXWRK = MAX( MAXWRK, M*M+4*M+2*M*
324: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
325: MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS*
326: $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, M, -1 ) )
327: MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )*
328: $ ILAENV( 1, 'DORMBR', 'PLN', M, NRHS, M, -1 ) )
329: IF( NRHS.GT.1 ) THEN
330: MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS )
331: ELSE
332: MAXWRK = MAX( MAXWRK, M*M+2*M )
333: END IF
334: MAXWRK = MAX( MAXWRK, M+NRHS*
335: $ ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, -1 ) )
336: MAXWRK = MAX( MAXWRK, M*M+4*M+WLALSD )
337: ! XXX: Ensure the Path 2a case below is triggered. The workspace
338: ! calculation should use queries for all routines eventually.
339: MAXWRK = MAX( MAXWRK,
340: $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) )
341: ELSE
342: *
343: * Path 2 - remaining underdetermined cases.
344: *
345: MAXWRK = 3*M + ( N+M )*ILAENV( 1, 'DGEBRD', ' ', M, N,
346: $ -1, -1 )
347: MAXWRK = MAX( MAXWRK, 3*M+NRHS*
348: $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, N, -1 ) )
349: MAXWRK = MAX( MAXWRK, 3*M+M*
350: $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, M, -1 ) )
351: MAXWRK = MAX( MAXWRK, 3*M+WLALSD )
352: END IF
353: MINWRK = MAX( 3*M+NRHS, 3*M+M, 3*M+WLALSD )
354: END IF
355: MINWRK = MIN( MINWRK, MAXWRK )
356: WORK( 1 ) = MAXWRK
1.5 bertrand 357: IWORK( 1 ) = LIWORK
358:
1.1 bertrand 359: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
360: INFO = -12
361: END IF
362: END IF
363: *
364: IF( INFO.NE.0 ) THEN
365: CALL XERBLA( 'DGELSD', -INFO )
366: RETURN
367: ELSE IF( LQUERY ) THEN
368: GO TO 10
369: END IF
370: *
371: * Quick return if possible.
372: *
373: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
374: RANK = 0
375: RETURN
376: END IF
377: *
378: * Get machine parameters.
379: *
380: EPS = DLAMCH( 'P' )
381: SFMIN = DLAMCH( 'S' )
382: SMLNUM = SFMIN / EPS
383: BIGNUM = ONE / SMLNUM
384: CALL DLABAD( SMLNUM, BIGNUM )
385: *
386: * Scale A if max entry outside range [SMLNUM,BIGNUM].
387: *
388: ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
389: IASCL = 0
390: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
391: *
392: * Scale matrix norm up to SMLNUM.
393: *
394: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
395: IASCL = 1
396: ELSE IF( ANRM.GT.BIGNUM ) THEN
397: *
398: * Scale matrix norm down to BIGNUM.
399: *
400: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
401: IASCL = 2
402: ELSE IF( ANRM.EQ.ZERO ) THEN
403: *
404: * Matrix all zero. Return zero solution.
405: *
406: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
407: CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
408: RANK = 0
409: GO TO 10
410: END IF
411: *
412: * Scale B if max entry outside range [SMLNUM,BIGNUM].
413: *
414: BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
415: IBSCL = 0
416: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
417: *
418: * Scale matrix norm up to SMLNUM.
419: *
420: CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
421: IBSCL = 1
422: ELSE IF( BNRM.GT.BIGNUM ) THEN
423: *
424: * Scale matrix norm down to BIGNUM.
425: *
426: CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
427: IBSCL = 2
428: END IF
429: *
430: * If M < N make sure certain entries of B are zero.
431: *
432: IF( M.LT.N )
433: $ CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
434: *
435: * Overdetermined case.
436: *
437: IF( M.GE.N ) THEN
438: *
439: * Path 1 - overdetermined or exactly determined.
440: *
441: MM = M
442: IF( M.GE.MNTHR ) THEN
443: *
444: * Path 1a - overdetermined, with many more rows than columns.
445: *
446: MM = N
447: ITAU = 1
448: NWORK = ITAU + N
449: *
450: * Compute A=Q*R.
451: * (Workspace: need 2*N, prefer N+N*NB)
452: *
453: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
454: $ LWORK-NWORK+1, INFO )
455: *
456: * Multiply B by transpose(Q).
457: * (Workspace: need N+NRHS, prefer N+NRHS*NB)
458: *
459: CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B,
460: $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
461: *
462: * Zero out below R.
463: *
464: IF( N.GT.1 ) THEN
465: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
466: END IF
467: END IF
468: *
469: IE = 1
470: ITAUQ = IE + N
471: ITAUP = ITAUQ + N
472: NWORK = ITAUP + N
473: *
474: * Bidiagonalize R in A.
475: * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
476: *
477: CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
478: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
479: $ INFO )
480: *
481: * Multiply B by transpose of left bidiagonalizing vectors of R.
482: * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
483: *
484: CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
485: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
486: *
487: * Solve the bidiagonal least squares problem.
488: *
489: CALL DLALSD( 'U', SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB,
490: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
491: IF( INFO.NE.0 ) THEN
492: GO TO 10
493: END IF
494: *
495: * Multiply B by right bidiagonalizing vectors of R.
496: *
497: CALL DORMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ),
498: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
499: *
500: ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
501: $ MAX( M, 2*M-4, NRHS, N-3*M, WLALSD ) ) THEN
502: *
503: * Path 2a - underdetermined, with many more columns than rows
504: * and sufficient workspace for an efficient algorithm.
505: *
506: LDWORK = M
507: IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
508: $ M*LDA+M+M*NRHS, 4*M+M*LDA+WLALSD ) )LDWORK = LDA
509: ITAU = 1
510: NWORK = M + 1
511: *
512: * Compute A=L*Q.
513: * (Workspace: need 2*M, prefer M+M*NB)
514: *
515: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
516: $ LWORK-NWORK+1, INFO )
517: IL = NWORK
518: *
519: * Copy L to WORK(IL), zeroing out above its diagonal.
520: *
521: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
522: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ),
523: $ LDWORK )
524: IE = IL + LDWORK*M
525: ITAUQ = IE + M
526: ITAUP = ITAUQ + M
527: NWORK = ITAUP + M
528: *
529: * Bidiagonalize L in WORK(IL).
530: * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
531: *
532: CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ),
533: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
534: $ LWORK-NWORK+1, INFO )
535: *
536: * Multiply B by transpose of left bidiagonalizing vectors of L.
537: * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
538: *
539: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK,
540: $ WORK( ITAUQ ), B, LDB, WORK( NWORK ),
541: $ LWORK-NWORK+1, INFO )
542: *
543: * Solve the bidiagonal least squares problem.
544: *
545: CALL DLALSD( 'U', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB,
546: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
547: IF( INFO.NE.0 ) THEN
548: GO TO 10
549: END IF
550: *
551: * Multiply B by right bidiagonalizing vectors of L.
552: *
553: CALL DORMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK,
554: $ WORK( ITAUP ), B, LDB, WORK( NWORK ),
555: $ LWORK-NWORK+1, INFO )
556: *
557: * Zero out below first M rows of B.
558: *
559: CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
560: NWORK = ITAU + M
561: *
562: * Multiply transpose(Q) by B.
563: * (Workspace: need M+NRHS, prefer M+NRHS*NB)
564: *
565: CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B,
566: $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
567: *
568: ELSE
569: *
570: * Path 2 - remaining underdetermined cases.
571: *
572: IE = 1
573: ITAUQ = IE + M
574: ITAUP = ITAUQ + M
575: NWORK = ITAUP + M
576: *
577: * Bidiagonalize A.
578: * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
579: *
580: CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
581: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
582: $ INFO )
583: *
584: * Multiply B by transpose of left bidiagonalizing vectors.
585: * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
586: *
587: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ),
588: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
589: *
590: * Solve the bidiagonal least squares problem.
591: *
592: CALL DLALSD( 'L', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB,
593: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
594: IF( INFO.NE.0 ) THEN
595: GO TO 10
596: END IF
597: *
598: * Multiply B by right bidiagonalizing vectors of A.
599: *
600: CALL DORMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ),
601: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
602: *
603: END IF
604: *
605: * Undo scaling.
606: *
607: IF( IASCL.EQ.1 ) THEN
608: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
609: CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
610: $ INFO )
611: ELSE IF( IASCL.EQ.2 ) THEN
612: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
613: CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
614: $ INFO )
615: END IF
616: IF( IBSCL.EQ.1 ) THEN
617: CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
618: ELSE IF( IBSCL.EQ.2 ) THEN
619: CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
620: END IF
621: *
622: 10 CONTINUE
623: WORK( 1 ) = MAXWRK
1.5 bertrand 624: IWORK( 1 ) = LIWORK
1.1 bertrand 625: RETURN
626: *
627: * End of DGELSD
628: *
629: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgelsd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack