src/gmxlib/gmx_lapack/slasr.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
  37 #include <types/simple.h>
  38 #include "gmx_lapack.h"
  39
  40 void
  41 F77_FUNC(slasr,SLASR)(const char *side,
  42        const char *pivot,
  43        const char *direct,
  44        int *m,
  45        int *n,
  46        float *c__,
  47        float *s,
  48        float *a,
  49        int *lda)
  50 {
  51     /* System generated locals */
  52     int a_dim1, a_offset, i__1, i__2;
  53
  54     /* Local variables */
  55     int i__, j, info;
  56     float temp;
  57     float ctemp, stemp;
  58
  59     --c__;
  60     --s;
  61     a_dim1 = *lda;
  62     a_offset = 1 + a_dim1;
  63     a -= a_offset;
  64
  65     /* Function Body */
  66     info = 0;
  67
  68     if (*m == 0 || *n == 0) {
  69         return;
  70     }
  71     if (*side=='L' || *side=='l') {
  72
  73         if (*pivot=='V' || *pivot=='v') {
  74             if (*direct=='F' || *direct=='f') {
  75                 i__1 = *m - 1;
  76                 for (j = 1; j <= i__1; ++j) {
  77                     ctemp = c__[j];
  78                     stemp = s[j];
  79                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
  80                         i__2 = *n;
  81                         for (i__ = 1; i__ <= i__2; ++i__) {
  82                             temp = a[j + 1 + i__ * a_dim1];
  83                             a[j + 1 + i__ * a_dim1] = ctemp * temp - stemp *
  84                                     a[j + i__ * a_dim1];
  85                             a[j + i__ * a_dim1] = stemp * temp + ctemp * a[j
  86                                     + i__ * a_dim1];
  87                         }
  88                     }
  89                 }
  90             } else if (*direct=='B' || *direct=='b') {
  91                 for (j = *m - 1; j >= 1; --j) {
  92                     ctemp = c__[j];
  93                     stemp = s[j];
  94                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
  95                         i__1 = *n;
  96                         for (i__ = 1; i__ <= i__1; ++i__) {
  97                             temp = a[j + 1 + i__ * a_dim1];
  98                             a[j + 1 + i__ * a_dim1] = ctemp * temp - stemp *
  99                                     a[j + i__ * a_dim1];
 100                             a[j + i__ * a_dim1] = stemp * temp + ctemp * a[j
 101                                     + i__ * a_dim1];
 102                         }
 103                     }
 104                 }
 105             }
 106         } else if (*pivot=='T' || *pivot=='t') {
 107             if (*direct=='F' || *direct=='f') {
 108                 i__1 = *m;
 109                 for (j = 2; j <= i__1; ++j) {
 110                     ctemp = c__[j - 1];
 111                     stemp = s[j - 1];
 112                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
 113                         i__2 = *n;
 114                         for (i__ = 1; i__ <= i__2; ++i__) {
 115                             temp = a[j + i__ * a_dim1];
 116                             a[j + i__ * a_dim1] = ctemp * temp - stemp * a[
 117                                     i__ * a_dim1 + 1];
 118                             a[i__ * a_dim1 + 1] = stemp * temp + ctemp * a[
 119                                     i__ * a_dim1 + 1];
 120                         }
 121                     }
 122                 }
 123             } else if (*direct=='B' || *direct=='b') {
 124                 for (j = *m; j >= 2; --j) {
 125                     ctemp = c__[j - 1];
 126                     stemp = s[j - 1];
 127                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
 128                         i__1 = *n;
 129                         for (i__ = 1; i__ <= i__1; ++i__) {
 130                             temp = a[j + i__ * a_dim1];
 131                             a[j + i__ * a_dim1] = ctemp * temp - stemp * a[
 132                                     i__ * a_dim1 + 1];
 133                             a[i__ * a_dim1 + 1] = stemp * temp + ctemp * a[
 134                                     i__ * a_dim1 + 1];
 135                         }
 136                     }
 137                 }
 138             }
 139         } else if (*pivot=='B' || *pivot=='b') {
 140             if (*direct=='F' || *direct=='f') {
 141                 i__1 = *m - 1;
 142                 for (j = 1; j <= i__1; ++j) {
 143                     ctemp = c__[j];
 144                     stemp = s[j];
 145                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
 146                         i__2 = *n;
 147                         for (i__ = 1; i__ <= i__2; ++i__) {
 148                             temp = a[j + i__ * a_dim1];
 149                             a[j + i__ * a_dim1] = stemp * a[*m + i__ * a_dim1]
 150                                      + ctemp * temp;
 151                             a[*m + i__ * a_dim1] = ctemp * a[*m + i__ *
 152                                     a_dim1] - stemp * temp;
 153                         }
 154                     }
 155                 }
 156             } else if (*direct=='B' || *direct=='b') {
 157                 for (j = *m - 1; j >= 1; --j) {
 158                     ctemp = c__[j];
 159                     stemp = s[j];
 160                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
 161                         i__1 = *n;
 162                         for (i__ = 1; i__ <= i__1; ++i__) {
 163                             temp = a[j + i__ * a_dim1];
 164                             a[j + i__ * a_dim1] = stemp * a[*m + i__ * a_dim1]
 165                                      + ctemp * temp;
 166                             a[*m + i__ * a_dim1] = ctemp * a[*m + i__ *
 167                                     a_dim1] - stemp * temp;
 168                         }
 169                     }
 170                 }
 171             }
 172         }
 173     } else if (*side=='R' || *side=='r') {
 174
 175         if (*pivot=='V' || *pivot=='v') {
 176             if (*direct=='F' || *direct=='f') {
 177                 i__1 = *n - 1;
 178                 for (j = 1; j <= i__1; ++j) {
 179                     ctemp = c__[j];
 180                     stemp = s[j];
 181                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
 182                         i__2 = *m;
 183                         for (i__ = 1; i__ <= i__2; ++i__) {
 184                             temp = a[i__ + (j + 1) * a_dim1];
 185                             a[i__ + (j + 1) * a_dim1] = ctemp * temp - stemp *
 186                                      a[i__ + j * a_dim1];
 187                             a[i__ + j * a_dim1] = stemp * temp + ctemp * a[
 188                                     i__ + j * a_dim1];
 189                         }
 190                     }
 191                 }
 192             } else if (*direct=='B' || *direct=='b') {
 193                 for (j = *n - 1; j >= 1; --j) {
 194                     ctemp = c__[j];
 195                     stemp = s[j];
 196                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
 197                         i__1 = *m;
 198                         for (i__ = 1; i__ <= i__1; ++i__) {
 199                             temp = a[i__ + (j + 1) * a_dim1];
 200                             a[i__ + (j + 1) * a_dim1] = ctemp * temp - stemp *
 201                                      a[i__ + j * a_dim1];
 202                             a[i__ + j * a_dim1] = stemp * temp + ctemp * a[
 203                                     i__ + j * a_dim1];
 204                         }
 205                     }
 206                 }
 207             }
 208         } else if (*pivot=='T' || *pivot=='t') {
 209             if (*direct=='F' || *direct=='f') {
 210                 i__1 = *n;
 211                 for (j = 2; j <= i__1; ++j) {
 212                     ctemp = c__[j - 1];
 213                     stemp = s[j - 1];
 214                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
 215                         i__2 = *m;
 216                         for (i__ = 1; i__ <= i__2; ++i__) {
 217                             temp = a[i__ + j * a_dim1];
 218                             a[i__ + j * a_dim1] = ctemp * temp - stemp * a[
 219                                     i__ + a_dim1];
 220                             a[i__ + a_dim1] = stemp * temp + ctemp * a[i__ +
 221                                     a_dim1];
 222                         }
 223                     }
 224                 }
 225             } else if (*direct=='B' || *direct=='b') {
 226                 for (j = *n; j >= 2; --j) {
 227                     ctemp = c__[j - 1];
 228                     stemp = s[j - 1];
 229                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
 230                         i__1 = *m;
 231                         for (i__ = 1; i__ <= i__1; ++i__) {
 232                             temp = a[i__ + j * a_dim1];
 233                             a[i__ + j * a_dim1] = ctemp * temp - stemp * a[
 234                                     i__ + a_dim1];
 235                             a[i__ + a_dim1] = stemp * temp + ctemp * a[i__ +
 236                                     a_dim1];
 237                         }
 238                     }
 239                 }
 240             }
 241         } else if (*pivot=='B' || *pivot=='b') {
 242             if (*direct=='F' || *direct=='f') {
 243                 i__1 = *n - 1;
 244                 for (j = 1; j <= i__1; ++j) {
 245                     ctemp = c__[j];
 246                     stemp = s[j];
 247                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
 248                         i__2 = *m;
 249                         for (i__ = 1; i__ <= i__2; ++i__) {
 250                             temp = a[i__ + j * a_dim1];
 251                             a[i__ + j * a_dim1] = stemp * a[i__ + *n * a_dim1]
 252                                      + ctemp * temp;
 253                             a[i__ + *n * a_dim1] = ctemp * a[i__ + *n *
 254                                     a_dim1] - stemp * temp;
 255                         }
 256                     }
 257                 }
 258             } else if (*direct=='B' || *direct=='b') {
 259                 for (j = *n - 1; j >= 1; --j) {
 260                     ctemp = c__[j];
 261                     stemp = s[j];
 262                     if (fabs(ctemp-1.0)>GMX_FLOAT_EPS || fabs(stemp)>GMX_FLOAT_MIN) {
 263                         i__1 = *m;
 264                         for (i__ = 1; i__ <= i__1; ++i__) {
 265                             temp = a[i__ + j * a_dim1];
 266                             a[i__ + j * a_dim1] = stemp * a[i__ + *n * a_dim1]
 267                                      + ctemp * temp;
 268                             a[i__ + *n * a_dim1] = ctemp * a[i__ + *n *
 269                                     a_dim1] - stemp * temp;
 270                         }
 271                     }
 272                 }
 273             }
 274         }
 275     }
 276
 277     return;
 278
 279 }
 280
 281