src/gmxlib/gmx_lapack/slasda.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_blas.h"
  36 #include "gmx_lapack.h"
  37
  38 void
  39 F77_FUNC(slasda,SLASDA)(int *icompq,
  40         int *smlsiz,
  41         int *n,
  42         int *sqre,
  43         float *d__,
  44         float *e,
  45         float *u,
  46         int *ldu,
  47         float *vt,
  48         int *k,
  49         float *difl,
  50         float *difr,
  51         float *z__,
  52         float *poles,
  53         int *givptr,
  54         int *givcol,
  55         int *ldgcol,
  56         int *perm,
  57         float *givnum,
  58         float *c__,
  59         float *s,
  60         float *work,
  61         int *iwork,
  62         int *info)
  63 {
  64     int givcol_dim1, givcol_offset, perm_dim1, perm_offset, difl_dim1,
  65             difl_offset, difr_dim1, difr_offset, givnum_dim1, givnum_offset,
  66             poles_dim1, poles_offset, u_dim1, u_offset, vt_dim1, vt_offset,
  67             z_dim1, z_offset, i__1, i__2;
  68
  69     int i__, j, m, i1, ic, lf, nd, ll, nl, vf, nr, vl, im1, ncc,
  70             nlf, nrf, vfi, iwk, vli, lvl, nru, ndb1, nlp1, lvl2, nrp1;
  71     float beta;
  72     int idxq, nlvl;
  73     float alpha;
  74     int inode, ndiml, ndimr, idxqi, itemp;
  75     int sqrei;
  76     int nwork1, nwork2;
  77     int smlszp;
  78     int c__0 = 0;
  79     float zero = 0.0;
  80     float one = 1.;
  81     int c__1 = 1;
  82     --d__;
  83     --e;
  84     givnum_dim1 = *ldu;
  85     givnum_offset = 1 + givnum_dim1;
  86     givnum -= givnum_offset;
  87     poles_dim1 = *ldu;
  88     poles_offset = 1 + poles_dim1;
  89     poles -= poles_offset;
  90     z_dim1 = *ldu;
  91     z_offset = 1 + z_dim1;
  92     z__ -= z_offset;
  93     difr_dim1 = *ldu;
  94     difr_offset = 1 + difr_dim1;
  95     difr -= difr_offset;
  96     difl_dim1 = *ldu;
  97     difl_offset = 1 + difl_dim1;
  98     difl -= difl_offset;
  99     vt_dim1 = *ldu;
 100     vt_offset = 1 + vt_dim1;
 101     vt -= vt_offset;
 102     u_dim1 = *ldu;
 103     u_offset = 1 + u_dim1;
 104     u -= u_offset;
 105     --k;
 106     --givptr;
 107     perm_dim1 = *ldgcol;
 108     perm_offset = 1 + perm_dim1;
 109     perm -= perm_offset;
 110     givcol_dim1 = *ldgcol;
 111     givcol_offset = 1 + givcol_dim1;
 112     givcol -= givcol_offset;
 113     --c__;
 114     --s;
 115     --work;
 116     --iwork;
 117     *info = 0;
 118
 119     m = *n + *sqre;
 120
 121     if (*n <= *smlsiz) {
 122         if (*icompq == 0) {
 123             F77_FUNC(slasdq,SLASDQ)("U", sqre, n, &c__0, &c__0, &c__0, &d__[1], &e[1], &vt[
 124                     vt_offset], ldu, &u[u_offset], ldu, &u[u_offset], ldu, &
 125                     work[1], info);
 126         } else {
 127             F77_FUNC(slasdq,SLASDQ)("U", sqre, n, &m, n, &c__0, &d__[1], &e[1], &vt[vt_offset]
 128                     , ldu, &u[u_offset], ldu, &u[u_offset], ldu, &work[1],
 129                     info);
 130         }
 131         return;
 132     }
 133
 134     inode = 1;
 135     ndiml = inode + *n;
 136     ndimr = ndiml + *n;
 137     idxq = ndimr + *n;
 138     iwk = idxq + *n;
 139
 140     ncc = 0;
 141     nru = 0;
 142
 143     smlszp = *smlsiz + 1;
 144     vf = 1;
 145     vl = vf + m;
 146     nwork1 = vl + m;
 147     nwork2 = nwork1 + smlszp * smlszp;
 148
 149     F77_FUNC(slasdt,SLASDT)(n, &nlvl, &nd, &iwork[inode], &iwork[ndiml], &iwork[ndimr],
 150             smlsiz);
 151
 152     ndb1 = (nd + 1) / 2;
 153     i__1 = nd;
 154     for (i__ = ndb1; i__ <= i__1; ++i__) {
 155         i1 = i__ - 1;
 156         ic = iwork[inode + i1];
 157         nl = iwork[ndiml + i1];
 158         nlp1 = nl + 1;
 159         nr = iwork[ndimr + i1];
 160         nlf = ic - nl;
 161         nrf = ic + 1;
 162         idxqi = idxq + nlf - 2;
 163         vfi = vf + nlf - 1;
 164         vli = vl + nlf - 1;
 165         sqrei = 1;
 166         if (*icompq == 0) {
 167             F77_FUNC(slaset,SLASET)("A", &nlp1, &nlp1, &zero, &one, &work[nwork1], &smlszp);
 168             F77_FUNC(slasdq,SLASDQ)("U", &sqrei, &nl, &nlp1, &nru, &ncc, &d__[nlf], &e[nlf], &
 169                     work[nwork1], &smlszp, &work[nwork2], &nl, &work[nwork2],
 170                     &nl, &work[nwork2], info);
 171             itemp = nwork1 + nl * smlszp;
 172             F77_FUNC(scopy,SCOPY)(&nlp1, &work[nwork1], &c__1, &work[vfi], &c__1);
 173             F77_FUNC(scopy,SCOPY)(&nlp1, &work[itemp], &c__1, &work[vli], &c__1);
 174         } else {
 175             F77_FUNC(slaset,SLASET)("A", &nl, &nl, &zero, &one, &u[nlf + u_dim1], ldu);
 176             F77_FUNC(slaset,SLASET)("A", &nlp1, &nlp1, &zero, &one, &vt[nlf + vt_dim1],
 177                     ldu);
 178             F77_FUNC(slasdq,SLASDQ)("U", &sqrei, &nl, &nlp1, &nl, &ncc, &d__[nlf], &e[nlf], &
 179                     vt[nlf + vt_dim1], ldu, &u[nlf + u_dim1], ldu, &u[nlf +
 180                     u_dim1], ldu, &work[nwork1], info);
 181             F77_FUNC(scopy,SCOPY)(&nlp1, &vt[nlf + vt_dim1], &c__1, &work[vfi], &c__1);
 182             F77_FUNC(scopy,SCOPY)(&nlp1, &vt[nlf + nlp1 * vt_dim1], &c__1, &work[vli], &c__1)
 183                     ;
 184         }
 185         if (*info != 0) {
 186             return;
 187         }
 188         i__2 = nl;
 189         for (j = 1; j <= i__2; ++j) {
 190             iwork[idxqi + j] = j;
 191         }
 192         if (i__ == nd && *sqre == 0) {
 193             sqrei = 0;
 194         } else {
 195             sqrei = 1;
 196         }
 197         idxqi += nlp1;
 198         vfi += nlp1;
 199         vli += nlp1;
 200         nrp1 = nr + sqrei;
 201         if (*icompq == 0) {
 202             F77_FUNC(slaset,SLASET)("A", &nrp1, &nrp1, &zero, &one, &work[nwork1], &smlszp);
 203             F77_FUNC(slasdq,SLASDQ)("U", &sqrei, &nr, &nrp1, &nru, &ncc, &d__[nrf], &e[nrf], &
 204                     work[nwork1], &smlszp, &work[nwork2], &nr, &work[nwork2],
 205                     &nr, &work[nwork2], info);
 206             itemp = nwork1 + (nrp1 - 1) * smlszp;
 207             F77_FUNC(scopy,SCOPY)(&nrp1, &work[nwork1], &c__1, &work[vfi], &c__1);
 208             F77_FUNC(scopy,SCOPY)(&nrp1, &work[itemp], &c__1, &work[vli], &c__1);
 209         } else {
 210             F77_FUNC(slaset,SLASET)("A", &nr, &nr, &zero, &one, &u[nrf + u_dim1], ldu);
 211             F77_FUNC(slaset,SLASET)("A", &nrp1, &nrp1, &zero, &one, &vt[nrf + vt_dim1],
 212                     ldu);
 213             F77_FUNC(slasdq,SLASDQ)("U", &sqrei, &nr, &nrp1, &nr, &ncc, &d__[nrf], &e[nrf], &
 214                     vt[nrf + vt_dim1], ldu, &u[nrf + u_dim1], ldu, &u[nrf +
 215                     u_dim1], ldu, &work[nwork1], info);
 216             F77_FUNC(scopy,SCOPY)(&nrp1, &vt[nrf + vt_dim1], &c__1, &work[vfi], &c__1);
 217             F77_FUNC(scopy,SCOPY)(&nrp1, &vt[nrf + nrp1 * vt_dim1], &c__1, &work[vli], &c__1)
 218                     ;
 219         }
 220         if (*info != 0) {
 221             return;
 222         }
 223         i__2 = nr;
 224         for (j = 1; j <= i__2; ++j) {
 225             iwork[idxqi + j] = j;
 226         }
 227     }
 228
 229     j = (1 << nlvl);
 230
 231     for (lvl = nlvl; lvl >= 1; --lvl) {
 232         lvl2 = (lvl << 1) - 1;
 233
 234         if (lvl == 1) {
 235             lf = 1;
 236             ll = 1;
 237         } else {
 238             i__1 = lvl - 1;
 239             lf = (1 << (lvl-1));
 240             ll = (lf << 1) - 1;
 241         }
 242         i__1 = ll;
 243         for (i__ = lf; i__ <= i__1; ++i__) {
 244             im1 = i__ - 1;
 245             ic = iwork[inode + im1];
 246             nl = iwork[ndiml + im1];
 247             nr = iwork[ndimr + im1];
 248             nlf = ic - nl;
 249             nrf = ic + 1;
 250             if (i__ == ll) {
 251                 sqrei = *sqre;
 252             } else {
 253                 sqrei = 1;
 254             }
 255             vfi = vf + nlf - 1;
 256             vli = vl + nlf - 1;
 257             idxqi = idxq + nlf - 1;
 258             alpha = d__[ic];
 259             beta = e[ic];
 260             if (*icompq == 0) {
 261                 F77_FUNC(slasd6,SLASD6)(icompq, &nl, &nr, &sqrei, &d__[nlf], &work[vfi], &
 262                         work[vli], &alpha, &beta, &iwork[idxqi], &perm[
 263                         perm_offset], &givptr[1], &givcol[givcol_offset],
 264                         ldgcol, &givnum[givnum_offset], ldu, &poles[
 265                         poles_offset], &difl[difl_offset], &difr[difr_offset],
 266                          &z__[z_offset], &k[1], &c__[1], &s[1], &work[nwork1],
 267                          &iwork[iwk], info);
 268             } else {
 269                 --j;
 270                 F77_FUNC(slasd6,SLASD6)(icompq, &nl, &nr, &sqrei, &d__[nlf], &work[vfi], &
 271                         work[vli], &alpha, &beta, &iwork[idxqi], &perm[nlf +
 272                         lvl * perm_dim1], &givptr[j], &givcol[nlf + lvl2 *
 273                         givcol_dim1], ldgcol, &givnum[nlf + lvl2 *
 274                         givnum_dim1], ldu, &poles[nlf + lvl2 * poles_dim1], &
 275                         difl[nlf + lvl * difl_dim1], &difr[nlf + lvl2 *
 276                         difr_dim1], &z__[nlf + lvl * z_dim1], &k[j], &c__[j],
 277                         &s[j], &work[nwork1], &iwork[iwk], info);
 278             }
 279             if (*info != 0) {
 280                 return;
 281             }
 282         }
 283     }
 284
 285     return;
 286
 287 }
 288
 289