LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sgehrd.f
Go to the documentation of this file.
00001 *> \brief \b SGEHRD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SGEHRD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgehrd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgehrd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgehrd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            IHI, ILO, INFO, LDA, LWORK, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       REAL              A( LDA, * ), TAU( * ), WORK( * )
00028 *       ..
00029 *  
00030 *
00031 *> \par Purpose:
00032 *  =============
00033 *>
00034 *> \verbatim
00035 *>
00036 *> SGEHRD reduces a real general matrix A to upper Hessenberg form H by
00037 *> an orthogonal similarity transformation:  Q**T * A * Q = H .
00038 *> \endverbatim
00039 *
00040 *  Arguments:
00041 *  ==========
00042 *
00043 *> \param[in] N
00044 *> \verbatim
00045 *>          N is INTEGER
00046 *>          The order of the matrix A.  N >= 0.
00047 *> \endverbatim
00048 *>
00049 *> \param[in] ILO
00050 *> \verbatim
00051 *>          ILO is INTEGER
00052 *> \endverbatim
00053 *>
00054 *> \param[in] IHI
00055 *> \verbatim
00056 *>          IHI is INTEGER
00057 *>
00058 *>          It is assumed that A is already upper triangular in rows
00059 *>          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
00060 *>          set by a previous call to SGEBAL; otherwise they should be
00061 *>          set to 1 and N respectively. See Further Details.
00062 *>          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
00063 *> \endverbatim
00064 *>
00065 *> \param[in,out] A
00066 *> \verbatim
00067 *>          A is REAL array, dimension (LDA,N)
00068 *>          On entry, the N-by-N general matrix to be reduced.
00069 *>          On exit, the upper triangle and the first subdiagonal of A
00070 *>          are overwritten with the upper Hessenberg matrix H, and the
00071 *>          elements below the first subdiagonal, with the array TAU,
00072 *>          represent the orthogonal matrix Q as a product of elementary
00073 *>          reflectors. See Further Details.
00074 *> \endverbatim
00075 *>
00076 *> \param[in] LDA
00077 *> \verbatim
00078 *>          LDA is INTEGER
00079 *>          The leading dimension of the array A.  LDA >= max(1,N).
00080 *> \endverbatim
00081 *>
00082 *> \param[out] TAU
00083 *> \verbatim
00084 *>          TAU is REAL array, dimension (N-1)
00085 *>          The scalar factors of the elementary reflectors (see Further
00086 *>          Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to
00087 *>          zero.
00088 *> \endverbatim
00089 *>
00090 *> \param[out] WORK
00091 *> \verbatim
00092 *>          WORK is REAL array, dimension (LWORK)
00093 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00094 *> \endverbatim
00095 *>
00096 *> \param[in] LWORK
00097 *> \verbatim
00098 *>          LWORK is INTEGER
00099 *>          The length of the array WORK.  LWORK >= max(1,N).
00100 *>          For optimum performance LWORK >= N*NB, where NB is the
00101 *>          optimal blocksize.
00102 *>
00103 *>          If LWORK = -1, then a workspace query is assumed; the routine
00104 *>          only calculates the optimal size of the WORK array, returns
00105 *>          this value as the first entry of the WORK array, and no error
00106 *>          message related to LWORK is issued by XERBLA.
00107 *> \endverbatim
00108 *>
00109 *> \param[out] INFO
00110 *> \verbatim
00111 *>          INFO is INTEGER
00112 *>          = 0:  successful exit
00113 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00114 *> \endverbatim
00115 *
00116 *  Authors:
00117 *  ========
00118 *
00119 *> \author Univ. of Tennessee 
00120 *> \author Univ. of California Berkeley 
00121 *> \author Univ. of Colorado Denver 
00122 *> \author NAG Ltd. 
00123 *
00124 *> \date November 2011
00125 *
00126 *> \ingroup realGEcomputational
00127 *
00128 *> \par Further Details:
00129 *  =====================
00130 *>
00131 *> \verbatim
00132 *>
00133 *>  The matrix Q is represented as a product of (ihi-ilo) elementary
00134 *>  reflectors
00135 *>
00136 *>     Q = H(ilo) H(ilo+1) . . . H(ihi-1).
00137 *>
00138 *>  Each H(i) has the form
00139 *>
00140 *>     H(i) = I - tau * v * v**T
00141 *>
00142 *>  where tau is a real scalar, and v is a real vector with
00143 *>  v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
00144 *>  exit in A(i+2:ihi,i), and tau in TAU(i).
00145 *>
00146 *>  The contents of A are illustrated by the following example, with
00147 *>  n = 7, ilo = 2 and ihi = 6:
00148 *>
00149 *>  on entry,                        on exit,
00150 *>
00151 *>  ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
00152 *>  (     a   a   a   a   a   a )    (      a   h   h   h   h   a )
00153 *>  (     a   a   a   a   a   a )    (      h   h   h   h   h   h )
00154 *>  (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
00155 *>  (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
00156 *>  (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
00157 *>  (                         a )    (                          a )
00158 *>
00159 *>  where a denotes an element of the original matrix A, h denotes a
00160 *>  modified element of the upper Hessenberg matrix H, and vi denotes an
00161 *>  element of the vector defining H(i).
00162 *>
00163 *>  This file is a slight modification of LAPACK-3.0's DGEHRD
00164 *>  subroutine incorporating improvements proposed by Quintana-Orti and
00165 *>  Van de Geijn (2006). (See DLAHR2.)
00166 *> \endverbatim
00167 *>
00168 *  =====================================================================
00169       SUBROUTINE SGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
00170 *
00171 *  -- LAPACK computational routine (version 3.4.0) --
00172 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00173 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00174 *     November 2011
00175 *
00176 *     .. Scalar Arguments ..
00177       INTEGER            IHI, ILO, INFO, LDA, LWORK, N
00178 *     ..
00179 *     .. Array Arguments ..
00180       REAL              A( LDA, * ), TAU( * ), WORK( * )
00181 *     ..
00182 *
00183 *  =====================================================================
00184 *
00185 *     .. Parameters ..
00186       INTEGER            NBMAX, LDT
00187       PARAMETER          ( NBMAX = 64, LDT = NBMAX+1 )
00188       REAL              ZERO, ONE
00189       PARAMETER          ( ZERO = 0.0E+0, 
00190      $                     ONE = 1.0E+0 )
00191 *     ..
00192 *     .. Local Scalars ..
00193       LOGICAL            LQUERY
00194       INTEGER            I, IB, IINFO, IWS, J, LDWORK, LWKOPT, NB,
00195      $                   NBMIN, NH, NX
00196       REAL              EI
00197 *     ..
00198 *     .. Local Arrays ..
00199       REAL              T( LDT, NBMAX )
00200 *     ..
00201 *     .. External Subroutines ..
00202       EXTERNAL           SAXPY, SGEHD2, SGEMM, SLAHR2, SLARFB, STRMM,
00203      $                   XERBLA
00204 *     ..
00205 *     .. Intrinsic Functions ..
00206       INTRINSIC          MAX, MIN
00207 *     ..
00208 *     .. External Functions ..
00209       INTEGER            ILAENV
00210       EXTERNAL           ILAENV
00211 *     ..
00212 *     .. Executable Statements ..
00213 *
00214 *     Test the input parameters
00215 *
00216       INFO = 0
00217       NB = MIN( NBMAX, ILAENV( 1, 'SGEHRD', ' ', N, ILO, IHI, -1 ) )
00218       LWKOPT = N*NB
00219       WORK( 1 ) = LWKOPT
00220       LQUERY = ( LWORK.EQ.-1 )
00221       IF( N.LT.0 ) THEN
00222          INFO = -1
00223       ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
00224          INFO = -2
00225       ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
00226          INFO = -3
00227       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00228          INFO = -5
00229       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
00230          INFO = -8
00231       END IF
00232       IF( INFO.NE.0 ) THEN
00233          CALL XERBLA( 'SGEHRD', -INFO )
00234          RETURN
00235       ELSE IF( LQUERY ) THEN
00236          RETURN
00237       END IF
00238 *
00239 *     Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
00240 *
00241       DO 10 I = 1, ILO - 1
00242          TAU( I ) = ZERO
00243    10 CONTINUE
00244       DO 20 I = MAX( 1, IHI ), N - 1
00245          TAU( I ) = ZERO
00246    20 CONTINUE
00247 *
00248 *     Quick return if possible
00249 *
00250       NH = IHI - ILO + 1
00251       IF( NH.LE.1 ) THEN
00252          WORK( 1 ) = 1
00253          RETURN
00254       END IF
00255 *
00256 *     Determine the block size
00257 *
00258       NB = MIN( NBMAX, ILAENV( 1, 'SGEHRD', ' ', N, ILO, IHI, -1 ) )
00259       NBMIN = 2
00260       IWS = 1
00261       IF( NB.GT.1 .AND. NB.LT.NH ) THEN
00262 *
00263 *        Determine when to cross over from blocked to unblocked code
00264 *        (last block is always handled by unblocked code)
00265 *
00266          NX = MAX( NB, ILAENV( 3, 'SGEHRD', ' ', N, ILO, IHI, -1 ) )
00267          IF( NX.LT.NH ) THEN
00268 *
00269 *           Determine if workspace is large enough for blocked code
00270 *
00271             IWS = N*NB
00272             IF( LWORK.LT.IWS ) THEN
00273 *
00274 *              Not enough workspace to use optimal NB:  determine the
00275 *              minimum value of NB, and reduce NB or force use of
00276 *              unblocked code
00277 *
00278                NBMIN = MAX( 2, ILAENV( 2, 'SGEHRD', ' ', N, ILO, IHI,
00279      $                 -1 ) )
00280                IF( LWORK.GE.N*NBMIN ) THEN
00281                   NB = LWORK / N
00282                ELSE
00283                   NB = 1
00284                END IF
00285             END IF
00286          END IF
00287       END IF
00288       LDWORK = N
00289 *
00290       IF( NB.LT.NBMIN .OR. NB.GE.NH ) THEN
00291 *
00292 *        Use unblocked code below
00293 *
00294          I = ILO
00295 *
00296       ELSE
00297 *
00298 *        Use blocked code
00299 *
00300          DO 40 I = ILO, IHI - 1 - NX, NB
00301             IB = MIN( NB, IHI-I )
00302 *
00303 *           Reduce columns i:i+ib-1 to Hessenberg form, returning the
00304 *           matrices V and T of the block reflector H = I - V*T*V**T
00305 *           which performs the reduction, and also the matrix Y = A*V*T
00306 *
00307             CALL SLAHR2( IHI, I, IB, A( 1, I ), LDA, TAU( I ), T, LDT,
00308      $                   WORK, LDWORK )
00309 *
00310 *           Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
00311 *           right, computing  A := A - Y * V**T. V(i+ib,ib-1) must be set
00312 *           to 1
00313 *
00314             EI = A( I+IB, I+IB-1 )
00315             A( I+IB, I+IB-1 ) = ONE
00316             CALL SGEMM( 'No transpose', 'Transpose', 
00317      $                  IHI, IHI-I-IB+1,
00318      $                  IB, -ONE, WORK, LDWORK, A( I+IB, I ), LDA, ONE,
00319      $                  A( 1, I+IB ), LDA )
00320             A( I+IB, I+IB-1 ) = EI
00321 *
00322 *           Apply the block reflector H to A(1:i,i+1:i+ib-1) from the
00323 *           right
00324 *
00325             CALL STRMM( 'Right', 'Lower', 'Transpose',
00326      $                  'Unit', I, IB-1,
00327      $                  ONE, A( I+1, I ), LDA, WORK, LDWORK )
00328             DO 30 J = 0, IB-2
00329                CALL SAXPY( I, -ONE, WORK( LDWORK*J+1 ), 1,
00330      $                     A( 1, I+J+1 ), 1 )
00331    30       CONTINUE
00332 *
00333 *           Apply the block reflector H to A(i+1:ihi,i+ib:n) from the
00334 *           left
00335 *
00336             CALL SLARFB( 'Left', 'Transpose', 'Forward',
00337      $                   'Columnwise',
00338      $                   IHI-I, N-I-IB+1, IB, A( I+1, I ), LDA, T, LDT,
00339      $                   A( I+1, I+IB ), LDA, WORK, LDWORK )
00340    40    CONTINUE
00341       END IF
00342 *
00343 *     Use unblocked code to reduce the rest of the matrix
00344 *
00345       CALL SGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO )
00346       WORK( 1 ) = IWS
00347 *
00348       RETURN
00349 *
00350 *     End of SGEHRD
00351 *
00352       END
 All Files Functions