Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /trunk/src/compiler/cl-target/fragments/eigen2x2.in
ViewVC logotype

View of /trunk/src/compiler/cl-target/fragments/eigen2x2.in

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1875 - (download) (annotate)
Wed May 16 23:07:38 2012 UTC (7 years, 2 months ago) by jhr
File size: 2907 byte(s)
  Change "eval->" to "(*eval)." to satisfy the AMD OpenCL compiler.
#define ROOT_DOUBLE 1
#define ROOT_TWO    2

/*
** Eigensolver for symmetric 2x2 matrix:
**
**  M00  M01
**  M01  M11
**
**
** Return value indicates something about the eigenvalue solution to
** the quadratic characteristic equation; see ROOT_ #defines above
** 
** HEY: the numerical precision issues here merit some more scrutiny.
*/
int Diderot_evals2x2 (
    float2 *eval,
    const Diderot_real_t _M00, const Diderot_real_t _M01,
    const Diderot_real_t _M11)
{
    Diderot_real_t mean, Q, D, M00, M01, M11;
    int roots;

  /* copy the given matrix elements */
    M00 = _M00;
    M01 = _M01;
    M11 = _M11;

  /*
  ** subtract out the eigenvalue mean (will add back to evals later);
  ** helps with numerical stability
  */
    mean = (M00 + M11)/2.0;
    M00 -= mean;
    M11 -= mean;
  
    Q = M00 - M11;
    D = 4.0*M01*M01 + Q*Q;
    if (D > EPSILON) {
      /* two distinct roots */
        Diderot_real_t vv;
        vv = sqrt(D)/2.0;
        (*eval)s0 = vv;
        (*eval)s1= -vv;
        roots = ROOT_TWO;
    }
    else {
      /* double root */
      (*eval)s0 = (*eval)s1 = 0.0;
      roots = ROOT_DOUBLE;
    }

  /* add back in the eigenvalue mean */
    (*eval)s0 += mean;
    (*eval)s1 += mean;

    return roots;
}

int Diderot_evecs2x2 (
    float2 * eval, float2 * evec,
    const Diderot_real_t _M00, const Diderot_real_t _M01,
    const Diderot_real_t _M11)
{
    Diderot_real_t mean, Q, D, M00, M01, M11;
    int roots;

  /* copy the given matrix elements */
    M00 = _M00;
    M01 = _M01;
    M11 = _M11;

  /*
  ** subtract out the eigenvalue mean (will add back to evals later);
  ** helps with numerical stability
  */
    mean = (M00 + M11)/2.0;
    M00 -= mean;
    M11 -= mean;
  
    Q = M00 - M11;
    D = 4.0*M01*M01 + Q*Q;
    if (D > EPSILON) {
      /* two distinct roots */
        Diderot_real_t vv;
        Diderot_vec2_t r1, r2;
        vv = sqrt(D)/2.0;
        (*eval)s0 = vv;
        (*eval)s1 = -vv;
      /* null space of T = M - evec[0]*I == 
         [M00 - vv      M01  ]
         [  M01      M11 - vv]
         is evec[0], but we know evec[0] and evec[1] are orthogonal,
         so row span of T is evec[1]
      */
        r1 = VEC2(M00 - vv, M01);
        r2 = VEC2(M01, M11 - vv);
        if (dot(r1,r2) > 0.0) {
            evec[1] = VEC2(r1.s0+r2.s0, r1.s1+r2.s1);
        }
        else {
            evec[1] = VEC2(r1.s0-r2.s0, r1.s1-r2.s1);
        }
        evec[1] = normalize(evec[1]);
        evec[0] = VEC2(evec[1].s1, -evec[1].s0);
        evec[0] = normalize(evec[0]);
        roots = ROOT_TWO;
    }
    else {
      /* double root */
        (*eval)s0 = (*eval)s1 = 0.0;
      /* use any basis for eigenvectors */
        evec[0] = VEC2(1.0, 0.0);
        evec[1] = VEC2(0.0, 1.0);
        roots = ROOT_DOUBLE;
    }

    /* add back in the eigenvalue mean */
    (*eval)s0 += mean;
    (*eval)s1 += mean;

    return roots;
}

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0