src/gromacs/linearalgebra/gmx_lapack/dlascl.c

   1 #include <ctype.h>
   2 #include <math.h>
   3 #include "gromacs/utility/real.h"
   4
   5 #include "../gmx_lapack.h"
   6 #include "lapack_limits.h"
   7
   8
   9 void
  10 F77_FUNC(dlascl,DLASCL)(const char *type,
  11                         int *kl,
  12                         int *ku,
  13                         double *cfrom,
  14                         double *cto,
  15                         int *m,
  16                         int *n,
  17                         double *a,
  18                         int *lda,
  19                         int *info)
  20 {
  21   const char ch=toupper(*type);
  22   int i,j,k,l,k1,k2,k3,k4;
  23   int done=0;
  24   double minval,smlnum,bignum;
  25   double cfromc, ctoc, cfrom1, cto1, mul;
  26
  27   if(*n<=0 || *m<=0)
  28     return;
  29
  30   minval = GMX_DOUBLE_MIN;
  31   smlnum = minval / GMX_DOUBLE_EPS;
  32   bignum = 1.0 / smlnum;
  33
  34   cfromc = *cfrom;
  35   ctoc   = *cto;
  36
  37   while(!done) {
  38
  39     cfrom1 = cfromc * smlnum;
  40     cto1   = ctoc / bignum;
  41
  42     if(fabs(cfrom1)>fabs(ctoc) && fabs(ctoc)>GMX_DOUBLE_MIN) {
  43       mul = smlnum;
  44       done = 0;
  45       cfromc = cfrom1;
  46     } else if(fabs(cto1)>fabs(cfromc)) {
  47       mul = bignum;
  48       done = 0;
  49       ctoc = cto1;
  50     } else {
  51       mul = ctoc / cfromc;
  52       done = 1;
  53     }
  54
  55     switch(ch) {
  56     case 'G':
  57       /* Full matrix */
  58       for(j=0;j<*n;j++)
  59         for(i=0;i<*m;i++)
  60           a[j*(*lda)+i] *= mul;
  61       break;
  62
  63     case 'L':
  64       /* Lower triangular matrix */
  65       for(j=0;j<*n;j++)
  66         for(i=j;i<*m;i++)
  67           a[j*(*lda)+i] *= mul;
  68       break;
  69
  70     case 'U':
  71       /* Upper triangular matrix */
  72       for(j=0;j<*n;j++) {
  73         k = (j < (*m-1)) ? j : (*m-1);
  74         for(i=0;i<=k;i++)
  75           a[j*(*lda)+i] *= mul;
  76       }
  77       break;
  78
  79     case 'H':
  80       /* Upper Hessenberg matrix */
  81       for(j=0;j<*n;j++) {
  82         k = ((j+1) < (*m-1)) ? (j+1) : (*m-1);
  83         for(i=0;i<=k;i++)
  84           a[j*(*lda)+i] *= mul;
  85       }
  86       break;
  87
  88     case 'B':
  89       /* Symmetric band matrix, lower bandwidth KL, upper KU,
  90        * only the lower half stored.
  91        */
  92       k3 = *kl;
  93       k4 = *n - 1;
  94       for(j=0;j<*n;j++) {
  95         k = (k3 < (k4-j)) ? k3 : (k4-j);
  96         for(i=0;i<=k;i++)
  97           a[j*(*lda)+i] *= mul;
  98       }
  99       break;
 100
 101     case 'Q':
 102       /* Symmetric band matrix, lower bandwidth KL, upper KU,
 103        * only the upper half stored.
 104        */
 105       k1 = *ku;
 106       k3 = *ku;
 107       for(j=0;j<*n;j++) {
 108         k = ((k1-j) > 0) ? (k1-j) : 0;
 109         for(i=k;i<=k3;i++)
 110           a[j*(*lda)+i] *= mul;
 111       }
 112       break;
 113
 114     case 'Z':
 115       /* Band matrix, lower bandwidth KL, upper KU. */
 116
 117       k1 = *kl + *ku;
 118       k2 = *kl;
 119       k3 = 2*(*kl) + *ku;
 120       k4 = *kl + *ku - 1 + *m;
 121       for(j=0;j<*n;j++) {
 122         k = ((k1-j) > k2) ? (k1-j) : k2;
 123         l = (k3 < (k4-j)) ? k3 : (k4-j);
 124         for(i=k;i<=l;i++)
 125           a[j*(*lda)+i] *= mul;
 126       }
 127       break;
 128
 129     default:
 130       *info = -1;
 131       return;
 132     }
 133   } /* finished */
 134
 135   *info = 0;
 136   return;
 137 }