clalsd.f.html

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

HTML
621 行
字号
123
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>clalsd.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="CLALSD.1"></a><a href="clalsd.f.html#CLALSD.1">CLALSD</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
      REAL               RCOND
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IWORK( * )
      REAL               D( * ), E( * ), RWORK( * )
      COMPLEX            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="CLALSD.22"></a><a href="clalsd.f.html#CLALSD.1">CLALSD</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) REAL 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) REAL 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 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) REAL
</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 array, dimension (N * NRHS).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RWORK  (workspace) REAL 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 (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>      REAL               ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
      COMPLEX            CZERO
      PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ) )
<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
      REAL               CS, EPS, ORGNRM, R, RCND, SN, TOL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            ISAMAX
      REAL               <a name="SLAMCH.130"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANST.130"></a><a href="slanst.f.html#SLANST.1">SLANST</a>
      EXTERNAL           ISAMAX, <a name="SLAMCH.131"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANST.131"></a><a href="slanst.f.html#SLANST.1">SLANST</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CCOPY, <a name="CLACPY.134"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>, <a name="CLALSA.134"></a><a href="clalsa.f.html#CLALSA.1">CLALSA</a>, <a name="CLASCL.134"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>, <a name="CLASET.134"></a><a href="claset.f.html#CLASET.1">CLASET</a>, CSROT,
     $                   SGEMM, <a name="SLARTG.135"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>, <a name="SLASCL.135"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>, <a name="SLASDA.135"></a><a href="slasda.f.html#SLASDA.1">SLASDA</a>, <a name="SLASDQ.135"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>, <a name="SLASET.135"></a><a href="slaset.f.html#SLASET.1">SLASET</a>,
     $                   <a name="SLASRT.136"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>, <a name="XERBLA.136"></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">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, AIMAG, CMPLX, INT, LOG, REAL, 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.155"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CLALSD.155"></a><a href="clalsd.f.html#CLALSD.1">CLALSD</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.159"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</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="CLASET.177"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'A'</span>, 1, NRHS, CZERO, CZERO, B, LDB )
         ELSE
            RANK = 1
            CALL <a name="CLASCL.180"></a><a href="clascl.f.html#CLASCL.1">CLASCL</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

clalsd.f.html - 源码说明

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

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

⌨️ 快捷键说明

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