src/gmxlib/gmx_lapack/dlabrd.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 #include "gmx_blas.h"
  37 #include "gmx_lapack.h"
  38
  39
  40 void
  41 F77_FUNC(dlabrd,DLABRD)(int *m,
  42         int *n,
  43         int *nb,
  44         double *a,
  45         int *lda,
  46         double *d__,
  47         double *e,
  48         double *tauq,
  49         double *taup,
  50         double *x,
  51         int *ldx,
  52         double *y,
  53         int *ldy)
  54 {
  55     int a_dim1, a_offset, x_dim1, x_offset, y_dim1, y_offset;
  56     int i__1, i__2, i__3;
  57     double one = 1.0;
  58     double minusone = -1.0;
  59     double zero = 0.0;
  60     int c__1 = 1;
  61     int i__;
  62
  63     a_dim1 = *lda;
  64     a_offset = 1 + a_dim1;
  65     a -= a_offset;
  66     --d__;
  67     --e;
  68     --tauq;
  69     --taup;
  70     x_dim1 = *ldx;
  71     x_offset = 1 + x_dim1;
  72     x -= x_offset;
  73     y_dim1 = *ldy;
  74     y_offset = 1 + y_dim1;
  75     y -= y_offset;
  76
  77     if (*m <= 0 || *n <= 0) {
  78         return;
  79     }
  80
  81     if (*m >= *n) {
  82
  83         i__1 = *nb;
  84         for (i__ = 1; i__ <= i__1; ++i__) {
  85
  86             i__2 = *m - i__ + 1;
  87             i__3 = i__ - 1;
  88             F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &a[i__ + a_dim1], lda,
  89                      &y[i__ + y_dim1], ldy, &one, &a[i__ + i__ * a_dim1], &c__1);
  90             i__2 = *m - i__ + 1;
  91             i__3 = i__ - 1;
  92             F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &x[i__ + x_dim1], ldx,
  93                    &a[i__*a_dim1+1],&c__1,&one,&a[i__+i__*a_dim1],&c__1);
  94
  95             i__2 = *m - i__ + 1;
  96             i__3 = i__ + 1;
  97             if(*m<i__3)
  98               i__3 = *m;
  99             F77_FUNC(dlarfg,DLARFG)(&i__2, &a[i__ + i__ * a_dim1], &a[i__3 + i__ * a_dim1],
 100                     &c__1, &tauq[i__]);
 101             d__[i__] = a[i__ + i__ * a_dim1];
 102             if (i__ < *n) {
 103                 a[i__ + i__ * a_dim1] = 1.;
 104
 105                 i__2 = *m - i__ + 1;
 106                 i__3 = *n - i__;
 107                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &a[i__ + (i__ + 1) *
 108                         a_dim1], lda, &a[i__ + i__ * a_dim1], &c__1, &zero, &
 109                         y[i__ + 1 + i__ * y_dim1], &c__1);
 110                 i__2 = *m - i__ + 1;
 111                 i__3 = i__ - 1;
 112                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &a[i__ + a_dim1],
 113                         lda, &a[i__ + i__ * a_dim1], &c__1, &zero, &y[i__ *
 114                         y_dim1 + 1], &c__1);
 115                 i__2 = *n - i__;
 116                 i__3 = i__ - 1;
 117                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &y[i__ + 1 +
 118                         y_dim1], ldy, &y[i__ * y_dim1 + 1], &c__1, &one, &y[
 119                         i__ + 1 + i__ * y_dim1], &c__1);
 120                 i__2 = *m - i__ + 1;
 121                 i__3 = i__ - 1;
 122                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &x[i__ + x_dim1],
 123                         ldx, &a[i__ + i__ * a_dim1], &c__1, &zero, &y[i__ *
 124                         y_dim1 + 1], &c__1);
 125                 i__2 = i__ - 1;
 126                 i__3 = *n - i__;
 127                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &minusone, &a[(i__ + 1) *
 128                         a_dim1 + 1], lda, &y[i__ * y_dim1 + 1], &c__1, &one,
 129                         &y[i__ + 1 + i__ * y_dim1], &c__1);
 130                 i__2 = *n - i__;
 131                 F77_FUNC(dscal,DSCAL)(&i__2, &tauq[i__], &y[i__ + 1 + i__ * y_dim1], &c__1);
 132
 133                 i__2 = *n - i__;
 134                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__, &minusone, &y[i__ + 1 +
 135                         y_dim1], ldy, &a[i__ + a_dim1], lda, &one, &a[i__ + (
 136                         i__ + 1) * a_dim1], lda);
 137                 i__2 = i__ - 1;
 138                 i__3 = *n - i__;
 139                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &minusone, &a[(i__ + 1) *
 140                         a_dim1 + 1], lda, &x[i__ + x_dim1], ldx, &one, &a[
 141                         i__ + (i__ + 1) * a_dim1], lda);
 142
 143                 i__2 = *n - i__;
 144                 i__3 = i__ + 2;
 145                 if(*n<i__3)
 146                   i__3 = *n;
 147                 F77_FUNC(dlarfg,DLARFG)(&i__2, &a[i__ + (i__ + 1) * a_dim1],
 148                         &a[i__ + i__3 * a_dim1], lda, &taup[i__]);
 149                 e[i__] = a[i__ + (i__ + 1) * a_dim1];
 150                 a[i__ + (i__ + 1) * a_dim1] = 1.;
 151
 152                 i__2 = *m - i__;
 153                 i__3 = *n - i__;
 154                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &one, &a[i__ + 1 + (i__
 155                         + 1) * a_dim1], lda, &a[i__ + (i__ + 1) * a_dim1],
 156                         lda, &zero, &x[i__ + 1 + i__ * x_dim1], &c__1);
 157                 i__2 = *n - i__;
 158                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__, &one, &y[i__ + 1 + y_dim1],
 159                         ldy, &a[i__ + (i__ + 1) * a_dim1], lda, &zero, &x[
 160                         i__ * x_dim1 + 1], &c__1);
 161                 i__2 = *m - i__;
 162                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__, &minusone, &a[i__ + 1 +
 163                         a_dim1], lda, &x[i__ * x_dim1 + 1], &c__1, &one, &x[
 164                         i__ + 1 + i__ * x_dim1], &c__1);
 165                 i__2 = i__ - 1;
 166                 i__3 = *n - i__;
 167                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &one, &a[(i__ + 1) *
 168                         a_dim1 + 1], lda, &a[i__ + (i__ + 1) * a_dim1], lda, &
 169                         zero, &x[i__ * x_dim1 + 1], &c__1);
 170                 i__2 = *m - i__;
 171                 i__3 = i__ - 1;
 172                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &x[i__ + 1 +
 173                         x_dim1], ldx, &x[i__ * x_dim1 + 1], &c__1, &one, &x[
 174                         i__ + 1 + i__ * x_dim1], &c__1);
 175                 i__2 = *m - i__;
 176                 F77_FUNC(dscal,DSCAL)(&i__2, &taup[i__], &x[i__ + 1 + i__ * x_dim1], &c__1);
 177             }
 178         }
 179     } else {
 180
 181         i__1 = *nb;
 182         for (i__ = 1; i__ <= i__1; ++i__) {
 183
 184             i__2 = *n - i__ + 1;
 185             i__3 = i__ - 1;
 186             F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &y[i__ + y_dim1], ldy,
 187                      &a[i__ + a_dim1], lda, &one, &a[i__ + i__ * a_dim1],lda);
 188             i__2 = i__ - 1;
 189             i__3 = *n - i__ + 1;
 190             F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &minusone, &a[i__ * a_dim1 + 1],
 191                     lda, &x[i__ + x_dim1], ldx, &one,&a[i__+i__*a_dim1],lda);
 192
 193             i__2 = *n - i__ + 1;
 194             i__3 = i__ + 1;
 195             if(*n<i__3)
 196               i__3 = *n;
 197             F77_FUNC(dlarfg,DLARFG)(&i__2, &a[i__ + i__ * a_dim1],
 198                     &a[i__ + i__3 * a_dim1], lda, &taup[i__]);
 199             d__[i__] = a[i__ + i__ * a_dim1];
 200             if (i__ < *m) {
 201                 a[i__ + i__ * a_dim1] = 1.;
 202
 203                 i__2 = *m - i__;
 204                 i__3 = *n - i__ + 1;
 205                 F77_FUNC(dgemv,DGEMV)("No transpose",&i__2,&i__3,&one,&a[i__+1+i__*a_dim1],
 206                        lda, &a[i__ + i__ * a_dim1], lda, &zero,
 207                        &x[i__ + 1 + i__ * x_dim1], &c__1);
 208                 i__2 = *n - i__ + 1;
 209                 i__3 = i__ - 1;
 210                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &y[i__ + y_dim1],
 211                         ldy, &a[i__ + i__ * a_dim1], lda, &zero, &x[i__ *
 212                         x_dim1 + 1], &c__1);
 213                 i__2 = *m - i__;
 214                 i__3 = i__ - 1;
 215                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &a[i__ + 1 +
 216                         a_dim1], lda, &x[i__ * x_dim1 + 1], &c__1, &one, &x[
 217                         i__ + 1 + i__ * x_dim1], &c__1);
 218                 i__2 = i__ - 1;
 219                 i__3 = *n - i__ + 1;
 220                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &one, &a[i__ * a_dim1 +
 221                         1], lda, &a[i__ + i__ * a_dim1], lda, &zero, &x[i__ *
 222                          x_dim1 + 1], &c__1);
 223                 i__2 = *m - i__;
 224                 i__3 = i__ - 1;
 225                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &x[i__ + 1 +
 226                         x_dim1], ldx, &x[i__ * x_dim1 + 1], &c__1, &one, &x[
 227                         i__ + 1 + i__ * x_dim1], &c__1);
 228                 i__2 = *m - i__;
 229                 F77_FUNC(dscal,DSCAL)(&i__2, &taup[i__], &x[i__ + 1 + i__ * x_dim1], &c__1);
 230
 231                 i__2 = *m - i__;
 232                 i__3 = i__ - 1;
 233                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &a[i__ + 1 +
 234                         a_dim1], lda, &y[i__ + y_dim1], ldy, &one, &a[i__ +
 235                         1 + i__ * a_dim1], &c__1);
 236                 i__2 = *m - i__;
 237                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__, &minusone, &x[i__ + 1 +
 238                         x_dim1], ldx, &a[i__ * a_dim1 + 1], &c__1, &one, &a[
 239                         i__ + 1 + i__ * a_dim1], &c__1);
 240
 241                 i__2 = *m - i__;
 242                 i__3 = i__ + 2;
 243                 if(*m<i__3)
 244                   i__3 = *m;
 245                 F77_FUNC(dlarfg,DLARFG)(&i__2, &a[i__ + 1 + i__ * a_dim1],
 246                         &a[i__3 + i__ * a_dim1], &c__1, &tauq[i__]);
 247                 e[i__] = a[i__ + 1 + i__ * a_dim1];
 248                 a[i__ + 1 + i__ * a_dim1] = 1.;
 249
 250                 i__2 = *m - i__;
 251                 i__3 = *n - i__;
 252                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &a[i__ + 1 + (i__ +
 253                         1) * a_dim1], lda, &a[i__ + 1 + i__ * a_dim1], &c__1,
 254                         &zero, &y[i__ + 1 + i__ * y_dim1], &c__1);
 255                 i__2 = *m - i__;
 256                 i__3 = i__ - 1;
 257                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &a[i__ + 1 + a_dim1],
 258                          lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &zero, &y[
 259                         i__ * y_dim1 + 1], &c__1);
 260                 i__2 = *n - i__;
 261                 i__3 = i__ - 1;
 262                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &y[i__ + 1 +
 263                         y_dim1], ldy, &y[i__ * y_dim1 + 1], &c__1, &one, &y[
 264                         i__ + 1 + i__ * y_dim1], &c__1);
 265                 i__2 = *m - i__;
 266                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__, &one, &x[i__ + 1 + x_dim1],
 267                         ldx, &a[i__ + 1 + i__ * a_dim1], &c__1, &zero, &y[
 268                         i__ * y_dim1 + 1], &c__1);
 269                 i__2 = *n - i__;
 270                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__, &i__2, &minusone, &a[(i__ + 1) * a_dim1
 271                         + 1], lda, &y[i__ * y_dim1 + 1], &c__1, &one, &y[i__
 272                         + 1 + i__ * y_dim1], &c__1);
 273                 i__2 = *n - i__;
 274                 F77_FUNC(dscal,DSCAL)(&i__2, &tauq[i__], &y[i__ + 1 + i__ * y_dim1], &c__1);
 275             }
 276         }
 277     }
 278     return;
 279 }
 280