Annotation of rpl/lapack/lapack/zhetrd_2stage.f, revision 1.6

1.1       bertrand    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
1.5       bertrand  126: *>          HOUS2 is COMPLEX*16 array, dimension (LHOUS2)
                    127: *>          Stores the Householder representation of the stage2
1.1       bertrand  128: *>          band to tridiagonal.
                    129: *> \endverbatim
                    130: *>
                    131: *> \param[in] LHOUS2
                    132: *> \verbatim
                    133: *>          LHOUS2 is INTEGER
1.5       bertrand  134: *>          The dimension of the array HOUS2.
                    135: *>          If LWORK = -1, or LHOUS2 = -1,
1.1       bertrand  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.
1.5       bertrand  140: *>          If VECT='N', LHOUS2 = max(1, 4*n);
                    141: *>          if VECT='V', option not yet available.
1.1       bertrand  142: *> \endverbatim
                    143: *>
                    144: *> \param[out] WORK
                    145: *> \verbatim
1.3       bertrand  146: *>          WORK is COMPLEX*16 array, dimension (LWORK)
1.1       bertrand  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: *
1.6     ! bertrand  227: *  -- LAPACK computational routine --
1.1       bertrand  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
1.3       bertrand  252:       INTEGER            ILAENV2STAGE
                    253:       EXTERNAL           LSAME, ILAENV2STAGE
1.1       bertrand  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: *
1.3       bertrand  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 )
1.1       bertrand  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>