src/gmxlib/gmx_lapack/slae2.c

   1 #include <math.h>
   2 #include "gmx_lapack.h"
   3
   4
   5 void
   6 F77_FUNC(slae2,SLAE2)(float *a,
   7        float *b,
   8        float *c__,
   9        float *rt1,
  10        float *rt2)
  11 {
  12     float d__1;
  13     float ab, df, tb, sm, rt, adf, acmn, acmx;
  14
  15
  16     sm = *a + *c__;
  17     df = *a - *c__;
  18     adf = fabs(df);
  19     tb = *b + *b;
  20     ab = fabs(tb);
  21     if (fabs(*a) > fabs(*c__)) {
  22         acmx = *a;
  23         acmn = *c__;
  24     } else {
  25         acmx = *c__;
  26         acmn = *a;
  27     }
  28     if (adf > ab) {
  29         d__1 = ab / adf;
  30         rt = adf * sqrt(d__1 * d__1 + 1.);
  31     } else if (adf < ab) {
  32         d__1 = adf / ab;
  33         rt = ab * sqrt(d__1 * d__1 + 1.);
  34     } else {
  35
  36         rt = ab * sqrt(2.);
  37     }
  38     if (sm < 0.) {
  39         *rt1 = (sm - rt) * .5;
  40         *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
  41     } else if (sm > 0.) {
  42         *rt1 = (sm + rt) * .5;
  43         *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
  44     } else {
  45         *rt1 = rt * .5;
  46         *rt2 = rt * -.5;
  47     }
  48     return;
  49
  50 }
  51
  52