LAPACK 3.3.0

cgebrd.f

Go to the documentation of this file.
00001       SUBROUTINE CGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
00002      $                   INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            INFO, LDA, LWORK, M, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       REAL               D( * ), E( * )
00014       COMPLEX            A( LDA, * ), TAUP( * ), TAUQ( * ),
00015      $                   WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CGEBRD reduces a general complex M-by-N matrix A to upper or lower
00022 *  bidiagonal form B by a unitary transformation: Q**H * A * P = B.
00023 *
00024 *  If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  M       (input) INTEGER
00030 *          The number of rows in the matrix A.  M >= 0.
00031 *
00032 *  N       (input) INTEGER
00033 *          The number of columns in the matrix A.  N >= 0.
00034 *
00035 *  A       (input/output) COMPLEX array, dimension (LDA,N)
00036 *          On entry, the M-by-N general matrix to be reduced.
00037 *          On exit,
00038 *          if m >= n, the diagonal and the first superdiagonal are
00039 *            overwritten with the upper bidiagonal matrix B; the
00040 *            elements below the diagonal, with the array TAUQ, represent
00041 *            the unitary matrix Q as a product of elementary
00042 *            reflectors, and the elements above the first superdiagonal,
00043 *            with the array TAUP, represent the unitary matrix P as
00044 *            a product of elementary reflectors;
00045 *          if m < n, the diagonal and the first subdiagonal are
00046 *            overwritten with the lower bidiagonal matrix B; the
00047 *            elements below the first subdiagonal, with the array TAUQ,
00048 *            represent the unitary matrix Q as a product of
00049 *            elementary reflectors, and the elements above the diagonal,
00050 *            with the array TAUP, represent the unitary matrix P as
00051 *            a product of elementary reflectors.
00052 *          See Further Details.
00053 *
00054 *  LDA     (input) INTEGER
00055 *          The leading dimension of the array A.  LDA >= max(1,M).
00056 *
00057 *  D       (output) REAL array, dimension (min(M,N))
00058 *          The diagonal elements of the bidiagonal matrix B:
00059 *          D(i) = A(i,i).
00060 *
00061 *  E       (output) REAL array, dimension (min(M,N)-1)
00062 *          The off-diagonal elements of the bidiagonal matrix B:
00063 *          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
00064 *          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
00065 *
00066 *  TAUQ    (output) COMPLEX array dimension (min(M,N))
00067 *          The scalar factors of the elementary reflectors which
00068 *          represent the unitary matrix Q. See Further Details.
00069 *
00070 *  TAUP    (output) COMPLEX array, dimension (min(M,N))
00071 *          The scalar factors of the elementary reflectors which
00072 *          represent the unitary matrix P. See Further Details.
00073 *
00074 *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
00075 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00076 *
00077 *  LWORK   (input) INTEGER
00078 *          The length of the array WORK.  LWORK >= max(1,M,N).
00079 *          For optimum performance LWORK >= (M+N)*NB, where NB
00080 *          is the optimal blocksize.
00081 *
00082 *          If LWORK = -1, then a workspace query is assumed; the routine
00083 *          only calculates the optimal size of the WORK array, returns
00084 *          this value as the first entry of the WORK array, and no error
00085 *          message related to LWORK is issued by XERBLA.
00086 *
00087 *  INFO    (output) INTEGER
00088 *          = 0:  successful exit.
00089 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00090 *
00091 *  Further Details
00092 *  ===============
00093 *
00094 *  The matrices Q and P are represented as products of elementary
00095 *  reflectors:
00096 *
00097 *  If m >= n,
00098 *
00099 *     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)
00100 *
00101 *  Each H(i) and G(i) has the form:
00102 *
00103 *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
00104 *
00105 *  where tauq and taup are complex scalars, and v and u are complex
00106 *  vectors; v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in
00107 *  A(i+1:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in
00108 *  A(i,i+2:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
00109 *
00110 *  If m < n,
00111 *
00112 *     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)
00113 *
00114 *  Each H(i) and G(i) has the form:
00115 *
00116 *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
00117 *
00118 *  where tauq and taup are complex scalars, and v and u are complex
00119 *  vectors; v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in
00120 *  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in
00121 *  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
00122 *
00123 *  The contents of A on exit are illustrated by the following examples:
00124 *
00125 *  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
00126 *
00127 *    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
00128 *    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
00129 *    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
00130 *    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
00131 *    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
00132 *    (  v1  v2  v3  v4  v5 )
00133 *
00134 *  where d and e denote diagonal and off-diagonal elements of B, vi
00135 *  denotes an element of the vector defining H(i), and ui an element of
00136 *  the vector defining G(i).
00137 *
00138 *  =====================================================================
00139 *
00140 *     .. Parameters ..
00141       COMPLEX            ONE
00142       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
00143 *     ..
00144 *     .. Local Scalars ..
00145       LOGICAL            LQUERY
00146       INTEGER            I, IINFO, J, LDWRKX, LDWRKY, LWKOPT, MINMN, NB,
00147      $                   NBMIN, NX
00148       REAL               WS
00149 *     ..
00150 *     .. External Subroutines ..
00151       EXTERNAL           CGEBD2, CGEMM, CLABRD, XERBLA
00152 *     ..
00153 *     .. Intrinsic Functions ..
00154       INTRINSIC          MAX, MIN, REAL
00155 *     ..
00156 *     .. External Functions ..
00157       INTEGER            ILAENV
00158       EXTERNAL           ILAENV
00159 *     ..
00160 *     .. Executable Statements ..
00161 *
00162 *     Test the input parameters
00163 *
00164       INFO = 0
00165       NB = MAX( 1, ILAENV( 1, 'CGEBRD', ' ', M, N, -1, -1 ) )
00166       LWKOPT = ( M+N )*NB
00167       WORK( 1 ) = REAL( LWKOPT )
00168       LQUERY = ( LWORK.EQ.-1 )
00169       IF( M.LT.0 ) THEN
00170          INFO = -1
00171       ELSE IF( N.LT.0 ) THEN
00172          INFO = -2
00173       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00174          INFO = -4
00175       ELSE IF( LWORK.LT.MAX( 1, M, N ) .AND. .NOT.LQUERY ) THEN
00176          INFO = -10
00177       END IF
00178       IF( INFO.LT.0 ) THEN
00179          CALL XERBLA( 'CGEBRD', -INFO )
00180          RETURN
00181       ELSE IF( LQUERY ) THEN
00182          RETURN
00183       END IF
00184 *
00185 *     Quick return if possible
00186 *
00187       MINMN = MIN( M, N )
00188       IF( MINMN.EQ.0 ) THEN
00189          WORK( 1 ) = 1
00190          RETURN
00191       END IF
00192 *
00193       WS = MAX( M, N )
00194       LDWRKX = M
00195       LDWRKY = N
00196 *
00197       IF( NB.GT.1 .AND. NB.LT.MINMN ) THEN
00198 *
00199 *        Set the crossover point NX.
00200 *
00201          NX = MAX( NB, ILAENV( 3, 'CGEBRD', ' ', M, N, -1, -1 ) )
00202 *
00203 *        Determine when to switch from blocked to unblocked code.
00204 *
00205          IF( NX.LT.MINMN ) THEN
00206             WS = ( M+N )*NB
00207             IF( LWORK.LT.WS ) THEN
00208 *
00209 *              Not enough work space for the optimal NB, consider using
00210 *              a smaller block size.
00211 *
00212                NBMIN = ILAENV( 2, 'CGEBRD', ' ', M, N, -1, -1 )
00213                IF( LWORK.GE.( M+N )*NBMIN ) THEN
00214                   NB = LWORK / ( M+N )
00215                ELSE
00216                   NB = 1
00217                   NX = MINMN
00218                END IF
00219             END IF
00220          END IF
00221       ELSE
00222          NX = MINMN
00223       END IF
00224 *
00225       DO 30 I = 1, MINMN - NX, NB
00226 *
00227 *        Reduce rows and columns i:i+ib-1 to bidiagonal form and return
00228 *        the matrices X and Y which are needed to update the unreduced
00229 *        part of the matrix
00230 *
00231          CALL CLABRD( M-I+1, N-I+1, NB, A( I, I ), LDA, D( I ), E( I ),
00232      $                TAUQ( I ), TAUP( I ), WORK, LDWRKX,
00233      $                WORK( LDWRKX*NB+1 ), LDWRKY )
00234 *
00235 *        Update the trailing submatrix A(i+ib:m,i+ib:n), using
00236 *        an update of the form  A := A - V*Y' - X*U'
00237 *
00238          CALL CGEMM( 'No transpose', 'Conjugate transpose', M-I-NB+1,
00239      $               N-I-NB+1, NB, -ONE, A( I+NB, I ), LDA,
00240      $               WORK( LDWRKX*NB+NB+1 ), LDWRKY, ONE,
00241      $               A( I+NB, I+NB ), LDA )
00242          CALL CGEMM( 'No transpose', 'No transpose', M-I-NB+1, N-I-NB+1,
00243      $               NB, -ONE, WORK( NB+1 ), LDWRKX, A( I, I+NB ), LDA,
00244      $               ONE, A( I+NB, I+NB ), LDA )
00245 *
00246 *        Copy diagonal and off-diagonal elements of B back into A
00247 *
00248          IF( M.GE.N ) THEN
00249             DO 10 J = I, I + NB - 1
00250                A( J, J ) = D( J )
00251                A( J, J+1 ) = E( J )
00252    10       CONTINUE
00253          ELSE
00254             DO 20 J = I, I + NB - 1
00255                A( J, J ) = D( J )
00256                A( J+1, J ) = E( J )
00257    20       CONTINUE
00258          END IF
00259    30 CONTINUE
00260 *
00261 *     Use unblocked code to reduce the remainder of the matrix
00262 *
00263       CALL CGEBD2( M-I+1, N-I+1, A( I, I ), LDA, D( I ), E( I ),
00264      $             TAUQ( I ), TAUP( I ), WORK, IINFO )
00265       WORK( 1 ) = WS
00266       RETURN
00267 *
00268 *     End of CGEBRD
00269 *
00270       END
 All Files Functions