src/gmxlib/gmx_lapack/dlaruv.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_lapack.h"
  37
  38 void
  39 F77_FUNC(dlaruv,DLARUV)(int *iseed, int *n, double *x)
  40 {
  41   const int
  42     mm[512] = {
  43       494,2637,255,2008,1253,
  44       3344,4084,1739,3143,3468,688,1657,1238,3166,1292,3422,1270,2016,
  45       154,2862,697,1706,491,931,1444,444,3577,3944,2184,1661,3482,657,
  46       3023,3618,1267,1828,164,3798,3087,2400,2870,3876,1905,1593,1797,
  47       1234,3460,328,2861,1950,617,2070,3331,769,1558,2412,2800,189,287,
  48       2045,1227,2838,209,2770,3654,3993,192,2253,3491,2889,2857,2094,
  49       1818,688,1407,634,3231,815,3524,1914,516,164,303,2144,3480,119,
  50       3357,837,2826,2332,2089,3780,1700,3712,150,2000,3375,1621,3090,
  51       3765,1149,3146,33,3082,2741,359,3316,1749,185,2784,2202,2199,1364,
  52       1244,2020,3160,2785,2772,1217,1822,1245,2252,3904,2774,997,2573,
  53       1148,545,322,789,1440,752,2859,123,1848,643,2405,2638,2344,46,
  54       3814,913,3649,339,3808,822,2832,3078,3633,2970,637,2249,2081,4019,
  55       1478,242,481,2075,4058,622,3376,812,234,641,4005,1122,3135,2640,
  56       2302,40,1832,2247,2034,2637,1287,1691,496,1597,2394,2584,1843,336,
  57       1472,2407,433,2096,1761,2810,566,442,41,1238,1086,603,840,3168,
  58       1499,1084,3438,2408,1589,2391,288,26,512,1456,171,1677,2657,2270,
  59       2587,2961,1970,1817,676,1410,3723,2803,3185,184,663,499,3784,1631,
  60       1925,3912,1398,1349,1441,2224,2411,1907,3192,2786,382,37,759,2948,
  61       1862,3802,2423,2051,2295,1332,1832,2405,3638,3661,327,3660,716,
  62       1842,3987,1368,1848,2366,2508,3754,1766,3572,2893,307,1297,3966,
  63       758,2598,3406,2922,1038,2934,2091,2451,1580,1958,2055,1507,1078,
  64       3273,17,854,2916,3971,2889,3831,2621,1541,893,736,3992,787,2125,
  65       2364,2460,257,1574,3912,1216,3248,3401,2124,2762,149,2245,166,466,
  66       4018,1399,190,2879,153,2320,18,712,2159,2318,2091,3443,1510,449,
  67       1956,2201,3137,3399,1321,2271,3667,2703,629,2365,2431,1113,3922,
  68       2554,184,2099,3228,4012,1921,3452,3901,572,3309,3171,817,3039,
  69       1696,1256,3715,2077,3019,1497,1101,717,51,981,1978,1813,3881,76,
  70       3846,3694,1682,124,1660,3997,479,1141,886,3514,1301,3604,1888,
  71       1836,1990,2058,692,1194,20,3285,2046,2107,3508,3525,3801,2549,
  72       1145,2253,305,3301,1065,3133,2913,3285,1241,1197,3729,2501,1673,
  73       541,2753,949,2361,1165,4081,2725,3305,3069,3617,3733,409,2157,
  74       1361,3973,1865,2525,1409,3445,3577,77,3761,2149,1449,3005,225,85,
  75       3673,3117,3089,1349,2057,413,65,1845,697,3085,3441,1573,3689,2941,
  76       929,533,2841,4077,721,2821,2249,2397,2817,245,1913,1997,3121,997,
  77       1833,2877,1633,981,2009,941,2449,197,2441,285,1473,2741,3129,909,
  78       2801,421,4073,2813,2337,1429,1177,1901,81,1669,2633,2269,129,1141,
  79       249,3917,2481,3941,2217,2749,3041,1877,345,2861,1809,3141,2825,
  80       157,2881,3637,1465,2829,2161,3365,361,2685,3745,2325,3609,3821,
  81       3537,517,3017,2141,1537
  82     };
  83
  84     int i__1;
  85
  86     int i__, i1, i2, i3, i4, it1, it2, it3, it4;
  87
  88
  89     --iseed;
  90     --x;
  91
  92     it1 = it2 = it3 = it4 = 0;
  93
  94     i1 = iseed[1];
  95     i2 = iseed[2];
  96     i3 = iseed[3];
  97     i4 = iseed[4];
  98
  99     i__1 = (*n<128) ? *n : 128;
 100     for (i__ = 1; i__ <= i__1; ++i__) {
 101
 102         it4 = i4 * mm[i__ + 383];
 103         it3 = it4 / 4096;
 104         it4 -= it3 << 12;
 105         it3 = it3 + i3 * mm[i__ + 383] + i4 * mm[i__ + 255];
 106         it2 = it3 / 4096;
 107         it3 -= it2 << 12;
 108         it2 = it2 + i2 * mm[i__ + 383] + i3 * mm[i__ + 255] +
 109           i4 * mm[i__ + 127];
 110         it1 = it2 / 4096;
 111         it2 -= it1 << 12;
 112         it1 = it1 + i1 * mm[i__ + 383] + i2 * mm[i__ + 255] +
 113           i3 * mm[i__ + 127] + i4 * mm[i__ - 1];
 114         it1 %= 4096;
 115
 116         x[i__] = ((double) it1 + ((double) it2 + ((double) it3 + (
 117                 double) it4 * 2.44140625e-4) * 2.44140625e-4) *
 118                 2.44140625e-4) * 2.44140625e-4;
 119     }
 120
 121     iseed[1] = it1;
 122     iseed[2] = it2;
 123     iseed[3] = it3;
 124     iseed[4] = it4;
 125     return;
 126
 127 }