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

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