src/gromacs/linearalgebra/gmx_lapack/dlasr.c

   1 #include <math.h>
   2
   3 #include "types/simple.h"
   4 #include "../gmx_lapack.h"
   5
   6 void
   7 F77_FUNC(dlasr,DLASR)(const char *side,
   8        const char *pivot,
   9        const char *direct,
  10        int *m,
  11        int *n,
  12        double *c__,
  13        double *s,
  14        double *a,
  15        int *lda)
  16 {
  17     /* System generated locals */
  18     int a_dim1, a_offset, i__1, i__2;
  19
  20     /* Local variables */
  21     int i__, j, info;
  22     double temp;
  23     double ctemp, stemp;
  24
  25     --c__;
  26     --s;
  27     a_dim1 = *lda;
  28     a_offset = 1 + a_dim1;
  29     a -= a_offset;
  30
  31     /* Function Body */
  32     info = 0;
  33
  34     if (*m == 0 || *n == 0) {
  35         return;
  36     }
  37     if (*side=='L' || *side=='l') {
  38
  39         if (*pivot=='V' || *pivot=='v') {
  40             if (*direct=='F' || *direct=='f') {
  41                 i__1 = *m - 1;
  42                 for (j = 1; j <= i__1; ++j) {
  43                     ctemp = c__[j];
  44                     stemp = s[j];
  45                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
  46                         i__2 = *n;
  47                         for (i__ = 1; i__ <= i__2; ++i__) {
  48                             temp = a[j + 1 + i__ * a_dim1];
  49                             a[j + 1 + i__ * a_dim1] = ctemp * temp - stemp *
  50                                     a[j + i__ * a_dim1];
  51                             a[j + i__ * a_dim1] = stemp * temp + ctemp * a[j
  52                                     + i__ * a_dim1];
  53                         }
  54                     }
  55                 }
  56             } else if (*direct=='B' || *direct=='b') {
  57                 for (j = *m - 1; j >= 1; --j) {
  58                     ctemp = c__[j];
  59                     stemp = s[j];
  60                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
  61                         i__1 = *n;
  62                         for (i__ = 1; i__ <= i__1; ++i__) {
  63                             temp = a[j + 1 + i__ * a_dim1];
  64                             a[j + 1 + i__ * a_dim1] = ctemp * temp - stemp *
  65                                     a[j + i__ * a_dim1];
  66                             a[j + i__ * a_dim1] = stemp * temp + ctemp * a[j
  67                                     + i__ * a_dim1];
  68                         }
  69                     }
  70                 }
  71             }
  72         } else if (*pivot=='T' || *pivot=='t') {
  73             if (*direct=='F' || *direct=='f') {
  74                 i__1 = *m;
  75                 for (j = 2; j <= i__1; ++j) {
  76                     ctemp = c__[j - 1];
  77                     stemp = s[j - 1];
  78                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
  79                         i__2 = *n;
  80                         for (i__ = 1; i__ <= i__2; ++i__) {
  81                             temp = a[j + i__ * a_dim1];
  82                             a[j + i__ * a_dim1] = ctemp * temp - stemp * a[
  83                                     i__ * a_dim1 + 1];
  84                             a[i__ * a_dim1 + 1] = stemp * temp + ctemp * a[
  85                                     i__ * a_dim1 + 1];
  86                         }
  87                     }
  88                 }
  89             } else if (*direct=='B' || *direct=='b') {
  90                 for (j = *m; j >= 2; --j) {
  91                     ctemp = c__[j - 1];
  92                     stemp = s[j - 1];
  93                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
  94                         i__1 = *n;
  95                         for (i__ = 1; i__ <= i__1; ++i__) {
  96                             temp = a[j + i__ * a_dim1];
  97                             a[j + i__ * a_dim1] = ctemp * temp - stemp * a[
  98                                     i__ * a_dim1 + 1];
  99                             a[i__ * a_dim1 + 1] = stemp * temp + ctemp * a[
 100                                     i__ * a_dim1 + 1];
 101                         }
 102                     }
 103                 }
 104             }
 105         } else if (*pivot=='B' || *pivot=='b') {
 106             if (*direct=='F' || *direct=='f') {
 107                 i__1 = *m - 1;
 108                 for (j = 1; j <= i__1; ++j) {
 109                     ctemp = c__[j];
 110                     stemp = s[j];
 111                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
 112                         i__2 = *n;
 113                         for (i__ = 1; i__ <= i__2; ++i__) {
 114                             temp = a[j + i__ * a_dim1];
 115                             a[j + i__ * a_dim1] = stemp * a[*m + i__ * a_dim1]
 116                                      + ctemp * temp;
 117                             a[*m + i__ * a_dim1] = ctemp * a[*m + i__ *
 118                                     a_dim1] - stemp * temp;
 119                         }
 120                     }
 121                 }
 122             } else if (*direct=='B' || *direct=='b') {
 123                 for (j = *m - 1; j >= 1; --j) {
 124                     ctemp = c__[j];
 125                     stemp = s[j];
 126                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
 127                         i__1 = *n;
 128                         for (i__ = 1; i__ <= i__1; ++i__) {
 129                             temp = a[j + i__ * a_dim1];
 130                             a[j + i__ * a_dim1] = stemp * a[*m + i__ * a_dim1]
 131                                      + ctemp * temp;
 132                             a[*m + i__ * a_dim1] = ctemp * a[*m + i__ *
 133                                     a_dim1] - stemp * temp;
 134                         }
 135                     }
 136                 }
 137             }
 138         }
 139     } else if (*side=='R' || *side=='r') {
 140
 141         if (*pivot=='V' || *pivot=='v') {
 142             if (*direct=='F' || *direct=='f') {
 143                 i__1 = *n - 1;
 144                 for (j = 1; j <= i__1; ++j) {
 145                     ctemp = c__[j];
 146                     stemp = s[j];
 147                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
 148                         i__2 = *m;
 149                         for (i__ = 1; i__ <= i__2; ++i__) {
 150                             temp = a[i__ + (j + 1) * a_dim1];
 151                             a[i__ + (j + 1) * a_dim1] = ctemp * temp - stemp *
 152                                      a[i__ + j * a_dim1];
 153                             a[i__ + j * a_dim1] = stemp * temp + ctemp * a[
 154                                     i__ + j * a_dim1];
 155                         }
 156                     }
 157                 }
 158             } else if (*direct=='B' || *direct=='b') {
 159                 for (j = *n - 1; j >= 1; --j) {
 160                     ctemp = c__[j];
 161                     stemp = s[j];
 162                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
 163                         i__1 = *m;
 164                         for (i__ = 1; i__ <= i__1; ++i__) {
 165                             temp = a[i__ + (j + 1) * a_dim1];
 166                             a[i__ + (j + 1) * a_dim1] = ctemp * temp - stemp *
 167                                      a[i__ + j * a_dim1];
 168                             a[i__ + j * a_dim1] = stemp * temp + ctemp * a[
 169                                     i__ + j * a_dim1];
 170                         }
 171                     }
 172                 }
 173             }
 174         } else if (*pivot=='T' || *pivot=='t') {
 175             if (*direct=='F' || *direct=='f') {
 176                 i__1 = *n;
 177                 for (j = 2; j <= i__1; ++j) {
 178                     ctemp = c__[j - 1];
 179                     stemp = s[j - 1];
 180                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
 181                         i__2 = *m;
 182                         for (i__ = 1; i__ <= i__2; ++i__) {
 183                             temp = a[i__ + j * a_dim1];
 184                             a[i__ + j * a_dim1] = ctemp * temp - stemp * a[
 185                                     i__ + a_dim1];
 186                             a[i__ + a_dim1] = stemp * temp + ctemp * a[i__ +
 187                                     a_dim1];
 188                         }
 189                     }
 190                 }
 191             } else if (*direct=='B' || *direct=='b') {
 192                 for (j = *n; j >= 2; --j) {
 193                     ctemp = c__[j - 1];
 194                     stemp = s[j - 1];
 195                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
 196                         i__1 = *m;
 197                         for (i__ = 1; i__ <= i__1; ++i__) {
 198                             temp = a[i__ + j * a_dim1];
 199                             a[i__ + j * a_dim1] = ctemp * temp - stemp * a[
 200                                     i__ + a_dim1];
 201                             a[i__ + a_dim1] = stemp * temp + ctemp * a[i__ +
 202                                     a_dim1];
 203                         }
 204                     }
 205                 }
 206             }
 207         } else if (*pivot=='B' || *pivot=='b') {
 208             if (*direct=='F' || *direct=='f') {
 209                 i__1 = *n - 1;
 210                 for (j = 1; j <= i__1; ++j) {
 211                     ctemp = c__[j];
 212                     stemp = s[j];
 213                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
 214                         i__2 = *m;
 215                         for (i__ = 1; i__ <= i__2; ++i__) {
 216                             temp = a[i__ + j * a_dim1];
 217                             a[i__ + j * a_dim1] = stemp * a[i__ + *n * a_dim1]
 218                                      + ctemp * temp;
 219                             a[i__ + *n * a_dim1] = ctemp * a[i__ + *n *
 220                                     a_dim1] - stemp * temp;
 221                         }
 222                     }
 223                 }
 224             } else if (*direct=='B' || *direct=='b') {
 225                 for (j = *n - 1; j >= 1; --j) {
 226                     ctemp = c__[j];
 227                     stemp = s[j];
 228                     if (fabs(ctemp-1.0)>GMX_DOUBLE_EPS || fabs(stemp)>GMX_DOUBLE_MIN) {
 229                         i__1 = *m;
 230                         for (i__ = 1; i__ <= i__1; ++i__) {
 231                             temp = a[i__ + j * a_dim1];
 232                             a[i__ + j * a_dim1] = stemp * a[i__ + *n * a_dim1]
 233                                      + ctemp * temp;
 234                             a[i__ + *n * a_dim1] = ctemp * a[i__ + *n *
 235                                     a_dim1] - stemp * temp;
 236                         }
 237                     }
 238                 }
 239             }
 240         }
 241     }
 242
 243     return;
 244
 245 }
 246
 247