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