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
334