src/gromacs/linearalgebra/gmx_lapack/dlaev2.c

   1 #include <math.h>
   2 #include "gromacs/utility/real.h"
   3
   4 #include "../gmx_lapack.h"
   5
   6
   7 void
   8 F77_FUNC(dlaev2,DLAEV2)(double *   a,
   9         double *   b,
  10         double *   c__,
  11         double *   rt1,
  12         double *   rt2,
  13         double *   cs1,
  14         double *   sn1)
  15 {
  16     double d__1;
  17
  18     double ab, df, cs, ct, tb, sm, tn, rt, adf, acs;
  19     int sgn1, sgn2;
  20     double acmn, acmx;
  21
  22     sm = *a + *c__;
  23     df = *a - *c__;
  24     adf = fabs(df);
  25     tb = *b + *b;
  26     ab = fabs(tb);
  27     if (fabs(*a) > fabs(*c__)) {
  28         acmx = *a;
  29         acmn = *c__;
  30     } else {
  31         acmx = *c__;
  32         acmn = *a;
  33     }
  34     if (adf > ab) {
  35         d__1 = ab / adf;
  36         rt = adf * sqrt(d__1 * d__1 + 1.);
  37     } else if (adf < ab) {
  38         d__1 = adf / ab;
  39         rt = ab * sqrt(d__1 * d__1 + 1.);
  40     } else {
  41
  42         rt = ab * sqrt(2.);
  43     }
  44     if (sm < 0.) {
  45         *rt1 = (sm - rt) * .5;
  46         sgn1 = -1;
  47
  48         *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
  49     } else if (sm > 0.) {
  50         *rt1 = (sm + rt) * .5;
  51         sgn1 = 1;
  52         *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
  53     } else {
  54         *rt1 = rt * .5;
  55         *rt2 = rt * -.5;
  56         sgn1 = 1;
  57     }
  58     if (df >= 0.) {
  59         cs = df + rt;
  60         sgn2 = 1;
  61     } else {
  62         cs = df - rt;
  63         sgn2 = -1;
  64     }
  65     acs = fabs(cs);
  66     if (acs > ab) {
  67         ct = -tb / cs;
  68         *sn1 = 1. / sqrt(ct * ct + 1.);
  69         *cs1 = ct * *sn1;
  70     } else {
  71         if (fabs(ab)<GMX_DOUBLE_MIN) {
  72             *cs1 = 1.;
  73             *sn1 = 0.;
  74         } else {
  75             tn = -cs / tb;
  76             *cs1 = 1. / sqrt(tn * tn + 1.);
  77             *sn1 = tn * *cs1;
  78         }
  79     }
  80     if (sgn1 == sgn2) {
  81         tn = *cs1;
  82         *cs1 = -(*sn1);
  83         *sn1 = tn;
  84     }
  85     return;
  86
  87 }
  88
  89