zlalsd.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 625 行 · 第 1/3 页

HTML
625 行
字号
123
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>zlalsd.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="ZLALSD.1"></a><a href="zlalsd.f.html#ZLALSD.1">ZLALSD</a>( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
     $                   RANK, WORK, RWORK, IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK 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>      CHARACTER          UPLO
      INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
      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   D( * ), E( * ), RWORK( * )
      COMPLEX*16         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="ZLALSD.22"></a><a href="zlalsd.f.html#ZLALSD.1">ZLALSD</a> uses the singular value decomposition of A to solve the least
</span><span class="comment">*</span><span class="comment">  squares problem of finding X to minimize the Euclidean norm of each
</span><span class="comment">*</span><span class="comment">  column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
</span><span class="comment">*</span><span class="comment">  are N-by-NRHS. The solution X overwrites B.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The singular values of A smaller than RCOND times the largest
</span><span class="comment">*</span><span class="comment">  singular value are treated as zero in solving the least squares
</span><span class="comment">*</span><span class="comment">  problem; in this case a minimum norm solution is returned.
</span><span class="comment">*</span><span class="comment">  The actual singular values are returned in D in ascending order.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  This code makes very mild assumptions about floating point
</span><span class="comment">*</span><span class="comment">  arithmetic. It will work on machines with a guard digit in
</span><span class="comment">*</span><span class="comment">  add/subtract, or on those binary machines without guard digits
</span><span class="comment">*</span><span class="comment">  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
</span><span class="comment">*</span><span class="comment">  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">  UPLO   (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">         = 'U': D and E define an upper bidiagonal matrix.
</span><span class="comment">*</span><span class="comment">         = 'L': D and E define a  lower bidiagonal matrix.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  SMLSIZ (input) INTEGER
</span><span class="comment">*</span><span class="comment">         The maximum size of the subproblems at the bottom of the
</span><span class="comment">*</span><span class="comment">         computation tree.
</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 dimension of the  bidiagonal matrix.  N &gt;= 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 columns of B. NRHS must be at least 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D      (input/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">         On entry D contains the main diagonal of the bidiagonal
</span><span class="comment">*</span><span class="comment">         matrix. On exit, if INFO = 0, D contains its singular values.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E      (input/output) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">         Contains the super-diagonal entries of the bidiagonal matrix.
</span><span class="comment">*</span><span class="comment">         On exit, E has been destroyed.
</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 input, B contains the right hand sides of the least
</span><span class="comment">*</span><span class="comment">         squares problem. On output, B contains the solution X.
</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 B in the calling subprogram.
</span><span class="comment">*</span><span class="comment">         LDB must be at least max(1,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">         The singular values of A less than or equal to RCOND times
</span><span class="comment">*</span><span class="comment">         the largest singular value are treated as zero in solving
</span><span class="comment">*</span><span class="comment">         the least squares problem. If RCOND is negative,
</span><span class="comment">*</span><span class="comment">         machine precision is used instead.
</span><span class="comment">*</span><span class="comment">         For example, if diag(S)*X=B were the least squares problem,
</span><span class="comment">*</span><span class="comment">         where diag(S) is a diagonal matrix of singular values, the
</span><span class="comment">*</span><span class="comment">         solution would be X(i) = B(i) / S(i) if S(i) is greater than
</span><span class="comment">*</span><span class="comment">         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
</span><span class="comment">*</span><span class="comment">         RCOND*max(S).
</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 number of singular values of A greater than RCOND times
</span><span class="comment">*</span><span class="comment">         the largest singular value.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK   (workspace) COMPLEX*16 array, dimension at least
</span><span class="comment">*</span><span class="comment">         (N * NRHS).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RWORK  (workspace) DOUBLE PRECISION array, dimension at least
</span><span class="comment">*</span><span class="comment">         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + (SMLSIZ+1)**2),
</span><span class="comment">*</span><span class="comment">         where
</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">
</span><span class="comment">*</span><span class="comment">  IWORK  (workspace) INTEGER array, dimension at least
</span><span class="comment">*</span><span class="comment">         (3*N*NLVL + 11*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">         &lt; 0:  if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">         &gt; 0:  The algorithm failed to compute an singular value while
</span><span class="comment">*</span><span class="comment">               working on the submatrix lying in rows and columns
</span><span class="comment">*</span><span class="comment">               INFO/(N+1) through MOD(INFO,N+1).
</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 )
      COMPLEX*16         CZERO
      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
     $                   GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
     $                   IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
     $                   JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
     $                   PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
     $                   U, VT, Z
      DOUBLE PRECISION   CS, EPS, ORGNRM, RCND, R, SN, TOL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            IDAMAX
      DOUBLE PRECISION   <a name="DLAMCH.132"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.132"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>
      EXTERNAL           IDAMAX, <a name="DLAMCH.133"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.133"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DGEMM, <a name="DLARTG.136"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>, <a name="DLASCL.136"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>, <a name="DLASDA.136"></a><a href="dlasda.f.html#DLASDA.1">DLASDA</a>, <a name="DLASDQ.136"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>, <a name="DLASET.136"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>,
     $                   <a name="DLASRT.137"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</a>, <a name="XERBLA.137"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, ZCOPY, ZDROT, <a name="ZLACPY.137"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZLALSA.137"></a><a href="zlalsa.f.html#ZLALSA.1">ZLALSA</a>,
     $                   <a name="ZLASCL.138"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>, <a name="ZLASET.138"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN
<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 parameters.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span>      IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( NRHS.LT.1 ) THEN
         INFO = -4
      ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
         INFO = -8
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.157"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZLALSD.157"></a><a href="zlalsd.f.html#ZLALSD.1">ZLALSD</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.161"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set up the tolerance.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
         RCND = EPS
      ELSE
         RCND = RCOND
      END IF
<span class="comment">*</span><span class="comment">
</span>      RANK = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible.
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 ) THEN
         RETURN
      ELSE IF( N.EQ.1 ) THEN
         IF( D( 1 ).EQ.ZERO ) THEN
            CALL <a name="ZLASET.179"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'A'</span>, 1, NRHS, CZERO, CZERO, B, LDB )
         ELSE
            RANK = 1
            CALL <a name="ZLASCL.182"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
            D( 1 ) = ABS( D( 1 ) )
         END IF
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Rotate the matrix if it is lower bidiagonal.
</span><span class="comment">*</span><span class="comment">
</span>      IF( UPLO.EQ.<span class="string">'L'</span> ) THEN
123

zlalsd.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 zlalsd.f.html 源码文件，采用 HTML 编程语言编写，共 625 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?