src/gmxlib/gmx_lapack/slartg.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include <math.h>
  36 #include "gmx_lapack.h"
  37 #include "lapack_limits.h"
  38
  39 #include <types/simple.h>
  40
  41 void
  42 F77_FUNC(slartg,SLARTG)(float *f,
  43         float *g,
  44         float *cs,
  45         float *sn,
  46         float *r)
  47 {
  48   float minval,safemin, safemin2, safemx2, eps;
  49   float f1,g1,f1a,g1a,scale;
  50   int i,n,count;
  51
  52   eps = GMX_FLOAT_EPS;
  53   minval = GMX_FLOAT_MIN;
  54   safemin = minval*(1.0+eps);
  55   n = 0.5*log( safemin/eps ) / log(2);
  56   safemin2 = pow(2,n);
  57
  58   safemx2 = 1.0 / safemin2;
  59
  60   if(fabs(*g)<GMX_FLOAT_MIN) {
  61     *cs = 1.0;
  62     *sn = 0.0;
  63     *r = *f;
  64   } else if (fabs(*f)<GMX_FLOAT_MIN) {
  65     *cs = 0.0;
  66     *sn = 1.0;
  67     *r = *g;
  68   } else {
  69     f1 = *f;
  70     g1 = *g;
  71     f1a = fabs(f1);
  72     g1a = fabs(g1);
  73     scale = (f1a > g1a) ? f1a : g1a;
  74     if(scale >= safemx2) {
  75       count = 0;
  76       while(scale >= safemx2) {
  77         count++;
  78         f1 *= safemin2;
  79         g1 *= safemin2;
  80         f1a = fabs(f1);
  81         g1a = fabs(g1);
  82         scale = (f1a > g1a) ? f1a : g1a;
  83       }
  84       *r = sqrt(f1*f1 + g1*g1);
  85       *cs = f1 / *r;
  86       *sn = g1 / *r;
  87       for(i=0;i<count;i++)
  88         *r *= safemx2;
  89     } else if (scale<=safemin2) {
  90       count = 0;
  91       while(scale <= safemin2) {
  92         count++;
  93         f1 *= safemx2;
  94         g1 *= safemx2;
  95         f1a = fabs(f1);
  96         g1a = fabs(g1);
  97         scale = (f1a > g1a) ? f1a : g1a;
  98       }
  99       *r = sqrt(f1*f1 + g1*g1);
 100       *cs = f1 / *r;
 101       *sn = g1 / *r;
 102       for(i=0;i<count;i++)
 103         *r *= safemin2;
 104     } else {
 105       *r = sqrt(f1*f1 + g1*g1);
 106       *cs = f1 / *r;
 107       *sn = g1 / *r;
 108     }
 109     if(fabs(*f)>fabs(*g) && *cs<0.0) {
 110       *cs *= -1.0;
 111       *sn *= -1.0;
 112       *r  *= -1.0;
 113     }
 114   }
 115   return;
 116 }
 117