LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
clabrd.f
Go to the documentation of this file.
00001 *> \brief \b CLABRD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CLABRD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clabrd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clabrd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clabrd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y,
00022 *                          LDY )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            LDA, LDX, LDY, M, N, NB
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       REAL               D( * ), E( * )
00029 *       COMPLEX            A( LDA, * ), TAUP( * ), TAUQ( * ), X( LDX, * ),
00030 *      $                   Y( LDY, * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> CLABRD reduces the first NB rows and columns of a complex general
00040 *> m by n matrix A to upper or lower real bidiagonal form by a unitary
00041 *> transformation Q**H * A * P, and returns the matrices X and Y which
00042 *> are needed to apply the transformation to the unreduced part of A.
00043 *>
00044 *> If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower
00045 *> bidiagonal form.
00046 *>
00047 *> This is an auxiliary routine called by CGEBRD
00048 *> \endverbatim
00049 *
00050 *  Arguments:
00051 *  ==========
00052 *
00053 *> \param[in] M
00054 *> \verbatim
00055 *>          M is INTEGER
00056 *>          The number of rows in the matrix A.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] N
00060 *> \verbatim
00061 *>          N is INTEGER
00062 *>          The number of columns in the matrix A.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] NB
00066 *> \verbatim
00067 *>          NB is INTEGER
00068 *>          The number of leading rows and columns of A to be reduced.
00069 *> \endverbatim
00070 *>
00071 *> \param[in,out] A
00072 *> \verbatim
00073 *>          A is COMPLEX array, dimension (LDA,N)
00074 *>          On entry, the m by n general matrix to be reduced.
00075 *>          On exit, the first NB rows and columns of the matrix are
00076 *>          overwritten; the rest of the array is unchanged.
00077 *>          If m >= n, elements on and below the diagonal in the first NB
00078 *>            columns, with the array TAUQ, represent the unitary
00079 *>            matrix Q as a product of elementary reflectors; and
00080 *>            elements above the diagonal in the first NB rows, with the
00081 *>            array TAUP, represent the unitary matrix P as a product
00082 *>            of elementary reflectors.
00083 *>          If m < n, elements below the diagonal in the first NB
00084 *>            columns, with the array TAUQ, represent the unitary
00085 *>            matrix Q as a product of elementary reflectors, and
00086 *>            elements on and above the diagonal in the first NB rows,
00087 *>            with the array TAUP, represent the unitary matrix P as
00088 *>            a product of elementary reflectors.
00089 *>          See Further Details.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] LDA
00093 *> \verbatim
00094 *>          LDA is INTEGER
00095 *>          The leading dimension of the array A.  LDA >= max(1,M).
00096 *> \endverbatim
00097 *>
00098 *> \param[out] D
00099 *> \verbatim
00100 *>          D is REAL array, dimension (NB)
00101 *>          The diagonal elements of the first NB rows and columns of
00102 *>          the reduced matrix.  D(i) = A(i,i).
00103 *> \endverbatim
00104 *>
00105 *> \param[out] E
00106 *> \verbatim
00107 *>          E is REAL array, dimension (NB)
00108 *>          The off-diagonal elements of the first NB rows and columns of
00109 *>          the reduced matrix.
00110 *> \endverbatim
00111 *>
00112 *> \param[out] TAUQ
00113 *> \verbatim
00114 *>          TAUQ is COMPLEX array dimension (NB)
00115 *>          The scalar factors of the elementary reflectors which
00116 *>          represent the unitary matrix Q. See Further Details.
00117 *> \endverbatim
00118 *>
00119 *> \param[out] TAUP
00120 *> \verbatim
00121 *>          TAUP is COMPLEX array, dimension (NB)
00122 *>          The scalar factors of the elementary reflectors which
00123 *>          represent the unitary matrix P. See Further Details.
00124 *> \endverbatim
00125 *>
00126 *> \param[out] X
00127 *> \verbatim
00128 *>          X is COMPLEX array, dimension (LDX,NB)
00129 *>          The m-by-nb matrix X required to update the unreduced part
00130 *>          of A.
00131 *> \endverbatim
00132 *>
00133 *> \param[in] LDX
00134 *> \verbatim
00135 *>          LDX is INTEGER
00136 *>          The leading dimension of the array X. LDX >= max(1,M).
00137 *> \endverbatim
00138 *>
00139 *> \param[out] Y
00140 *> \verbatim
00141 *>          Y is COMPLEX array, dimension (LDY,NB)
00142 *>          The n-by-nb matrix Y required to update the unreduced part
00143 *>          of A.
00144 *> \endverbatim
00145 *>
00146 *> \param[in] LDY
00147 *> \verbatim
00148 *>          LDY is INTEGER
00149 *>          The leading dimension of the array Y. LDY >= max(1,N).
00150 *> \endverbatim
00151 *
00152 *  Authors:
00153 *  ========
00154 *
00155 *> \author Univ. of Tennessee 
00156 *> \author Univ. of California Berkeley 
00157 *> \author Univ. of Colorado Denver 
00158 *> \author NAG Ltd. 
00159 *
00160 *> \date November 2011
00161 *
00162 *> \ingroup complexOTHERauxiliary
00163 *
00164 *> \par Further Details:
00165 *  =====================
00166 *>
00167 *> \verbatim
00168 *>
00169 *>  The matrices Q and P are represented as products of elementary
00170 *>  reflectors:
00171 *>
00172 *>     Q = H(1) H(2) . . . H(nb)  and  P = G(1) G(2) . . . G(nb)
00173 *>
00174 *>  Each H(i) and G(i) has the form:
00175 *>
00176 *>     H(i) = I - tauq * v * v**H  and G(i) = I - taup * u * u**H
00177 *>
00178 *>  where tauq and taup are complex scalars, and v and u are complex
00179 *>  vectors.
00180 *>
00181 *>  If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in
00182 *>  A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in
00183 *>  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
00184 *>
00185 *>  If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in
00186 *>  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in
00187 *>  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
00188 *>
00189 *>  The elements of the vectors v and u together form the m-by-nb matrix
00190 *>  V and the nb-by-n matrix U**H which are needed, with X and Y, to apply
00191 *>  the transformation to the unreduced part of the matrix, using a block
00192 *>  update of the form:  A := A - V*Y**H - X*U**H.
00193 *>
00194 *>  The contents of A on exit are illustrated by the following examples
00195 *>  with nb = 2:
00196 *>
00197 *>  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
00198 *>
00199 *>    (  1   1   u1  u1  u1 )           (  1   u1  u1  u1  u1  u1 )
00200 *>    (  v1  1   1   u2  u2 )           (  1   1   u2  u2  u2  u2 )
00201 *>    (  v1  v2  a   a   a  )           (  v1  1   a   a   a   a  )
00202 *>    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
00203 *>    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
00204 *>    (  v1  v2  a   a   a  )
00205 *>
00206 *>  where a denotes an element of the original matrix which is unchanged,
00207 *>  vi denotes an element of the vector defining H(i), and ui an element
00208 *>  of the vector defining G(i).
00209 *> \endverbatim
00210 *>
00211 *  =====================================================================
00212       SUBROUTINE CLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y,
00213      $                   LDY )
00214 *
00215 *  -- LAPACK auxiliary routine (version 3.4.0) --
00216 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00217 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00218 *     November 2011
00219 *
00220 *     .. Scalar Arguments ..
00221       INTEGER            LDA, LDX, LDY, M, N, NB
00222 *     ..
00223 *     .. Array Arguments ..
00224       REAL               D( * ), E( * )
00225       COMPLEX            A( LDA, * ), TAUP( * ), TAUQ( * ), X( LDX, * ),
00226      $                   Y( LDY, * )
00227 *     ..
00228 *
00229 *  =====================================================================
00230 *
00231 *     .. Parameters ..
00232       COMPLEX            ZERO, ONE
00233       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ),
00234      $                   ONE = ( 1.0E+0, 0.0E+0 ) )
00235 *     ..
00236 *     .. Local Scalars ..
00237       INTEGER            I
00238       COMPLEX            ALPHA
00239 *     ..
00240 *     .. External Subroutines ..
00241       EXTERNAL           CGEMV, CLACGV, CLARFG, CSCAL
00242 *     ..
00243 *     .. Intrinsic Functions ..
00244       INTRINSIC          MIN
00245 *     ..
00246 *     .. Executable Statements ..
00247 *
00248 *     Quick return if possible
00249 *
00250       IF( M.LE.0 .OR. N.LE.0 )
00251      $   RETURN
00252 *
00253       IF( M.GE.N ) THEN
00254 *
00255 *        Reduce to upper bidiagonal form
00256 *
00257          DO 10 I = 1, NB
00258 *
00259 *           Update A(i:m,i)
00260 *
00261             CALL CLACGV( I-1, Y( I, 1 ), LDY )
00262             CALL CGEMV( 'No transpose', M-I+1, I-1, -ONE, A( I, 1 ),
00263      $                  LDA, Y( I, 1 ), LDY, ONE, A( I, I ), 1 )
00264             CALL CLACGV( I-1, Y( I, 1 ), LDY )
00265             CALL CGEMV( 'No transpose', M-I+1, I-1, -ONE, X( I, 1 ),
00266      $                  LDX, A( 1, I ), 1, ONE, A( I, I ), 1 )
00267 *
00268 *           Generate reflection Q(i) to annihilate A(i+1:m,i)
00269 *
00270             ALPHA = A( I, I )
00271             CALL CLARFG( M-I+1, ALPHA, A( MIN( I+1, M ), I ), 1,
00272      $                   TAUQ( I ) )
00273             D( I ) = ALPHA
00274             IF( I.LT.N ) THEN
00275                A( I, I ) = ONE
00276 *
00277 *              Compute Y(i+1:n,i)
00278 *
00279                CALL CGEMV( 'Conjugate transpose', M-I+1, N-I, ONE,
00280      $                     A( I, I+1 ), LDA, A( I, I ), 1, ZERO,
00281      $                     Y( I+1, I ), 1 )
00282                CALL CGEMV( 'Conjugate transpose', M-I+1, I-1, ONE,
00283      $                     A( I, 1 ), LDA, A( I, I ), 1, ZERO,
00284      $                     Y( 1, I ), 1 )
00285                CALL CGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
00286      $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
00287                CALL CGEMV( 'Conjugate transpose', M-I+1, I-1, ONE,
00288      $                     X( I, 1 ), LDX, A( I, I ), 1, ZERO,
00289      $                     Y( 1, I ), 1 )
00290                CALL CGEMV( 'Conjugate transpose', I-1, N-I, -ONE,
00291      $                     A( 1, I+1 ), LDA, Y( 1, I ), 1, ONE,
00292      $                     Y( I+1, I ), 1 )
00293                CALL CSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
00294 *
00295 *              Update A(i,i+1:n)
00296 *
00297                CALL CLACGV( N-I, A( I, I+1 ), LDA )
00298                CALL CLACGV( I, A( I, 1 ), LDA )
00299                CALL CGEMV( 'No transpose', N-I, I, -ONE, Y( I+1, 1 ),
00300      $                     LDY, A( I, 1 ), LDA, ONE, A( I, I+1 ), LDA )
00301                CALL CLACGV( I, A( I, 1 ), LDA )
00302                CALL CLACGV( I-1, X( I, 1 ), LDX )
00303                CALL CGEMV( 'Conjugate transpose', I-1, N-I, -ONE,
00304      $                     A( 1, I+1 ), LDA, X( I, 1 ), LDX, ONE,
00305      $                     A( I, I+1 ), LDA )
00306                CALL CLACGV( I-1, X( I, 1 ), LDX )
00307 *
00308 *              Generate reflection P(i) to annihilate A(i,i+2:n)
00309 *
00310                ALPHA = A( I, I+1 )
00311                CALL CLARFG( N-I, ALPHA, A( I, MIN( I+2, N ) ),
00312      $                      LDA, TAUP( I ) )
00313                E( I ) = ALPHA
00314                A( I, I+1 ) = ONE
00315 *
00316 *              Compute X(i+1:m,i)
00317 *
00318                CALL CGEMV( 'No transpose', M-I, N-I, ONE, A( I+1, I+1 ),
00319      $                     LDA, A( I, I+1 ), LDA, ZERO, X( I+1, I ), 1 )
00320                CALL CGEMV( 'Conjugate transpose', N-I, I, ONE,
00321      $                     Y( I+1, 1 ), LDY, A( I, I+1 ), LDA, ZERO,
00322      $                     X( 1, I ), 1 )
00323                CALL CGEMV( 'No transpose', M-I, I, -ONE, A( I+1, 1 ),
00324      $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
00325                CALL CGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ),
00326      $                     LDA, A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
00327                CALL CGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
00328      $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
00329                CALL CSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
00330                CALL CLACGV( N-I, A( I, I+1 ), LDA )
00331             END IF
00332    10    CONTINUE
00333       ELSE
00334 *
00335 *        Reduce to lower bidiagonal form
00336 *
00337          DO 20 I = 1, NB
00338 *
00339 *           Update A(i,i:n)
00340 *
00341             CALL CLACGV( N-I+1, A( I, I ), LDA )
00342             CALL CLACGV( I-1, A( I, 1 ), LDA )
00343             CALL CGEMV( 'No transpose', N-I+1, I-1, -ONE, Y( I, 1 ),
00344      $                  LDY, A( I, 1 ), LDA, ONE, A( I, I ), LDA )
00345             CALL CLACGV( I-1, A( I, 1 ), LDA )
00346             CALL CLACGV( I-1, X( I, 1 ), LDX )
00347             CALL CGEMV( 'Conjugate transpose', I-1, N-I+1, -ONE,
00348      $                  A( 1, I ), LDA, X( I, 1 ), LDX, ONE, A( I, I ),
00349      $                  LDA )
00350             CALL CLACGV( I-1, X( I, 1 ), LDX )
00351 *
00352 *           Generate reflection P(i) to annihilate A(i,i+1:n)
00353 *
00354             ALPHA = A( I, I )
00355             CALL CLARFG( N-I+1, ALPHA, A( I, MIN( I+1, N ) ), LDA,
00356      $                   TAUP( I ) )
00357             D( I ) = ALPHA
00358             IF( I.LT.M ) THEN
00359                A( I, I ) = ONE
00360 *
00361 *              Compute X(i+1:m,i)
00362 *
00363                CALL CGEMV( 'No transpose', M-I, N-I+1, ONE, A( I+1, I ),
00364      $                     LDA, A( I, I ), LDA, ZERO, X( I+1, I ), 1 )
00365                CALL CGEMV( 'Conjugate transpose', N-I+1, I-1, ONE,
00366      $                     Y( I, 1 ), LDY, A( I, I ), LDA, ZERO,
00367      $                     X( 1, I ), 1 )
00368                CALL CGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
00369      $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
00370                CALL CGEMV( 'No transpose', I-1, N-I+1, ONE, A( 1, I ),
00371      $                     LDA, A( I, I ), LDA, ZERO, X( 1, I ), 1 )
00372                CALL CGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
00373      $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
00374                CALL CSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
00375                CALL CLACGV( N-I+1, A( I, I ), LDA )
00376 *
00377 *              Update A(i+1:m,i)
00378 *
00379                CALL CLACGV( I-1, Y( I, 1 ), LDY )
00380                CALL CGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
00381      $                     LDA, Y( I, 1 ), LDY, ONE, A( I+1, I ), 1 )
00382                CALL CLACGV( I-1, Y( I, 1 ), LDY )
00383                CALL CGEMV( 'No transpose', M-I, I, -ONE, X( I+1, 1 ),
00384      $                     LDX, A( 1, I ), 1, ONE, A( I+1, I ), 1 )
00385 *
00386 *              Generate reflection Q(i) to annihilate A(i+2:m,i)
00387 *
00388                ALPHA = A( I+1, I )
00389                CALL CLARFG( M-I, ALPHA, A( MIN( I+2, M ), I ), 1,
00390      $                      TAUQ( I ) )
00391                E( I ) = ALPHA
00392                A( I+1, I ) = ONE
00393 *
00394 *              Compute Y(i+1:n,i)
00395 *
00396                CALL CGEMV( 'Conjugate transpose', M-I, N-I, ONE,
00397      $                     A( I+1, I+1 ), LDA, A( I+1, I ), 1, ZERO,
00398      $                     Y( I+1, I ), 1 )
00399                CALL CGEMV( 'Conjugate transpose', M-I, I-1, ONE,
00400      $                     A( I+1, 1 ), LDA, A( I+1, I ), 1, ZERO,
00401      $                     Y( 1, I ), 1 )
00402                CALL CGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
00403      $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
00404                CALL CGEMV( 'Conjugate transpose', M-I, I, ONE,
00405      $                     X( I+1, 1 ), LDX, A( I+1, I ), 1, ZERO,
00406      $                     Y( 1, I ), 1 )
00407                CALL CGEMV( 'Conjugate transpose', I, N-I, -ONE,
00408      $                     A( 1, I+1 ), LDA, Y( 1, I ), 1, ONE,
00409      $                     Y( I+1, I ), 1 )
00410                CALL CSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
00411             ELSE
00412                CALL CLACGV( N-I+1, A( I, I ), LDA )
00413             END IF
00414    20    CONTINUE
00415       END IF
00416       RETURN
00417 *
00418 *     End of CLABRD
00419 *
00420       END
 All Files Functions