sbdsqr.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 767 行 · 第 1/4 页
HTML
767 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>sbdsqr.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="SBDSQR.1"></a><a href="sbdsqr.f.html#SBDSQR.1">SBDSQR</a>( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
$ LDU, C, LDC, WORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK routine (version 3.1.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"> January 2007
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Scalar Arguments ..
</span> CHARACTER UPLO
INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> REAL C( LDC, * ), D( * ), E( * ), U( LDU, * ),
$ VT( LDVT, * ), 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="SBDSQR.20"></a><a href="sbdsqr.f.html#SBDSQR.1">SBDSQR</a> computes the singular values and, optionally, the right and/or
</span><span class="comment">*</span><span class="comment"> left singular vectors from the singular value decomposition (SVD) of
</span><span class="comment">*</span><span class="comment"> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
</span><span class="comment">*</span><span class="comment"> zero-shift QR algorithm. The SVD of B has the form
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> B = Q * S * P**T
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> where S is the diagonal matrix of singular values, Q is an orthogonal
</span><span class="comment">*</span><span class="comment"> matrix of left singular vectors, and P is an orthogonal matrix of
</span><span class="comment">*</span><span class="comment"> right singular vectors. If left singular vectors are requested, this
</span><span class="comment">*</span><span class="comment"> subroutine actually returns U*Q instead of Q, and, if right singular
</span><span class="comment">*</span><span class="comment"> vectors are requested, this subroutine returns P**T*VT instead of
</span><span class="comment">*</span><span class="comment"> P**T, for given real input matrices U and VT. When U and VT are the
</span><span class="comment">*</span><span class="comment"> orthogonal matrices that reduce a general matrix A to bidiagonal
</span><span class="comment">*</span><span class="comment"> form: A = U*B*VT, as computed by <a name="SGEBRD.34"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>, then
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> A = (U*Q) * S * (P**T*VT)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> is the SVD of A. Optionally, the subroutine may also compute Q**T*C
</span><span class="comment">*</span><span class="comment"> for a given real input matrix C.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> See "Computing Small Singular Values of Bidiagonal Matrices With
</span><span class="comment">*</span><span class="comment"> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
</span><span class="comment">*</span><span class="comment"> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
</span><span class="comment">*</span><span class="comment"> no. 5, pp. 873-912, Sept 1990) and
</span><span class="comment">*</span><span class="comment"> "Accurate singular values and differential qd algorithms," by
</span><span class="comment">*</span><span class="comment"> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
</span><span class="comment">*</span><span class="comment"> Department, University of California at Berkeley, July 1992
</span><span class="comment">*</span><span class="comment"> for a detailed description of the algorithm.
</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"> UPLO (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment"> = 'U': B is upper bidiagonal;
</span><span class="comment">*</span><span class="comment"> = 'L': B is lower bidiagonal.
</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 order of the matrix B. N >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> NCVT (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of columns of the matrix VT. NCVT >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> NRU (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of rows of the matrix U. NRU >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> NCC (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of columns of the matrix C. NCC >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> D (input/output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment"> On entry, the n diagonal elements of the bidiagonal matrix B.
</span><span class="comment">*</span><span class="comment"> On exit, if INFO=0, the singular values of B in decreasing
</span><span class="comment">*</span><span class="comment"> order.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> E (input/output) REAL array, dimension (N-1)
</span><span class="comment">*</span><span class="comment"> On entry, the N-1 offdiagonal elements of the bidiagonal
</span><span class="comment">*</span><span class="comment"> matrix B.
</span><span class="comment">*</span><span class="comment"> On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
</span><span class="comment">*</span><span class="comment"> will contain the diagonal and superdiagonal elements of a
</span><span class="comment">*</span><span class="comment"> bidiagonal matrix orthogonally equivalent to the one given
</span><span class="comment">*</span><span class="comment"> as input.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> VT (input/output) REAL array, dimension (LDVT, NCVT)
</span><span class="comment">*</span><span class="comment"> On entry, an N-by-NCVT matrix VT.
</span><span class="comment">*</span><span class="comment"> On exit, VT is overwritten by P**T * VT.
</span><span class="comment">*</span><span class="comment"> Not referenced if NCVT = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDVT (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array VT.
</span><span class="comment">*</span><span class="comment"> LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> U (input/output) REAL array, dimension (LDU, N)
</span><span class="comment">*</span><span class="comment"> On entry, an NRU-by-N matrix U.
</span><span class="comment">*</span><span class="comment"> On exit, U is overwritten by U * Q.
</span><span class="comment">*</span><span class="comment"> Not referenced if NRU = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDU (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array U. LDU >= max(1,NRU).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> C (input/output) REAL array, dimension (LDC, NCC)
</span><span class="comment">*</span><span class="comment"> On entry, an N-by-NCC matrix C.
</span><span class="comment">*</span><span class="comment"> On exit, C is overwritten by Q**T * C.
</span><span class="comment">*</span><span class="comment"> Not referenced if NCC = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDC (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array C.
</span><span class="comment">*</span><span class="comment"> LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> WORK (workspace) REAL array, dimension (2*N)
</span><span class="comment">*</span><span class="comment"> if NCVT = NRU = NCC = 0, (max(1, 4*N)) otherwise
</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 did not converge; D and E contain the
</span><span class="comment">*</span><span class="comment"> elements of a bidiagonal matrix which is orthogonally
</span><span class="comment">*</span><span class="comment"> similar to the input matrix B; if INFO = i, i
</span><span class="comment">*</span><span class="comment"> elements of E have not converged to zero.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Internal Parameters
</span><span class="comment">*</span><span class="comment"> ===================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> TOLMUL REAL, default = max(10,min(100,EPS**(-1/8)))
</span><span class="comment">*</span><span class="comment"> TOLMUL controls the convergence criterion of the QR loop.
</span><span class="comment">*</span><span class="comment"> If it is positive, TOLMUL*EPS is the desired relative
</span><span class="comment">*</span><span class="comment"> precision in the computed singular values.
</span><span class="comment">*</span><span class="comment"> If it is negative, abs(TOLMUL*EPS*sigma_max) is the
</span><span class="comment">*</span><span class="comment"> desired absolute accuracy in the computed singular
</span><span class="comment">*</span><span class="comment"> values (corresponds to relative accuracy
</span><span class="comment">*</span><span class="comment"> abs(TOLMUL*EPS) in the largest singular value.
</span><span class="comment">*</span><span class="comment"> abs(TOLMUL) should be between 1 and 1/EPS, and preferably
</span><span class="comment">*</span><span class="comment"> between 10 (for fast convergence) and .1/EPS
</span><span class="comment">*</span><span class="comment"> (for there to be some accuracy in the results).
</span><span class="comment">*</span><span class="comment"> Default is to lose at either one eighth or 2 of the
</span><span class="comment">*</span><span class="comment"> available decimal digits in each computed singular value
</span><span class="comment">*</span><span class="comment"> (whichever is smaller).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> MAXITR INTEGER, default = 6
</span><span class="comment">*</span><span class="comment"> MAXITR controls the maximum number of passes of the
</span><span class="comment">*</span><span class="comment"> algorithm through its inner loop. The algorithms stops
</span><span class="comment">*</span><span class="comment"> (and so fails to converge) if the number of passes
</span><span class="comment">*</span><span class="comment"> through the inner loop exceeds MAXITR*N**2.
</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> REAL ZERO
PARAMETER ( ZERO = 0.0E0 )
REAL ONE
PARAMETER ( ONE = 1.0E0 )
REAL NEGONE
PARAMETER ( NEGONE = -1.0E0 )
REAL HNDRTH
PARAMETER ( HNDRTH = 0.01E0 )
REAL TEN
PARAMETER ( TEN = 10.0E0 )
REAL HNDRD
PARAMETER ( HNDRD = 100.0E0 )
REAL MEIGTH
PARAMETER ( MEIGTH = -0.125E0 )
INTEGER MAXITR
PARAMETER ( MAXITR = 6 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL LOWER, ROTATE
INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
$ NM12, NM13, OLDLL, OLDM
REAL ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
$ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
$ SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
$ SN, THRESH, TOL, TOLMUL, UNFL
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> LOGICAL <a name="LSAME.173"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?