src/gmxlib/gmx_lapack/slascl.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 <ctype.h>
  36 #include <math.h>
  37 #include <types/simple.h>
  38
  39 #include "gmx_lapack.h"
  40 #include "lapack_limits.h"
  41
  42
  43 void
  44 F77_FUNC(slascl,SLASCL)(const char *type,
  45                         int *kl,
  46                         int *ku,
  47                         float *cfrom,
  48                         float *cto,
  49                         int *m,
  50                         int *n,
  51                         float *a,
  52                         int *lda,
  53                         int *info)
  54 {
  55   const char ch=toupper(*type);
  56   int i,j,k,l,k1,k2,k3,k4;
  57   int done=0;
  58   float minval,smlnum,bignum;
  59   float cfromc, ctoc, cfrom1, cto1, mul;
  60
  61   if(*n<=0 || *m<=0)
  62     return;
  63
  64   minval = GMX_FLOAT_MIN;
  65   smlnum = minval / GMX_FLOAT_EPS;
  66   bignum = 1.0 / smlnum;
  67
  68   cfromc = *cfrom;
  69   ctoc   = *cto;
  70
  71   while(!done) {
  72
  73     cfrom1 = cfromc * smlnum;
  74     cto1   = ctoc / bignum;
  75
  76     if(fabs(cfrom1)>fabs(ctoc) && fabs(ctoc)>GMX_FLOAT_MIN) {
  77       mul = smlnum;
  78       done = 0;
  79       cfromc = cfrom1;
  80     } else if(fabs(cto1)>fabs(cfromc)) {
  81       mul = bignum;
  82       done = 0;
  83       ctoc = cto1;
  84     } else {
  85       mul = ctoc / cfromc;
  86       done = 1;
  87     }
  88
  89     switch(ch) {
  90     case 'G':
  91       /* Full matrix */
  92       for(j=0;j<*n;j++)
  93         for(i=0;i<*m;i++)
  94           a[j*(*lda)+i] *= mul;
  95       break;
  96
  97     case 'L':
  98       /* Lower triangular matrix */
  99       for(j=0;j<*n;j++)
 100         for(i=j;i<*m;i++)
 101           a[j*(*lda)+i] *= mul;
 102       break;
 103
 104     case 'U':
 105       /* Upper triangular matrix */
 106       for(j=0;j<*n;j++) {
 107         k = (j < (*m-1)) ? j : (*m-1);
 108         for(i=0;i<=k;i++)
 109           a[j*(*lda)+i] *= mul;
 110       }
 111       break;
 112
 113     case 'H':
 114       /* Upper Hessenberg matrix */
 115       for(j=0;j<*n;j++) {
 116         k = ((j+1) < (*m-1)) ? (j+1) : (*m-1);
 117         for(i=0;i<=k;i++)
 118           a[j*(*lda)+i] *= mul;
 119       }
 120       break;
 121
 122     case 'B':
 123       /* Symmetric band matrix, lower bandwidth KL, upper KU,
 124        * only the lower half stored.
 125        */
 126       k3 = *kl;
 127       k4 = *n - 1;
 128       for(j=0;j<*n;j++) {
 129         k = (k3 < (k4-j)) ? k3 : (k4-j);
 130         for(i=0;i<=k;i++)
 131           a[j*(*lda)+i] *= mul;
 132       }
 133       break;
 134
 135     case 'Q':
 136       /* Symmetric band matrix, lower bandwidth KL, upper KU,
 137        * only the upper half stored.
 138        */
 139       k1 = *ku;
 140       k3 = *ku;
 141       for(j=0;j<*n;j++) {
 142         k = ((k1-j) > 0) ? (k1-j) : 0;
 143         for(i=k;i<=k3;i++)
 144           a[j*(*lda)+i] *= mul;
 145       }
 146       break;
 147
 148     case 'Z':
 149       /* Band matrix, lower bandwidth KL, upper KU. */
 150
 151       k1 = *kl + *ku;
 152       k2 = *kl;
 153       k3 = 2*(*kl) + *ku;
 154       k4 = *kl + *ku - 1 + *m;
 155       for(j=0;j<*n;j++) {
 156         k = ((k1-j) > k2) ? (k1-j) : k2;
 157         l = (k3 < (k4-j)) ? k3 : (k4-j);
 158         for(i=k;i<=l;i++)
 159           a[j*(*lda)+i] *= mul;
 160       }
 161       break;
 162
 163     default:
 164       *info = -1;
 165       return;
 166     }
 167   } /* finished */
 168
 169   *info = 0;
 170   return;
 171 }