LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
chetrd_hb2st.F
Go to the documentation of this file.
1*> \brief \b CHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric tridiagonal form T
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHETRD_HB2ST + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetrd_hb2st.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetrd_hb2st.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetrd_hb2st.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CHETRD_HB2ST( STAGE1, VECT, UPLO, N, KD, AB, LDAB,
22* D, E, HOUS, LHOUS, WORK, LWORK, INFO )
23*
24* #if defined(_OPENMP)
25* use omp_lib
26* #endif
27*
28* IMPLICIT NONE
29*
30* .. Scalar Arguments ..
31* CHARACTER STAGE1, UPLO, VECT
32* INTEGER N, KD, IB, LDAB, LHOUS, LWORK, INFO
33* ..
34* .. Array Arguments ..
35* REAL D( * ), E( * )
36* COMPLEX AB( LDAB, * ), HOUS( * ), WORK( * )
37* ..
38*
39*
40*> \par Purpose:
41* =============
42*>
43*> \verbatim
44*>
45*> CHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric
46*> tridiagonal form T by a unitary similarity transformation:
47*> Q**H * A * Q = T.
48*> \endverbatim
49*
50* Arguments:
51* ==========
52*
53*> \param[in] STAGE1
54*> \verbatim
55*> STAGE1 is CHARACTER*1
56*> = 'N': "No": to mention that the stage 1 of the reduction
57*> from dense to band using the chetrd_he2hb routine
58*> was not called before this routine to reproduce AB.
59*> In other term this routine is called as standalone.
60*> = 'Y': "Yes": to mention that the stage 1 of the
61*> reduction from dense to band using the chetrd_he2hb
62*> routine has been called to produce AB (e.g., AB is
63*> the output of chetrd_he2hb.
64*> \endverbatim
65*>
66*> \param[in] VECT
67*> \verbatim
68*> VECT is CHARACTER*1
69*> = 'N': No need for the Housholder representation,
70*> and thus LHOUS is of size max(1, 4*N);
71*> = 'V': the Householder representation is needed to
72*> either generate or to apply Q later on,
73*> then LHOUS is to be queried and computed.
74*> (NOT AVAILABLE IN THIS RELEASE).
75*> \endverbatim
76*>
77*> \param[in] UPLO
78*> \verbatim
79*> UPLO is CHARACTER*1
80*> = 'U': Upper triangle of A is stored;
81*> = 'L': Lower triangle of A is stored.
82*> \endverbatim
83*>
84*> \param[in] N
85*> \verbatim
86*> N is INTEGER
87*> The order of the matrix A. N >= 0.
88*> \endverbatim
89*>
90*> \param[in] KD
91*> \verbatim
92*> KD is INTEGER
93*> The number of superdiagonals of the matrix A if UPLO = 'U',
94*> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
95*> \endverbatim
96*>
97*> \param[in,out] AB
98*> \verbatim
99*> AB is COMPLEX array, dimension (LDAB,N)
100*> On entry, the upper or lower triangle of the Hermitian band
101*> matrix A, stored in the first KD+1 rows of the array. The
102*> j-th column of A is stored in the j-th column of the array AB
103*> as follows:
104*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
105*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
106*> On exit, the diagonal elements of AB are overwritten by the
107*> diagonal elements of the tridiagonal matrix T; if KD > 0, the
108*> elements on the first superdiagonal (if UPLO = 'U') or the
109*> first subdiagonal (if UPLO = 'L') are overwritten by the
110*> off-diagonal elements of T; the rest of AB is overwritten by
111*> values generated during the reduction.
112*> \endverbatim
113*>
114*> \param[in] LDAB
115*> \verbatim
116*> LDAB is INTEGER
117*> The leading dimension of the array AB. LDAB >= KD+1.
118*> \endverbatim
119*>
120*> \param[out] D
121*> \verbatim
122*> D is REAL array, dimension (N)
123*> The diagonal elements of the tridiagonal matrix T.
124*> \endverbatim
125*>
126*> \param[out] E
127*> \verbatim
128*> E is REAL array, dimension (N-1)
129*> The off-diagonal elements of the tridiagonal matrix T:
130*> E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
131*> \endverbatim
132*>
133*> \param[out] HOUS
134*> \verbatim
135*> HOUS is COMPLEX array, dimension LHOUS, that
136*> store the Householder representation.
137*> \endverbatim
138*>
139*> \param[in] LHOUS
140*> \verbatim
141*> LHOUS is INTEGER
142*> The dimension of the array HOUS. LHOUS = MAX(1, dimension)
143*> If LWORK = -1, or LHOUS=-1,
144*> then a query is assumed; the routine
145*> only calculates the optimal size of the HOUS array, returns
146*> this value as the first entry of the HOUS array, and no error
147*> message related to LHOUS is issued by XERBLA.
148*> LHOUS = MAX(1, dimension) where
149*> dimension = 4*N if VECT='N'
150*> not available now if VECT='H'
151*> \endverbatim
152*>
153*> \param[out] WORK
154*> \verbatim
155*> WORK is COMPLEX array, dimension LWORK.
156*> \endverbatim
157*>
158*> \param[in] LWORK
159*> \verbatim
160*> LWORK is INTEGER
161*> The dimension of the array WORK. LWORK = MAX(1, dimension)
162*> If LWORK = -1, or LHOUS=-1,
163*> then a workspace query is assumed; the routine
164*> only calculates the optimal size of the WORK array, returns
165*> this value as the first entry of the WORK array, and no error
166*> message related to LWORK is issued by XERBLA.
167*> LWORK = MAX(1, dimension) where
168*> dimension = (2KD+1)*N + KD*NTHREADS
169*> where KD is the blocking size of the reduction,
170*> FACTOPTNB is the blocking used by the QR or LQ
171*> algorithm, usually FACTOPTNB=128 is a good choice
172*> NTHREADS is the number of threads used when
173*> openMP compilation is enabled, otherwise =1.
174*> \endverbatim
175*>
176*> \param[out] INFO
177*> \verbatim
178*> INFO is INTEGER
179*> = 0: successful exit
180*> < 0: if INFO = -i, the i-th argument had an illegal value
181*> \endverbatim
182*
183* Authors:
184* ========
185*
186*> \author Univ. of Tennessee
187*> \author Univ. of California Berkeley
188*> \author Univ. of Colorado Denver
189*> \author NAG Ltd.
190*
191*> \ingroup hetrd_hb2st
192*
193*> \par Further Details:
194* =====================
195*>
196*> \verbatim
197*>
198*> Implemented by Azzam Haidar.
199*>
200*> All details are available on technical report, SC11, SC13 papers.
201*>
202*> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
203*> Parallel reduction to condensed forms for symmetric eigenvalue problems
204*> using aggregated fine-grained and memory-aware kernels. In Proceedings
205*> of 2011 International Conference for High Performance Computing,
206*> Networking, Storage and Analysis (SC '11), New York, NY, USA,
207*> Article 8 , 11 pages.
208*> http://doi.acm.org/10.1145/2063384.2063394
209*>
210*> A. Haidar, J. Kurzak, P. Luszczek, 2013.
211*> An improved parallel singular value algorithm and its implementation
212*> for multicore hardware, In Proceedings of 2013 International Conference
213*> for High Performance Computing, Networking, Storage and Analysis (SC '13).
214*> Denver, Colorado, USA, 2013.
215*> Article 90, 12 pages.
216*> http://doi.acm.org/10.1145/2503210.2503292
217*>
218*> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
219*> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
220*> calculations based on fine-grained memory aware tasks.
221*> International Journal of High Performance Computing Applications.
222*> Volume 28 Issue 2, Pages 196-209, May 2014.
223*> http://hpc.sagepub.com/content/28/2/196
224*>
225*> \endverbatim
226*>
227* =====================================================================
228 SUBROUTINE chetrd_hb2st( STAGE1, VECT, UPLO, N, KD, AB, LDAB,
229 $ D, E, HOUS, LHOUS, WORK, LWORK, INFO )
230*
231*
232#if defined(_OPENMP)
233 use omp_lib
234#endif
235*
236 IMPLICIT NONE
237*
238* -- LAPACK computational routine --
239* -- LAPACK is a software package provided by Univ. of Tennessee, --
240* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
241*
242* .. Scalar Arguments ..
243 CHARACTER STAGE1, UPLO, VECT
244 INTEGER N, KD, LDAB, LHOUS, LWORK, INFO
245* ..
246* .. Array Arguments ..
247 REAL D( * ), E( * )
248 COMPLEX AB( LDAB, * ), HOUS( * ), WORK( * )
249* ..
250*
251* =====================================================================
252*
253* .. Parameters ..
254 REAL RZERO
255 COMPLEX ZERO, ONE
256 parameter( rzero = 0.0e+0,
257 $ zero = ( 0.0e+0, 0.0e+0 ),
258 $ one = ( 1.0e+0, 0.0e+0 ) )
259* ..
260* .. Local Scalars ..
261 LOGICAL LQUERY, WANTQ, UPPER, AFTERS1
262 INTEGER I, M, K, IB, SWEEPID, MYID, SHIFT, STT, ST,
263 $ ed, stind, edind, blklastind, colpt, thed,
264 $ stepercol, grsiz, thgrsiz, thgrnb, thgrid,
265 $ nbtiles, ttype, tid, nthreads, debug,
266 $ abdpos, abofdpos, dpos, ofdpos, awpos,
267 $ inda, indw, apos, sizea, lda, indv, indtau,
268 $ sicev, sizetau, ldv, lhmin, lwmin
269 REAL ABSTMP
270 COMPLEX TMP
271* ..
272* .. External Subroutines ..
274* ..
275* .. Intrinsic Functions ..
276 INTRINSIC min, max, ceiling, real
277* ..
278* .. External Functions ..
279 LOGICAL LSAME
280 INTEGER ILAENV2STAGE
281 REAL SROUNDUP_LWORK
282 EXTERNAL lsame, ilaenv2stage, sroundup_lwork
283* ..
284* .. Executable Statements ..
285*
286* Determine the minimal workspace size required.
287* Test the input parameters
288*
289 debug = 0
290 info = 0
291 afters1 = lsame( stage1, 'Y' )
292 wantq = lsame( vect, 'V' )
293 upper = lsame( uplo, 'U' )
294 lquery = ( lwork.EQ.-1 ) .OR. ( lhous.EQ.-1 )
295*
296* Determine the block size, the workspace size and the hous size.
297*
298 ib = ilaenv2stage( 2, 'CHETRD_HB2ST', vect, n, kd, -1, -1 )
299 lhmin = ilaenv2stage( 3, 'CHETRD_HB2ST', vect, n, kd, ib, -1 )
300 lwmin = ilaenv2stage( 4, 'CHETRD_HB2ST', vect, n, kd, ib, -1 )
301*
302 IF( .NOT.afters1 .AND. .NOT.lsame( stage1, 'N' ) ) THEN
303 info = -1
304 ELSE IF( .NOT.lsame( vect, 'N' ) ) THEN
305 info = -2
306 ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
307 info = -3
308 ELSE IF( n.LT.0 ) THEN
309 info = -4
310 ELSE IF( kd.LT.0 ) THEN
311 info = -5
312 ELSE IF( ldab.LT.(kd+1) ) THEN
313 info = -7
314 ELSE IF( lhous.LT.lhmin .AND. .NOT.lquery ) THEN
315 info = -11
316 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
317 info = -13
318 END IF
319*
320 IF( info.EQ.0 ) THEN
321 hous( 1 ) = lhmin
322 work( 1 ) = sroundup_lwork(lwmin)
323 END IF
324*
325 IF( info.NE.0 ) THEN
326 CALL xerbla( 'CHETRD_HB2ST', -info )
327 RETURN
328 ELSE IF( lquery ) THEN
329 RETURN
330 END IF
331*
332* Quick return if possible
333*
334 IF( n.EQ.0 ) THEN
335 hous( 1 ) = 1
336 work( 1 ) = 1
337 RETURN
338 END IF
339*
340* Determine pointer position
341*
342 ldv = kd + ib
343 sizetau = 2 * n
344 sicev = 2 * n
345 indtau = 1
346 indv = indtau + sizetau
347 lda = 2 * kd + 1
348 sizea = lda * n
349 inda = 1
350 indw = inda + sizea
351 nthreads = 1
352 tid = 0
353*
354 IF( upper ) THEN
355 apos = inda + kd
356 awpos = inda
357 dpos = apos + kd
358 ofdpos = dpos - 1
359 abdpos = kd + 1
360 abofdpos = kd
361 ELSE
362 apos = inda
363 awpos = inda + kd + 1
364 dpos = apos
365 ofdpos = dpos + 1
366 abdpos = 1
367 abofdpos = 2
368
369 ENDIF
370*
371* Case KD=0:
372* The matrix is diagonal. We just copy it (convert to "real" for
373* complex because D is double and the imaginary part should be 0)
374* and store it in D. A sequential code here is better or
375* in a parallel environment it might need two cores for D and E
376*
377 IF( kd.EQ.0 ) THEN
378 DO 30 i = 1, n
379 d( i ) = real( ab( abdpos, i ) )
380 30 CONTINUE
381 DO 40 i = 1, n-1
382 e( i ) = rzero
383 40 CONTINUE
384*
385 hous( 1 ) = 1
386 work( 1 ) = 1
387 RETURN
388 END IF
389*
390* Case KD=1:
391* The matrix is already Tridiagonal. We have to make diagonal
392* and offdiagonal elements real, and store them in D and E.
393* For that, for real precision just copy the diag and offdiag
394* to D and E while for the COMPLEX case the bulge chasing is
395* performed to convert the hermetian tridiagonal to symmetric
396* tridiagonal. A simpler conversion formula might be used, but then
397* updating the Q matrix will be required and based if Q is generated
398* or not this might complicate the story.
399*
400 IF( kd.EQ.1 ) THEN
401 DO 50 i = 1, n
402 d( i ) = real( ab( abdpos, i ) )
403 50 CONTINUE
404*
405* make off-diagonal elements real and copy them to E
406*
407 IF( upper ) THEN
408 DO 60 i = 1, n - 1
409 tmp = ab( abofdpos, i+1 )
410 abstmp = abs( tmp )
411 ab( abofdpos, i+1 ) = abstmp
412 e( i ) = abstmp
413 IF( abstmp.NE.rzero ) THEN
414 tmp = tmp / abstmp
415 ELSE
416 tmp = one
417 END IF
418 IF( i.LT.n-1 )
419 $ ab( abofdpos, i+2 ) = ab( abofdpos, i+2 )*tmp
420C IF( WANTZ ) THEN
421C CALL CSCAL( N, CONJG( TMP ), Q( 1, I+1 ), 1 )
422C END IF
423 60 CONTINUE
424 ELSE
425 DO 70 i = 1, n - 1
426 tmp = ab( abofdpos, i )
427 abstmp = abs( tmp )
428 ab( abofdpos, i ) = abstmp
429 e( i ) = abstmp
430 IF( abstmp.NE.rzero ) THEN
431 tmp = tmp / abstmp
432 ELSE
433 tmp = one
434 END IF
435 IF( i.LT.n-1 )
436 $ ab( abofdpos, i+1 ) = ab( abofdpos, i+1 )*tmp
437C IF( WANTQ ) THEN
438C CALL CSCAL( N, TMP, Q( 1, I+1 ), 1 )
439C END IF
440 70 CONTINUE
441 ENDIF
442*
443 hous( 1 ) = 1
444 work( 1 ) = 1
445 RETURN
446 END IF
447*
448* Main code start here.
449* Reduce the hermitian band of A to a tridiagonal matrix.
450*
451 thgrsiz = n
452 grsiz = 1
453 shift = 3
454 nbtiles = ceiling( real(n)/real(kd) )
455 stepercol = ceiling( real(shift)/real(grsiz) )
456 thgrnb = ceiling( real(n-1)/real(thgrsiz) )
457*
458 CALL clacpy( "A", kd+1, n, ab, ldab, work( apos ), lda )
459 CALL claset( "A", kd, n, zero, zero, work( awpos ), lda )
460*
461*
462* openMP parallelisation start here
463*
464#if defined(_OPENMP)
465!$OMP PARALLEL PRIVATE( TID, THGRID, BLKLASTIND )
466!$OMP$ PRIVATE( THED, I, M, K, ST, ED, STT, SWEEPID )
467!$OMP$ PRIVATE( MYID, TTYPE, COLPT, STIND, EDIND )
468!$OMP$ SHARED ( UPLO, WANTQ, INDV, INDTAU, HOUS, WORK)
469!$OMP$ SHARED ( N, KD, IB, NBTILES, LDA, LDV, INDA )
470!$OMP$ SHARED ( STEPERCOL, THGRNB, THGRSIZ, GRSIZ, SHIFT )
471!$OMP MASTER
472#endif
473*
474* main bulge chasing loop
475*
476 DO 100 thgrid = 1, thgrnb
477 stt = (thgrid-1)*thgrsiz+1
478 thed = min( (stt + thgrsiz -1), (n-1))
479 DO 110 i = stt, n-1
480 ed = min( i, thed )
481 IF( stt.GT.ed ) EXIT
482 DO 120 m = 1, stepercol
483 st = stt
484 DO 130 sweepid = st, ed
485 DO 140 k = 1, grsiz
486 myid = (i-sweepid)*(stepercol*grsiz)
487 $ + (m-1)*grsiz + k
488 IF ( myid.EQ.1 ) THEN
489 ttype = 1
490 ELSE
491 ttype = mod( myid, 2 ) + 2
492 ENDIF
493
494 IF( ttype.EQ.2 ) THEN
495 colpt = (myid/2)*kd + sweepid
496 stind = colpt-kd+1
497 edind = min(colpt,n)
498 blklastind = colpt
499 ELSE
500 colpt = ((myid+1)/2)*kd + sweepid
501 stind = colpt-kd+1
502 edind = min(colpt,n)
503 IF( ( stind.GE.edind-1 ).AND.
504 $ ( edind.EQ.n ) ) THEN
505 blklastind = n
506 ELSE
507 blklastind = 0
508 ENDIF
509 ENDIF
510*
511* Call the kernel
512*
513#if defined(_OPENMP) && _OPENMP >= 201307
514 IF( ttype.NE.1 ) THEN
515!$OMP TASK DEPEND(in:WORK(MYID+SHIFT-1))
516!$OMP$ DEPEND(in:WORK(MYID-1))
517!$OMP$ DEPEND(out:WORK(MYID))
518 tid = omp_get_thread_num()
519 CALL chb2st_kernels( uplo, wantq, ttype,
520 $ stind, edind, sweepid, n, kd, ib,
521 $ work( inda ), lda,
522 $ hous( indv ), hous( indtau ), ldv,
523 $ work( indw + tid*kd ) )
524!$OMP END TASK
525 ELSE
526!$OMP TASK DEPEND(in:WORK(MYID+SHIFT-1))
527!$OMP$ DEPEND(out:WORK(MYID))
528 tid = omp_get_thread_num()
529 CALL chb2st_kernels( uplo, wantq, ttype,
530 $ stind, edind, sweepid, n, kd, ib,
531 $ work( inda ), lda,
532 $ hous( indv ), hous( indtau ), ldv,
533 $ work( indw + tid*kd ) )
534!$OMP END TASK
535 ENDIF
536#else
537 CALL chb2st_kernels( uplo, wantq, ttype,
538 $ stind, edind, sweepid, n, kd, ib,
539 $ work( inda ), lda,
540 $ hous( indv ), hous( indtau ), ldv,
541 $ work( indw ) )
542#endif
543 IF ( blklastind.GE.(n-1) ) THEN
544 stt = stt + 1
545 EXIT
546 ENDIF
547 140 CONTINUE
548 130 CONTINUE
549 120 CONTINUE
550 110 CONTINUE
551 100 CONTINUE
552*
553#if defined(_OPENMP)
554!$OMP END MASTER
555!$OMP END PARALLEL
556#endif
557*
558* Copy the diagonal from A to D. Note that D is REAL thus only
559* the Real part is needed, the imaginary part should be zero.
560*
561 DO 150 i = 1, n
562 d( i ) = real( work( dpos+(i-1)*lda ) )
563 150 CONTINUE
564*
565* Copy the off diagonal from A to E. Note that E is REAL thus only
566* the Real part is needed, the imaginary part should be zero.
567*
568 IF( upper ) THEN
569 DO 160 i = 1, n-1
570 e( i ) = real( work( ofdpos+i*lda ) )
571 160 CONTINUE
572 ELSE
573 DO 170 i = 1, n-1
574 e( i ) = real( work( ofdpos+(i-1)*lda ) )
575 170 CONTINUE
576 ENDIF
577*
578 hous( 1 ) = lhmin
579 work( 1 ) = sroundup_lwork(lwmin)
580 RETURN
581*
582* End of CHETRD_HB2ST
583*
584 END
585
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine chb2st_kernels(uplo, wantz, ttype, st, ed, sweep, n, nb, ib, a, lda, v, tau, ldvt, work)
CHB2ST_KERNELS
subroutine chetrd_hb2st(stage1, vect, uplo, n, kd, ab, ldab, d, e, hous, lhous, work, lwork, info)
CHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric tridiagonal form T
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106