src/gmxlib/gmx_lapack/slasd7.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 <types/simple.h>
  37
  38 #include "gmx_blas.h"
  39 #include "gmx_lapack.h"
  40 #include "lapack_limits.h"
  41
  42 void
  43 F77_FUNC(slasd7,SLASD7)(int *icompq,
  44         int *nl,
  45         int *nr,
  46         int *sqre,
  47         int *k,
  48         float *d__,
  49         float *z__,
  50         float *zw,
  51         float *vf,
  52         float *vfw,
  53         float *vl,
  54         float *vlw,
  55         float *alpha,
  56         float *beta,
  57         float *dsigma,
  58         int *idx,
  59         int *idxp,
  60         int *idxq,
  61         int *perm,
  62         int *givptr,
  63         int *givcol,
  64         int *ldgcol,
  65         float *givnum,
  66         int *ldgnum,
  67         float *c__,
  68         float *s,
  69         int *info)
  70 {
  71     int givcol_dim1, givcol_offset, givnum_dim1, givnum_offset, i__1;
  72     float d__1, d__2;
  73
  74     int i__, j, m, n, k2;
  75     float z1;
  76     int jp;
  77     float eps, tau, tol;
  78     int nlp1, nlp2, idxi, idxj;
  79     int idxjp;
  80     int jprev = 0;
  81     float hlftol;
  82     int c__1 = 1;
  83
  84     --d__;
  85     --z__;
  86     --zw;
  87     --vf;
  88     --vfw;
  89     --vl;
  90     --vlw;
  91     --dsigma;
  92     --idx;
  93     --idxp;
  94     --idxq;
  95     --perm;
  96     givcol_dim1 = *ldgcol;
  97     givcol_offset = 1 + givcol_dim1;
  98     givcol -= givcol_offset;
  99     givnum_dim1 = *ldgnum;
 100     givnum_offset = 1 + givnum_dim1;
 101     givnum -= givnum_offset;
 102
 103     *info = 0;
 104     n = *nl + *nr + 1;
 105     m = n + *sqre;
 106
 107     nlp1 = *nl + 1;
 108     nlp2 = *nl + 2;
 109     if (*icompq == 1) {
 110         *givptr = 0;
 111     }
 112
 113     z1 = *alpha * vl[nlp1];
 114     vl[nlp1] = 0.;
 115     tau = vf[nlp1];
 116     for (i__ = *nl; i__ >= 1; --i__) {
 117         z__[i__ + 1] = *alpha * vl[i__];
 118         vl[i__] = 0.;
 119         vf[i__ + 1] = vf[i__];
 120         d__[i__ + 1] = d__[i__];
 121         idxq[i__ + 1] = idxq[i__] + 1;
 122     }
 123     vf[1] = tau;
 124
 125     i__1 = m;
 126     for (i__ = nlp2; i__ <= i__1; ++i__) {
 127         z__[i__] = *beta * vf[i__];
 128         vf[i__] = 0.;
 129     }
 130     i__1 = n;
 131     for (i__ = nlp2; i__ <= i__1; ++i__) {
 132         idxq[i__] += nlp1;
 133     }
 134
 135     i__1 = n;
 136     for (i__ = 2; i__ <= i__1; ++i__) {
 137         dsigma[i__] = d__[idxq[i__]];
 138         zw[i__] = z__[idxq[i__]];
 139         vfw[i__] = vf[idxq[i__]];
 140         vlw[i__] = vl[idxq[i__]];
 141     }
 142
 143     F77_FUNC(slamrg,SLAMRG)(nl, nr, &dsigma[2], &c__1, &c__1, &idx[2]);
 144
 145     i__1 = n;
 146     for (i__ = 2; i__ <= i__1; ++i__) {
 147         idxi = idx[i__] + 1;
 148         d__[i__] = dsigma[idxi];
 149         z__[i__] = zw[idxi];
 150         vf[i__] = vfw[idxi];
 151         vl[i__] = vlw[idxi];
 152     }
 153
 154     eps = GMX_FLOAT_EPS;
 155
 156     d__1 = fabs(*alpha);
 157     d__2 = fabs(*beta);
 158     tol = (d__1>d__2) ? d__1 : d__2;
 159     d__2 = fabs(d__[n]);
 160     tol = eps * 64. * ((d__2>tol) ? d__2 : tol);
 161
 162     *k = 1;
 163     k2 = n + 1;
 164     i__1 = n;
 165     for (j = 2; j <= i__1; ++j) {
 166         if (fabs(z__[j]) <= tol) {
 167
 168             --k2;
 169             idxp[k2] = j;
 170             if (j == n) {
 171                 goto L100;
 172             }
 173         } else {
 174             jprev = j;
 175             goto L70;
 176         }
 177     }
 178 L70:
 179     j = jprev;
 180 L80:
 181     ++j;
 182     if (j > n) {
 183         goto L90;
 184     }
 185     if (fabs(z__[j]) <= tol) {
 186
 187         --k2;
 188         idxp[k2] = j;
 189     } else {
 190
 191         if (fabs(d__[j] - d__[jprev]) <= tol) {
 192
 193             *s = z__[jprev];
 194             *c__ = z__[j];
 195
 196             tau = F77_FUNC(slapy2,SLAPY2)(c__, s);
 197             z__[j] = tau;
 198             z__[jprev] = 0.;
 199             *c__ /= tau;
 200             *s = -(*s) / tau;
 201
 202
 203             if (*icompq == 1) {
 204                 ++(*givptr);
 205                 idxjp = idxq[idx[jprev] + 1];
 206                 idxj = idxq[idx[j] + 1];
 207                 if (idxjp <= nlp1) {
 208                     --idxjp;
 209                 }
 210                 if (idxj <= nlp1) {
 211                     --idxj;
 212                 }
 213                 givcol[*givptr + (givcol_dim1 << 1)] = idxjp;
 214                 givcol[*givptr + givcol_dim1] = idxj;
 215                 givnum[*givptr + (givnum_dim1 << 1)] = *c__;
 216                 givnum[*givptr + givnum_dim1] = *s;
 217             }
 218             F77_FUNC(srot,SROT)(&c__1, &vf[jprev], &c__1, &vf[j], &c__1, c__, s);
 219             F77_FUNC(srot,SROT)(&c__1, &vl[jprev], &c__1, &vl[j], &c__1, c__, s);
 220             --k2;
 221             idxp[k2] = jprev;
 222             jprev = j;
 223         } else {
 224             ++(*k);
 225             zw[*k] = z__[jprev];
 226             dsigma[*k] = d__[jprev];
 227             idxp[*k] = jprev;
 228             jprev = j;
 229         }
 230     }
 231     goto L80;
 232 L90:
 233
 234     ++(*k);
 235     zw[*k] = z__[jprev];
 236     dsigma[*k] = d__[jprev];
 237     idxp[*k] = jprev;
 238
 239 L100:
 240
 241     i__1 = n;
 242     for (j = 2; j <= i__1; ++j) {
 243         jp = idxp[j];
 244         dsigma[j] = d__[jp];
 245         vfw[j] = vf[jp];
 246         vlw[j] = vl[jp];
 247     }
 248     if (*icompq == 1) {
 249         i__1 = n;
 250         for (j = 2; j <= i__1; ++j) {
 251             jp = idxp[j];
 252             perm[j] = idxq[idx[jp] + 1];
 253             if (perm[j] <= nlp1) {
 254                 --perm[j];
 255             }
 256         }
 257     }
 258     i__1 = n - *k;
 259     F77_FUNC(scopy,SCOPY)(&i__1, &dsigma[*k + 1], &c__1, &d__[*k + 1], &c__1);
 260
 261     dsigma[1] = 0.;
 262     hlftol = tol / 2.;
 263     if (fabs(dsigma[2]) <= hlftol) {
 264         dsigma[2] = hlftol;
 265     }
 266     if (m > n) {
 267         z__[1] = F77_FUNC(slapy2,SLAPY2)(&z1, &z__[m]);
 268         if (z__[1] <= tol) {
 269             *c__ = 1.;
 270             *s = 0.;
 271             z__[1] = tol;
 272         } else {
 273             *c__ = z1 / z__[1];
 274             *s = -z__[m] / z__[1];
 275         }
 276         F77_FUNC(srot,SROT)(&c__1, &vf[m], &c__1, &vf[1], &c__1, c__, s);
 277         F77_FUNC(srot,SROT)(&c__1, &vl[m], &c__1, &vl[1], &c__1, c__, s);
 278     } else {
 279         if (fabs(z1) <= tol) {
 280             z__[1] = tol;
 281         } else {
 282             z__[1] = z1;
 283         }
 284     }
 285
 286     i__1 = *k - 1;
 287     F77_FUNC(scopy,SCOPY)(&i__1, &zw[2], &c__1, &z__[2], &c__1);
 288     i__1 = n - 1;
 289     F77_FUNC(scopy,SCOPY)(&i__1, &vfw[2], &c__1, &vf[2], &c__1);
 290     i__1 = n - 1;
 291     F77_FUNC(scopy,SCOPY)(&i__1, &vlw[2], &c__1, &vl[2], &c__1);
 292
 293     return;
 294
 295 }
 296
 297