LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ dsytrd_sb2st()

 subroutine dsytrd_sb2st ( character STAGE1, character VECT, character UPLO, integer N, integer KD, double precision, dimension( ldab, * ) AB, integer LDAB, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( * ) HOUS, integer LHOUS, double precision, dimension( * ) WORK, integer LWORK, integer INFO )

DSYTRD_SB2ST reduces a real symmetric band matrix A to real symmetric tridiagonal form T

Purpose:
``` DSYTRD_SB2ST reduces a real symmetric band matrix A to real symmetric
tridiagonal form T by a orthogonal similarity transformation:
Q**T * A * Q = T.```
Parameters
 [in] STAGE1 ``` STAGE1 is CHARACTER*1 = 'N': "No": to mention that the stage 1 of the reduction from dense to band using the dsytrd_sy2sb routine was not called before this routine to reproduce AB. In other term this routine is called as standalone. = 'Y': "Yes": to mention that the stage 1 of the reduction from dense to band using the dsytrd_sy2sb routine has been called to produce AB (e.g., AB is the output of dsytrd_sy2sb.``` [in] VECT ``` VECT is CHARACTER*1 = 'N': No need for the Housholder representation, and thus LHOUS is of size max(1, 4*N); = 'V': the Householder representation is needed to either generate or to apply Q later on, then LHOUS is to be queried and computed. (NOT AVAILABLE IN THIS RELEASE).``` [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] KD ``` KD is INTEGER The number of superdiagonals of the matrix A if UPLO = 'U', or the number of subdiagonals if UPLO = 'L'. KD >= 0.``` [in,out] AB ``` AB is DOUBLE PRECISION array, dimension (LDAB,N) On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). On exit, the diagonal elements of AB are overwritten by the diagonal elements of the tridiagonal matrix T; if KD > 0, the elements on the first superdiagonal (if UPLO = 'U') or the first subdiagonal (if UPLO = 'L') are overwritten by the off-diagonal elements of T; the rest of AB is overwritten by values generated during the reduction.``` [in] LDAB ``` LDAB is INTEGER The leading dimension of the array AB. LDAB >= KD+1.``` [out] D ``` D is DOUBLE PRECISION array, dimension (N) The diagonal elements of the tridiagonal matrix T.``` [out] E ``` E is DOUBLE PRECISION array, dimension (N-1) The off-diagonal elements of the tridiagonal matrix T: E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.``` [out] HOUS ``` HOUS is DOUBLE PRECISION array, dimension LHOUS, that store the Householder representation.``` [in] LHOUS ``` LHOUS is INTEGER The dimension of the array HOUS. LHOUS = MAX(1, dimension) If LWORK = -1, or LHOUS=-1, then a query is assumed; the routine only calculates the optimal size of the HOUS array, returns this value as the first entry of the HOUS array, and no error message related to LHOUS is issued by XERBLA. LHOUS = MAX(1, dimension) where dimension = 4*N if VECT='N' not available now if VECT='H' ``` [out] WORK ` WORK is DOUBLE PRECISION array, dimension LWORK.` [in] LWORK ``` LWORK is INTEGER The dimension of the array WORK. LWORK = MAX(1, dimension) If LWORK = -1, or LHOUS=-1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA. LWORK = MAX(1, dimension) where dimension = (2KD+1)*N + KD*NTHREADS where KD is the blocking size of the reduction, FACTOPTNB is the blocking used by the QR or LQ algorithm, usually FACTOPTNB=128 is a good choice NTHREADS is the number of threads used when openMP compilation is enabled, otherwise =1.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
Further Details:
```  Implemented by Azzam Haidar.

