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