zgelsd.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 591 行 · 第 1/4 页
HTML
591 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>zgelsd.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="ZGELSD.1"></a><a href="zgelsd.f.html#ZGELSD.1">ZGELSD</a>( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
$ WORK, LWORK, RWORK, 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 RWORK( * ), S( * )
COMPLEX*16 A( LDA, * ), B( LDB, * ), 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="ZGELSD.21"></a><a href="zgelsd.f.html#ZGELSD.1">ZGELSD</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 tranformations, 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 the matrix 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 the matrix 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) COMPLEX*16 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) COMPLEX*16 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 matrix X.
</span><span class="comment">*</span><span class="comment"> If m >= n and RANK = n, the residual sum-of-squares for
</span><span class="comment">*</span><span class="comment"> the solution in the i-th column is given by the sum of
</span><span class="comment">*</span><span class="comment"> squares of the modulus 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,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) COMPLEX*16 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"> 2*N + N*NRHS
</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"> 2*M + M*NRHS
</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"> 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 array WORK and the
</span><span class="comment">*</span><span class="comment"> minimum sizes of the arrays RWORK and IWORK, and returns
</span><span class="comment">*</span><span class="comment"> these values as the first entries of the WORK, RWORK and
</span><span class="comment">*</span><span class="comment"> IWORK arrays, and no error message related to LWORK is issued
</span><span class="comment">*</span><span class="comment"> by <a name="XERBLA.112"></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"> RWORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
</span><span class="comment">*</span><span class="comment"> LRWORK >=
</span><span class="comment">*</span><span class="comment"> 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
</span><span class="comment">*</span><span class="comment"> (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"> 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
</span><span class="comment">*</span><span class="comment"> (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.122"></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"> On exit, if INFO = 0, RWORK(1) returns the minimum LRWORK.
</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 >= max(1, 3*MINMN*NLVL + 11*MINMN),
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?