src/gromacs/linearalgebra/gmx_lapack/dstein.c

   1 #include <math.h>
   2 #include "../gmx_blas.h"
   3 #include "../gmx_lapack.h"
   4 #include "lapack_limits.h"
   5
   6 #include "gromacs/utility/real.h"
   7
   8 void
   9 F77_FUNC(dstein,DSTEIN)(int *n,
  10         double *d__,
  11         double *e,
  12         int *m,
  13         double *w,
  14         int *iblock,
  15         int *isplit,
  16         double *z__,
  17         int *ldz,
  18         double *work,
  19         int *iwork,
  20         int *ifail,
  21         int *info)
  22 {
  23     int z_dim1, z_offset, i__1, i__2, i__3;
  24     double d__2, d__3, d__4, d__5;
  25
  26     int i__, j, b1, j1, bn;
  27     double xj, scl, eps, sep, nrm, tol;
  28     int its;
  29     double xjm, ztr, eps1;
  30     int jblk, nblk;
  31     int jmax;
  32
  33     int iseed[4], gpind, iinfo;
  34     double ortol;
  35     int indrv1, indrv2, indrv3, indrv4, indrv5;
  36     int nrmchk;
  37     int blksiz;
  38     double onenrm, dtpcrt, pertol;
  39     int c__2 = 2;
  40     int c__1 = 1;
  41     int c_n1 = -1;
  42
  43     --d__;
  44     --e;
  45     --w;
  46     --iblock;
  47     --isplit;
  48     z_dim1 = *ldz;
  49     z_offset = 1 + z_dim1;
  50     z__ -= z_offset;
  51     --work;
  52     --iwork;
  53     --ifail;
  54
  55     *info = 0;
  56
  57     xjm = 0.0;
  58     i__1 = *m;
  59     for (i__ = 1; i__ <= i__1; ++i__) {
  60         ifail[i__] = 0;
  61     }
  62
  63     if (*n < 0) {
  64         *info = -1;
  65     } else if (*m < 0 || *m > *n) {
  66         *info = -4;
  67     } else if (*ldz < (*n)) {
  68         *info = -9;
  69     } else {
  70         i__1 = *m;
  71         for (j = 2; j <= i__1; ++j) {
  72             if (iblock[j] < iblock[j - 1]) {
  73                 *info = -6;
  74                 break;
  75             }
  76             if (iblock[j] == iblock[j - 1] && w[j] < w[j - 1]) {
  77                 *info = -5;
  78                 break;
  79             }
  80         }
  81     }
  82
  83     if (*info != 0) {
  84         i__1 = -(*info);
  85         return;
  86     }
  87
  88     if (*n == 0 || *m == 0) {
  89         return;
  90     } else if (*n == 1) {
  91         z__[z_dim1 + 1] = 1.;
  92         return;
  93     }
  94
  95     eps = GMX_DOUBLE_EPS;
  96
  97     for (i__ = 1; i__ <= 4; ++i__) {
  98         iseed[i__ - 1] = 1;
  99     }
 100
 101     indrv1 = 0;
 102     indrv2 = indrv1 + *n;
 103     indrv3 = indrv2 + *n;
 104     indrv4 = indrv3 + *n;
 105     indrv5 = indrv4 + *n;
 106
 107     j1 = 1;
 108     i__1 = iblock[*m];
 109     for (nblk = 1; nblk <= i__1; ++nblk) {
 110
 111         if (nblk == 1) {
 112             b1 = 1;
 113         } else {
 114             b1 = isplit[nblk - 1] + 1;
 115         }
 116         bn = isplit[nblk];
 117         blksiz = bn - b1 + 1;
 118         if (blksiz == 1) {
 119             continue;
 120         }
 121         gpind = b1;
 122
 123         onenrm = fabs(d__[b1]) + fabs(e[b1]);
 124         d__3 = onenrm;
 125         d__4 = fabs(d__[bn]) + fabs(e[bn - 1]);
 126         onenrm = (d__3>d__4) ? d__3 : d__4;
 127         i__2 = bn - 1;
 128         for (i__ = b1 + 1; i__ <= i__2; ++i__) {
 129           d__4 = onenrm;
 130           d__5 = fabs(d__[i__]) + fabs(e[i__ - 1]) + fabs(e[i__]);
 131             onenrm = (d__4>d__5) ? d__4 : d__5;
 132         }
 133         ortol = onenrm * .001;
 134
 135         dtpcrt = sqrt(.1 / blksiz);
 136
 137         jblk = 0;
 138         i__2 = *m;
 139         for (j = j1; j <= i__2; ++j) {
 140             if (iblock[j] != nblk) {
 141                 j1 = j;
 142                 break;
 143             }
 144             ++jblk;
 145             xj = w[j];
 146
 147             if (blksiz == 1) {
 148                 work[indrv1 + 1] = 1.;
 149                 goto L120;
 150             }
 151
 152             if (jblk > 1) {
 153                 eps1 = fabs(eps * xj);
 154                 pertol = eps1 * 10.;
 155                 sep = xj - xjm;
 156                 if (sep < pertol) {
 157                     xj = xjm + pertol;
 158                 }
 159             }
 160
 161             its = 0;
 162             nrmchk = 0;
 163
 164             F77_FUNC(dlarnv,DLARNV)(&c__2, iseed, &blksiz, &work[indrv1 + 1]);
 165
 166             F77_FUNC(dcopy,DCOPY)(&blksiz, &d__[b1], &c__1, &work[indrv4 + 1], &c__1);
 167             i__3 = blksiz - 1;
 168             F77_FUNC(dcopy,DCOPY)(&i__3, &e[b1], &c__1, &work[indrv2 + 2], &c__1);
 169             i__3 = blksiz - 1;
 170             F77_FUNC(dcopy,DCOPY)(&i__3, &e[b1], &c__1, &work[indrv3 + 1], &c__1);
 171
 172             tol = 0.;
 173             F77_FUNC(dlagtf,DLAGTF)(&blksiz, &work[indrv4 + 1], &xj, &work[indrv2 + 2], &work[
 174                     indrv3 + 1], &tol, &work[indrv5 + 1], &iwork[1], &iinfo);
 175
 176 L70:
 177             ++its;
 178             if (its > 5) {
 179                 goto L100;
 180             }
 181
 182             d__2 = eps;
 183             d__3 = fabs(work[indrv4 + blksiz]);
 184             scl = blksiz * onenrm * ((d__2>d__3) ? d__2 : d__3) / F77_FUNC(dasum,DASUM)(&blksiz, &work[
 185                     indrv1 + 1], &c__1);
 186             F77_FUNC(dscal,DSCAL)(&blksiz, &scl, &work[indrv1 + 1], &c__1);
 187
 188             F77_FUNC(dlagts,DLAGTS)(&c_n1, &blksiz, &work[indrv4 + 1], &work[indrv2 + 2], &
 189                     work[indrv3 + 1], &work[indrv5 + 1], &iwork[1], &work[
 190                     indrv1 + 1], &tol, &iinfo);
 191
 192             if (jblk == 1) {
 193                 goto L90;
 194             }
 195             if (fabs(xj - xjm) > ortol) {
 196                 gpind = j;
 197             }
 198             if (gpind != j) {
 199                 i__3 = j - 1;
 200                 for (i__ = gpind; i__ <= i__3; ++i__) {
 201                     ztr = -F77_FUNC(ddot,DDOT)(&blksiz, &work[indrv1 + 1], &c__1, &z__[b1 +
 202                             i__ * z_dim1], &c__1);
 203                     F77_FUNC(daxpy,DAXPY)(&blksiz, &ztr, &z__[b1 + i__ * z_dim1], &c__1, &
 204                             work[indrv1 + 1], &c__1);
 205                 }
 206             }
 207
 208 L90:
 209             jmax = F77_FUNC(idamax,IDAMAX)(&blksiz, &work[indrv1 + 1], &c__1);
 210             nrm = fabs(work[indrv1 + jmax]);
 211
 212             if (nrm < dtpcrt) {
 213                 goto L70;
 214             }
 215             ++nrmchk;
 216             if (nrmchk < 3) {
 217                 goto L70;
 218             }
 219
 220             goto L110;
 221
 222 L100:
 223             ++(*info);
 224             ifail[*info] = j;
 225
 226 L110:
 227             scl = 1. / F77_FUNC(dnrm2,DNRM2)(&blksiz, &work[indrv1 + 1], &c__1);
 228             jmax = F77_FUNC(idamax,IDAMAX)(&blksiz, &work[indrv1 + 1], &c__1);
 229             if (work[indrv1 + jmax] < 0.) {
 230                 scl = -scl;
 231             }
 232             F77_FUNC(dscal,DSCAL)(&blksiz, &scl, &work[indrv1 + 1], &c__1);
 233 L120:
 234             i__3 = *n;
 235             for (i__ = 1; i__ <= i__3; ++i__) {
 236                 z__[i__ + j * z_dim1] = 0.;
 237             }
 238             i__3 = blksiz;
 239             for (i__ = 1; i__ <= i__3; ++i__) {
 240                 z__[b1 + i__ - 1 + j * z_dim1] = work[indrv1 + i__];
 241             }
 242
 243             xjm = xj;
 244         }
 245     }
 246
 247     return;
 248
 249 }
 250
 251