SRC/dsytrd_2stage.f

   1 *> \brief \b DSYTRD_2STAGE
   2 *
   3 *  @generated from zhetrd_2stage.f, fortran z -> d, Sun Nov  6 19:34:06 2016
   4 *
   5 *  =========== DOCUMENTATION ===========
   6 *
   7 * Online html documentation available at
   8 *            http://www.netlib.org/lapack/explore-html/
   9 *
  10 *> \htmlonly
  11 *> Download DSYTRD_2STAGE + dependencies
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrd_2stage.f">
  13 *> [TGZ]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrd_2stage.f">
  15 *> [ZIP]</a>
  16 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrd_2stage.f">
  17 *> [TXT]</a>
  18 *> \endhtmlonly
  19 *
  20 *  Definition:
  21 *  ===========
  22 *
  23 *       SUBROUTINE DSYTRD_2STAGE( VECT, UPLO, N, A, LDA, D, E, TAU,
  24 *                                 HOUS2, LHOUS2, WORK, LWORK, INFO )
  25 *
  26 *       IMPLICIT NONE
  27 *
  28 *      .. Scalar Arguments ..
  29 *       CHARACTER          VECT, UPLO
  30 *       INTEGER            N, LDA, LWORK, LHOUS2, INFO
  31 *      ..
  32 *      .. Array Arguments ..
  33 *       DOUBLE PRECISION   D( * ), E( * )
  34 *       DOUBLE PRECISION   A( LDA, * ), TAU( * ),
  35 *                          HOUS2( * ), WORK( * )
  36 *       ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *> DSYTRD_2STAGE reduces a real symmetric matrix A to real symmetric
  45 *> tridiagonal form T by a orthogonal similarity transformation:
  46 *> Q1**T Q2**T* A * Q2 * Q1 = T.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] VECT
  53 *> \verbatim
  54 *>          VECT is CHARACTER*1
  55 *>          = 'N':  No need for the Housholder representation,
  56 *>                  in particular for the second stage (Band to
  57 *>                  tridiagonal) and thus LHOUS2 is of size max(1, 4*N);
  58 *>          = 'V':  the Householder representation is needed to
  59 *>                  either generate Q1 Q2 or to apply Q1 Q2,
  60 *>                  then LHOUS2 is to be queried and computed.
  61 *>                  (NOT AVAILABLE IN THIS RELEASE).
  62 *> \endverbatim
  63 *>
  64 *> \param[in] UPLO
  65 *> \verbatim
  66 *>          UPLO is CHARACTER*1
  67 *>          = 'U':  Upper triangle of A is stored;
  68 *>          = 'L':  Lower triangle of A is stored.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] N
  72 *> \verbatim
  73 *>          N is INTEGER
  74 *>          The order of the matrix A.  N >= 0.
  75 *> \endverbatim
  76 *>
  77 *> \param[in,out] A
  78 *> \verbatim
  79 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  80 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
  81 *>          N-by-N upper triangular part of A contains the upper
  82 *>          triangular part of the matrix A, and the strictly lower
  83 *>          triangular part of A is not referenced.  If UPLO = 'L', the
  84 *>          leading N-by-N lower triangular part of A contains the lower
  85 *>          triangular part of the matrix A, and the strictly upper
  86 *>          triangular part of A is not referenced.
  87 *>          On exit, if UPLO = 'U', the band superdiagonal
  88 *>          of A are overwritten by the corresponding elements of the
  89 *>          internal band-diagonal matrix AB, and the elements above
  90 *>          the KD superdiagonal, with the array TAU, represent the orthogonal
  91 *>          matrix Q1 as a product of elementary reflectors; if UPLO
  92 *>          = 'L', the diagonal and band subdiagonal of A are over-
  93 *>          written by the corresponding elements of the internal band-diagonal
  94 *>          matrix AB, and the elements below the KD subdiagonal, with
  95 *>          the array TAU, represent the orthogonal matrix Q1 as a product
  96 *>          of elementary reflectors. See Further Details.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] LDA
 100 *> \verbatim
 101 *>          LDA is INTEGER
 102 *>          The leading dimension of the array A.  LDA >= max(1,N).
 103 *> \endverbatim
 104 *>
 105 *> \param[out] D
 106 *> \verbatim
 107 *>          D is DOUBLE PRECISION array, dimension (N)
 108 *>          The diagonal elements of the tridiagonal matrix T.
 109 *> \endverbatim
 110 *>
 111 *> \param[out] E
 112 *> \verbatim
 113 *>          E is DOUBLE PRECISION array, dimension (N-1)
 114 *>          The off-diagonal elements of the tridiagonal matrix T.
 115 *> \endverbatim
 116 *>
 117 *> \param[out] TAU
 118 *> \verbatim
 119 *>          TAU is DOUBLE PRECISION array, dimension (N-KD)
 120 *>          The scalar factors of the elementary reflectors of
 121 *>          the first stage (see Further Details).
 122 *> \endverbatim
 123 *>
 124 *> \param[out] HOUS2
 125 *> \verbatim
 126 *>          HOUS2 is DOUBLE PRECISION array, dimension LHOUS2, that
 127 *>          store the Householder representation of the stage2
 128 *>          band to tridiagonal.
 129 *> \endverbatim
 130 *>
 131 *> \param[in] LHOUS2
 132 *> \verbatim
 133 *>          LHOUS2 is INTEGER
 134 *>          The dimension of the array HOUS2. LHOUS2 = MAX(1, dimension)
 135 *>          If LWORK = -1, or LHOUS2=-1,
 136 *>          then a query is assumed; the routine
 137 *>          only calculates the optimal size of the HOUS2 array, returns
 138 *>          this value as the first entry of the HOUS2 array, and no error
 139 *>          message related to LHOUS2 is issued by XERBLA.
 140 *>          LHOUS2 = MAX(1, dimension) where
 141 *>          dimension = 4*N if VECT='N'
 142 *>          not available now if VECT='H'
 143 *> \endverbatim
 144 *>
 145 *> \param[out] WORK
 146 *> \verbatim
 147 *>          WORK is DOUBLE PRECISION array, dimension LWORK.
 148 *> \endverbatim
 149 *>
 150 *> \param[in] LWORK
 151 *> \verbatim
 152 *>          LWORK is INTEGER
 153 *>          The dimension of the array WORK. LWORK = MAX(1, dimension)
 154 *>          If LWORK = -1, or LHOUS2=-1,
 155 *>          then a workspace query is assumed; the routine
 156 *>          only calculates the optimal size of the WORK array, returns
 157 *>          this value as the first entry of the WORK array, and no error
 158 *>          message related to LWORK is issued by XERBLA.
 159 *>          LWORK = MAX(1, dimension) where
 160 *>          dimension   = max(stage1,stage2) + (KD+1)*N
 161 *>                      = N*KD + N*max(KD+1,FACTOPTNB)
 162 *>                        + max(2*KD*KD, KD*NTHREADS)
 163 *>                        + (KD+1)*N
 164 *>          where KD is the blocking size of the reduction,
 165 *>          FACTOPTNB is the blocking used by the QR or LQ
 166 *>          algorithm, usually FACTOPTNB=128 is a good choice
 167 *>          NTHREADS is the number of threads used when
 168 *>          openMP compilation is enabled, otherwise =1.
 169 *> \endverbatim
 170 *>
 171 *> \param[out] INFO
 172 *> \verbatim
 173 *>          INFO is INTEGER
 174 *>          = 0:  successful exit
 175 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 176 *> \endverbatim
 177 *
 178 *  Authors:
 179 *  ========
 180 *
 181 *> \author Univ. of Tennessee
 182 *> \author Univ. of California Berkeley
 183 *> \author Univ. of Colorado Denver
 184 *> \author NAG Ltd.
 185 *
 186 *> \date November 2016
 187 *
 188 *> \ingroup doubleSYcomputational
 189 *
 190 *> \par Further Details:
 191 *  =====================
 192 *>
 193 *> \verbatim
 194 *>
 195 *>  Implemented by Azzam Haidar.
 196 *>
 197 *>  All details are available on technical report, SC11, SC13 papers.
 198 *>
 199 *>  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
 200 *>  Parallel reduction to condensed forms for symmetric eigenvalue problems
 201 *>  using aggregated fine-grained and memory-aware kernels. In Proceedings
 202 *>  of 2011 International Conference for High Performance Computing,
 203 *>  Networking, Storage and Analysis (SC '11), New York, NY, USA,
 204 *>  Article 8 , 11 pages.
 205 *>  http://doi.acm.org/10.1145/2063384.2063394
 206 *>
 207 *>  A. Haidar, J. Kurzak, P. Luszczek, 2013.
 208 *>  An improved parallel singular value algorithm and its implementation
 209 *>  for multicore hardware, In Proceedings of 2013 International Conference
 210 *>  for High Performance Computing, Networking, Storage and Analysis (SC '13).
 211 *>  Denver, Colorado, USA, 2013.
 212 *>  Article 90, 12 pages.
 213 *>  http://doi.acm.org/10.1145/2503210.2503292
 214 *>
 215 *>  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
 216 *>  A novel hybrid CPU-GPU generalized eigensolver for electronic structure
 217 *>  calculations based on fine-grained memory aware tasks.
 218 *>  International Journal of High Performance Computing Applications.
 219 *>  Volume 28 Issue 2, Pages 196-209, May 2014.
 220 *>  http://hpc.sagepub.com/content/28/2/196
 221 *>
 222 *> \endverbatim
 223 *>
 224 *  =====================================================================
 225       SUBROUTINE DSYTRD_2STAGE( VECT, UPLO, N, A, LDA, D, E, TAU,
 226      $                          HOUS2, LHOUS2, WORK, LWORK, INFO )
 227 *
 228       IMPLICIT NONE
 229 *
 230 *  -- LAPACK computational routine (version 3.4.0) --
 231 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 232 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 233 *     November 2016
 234 *
 235 *     .. Scalar Arguments ..
 236       CHARACTER          VECT, UPLO
 237       INTEGER            N, LDA, LWORK, LHOUS2, INFO
 238 *     ..
 239 *     .. Array Arguments ..
 240       DOUBLE PRECISION   D( * ), E( * )
 241       DOUBLE PRECISION   A( LDA, * ), TAU( * ),
 242      $                   HOUS2( * ), WORK( * )
 243 *     ..
 244 *
 245 *  =====================================================================
 246 *     ..
 247 *     .. Local Scalars ..
 248       LOGICAL            LQUERY, UPPER, WANTQ
 249       INTEGER            KD, IB, LWMIN, LHMIN, LWRK, LDAB, WPOS, ABPOS
 250 *     ..
 251 *     .. External Subroutines ..
 252       EXTERNAL           XERBLA, DSYTRD_SY2SB, DSYTRD_SB2ST
 253 *     ..
 254 *     .. External Functions ..
 255       LOGICAL            LSAME
 256       INTEGER            ILAENV
 257       EXTERNAL           LSAME, ILAENV
 258 *     ..
 259 *     .. Executable Statements ..
 260 *
 261 *     Test the input parameters
 262 *
 263       INFO   = 0
 264       WANTQ  = LSAME( VECT, 'V' )
 265       UPPER  = LSAME( UPLO, 'U' )
 266       LQUERY = ( LWORK.EQ.-1 ) .OR. ( LHOUS2.EQ.-1 )
 267 *
 268 *     Determine the block size, the workspace size and the hous size.
 269 *
 270       KD     = ILAENV( 17, 'DSYTRD_2STAGE', VECT, N, -1, -1, -1 )
 271       IB     = ILAENV( 18, 'DSYTRD_2STAGE', VECT, N, KD, -1, -1 )
 272       LHMIN  = ILAENV( 19, 'DSYTRD_2STAGE', VECT, N, KD, IB, -1 )
 273       LWMIN  = ILAENV( 20, 'DSYTRD_2STAGE', VECT, N, KD, IB, -1 )
 274 *      WRITE(*,*),'DSYTRD_2STAGE N KD UPLO LHMIN LWMIN ',N, KD, UPLO,
 275 *     $            LHMIN, LWMIN
 276 *
 277       IF( .NOT.LSAME( VECT, 'N' ) ) THEN
 278          INFO = -1
 279       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 280          INFO = -2
 281       ELSE IF( N.LT.0 ) THEN
 282          INFO = -3
 283       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 284          INFO = -5
 285       ELSE IF( LHOUS2.LT.LHMIN .AND. .NOT.LQUERY ) THEN
 286          INFO = -10
 287       ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 288          INFO = -12
 289       END IF
 290 *
 291       IF( INFO.EQ.0 ) THEN
 292          HOUS2( 1 ) = LHMIN
 293          WORK( 1 )  = LWMIN
 294       END IF
 295 *
 296       IF( INFO.NE.0 ) THEN
 297          CALL XERBLA( 'DSYTRD_2STAGE', -INFO )
 298          RETURN
 299       ELSE IF( LQUERY ) THEN
 300          RETURN
 301       END IF
 302 *
 303 *     Quick return if possible
 304 *
 305       IF( N.EQ.0 ) THEN
 306          WORK( 1 ) = 1
 307          RETURN
 308       END IF
 309 *
 310 *     Determine pointer position
 311 *
 312       LDAB  = KD+1
 313       LWRK  = LWORK-LDAB*N
 314       ABPOS = 1
 315       WPOS  = ABPOS + LDAB*N
 316       CALL DSYTRD_SY2SB( UPLO, N, KD, A, LDA, WORK( ABPOS ), LDAB,
 317      $                   TAU, WORK( WPOS ), LWRK, INFO )
 318       IF( INFO.NE.0 ) THEN
 319          CALL XERBLA( 'DSYTRD_SY2SB', -INFO )
 320          RETURN
 321       END IF
 322       CALL DSYTRD_SB2ST( 'Y', VECT, UPLO, N, KD,
 323      $                   WORK( ABPOS ), LDAB, D, E,
 324      $                   HOUS2, LHOUS2, WORK( WPOS ), LWRK, INFO )
 325       IF( INFO.NE.0 ) THEN
 326          CALL XERBLA( 'DSYTRD_SB2ST', -INFO )
 327          RETURN
 328       END IF
 329 *
 330 *
 331       HOUS2( 1 ) = LHMIN
 332       WORK( 1 )  = LWMIN
 333       RETURN
 334 *
 335 *     End of DSYTRD_2STAGE
 336 *
 337       END