src/gmxlib/gmx_lapack/slarfb.c

   1 #include "gmx_blas.h"
   2 #include "gmx_lapack.h"
   3
   4
   5 void
   6 F77_FUNC(slarfb,SLARFB)(const char *side,
   7         const char *trans,
   8         const char *direct,
   9         const char *storev,
  10         int *m,
  11         int *n,
  12         int *k,
  13         float *v,
  14         int *ldv,
  15         float *t,
  16         int *ldt,
  17         float *c__,
  18         int *ldc,
  19         float *work,
  20         int *ldwork)
  21 {
  22     int c_dim1, c_offset, t_dim1, t_offset, v_dim1, v_offset, work_dim1,
  23             work_offset, i__1, i__2;
  24
  25     int i__, j;
  26     char transt[1];
  27     int c__1 = 1;
  28     float one = 1.0;
  29     float minusone = -1.0;
  30
  31     v_dim1 = *ldv;
  32     v_offset = 1 + v_dim1;
  33     v -= v_offset;
  34     t_dim1 = *ldt;
  35     t_offset = 1 + t_dim1;
  36     t -= t_offset;
  37     c_dim1 = *ldc;
  38     c_offset = 1 + c_dim1;
  39     c__ -= c_offset;
  40     work_dim1 = *ldwork;
  41     work_offset = 1 + work_dim1;
  42     work -= work_offset;
  43
  44     if (*m <= 0 || *n <= 0) {
  45         return;
  46     }
  47     if (*trans=='N' || *trans=='n') {
  48       *(unsigned char *)transt = 'T';
  49     } else {
  50         *(unsigned char *)transt = 'N';
  51     }
  52
  53     if (*storev=='C' || *storev=='c') {
  54
  55         if (*direct=='F' || *direct=='f') {
  56           if (*side=='l' || *side=='L') {
  57
  58                 i__1 = *k;
  59                 for (j = 1; j <= i__1; ++j) {
  60                     F77_FUNC(scopy,SCOPY)(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1],
  61                              &c__1);
  62                 }
  63
  64                 F77_FUNC(strmm,STRMM)("Right", "Lower", "No transpose", "Unit", n, k, &one,
  65                          &v[v_offset], ldv, &work[work_offset], ldwork);
  66                 if (*m > *k) {
  67
  68                     i__1 = *m - *k;
  69                     F77_FUNC(sgemm,SGEMM)("Transpose", "No transpose", n, k, &i__1, &one, &
  70                             c__[*k + 1 + c_dim1], ldc, &v[*k + 1 + v_dim1],
  71                             ldv, &one, &work[work_offset], ldwork);
  72                 }
  73
  74                 F77_FUNC(strmm,STRMM)("Right", "Upper", transt, "Non-unit", n, k, &one, &t[
  75                         t_offset], ldt, &work[work_offset], ldwork);
  76
  77                 if (*m > *k) {
  78                     i__1 = *m - *k;
  79                     F77_FUNC(sgemm,SGEMM)("No transpose", "Transpose", &i__1, n, k, &minusone, &
  80                             v[*k + 1 + v_dim1], ldv, &work[work_offset],
  81                             ldwork, &one, &c__[*k + 1 + c_dim1], ldc);
  82                 }
  83
  84                 F77_FUNC(strmm,STRMM)("Right", "Lower", "Transpose", "Unit", n, k, &one, &
  85                         v[v_offset], ldv, &work[work_offset], ldwork);
  86
  87                 i__1 = *k;
  88                 for (j = 1; j <= i__1; ++j) {
  89                     i__2 = *n;
  90                     for (i__ = 1; i__ <= i__2; ++i__) {
  91                         c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
  92                     }
  93                 }
  94
  95             } else if (*side=='r' || *side=='R') {
  96
  97                 i__1 = *k;
  98                 for (j = 1; j <= i__1; ++j) {
  99                     F77_FUNC(scopy,SCOPY)(m, &c__[j * c_dim1 + 1], &c__1, &work[j *
 100                             work_dim1 + 1], &c__1);
 101                 }
 102
 103                 F77_FUNC(strmm,STRMM)("Right", "Lower", "No transpose", "Unit", m, k, &one,
 104                          &v[v_offset], ldv, &work[work_offset], ldwork);
 105                 if (*n > *k) {
 106
 107                     i__1 = *n - *k;
 108                     F77_FUNC(sgemm,SGEMM)("No transpose", "No transpose", m, k, &i__1, &
 109                             one, &c__[(*k + 1) * c_dim1 + 1], ldc, &v[*k +
 110                             1 + v_dim1], ldv, &one, &work[work_offset],
 111                             ldwork);
 112                 }
 113
 114                 F77_FUNC(strmm,STRMM)("Right", "Upper", trans, "Non-unit", m, k, &one, &t[
 115                         t_offset], ldt, &work[work_offset], ldwork);
 116
 117                 if (*n > *k) {
 118                     i__1 = *n - *k;
 119                     F77_FUNC(sgemm,SGEMM)("No transpose", "Transpose", m, &i__1, k, &minusone, &
 120                             work[work_offset], ldwork, &v[*k + 1 + v_dim1],
 121                             ldv, &one, &c__[(*k + 1) * c_dim1 + 1], ldc);
 122                 }
 123
 124                 F77_FUNC(strmm,STRMM)("Right", "Lower", "Transpose", "Unit", m, k, &one, &
 125                         v[v_offset], ldv, &work[work_offset], ldwork);
 126
 127                 i__1 = *k;
 128                 for (j = 1; j <= i__1; ++j) {
 129                     i__2 = *m;
 130                     for (i__ = 1; i__ <= i__2; ++i__) {
 131                         c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
 132                     }
 133                 }
 134             }
 135
 136         } else {
 137
 138           if (*side=='l' || *side=='L') {
 139                 i__1 = *k;
 140                 for (j = 1; j <= i__1; ++j) {
 141                     F77_FUNC(scopy,SCOPY)(n, &c__[*m - *k + j + c_dim1], ldc, &work[j *
 142                             work_dim1 + 1], &c__1);
 143                 }
 144
 145                 F77_FUNC(strmm,STRMM)("Right", "Upper", "No transpose", "Unit", n, k, &one,
 146                          &v[*m - *k + 1 + v_dim1], ldv, &work[work_offset],
 147                         ldwork);
 148                 if (*m > *k) {
 149                     i__1 = *m - *k;
 150                     F77_FUNC(sgemm,SGEMM)("Transpose", "No transpose", n, k, &i__1, &one, &
 151                             c__[c_offset], ldc, &v[v_offset], ldv, &one, &
 152                             work[work_offset], ldwork);
 153                 }
 154
 155                 F77_FUNC(strmm,STRMM)("Right", "Lower", transt, "Non-unit", n, k, &one, &t[
 156                         t_offset], ldt, &work[work_offset], ldwork);
 157
 158                 if (*m > *k) {
 159
 160                     i__1 = *m - *k;
 161                     F77_FUNC(sgemm,SGEMM)("No transpose", "Transpose", &i__1, n, k, &minusone, &
 162                             v[v_offset], ldv, &work[work_offset], ldwork, &
 163                             one, &c__[c_offset], ldc)
 164                             ;
 165                 }
 166
 167                 F77_FUNC(strmm,STRMM)("Right", "Upper", "Transpose", "Unit", n, k, &one, &
 168                         v[*m - *k + 1 + v_dim1], ldv, &work[work_offset],
 169                         ldwork);
 170
 171                 i__1 = *k;
 172                 for (j = 1; j <= i__1; ++j) {
 173                     i__2 = *n;
 174                     for (i__ = 1; i__ <= i__2; ++i__) {
 175                         c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j *
 176                                 work_dim1];
 177                     }
 178                 }
 179
 180             } else if (*side=='r' || *side=='R') {
 181                 i__1 = *k;
 182                 for (j = 1; j <= i__1; ++j) {
 183                     F77_FUNC(scopy,SCOPY)(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
 184                             j * work_dim1 + 1], &c__1);
 185                 }
 186
 187                 F77_FUNC(strmm,STRMM)("Right", "Upper", "No transpose", "Unit", m, k, &one,
 188                          &v[*n - *k + 1 + v_dim1], ldv, &work[work_offset],
 189                         ldwork);
 190                 if (*n > *k) {
 191                     i__1 = *n - *k;
 192                     F77_FUNC(sgemm,SGEMM)("No transpose", "No transpose", m, k, &i__1, &
 193                             one, &c__[c_offset], ldc, &v[v_offset], ldv, &
 194                             one, &work[work_offset], ldwork);
 195                 }
 196                 F77_FUNC(strmm,STRMM)("Right", "Lower", trans, "Non-unit", m, k, &one, &t[
 197                         t_offset], ldt, &work[work_offset], ldwork);
 198                 if (*n > *k) {
 199                     i__1 = *n - *k;
 200                     F77_FUNC(sgemm,SGEMM)("No transpose", "Transpose", m, &i__1, k, &minusone, &
 201                             work[work_offset], ldwork, &v[v_offset], ldv, &
 202                             one, &c__[c_offset], ldc)
 203                             ;
 204                 }
 205                 F77_FUNC(strmm,STRMM)("Right", "Upper", "Transpose", "Unit", m, k, &one, &
 206                         v[*n - *k + 1 + v_dim1], ldv, &work[work_offset],
 207                         ldwork);
 208                 i__1 = *k;
 209                 for (j = 1; j <= i__1; ++j) {
 210                     i__2 = *m;
 211                     for (i__ = 1; i__ <= i__2; ++i__) {
 212                         c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j *
 213                                 work_dim1];
 214                     }
 215                 }
 216             }
 217         }
 218
 219     } else  if (*storev=='r' || *storev=='R') {
 220       if (*direct=='F' || *direct=='f') {
 221           if (*side=='l' || *side=='L') {
 222                 i__1 = *k;
 223                 for (j = 1; j <= i__1; ++j) {
 224                     F77_FUNC(scopy,SCOPY)(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1],
 225                              &c__1);
 226                 }
 227                 F77_FUNC(strmm,STRMM)("Right", "Upper", "Transpose", "Unit", n, k, &one, &
 228                         v[v_offset], ldv, &work[work_offset], ldwork);
 229                 if (*m > *k) {
 230                     i__1 = *m - *k;
 231                     F77_FUNC(sgemm,SGEMM)("Transpose", "Transpose", n, k, &i__1, &one, &
 232                             c__[*k + 1 + c_dim1], ldc, &v[(*k + 1) * v_dim1 +
 233                             1], ldv, &one, &work[work_offset], ldwork);
 234                 }
 235
 236                 F77_FUNC(strmm,STRMM)("Right", "Upper", transt, "Non-unit", n, k, &one, &t[
 237                         t_offset], ldt, &work[work_offset], ldwork);
 238                 if (*m > *k) {
 239
 240                     i__1 = *m - *k;
 241                     F77_FUNC(sgemm,SGEMM)("Transpose", "Transpose", &i__1, n, k, &minusone, &v[(
 242                             *k + 1) * v_dim1 + 1], ldv, &work[work_offset],
 243                             ldwork, &one, &c__[*k + 1 + c_dim1], ldc);
 244                 }
 245
 246                 F77_FUNC(strmm,STRMM)("Right", "Upper", "No transpose", "Unit", n, k, &one,
 247                          &v[v_offset], ldv, &work[work_offset], ldwork);
 248
 249                 i__1 = *k;
 250                 for (j = 1; j <= i__1; ++j) {
 251                     i__2 = *n;
 252                     for (i__ = 1; i__ <= i__2; ++i__) {
 253                         c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
 254                     }
 255                 }
 256
 257             } else if (*side=='r' || *side=='R') {
 258
 259                 i__1 = *k;
 260                 for (j = 1; j <= i__1; ++j) {
 261                     F77_FUNC(scopy,SCOPY)(m, &c__[j * c_dim1 + 1], &c__1, &work[j *
 262                             work_dim1 + 1], &c__1);
 263                 }
 264
 265                 F77_FUNC(strmm,STRMM)("Right", "Upper", "Transpose", "Unit", m, k, &one, &
 266                         v[v_offset], ldv, &work[work_offset], ldwork);
 267                 if (*n > *k) {
 268
 269                     i__1 = *n - *k;
 270                     F77_FUNC(sgemm,SGEMM)("No transpose", "Transpose", m, k, &i__1, &one, &
 271                             c__[(*k + 1) * c_dim1 + 1], ldc, &v[(*k + 1) *
 272                             v_dim1 + 1], ldv, &one, &work[work_offset],
 273                             ldwork);
 274                 }
 275
 276                 F77_FUNC(strmm,STRMM)("Right", "Upper", trans, "Non-unit", m, k, &one, &t[
 277                         t_offset], ldt, &work[work_offset], ldwork);
 278
 279                 if (*n > *k) {
 280
 281                     i__1 = *n - *k;
 282                     F77_FUNC(sgemm,SGEMM)("No transpose", "No transpose", m, &i__1, k, &
 283                             minusone, &work[work_offset], ldwork, &v[(*k + 1) *
 284                             v_dim1 + 1], ldv, &one, &c__[(*k + 1) * c_dim1
 285                             + 1], ldc);
 286                 }
 287                 F77_FUNC(strmm,STRMM)("Right", "Upper", "No transpose", "Unit", m, k, &one,
 288                          &v[v_offset], ldv, &work[work_offset], ldwork);
 289                 i__1 = *k;
 290                 for (j = 1; j <= i__1; ++j) {
 291                     i__2 = *m;
 292                     for (i__ = 1; i__ <= i__2; ++i__) {
 293                         c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
 294                     }
 295                 }
 296
 297             }
 298
 299         } else {
 300
 301             if (*side=='l' || *side=='L') {
 302
 303                 i__1 = *k;
 304                 for (j = 1; j <= i__1; ++j) {
 305                     F77_FUNC(scopy,SCOPY)(n, &c__[*m - *k + j + c_dim1], ldc, &work[j *
 306                             work_dim1 + 1], &c__1);
 307                 }
 308
 309                 F77_FUNC(strmm,STRMM)("Right", "Lower", "Transpose", "Unit", n, k, &one, &
 310                         v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
 311                         , ldwork);
 312                 if (*m > *k) {
 313
 314                     i__1 = *m - *k;
 315                     F77_FUNC(sgemm,SGEMM)("Transpose", "Transpose", n, k, &i__1, &one, &
 316                             c__[c_offset], ldc, &v[v_offset], ldv, &one, &
 317                             work[work_offset], ldwork);
 318                 }
 319
 320                 F77_FUNC(strmm,STRMM)("Right", "Lower", transt, "Non-unit", n, k, &one, &t[
 321                         t_offset], ldt, &work[work_offset], ldwork);
 322
 323                 if (*m > *k) {
 324
 325                     i__1 = *m - *k;
 326                     F77_FUNC(sgemm,SGEMM)("Transpose", "Transpose", &i__1, n, k, &minusone, &v[
 327                             v_offset], ldv, &work[work_offset], ldwork, &
 328                             one, &c__[c_offset], ldc);
 329                 }
 330
 331                 F77_FUNC(strmm,STRMM)("Right", "Lower", "No transpose", "Unit", n, k, &one,
 332                          &v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[
 333                         work_offset], ldwork);
 334
 335                 i__1 = *k;
 336                 for (j = 1; j <= i__1; ++j) {
 337                     i__2 = *n;
 338                     for (i__ = 1; i__ <= i__2; ++i__) {
 339                         c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j *
 340                                 work_dim1];
 341                     }
 342                 }
 343
 344             } else if (*side=='r' || *side=='R') {
 345
 346                 i__1 = *k;
 347                 for (j = 1; j <= i__1; ++j) {
 348                     F77_FUNC(scopy,SCOPY)(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
 349                             j * work_dim1 + 1], &c__1);
 350                 }
 351
 352                 F77_FUNC(strmm,STRMM)("Right", "Lower", "Transpose", "Unit", m, k, &one, &
 353                         v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
 354                         , ldwork);
 355                 if (*n > *k) {
 356
 357                     i__1 = *n - *k;
 358                     F77_FUNC(sgemm,SGEMM)("No transpose", "Transpose", m, k, &i__1, &one, &
 359                             c__[c_offset], ldc, &v[v_offset], ldv, &one, &
 360                             work[work_offset], ldwork);
 361                 }
 362
 363                 F77_FUNC(strmm,STRMM)("Right", "Lower", trans, "Non-unit", m, k, &one, &t[
 364                         t_offset], ldt, &work[work_offset], ldwork);
 365
 366                 if (*n > *k) {
 367
 368                     i__1 = *n - *k;
 369                     F77_FUNC(sgemm,SGEMM)("No transpose", "No transpose", m, &i__1, k, &
 370                             minusone, &work[work_offset], ldwork, &v[v_offset],
 371                             ldv, &one, &c__[c_offset], ldc);
 372                 }
 373
 374                 F77_FUNC(strmm,STRMM)("Right", "Lower", "No transpose", "Unit", m, k, &one,
 375                          &v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[
 376                         work_offset], ldwork);
 377
 378                 i__1 = *k;
 379                 for (j = 1; j <= i__1; ++j) {
 380                     i__2 = *m;
 381                     for (i__ = 1; i__ <= i__2; ++i__) {
 382                         c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j *
 383                                 work_dim1];
 384                     }
 385                 }
 386
 387             }
 388
 389         }
 390     }
 391
 392     return;
 393
 394
 395 }
 396