src/gromacs/linearalgebra/gmx_lapack/dlasd2.c

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