src/gmxlib/gmx_lapack/dlasd2.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 <math.h>
  36 #include "gmx_blas.h"
  37 #include "gmx_lapack.h"
  38 #include "lapack_limits.h"
  39
  40 #include <types/simple.h>
  41
  42 void
  43 F77_FUNC(dlasd2,DLASD2)(int *nl,
  44                         int *nr,
  45                         int *sqre,
  46                         int *k,
  47                         double *d__,
  48                         double *z__,
  49                         double *alpha,
  50                         double *beta,
  51                         double *u,
  52                         int *ldu,
  53                         double *vt,
  54                         int *ldvt,
  55                         double *dsigma,
  56                         double *u2,
  57                         int *ldu2,
  58                         double *vt2,
  59                         int *ldvt2,
  60                         int *idxp,
  61                         int *idx,
  62                         int *idxc,
  63                         int *idxq,
  64                         int *coltyp,
  65                         int *info)
  66 {
  67     int u_dim1, u_offset, u2_dim1, u2_offset, vt_dim1, vt_offset;
  68     int vt2_dim1, vt2_offset, i__1;
  69     double d__1, d__2;
  70
  71     double c__;
  72     int i__, j, m, n;
  73     double s;
  74     int k2;
  75     double z1;
  76     int ct, jp;
  77     double eps, tau, tol;
  78     int psm[4], nlp1, nlp2, idxi, idxj;
  79     int ctot[4], idxjp;
  80     int jprev = 0;
  81     double hlftol;
  82     double zero = 0.0;
  83     int c__1 = 1;
  84
  85
  86     --d__;
  87     --z__;
  88     u_dim1 = *ldu;
  89     u_offset = 1 + u_dim1;
  90     u -= u_offset;
  91     vt_dim1 = *ldvt;
  92     vt_offset = 1 + vt_dim1;
  93     vt -= vt_offset;
  94     --dsigma;
  95     u2_dim1 = *ldu2;
  96     u2_offset = 1 + u2_dim1;
  97     u2 -= u2_offset;
  98     vt2_dim1 = *ldvt2;
  99     vt2_offset = 1 + vt2_dim1;
 100     vt2 -= vt2_offset;
 101     --idxp;
 102     --idx;
 103     --idxc;
 104     --idxq;
 105     --coltyp;
 106
 107     *info = 0;
 108
 109     n = *nl + *nr + 1;
 110     m = n + *sqre;
 111
 112     nlp1 = *nl + 1;
 113     nlp2 = *nl + 2;
 114
 115     z1 = *alpha * vt[nlp1 + nlp1 * vt_dim1];
 116     z__[1] = z1;
 117     for (i__ = *nl; i__ >= 1; --i__) {
 118         z__[i__ + 1] = *alpha * vt[i__ + nlp1 * vt_dim1];
 119         d__[i__ + 1] = d__[i__];
 120         idxq[i__ + 1] = idxq[i__] + 1;
 121     }
 122
 123     i__1 = m;
 124     for (i__ = nlp2; i__ <= i__1; ++i__) {
 125         z__[i__] = *beta * vt[i__ + nlp2 * vt_dim1];
 126     }
 127
 128     i__1 = nlp1;
 129     for (i__ = 2; i__ <= i__1; ++i__) {
 130         coltyp[i__] = 1;
 131     }
 132     i__1 = n;
 133     for (i__ = nlp2; i__ <= i__1; ++i__) {
 134         coltyp[i__] = 2;
 135     }
 136
 137     i__1 = n;
 138     for (i__ = nlp2; i__ <= i__1; ++i__) {
 139         idxq[i__] += nlp1;
 140     }
 141
 142     i__1 = n;
 143     for (i__ = 2; i__ <= i__1; ++i__) {
 144         dsigma[i__] = d__[idxq[i__]];
 145         u2[i__ + u2_dim1] = z__[idxq[i__]];
 146         idxc[i__] = coltyp[idxq[i__]];
 147     }
 148
 149     F77_FUNC(dlamrg,DLAMRG)(nl, nr, &dsigma[2], &c__1, &c__1, &idx[2]);
 150
 151     i__1 = n;
 152     for (i__ = 2; i__ <= i__1; ++i__) {
 153         idxi = idx[i__] + 1;
 154         d__[i__] = dsigma[idxi];
 155         z__[i__] = u2[idxi + u2_dim1];
 156         coltyp[i__] = idxc[idxi];
 157     }
 158
 159     eps = GMX_DOUBLE_EPS;
 160     d__1 = fabs(*alpha), d__2 = fabs(*beta);
 161     tol = (d__1 > d__2) ? d__1 : d__2;
 162     d__2 = fabs(d__[n]);
 163     tol = eps * 8. * ((d__2 > tol) ? d__2 : tol);
 164
 165     *k = 1;
 166     k2 = n + 1;
 167     i__1 = n;
 168     for (j = 2; j <= i__1; ++j) {
 169         if (fabs(z__[j]) <= tol) {
 170
 171             --k2;
 172             idxp[k2] = j;
 173             coltyp[j] = 4;
 174             if (j == n) {
 175                 goto L120;
 176             }
 177         } else {
 178             jprev = j;
 179             goto L90;
 180         }
 181     }
 182 L90:
 183     j = jprev;
 184 L100:
 185     ++j;
 186     if (j > n) {
 187         goto L110;
 188     }
 189     if (fabs(z__[j]) <= tol) {
 190
 191         --k2;
 192         idxp[k2] = j;
 193         coltyp[j] = 4;
 194     } else {
 195
 196         if (fabs(d__[j] - d__[jprev]) <= tol) {
 197
 198             s = z__[jprev];
 199             c__ = z__[j];
 200
 201             tau = F77_FUNC(dlapy2,DLAPY2)(&c__, &s);
 202             c__ /= tau;
 203             s = -s / tau;
 204             z__[j] = tau;
 205             z__[jprev] = 0.;
 206
 207             idxjp = idxq[idx[jprev] + 1];
 208             idxj = idxq[idx[j] + 1];
 209             if (idxjp <= nlp1) {
 210                 --idxjp;
 211             }
 212             if (idxj <= nlp1) {
 213                 --idxj;
 214             }
 215             F77_FUNC(drot,DROT)(&n, &u[idxjp * u_dim1 + 1], &c__1, &u[idxj * u_dim1 + 1], &
 216                     c__1, &c__, &s);
 217             F77_FUNC(drot,DROT)(&m, &vt[idxjp + vt_dim1], ldvt, &vt[idxj + vt_dim1], ldvt, &
 218                     c__, &s);
 219             if (coltyp[j] != coltyp[jprev]) {
 220                 coltyp[j] = 3;
 221             }
 222             coltyp[jprev] = 4;
 223             --k2;
 224             idxp[k2] = jprev;
 225             jprev = j;
 226         } else {
 227             ++(*k);
 228             u2[*k + u2_dim1] = z__[jprev];
 229             dsigma[*k] = d__[jprev];
 230             idxp[*k] = jprev;
 231             jprev = j;
 232         }
 233     }
 234     goto L100;
 235 L110:
 236
 237     ++(*k);
 238     u2[*k + u2_dim1] = z__[jprev];
 239     dsigma[*k] = d__[jprev];
 240     idxp[*k] = jprev;
 241
 242 L120:
 243
 244     for (j = 1; j <= 4; ++j) {
 245         ctot[j - 1] = 0;
 246     }
 247     i__1 = n;
 248     for (j = 2; j <= i__1; ++j) {
 249         ct = coltyp[j];
 250         ++ctot[ct - 1];
 251     }
 252
 253     psm[0] = 2;
 254     psm[1] = ctot[0] + 2;
 255     psm[2] = psm[1] + ctot[1];
 256     psm[3] = psm[2] + ctot[2];
 257
 258     i__1 = n;
 259     for (j = 2; j <= i__1; ++j) {
 260         jp = idxp[j];
 261         ct = coltyp[jp];
 262         idxc[psm[ct - 1]] = j;
 263         ++psm[ct - 1];
 264     }
 265
 266     i__1 = n;
 267     for (j = 2; j <= i__1; ++j) {
 268         jp = idxp[j];
 269         dsigma[j] = d__[jp];
 270         idxj = idxq[idx[idxp[idxc[j]]] + 1];
 271         if (idxj <= nlp1) {
 272             --idxj;
 273         }
 274         F77_FUNC(dcopy,DCOPY)(&n, &u[idxj * u_dim1 + 1], &c__1, &u2[j * u2_dim1 + 1], &c__1);
 275         F77_FUNC(dcopy,DCOPY)(&m, &vt[idxj + vt_dim1], ldvt, &vt2[j + vt2_dim1], ldvt2);
 276     }
 277
 278     dsigma[1] = 0.;
 279     hlftol = tol / 2.;
 280     if (fabs(dsigma[2]) <= hlftol) {
 281         dsigma[2] = hlftol;
 282     }
 283     if (m > n) {
 284         z__[1] = F77_FUNC(dlapy2,DLAPY2)(&z1, &z__[m]);
 285         if (z__[1] <= tol) {
 286             c__ = 1.;
 287             s = 0.;
 288             z__[1] = tol;
 289         } else {
 290             c__ = z1 / z__[1];
 291             s = z__[m] / z__[1];
 292         }
 293     } else {
 294         if (fabs(z1) <= tol) {
 295             z__[1] = tol;
 296         } else {
 297             z__[1] = z1;
 298         }
 299     }
 300
 301     i__1 = *k - 1;
 302     F77_FUNC(dcopy,DCOPY)(&i__1, &u2[u2_dim1 + 2], &c__1, &z__[2], &c__1);
 303
 304     F77_FUNC(dlaset,DLASET)("A", &n, &c__1, &zero, &zero, &u2[u2_offset], ldu2);
 305     u2[nlp1 + u2_dim1] = 1.;
 306     if (m > n) {
 307         i__1 = nlp1;
 308         for (i__ = 1; i__ <= i__1; ++i__) {
 309             vt[m + i__ * vt_dim1] = -s * vt[nlp1 + i__ * vt_dim1];
 310             vt2[i__ * vt2_dim1 + 1] = c__ * vt[nlp1 + i__ * vt_dim1];
 311         }
 312         i__1 = m;
 313         for (i__ = nlp2; i__ <= i__1; ++i__) {
 314             vt2[i__ * vt2_dim1 + 1] = s * vt[m + i__ * vt_dim1];
 315             vt[m + i__ * vt_dim1] = c__ * vt[m + i__ * vt_dim1];
 316         }
 317     } else {
 318         F77_FUNC(dcopy,DCOPY)(&m, &vt[nlp1 + vt_dim1], ldvt, &vt2[vt2_dim1 + 1], ldvt2);
 319     }
 320     if (m > n) {
 321         F77_FUNC(dcopy,DCOPY)(&m, &vt[m + vt_dim1], ldvt, &vt2[m + vt2_dim1], ldvt2);
 322     }
 323
 324     if (n > *k) {
 325         i__1 = n - *k;
 326         F77_FUNC(dcopy,DCOPY)(&i__1, &dsigma[*k + 1], &c__1, &d__[*k + 1], &c__1);
 327         i__1 = n - *k;
 328         F77_FUNC(dlacpy,DLACPY)("A", &n, &i__1, &u2[(*k + 1) * u2_dim1 + 1], ldu2, &u[(*k + 1)
 329                  * u_dim1 + 1], ldu);
 330         i__1 = n - *k;
 331         F77_FUNC(dlacpy,DLACPY)("A", &i__1, &m, &vt2[*k + 1 + vt2_dim1], ldvt2, &vt[*k + 1 +
 332                 vt_dim1], ldvt);
 333     }
 334     for (j = 1; j <= 4; ++j) {
 335         coltyp[j] = ctot[j - 1];
 336     }
 337
 338     return;
 339
 340 }
 341
 342