text: /* ../netlib/dsyev.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
   on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
  #include "FLA_f2c.h" /* Table of constant values */
  static integer c__1 = 1;
  static integer c_n1 = -1;
  static integer c__0 = 0;
  static doublereal c_b17 = 1.;
  /* > \brief <b> DSYEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matr ices</b> */
  /* =========== DOCUMENTATION =========== */
  /* Online html documentation available at */
  /* http://www.netlib.org/lapack/explore-html/ */
  /* > \htmlonly */
  /* > Download DSYEV + dependencies */
  /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyev.f "> */
  /* > [TGZ]</a> */
  /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyev.f "> */
  /* > [ZIP]</a> */
  /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyev.f "> */
  /* > [TXT]</a> */
  /* > \endhtmlonly */
  /* Definition: */
  /* =========== */
  /* SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) */
  /* .. Scalar Arguments .. */
  /* CHARACTER JOBZ, UPLO */
  /* INTEGER INFO, LDA, LWORK, N */
  /* .. */
  /* .. Array Arguments .. */
  /* DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ) */
  /* .. */
  /* > \par Purpose: */
  /* ============= */
  /* > */
  /* > \verbatim */
  /* > */
  /* > DSYEV computes all eigenvalues and, optionally, eigenvectors of a */
  /* > real symmetric matrix A. */
  /* > \endverbatim */
  /* Arguments: */
  /* ========== */
  /* > \param[in] JOBZ */
  /* > \verbatim */
  /* > JOBZ is CHARACTER*1 */
  /* > = 'N': Compute eigenvalues only;
  */
  /* > = 'V': Compute eigenvalues and eigenvectors. */
  /* > \endverbatim */
  /* > */
  /* > \param[in] UPLO */
  /* > \verbatim */
  /* > UPLO is CHARACTER*1 */
  /* > = 'U': Upper triangle of A is stored;
  */
  /* > = 'L': Lower triangle of A is stored. */
  /* > \endverbatim */
  /* > */
  /* > \param[in] N */
  /* > \verbatim */
  /* > N is INTEGER */
  /* > The order of the matrix A. N >= 0. */
  /* > \endverbatim */
  /* > */
  /* > \param[in,out] A */
  /* > \verbatim */
  /* > A is DOUBLE PRECISION array, dimension (LDA, N) */
  /* > On entry, the symmetric matrix A. If UPLO = 'U', the */
  /* > leading N-by-N upper triangular part of A contains the */
  /* > upper triangular part of the matrix A. If UPLO = 'L', */
  /* > the leading N-by-N lower triangular part of A contains */
  /* > the lower triangular part of the matrix A. */
  /* > On exit, if JOBZ = 'V', then if INFO = 0, A contains the */
  /* > orthonormal eigenvectors of the matrix A. */
  /* > If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') */
  /* > or the upper triangle (if UPLO='U') of A, including the */
  /* > diagonal, is destroyed. */
  /* > \endverbatim */
  /* > */
  /* > \param[in] LDA */
  /* > \verbatim */
  /* > LDA is INTEGER */
  /* > The leading dimension of the array A. LDA >= max(1,N). */
  /* > \endverbatim */
  /* > */
  /* > \param[out] W */
  /* > \verbatim */
  /* > W is DOUBLE PRECISION array, dimension (N) */
  /* > If INFO = 0, the eigenvalues in ascending order. */
  /* > \endverbatim */
  /* > */
  /* > \param[out] WORK */
  /* > \verbatim */
  /* > WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
  /* > On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
  /* > \endverbatim */
  /* > */
  /* > \param[in] LWORK */
  /* > \verbatim */
  /* > LWORK is INTEGER */
  /* > The length of the array WORK. LWORK >= max(1,3*N-1). */
  /* > For optimal efficiency, LWORK >= (NB+2)*N, */
  /* > where NB is the blocksize for DSYTRD returned by ILAENV. */
  /* > */
  /* > If LWORK = -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. */
  /* > \endverbatim */
  /* > */
  /* > \param[out] INFO */
  /* > \verbatim */
  /* > INFO is INTEGER */
  /* > = 0: successful exit */
  /* > < 0: if INFO = -i, the i-th argument had an illegal value */
  /* > > 0: if INFO = i, the algorithm failed to converge;
  i */
  /* > off-diagonal elements of an intermediate tridiagonal */
  /* > form did not converge to zero. */
  /* > \endverbatim */
  /* Authors: */
  /* ======== */
  /* > \author Univ. of Tennessee */
  /* > \author Univ. of California Berkeley */
  /* > \author Univ. of Colorado Denver */
  /* > \author NAG Ltd. */
  /* > \date November 2011 */
  /* > \ingroup doubleSYeigen */
  /* ===================================================================== */
  /* Subroutine */
  int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a, integer *lda, doublereal *w, doublereal *work, integer *lwork, integer *info)
  {
      /* System generated locals */
      integer a_dim1, a_offset, i__1, i__2;
      doublereal d__1;
      /* Builtin functions */
      double sqrt(doublereal);
      /* Local variables */
      integer nb;
      doublereal eps;
      integer inde;
      doublereal anrm;
      integer imax;
      doublereal rmin, rmax;
      extern /* Subroutine */
      int dscal_(integer *, doublereal *, doublereal *, integer *);
      doublereal sigma;
      extern logical lsame_(char *, char *);
      integer iinfo;
      logical lower, wantz;
      extern doublereal dlamch_(char *);
      integer iscale;
      extern /* Subroutine */
      int dlascl_(char *, integer *, integer *, doublereal *, doublereal *, integer *, integer *, doublereal *, integer *, integer *);
      doublereal safmin;
      extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *);
      extern /* Subroutine */
      int xerbla_(char *, integer *);
      doublereal bignum;
      integer indtau;
      extern /* Subroutine */
      int dsterf_(integer *, doublereal *, doublereal *, integer *);
      extern doublereal dlansy_(char *, char *, integer *, doublereal *, integer *, doublereal *);
      integer indwrk;
      extern /* Subroutine */
      int dorgtr_(char *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *), dsteqr_(char *, integer *, doublereal *, doublereal *, doublereal *, integer *, doublereal *, integer *), dsytrd_(char *, integer *, doublereal *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, integer *);
      integer llwork;
      doublereal smlnum;
      integer lwkopt;
      logical lquery;
      /* -- LAPACK driver routine (version 3.4.0) -- */
      /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
      /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
      /* November 2011 */
      /* .. Scalar Arguments .. */
      /* .. */
      /* .. Array Arguments .. */
      /* .. */
      /* ===================================================================== */
      /* .. Parameters .. */
      /* .. */
      /* .. Local Scalars .. */
      /* .. */
      /* .. External Functions .. */
      /* .. */
      /* .. External Subroutines .. */
      /* .. */
      /* .. Intrinsic Functions .. */
      /* .. */
      /* .. Executable Statements .. */
      /* Test the input parameters. */
      /* Parameter adjustments */
      a_dim1 = *lda;
      a_offset = 1 + a_dim1;
      a -= a_offset;
      --w;
      --work;
      /* Function Body */
      wantz = lsame_(jobz, "V");
      lower = lsame_(uplo, "L");
      lquery = *lwork == -1;
      *info = 0;
      if (! (wantz || lsame_(jobz, "N")))
      {
          *info = -1;
      }
      else if (! (lower || lsame_(uplo, "U")))
      {
          *info = -2;
      }
      else if (*n < 0)
      {
          *info = -3;
      }
      else if (*lda < max(1,*n))
      {
          *info = -5;
      }
      if (*info == 0)
      {
          nb = ilaenv_(&c__1, "DSYTRD", uplo, n, &c_n1, &c_n1, &c_n1);
          /* Computing MAX */
          i__1 = 1;
          i__2 = (nb + 2) * *n; // , expr subst
          lwkopt = max(i__1,i__2);
          work[1] = (doublereal) lwkopt;
          /* Computing MAX */
          i__1 = 1;
          i__2 = *n * 3 - 1; // , expr subst
          if (*lwork < max(i__1,i__2) && ! lquery)
          {
              *info = -8;
          }
      }
      if (*info != 0)
      {
          i__1 = -(*info);
          xerbla_("DSYEV ", &i__1);
          return 0;
      }
      else if (lquery)
      {
          return 0;
      }
      /* Quick return if possible */
      if (*n == 0)
      {
          return 0;
      }
      if (*n == 1)
      {
          w[1] = a[a_dim1 + 1];
          work[1] = 2.;
          if (wantz)
          {
              a[a_dim1 + 1] = 1.;
          }
          return 0;
      }
      /* Get machine constants. */
      safmin = dlamch_("Safe minimum");
      eps = dlamch_("Precision");
      smlnum = safmin / eps;
      bignum = 1. / smlnum;
      rmin = sqrt(smlnum);
      rmax = sqrt(bignum);
      /* Scale matrix to allowable range, if necessary. */
      anrm = dlansy_("M", uplo, n, &a[a_offset], lda, &work[1]);
      iscale = 0;
      if (anrm > 0. && anrm < rmin)
      {
          iscale = 1;
          sigma = rmin / anrm;
      }
      else if (anrm > rmax)
      {
          iscale = 1;
          sigma = rmax / anrm;
      }
      if (iscale == 1)
      {
          dlascl_(uplo, &c__0, &c__0, &c_b17, &sigma, n, n, &a[a_offset], lda, info);
      }
      /* Call DSYTRD to reduce symmetric matrix to tridiagonal form. */
      inde = 1;
      indtau = inde + *n;
      indwrk = indtau + *n;
      llwork = *lwork - indwrk + 1;
      dsytrd_(uplo, n, &a[a_offset], lda, &w[1], &work[inde], &work[indtau], & work[indwrk], &llwork, &iinfo);
      /* For eigenvalues only, call DSTERF. For eigenvectors, first call */
      /* DORGTR to generate the orthogonal matrix, then call DSTEQR. */
      if (! wantz)
      {
          dsterf_(n, &w[1], &work[inde], info);
      }
      else
      {
          dorgtr_(uplo, n, &a[a_offset], lda, &work[indtau], &work[indwrk], & llwork, &iinfo);
          dsteqr_(jobz, n, &w[1], &work[inde], &a[a_offset], lda, &work[indtau], info);
      }
      /* If matrix was scaled, then rescale eigenvalues appropriately. */
      if (iscale == 1)
      {
          if (*info == 0)
          {
              imax = *n;
          }
          else
          {
              imax = *info - 1;
          }
          d__1 = 1. / sigma;
          dscal_(&imax, &d__1, &w[1], &c__1);
      }
      /* Set WORK(1) to optimal workspace size. */
      work[1] = (doublereal) lwkopt;
      return 0;
      /* End of DSYEV */
  }
  /* dsyev_ */
  