Module_Random¶
Copyright 2020 IRD
This file is part of statpack.
statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
statpack is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.
THIS MODULE REPLACES THE FORTRAN 90 INTRINSICS random_number AND random_seed BY SEVERAL IMPLEMENTATIONS OF THE KISS (Keep It Simple Stupid), L’Ecuyer’s LFSR113, MERSENNE TWISTER MT19937 AND MEMT19937-II RANDOM NUMBER GENERATORS.
IN ADDITION TO 10 DIFFERENT UNIFORM RANDOM GENERATORS, GAUSSIAN RANDOM GENERATORS, SHUFFLING AND SAMPLING SUBROUTINES ARE ALSO PROVIDED, AS WELL AS SUBROUTINES FOR GENERATING PSEUDO-RANDOM ORTHOGONAL MATRICES FOLLOWING THE HAAR DISTRIBUTION OVER THE GROUP OF ORTHOGONAL MATRICES, PSEUDO_RANDOM SYMMETRIC MATRICES WITH A PRESCRIBED SPECTRUM OR PSEUDO_RANDOM MATRICES WITH A PRESCRIBED SINGULAR VALUE DISTRIBUTION.
MANY PARTS OF THIS MODULE ARE ADAPTED FROM :
- Hennecke, M., 1995: A Fortran90 interface to random
- number generation. Computer Physics Communications. Volume 90, Number 1, 117-120
LATEST REVISION : 20/10/2020
function rand_number ()
¶
purpose¶
This function returns a uniformly distributed random number between 0 and 1, exclusive of the two endpoints 0 and 1.
Arguments¶
None
Further Details¶
If the CPP macro _RANDOM_WITH0 is used during compilation, this routine may return 0 value.
subroutine random_number_ ( harvest )
¶
purpose¶
This subroutine returns a uniformly distributed random number HARVEST between 0 and 1, exclusive of the two endpoints 0 and 1.
Arguments¶
- HARVEST (OUTPUT) real(stnd)
- A uniformly distributed random real number between 0 and 1.
Further Details¶
If the CPP macro _RANDOM_WITH0 is used during compilation, this routine may return 0 value.
subroutine random_number_ ( harvest )
¶
purpose¶
This subroutine returns a uniformly distributed random vector HARVEST between 0 and 1, exclusive of the two endpoints 0 and 1.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:)
- A uniformly distributed random real vector between 0 and 1.
Further Details¶
If the CPP macro _RANDOM_WITH0 is used during compilation, this routine may return 0 value.
subroutine random_number_ ( harvest )
¶
purpose¶
This subroutine returns a uniformly distributed random matrix HARVEST between 0 and 1, exclusive of the two endpoints 0 and 1.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:,:)
- A uniformly distributed random real matrix between 0 and 1.
Further Details¶
If the CPP macro _RANDOM_WITH0 is used during compilation, this routine may return 0 value.
subroutine random_number_ ( harvest )
¶
purpose¶
This subroutine returns a uniformly distributed random array of dimension 3 HARVEST between 0 and 1, exclusive of the two endpoints 0 and 1.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:,:,:)
- A uniformly distributed random real array of dimension 3 between 0 and 1.
Further Details¶
If the CPP macro _RANDOM_WITH0 is used during compilation, this routine may return 0 value.
subroutine random_number_ ( harvest )
¶
purpose¶
This subroutine returns a uniformly distributed random array of dimension 4 HARVEST between 0 and 1, exclusive of the two endpoints 0 and 1.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:,:,:,:)
- A uniformly distributed random real array of dimension 4 between 0 and 1.
Further Details¶
If the CPP macro _RANDOM_WITH0 is used during compilation, this routine may return 0 value.
subroutine random_number_ ( harvest )
¶
purpose¶
This subroutine returns a uniformly distributed random array of dimension 5 HARVEST between 0 and 1, exclusive of the two endpoints 0 and 1.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:,:,:,:,:)
- A uniformly distributed random real array of dimension 5 between 0 and 1.
Further Details¶
If the CPP macro _RANDOM_WITH0 is used during compilation, this routine may return 0 value.
subroutine random_number_ ( harvest )
¶
purpose¶
This subroutine returns a uniformly distributed random array of dimension 6 HARVEST between 0 and 1, exclusive of the two endpoints 0 and 1.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:,:,:,:,:,:)
- A uniformly distributed random real array of dimension 6 between 0 and 1.
Further Details¶
If the CPP macro _RANDOM_WITH0 is used during compilation, this routine may return 0 value.
subroutine random_number_ ( harvest )
¶
purpose¶
This subroutine returns a uniformly distributed random array of dimension 7 HARVEST between 0 and 1, exclusive of the two endpoints 0 and 1.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:,:,:,:,:,:,:)
- A uniformly distributed random real array of dimension 7 between 0 and 1.
Further Details¶
If the CPP macro _RANDOM_WITH0 is used during compilation, this routine may return 0 value.
function rand_integer32 ()
¶
purpose¶
This function returns a random integer in the interval (-2147483648,2147483647) inclusive of the two endpoints. The returned integer is equivalent to a signed 32-bit integer.
Arguments¶
None
subroutine random_integer32_ ( harvest )
¶
purpose¶
This subroutine returns a random integer in the interval (-2147483648,2147483647) inclusive of the two endpoints. The returned integer is equivalent to a signed 32-bit integer.
Arguments¶
- HARVEST (OUTPUT) integer(i4b)
- A random integer in the interval (-2147483648,2147483647).
subroutine random_integer32_ ( harvest )
¶
purpose¶
This subroutine returns a vector of random integers in the interval (-2147483648,2147483647) inclusive of the two endpoints. The returned integers are equivalent to signed 32-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:)
- A vector of random integers in the interval (-2147483648,2147483647).
subroutine random_integer32_ ( harvest )
¶
purpose¶
This subroutine returns a matrix of random integers in the interval (-2147483648,2147483647) inclusive of the two endpoints. The returned integers are equivalent to signed 32-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:)
- A matrix of random integers in the interval (-2147483648,2147483647).
subroutine random_integer32_ ( harvest )
¶
purpose¶
This subroutine returns an array of random integers in the interval (-2147483648,2147483647) inclusive of the two endpoints. The returned integers are equivalent to signed 32-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:,:)
- An array of random integers in the interval (-2147483648,2147483647).
subroutine random_integer32_ ( harvest )
¶
purpose¶
This subroutine returns an array of random integers in the interval (-2147483648,2147483647) inclusive of the two endpoints. The returned integers are equivalent to signed 32-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:,:,:)
- An array of random integers in the interval (-2147483648,2147483647).
subroutine random_integer32_ ( harvest )
¶
purpose¶
This subroutine returns an array of random integers in the interval (-2147483648,2147483647) inclusive of the two endpoints. The returned integers are equivalent to signed 32-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:,:,:,:)
- An array of random integers in the interval (-2147483648,2147483647).
subroutine random_integer32_ ( harvest )
¶
purpose¶
This subroutine returns an array of random integers in the interval (-2147483648,2147483647) inclusive of the two endpoints. The returned integers are equivalent to signed 32-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:,:,:,:,:)
- An array of random integers in the interval (-2147483648,2147483647).
subroutine random_integer32_ ( harvest )
¶
purpose¶
This subroutine returns an array of random integers in the interval (-2147483648,2147483647) inclusive of the two endpoints. The returned integers are equivalent to signed 32-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:,:,:,:,:,:)
- An array of random integers in the interval (-2147483648,2147483647).
function rand_integer31 ()
¶
purpose¶
This function returns a random integer in the interval (0,2147483647) inclusive of the two endpoints. The returned integer is equivalent to an unsigned 31-bit integer.
Arguments¶
None
subroutine random_integer31_ ( harvest )
¶
purpose¶
This subroutine returns a random integer in the interval (0,2147483647) inclusive of the two endpoints. The returned integer is equivalent to an unsigned 31-bit integer.
Arguments¶
- HARVEST (OUTPUT) integer(i4b)
- A random integer in the interval (0,2147483647).
subroutine random_integer31_ ( harvest )
¶
purpose¶
This subroutine returns a vector of random integers in the interval (0,2147483647) inclusive of the two endpoints. The returned integers are equivalent to unsigned 31-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:)
- A vector of random integers in the interval (0,2147483647).
subroutine random_integer31_ ( harvest )
¶
purpose¶
This subroutine returns a matrix of random integers in the interval (0,2147483647) inclusive of the two endpoints. The returned integers are equivalent to unsigned 31-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:)
- A matrix of random integers in the interval (0,2147483647).
subroutine random_integer31_ ( harvest )
¶
purpose¶
This subroutine returns an array of random integers in the interval (0,2147483647) inclusive of the two endpoints. The returned integers are equivalent to unsigned 31-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:,:)
- An array of random integers in the interval (0,2147483647).
subroutine random_integer31_ ( harvest )
¶
purpose¶
This subroutine returns an array of random integers in the interval (0,2147483647) inclusive of the two endpoints. The returned integers are equivalent to unsigned 31-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:,:,:)
- An array of random integers in the interval (0,2147483647).
subroutine random_integer31_ ( harvest )
¶
purpose¶
This subroutine returns an array of random integers in the interval (0,2147483647) inclusive of the two endpoints. The returned integers are equivalent to unsigned 31-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:,:,:,:)
- An array of random integers in the interval (0,2147483647).
subroutine random_integer31_ ( harvest )
¶
purpose¶
This subroutine returns an array of random integers in the interval (0,2147483647) inclusive of the two endpoints. The returned integers are equivalent to unsigned 31-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:,:,:,:,:)
- An array of random integers in the interval (0,2147483647).
subroutine random_integer31_ ( harvest )
¶
purpose¶
This subroutine returns an array of random integers in the interval (0,2147483647) inclusive of the two endpoints. The returned integers are equivalent to unsigned 31-bit integers.
Arguments¶
- HARVEST (OUTPUT) integer(i4b), dimension(:,:,:,:,:,:,:)
- An array of random integers in the interval (0,2147483647).
subroutine init_mt19937 ( seed )
¶
purpose¶
User interface subroutine for initializing the state of the MT19937 Random Number Generator (RNG) with a scalar seed, directly, without using the subroutine RANDOM_SEED_ and its interface.
Arguments¶
- SEED (INPUT) integer(i4b)
- On entry, a scalar integer that will be used to initialize the MT19937 RNG.
Further Details¶
Only the first 32 bits of the scalar SEED will be used.
For more informations on the MT19937 RNG, see:
- Matsumoto, M., and Nishimura, T., 1998:
- Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator. ACM Trans. on Modeling and Computer Simulation, Vol. 8, No. 1, January pp.3-30
subroutine init_mt19937 ( seed )
¶
purpose¶
User interface subroutine for initializing the state of the MT19937 Random Number Generator (RNG) with an array of seeds, directly, without using the subroutine RANDOM_SEED_ and its interface.
Arguments¶
- SEED (INPUT) integer(i4b), dimension(:)
- On entry, a vector of integers that will be used to initialize the MT19937 RNG. If size(SEED) is zero, a default scalar seed will be used instead.
Further Details¶
Only the first 32 bits of each element of the array SEED will be used.
For more informations on the MT19937 RNG, see:
- Matsumoto, M., and Nishimura, T., 1998:
- Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator. ACM Trans. on Modeling and Computer Simulation, Vol. 8, No. 1, January pp.3-30
subroutine init_memt19937 ( seed )
¶
purpose¶
User interface subroutine for initializing the state of the MEMT19937-II Random Number Generator (RNG) with a seed, directly, without using the subroutine RANDOM_SEED_ and its interface.
Arguments¶
- SEED (INPUT) integer(i4b)
- On entry, a scalar integer that will be used to initialize the MEMT19937-II RNG.
Further Details¶
Only the first 32 bits of the scalar SEED will be used.
For more informations on the MEMT19937-II RNG, see:
- Harase, S., 2014:
- On the F2-linear relations of Mersenne Twister pseudorandom number generators. Mathematics and Computers in Simulation, Volume 100, Pages 103-113.
subroutine init_memt19937 ( seed )
¶
purpose¶
User interface subroutine for initializing the state of the MEMT19937-II Random Number Generator (RNG) with an array of seeds, directly, without using the subroutine RANDOM_SEED_ and its interface.
Arguments¶
- SEED (INPUT) integer(i4b), dimension(:)
- On entry, a vector of integers that will be used to initialize the MT19937 RNG. If size(SEED) is zero, a default scalar seed will be used instead.
Further Details¶
Only the first 32 bits of each element of the array SEED will be used.
For more informations on the MEMT19937-II RNG, see:
- Harase, S., 2014:
- On the F2-linear relations of Mersenne Twister pseudorandom number generators. Mathematics and Computers in Simulation, Volume 100, Pages 103-113.
subroutine random_seed_ ( alg, size, put, get )
¶
purpose¶
User interface for seeding the random number routines in module RANDOM.
Syntax is like RANDOM_SEED intrinsic subroutine and a call to RANDOM_SEED_() without arguments initiates a non-repeatable reset of the seeds used by the random number subroutines in module RANDOM.
As for RANDOM_SEED intrinsic subroutine, no more than one argument may be specified in a call to RANDOM_SEED_ .
Arguments¶
- ALG (INPUT, OPTIONAL) integer
On entry, a scalar default integer to select the random number generator used in subsequent calls to subroutines RANDOM_NUMBER_ , RANDOM_INTEGER32_ , RANDOM_INTEGER31_ and functions RAND_NUMBER, RAND_INTEGER32 and RAND_INTEGER31. The possible values are:
- ALG=1 : selects the Marsaglia’s KISS random number generator;
- ALG=2 : selects the fast Marsaglia’s KISS random number generator;
- ALG=3 : selects the L’Ecuyer’s LFSR113 random number generator;
- ALG=4 : selects the Mersenne Twister random number generator;
- ALG=5 : selects the maximally equidistributed Mersenne Twister random number generator;
- ALG=6 : selects the extended precision of the Marsaglia’s KISS random number generator;
- ALG=7 : selects the extended precision of the fast Marsaglia’s KISS random number generator;
- ALG=8 : selects the extended precision of the L’Ecuyer’s LFSR113 random number generator.
- ALG=9 : selects the extended precision of Mersenne Twister random number generator;
- ALG=10 : selects the extended precision of maximally equidistributed Mersenne Twister random number generator;
For other values the random number generator is not changed. The default value is the L’Ecuyer’s LFSR113 random number generator (e.g. ALG=3).
- SIZE (OUTPUT, OPTIONAL) integer
- On exit, the size of the seed array used by the random number generators.
- PUT (INPUT, OPTIONAL) integer, dimension(:)
- On entry, a scalar default integer vector of size at least equal to the size of the seed array (as returned by a call to RANDOM_SEED_ with argument SIZE) that will be used to reset the seed for subsequent calls to subroutine RANDOM_NUMBER_.
- GET (OUTPUT, OPTIONAL) integer, dimension(:)
- On exit, a scalar default integer vector, which is the current value of the seed array. Argument GET must be of size at least equal to the size of the seed array (as returned by a call to RANDOM_SEED_ with argument SIZE).
Further Details¶
This subroutine is not thread-safe and must not be called in parallel when OPENMP is used. On the other hand, the associated routines RAND_NUMBER, RANDOM_NUMBER_, RAND_INTEGER32, RANDOM_INTEGER32_, RAND_INTEGER31 and RANDOM_INTEGER31_ are thread-safe if used with OpenMP directives and their states will be consistent while called from multiple OpenMP threads.
The Marsaglia’s KISS (Keep It Simple Stupid) random number generator combines:
- The congruential generator x(n) = 69069 cdot x(n-1) + 1327217885 with a period of 2^32;
- A 3-shift shift-register generator with a period of 2^32 - 1;
- Two 16-bit multiply-with-carry generators with a period of 597273182964842497 > 2^59.
The overall period of this KISS random number generator exceeds 2^123. More details on this Marsaglia’s KISS random number generator are available in the references (3) and (4). This generator is also the one used by the intrinsic subroutine RANDOM_NUMBER as implemented in the GNU gfortran compiler.
The “fast” version of the Marsaglia’s KISS random number generator uses only add, shift, exclusive-or and ‘and’ operations to produces exactly the same 32-bit integer output, which C views as unsigned and Fortran views as signed integers. This version avoids multiplication and is probably faster. More details are available in the reference (5).
The LFSR113 random number generator is described in the reference (2). This random number generator has a period length of about 2^113.
The MT19937 Mersenne Twister random number generator is described in the reference (7). This random number generator has a period length of about 2^19937-1, and 623-dimensional equidistribution property is assured.
The MEMT19937-II Mersenne Twister random number generator is described in the reference (8). This random number generator has also a period length of about 2^19937-1, and a new set of parameters is introduced in the tempering phase of MT19937, which gives a maximally equidistributed Mersenne Twister random number generator.
Note that the size of the seed array varies according to the selected random number generator.
For all the random number generators described above, extended precision versions are also available to generate full precision random real numbers of kind STND (up to 63-bit precision), using the method described in the reference (6).
The FORTRAN versions of these random number generators as implemented here require that 32-bits integer type is available on your computer and that 32-bits integers are represented in base 2 with two’s complement notation.
However, the LFSR113, MT19937 and MEMT19937-II Mersenne Twister random number generators will also work if only 64-bits integer type is available on your system, but in that case you must specify the CPP macro _RANDOM_NOINT32 at compilation. Note, however that the other random number generators will not work properly with 64-bits integer type so they cannot be used on such system.
The KISS random number generators also assumed that integer overflows do not cause crashes. These assumptions are checked before using these random number generators.
On the other hand, the LFSR113, MT19937 and MEMT19937-II random number generators do not use integer arithmetic and are free of such assumptions.
This subroutine is adapted from:
- Hennecke, M., 1995:
- A Fortran90 interface to random number generation. Computer Physics Communications, Volume 90, Number 1, 117-120
- L’Ecuyer, P., 1999:
- Tables of Maximally-Equidistributed Combined LFSR Generators. Mathematics of Computation, 68, 225, 261-269.
- Marsaglia, G., 1999:
- Random number generators for Fortran. Posted to the computer-programming-forum. See: http://computer-programming-forum.com/49-fortran/b89977aa62f72ee8.htm
- Marsaglia, G., 2005:
- Double precision RNGs. Posted to the electronic billboard sci.math.num-analysis. See: http://sci.tech-archive.net/Archive/sci.math.num-analysis/2005-11/msg00352.html
- Marsaglia, G., 2007:
- Fortran and C: United with a KISS. Posted to the Google comp.lang.forum . See: http://groups.google.co.uk/group/comp.lang.fortran/msg/6edb8ad6ec5421a5
- Doornik, J.A, 2007:
- Conversion of high-period random number to floating point. ACM Transactions on Modeling and Computer Simulation, Volume 17, Issue 1, Article No. 3.
- Matsumoto, M., and Nishimura, T., 1998:
- Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator. ACM Trans. on Modeling and Computer Simulation, Vol. 8, No. 1, January pp.3-30
- Harase, S., 2014:
- On the F2-linear relations of Mersenne Twister pseudorandom number generators. Mathematics and Computers in Simulation, Volume 100, Pages 103-113.
function normal_rand_number ()
¶
purpose¶
This function returns a Gaussian distributed random real number.
Arguments¶
None
Further Details¶
This function uses the Cumulative Density Function (CDF) inversion method to generate a Gaussian random real number. Starting with a random number produced by the STATPACK uniform random number generator that can produce random numbers with the uniform distribution over the continuous range (0, 1) (denoted U(0, 1)), the CDF method simply inverts the CDF of a standard Gaussian distribution to produce a standard Gaussian (e.g. a Gaussian distribution with mean zero and standard deviation one) random real number.
The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 2 and 3) by the subroutine PPND7 described in the reference (1). This method gives 7 decimal digits of accuracy in the range [10**(-316), 1-10**(-316)] if computations are done in double or higher precision.
For more details on Uniform and Gaussian random number generators or the approximation of the inverse Gaussian CDF used here, see:
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Wichura, M.J., 1988:
- Algorithm AS 241: The percentage points of the normal distribution. Appl. Statis. 37, 3, 477-484.
subroutine normal_random_number_ ( harvest )
¶
purpose¶
This subroutine returns a random real number HARVEST following the standard Gaussian distribution.
Arguments¶
- HARVEST (OUTPUT) real(stnd)
- A Gaussian distributed random real number.
Further Details¶
This subroutine uses the Cumulative Density Function (CDF) inversion method to generate a Gaussian random real number. Starting with a random number produced by the STATPACK uniform random number generator that can produce random numbers with the uniform distribution over the continuous range (0, 1) (denoted U(0, 1)), the CDF method simply inverts the CDF of a standard Gaussian distribution to produce a standard Gaussian (e.g. a Gaussian distribution with mean zero and standard deviation one) random real number.
The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 2 and 3) by the subroutine PPND7 described in the reference (1). This method gives 7 decimal digits of accuracy in the range [10**(-316), 1-10**(-316)] if computations are done in double or higher precision.
For more details on Uniform and Gaussian random number generators or the approximation of the inverse Gaussian CDF used here, see:
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Wichura, M.J., 1988:
- Algorithm AS 241: The percentage points of the normal distribution. Appl. Statis. 37, 3, 477-484.
subroutine normal_random_number_ ( harvest )
¶
purpose¶
This subroutine returns a random vector HARVEST following the standard normal (Gaussian) distribution.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:)
- A Gaussian distributed random real vector.
Further Details¶
This subroutines uses the Cumulative Density Function (CDF) inversion method to generate Gaussian random numbers. Starting with random numbers produced by the STATPACK uniform random number generator that can produce random numbers with the uniform distribution over the continuous range (0, 1) (denoted U(0, 1)), the CDF method simply inverts the CDF of a standard Gaussian distribution to produce standard Gaussian (e.g. a Gaussian distribution with mean zero and standard deviation one) random numbers.
The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 2 and 3) by the subroutine PPND7 described in the reference (1). This method gives 7 decimal digits of accuracy in the range [10**(-316), 1-10**(-316)] if computations are done in double or higher precision.
For more details on Uniform and Gaussian random number generators or the approximation of the inverse Gaussian CDF used here, see:
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Wichura, M.J., 1988:
- Algorithm AS 241: The percentage points of the normal distribution. Appl. Statis. 37, 3, 477-484.
subroutine normal_random_number_ ( harvest )
¶
purpose¶
This subroutine returns a random matrix HARVEST following the standard normal (Gaussian) distribution.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:,:)
- A Gaussian distributed random real matrix.
Further Details¶
This subroutines uses the Cumulative Density Function (CDF) inversion method to generate Gaussian random real numbers. Starting with random numbers produced by the STATPACK uniform random number generator that can produce random numbers with the uniform distribution over the continuous range (0, 1) (denoted U(0, 1)), the CDF method simply inverts the CDF of a standard Gaussian distribution to produce standard Gaussian (e.g. a Gaussian distribution with mean zero and standard deviation one) random real numbers.
The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 2 and 3) by the subroutine PPND7 described in the reference (1). This method gives 7 decimal digits of accuracy in the range [10**(-316), 1-10**(-316)] if computations are done in double or higher precision.
For more details on Uniform and Gaussian random number generators or the approximation of the inverse Gaussian CDF used here, see :
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Wichura, M.J., 1988:
- Algorithm AS 241: The percentage points of the normal distribution. Appl. Statis. 37, 3, 477-484.
function normal_rand_number2 ()
¶
purpose¶
This function returns a Gaussian distributed real random number of kind extd.
Arguments¶
None
Further Details¶
This function uses the Cumulative Density Function (CDF) inversion method to generate a Gaussian random real number of kind extd. Starting with a random number produced by the STATPACK uniform random number generator that can produce random numbers with the uniform distribution over the continuous range (0, 1) (denoted U(0, 1)), the CDF method simply inverts the CDF of a standard Gaussian distribution to produce a standard Gaussian (e.g. a Gaussian distribution with mean zero and standard deviation one) random real number of kind extd.
The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 7) by the subroutine PPND16 described in the reference (1). This method gives about 16 decimal digits of accuracy in the range [10**(-316), 1-10**(-316)] if computations are done in double or higher precision.
This function is more accurate than NORMAL_RAND_NUMBER function, but it is slower.
For more details on Uniform and Gaussian random number generators or the approximation of the inverse Gaussian CDF used here, see :
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Wichura, M.J., 1988:
- Algorithm AS 241: The percentage points of the normal distribution. Appl. Statis. 37, 3, 477-484.
subroutine normal_random_number2_ ( harvest )
¶
purpose¶
This subroutine returns a Gaussian distributed real random number of kind extd.
Arguments¶
- HARVEST (OUTPUT) real(extd)
- A Gaussian distributed random real number of kind extd.
Further Details¶
This subroutine uses the Cumulative Density Function (CDF) inversion method to generate a Gaussian random ral number of kind extd. Starting with a random number produced by the STATPACK uniform random number generator that can produce random numbers with the uniform distribution over the continuous range (0, 1) (denoted U(0, 1)), the CDF method simply inverts the CDF of a standard Gaussian distribution to produce a standard Gaussian (e.g. a Gaussian distribution with mean zero and standard deviation one) random real number of kind extd.
The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 7) by the subroutine PPND16 described in the reference (1). This method gives about 16 decimal digits of accuracy in the range [10**(-316), 1-10**(-316)] if computations are done in double or higher precision.
For more details on Uniform and Gaussian random number generators or the approximation of the inverse Gaussian CDF used here, see :
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Wichura, M.J., 1988:
- Algorithm AS 241: The percentage points of the normal distribution. Appl. Statis. 37, 3, 477-484.
subroutine normal_random_number2_ ( harvest )
¶
purpose¶
This subroutine returns a random real vector of kind extd following the standard normal (Gaussian) distribution.
Arguments¶
- HARVEST (OUTPUT) real(extd), dimension(:)
- A Gaussian distributed random real vector of kind extd.
Further Details¶
This subroutines uses the Cumulative Density Function (CDF) inversion method to generate Gaussian random real numbers of kind extd. Starting with random numbers produced by the STATPACK uniform random number generator that can produce random numbers with the uniform distribution over the continuous range (0, 1) (denoted U(0, 1)), the CDF method simply inverts the CDF of a standard Gaussian distribution to produce standard Gaussian (e.g. a Gaussian distribution with mean zero and standard deviation one) random real numbers of kind extd.
The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 7) by the subroutine PPND16 described in the reference (1). This method gives about 16 decimal digits of accuracy in the range [10**(-316), 1-10**(-316)] if computations are done in double or higher precision.
For more details on Uniform and Gaussian random number generators or the approximation of the inverse Gaussian CDF used here, see :
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Wichura, M.J., 1988:
- Algorithm AS 241: The percentage points of the normal distribution. Appl. Statis. 37, 3, 477-484.
subroutine normal_random_number2_ ( harvest )
¶
purpose¶
This subroutine returns a random real matrix of kind extd following the standard normal (Gaussian) distribution.
Arguments¶
- HARVEST (OUTPUT) real(extd), dimension(:,:)
- A Gaussian distributed random real matrix of kind extd.
Further Details¶
This subroutines uses the Cumulative Density Function (CDF) inversion method to generate Gaussian random real numbers of kind extd. Starting with random numbers produced by the STATPACK uniform random number generator that can produce random numbers with the uniform distribution over the continuous range (0, 1) (denoted U(0, 1)), the CDF method simply inverts the CDF of a standard Gaussian distribution to produce standard Gaussian (e.g. a Gaussian distribution with mean zero and standard deviation one) random real numbers of kind extd.
The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 7) by the subroutine PPND16 described in the reference (1). This method gives about 16 decimal digits of accuracy in the range [10**(-316), 1-10**(-316)] if computations are done in double or higher precision.
For more details on Uniform and Gaussian random number generators or the approximation of the inverse Gaussian CDF used here, see :
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Wichura, M.J., 1988:
- Algorithm AS 241: The percentage points of the normal distribution. Appl. Statis. 37, 3, 477-484.
function normal_rand_number3 ()
¶
purpose¶
This function returns a Gaussian distributed random real number.
Arguments¶
None
Further Details¶
This function uses the classical Box-Muller method to generate a Gaussian random real number.
For more details on Uniform and Gaussian random number generators or the Box-Muller method used here, see :
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Brent, R.P., 1993:
- Fast Normal Random Generators for vector processors. Report TR-CS-93-04, Computer Sciences Laboratory, Australian National University.
subroutine normal_random_number3_ ( harvest )
¶
purpose¶
This subroutine returns a random real number HARVEST following the standard Gaussian distribution.
Arguments¶
- HARVEST (OUTPUT) real(stnd)
- A Gaussian distributed random real number.
Further Details¶
This subroutine uses the classical Box-Muller method to generate a Gaussian random real number.
For more details on Uniform and Gaussian random number generators or the Box-Muller method used here, see :
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Brent, R.P., 1993:
- Fast Normal Random Generators for vector processors. Report TR-CS-93-04, Computer Sciences Laboratory, Australian National University.
subroutine normal_random_number3_ ( harvest )
¶
purpose¶
This subroutine returns a random real vector HARVEST following the standard normal (Gaussian) distribution.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:)
- A Gaussian distributed random real vector.
Further Details¶
This subroutine uses the classical Box-Muller method to generate Gaussian random real numbers. The computations are parallelized if OPENMP is used.
For more details on Uniform and Gaussian random number generators or the Box-Muller method used here, see :
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Brent, R.P., 1993:
- Fast Normal Random Generators for vector processors. Report TR-CS-93-04, Computer Sciences Laboratory, Australian National University.
subroutine normal_random_number3_ ( harvest )
¶
purpose¶
This subroutine returns a random matrix HARVEST following the standard normal (Gaussian) distribution.
Arguments¶
- HARVEST (OUTPUT) real(stnd), dimension(:,:)
- A Gaussian distributed random real matrix.
Further Details¶
This subroutine uses the classical Box-Muller method to generate Gaussian random real numbers. The computations are parallelized if OPENMP is used.
For more details on Uniform and Gaussian random number generators or the Box-Muller method used here, see :
- Devroye, L., 1986:
- Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.
- Thomas, D.B., Luk, W., Leong, P.H.W., and Villasenor, J.D., 2007:
- Gaussian random number generators. ACM Comput. Surv. 39, 4, Article 11 (October 2007), 38 pages, DOI = 10.1145/1287620.1287622 (http://doi.acm.org/10.1145/1287620.1287622)
- Brent, R.P., 1993:
- Fast Normal Random Generators for vector processors. Report TR-CS-93-04, Computer Sciences Laboratory, Australian National University.
subroutine random_qr_cmp ( mat, diagr, beta, fillr, initseed )
¶
Purpose¶
RANDOM_QR_CMP generates the first k columns of a pseudo-random QR factorization of a hypothetical real n-by-n matrix MAT, whose elements follow independently a Laplace_Gauss(0;1) distribution (e.g. the standard normal distribution):
MAT = Q * Rwhere Q is a pseudo-random orthogonal matrix following the Haar distribution from the group of orthogonal matrices and R is upper triangular. The upper-diagonal elements of R follow a Laplace_Gauss(0;1) distribution (e.g. the standard normal distribution) and the squares of the diagonal elements of R, (R(i,i))**(2) follow a chi-squared distribution with n-i+1 degrees of freedom.
Arguments¶
- MAT (OUTPUT) real(stnd), dimension(:,:)
On exit, the elements above the diagonal of the array MAT contain the corresponding elements of R if FILLR is present and set to true; otherwise the upper-diagonal elements of MAT are not modified. The elements on and below the diagonal, with the arrays BETA and DIAGR, represent the first k columns (with k=size(MAT,2) and k<=n ) of a pseudo-random orthogonal matrix Q following the Haar distribution as a product of elementary reflectors and a diagonal matrix, see Further Details.
- The shape of MAT must verify :
- size( MAT, 2 ) <= size( MAT, 1 ).
- DIAGR (OUTPUT) real(stnd), dimension(:)
On exit, the diagonal elements of the matrix R, see Further Details.
- The size of DIAGR must verify :
- size( DIAGR ) = size( MAT, 2 ) <= size( MAT, 1 ).
- BETA (OUTPUT) real(stnd), dimension(:)
On exit, the scalars factors of the elementary reflectors, see Further Details.
- The size of BETA must verify :
- size( BETA ) = size( DIAGR ) = size( MAT, 2 ) <= size( MAT, 1 ).
- FILLR (INPUT, OPTIONAL) logical(lgl)
on entry, if FILLR is set to true, the super-diagonal elements of R are generated in MAT on exit. If FILLR is set to false, the super-diagonal elements of MAT are not modified or used.
The default is FILLR = false.
- INITSEED (INPUT, OPTIONAL) logical(lgl)
On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiates a non-repeatable reset of the seed used by the STATPACK random generator.
The default is INITSEED = false.
Further Details¶
RANDOM_QR_CMP uses the method described in the reference (1), based on Householder transformations, for generating the first k columns of a n-by-n pseudo-random orthogonal matrices Q distributed according to the Haar measure over the orthogonal group of order n, in a factored form.
The pseudo-random orthogonal matrix Q is represented as a product of n elementary reflectors and of a diagonal matrix
Q = H(1) * H(2) * … * H(n) * diag(sign(DIAGR))Each H(i) has the form
H(i) = I + beta * ( v * v’ ) ,where beta is a real scalar and v is a real n-element vector with v(1:i-1) = 0. v(i:n) is stored on exit in MAT(i:n,i) and beta in BETA(i).
diag(sign(DIAGR)) is the n-by-n diagonal matrix with diagonal elements equal to sign( one, DIAGR) (e.g. its ith diagonal element equals to one if DIAGR(i) is positive and -one otherwise).
It is possible to compute only the first k columns of Q and R by restricting the number of columns of MAT and the sizes of DIAGR and BETA on entry of the subroutine.
Finally, note that the computations are parallelized if OPENMP is used.
Q can be generated with the help of subroutine ORTHO_GEN_RANDOM_QR or can be applied to a vector or a matrix with the help of subroutine APPLY_Q_QR.
For further details on the QR factorization and its use or pseudo-random orthogonal matrices distributed according to the Haar measure over the orthogonal group, see:
- Stewart, G.W., 1980:
- The efficient generation of random orthogonal matrices with an application to condition estimators. SIAM J. Numer. Anal., 17, 403-409
- Lawson, C.L., and Hanson, R.J., 1974:
- Solving least square problems. Prentice-Hall.
- Golub, G.H., and Van Loan, C.F., 1996:
- Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
subroutine ortho_gen_random_qr ( mat, diagr, beta )
¶
Purpose¶
ORTHO_GEN_RANDOM_QR generates a n-by-n real pseudo-random orthogonal matrix following the Haar distribution, which is defined as the product of n elementary reflectors of order n and of a n-by-n diagonal matrix with diagonal elements equal to sign( one, DIAGR):
Q = H(1) * H(2) * … * H(n) * diag(sign(DIAGR))as returned by RANDOM_QR_CMP.
Optionnally, it is possible to generate only the first k columns of Q by restricting arguments MAT, BETA and DIAGR to the first k columns or elements of the corresponding arguments as returned by RANDOM_QR_CMP.
Arguments¶
- MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,…,k, ( with k<=n and k=size(MAT,2) ) as returned by RANDOM_QR_CMP in its array argument MAT. On exit, the first k columns of the pseudo-random n-by-n orthogonal matrix Q.
The shape of MAT must verify : size( MAT, 2 ) <= size( MAT, 1 ).
- DIAGR (INPUT) real(stnd), dimension(:)
On entry, the diagonal elements of the matrix R, as returned by RANDOM_QR_CMP in its argument DIAGR.
- The size of DIAGR must verify :
- size( DIAGR ) = size( MAT, 2 ).
- BETA (INPUT) real(stnd), dimension(:)
On entry, BETA(i) must contain the scalar factor of the elementary reflector H(i), as returned by RANDOM_QR_CMP in its argument BETA.
- The size of BETA must verify :
- size( BETA ) = size( DIAGR ) = size( MAT, 2 ) .
Further Details¶
This subroutine used a blocked algorithm for agregating the Householder transformations (e.g. the elementary reflectors) stored in the lower triangle of MAT and generating the pseudo-random orthogonal matrix Q of the QR factorization returned by subroutine RANDOM_QR_CMP.
Furthermore, the computations are parallelized if OPENMP is used.
For further details on the QR factorization and its use or the blocked algorithm used here, see:
- Lawson, C.L., and Hanson, R.J., 1974:
- Solving least square problems. Prentice-Hall.
- Golub, G.H., and Van Loan, C.F., 1996:
- Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
- Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:
- Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.
- Walker, H.F., 1988:
- Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.
subroutine gen_random_sym_mat ( eigval, mat, eigvec, initseed )
¶
Purpose¶
GEN_RANDOM_SYM_MAT generates a pseudo-random n-by-n real symmetric matrix with prescribed eigenvalues.
Optionally, the corresponding eigenvectors of the generated pseudo-random n-by-n real symmetric matrix can be output if required.
Arguments¶
- EIGVAL (INPUT) real(stnd), dimension(:)
On entry, the prescribed eigenvalues of the pseudo-random n-by-n real symmetric matrix. IF size(EIGVAL)<n, the other eigenvalues are assumed to be zero.
- The size of EIGVAL must verify :
- size( EIGVAL ) <= size( MAT, 1 ) = size( MAT, 2 ) .
- MAT (OUTPUT) real(stnd), dimension(:,:)
On exit, the pseudo-random n-by-n real symmetric matrix.
The shape of MAT must verify size( MAT, 2 ) = size( MAT, 1 ).
- EIGVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)
On exit, the pseudo-random eigenvectors corresponding to the eigenvalues prescribed in EIGVAL. The eigenvectors are returned columnwise.
The shape of EIGVEC must verify:
- size( EIGVEC, 1 ) = size( MAT, 1 ) = size( MAT, 2 )
- size( EIGVEC, 2 ) = size( EIGVAL ) .
- INITSEED (INPUT, OPTIONAL) logical(lgl)
On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiates a non-repeatable reset of the seed used by the STATPACK random generator.
The default is INITSEED = false.
Further Details¶
Pseudo-random eigenvectors are generated as a pseudo-random orthogonal matrix following the Haar distribution from the group of orthogonal matrices and are computed with the help of subroutines RANDOM_QR_CMP and ORTHO_GEN_RANDOM_QR.
These computed pseudo-random eigenvectors and the corresponding prescribed eigenvalues are then used to generate a pseudo-random n-by-n symmetric matrix whose eigenvalues are the prescribed eigenvalues.
These computations are parallelized if OPENMP is used.
For further details, see:
- Stewart, G.W., 1980:
- The efficient generation of random orthogonal matrices with an application to condition estimators. SIAM J. Numer. Anal., 17, 403-409
subroutine gen_random_mat ( singval, mat, leftvec, rightvec, initseed )
¶
Purpose¶
GEN_RANDOM_MAT generates a pseudo-random m-by-n real matrix with prescribed singular values.
Optionnally, the corresponding singular vectors of the generated pseudo-random m-by-n real matrix can be output if required.
Arguments¶
- SINGVAL (INPUT) real(stnd), dimension(:)
On entry , the prescribed singular values of the pseudo-random m-by-n real matrix. IF size(SINGVAL)<min(m,n), the other eigenvalues are assumed to be zero.
- The size of SINGVAL must verify :
- size( SINGVAL ) <= min( size( MAT, 1 ), size( MAT, 2 ) ) .
- MAT (OUTPUT) real(stnd), dimension(:,:)
- On exit, the pseudo-random m-by-n real matrix.
- LEFTVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)
On exit, the pseudo-random left singular vectors corresponding to the singular values prescribed in SINGVAL. The left singular vectors are returned columnwise.
The shape of LEFTVEC must verify:
- size( LEFTVEC, 1 ) = size( MAT, 1 )
- size( LEFTVEC, 2 ) = size( SINGVAL ) .
- RIGHTVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)
On exit, the pseudo-random right singular vectors corresponding to the singular values prescribed in SINGVAL. The right singular vectors are returned columnwise.
The shape of RIGHTVEC must verify:
- size( RIGHTVEC, 1 ) = size( MAT, 2 )
- size( RIGHTVEC, 2 ) = size( SINGVAL ) .
- INITSEED (INPUT, OPTIONAL) logical(lgl)
On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiates a non-repeatable reset of the seed used by the STATPACK random generator.
The default is INITSEED = false.
Further Details¶
Pseudo-random singular vectors are generated as pseudo-random orthogonal matrices following the Haar distribution from the group of orthogonal matrices and are computed with the help of subroutines RANDOM_QR_CMP and ORTHO_GEN_RANDOM_QR.
These computed pseudo-random singular vectors and the corresponding prescribed singular values are then used to generate a pseudo-random m-by-n real matrix whose singular values are the prescribed singular values.
These computations are parallelized if OPENMP is used.
For further details, see:
- Stewart, G.W., 1980:
- The efficient generation of random orthogonal matrices with an application to condition estimators. SIAM J. Numer. Anal., 17, 403-409
subroutine simple_shuffle ( vec )
¶
purpose¶
This subroutine shuffles all the elements of the real vector VEC.
Arguments¶
- VEC (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the real vector to be shuffled.
On exit, the permuted real vector.
Further Details¶
For more details and algorithm, see:
- Noreen, E.W., 1989:
- Computer-intensive methods for testing hypotheses: an introduction. Wiley and Sons, New York, USA, ISBN:978-0-471-61136-3
subroutine simple_shuffle ( vec )
¶
purpose¶
This subroutine shuffles all the elements of the complex vector VEC.
Arguments¶
- VEC (INPUT/OUTPUT) complex(stnd), dimension(:)
On entry, the complex vector to be shuffled.
On exit, the permuted complex vector.
Further Details¶
For more details and algorithm, see:
- Noreen, E.W., 1989:
- Computer-intensive methods for testing hypotheses: an introduction. Wiley and Sons, New York, USA, ISBN:978-0-471-61136-3
subroutine simple_shuffle ( vec )
¶
purpose¶
This subroutine shuffles all the elements of the integer vector VEC.
Arguments¶
- VEC (INPUT/OUTPUT) intger(i4b), dimension(:)
On entry, the integer vector to be shuffled.
On exit, the permuted integer vector.
Further Details¶
For more details and algorithm, see:
- Noreen, E.W., 1989:
- Computer-intensive methods for testing hypotheses: an introduction. Wiley and Sons, New York, USA, ISBN:978-0-471-61136-3
subroutine drawsample ( nsample, pop )
¶
purpose¶
This subroutine may be used to draw a sample, without replacement of size NSAMPLE from a population of size SIZE(POP). On output, the integer vector POP(1:NSAMPLE) indicates which observations are included in the sample.
The integer vector POP must be dimensioned at least as large as NSAMPLE in the calling program.
Arguments¶
- NSAMPLE (INPUT) intger(i4b)
- On entry, the size of the sample.
- POP (OUTPUT) integer(i4b), dimension(:)
On exit, the indices of the observations belonging to the sample are in POP(1:NSAMPLE) and the indices of the observations, which are not in the sample are in POP(NSAMPLE+1:).
The size of POP must greater or equal to NSAMPLE. If this condition is not meet POP(:) is set to -1.
Further Details¶
For more details and algorithm, see:
- Noreen, E.W., 1989:
- Computer-intensive methods for testing hypotheses: an introduction. Wiley and Sons, New York, USA, ISBN:978-0-471-61136-3
subroutine drawbootsample ( npop, sample )
¶
purpose¶
This subroutine may be used to draw a bootstrap random sample of size SIZE(SAMPLE) from a population of size NPOP. On output, the integer vector SAMPLE indicates which observations are included in the bootstrap sample.
Arguments¶
- NPOP (INPUT) intger(i4b)
- On entry, the size of the population.
- SAMPLE (OUTPUT) integer(i4b), dimension(:)
- On exit, the indices of the observations belonging to the sample.
Further Details¶
The sampling is done with replacement, meaning that the sample may contain duplicate observations.
For more details and algorithm, see:
- Noreen, E.W., 1989:
- Computer-intensive methods for testing hypotheses: an introduction. Wiley and Sons, New York, USA, ISBN:978-0-471-61136-3