own_lapack_routines.cpp

Go to the documentation of this file.

/* ev.f -- translated by f2c (version 19940714).
   You must link the resulting object file with the libraries:
        -lf2c -lm   (in that order)
*/

#ifdef __cplusplus
extern "C" {
#endif
#include "lapack/f2c.h"
#include <stdio.h>

#ifdef _DOXYGEN
# define integer long int
# define doublereal double
# define doublecomplex __complex__ double
# define real float
# define complex __complex__ float
#endif

/* Table of constant values */

static integer c__9 = 9;
static integer c__1 = 1;
static doublereal c_b14 = 1.;
static doublecomplex c_b26 = {0.,0.};
static doublecomplex c_b27 = {1.,0.};

integer own_ew_(integer *job, integer *up, integer *n, integer *m, doublecomplex *
        ab, integer *ldab, doublecomplex *bb, integer *ldbb, doublereal *vl, 
        doublereal *vu, doublereal *w, doublecomplex *q, integer *ldq, 
        doublecomplex *work, doublereal *rwork, integer *iwork)
{
    /* System generated locals */
    integer ab_dim1, ab_offset, bb_dim1, bb_offset, q_dim1, q_offset, ret_val,
             i__1;
    doublereal d__1, d__2;

    /* Builtin functions */
    double sqrt(doublereal);
    integer s_wsle(cilist *), do_lio(integer *, integer *, const char *, ftnlen), 
            e_wsle();

    /* Local variables */
    static integer indd, info;
    static doublereal anrm;
    static char vect[1], jobz[1];
    static doublereal rmin, rmax;
    static char uplo[1];
    static integer indnd;
    static doublereal sigma;
    static char order[1];
    static logical upper, wantz;
    static integer ka, kb, il;
    extern doublereal dlamch_(const char *, ftnlen);
    static integer iu, iscale, indibl;
    static logical valeig;
    static doublereal safmin;
    extern doublereal zlanhb_(const char *, const char *, integer *, integer *, 
                doublecomplex *, integer *, doublereal *, ftnlen, ftnlen);
    extern /* Subroutine */ int xerbla_(const char *, integer *, ftnlen);
    static doublereal abstll, bignum, abstol;
    static integer indiwk, indisp;
    extern /* Subroutine */ int zlascl_(const char *, integer *, integer *, 
                doublereal *, doublereal *, integer *, integer *,
                doublecomplex *, integer *, integer *, ftnlen), 
            dstebz_(const char *, const char *, integer *, doublereal *, 
                doublereal *, integer *, integer *, doublereal *, doublereal *,
                doublereal *, integer *, integer *, doublereal *, integer *, 
                integer *, doublereal *, integer *, integer *, ftnlen, ftnlen),
            zhbtrd_(const char *, const char *, integer *, integer *, 
                doublecomplex *, integer *, doublereal *, doublereal *, 
                doublecomplex *, integer *, doublecomplex *, integer *, ftnlen,
                ftnlen);
    static integer indrwk;
    extern /* Subroutine */ int zhbgst_(const char *, const char *, integer *,
                integer *, integer *, doublecomplex *, integer *, 
                doublecomplex *, integer *, doublecomplex *, integer *, 
                doublecomplex *, doublereal *, integer *, ftnlen, ftnlen), 
           zpbstf_(const char *, integer *, integer *, 
                doublecomplex *, integer *, integer *, ftnlen);
    static integer nsplit;
    static doublereal smlnum, eps, vll, vuu;

    /* Fortran I/O blocks */
    static cilist io___22 = { 0, 6, 0, 0, 0 };



/*     .. Executable Statements .. */

/*     Test the input parameters. */

    /* Parameter adjustments */
    ab_dim1 = *ldab;
    ab_offset = ab_dim1 + 1;
    ab -= ab_offset;
    bb_dim1 = *ldbb;
    bb_offset = bb_dim1 + 1;
    bb -= bb_offset;
    --w;
    q_dim1 = *ldq;
    q_offset = q_dim1 + 1;
    q -= q_offset;
    --work;
    --rwork;
    --iwork;

    /* Function Body */
    wantz = *job == 1;
    if (wantz) {
        *jobz = 'V';
    } else {
        *jobz = 'N';
    }
    upper = *up == 1;
    if (upper) {
        *uplo = 'U';
    } else {
        *uplo = 'L';
    }
    valeig = TRUE_;
/*     Werte setzen */
    ka = *ldab - 1;
    kb = *ldbb - 1;

    info = 0;
    if (*n < 0) {
        info = -3;
    } else if (ka < 0) {
        info = -4;
    } else if (kb < 0 || kb > ka) {
        info = -5;
    } else if (*ldab < ka + 1) {
        info = -7;
    } else if (*ldbb < kb + 1) {
        info = -9;
    } else if (*n <= 0) {
        info = -110;
    }
    if (info != 0) {
        i__1 = -info;
        xerbla_("ZHBGVX", &i__1, 6L);
        ret_val = -info;
        return ret_val;
    }

/*     Get machine constants. */

    safmin = dlamch_("Safe minimum", 12L);
    eps = dlamch_("Precision", 9L);
    smlnum = safmin / eps;
    bignum = 1. / smlnum;
    rmin = sqrt(smlnum);
/* Computing MIN */
    d__1 = sqrt(bignum), d__2 = 1. / sqrt(sqrt(safmin));
    rmax = min(d__1,d__2);

/*     Fuer maximale Genauigkeit der Eigenwerte wird gesetzt: */
    abstol = safmin * 2.;

/*     Scale both matrices to allowable range, if necessary. */

    iscale = 0;
    abstll = abstol;
    vll = *vl;
    vuu = *vu;
    anrm = zlanhb_("M", uplo, n, &ka, &ab[ab_offset], ldab, &rwork[1], 1L, 1L)
            ;
    if (anrm > 0. && anrm < rmin) {
        iscale = 1;
        sigma = rmin / anrm;
    } else if (anrm > rmax) {
        iscale = 1;
        sigma = rmax / anrm;
    }
    if (iscale == 1) {
        s_wsle(&io___22);
        do_lio(&c__9, &c__1, "Matrix in bad condition, scale matrix ...", 49L);
        e_wsle();
        if (! upper) {
            zlascl_("B", &ka, &ka, &c_b14, &sigma, n, n, &ab[ab_offset], ldab,
                     &info, 1L);
            zlascl_("B", &kb, &kb, &c_b14, &sigma, n, n, &bb[bb_offset], ldbb,
                     &info, 1L);
        } else {
            zlascl_("Q", &ka, &ka, &c_b14, &sigma, n, n, &ab[ab_offset], ldab,
                     &info, 1L);
            zlascl_("Q", &kb, &kb, &c_b14, &sigma, n, n, &bb[bb_offset], ldbb,
                     &info, 1L);
        }
        if (abstol > 0.) {
            abstll = abstol * sigma;
        }
        vll = *vl * sigma;
        vuu = *vu * sigma;
    }

/*     Form a split Cholesky factorization of B. */

    zpbstf_(uplo, n, &kb, &bb[bb_offset], ldbb, &info, 1L);

        //printf("Cholesky factorization performed..\n");
    if (info != 0) {
        info = *n + info;
        return info;
    }


/*     Transform problem to standard eigenvalue problem. */

    indd = 1;
    indnd = indd + *n;
    indrwk = *n << 1;
    zhbgst_(jobz, uplo, n, &ka, &kb, &ab[ab_offset], ldab, &bb[bb_offset], 
      ldbb, &q[q_offset], ldq, &work[1], &rwork[indrwk], &info, 1L, 1L);
        //printf("Transformed to standard eigenvalueform..\n");

    if (info != 0)      return info;




/*     Reduce to tridiagonal form. */

    if (wantz) {
        *vect = 'U';
    } else {
        *vect = 'N';
    }
    zhbtrd_(vect, uplo, n, &ka, &ab[ab_offset], ldab, &rwork[indd], &rwork[
            indnd], &q[q_offset], ldq, &work[1], &info, 1L, 1L);


        //printf("Reduced to tridiagonal form..\n");

        
/*     Bei der Rueckgabe enthält nach dieser Routine RWORK(IndD) */
/*     die Diagonalelemente,RWORK(IndND) die Nebendiagonale der */
/*     Tridiagonalmatrix */

/*     Call of the Eigenvaluesolver DSTEBZ */

    if (wantz) {
        *order = 'B';
    } else {
        *order = 'E';
    }

    indibl = 1;
    indisp = indibl + *n;
    indiwk = indisp + *n;
    dstebz_("V", order, n, &vll, &vuu, &il, &iu, &abstll, &rwork[indd], &
            rwork[indnd], m, &nsplit, &w[1], &iwork[indibl], &iwork[indisp], &
            rwork[indrwk], &iwork[indiwk], &info, 1L, 1L);

        //printf("Eigenvaluesolver DSTEBZ called..\n");


/*     End of INTEGER FUNCTION ew */
    ret_val = info;
    return ret_val;
} /* ew_ */

/* ********************************************************************** */
integer own_ev_(integer *n, integer *m, doublereal *w, doublecomplex *z, integer *
        ldz, doublecomplex *q, integer *ldq, doublecomplex *work, doublereal *
        rwork, integer *iwork, integer *ifail)
{
    /* System generated locals */
    integer z_dim1, z_offset, q_dim1, q_offset, ret_val, i__1, i__2;

    /* Local variables */
    static integer indd, info, itmp1, i, j, indnd;
    extern /* Subroutine */ int zgemv_(const char *, integer *, integer *, 
                doublecomplex *, doublecomplex *, integer *, doublecomplex *,
                integer *, doublecomplex *, doublecomplex *, integer *, ftnlen),
            zcopy_(integer *, doublecomplex *, integer *, doublecomplex *,
                integer *), 
            zswap_(integer *, doublecomplex *, integer *, doublecomplex *, 
                integer *);
    static integer jj, indibl, indiwk, indisp, indrwk;
    extern /* Subroutine */ int zstein_(integer *, doublereal *, doublereal *,
                integer *, doublereal *, integer *, integer *, doublecomplex *,
                integer *, doublereal *, integer *, integer *, integer *);
    static doublereal tmp1;



/*  ===================================================================== 
*/
/*     .. Local Scalars .. */


/* ===================================================================== 
*/
/*     Executable statements */

    /* Parameter adjustments */
    --w;
    z_dim1 = *ldz;
    z_offset = z_dim1 + 1;
    z -= z_offset;
    q_dim1 = *ldq;
    q_offset = q_dim1 + 1;
    q -= q_offset;
    --work;
    --rwork;
    --iwork;
    --ifail;

    /* Function Body */
    indd = 1;
    indnd = indd + *n;
    indrwk = *n << 1;
    indibl = 1;
    indisp = indibl + *n;
    indiwk = indisp + *n;

    zstein_(n, &rwork[indd], &rwork[indnd], m, &w[1], &iwork[indibl], &iwork[
            indisp], &z[z_offset], ldz, &rwork[indrwk], &iwork[indiwk], &
            ifail[1], &info);

/*      Apply unitary matrix used in reduction to tridiagonal */
/*      form to eigenvectors returned by ZSTEIN. */

    i__1 = *m;
    for (j = 1; j <= i__1; ++j) {
        zcopy_(n, &z[j * z_dim1 + 1], &c__1, &work[1], &c__1);
        zgemv_("N", n, n, &c_b27, &q[q_offset], ldq, &work[1], &c__1, &c_b26, 
                &z[j * z_dim1 + 1], &c__1, 1L);
    }

/*     If eigenvalues are not in order, then sort them, along with */
/*     eigenvectors. */

    i__1 = *m - 1;
    for (j = 1; j <= i__1; ++j) {
        i = 0;
        tmp1 = w[j];
        i__2 = *m;
        for (jj = j + 1; jj <= i__2; ++jj) {
            if (w[jj] < tmp1) {
                i = jj;
                tmp1 = w[jj];
            }
        }

        if (i != 0) {
            itmp1 = iwork[indibl + i - 1];
            w[i] = w[j];
            iwork[indibl + i - 1] = iwork[indibl + j - 1];
            w[j] = tmp1;
            iwork[indibl + j - 1] = itmp1;
            zswap_(n, &z[i * z_dim1 + 1], &c__1, &z[j * z_dim1 + 1], &c__1);
            if (info != 0) {
                itmp1 = ifail[i];
                ifail[i] = ifail[j];
                ifail[j] = itmp1;
            }
        }
    }

/*     End of EV */

    ret_val = info;
    return ret_val;
} /* ev_ */

#ifdef __cplusplus
        }
#endif

---