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