zlalsd.f

计算矩阵的经典开源库．全世界都在用它．相信你也不能例外．

      SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
     $                   RANK, WORK, RWORK, IWORK, INFO )
*
*  -- LAPACK routine (instrumented to count ops, version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1999
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
      DOUBLE PRECISION   RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
      COMPLEX*16         B( LDB, * ), WORK( * )
*     ..
*     .. Common block to return operation count ..
      COMMON             / LATIME / OPS, ITCNT
*     ..
*     .. Scalars in Common ..
      DOUBLE PRECISION   ITCNT, OPS
*     ..
*
*  Purpose
*  =======
*
*  ZLALSD uses the singular value decomposition of A to solve the least
*  squares problem of finding X to minimize the Euclidean norm of each
*  column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
*  are N-by-NRHS. The solution X overwrites B.
*
*  The singular values of A smaller than RCOND times the largest
*  singular value are treated as zero in solving the least squares
*  problem; in this case a minimum norm solution is returned.
*  The actual singular values are returned in D in ascending order.
*
*  This code makes very mild assumptions about floating point
*  arithmetic. It will work on machines with a guard digit in
*  add/subtract, or on those binary machines without guard digits
*  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
*  It could conceivably fail on hexadecimal or decimal machines
*  without guard digits, but we know of none.
*
*  Arguments
*  =========
*
*  UPLO   (input) CHARACTER*1
*         = 'U': D and E define an upper bidiagonal matrix.
*         = 'L': D and E define a  lower bidiagonal matrix.
*
*  SMLSIZ (input) INTEGER
*         The maximum size of the subproblems at the bottom of the
*         computation tree.
*
*  N      (input) INTEGER
*         The dimension of the  bidiagonal matrix.  N >= 0.
*
*  NRHS   (input) INTEGER
*         The number of columns of B. NRHS must be at least 1.
*
*  D      (input/output) DOUBLE PRECISION array, dimension (N)
*         On entry D contains the main diagonal of the bidiagonal
*         matrix. On exit, if INFO = 0, D contains its singular values.
*
*  E      (input) DOUBLE PRECISION array, dimension (N-1)
*         Contains the super-diagonal entries of the bidiagonal matrix.
*         On exit, E has been destroyed.
*
*  B      (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
*         On input, B contains the right hand sides of the least
*         squares problem. On output, B contains the solution X.
*
*  LDB    (input) INTEGER
*         The leading dimension of B in the calling subprogram.
*         LDB must be at least max(1,N).
*
*  RCOND  (input) DOUBLE PRECISION
*         The singular values of A less than or equal to RCOND times
*         the largest singular value are treated as zero in solving
*         the least squares problem. If RCOND is negative,
*         machine precision is used instead.
*         For example, if diag(S)*X=B were the least squares problem,
*         where diag(S) is a diagonal matrix of singular values, the
*         solution would be X(i) = B(i) / S(i) if S(i) is greater than
*         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
*         RCOND*max(S).
*
*  RANK   (output) INTEGER
*         The number of singular values of A greater than RCOND times
*         the largest singular value.
*
*  WORK   (workspace) COMPLEX*16 array, dimension at least
*         (N * NRHS).
*
*  RWORK  (workspace) DOUBLE PRECISION array, dimension at least
*         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + (SMLSIZ+1)**2),
*         where
*         NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
*
*  IWORK  (workspace) INTEGER array, dimension at least
*         (3*N*NLVL + 11*N).
*
*  INFO   (output) INTEGER
*         = 0:  successful exit.
*         < 0:  if INFO = -i, the i-th argument had an illegal value.
*         > 0:  The algorithm failed to compute an singular value while
*               working on the submatrix lying in rows and columns
*               INFO/(N+1) through MOD(INFO,N+1).
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
      COMPLEX*16         CZERO
      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
*     ..
*     .. Local Scalars ..
      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, R, SN, TOL
*     ..
*     .. External Subroutines ..
      EXTERNAL           ZDROT, ZCOPY, ZLACPY,
     $                   ZLALSA, ZLASCL, ZLASET, DGEMM,
     $                   DLARTG, DLASCL, DLASDA, DLASDQ,
     $                   DLASET, DLASRT, XERBLA
*     ..
*     .. External Functions ..
      INTEGER            IDAMAX
      DOUBLE PRECISION   DLAMCH, DLANST, DOPBL3
      EXTERNAL           IDAMAX, DLAMCH, DLANST, DOPBL3
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          DCMPLX, DBLE, DIMAG, ABS, INT, LOG, SIGN
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      INFO = 0
*
      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 XERBLA( 'ZLALSD', -INFO )
         RETURN
      END IF
*
      EPS = DLAMCH( 'Epsilon' )
*
*     Set up the tolerance.
*
      IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
         RCOND = EPS
      END IF
*
      RANK = 0
*
*     Quick return if possible.
*
      IF( N.EQ.0 ) THEN
         RETURN
      ELSE IF( N.EQ.1 ) THEN
         IF( D( 1 ).EQ.ZERO ) THEN
            CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
         ELSE
            RANK = 1
            OPS = OPS + DBLE( 2*NRHS )
            CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
            D( 1 ) = ABS( D( 1 ) )
         END IF
         RETURN
      END IF
*
*     Rotate the matrix if it is lower bidiagonal.
*
      IF( UPLO.EQ.'L' ) THEN
         OPS = OPS + DBLE( 6*( N-1 ) )
         DO 10 I = 1, N - 1
            CALL DLARTG( D( I ), E( I ), CS, SN, R )
            D( I ) = R
            E( I ) = SN*D( I+1 )
            D( I+1 ) = CS*D( I+1 )
            IF( NRHS.EQ.1 ) THEN
               OPS = OPS + DBLE( 12 )
               CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
            ELSE
               RWORK( I*2-1 ) = CS
               RWORK( I*2 ) = SN
            END IF
   10    CONTINUE
         IF( NRHS.GT.1 ) THEN
            OPS = OPS + DBLE( 12*( N-1 )*NRHS )
            DO 30 I = 1, NRHS
               DO 20 J = 1, N - 1
                  CS = RWORK( J*2-1 )
                  SN = RWORK( J*2 )
                  CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
   20          CONTINUE
   30       CONTINUE
         END IF
      END IF
*
*     Scale.
*
      NM1 = N - 1
      ORGNRM = DLANST( 'M', N, D, E )
      IF( ORGNRM.EQ.ZERO ) THEN
         CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
         RETURN
      END IF
*
      OPS = OPS + DBLE( N + NM1 )
      CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
      CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
*
*     If N is smaller than the minimum divide size SMLSIZ, then solve
*     the problem with another solver.
*
      IF( N.LE.SMLSIZ ) THEN
         IRWU = 1
         IRWVT = IRWU + N*N
         IRWWRK = IRWVT + N*N
         IRWRB = IRWWRK
         IRWIB = IRWRB + N*NRHS
         IRWB = IRWIB + N*NRHS
         CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
         CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
         CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
     $                RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
     $                RWORK( IRWWRK ), INFO )
         IF( INFO.NE.0 ) THEN
            RETURN
         END IF
