util/X86MAC64/cuda/samples/6_Advanced/eigenvalues/gerschgorin.cpp

   1 /*
   2  * Copyright 1993-2013 NVIDIA Corporation.  All rights reserved.
   3  *
   4  * Please refer to the NVIDIA end user license agreement (EULA) associated
   5  * with this source code for terms and conditions that govern your use of
   6  * this software. Any use, reproduction, disclosure, or distribution of
   7  * this software and related documentation outside the terms of the EULA
   8  * is strictly prohibited.
   9  *
  10  */
  11
  12 /* Computation of Gerschgorin interval for symmetric, tridiagonal matrix */
  13
  14 #include <cstdio>
  15 #include <cstdlib>
  16 #include <cmath>
  17 #include <cfloat>
  18
  19 #include "util.h"
  20 #include "gerschgorin.h"
  21
  22 ////////////////////////////////////////////////////////////////////////////////
  23 //! Compute Gerschgorin interval for symmetric, tridiagonal matrix
  24 //! @param  d  diagonal elements
  25 //! @param  s  superdiagonal elements
  26 //! @param  n  size of matrix
  27 //! @param  lg  lower limit of Gerschgorin interval
  28 //! @param  ug  upper limit of Gerschgorin interval
  29 ////////////////////////////////////////////////////////////////////////////////
  30 void
  31 computeGerschgorin(float *d, float *s, unsigned int n, float &lg, float &ug)
  32 {
  33
  34     lg = FLT_MAX;
  35     ug = -FLT_MAX;
  36
  37     // compute bounds
  38     for (unsigned int i = 1; i < (n - 1); ++i)
  39     {
  40
  41         // sum over the absolute values of all elements of row i
  42         float sum_abs_ni = fabsf(s[i-1]) + fabsf(s[i]);
  43
  44         lg = min(lg, d[i] - sum_abs_ni);
  45         ug = max(ug, d[i] + sum_abs_ni);
  46     }
  47
  48     // first and last row, only one superdiagonal element
  49
  50     // first row
  51     lg = min(lg, d[0] - fabsf(s[0]));
  52     ug = max(ug, d[0] + fabsf(s[0]));
  53
  54     // last row
  55     lg = min(lg, d[n-1] - fabsf(s[n-2]));
  56     ug = max(ug, d[n-1] + fabsf(s[n-2]));
  57
  58     // increase interval to avoid side effects of fp arithmetic
  59     float bnorm = max(fabsf(ug), fabsf(lg));
  60
  61     // these values depend on the implmentation of floating count that is
  62     // employed in the following
  63     float psi_0 = 11 * FLT_EPSILON * bnorm;
  64     float psi_n = 11 * FLT_EPSILON * bnorm;
  65
  66     lg = lg - bnorm * 2 * n * FLT_EPSILON - psi_0;
  67     ug = ug + bnorm * 2 * n * FLT_EPSILON + psi_n;
  68
  69     ug = max(lg, ug);
  70 }
  71