src/gmxlib/gmx_lapack/slasd0.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 "gmx_lapack.h"
  36
  37 void
  38 F77_FUNC(slasd0,SLASD0)(int *n,
  39         int *sqre,
  40         float *d__,
  41         float *e,
  42         float *u,
  43         int *ldu,
  44         float *vt,
  45         int *ldvt,
  46         int *smlsiz,
  47         int *iwork,
  48         float *work,
  49         int *info)
  50 {
  51     int u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
  52
  53     int i__, j, m, i1, ic, lf, nd, ll, nl, nr, im1, ncc, nlf, nrf,
  54             iwk, lvl, ndb1, nlp1, nrp1;
  55     float beta;
  56     int idxq, nlvl;
  57     float alpha;
  58     int inode, ndiml, idxqc, ndimr, itemp, sqrei;
  59     int c__0 = 0;
  60
  61
  62     --d__;
  63     --e;
  64     u_dim1 = *ldu;
  65     u_offset = 1 + u_dim1;
  66     u -= u_offset;
  67     vt_dim1 = *ldvt;
  68     vt_offset = 1 + vt_dim1;
  69     vt -= vt_offset;
  70     --iwork;
  71     --work;
  72
  73     *info = 0;
  74
  75     if (*n < 0) {
  76         *info = -1;
  77     } else if (*sqre < 0 || *sqre > 1) {
  78         *info = -2;
  79     }
  80
  81     m = *n + *sqre;
  82
  83     if (*ldu < *n) {
  84         *info = -6;
  85     } else if (*ldvt < m) {
  86         *info = -8;
  87     } else if (*smlsiz < 3) {
  88         *info = -9;
  89     }
  90     if (*info != 0) {
  91         i__1 = -(*info);
  92         return;
  93     }
  94
  95     if (*n <= *smlsiz) {
  96         F77_FUNC(slasdq,SLASDQ)("U", sqre, n, &m, n, &c__0, &d__[1], &e[1], &vt[vt_offset],
  97                 ldvt, &u[u_offset], ldu, &u[u_offset], ldu, &work[1], info);
  98         return;
  99     }
 100
 101     inode = 1;
 102     ndiml = inode + *n;
 103     ndimr = ndiml + *n;
 104     idxq = ndimr + *n;
 105     iwk = idxq + *n;
 106     F77_FUNC(slasdt,SLASDT)(n, &nlvl, &nd, &iwork[inode], &iwork[ndiml], &iwork[ndimr],
 107             smlsiz);
 108
 109     ndb1 = (nd + 1) / 2;
 110     ncc = 0;
 111     i__1 = nd;
 112     for (i__ = ndb1; i__ <= i__1; ++i__) {
 113
 114         i1 = i__ - 1;
 115         ic = iwork[inode + i1];
 116         nl = iwork[ndiml + i1];
 117         nlp1 = nl + 1;
 118         nr = iwork[ndimr + i1];
 119         nrp1 = nr + 1;
 120         nlf = ic - nl;
 121         nrf = ic + 1;
 122         sqrei = 1;
 123         F77_FUNC(slasdq,SLASDQ)("U", &sqrei, &nl, &nlp1, &nl, &ncc, &d__[nlf], &e[nlf], &vt[
 124                 nlf + nlf * vt_dim1], ldvt, &u[nlf + nlf * u_dim1], ldu, &u[
 125                 nlf + nlf * u_dim1], ldu, &work[1], info);
 126         if (*info != 0) {
 127             return;
 128         }
 129         itemp = idxq + nlf - 2;
 130         i__2 = nl;
 131         for (j = 1; j <= i__2; ++j) {
 132             iwork[itemp + j] = j;
 133         }
 134         if (i__ == nd) {
 135             sqrei = *sqre;
 136         } else {
 137             sqrei = 1;
 138         }
 139         nrp1 = nr + sqrei;
 140         F77_FUNC(slasdq,SLASDQ)("U", &sqrei, &nr, &nrp1, &nr, &ncc, &d__[nrf], &e[nrf], &vt[
 141                 nrf + nrf * vt_dim1], ldvt, &u[nrf + nrf * u_dim1], ldu, &u[
 142                 nrf + nrf * u_dim1], ldu, &work[1], info);
 143         if (*info != 0) {
 144             return;
 145         }
 146         itemp = idxq + ic;
 147         i__2 = nr;
 148         for (j = 1; j <= i__2; ++j) {
 149             iwork[itemp + j - 1] = j;
 150         }
 151     }
 152
 153     for (lvl = nlvl; lvl >= 1; --lvl) {
 154
 155         if (lvl == 1) {
 156             lf = 1;
 157             ll = 1;
 158         } else {
 159             i__1 = lvl - 1;
 160             lf = (1 << i__1);
 161             ll = (lf << 1) - 1;
 162         }
 163         i__1 = ll;
 164         for (i__ = lf; i__ <= i__1; ++i__) {
 165             im1 = i__ - 1;
 166             ic = iwork[inode + im1];
 167             nl = iwork[ndiml + im1];
 168             nr = iwork[ndimr + im1];
 169             nlf = ic - nl;
 170             if (*sqre == 0 && i__ == ll) {
 171                 sqrei = *sqre;
 172             } else {
 173                 sqrei = 1;
 174             }
 175             idxqc = idxq + nlf - 1;
 176             alpha = d__[ic];
 177             beta = e[ic];
 178             F77_FUNC(slasd1,SLASD1)(&nl, &nr, &sqrei, &d__[nlf], &alpha, &beta, &u[nlf + nlf *
 179                      u_dim1], ldu, &vt[nlf + nlf * vt_dim1], ldvt, &iwork[
 180                     idxqc], &iwork[iwk], &work[1], info);
 181             if (*info != 0) {
 182                 return;
 183             }
 184         }
 185     }
 186
 187     return;
 188
 189 }