*
*        In the real version, B is passed to DLASDQ and multiplied
*        internally by Q'. Here B is complex and that product is
*        computed below in two steps (real and imaginary parts).
*
         J = IRWB - 1
         DO 50 JCOL = 1, NRHS
            DO 40 JROW = 1, N
               J = J + 1
               RWORK( J ) = DBLE( B( JROW, JCOL ) )
   40       CONTINUE
   50    CONTINUE
         OPS = OPS + DOPBL3( 'DGEMM ', N, NRHS, N ) 
         CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
     $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
         J = IRWB - 1
         DO 70 JCOL = 1, NRHS
            DO 60 JROW = 1, N
               J = J + 1
               RWORK( J ) = DIMAG( B( JROW, JCOL ) )
   60       CONTINUE
   70    CONTINUE
         OPS = OPS + DOPBL3( 'DGEMM ', N, NRHS, N ) 
         CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
     $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
         JREAL = IRWRB - 1
         JIMAG = IRWIB - 1
         DO 90 JCOL = 1, NRHS
            DO 80 JROW = 1, N
               JREAL = JREAL + 1
               JIMAG = JIMAG + 1
               B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
     $                           RWORK( JIMAG ) )
   80       CONTINUE
   90    CONTINUE
*
         OPS = OPS + DBLE( 1 )
         TOL = RCOND*ABS( D( IDAMAX( N, D, 1 ) ) )
         DO 100 I = 1, N
            IF( D( I ).LE.TOL ) THEN
               CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
            ELSE
               OPS = OPS + DBLE( 6*NRHS )
               CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
     $                      LDB, INFO )
               RANK = RANK + 1
            END IF
  100    CONTINUE
*
*        Since B is complex, the following call to DGEMM is performed
*        in two steps (real and imaginary parts). That is for V * B
*        (in the real version of the code V' is stored in WORK).
*
*        CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
*    $               WORK( NWORK ), N )
*
         J = IRWB - 1
         DO 120 JCOL = 1, NRHS
            DO 110 JROW = 1, N
               J = J + 1
               RWORK( J ) = DBLE( B( JROW, JCOL ) )
  110       CONTINUE
  120    CONTINUE