src/gmxlib/gmx_lapack/slasdq.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, 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 <ctype.h>
  36
  37 #include "gmx_blas.h"
  38 #include "gmx_lapack.h"
  39
  40
  41 void
  42 F77_FUNC(slasdq,SLASDQ)(const char *uplo,
  43                         int *sqre,
  44                         int *n,
  45                         int *ncvt,
  46                         int *nru,
  47                         int *ncc,
  48                         float *d__,
  49                         float *e,
  50                         float *vt,
  51                         int *ldvt,
  52                         float *u,
  53                         int *ldu,
  54                         float *c__,
  55                         int *ldc,
  56                         float *work,
  57                         int *info)
  58 {
  59     const char xuplo=toupper(*uplo);
  60     int c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
  61             i__2;
  62     int c__1 = 1;
  63     int itmp1,itmp2;
  64     int i__, j;
  65     float r__, cs, sn;
  66     int np1, isub;
  67     float smin;
  68     int sqre1;
  69     int iuplo;
  70     int rotate;
  71
  72     --d__;
  73     --e;
  74     vt_dim1 = *ldvt;
  75     vt_offset = 1 + vt_dim1;
  76     vt -= vt_offset;
  77     u_dim1 = *ldu;
  78     u_offset = 1 + u_dim1;
  79     u -= u_offset;
  80     c_dim1 = *ldc;
  81     c_offset = 1 + c_dim1;
  82     c__ -= c_offset;
  83     --work;
  84
  85     *info = 0;
  86     iuplo = 0;
  87     if (xuplo == 'U') {
  88         iuplo = 1;
  89     }
  90     if (xuplo == 'L') {
  91         iuplo = 2;
  92     }
  93
  94     itmp1 = (*n > 1) ? *n : 1;
  95     itmp2 = (*nru > 1) ? *nru : 1;
  96     if (iuplo == 0) {
  97         *info = -1;
  98     } else if (*sqre < 0 || *sqre > 1) {
  99         *info = -2;
 100     } else if (*n < 0) {
 101         *info = -3;
 102     } else if (*ncvt < 0) {
 103         *info = -4;
 104     } else if (*nru < 0) {
 105         *info = -5;
 106     } else if (*ncc < 0) {
 107         *info = -6;
 108     } else if ( (*ncvt == 0 && *ldvt < 1) || (*ncvt > 0 && *ldvt < itmp1)) {
 109         *info = -10;
 110     } else if (*ldu < itmp2) {
 111         *info = -12;
 112     } else if ((*ncc == 0 && *ldc < 1) || (*ncc > 0 && *ldc < itmp1)) {
 113         *info = -14;
 114     }
 115     if (*info != 0) {
 116         return;
 117     }
 118     if (*n == 0) {
 119         return;
 120     }
 121
 122     rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
 123     np1 = *n + 1;
 124     sqre1 = *sqre;
 125
 126     if (iuplo == 1 && sqre1 == 1) {
 127         i__1 = *n - 1;
 128         for (i__ = 1; i__ <= i__1; ++i__) {
 129             F77_FUNC(slartg,SLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
 130             d__[i__] = r__;
 131             e[i__] = sn * d__[i__ + 1];
 132             d__[i__ + 1] = cs * d__[i__ + 1];
 133             if (rotate) {
 134                 work[i__] = cs;
 135                 work[*n + i__] = sn;
 136             }
 137         }
 138         F77_FUNC(slartg,SLARTG)(&d__[*n], &e[*n], &cs, &sn, &r__);
 139         d__[*n] = r__;
 140         e[*n] = 0.f;
 141         if (rotate) {
 142             work[*n] = cs;
 143             work[*n + *n] = sn;
 144         }
 145         iuplo = 2;
 146         sqre1 = 0;
 147
 148         if (*ncvt > 0) {
 149             F77_FUNC(slasr,SLASR)("L", "V", "F", &np1, ncvt, &work[1], &work[np1], &vt[
 150                     vt_offset], ldvt);
 151         }
 152     }
 153     if (iuplo == 2) {
 154         i__1 = *n - 1;
 155         for (i__ = 1; i__ <= i__1; ++i__) {
 156             F77_FUNC(slartg,SLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
 157             d__[i__] = r__;
 158             e[i__] = sn * d__[i__ + 1];
 159             d__[i__ + 1] = cs * d__[i__ + 1];
 160             if (rotate) {
 161                 work[i__] = cs;
 162                 work[*n + i__] = sn;
 163             }
 164         }
 165
 166         if (sqre1 == 1) {
 167             F77_FUNC(slartg,SLARTG)(&d__[*n], &e[*n], &cs, &sn, &r__);
 168             d__[*n] = r__;
 169             if (rotate) {
 170                 work[*n] = cs;
 171                 work[*n + *n] = sn;
 172             }
 173         }
 174         if (*nru > 0) {
 175             if (sqre1 == 0) {
 176                 F77_FUNC(slasr,SLASR)("R", "V", "F", nru, n, &work[1], &work[np1], &u[
 177                         u_offset], ldu);
 178             } else {
 179                 F77_FUNC(slasr,SLASR)("R", "V", "F", nru, &np1, &work[1], &work[np1], &u[
 180                         u_offset], ldu);
 181             }
 182         }
 183         if (*ncc > 0) {
 184             if (sqre1 == 0) {
 185                 F77_FUNC(slasr,SLASR)("L", "V", "F", n, ncc, &work[1], &work[np1], &c__[
 186                         c_offset], ldc);
 187             } else {
 188                 F77_FUNC(slasr,SLASR)("L", "V", "F", &np1, ncc, &work[1], &work[np1], &c__[
 189                         c_offset], ldc);
 190             }
 191         }
 192     }
 193
 194     F77_FUNC(sbdsqr,SBDSQR)("U", n, ncvt, nru, ncc, &d__[1], &e[1], &vt[vt_offset], ldvt, &u[
 195             u_offset], ldu, &c__[c_offset], ldc, &work[1], info);
 196
 197     i__1 = *n;
 198     for (i__ = 1; i__ <= i__1; ++i__) {
 199
 200         isub = i__;
 201         smin = d__[i__];
 202         i__2 = *n;
 203         for (j = i__ + 1; j <= i__2; ++j) {
 204             if (d__[j] < smin) {
 205                 isub = j;
 206                 smin = d__[j];
 207             }
 208         }
 209         if (isub != i__) {
 210             d__[isub] = d__[i__];
 211             d__[i__] = smin;
 212             if (*ncvt > 0) {
 213                 F77_FUNC(sswap,SSWAP)(ncvt, &vt[isub + vt_dim1], ldvt, &vt[i__ + vt_dim1],
 214                         ldvt);
 215             }
 216             if (*nru > 0) {
 217                 F77_FUNC(sswap,SSWAP)(nru, &u[isub * u_dim1 + 1], &c__1, &u[i__ * u_dim1 + 1]
 218                         , &c__1);
 219             }
 220             if (*ncc > 0) {
 221                 F77_FUNC(sswap,SSWAP)(ncc, &c__[isub + c_dim1], ldc, &c__[i__ + c_dim1], ldc)
 222                         ;
 223             }
 224         }
 225     }
 226
 227     return;
 228 }
 229
 230