src/gromacs/linearalgebra/gmx_blas/dtrmm.c

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