All details are available on technical report, SC11, SC13 papers.

Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
Parallel reduction to condensed forms for symmetric eigenvalue problems
using aggregated fine-grained and memory-aware kernels. In Proceedings
of 2011 International Conference for High Performance Computing,
Networking, Storage and Analysis (SC '11), New York, NY, USA,
Article 8 , 11 pages.
http://doi.acm.org/10.1145/2063384.2063394

A. Haidar, J. Kurzak, P. Luszczek, 2013.
An improved parallel singular value algorithm and its implementation
for multicore hardware, In Proceedings of 2013 International Conference
for High Performance Computing, Networking, Storage and Analysis (SC '13).
Article 90, 12 pages.
http://doi.acm.org/10.1145/2503210.2503292

A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
A novel hybrid CPU-GPU generalized eigensolver for electronic structure
calculations based on fine-grained memory aware tasks.
International Journal of High Performance Computing Applications.
Volume 28 Issue 2, Pages 196-209, May 2014.
http://hpc.sagepub.com/content/28/2/196 ```

Definition at line 228 of file dsytrd_sb2st.F.

230 *
231 #if defined(_OPENMP)
232  use omp_lib
233 #endif
234 *
235  IMPLICIT NONE
236 *
237 * -- LAPACK computational routine --
238 * -- LAPACK is a software package provided by Univ. of Tennessee, --
239 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
240 *
241 * .. Scalar Arguments ..
242  CHARACTER STAGE1, UPLO, VECT
243  INTEGER N, KD, LDAB, LHOUS, LWORK, INFO
244 * ..
245 * .. Array Arguments ..
246  DOUBLE PRECISION D( * ), E( * )
247  DOUBLE PRECISION AB( LDAB, * ), HOUS( * ), WORK( * )
248 * ..
249 *
250 * =====================================================================
251 *
252 * .. Parameters ..
253  DOUBLE PRECISION RZERO
254  DOUBLE PRECISION ZERO, ONE
255  parameter( rzero = 0.0d+0,
256  \$ zero = 0.0d+0,
257  \$ one = 1.0d+0 )
258 * ..
259 * .. Local Scalars ..
260  LOGICAL LQUERY, WANTQ, UPPER, AFTERS1
261  INTEGER I, M, K, IB, SWEEPID, MYID, SHIFT, STT, ST,
262  \$ ED, STIND, EDIND, BLKLASTIND, COLPT, THED,
263  \$ STEPERCOL, GRSIZ, THGRSIZ, THGRNB, THGRID,
264  \$ NBTILES, TTYPE, TID, NTHREADS, DEBUG,
265  \$ ABDPOS, ABOFDPOS, DPOS, OFDPOS, AWPOS,
266  \$ INDA, INDW, APOS, SIZEA, LDA, INDV, INDTAU,
267  \$ SIDEV, SIZETAU, LDV, LHMIN, LWMIN
268 * ..
269 * .. External Subroutines ..
270  EXTERNAL dsb2st_kernels, dlacpy, dlaset, xerbla
271 * ..
272 * .. Intrinsic Functions ..
273  INTRINSIC min, max, ceiling, real
274 * ..
275 * .. External Functions ..
276  LOGICAL LSAME
277  INTEGER ILAENV2STAGE
278  EXTERNAL lsame, ilaenv2stage
279 * ..
280 * .. Executable Statements ..
281 *
282 * Determine the minimal workspace size required.
283 * Test the input parameters
284 *
285  debug = 0
286  info = 0
287  afters1 = lsame( stage1, 'Y' )
288  wantq = lsame( vect, 'V' )
289  upper = lsame( uplo, 'U' )
290  lquery = ( lwork.EQ.-1 ) .OR. ( lhous.EQ.-1 )
291 *
292 * Determine the block size, the workspace size and the hous size.
293 *
294  ib = ilaenv2stage( 2, 'DSYTRD_SB2ST', vect, n, kd, -1, -1 )
295  lhmin = ilaenv2stage( 3, 'DSYTRD_SB2ST', vect, n, kd, ib, -1 )
296  lwmin = ilaenv2stage( 4, 'DSYTRD_SB2ST', vect, n, kd, ib, -1 )
297 *
298  IF( .NOT.afters1 .AND. .NOT.lsame( stage1, 'N' ) ) THEN
299  info = -1
300  ELSE IF( .NOT.lsame( vect, 'N' ) ) THEN
301  info = -2
302  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
303  info = -3
304  ELSE IF( n.LT.0 ) THEN
305  info = -4
306  ELSE IF( kd.LT.0 ) THEN
307  info = -5
308  ELSE IF( ldab.LT.(kd+1) ) THEN
309  info = -7
310  ELSE IF( lhous.LT.lhmin .AND. .NOT.lquery ) THEN
311  info = -11
312  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
313  info = -13
314  END IF
315 *
316  IF( info.EQ.0 ) THEN
317  hous( 1 ) = lhmin
318  work( 1 ) = lwmin
319  END IF
320 *
321  IF( info.NE.0 ) THEN
322  CALL xerbla( 'DSYTRD_SB2ST', -info )
323  RETURN
324  ELSE IF( lquery ) THEN
325  RETURN
326  END IF
327 *
328 * Quick return if possible
329 *
330  IF( n.EQ.0 ) THEN
331  hous( 1 ) = 1
332  work( 1 ) = 1
333  RETURN
334  END IF
335 *
336 * Determine pointer position
337 *
338  ldv = kd + ib
339  sizetau = 2 * n
340  sidev = 2 * n
341  indtau = 1
342  indv = indtau + sizetau
343  lda = 2 * kd + 1
344  sizea = lda * n
345  inda = 1
346  indw = inda + sizea
348  tid = 0
349 *
350  IF( upper ) THEN
351  apos = inda + kd
352  awpos = inda
353  dpos = apos + kd
354  ofdpos = dpos - 1
355  abdpos = kd + 1
356  abofdpos = kd
357  ELSE
358  apos = inda
359  awpos = inda + kd + 1
360  dpos = apos
361  ofdpos = dpos + 1
362  abdpos = 1
363  abofdpos = 2
364
365  ENDIF
366 *
367 * Case KD=0:
368 * The matrix is diagonal. We just copy it (convert to "real" for
369 * real because D is double and the imaginary part should be 0)
370 * and store it in D. A sequential code here is better or
371 * in a parallel environment it might need two cores for D and E
372 *
373  IF( kd.EQ.0 ) THEN
374  DO 30 i = 1, n
375  d( i ) = ( ab( abdpos, i ) )
376  30 CONTINUE
377  DO 40 i = 1, n-1
378  e( i ) = rzero
379  40 CONTINUE
380 *
381  hous( 1 ) = 1
382  work( 1 ) = 1
383  RETURN
384  END IF
385 *
386 * Case KD=1:
387 * The matrix is already Tridiagonal. We have to make diagonal
388 * and offdiagonal elements real, and store them in D and E.
389 * For that, for real precision just copy the diag and offdiag
390 * to D and E while for the COMPLEX case the bulge chasing is
391 * performed to convert the hermetian tridiagonal to symmetric
392 * tridiagonal. A simpler conversion formula might be used, but then
393 * updating the Q matrix will be required and based if Q is generated
394 * or not this might complicate the story.
395 *
396  IF( kd.EQ.1 ) THEN
397  DO 50 i = 1, n
398  d( i ) = ( ab( abdpos, i ) )
399  50 CONTINUE
400 *
401  IF( upper ) THEN
402  DO 60 i = 1, n-1
403  e( i ) = ( ab( abofdpos, i+1 ) )
404  60 CONTINUE
405  ELSE
406  DO 70 i = 1, n-1
407  e( i ) = ( ab( abofdpos, i ) )
408  70 CONTINUE
409  ENDIF
410 *
411  hous( 1 ) = 1
412  work( 1 ) = 1
413  RETURN
414  END IF
415 *
416 * Main code start here.
417 * Reduce the symmetric band of A to a tridiagonal matrix.
418 *
419  thgrsiz = n
420  grsiz = 1
421  shift = 3
422  nbtiles = ceiling( real(n)/real(kd) )
423  stepercol = ceiling( real(shift)/real(grsiz) )
424  thgrnb = ceiling( real(n-1)/real(thgrsiz) )
425 *
426  CALL dlacpy( "A", kd+1, n, ab, ldab, work( apos ), lda )
427  CALL dlaset( "A", kd, n, zero, zero, work( awpos ), lda )
428 *
429 *
430 * openMP parallelisation start here
431 *
432 #if defined(_OPENMP)
433 !\$OMP PARALLEL PRIVATE( TID, THGRID, BLKLASTIND )
434 !\$OMP\$ PRIVATE( THED, I, M, K, ST, ED, STT, SWEEPID )
435 !\$OMP\$ PRIVATE( MYID, TTYPE, COLPT, STIND, EDIND )
436 !\$OMP\$ SHARED ( UPLO, WANTQ, INDV, INDTAU, HOUS, WORK)
437 !\$OMP\$ SHARED ( N, KD, IB, NBTILES, LDA, LDV, INDA )
438 !\$OMP\$ SHARED ( STEPERCOL, THGRNB, THGRSIZ, GRSIZ, SHIFT )
439 !\$OMP MASTER
440 #endif
441 *
442 * main bulge chasing loop
443 *
444  DO 100 thgrid = 1, thgrnb
445  stt = (thgrid-1)*thgrsiz+1
446  thed = min( (stt + thgrsiz -1), (n-1))
447  DO 110 i = stt, n-1
448  ed = min( i, thed )
449  IF( stt.GT.ed ) EXIT
450  DO 120 m = 1, stepercol
451  st = stt
452  DO 130 sweepid = st, ed
453  DO 140 k = 1, grsiz
454  myid = (i-sweepid)*(stepercol*grsiz)
455  \$ + (m-1)*grsiz + k
456  IF ( myid.EQ.1 ) THEN
457  ttype = 1
458  ELSE
459  ttype = mod( myid, 2 ) + 2
460  ENDIF
461
462  IF( ttype.EQ.2 ) THEN
463  colpt = (myid/2)*kd + sweepid
464  stind = colpt-kd+1
465  edind = min(colpt,n)
466  blklastind = colpt
467  ELSE
468  colpt = ((myid+1)/2)*kd + sweepid
469  stind = colpt-kd+1
470  edind = min(colpt,n)
471  IF( ( stind.GE.edind-1 ).AND.
472  \$ ( edind.EQ.n ) ) THEN
473  blklastind = n
474  ELSE
475  blklastind = 0
476  ENDIF
477  ENDIF
478 *
479 * Call the kernel
480 *
481 #if defined(_OPENMP) && _OPENMP >= 201307
482  IF( ttype.NE.1 ) THEN
484 !\$OMP\$ DEPEND(in:WORK(MYID-1))
485 !\$OMP\$ DEPEND(out:WORK(MYID))
487  CALL dsb2st_kernels( uplo, wantq, ttype,
488  \$ stind, edind, sweepid, n, kd, ib,
489  \$ work( inda ), lda,
490  \$ hous( indv ), hous( indtau ), ldv,
491  \$ work( indw + tid*kd ) )
493  ELSE
495 !\$OMP\$ DEPEND(out:WORK(MYID))
497  CALL dsb2st_kernels( uplo, wantq, ttype,
498  \$ stind, edind, sweepid, n, kd, ib,
499  \$ work( inda ), lda,
500  \$ hous( indv ), hous( indtau ), ldv,
501  \$ work( indw + tid*kd ) )
503  ENDIF
504 #else
505  CALL dsb2st_kernels( uplo, wantq, ttype,
506  \$ stind, edind, sweepid, n, kd, ib,
507  \$ work( inda ), lda,
508  \$ hous( indv ), hous( indtau ), ldv,
509  \$ work( indw + tid*kd ) )
510 #endif
511  IF ( blklastind.GE.(n-1) ) THEN
512  stt = stt + 1
513  EXIT
514  ENDIF
515  140 CONTINUE
516  130 CONTINUE
517  120 CONTINUE
518  110 CONTINUE
519  100 CONTINUE
520 *
521 #if defined(_OPENMP)
522 !\$OMP END MASTER
523 !\$OMP END PARALLEL
524 #endif
525 *
526 * Copy the diagonal from A to D. Note that D is REAL thus only
527 * the Real part is needed, the imaginary part should be zero.
528 *
529  DO 150 i = 1, n
530  d( i ) = ( work( dpos+(i-1)*lda ) )
531  150 CONTINUE
532 *
533 * Copy the off diagonal from A to E. Note that E is REAL thus only
534 * the Real part is needed, the imaginary part should be zero.
535 *
536  IF( upper ) THEN
537  DO 160 i = 1, n-1
538  e( i ) = ( work( ofdpos+i*lda ) )
539  160 CONTINUE
540  ELSE
541  DO 170 i = 1, n-1
542  e( i ) = ( work( ofdpos+(i-1)*lda ) )
543  170 CONTINUE
544  ENDIF
545 *
546  hous( 1 ) = lhmin
547  work( 1 ) = lwmin
548  RETURN
549 *
550 * End of DSYTRD_SB2ST
551 *
subroutine dsb2st_kernels(UPLO, WANTZ, TTYPE, ST, ED, SWEEP, N, NB, IB, A, LDA, V, TAU, LDVT, WORK)
DSB2ST_KERNELS
integer function ilaenv2stage(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV2STAGE
Definition: ilaenv2stage.f:149
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
Here is the call graph for this function:
Here is the caller graph for this function: