File:  [local] / rpl / lapack / lapack / zhetrd_2stage.f
Revision 1.6: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:24 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    1: *> \brief \b ZHETRD_2STAGE
    2: *
    3: *  @precisions fortran z -> s d c
    4: *
    5: *  =========== DOCUMENTATION ===========
    6: *
    7: * Online html documentation available at 
    8: *            http://www.netlib.org/lapack/explore-html/ 
    9: *
   10: *> \htmlonly
   11: *> Download ZHETRD_2STAGE + dependencies 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrd_2stage.f"> 
   13: *> [TGZ]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrd_2stage.f"> 
   15: *> [ZIP]</a> 
   16: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrd_2stage.f"> 
   17: *> [TXT]</a>
   18: *> \endhtmlonly 
   19: *
   20: *  Definition:
   21: *  ===========
   22: *
   23: *       SUBROUTINE ZHETRD_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: *       COMPLEX*16         A( LDA, * ), TAU( * ),
   35: *                          HOUS2( * ), WORK( * )
   36: *       ..
   37: *  
   38: *
   39: *> \par Purpose:
   40: *  =============
   41: *>
   42: *> \verbatim
   43: *>
   44: *> ZHETRD_2STAGE reduces a complex Hermitian matrix A to real symmetric
   45: *> tridiagonal form T by a unitary similarity transformation:
   46: *> Q1**H Q2**H* 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 COMPLEX*16 array, dimension (LDA,N)
   80: *>          On entry, the Hermitian 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 unitary
   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 unitary 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LHOUS2)
  127: *>          Stores 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.
  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: *>          If VECT='N', LHOUS2 = max(1, 4*n);
  141: *>          if VECT='V', option not yet available.
  142: *> \endverbatim
  143: *>
  144: *> \param[out] WORK
  145: *> \verbatim
  146: *>          WORK is COMPLEX*16 array, dimension (LWORK)
  147: *> \endverbatim
  148: *>
  149: *> \param[in] LWORK
  150: *> \verbatim
  151: *>          LWORK is INTEGER
  152: *>          The dimension of the array WORK. LWORK = MAX(1, dimension)
  153: *>          If LWORK = -1, or LHOUS2=-1,
  154: *>          then a workspace query is assumed; the routine
  155: *>          only calculates the optimal size of the WORK array, returns
  156: *>          this value as the first entry of the WORK array, and no error
  157: *>          message related to LWORK is issued by XERBLA.
  158: *>          LWORK = MAX(1, dimension) where
  159: *>          dimension   = max(stage1,stage2) + (KD+1)*N
  160: *>                      = N*KD + N*max(KD+1,FACTOPTNB) 
  161: *>                        + max(2*KD*KD, KD*NTHREADS) 
  162: *>                        + (KD+1)*N 
  163: *>          where KD is the blocking size of the reduction,
  164: *>          FACTOPTNB is the blocking used by the QR or LQ
  165: *>          algorithm, usually FACTOPTNB=128 is a good choice
  166: *>          NTHREADS is the number of threads used when
  167: *>          openMP compilation is enabled, otherwise =1.
  168: *> \endverbatim
  169: *>
  170: *> \param[out] INFO
  171: *> \verbatim
  172: *>          INFO is INTEGER
  173: *>          = 0:  successful exit
  174: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  175: *> \endverbatim
  176: *
  177: *  Authors:
  178: *  ========
  179: *
  180: *> \author Univ. of Tennessee 
  181: *> \author Univ. of California Berkeley 
  182: *> \author Univ. of Colorado Denver 
  183: *> \author NAG Ltd. 
  184: *
  185: *> \ingroup complex16HEcomputational
  186: *
  187: *> \par Further Details:
  188: *  =====================
  189: *>
  190: *> \verbatim
  191: *>
  192: *>  Implemented by Azzam Haidar.
  193: *>
  194: *>  All details are available on technical report, SC11, SC13 papers.
  195: *>
  196: *>  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
  197: *>  Parallel reduction to condensed forms for symmetric eigenvalue problems
  198: *>  using aggregated fine-grained and memory-aware kernels. In Proceedings
  199: *>  of 2011 International Conference for High Performance Computing,
  200: *>  Networking, Storage and Analysis (SC '11), New York, NY, USA,
  201: *>  Article 8 , 11 pages.
  202: *>  http://doi.acm.org/10.1145/2063384.2063394
  203: *>
  204: *>  A. Haidar, J. Kurzak, P. Luszczek, 2013.
  205: *>  An improved parallel singular value algorithm and its implementation 
  206: *>  for multicore hardware, In Proceedings of 2013 International Conference
  207: *>  for High Performance Computing, Networking, Storage and Analysis (SC '13).
  208: *>  Denver, Colorado, USA, 2013.
  209: *>  Article 90, 12 pages.
  210: *>  http://doi.acm.org/10.1145/2503210.2503292
  211: *>
  212: *>  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
  213: *>  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 
  214: *>  calculations based on fine-grained memory aware tasks.
  215: *>  International Journal of High Performance Computing Applications.
  216: *>  Volume 28 Issue 2, Pages 196-209, May 2014.
  217: *>  http://hpc.sagepub.com/content/28/2/196 
  218: *>
  219: *> \endverbatim
  220: *>
  221: *  =====================================================================
  222:       SUBROUTINE ZHETRD_2STAGE( VECT, UPLO, N, A, LDA, D, E, TAU, 
  223:      $                          HOUS2, LHOUS2, WORK, LWORK, INFO )
  224: *
  225:       IMPLICIT NONE
  226: *
  227: *  -- LAPACK computational routine --
  228: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  229: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  230: *
  231: *     .. Scalar Arguments ..
  232:       CHARACTER          VECT, UPLO
  233:       INTEGER            N, LDA, LWORK, LHOUS2, INFO
  234: *     ..
  235: *     .. Array Arguments ..
  236:       DOUBLE PRECISION   D( * ), E( * )
  237:       COMPLEX*16         A( LDA, * ), TAU( * ),
  238:      $                   HOUS2( * ), WORK( * )
  239: *     ..
  240: *
  241: *  =====================================================================
  242: *     ..
  243: *     .. Local Scalars ..
  244:       LOGICAL            LQUERY, UPPER, WANTQ
  245:       INTEGER            KD, IB, LWMIN, LHMIN, LWRK, LDAB, WPOS, ABPOS
  246: *     ..
  247: *     .. External Subroutines ..
  248:       EXTERNAL           XERBLA, ZHETRD_HE2HB, ZHETRD_HB2ST
  249: *     ..
  250: *     .. External Functions ..
  251:       LOGICAL            LSAME
  252:       INTEGER            ILAENV2STAGE
  253:       EXTERNAL           LSAME, ILAENV2STAGE
  254: *     ..
  255: *     .. Executable Statements ..
  256: *
  257: *     Test the input parameters
  258: *
  259:       INFO   = 0
  260:       WANTQ  = LSAME( VECT, 'V' )
  261:       UPPER  = LSAME( UPLO, 'U' )
  262:       LQUERY = ( LWORK.EQ.-1 ) .OR. ( LHOUS2.EQ.-1 )
  263: *
  264: *     Determine the block size, the workspace size and the hous size.
  265: *
  266:       KD     = ILAENV2STAGE( 1, 'ZHETRD_2STAGE', VECT, N, -1, -1, -1 )
  267:       IB     = ILAENV2STAGE( 2, 'ZHETRD_2STAGE', VECT, N, KD, -1, -1 )
  268:       LHMIN  = ILAENV2STAGE( 3, 'ZHETRD_2STAGE', VECT, N, KD, IB, -1 )
  269:       LWMIN  = ILAENV2STAGE( 4, 'ZHETRD_2STAGE', VECT, N, KD, IB, -1 )
  270: *      WRITE(*,*),'ZHETRD_2STAGE N KD UPLO LHMIN LWMIN ',N, KD, UPLO,
  271: *     $            LHMIN, LWMIN
  272: *
  273:       IF( .NOT.LSAME( VECT, 'N' ) ) THEN
  274:          INFO = -1
  275:       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  276:          INFO = -2
  277:       ELSE IF( N.LT.0 ) THEN
  278:          INFO = -3
  279:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  280:          INFO = -5
  281:       ELSE IF( LHOUS2.LT.LHMIN .AND. .NOT.LQUERY ) THEN
  282:          INFO = -10
  283:       ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
  284:          INFO = -12
  285:       END IF
  286: *
  287:       IF( INFO.EQ.0 ) THEN
  288:          HOUS2( 1 ) = LHMIN
  289:          WORK( 1 )  = LWMIN
  290:       END IF
  291: *
  292:       IF( INFO.NE.0 ) THEN
  293:          CALL XERBLA( 'ZHETRD_2STAGE', -INFO )
  294:          RETURN
  295:       ELSE IF( LQUERY ) THEN
  296:          RETURN
  297:       END IF
  298: *
  299: *     Quick return if possible
  300: *
  301:       IF( N.EQ.0 ) THEN
  302:          WORK( 1 ) = 1
  303:          RETURN
  304:       END IF
  305: *
  306: *     Determine pointer position
  307: *
  308:       LDAB  = KD+1
  309:       LWRK  = LWORK-LDAB*N
  310:       ABPOS = 1
  311:       WPOS  = ABPOS + LDAB*N
  312:       CALL ZHETRD_HE2HB( UPLO, N, KD, A, LDA, WORK( ABPOS ), LDAB, 
  313:      $                   TAU, WORK( WPOS ), LWRK, INFO )
  314:       IF( INFO.NE.0 ) THEN
  315:          CALL XERBLA( 'ZHETRD_HE2HB', -INFO )
  316:          RETURN
  317:       END IF
  318:       CALL ZHETRD_HB2ST( 'Y', VECT, UPLO, N, KD, 
  319:      $                   WORK( ABPOS ), LDAB, D, E, 
  320:      $                   HOUS2, LHOUS2, WORK( WPOS ), LWRK, INFO )
  321:       IF( INFO.NE.0 ) THEN
  322:          CALL XERBLA( 'ZHETRD_HB2ST', -INFO )
  323:          RETURN
  324:       END IF
  325: *
  326: *
  327:       HOUS2( 1 ) = LHMIN
  328:       WORK( 1 )  = LWMIN
  329:       RETURN
  330: *
  331: *     End of ZHETRD_2STAGE
  332: *
  333:       END

CVSweb interface <joel.bertrand@systella.fr>