Annotation of rpl/lapack/lapack/dsbev_2stage.f, revision 1.1

1.1     ! bertrand    1: *> \brief <b> DSBEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
        !             2: *
        !             3: *  @precisions fortran d -> s
        !             4: *
        !             5: *  =========== DOCUMENTATION ===========
        !             6: *
        !             7: * Online html documentation available at
        !             8: *            http://www.netlib.org/lapack/explore-html/
        !             9: *
        !            10: *> \htmlonly
        !            11: *> Download DSBEV_2STAGE + dependencies
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbev_2stage.f">
        !            13: *> [TGZ]</a>
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbev_2stage.f">
        !            15: *> [ZIP]</a>
        !            16: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbev_2stage.f">
        !            17: *> [TXT]</a>
        !            18: *> \endhtmlonly
        !            19: *
        !            20: *  Definition:
        !            21: *  ===========
        !            22: *
        !            23: *       SUBROUTINE DSBEV_2STAGE( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
        !            24: *                                WORK, LWORK, INFO )
        !            25: *
        !            26: *       IMPLICIT NONE
        !            27: *
        !            28: *       .. Scalar Arguments ..
        !            29: *       CHARACTER          JOBZ, UPLO
        !            30: *       INTEGER            INFO, KD, LDAB, LDZ, N, LWORK
        !            31: *       ..
        !            32: *       .. Array Arguments ..
        !            33: *       DOUBLE PRECISION   AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
        !            34: *       ..
        !            35: *
        !            36: *
        !            37: *> \par Purpose:
        !            38: *  =============
        !            39: *>
        !            40: *> \verbatim
        !            41: *>
        !            42: *> DSBEV_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
        !            43: *> a real symmetric band matrix A using the 2stage technique for
        !            44: *> the reduction to tridiagonal.
        !            45: *> \endverbatim
        !            46: *
        !            47: *  Arguments:
        !            48: *  ==========
        !            49: *
        !            50: *> \param[in] JOBZ
        !            51: *> \verbatim
        !            52: *>          JOBZ is CHARACTER*1
        !            53: *>          = 'N':  Compute eigenvalues only;
        !            54: *>          = 'V':  Compute eigenvalues and eigenvectors.
        !            55: *>                  Not available in this release.
        !            56: *> \endverbatim
        !            57: *>
        !            58: *> \param[in] UPLO
        !            59: *> \verbatim
        !            60: *>          UPLO is CHARACTER*1
        !            61: *>          = 'U':  Upper triangle of A is stored;
        !            62: *>          = 'L':  Lower triangle of A is stored.
        !            63: *> \endverbatim
        !            64: *>
        !            65: *> \param[in] N
        !            66: *> \verbatim
        !            67: *>          N is INTEGER
        !            68: *>          The order of the matrix A.  N >= 0.
        !            69: *> \endverbatim
        !            70: *>
        !            71: *> \param[in] KD
        !            72: *> \verbatim
        !            73: *>          KD is INTEGER
        !            74: *>          The number of superdiagonals of the matrix A if UPLO = 'U',
        !            75: *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
        !            76: *> \endverbatim
        !            77: *>
        !            78: *> \param[in,out] AB
        !            79: *> \verbatim
        !            80: *>          AB is DOUBLE PRECISION array, dimension (LDAB, N)
        !            81: *>          On entry, the upper or lower triangle of the symmetric band
        !            82: *>          matrix A, stored in the first KD+1 rows of the array.  The
        !            83: *>          j-th column of A is stored in the j-th column of the array AB
        !            84: *>          as follows:
        !            85: *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
        !            86: *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
        !            87: *>
        !            88: *>          On exit, AB is overwritten by values generated during the
        !            89: *>          reduction to tridiagonal form.  If UPLO = 'U', the first
        !            90: *>          superdiagonal and the diagonal of the tridiagonal matrix T
        !            91: *>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
        !            92: *>          the diagonal and first subdiagonal of T are returned in the
        !            93: *>          first two rows of AB.
        !            94: *> \endverbatim
        !            95: *>
        !            96: *> \param[in] LDAB
        !            97: *> \verbatim
        !            98: *>          LDAB is INTEGER
        !            99: *>          The leading dimension of the array AB.  LDAB >= KD + 1.
        !           100: *> \endverbatim
        !           101: *>
        !           102: *> \param[out] W
        !           103: *> \verbatim
        !           104: *>          W is DOUBLE PRECISION array, dimension (N)
        !           105: *>          If INFO = 0, the eigenvalues in ascending order.
        !           106: *> \endverbatim
        !           107: *>
        !           108: *> \param[out] Z
        !           109: *> \verbatim
        !           110: *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
        !           111: *>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
        !           112: *>          eigenvectors of the matrix A, with the i-th column of Z
        !           113: *>          holding the eigenvector associated with W(i).
        !           114: *>          If JOBZ = 'N', then Z is not referenced.
        !           115: *> \endverbatim
        !           116: *>
        !           117: *> \param[in] LDZ
        !           118: *> \verbatim
        !           119: *>          LDZ is INTEGER
        !           120: *>          The leading dimension of the array Z.  LDZ >= 1, and if
        !           121: *>          JOBZ = 'V', LDZ >= max(1,N).
        !           122: *> \endverbatim
        !           123: *>
        !           124: *> \param[out] WORK
        !           125: *> \verbatim
        !           126: *>          WORK is DOUBLE PRECISION array, dimension LWORK
        !           127: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
        !           128: *> \endverbatim
        !           129: *>
        !           130: *> \param[in] LWORK
        !           131: *> \verbatim
        !           132: *>          LWORK is INTEGER
        !           133: *>          The length of the array WORK. LWORK >= 1, when N <= 1;
        !           134: *>          otherwise  
        !           135: *>          If JOBZ = 'N' and N > 1, LWORK must be queried.
        !           136: *>                                   LWORK = MAX(1, dimension) where
        !           137: *>                                   dimension = (2KD+1)*N + KD*NTHREADS + N
        !           138: *>                                   where KD is the size of the band.
        !           139: *>                                   NTHREADS is the number of threads used when
        !           140: *>                                   openMP compilation is enabled, otherwise =1.
        !           141: *>          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available.
        !           142: *>
        !           143: *>          If LWORK = -1, then a workspace query is assumed; the routine
        !           144: *>          only calculates the optimal size of the WORK array, returns
        !           145: *>          this value as the first entry of the WORK array, and no error
        !           146: *>          message related to LWORK is issued by XERBLA.
        !           147: *> \endverbatim
        !           148: *>
        !           149: *> \param[out] INFO
        !           150: *> \verbatim
        !           151: *>          INFO is INTEGER
        !           152: *>          = 0:  successful exit
        !           153: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
        !           154: *>          > 0:  if INFO = i, the algorithm failed to converge; i
        !           155: *>                off-diagonal elements of an intermediate tridiagonal
        !           156: *>                form did not converge to zero.
        !           157: *> \endverbatim
        !           158: *
        !           159: *  Authors:
        !           160: *  ========
        !           161: *
        !           162: *> \author Univ. of Tennessee
        !           163: *> \author Univ. of California Berkeley
        !           164: *> \author Univ. of Colorado Denver
        !           165: *> \author NAG Ltd.
        !           166: *
        !           167: *> \date December 2016
        !           168: *
        !           169: *> \ingroup doubleOTHEReigen
        !           170: *
        !           171: *> \par Further Details:
        !           172: *  =====================
        !           173: *>
        !           174: *> \verbatim
        !           175: *>
        !           176: *>  All details about the 2stage techniques are available in:
        !           177: *>
        !           178: *>  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
        !           179: *>  Parallel reduction to condensed forms for symmetric eigenvalue problems
        !           180: *>  using aggregated fine-grained and memory-aware kernels. In Proceedings
        !           181: *>  of 2011 International Conference for High Performance Computing,
        !           182: *>  Networking, Storage and Analysis (SC '11), New York, NY, USA,
        !           183: *>  Article 8 , 11 pages.
        !           184: *>  http://doi.acm.org/10.1145/2063384.2063394
        !           185: *>
        !           186: *>  A. Haidar, J. Kurzak, P. Luszczek, 2013.
        !           187: *>  An improved parallel singular value algorithm and its implementation 
        !           188: *>  for multicore hardware, In Proceedings of 2013 International Conference
        !           189: *>  for High Performance Computing, Networking, Storage and Analysis (SC '13).
        !           190: *>  Denver, Colorado, USA, 2013.
        !           191: *>  Article 90, 12 pages.
        !           192: *>  http://doi.acm.org/10.1145/2503210.2503292
        !           193: *>
        !           194: *>  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
        !           195: *>  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 
        !           196: *>  calculations based on fine-grained memory aware tasks.
        !           197: *>  International Journal of High Performance Computing Applications.
        !           198: *>  Volume 28 Issue 2, Pages 196-209, May 2014.
        !           199: *>  http://hpc.sagepub.com/content/28/2/196 
        !           200: *>
        !           201: *> \endverbatim
        !           202: *
        !           203: *  =====================================================================
        !           204:       SUBROUTINE DSBEV_2STAGE( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
        !           205:      $                         WORK, LWORK, INFO )
        !           206: *
        !           207:       IMPLICIT NONE
        !           208: *
        !           209: *  -- LAPACK driver routine (version 3.7.0) --
        !           210: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
        !           211: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
        !           212: *     December 2016
        !           213: *
        !           214: *     .. Scalar Arguments ..
        !           215:       CHARACTER          JOBZ, UPLO
        !           216:       INTEGER            INFO, KD, LDAB, LDZ, N, LWORK
        !           217: *     ..
        !           218: *     .. Array Arguments ..
        !           219:       DOUBLE PRECISION   AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
        !           220: *     ..
        !           221: *
        !           222: *  =====================================================================
        !           223: *
        !           224: *     .. Parameters ..
        !           225:       DOUBLE PRECISION   ZERO, ONE
        !           226:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
        !           227: *     ..
        !           228: *     .. Local Scalars ..
        !           229:       LOGICAL            LOWER, WANTZ, LQUERY
        !           230:       INTEGER            IINFO, IMAX, INDE, INDWRK, ISCALE,
        !           231:      $                   LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS
        !           232:       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
        !           233:      $                   SMLNUM
        !           234: *     ..
        !           235: *     .. External Functions ..
        !           236:       LOGICAL            LSAME
        !           237:       INTEGER            ILAENV
        !           238:       DOUBLE PRECISION   DLAMCH, DLANSB
        !           239:       EXTERNAL           LSAME, DLAMCH, DLANSB, ILAENV
        !           240: *     ..
        !           241: *     .. External Subroutines ..
        !           242:       EXTERNAL           DLASCL, DSCAL, DSTEQR, DSTERF, XERBLA
        !           243:      $                   DSYTRD_SB2ST 
        !           244: *     ..
        !           245: *     .. Intrinsic Functions ..
        !           246:       INTRINSIC          SQRT
        !           247: *     ..
        !           248: *     .. Executable Statements ..
        !           249: *
        !           250: *     Test the input parameters.
        !           251: *
        !           252:       WANTZ = LSAME( JOBZ, 'V' )
        !           253:       LOWER = LSAME( UPLO, 'L' )
        !           254:       LQUERY = ( LWORK.EQ.-1 )
        !           255: *
        !           256:       INFO = 0
        !           257:       IF( .NOT.( LSAME( JOBZ, 'N' ) ) ) THEN
        !           258:          INFO = -1
        !           259:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
        !           260:          INFO = -2
        !           261:       ELSE IF( N.LT.0 ) THEN
        !           262:          INFO = -3
        !           263:       ELSE IF( KD.LT.0 ) THEN
        !           264:          INFO = -4
        !           265:       ELSE IF( LDAB.LT.KD+1 ) THEN
        !           266:          INFO = -6
        !           267:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
        !           268:          INFO = -9
        !           269:       END IF
        !           270: *
        !           271:       IF( INFO.EQ.0 ) THEN
        !           272:          IF( N.LE.1 ) THEN
        !           273:             LWMIN = 1
        !           274:             WORK( 1 ) = LWMIN
        !           275:          ELSE
        !           276:             IB    = ILAENV( 18, 'DSYTRD_SB2ST', JOBZ, N, KD, -1, -1 )
        !           277:             LHTRD = ILAENV( 19, 'DSYTRD_SB2ST', JOBZ, N, KD, IB, -1 )
        !           278:             LWTRD = ILAENV( 20, 'DSYTRD_SB2ST', JOBZ, N, KD, IB, -1 )
        !           279:             LWMIN = N + LHTRD + LWTRD
        !           280:             WORK( 1 )  = LWMIN
        !           281:          ENDIF
        !           282: *
        !           283:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY )
        !           284:      $      INFO = -11
        !           285:       END IF
        !           286: *
        !           287:       IF( INFO.NE.0 ) THEN
        !           288:          CALL XERBLA( 'DSBEV_2STAGE ', -INFO )
        !           289:          RETURN
        !           290:       ELSE IF( LQUERY ) THEN
        !           291:          RETURN
        !           292:       END IF
        !           293: *
        !           294: *     Quick return if possible
        !           295: *
        !           296:       IF( N.EQ.0 )
        !           297:      $   RETURN
        !           298: *
        !           299:       IF( N.EQ.1 ) THEN
        !           300:          IF( LOWER ) THEN
        !           301:             W( 1 ) = AB( 1, 1 )
        !           302:          ELSE
        !           303:             W( 1 ) = AB( KD+1, 1 )
        !           304:          END IF
        !           305:          IF( WANTZ )
        !           306:      $      Z( 1, 1 ) = ONE
        !           307:          RETURN
        !           308:       END IF
        !           309: *
        !           310: *     Get machine constants.
        !           311: *
        !           312:       SAFMIN = DLAMCH( 'Safe minimum' )
        !           313:       EPS    = DLAMCH( 'Precision' )
        !           314:       SMLNUM = SAFMIN / EPS
        !           315:       BIGNUM = ONE / SMLNUM
        !           316:       RMIN   = SQRT( SMLNUM )
        !           317:       RMAX   = SQRT( BIGNUM )
        !           318: *
        !           319: *     Scale matrix to allowable range, if necessary.
        !           320: *
        !           321:       ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
        !           322:       ISCALE = 0
        !           323:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
        !           324:          ISCALE = 1
        !           325:          SIGMA = RMIN / ANRM
        !           326:       ELSE IF( ANRM.GT.RMAX ) THEN
        !           327:          ISCALE = 1
        !           328:          SIGMA = RMAX / ANRM
        !           329:       END IF
        !           330:       IF( ISCALE.EQ.1 ) THEN
        !           331:          IF( LOWER ) THEN
        !           332:             CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
        !           333:          ELSE
        !           334:             CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
        !           335:          END IF
        !           336:       END IF
        !           337: *
        !           338: *     Call DSYTRD_SB2ST to reduce symmetric band matrix to tridiagonal form.
        !           339: *
        !           340:       INDE    = 1
        !           341:       INDHOUS = INDE + N
        !           342:       INDWRK  = INDHOUS + LHTRD
        !           343:       LLWORK  = LWORK - INDWRK + 1
        !           344: *
        !           345:       CALL DSYTRD_SB2ST( "N", JOBZ, UPLO, N, KD, AB, LDAB, W,
        !           346:      $                    WORK( INDE ), WORK( INDHOUS ), LHTRD, 
        !           347:      $                    WORK( INDWRK ), LLWORK, IINFO )
        !           348: *
        !           349: *     For eigenvalues only, call DSTERF.  For eigenvectors, call SSTEQR.
        !           350: *
        !           351:       IF( .NOT.WANTZ ) THEN
        !           352:          CALL DSTERF( N, W, WORK( INDE ), INFO )
        !           353:       ELSE
        !           354:          CALL DSTEQR( JOBZ, N, W, WORK( INDE ), Z, LDZ, WORK( INDWRK ),
        !           355:      $                INFO )
        !           356:       END IF
        !           357: *
        !           358: *     If matrix was scaled, then rescale eigenvalues appropriately.
        !           359: *
        !           360:       IF( ISCALE.EQ.1 ) THEN
        !           361:          IF( INFO.EQ.0 ) THEN
        !           362:             IMAX = N
        !           363:          ELSE
        !           364:             IMAX = INFO - 1
        !           365:          END IF
        !           366:          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
        !           367:       END IF
        !           368: *
        !           369: *     Set WORK(1) to optimal workspace size.
        !           370: *
        !           371:       WORK( 1 ) = LWMIN
        !           372: *
        !           373:       RETURN
        !           374: *
        !           375: *     End of DSBEV_2STAGE
        !           376: *
        !           377:       END

CVSweb interface <joel.bertrand@systella.fr>