src/gmxlib/gmx_lapack/ssytrd.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 #include "lapack_limits.h"
  38
  39 void
  40 F77_FUNC(ssytrd,SSYTRD)(const char *uplo, int *n, float *a, int *
  41         lda, float *d__, float *e, float *tau, float *
  42         work, int *lwork, int *info)
  43 {
  44     /* System generated locals */
  45     int a_dim1, a_offset, i__1, i__2, i__3;
  46
  47     /* Local variables */
  48     int i__, j, nb, kk, nx, iws;
  49     int nbmin, iinfo;
  50     int upper;
  51     int ldwork, lwkopt;
  52     int lquery;
  53     float c_b22 = -1.;
  54     float c_b23 = 1.;
  55
  56
  57     /* Parameter adjustments */
  58     a_dim1 = *lda;
  59     a_offset = 1 + a_dim1;
  60     a -= a_offset;
  61     --d__;
  62     --e;
  63     --tau;
  64     --work;
  65
  66     /* Function Body */
  67     *info = 0;
  68     upper = (*uplo=='U' || *uplo=='u');
  69     lquery = (*lwork == -1);
  70
  71     if (! upper && ! (*uplo=='L' || *uplo=='l')) {
  72         *info = -1;
  73     } else if (*n < 0) {
  74         *info = -2;
  75     } else if (*lda < ((1>*n) ? 1 : *n)) {
  76         *info = -4;
  77     } else if (*lwork < 1 && ! lquery) {
  78         *info = -9;
  79     }
  80
  81     if (*info == 0) {
  82
  83       nb = DSYTRD_BLOCKSIZE;
  84       lwkopt = *n * nb;
  85       work[1] = (float) lwkopt;
  86     } else
  87       return;
  88
  89     if (lquery)
  90       return;
  91
  92     if (*n == 0) {
  93         work[1] = 1.;
  94         return;
  95     }
  96
  97     nx = *n;
  98     iws = 1;
  99     if (nb > 1 && nb < *n) {
 100
 101         nx = DSYTRD_CROSSOVER;
 102         if (nx < *n) {
 103
 104             ldwork = *n;
 105             iws = ldwork * nb;
 106             if (*lwork < iws) {
 107
 108                 i__1 = *lwork / ldwork;
 109                 nb = (i__1>1) ? i__1 : 1;
 110                 nbmin = DSYTRD_MINBLOCKSIZE;
 111                 if (nb < nbmin) {
 112                     nx = *n;
 113                 }
 114             }
 115         } else {
 116             nx = *n;
 117         }
 118     } else {
 119         nb = 1;
 120     }
 121
 122     if (upper) {
 123
 124         kk = *n - (*n - nx + nb - 1) / nb * nb;
 125         i__1 = kk + 1;
 126         i__2 = -nb;
 127         for (i__ = *n - nb + 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
 128                 i__2) {
 129
 130             i__3 = i__ + nb - 1;
 131             F77_FUNC(slatrd,SLATRD)(uplo, &i__3, &nb, &a[a_offset], lda, &e[1], &tau[1], &
 132                     work[1], &ldwork);
 133
 134             i__3 = i__ - 1;
 135             F77_FUNC(ssyr2k,SSYR2K)(uplo, "No transpose", &i__3, &nb, &c_b22, &a[i__ * a_dim1
 136                     + 1], lda, &work[1], &ldwork, &c_b23, &a[a_offset], lda);
 137
 138             i__3 = i__ + nb - 1;
 139             for (j = i__; j <= i__3; ++j) {
 140                 a[j - 1 + j * a_dim1] = e[j - 1];
 141                 d__[j] = a[j + j * a_dim1];
 142
 143             }
 144
 145         }
 146
 147         F77_FUNC(ssytd2,SSYTD2)(uplo, &kk, &a[a_offset], lda, &d__[1], &e[1], &tau[1], &iinfo);
 148     } else {
 149
 150         i__2 = *n - nx;
 151         i__1 = nb;
 152         for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
 153
 154
 155             i__3 = *n - i__ + 1;
 156             F77_FUNC(slatrd,SLATRD)(uplo, &i__3, &nb, &a[i__ + i__ * a_dim1], lda, &e[i__], &
 157                     tau[i__], &work[1], &ldwork);
 158
 159             i__3 = *n - i__ - nb + 1;
 160             F77_FUNC(ssyr2k,SSYR2K)(uplo, "No transpose", &i__3, &nb, &c_b22, &a[i__ + nb +
 161                     i__ * a_dim1], lda, &work[nb + 1], &ldwork, &c_b23, &a[
 162                     i__ + nb + (i__ + nb) * a_dim1], lda);
 163
 164
 165             i__3 = i__ + nb - 1;
 166             for (j = i__; j <= i__3; ++j) {
 167                 a[j + 1 + j * a_dim1] = e[j];
 168                 d__[j] = a[j + j * a_dim1];
 169
 170             }
 171
 172         }
 173
 174
 175         i__1 = *n - i__ + 1;
 176         F77_FUNC(ssytd2,SSYTD2)(uplo, &i__1, &a[i__ + i__ * a_dim1], lda, &d__[i__], &e[i__],
 177                 &tau[i__], &iinfo);
 178     }
 179
 180     work[1] = (float) lwkopt;
 181     return;
 182
 183 }
 184
 185