src/gmxlib/gmx_lapack/dlasq1.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, 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_blas.h"
  37 #include "gmx_lapack.h"
  38 #include "lapack_limits.h"
  39
  40 #include <types/simple.h>
  41
  42 void
  43 F77_FUNC(dlasq1,DLASQ1)(int *n,
  44         double *d,
  45         double *e,
  46         double *work,
  47         int *info)
  48 {
  49   double sigmx = 0.0;
  50   int i,j,k,iinfo;
  51   double minval,safemin;
  52   double dtemp,scale;
  53   double eps;
  54
  55   eps = GMX_DOUBLE_EPS;
  56   minval = GMX_DOUBLE_MIN;
  57   safemin = minval*(1.0+GMX_DOUBLE_EPS);
  58   *info = 0;
  59
  60   if(*n<0) {
  61     *info = -2;
  62     return;
  63   }
  64
  65   for(i=0;i<*n-1;i++) {
  66     d[i] = fabs(d[i]);
  67     dtemp = fabs(e[i]);
  68     if(dtemp>sigmx)
  69       sigmx=dtemp;
  70   }
  71   d[*n-1] = fabs(d[*n-1]);
  72
  73   if(fabs(sigmx)<GMX_DOUBLE_MIN) {
  74     F77_FUNC(dlasrt,DLASRT)("D",n,d,&iinfo);
  75     return;
  76   }
  77
  78   for(i=0;i<*n;i++) {
  79     if(d[i]>sigmx)
  80       sigmx=d[i];
  81   }
  82
  83   /* Copy d and e into work (z format) and scale.
  84    * Squaring input data makes scaling by a power of the
  85    * radix pointless.
  86    */
  87   scale = sqrt(eps/safemin);
  88   i = 1;
  89   j = 2;
  90   F77_FUNC(dcopy,DCOPY)(n,d,&i,work,&j);
  91   k = *n-1;
  92   F77_FUNC(dcopy,DCOPY)(&k,e,&i,work+1,&j);
  93   i = 0;
  94   j = 2*(*n)-1;
  95   k = 1;
  96   F77_FUNC(dlascl,DLASCL)("G",&i,&i,&sigmx,&scale,&j,&k,work,&j,&iinfo);
  97
  98
  99   /* Compute q and e elements */
 100   for(i=0;i<2*(*n)-1;i++)
 101     work[i] = work[i]*work[i];
 102
 103   work[2*(*n)-1] = 0.0;
 104
 105   F77_FUNC(dlasq2,DLASQ2)(n,work,info);
 106
 107   j = 0;
 108   k = 1;
 109   if(*info==0) {
 110     for(i=0;i<*n;i++)
 111       d[i]=sqrt(work[i]);
 112     F77_FUNC(dlascl,DLASCL)("G",&j,&j,&scale,&sigmx,n,&k,d,n,&iinfo);
 113   }
 114   return;
 115 }