dgelsd.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 553 行 · 第 1/3 页
HTML
553 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>dgelsd.f</title>
<meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!--
body { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default a { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.string { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.string a { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.comment { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.comment a { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
--></style>
</head>
<body>
<pre>
SUBROUTINE <a name="DGELSD.1"></a><a href="dgelsd.f.html#DGELSD.1">DGELSD</a>( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
$ WORK, LWORK, IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK driver routine (version 3.1) --
</span><span class="comment">*</span><span class="comment"> Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment"> November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Scalar Arguments ..
</span> INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
DOUBLE PRECISION RCOND
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> INTEGER IWORK( * )
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Purpose
</span><span class="comment">*</span><span class="comment"> =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> <a name="DGELSD.20"></a><a href="dgelsd.f.html#DGELSD.1">DGELSD</a> computes the minimum-norm solution to a real linear least
</span><span class="comment">*</span><span class="comment"> squares problem:
</span><span class="comment">*</span><span class="comment"> minimize 2-norm(| b - A*x |)
</span><span class="comment">*</span><span class="comment"> using the singular value decomposition (SVD) of A. A is an M-by-N
</span><span class="comment">*</span><span class="comment"> matrix which may be rank-deficient.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Several right hand side vectors b and solution vectors x can be
</span><span class="comment">*</span><span class="comment"> handled in a single call; they are stored as the columns of the
</span><span class="comment">*</span><span class="comment"> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
</span><span class="comment">*</span><span class="comment"> matrix X.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The problem is solved in three steps:
</span><span class="comment">*</span><span class="comment"> (1) Reduce the coefficient matrix A to bidiagonal form with
</span><span class="comment">*</span><span class="comment"> Householder transformations, reducing the original problem
</span><span class="comment">*</span><span class="comment"> into a "bidiagonal least squares problem" (BLS)
</span><span class="comment">*</span><span class="comment"> (2) Solve the BLS using a divide and conquer approach.
</span><span class="comment">*</span><span class="comment"> (3) Apply back all the Householder tranformations to solve
</span><span class="comment">*</span><span class="comment"> the original least squares problem.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The effective rank of A is determined by treating as zero those
</span><span class="comment">*</span><span class="comment"> singular values which are less than RCOND times the largest singular
</span><span class="comment">*</span><span class="comment"> value.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The divide and conquer algorithm makes very mild assumptions about
</span><span class="comment">*</span><span class="comment"> floating point arithmetic. It will work on machines with a guard
</span><span class="comment">*</span><span class="comment"> digit in add/subtract, or on those binary machines without guard
</span><span class="comment">*</span><span class="comment"> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
</span><span class="comment">*</span><span class="comment"> Cray-2. It could conceivably fail on hexadecimal or decimal machines
</span><span class="comment">*</span><span class="comment"> without guard digits, but we know of none.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Arguments
</span><span class="comment">*</span><span class="comment"> =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> M (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of rows of A. M >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> N (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of columns of A. N >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> NRHS (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of right hand sides, i.e., the number of columns
</span><span class="comment">*</span><span class="comment"> of the matrices B and X. NRHS >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> A (input) DOUBLE PRECISION array, dimension (LDA,N)
</span><span class="comment">*</span><span class="comment"> On entry, the M-by-N matrix A.
</span><span class="comment">*</span><span class="comment"> On exit, A has been destroyed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDA (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array A. LDA >= max(1,M).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
</span><span class="comment">*</span><span class="comment"> On entry, the M-by-NRHS right hand side matrix B.
</span><span class="comment">*</span><span class="comment"> On exit, B is overwritten by the N-by-NRHS solution
</span><span class="comment">*</span><span class="comment"> matrix X. If m >= n and RANK = n, the residual
</span><span class="comment">*</span><span class="comment"> sum-of-squares for the solution in the i-th column is given
</span><span class="comment">*</span><span class="comment"> by the sum of squares of elements n+1:m in that column.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDB (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array B. LDB >= max(1,max(M,N)).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> S (output) DOUBLE PRECISION array, dimension (min(M,N))
</span><span class="comment">*</span><span class="comment"> The singular values of A in decreasing order.
</span><span class="comment">*</span><span class="comment"> The condition number of A in the 2-norm = S(1)/S(min(m,n)).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> RCOND (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> RCOND is used to determine the effective rank of A.
</span><span class="comment">*</span><span class="comment"> Singular values S(i) <= RCOND*S(1) are treated as zero.
</span><span class="comment">*</span><span class="comment"> If RCOND < 0, machine precision is used instead.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> RANK (output) INTEGER
</span><span class="comment">*</span><span class="comment"> The effective rank of A, i.e., the number of singular values
</span><span class="comment">*</span><span class="comment"> which are greater than RCOND*S(1).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
</span><span class="comment">*</span><span class="comment"> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LWORK (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The dimension of the array WORK. LWORK must be at least 1.
</span><span class="comment">*</span><span class="comment"> The exact minimum amount of workspace needed depends on M,
</span><span class="comment">*</span><span class="comment"> N and NRHS. As long as LWORK is at least
</span><span class="comment">*</span><span class="comment"> 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,
</span><span class="comment">*</span><span class="comment"> if M is greater than or equal to N or
</span><span class="comment">*</span><span class="comment"> 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,
</span><span class="comment">*</span><span class="comment"> if M is less than N, the code will execute correctly.
</span><span class="comment">*</span><span class="comment"> SMLSIZ is returned by <a name="ILAENV.104"></a><a href="hfy-index.html#ILAENV">ILAENV</a> and is equal to the maximum
</span><span class="comment">*</span><span class="comment"> size of the subproblems at the bottom of the computation
</span><span class="comment">*</span><span class="comment"> tree (usually about 25), and
</span><span class="comment">*</span><span class="comment"> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
</span><span class="comment">*</span><span class="comment"> For good performance, LWORK should generally be larger.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> If LWORK = -1, then a workspace query is assumed; the routine
</span><span class="comment">*</span><span class="comment"> only calculates the optimal size of the WORK array, returns
</span><span class="comment">*</span><span class="comment"> this value as the first entry of the WORK array, and no error
</span><span class="comment">*</span><span class="comment"> message related to LWORK is issued by <a name="XERBLA.113"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
</span><span class="comment">*</span><span class="comment"> LIWORK >= 3 * MINMN * NLVL + 11 * MINMN,
</span><span class="comment">*</span><span class="comment"> where MINMN = MIN( M,N ).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> INFO (output) INTEGER
</span><span class="comment">*</span><span class="comment"> = 0: successful exit
</span><span class="comment">*</span><span class="comment"> < 0: if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment"> > 0: the algorithm for computing the SVD failed to converge;
</span><span class="comment">*</span><span class="comment"> if INFO = i, i off-diagonal elements of an intermediate
</span><span class="comment">*</span><span class="comment"> bidiagonal form did not converge to zero.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Further Details
</span><span class="comment">*</span><span class="comment"> ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Based on contributions by
</span><span class="comment">*</span><span class="comment"> Ming Gu and Ren-Cang Li, Computer Science Division, University of
</span><span class="comment">*</span><span class="comment"> California at Berkeley, USA
</span><span class="comment">*</span><span class="comment"> Osni Marques, LBNL/NERSC, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Parameters ..
</span> DOUBLE PRECISION ZERO, ONE, TWO
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL LQUERY
INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
$ LDWORK, MAXMN, MAXWRK, MINMN, MINWRK, MM,
$ MNTHR, NLVL, NWORK, SMLSIZ, WLALSD
DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL <a name="DGEBRD.148"></a><a href="dgebrd.f.html#DGEBRD.1">DGEBRD</a>, <a name="DGELQF.148"></a><a href="dgelqf.f.html#DGELQF.1">DGELQF</a>, <a name="DGEQRF.148"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>, <a name="DLABAD.148"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>, <a name="DLACPY.148"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>, <a name="DLALSD.148"></a><a href="dlalsd.f.html#DLALSD.1">DLALSD</a>,
$ <a name="DLASCL.149"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>, <a name="DLASET.149"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>, <a name="DORMBR.149"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>, <a name="DORMLQ.149"></a><a href="dormlq.f.html#DORMLQ.1">DORMLQ</a>, <a name="DORMQR.149"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>, <a name="XERBLA.149"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> INTEGER <a name="ILAENV.152"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
DOUBLE PRECISION <a name="DLAMCH.153"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANGE.153"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
EXTERNAL <a name="ILAENV.154"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, <a name="DLAMCH.154"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANGE.154"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC DBLE, INT, LOG, MAX, MIN
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Test the input arguments.
</span><span class="comment">*</span><span class="comment">
</span> INFO = 0
MINMN = MIN( M, N )
MAXMN = MAX( M, N )
MNTHR = <a name="ILAENV.166"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 6, <span class="string">'<a name="DGELSD.166"></a><a href="dgelsd.f.html#DGELSD.1">DGELSD</a>'</span>, <span class="string">' '</span>, M, N, NRHS, -1 )
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?