src/gromacs/linearalgebra/gmx_lapack/slasd7.c

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