git.maemo.org Git - opencv/blob - 3rdparty/lapack/dlarra.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlarra_(integer *n, doublereal *d__, doublereal *e,
   4         doublereal *e2, doublereal *spltol, doublereal *tnrm, integer *nsplit,
   5          integer *isplit, integer *info)
   6 {
   7     /* System generated locals */
   8     integer i__1;
   9     doublereal d__1, d__2;
  10
  11     /* Builtin functions */
  12     double sqrt(doublereal);
  13
  14     /* Local variables */
  15     integer i__;
  16     doublereal tmp1, eabs;
  17
  18
  19 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  20 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  21 /*     November 2006 */
  22
  23 /*     .. Scalar Arguments .. */
  24 /*     .. */
  25 /*     .. Array Arguments .. */
  26 /*     .. */
  27
  28 /*  Purpose */
  29 /*  ======= */
  30
  31 /*  Compute the splitting points with threshold SPLTOL. */
  32 /*  DLARRA sets any "small" off-diagonal elements to zero. */
  33
  34 /*  Arguments */
  35 /*  ========= */
  36
  37 /*  N       (input) INTEGER */
  38 /*          The order of the matrix. N > 0. */
  39
  40 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
  41 /*          On entry, the N diagonal elements of the tridiagonal */
  42 /*          matrix T. */
  43
  44 /*  E       (input/output) DOUBLE PRECISION array, dimension (N) */
  45 /*          On entry, the first (N-1) entries contain the subdiagonal */
  46 /*          elements of the tridiagonal matrix T; E(N) need not be set. */
  47 /*          On exit, the entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT, */
  48 /*          are set to zero, the other entries of E are untouched. */
  49
  50 /*  E2      (input/output) DOUBLE PRECISION array, dimension (N) */
  51 /*          On entry, the first (N-1) entries contain the SQUARES of the */
  52 /*          subdiagonal elements of the tridiagonal matrix T; */
  53 /*          E2(N) need not be set. */
  54 /*          On exit, the entries E2( ISPLIT( I ) ), */
  55 /*          1 <= I <= NSPLIT, have been set to zero */
  56
  57 /*  SPLTOL (input) DOUBLE PRECISION */
  58 /*          The threshold for splitting. Two criteria can be used: */
  59 /*          SPLTOL<0 : criterion based on absolute off-diagonal value */
  60 /*          SPLTOL>0 : criterion that preserves relative accuracy */
  61
  62 /*  TNRM (input) DOUBLE PRECISION */
  63 /*          The norm of the matrix. */
  64
  65 /*  NSPLIT  (output) INTEGER */
  66 /*          The number of blocks T splits into. 1 <= NSPLIT <= N. */
  67
  68 /*  ISPLIT  (output) INTEGER array, dimension (N) */
  69 /*          The splitting points, at which T breaks up into blocks. */
  70 /*          The first block consists of rows/columns 1 to ISPLIT(1), */
  71 /*          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
  72 /*          etc., and the NSPLIT-th consists of rows/columns */
  73 /*          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
  74
  75
  76 /*  INFO    (output) INTEGER */
  77 /*          = 0:  successful exit */
  78
  79 /*  Further Details */
  80 /*  =============== */
  81
  82 /*  Based on contributions by */
  83 /*     Beresford Parlett, University of California, Berkeley, USA */
  84 /*     Jim Demmel, University of California, Berkeley, USA */
  85 /*     Inderjit Dhillon, University of Texas, Austin, USA */
  86 /*     Osni Marques, LBNL/NERSC, USA */
  87 /*     Christof Voemel, University of California, Berkeley, USA */
  88
  89 /*  ===================================================================== */
  90
  91 /*     .. Parameters .. */
  92 /*     .. */
  93 /*     .. Local Scalars .. */
  94 /*     .. */
  95 /*     .. Intrinsic Functions .. */
  96 /*     .. */
  97 /*     .. Executable Statements .. */
  98
  99     /* Parameter adjustments */
 100     --isplit;
 101     --e2;
 102     --e;
 103     --d__;
 104
 105     /* Function Body */
 106     *info = 0;
 107 /*     Compute splitting points */
 108     *nsplit = 1;
 109     if (*spltol < 0.) {
 110 /*        Criterion based on absolute off-diagonal value */
 111         tmp1 = abs(*spltol) * *tnrm;
 112         i__1 = *n - 1;
 113         for (i__ = 1; i__ <= i__1; ++i__) {
 114             eabs = (d__1 = e[i__], abs(d__1));
 115             if (eabs <= tmp1) {
 116                 e[i__] = 0.;
 117                 e2[i__] = 0.;
 118                 isplit[*nsplit] = i__;
 119                 ++(*nsplit);
 120             }
 121 /* L9: */
 122         }
 123     } else {
 124 /*        Criterion that guarantees relative accuracy */
 125         i__1 = *n - 1;
 126         for (i__ = 1; i__ <= i__1; ++i__) {
 127             eabs = (d__1 = e[i__], abs(d__1));
 128             if (eabs <= *spltol * sqrt((d__1 = d__[i__], abs(d__1))) * sqrt((
 129                     d__2 = d__[i__ + 1], abs(d__2)))) {
 130                 e[i__] = 0.;
 131                 e2[i__] = 0.;
 132                 isplit[*nsplit] = i__;
 133                 ++(*nsplit);
 134             }
 135 /* L10: */
 136         }
 137     }
 138     isplit[*nsplit] = *n;
 139     return 0;
 140
 141 /*     End of DLARRA */
 142
 143 } /* dlarra_ */