Module_Random

Copyright 2024 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.

Authors: Pascal Terray (LOCEAN/IPSL, Paris, France)

---

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. FINALLY, RANDOMIZED LINEAR ALGEBRA ROUTINES LIKE RANDOMIZED QB, QR, COLUMN INTERPOLATIVE, TWO_SIDED INTERPOLATIVE AND CUR DECOMPOSITIONS ARE ALSO INCLUDED.

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 : 12/03/2024

---

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:

1. 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:

1. 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:

1. 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:

1. 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:

1. The congruential generator x(n) = 69069 cdot x(n-1) + 1327217885 with a period of 2^32;
2. A 3-shift shift-register generator with a period of 2^32 - 1;
3. 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:

1. Hennecke, M., 1995:

   A Fortran90 interface to random number generation. Computer Physics Communications, Volume 90, Number 1, 117-120

2. L’Ecuyer, P., 1999:

   Tables of Maximally-Equidistributed Combined LFSR Generators. Mathematics of Computation, 68, 225, 261-269.

3. Marsaglia, G., 1999:

   Random number generators for Fortran. Posted to the computer-programming-forum. See: http://computer-programming-forum.com/49-fortran/b89977aa62f72ee8.htm

4. 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

5. 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

6. 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.

7. 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

8. 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:

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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:

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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:

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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 :

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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 :

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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 :

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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 :

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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 :

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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 :

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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 :

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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 :

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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 :

1. Devroye, L., 1986:

   Non-Uniform Random Variate Generation. Springer-Verlag, http://cg.scs.carleton.ca/~luc/rnbookindex.html, New York.

2. 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)

3. 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 * R

where 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=size(MAT,1) ) 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 initiate 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 in module QR_Procedures.

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:

1. Stewart, G.W., 1980:

   The efficient generation of random orthogonal matrices with an application to condition estimators. SIAM J. Numer. Anal., 17, 403-409

2. Lawson, C.L., and Hanson, R.J., 1974:

   Solving least square problems. Prentice-Hall.

3. 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.

Optionally, 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=size(MAT,1) 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:

1. Lawson, C.L., and Hanson, R.J., 1974:

   Solving least square problems. Prentice-Hall.

2. Golub, G.H., and Van Loan, C.F., 1996:

   Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

3. 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.

4. 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, hous )

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 remaining eigenvalues are assumed to be zero.

The size of EIGVAL must verify:

• size( EIGVAL ) <= size( MAT, 1 ) = size( MAT, 2 ) = n .

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 ) = n .

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 ) = n ;
• 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 initiate a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED = false.

HOUS (INPUT, OPTIONAL) logical(lgl)

On entry, if:

• HOUS=true, the pseudo-random eigenvectors are generated as Householder transformation matrices ;
• HOUS=false, the pseudo-random eigenvectors are generated as pseudo-random orthogonal matrices following the Haar distribution from the group of orthogonal matrices.

The default is HOUS = true.

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 if the logical optional argument HOUS has the value false. If the logical optional argument HOUS has the value true, they are generated as Householder transformation matrices with the help of function GEN_RANDOM_ORTHO_MAT.

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:

For further details on Householder matrices and pseudo-random orthogonal matrices following the Haar distribution, see:

1. Lawson, C.L., and Hanson, R.J., 1974:

   Solving least square problems. Prentice-Hall.

2. Golub, G.H., and Van Loan, C.F., 1996:

   Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

3. 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, hous )

Purpose

GEN_RANDOM_MAT generates a pseudo-random m-by-n real matrix with prescribed singular values.

Optionally, 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.

The prescribed singular values must be greater or equal to zero. Furthermore, iF size(SINGVAL)<min(m,n), the remaining singular values are assumed to be zero.

The size of SINGVAL must verify:

• size( SINGVAL ) <= min( size( MAT, 1 ), size( MAT, 2 ) ) = min( m, n) .

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 ) = m ;
• 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 ) = n ;
• 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 initiate a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED = false.

HOUS (INPUT, OPTIONAL) logical(lgl)

On entry, if:

• HOUS=true, the pseudo-random singular vectors are generated as Householder transformation matrices ;
• HOUS=false, the pseudo-random singular vectors are generated as pseudo-random orthogonal matrices following the Haar distribution from the group of orthogonal matrices.

The default is HOUS = true.

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 if the logical optional argument HOUS has the value false. If the logical optional argument HOUS has the value true, they are generated as Householder transformation matrices with the help of function GEN_RANDOM_ORTHO_MAT.

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 on Householder matrices and pseudo-random orthogonal matrices following the Haar distribution, see:

1. Lawson, C.L., and Hanson, R.J., 1974:

   Solving least square problems. Prentice-Hall.

2. Golub, G.H., and Van Loan, C.F., 1996:

   Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

3. Stewart, G.W., 1980:

   The efficient generation of random orthogonal matrices with an application to condition estimators. SIAM J. Numer. Anal., 17, 403-409

function gen_random_ortho_mat ( m, n, initseed )

Purpose

GEN_RANDOM_ORTHO_MAT generates a pseudo-random m-by-n real matrix with orthonormal columns.

Arguments

M (INPUT) integer(i4b)

On entry, the number of rows of the pseudo-random m-by-n real orthonormal matrix.

N (INPUT) integer(i4b)

On entry, the number of columns of the pseudo-random m-by-n real orthonormal matrix.

N must verify:

• N <= M .

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 initiate a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED = false.

Further Details

The pseudo-random m-by-n real matrix with orthonormal columns is generated from the first n columns of an m-by-m Householder matrix H (e.g. a Householder transformation), which is an orthogonal matrix defined as

H = I + tau * ( v * v’ ) ,

where tau is a real scalar and v is a real m-element vector. The scalar tau and vector v are computed to annilhate all the elements of a m-element pseudo-random uniform vector, excepted the first element.

These computations are parallelized if OPENMP is used.

For further details, on Householder matrices, see:

1. Lawson, C.L., and Hanson, R.J., 1974:

   Solving least square problems. Prentice-Hall.

2. Golub, G.H., and Van Loan, C.F., 1996:

   Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

subroutine partial_rqr_cmp ( mat, diagr, beta, ip, krank, tol, tau, rng_alg, blk_size, nover )

Purpose

PARTIAL_RQR_CMP computes a randomized (partial or complete) QR factorization with column pivoting or orthogonal factorization of a m-by-n matrix MAT:

MAT * P = Q * R

, where P is a n-by-n permutation matrix, R is an upper triangular or trapezoidal (i.e., if n>m) m-by-n matrix and Q is a m-by-m orthogonal matrix.

At the user option, the randomized QR factorization can be only partial, e.g., the subroutine ends when the numbers of columns of Q is equal to a predefined value equals to kpartial = size( DIAGR ) = size( BETA ).

This leads implicitly to the following partition of Q:

[ Q1 Q2 ]

where Q1 is a m-by-kpartial orthonormal matrix and Q2 is a m-by-(m-kpartial) orthonormal matrix orthogonal to Q1, and to the following corresponding partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R11 is a kpartial-by-kpartial triangular matrix, R21 is zero by construction, R12 is a full kpartial-by-(n-kpartial) matrix and R22 is a full (m-kpartial)-by-(n-kpartial) matrix.

Then, if the optional scalar argument TOL is present and:

• is in ]0,1[, the rank of R11 is determined by finding the submatrix of R11 which is defined as the largest leading submatrix whose estimated condition number, in the 1-norm, is less than 1/TOL.
• is equal to 0, the rank of R11, krank, is determined by finding the largest submatrix of R11 such that abs(R11[j,j])>0.

In both cases, the order of this submatrix, krank, is the effective rank of R11 (and MAT if krank is less than kpartial or if krank=kpartial=min(m,n)).

If TOL is absent or outside [0,1[, the rank of R11 is not checked and is assumed to be equal to kpartial.

If krank is less than kpartial, then MAT is not of full rank (i.e., certain columns of MAT(:m,:kpartial) are linear combinations of other columns of MAT(:m,:kpartial)) and krank is also an estimate of the rank of MAT.

This leads to a redefinition of the partition of Q = [ Q1 Q2 ], where Q1 and Q2 are now m-by-krank and m-by-(m-krank) orthonormal matrices, and a corresponding redefinition of the associated partition of R, where R11 is now a krank-by-krank triangular matrix, R21 is again zero by construction, R12 is a full krank-by-(n-krank) matrix and R22 is a (m-krank)-by-(n-krank) matrix.

In a final step, if TAU is present, R22 is considered to be negligible and R12 is annihilated by orthogonal transformations from the right, arriving at the partial or complete orthogonal factorization:

MAT * P = Q * T * Z

, where P is a n-by-n permutation matrix, Q is a m-by-m orthogonal matrix, Z is a n-by-n orthogonal matrix and T is a m-by-n matrix and has the form:

[ T11 T12 ]

[ T21 T22 ]

Here T21 (=R21) and T12 are all zero, T22 (=R22) is considered to be negligible and T11 is a krank-by-krank upper triangular matrix.

On exit, P is stored compactly in the integer vector argument IP and if:

• TAU is absent, PARTIAL_RQR_CMP computes Q and submatrices R11, R12 and R22. Submatrices R11, R12 and R22 are stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA.
• TAU is present, PARTIAL_RQR_CMP computes Q, Z and submatrices T11 and T22 (=R22). Submatrices T11 and T22 are stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA. Z is stored compactly in factored form in the array arguments MAT and TAU.

See Further Details for more information on how the (partial or complete) QR or orthogonal decomposition is stored in MAT on exit.

PARTIAL_RQR_CMP performs the same task as subroutine PARTIAL_QR_CMP in module QR_Procedures, but is much faster on large matrices because of the use of a randomized and blocked “BLAS3” algorithm instead of a standard “BLAS2” algorithm. Note, however, that PARTIAL_RQR_CMP is an effective and efficient way for computing a low-rank approximation of MAT, but is less effective to find the rank of MAT because of the use of randomization. As an illustration, the diagonal elements of R11 are not necessarily of decreasing magnitude when computed by PARTIAL_RQR_CMP, while this property is enforced with PARTIAL_QR_CMP.

On the other hand, PARTIAL_QR_CMP can be used for both tasks, but is much slower than PARTIAL_RQR_CMP.

Note that PARTIAL_RQR_CMP performs also exactly the same task as subroutine PARTIAL_RQR_CMP2 in module Random. The main difference between the two routines is that PARTIAL_RQR_CMP2 recomputes the Gaussian and compression matrices at each iteration of the randomized (partial or complete) QR algorithm, while PARTIAL_RQR_CMP uses an efficient updating formulae to recompute the compression matrix at each iteration. In other words, PARTIAL_RQR_CMP is usually faster than PARTIAL_RQR_CMP2, but may be sligthly less robust. See references (3), (4), (5) and (6) for further information.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real m-by-n matrix to be decomposed.

On exit, MAT has been overwritten by details of its (partial) QR or orthogonal factorization.

See Further Details.

DIAGR (OUTPUT) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R11 or T11.

See Further Details.

The size of DIAGR must verify:

• size( DIAGR ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).

BETA (OUTPUT) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q.

See Further Details.

The size of BETA must verify:

• size( BETA ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).

IP (OUTPUT) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

See Further Details.

The size of IP must be equal to size( MAT, 2 ) = n.

KRANK (OUTPUT) integer(i4b)

On exit, KRANK contains the effective rank of R11, i.e., krank, which is the order of this submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of MAT and is also the rank of MAT if krank is less than kpartial = size( BETA ).

If the computed pseudo-rank, krank, is less than kpartial = size( BETA ), BETA(krank+1:kpartial) and, eventually, TAU(krank+1:kpartial) are set to zero and MAT(krank+1:m,krank+1:n) (e.g., R22=T22) is updated on exit.

In other words, the subroutine outputs a partial QR factorization of rank krank instead of rank kpartial.

In all cases, norm(MAT(krank+1:m,krank+1:n)) gives the error of the associated matrix approximation in the Frobenius norm, on exit.

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, if TOL is present then:

• if TOL is in ]0,1[, the calculations to determine the condition number of R11 are performed. Then, TOL is used to determine the effective pseudo-rank of R11, which is defined as the order of the largest leading triangular submatrix in the partial QR factorization with column pivoting of MAT, whose estimated condition number in the 1-norm is less than 1/TOL. On exit, the reciprocal of the condition number is returned in TOL.
• if TOL=0 is specified, the calculations to determine the condition number of R11 are not performed and crude tests on R(j,j) are done to determine the numerical pseudo-rank of R11. On exit, TOL is not changed.

If TOL is not specified or is outside [0,1[, the calculations to determine the rank of R11 are not performed and this rank is assumed to be equal to kpartial.

TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a (partial) complete orthogonal factorization of MAT is computed. Otherwise, a simple QR factorization with column pivoting of MAT is computed.

On exit, the scalars factors of the elementary reflectors defining the orthogonal matrix Z in the (partial) orthogonal factorization of MAT.

See Further Details.

The size of TAU must verify:

• size( TAU ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).

RNG_ALG (INPUT, OPTIONAL) integer(i4b)

On entry, a scalar integer to select the random (uniform) number generator used to build the random gaussian matrix in the randomized partial QR algorithm.

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 current random number generator and its current state are not changed. Note further, that, on exit, the current random number generator is not reset to its previous value before the call to PARTIAL_RQR_CMP.

See the documentation of subroutine RANDOM_SEED_ in module Random for further information.

BLK_SIZE (INPUT, OPTIONAL) integer(i4b)

On entry, the block size used in the randomized partial QR factorization.

BLK_SIZE must be greater or equal to one and less than min(m,n) and must be set to a much smaller value than min(m,n) usually, depending also on the architecture of the computer.

By default, BLK_SIZE is set to min( BLKSZ_QR, min(m,n) ), where parameter BLKSZ_QR is the default block size for QR related algorithms specified in module Select_Parameters.

NOVER (INPUT, OPTIONAL) integer(i4b)

The oversampling size used in the randomized partial QR algorithm.

NOVER must be positive or null and verifies the relationship:

NOVER + BLK_SIZE <= size(MAT,1)

and is adjusted if necessary to verify this relationship in all cases.

See Further Details and the cited references for the meaning and usefulness of the oversampling size in the partial randomized QR algorithm.

By default, the oversampling size is set to 10.

Further Details

The matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(krank), where krank = size( BETA ) <= min( m , n ).

Each H(i) has the form

H(i) = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).

On exit of PARTIAL_RQR_CMP, the orthonormal matrix Q (or its first n columns) can be computed explicitly by a call to subroutine ORTHO_GEN_QR with arguments MAT and BETA.

The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

On exit, if the optional argument TAU is absent:

• The elements above the diagonal of the array MAT(:krank,:krank) contain the corresponding elements of the triangular matrix R11.
• The elements of the diagonal of R11 are stored in the array DIAGR.
• The submatrix R12 is stored in MAT(:krank,krank+1:n).
• The submatrix R22 is stored in MAT(krank+1:m,krank+1:n).
• krank is stored in the real argument KRANK.

If TAU is present, a partial or complete orthogonal factorization of MAT is computed. The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R (e.g., in R12), is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-krank) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.

The Z n-by-n orthogonal matrix is given by

Z = Z(1) * Z(2) * … * Z(krank)

On exit, if the optional argument TAU is present:

• The scalar tau defining T(k) is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,krank+1:n). The other elements of u(k) are not stored.
• The elements above the diagonal of the array section MAT(:krank,:krank) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array DIAGR.
• The submatrix T22 (=R22) is stored in MAT(krank+1:m,krank+1:n).
• krank is stored in the real argument KRANK.

In both cases, norm(MAT(krank+1:m,krank+1:n)) gives the error of the associated matrix approximation in the Frobenius norm, on exit.

If it is possible that MAT may not be of full rank (i.e., certain columns of MAT are linear combinations of other columns), then the eventual linearly dependent columns in the (partial or complete) QR decomposition of MAT, which is sought, can be determined by using TOL=relative precision of the elements in MAT. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.

The computations are parallelized if OPENMP is used. Note also that PARTIAL_RQR_CMP uses a randomized “BLAS3” algorithm described in the references (3), (4), (5) and (6), which has about the same efficiency of a “BLAS3” QR algorithm without column pivoting and is thus highly efficient on large matrices.

For further details, see:

1. Lawson, C.L., and Hanson, R.J., 1974:

   Solving least square problems. Prentice-Hall.

2. Golub, G.H., and Van Loan, C.F., 1996:

   Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

3. Duersch, J.A., and Gu, M., 2017:

   Randomized QR with column pivoting. SIAM J. Sci. Comput., Volume 39, Issue 4, C263-C291.

4. Martinsson, P.G., Quintana-Orti, G., Heavner, N., and Van de Geijn, R., 2017:

   Householder QR factorization with randomization for column pivoting (HQRRP). SIAM J. Sci. Comput., Volume 39, Issue 2, C96-C115.

5. Xiao, J., Gu, M., and Langou, J., 2017:

   Fast parallel randomized QR with column pivoting algorithms for reliable low-rank matrix approximations. IEEE 24th International Conference on High Performance Computing (HiPC), IEEE, 2017, 233-242.

6. Duersch, J.A., and Gu, M., 2020:

   Randomized projection for rank-revealing matrix factorizations and low-rank approximations. SIAM Review, Volume 62, Issue 3, 661-682.

subroutine partial_rqr_cmp2 ( mat, diagr, beta, ip, krank, tol, tau, rng_alg, blk_size, nover )

Purpose

PARTIAL_RQR_CMP2 computes a randomized (partial or complete) QR factorization with column pivoting or orthogonal factorization of a m-by-n matrix MAT:

MAT * P = Q * R

, where P is a n-by-n permutation matrix, R is an upper triangular or trapezoidal (i.e., if n>m) m-by-n matrix and Q is a m-by-m orthogonal matrix.

At the user option, the randomized QR factorization can be only partial, e.g., the subroutine ends when the numbers of columns of Q is equal to a predefined value equals to kpartial = size( DIAGR ) = size( BETA ).

This leads implicitly to the following partition of Q:

[ Q1 Q2 ]

where Q1 is a m-by-kpartial orthonormal matrix and Q2 is a m-by-(m-kpartial) orthonormal matrix orthogonal to Q1, and to the following corresponding partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R11 is a kpartial-by-kpartial triangular matrix, R21 is zero by construction, R12 is a full kpartial-by-(n-kpartial) matrix and R22 is a full (m-kpartial)-by-(n-kpartial) matrix.

Then, if the optional scalar argument TOL is present and:

• is in ]0,1[, the rank of R11 is determined by finding the submatrix of R11 which is defined as the largest leading submatrix whose estimated condition number, in the 1-norm, is less than 1/TOL.
• is equal to 0, the rank of R11, krank, is determined by finding the largest submatrix of R11 such that abs(R11[j,j])>0.

In both cases, the order of this submatrix, krank, is the effective rank of R11 (and MAT if krank is less than kpartial or if krank=kpartial=min(m,n)).

If TOL is absent or outside [0,1[, the rank of R11 is not checked and is assumed to be equal to kpartial.

If krank is less than kpartial, then MAT is not of full rank (i.e., certain columns of MAT(:m,:kpartial) are linear combinations of other columns of MAT(:m,:kpartial)) and krank is also an estimate of the rank of MAT.

This leads to a redefinition of the partition of Q = [ Q1 Q2 ], where Q1 and Q2 are now m-by-krank and m-by-(m-krank) orthonormal matrices, and a corresponding redefinition of the associated partition of R, where R11 is now a krank-by-krank triangular matrix, R21 is again zero by construction, R12 is a full krank-by-(n-krank) matrix and R22 is a (m-krank)-by-(n-krank) matrix.

In a final step, if TAU is present, R22 is considered to be negligible and R12 is annihilated by orthogonal transformations from the right, arriving at the partial or complete orthogonal factorization:

MAT * P = Q * T * Z

, where P is a n-by-n permutation matrix, Q is a m-by-m orthogonal matrix, Z is a n-by-n orthogonal matrix and T is a m-by-n matrix and has the form:

[ T11 T12 ]

[ T21 T22 ]

Here T21 (=R21) and T12 are all zero, T22 (=R22) is considered to be negligible and T11 is a krank-by-krank upper triangular matrix.

On exit, P is stored compactly in the vector argument IP and if:

• TAU is absent, PARTIAL_RQR_CMP2 computes Q and submatrices R11, R12 and R22. Submatrices R11, R12 and R22 are stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA.
• TAU is present, PARTIAL_RQR_CMP2 computes Q, Z and submatrices T11 and T22 (=R22). Submatrices T11 and T22 are stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA. Z is stored compactly in factored form in the array arguments MAT and TAU.

See Further Details for more information on how the (partial or complete) QR or orthogonal decomposition is stored in MAT on exit.

PARTIAL_RQR_CMP2 performs the same task as subroutine PARTIAL_QR_CMP in module QR_Procedures, but is much faster on large matrices because of the use of a randomized and blocked “BLAS3” algorithm instead of a standard “BLAS2” algorithm. Note, however, that PARTIAL_RQR_CMP2 is an effective and efficient way for computing a low-rank approximation of MAT, but is less effective to find the rank of MAT because of the use of randomization. As an illustration, the diagonal elements of R11 are not necessarily of decreasing magnitude when computed by PARTIAL_RQR_CMP2, while this property is enforced with PARTIAL_QR_CMP.

On the other hand, PARTIAL_QR_CMP can be used for both tasks, but is much slower than PARTIAL_RQR_CMP2.

Note that PARTIAL_RQR_CMP2 performs also exactly the same task as subroutine PARTIAL_RQR_CMP in module Random. The main difference between the two routines is that PARTIAL_RQR_CMP2 recomputes the Gaussian and compression matrices at each iteration of the randomized (partial or complete) QR algorithm, while PARTIAL_RQR_CMP uses an efficient updating formulae to recompute the compression matrix at each iteration. In other words, PARTIAL_RQR_CMP is usually faster than PARTIAL_RQR_CMP2, but may be sligthly less robust. See the references (3), (4), (5) and (6) for further information.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real m-by-n matrix to be decomposed.

On exit, MAT has been overwritten by details of its (partial) QR or orthogonal factorization.

See Further Details.

DIAGR (OUTPUT) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R11 or T11.

See Further Details.

The size of DIAGR must verify:

• size( DIAGR ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).

BETA (OUTPUT) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q.

See Further Details.

The size of BETA must verify:

• size( BETA ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).

IP (OUTPUT) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

See Further Details.

The size of IP must be equal to size( MAT, 2 ) = n.

KRANK (OUTPUT) integer(i4b)

On exit, KRANK contains the effective rank of R11, i.e., krank, which is the order of this submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of MAT and is also the rank of MAT if krank is less than kpartial = size( BETA ).

If the computed pseudo-rank, krank, is less than kpartial = size( BETA ), BETA(krank+1:kpartial) and, eventually, TAU(krank+1:kpartial) are set to zero and MAT(krank+1:m,krank+1:n) (e.g., R22=T22) is updated on exit. In other words, the subroutine outputs a partial QR factorization of rank krank instead of rank kpartial.

In all cases, norm(MAT(krank+1:m,krank+1:n)) gives the error of the associated matrix approximation in the Frobenius norm, on exit.

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, if TOL is present then:

• if TOL is in ]0,1[, the calculations to determine the condition number of R11 are performed. Then, TOL is used to determine the effective pseudo-rank of R11, which is defined as the order of the largest leading triangular submatrix in the partial QR factorization with column pivoting of MAT, whose estimated condition number in the 1-norm is less than 1/TOL. On exit, the reciprocal of the condition number is returned in TOL.
• if TOL=0 is specified, the calculations to determine the condition number of R11 are not performed and crude tests on R(j,j) are done to determine the numerical pseudo-rank of R11. On exit, TOL is not changed.

If TOL is not specified or is outside [0,1[, the calculations to determine the rank of R11 are not performed and this rank is assumed to be equal to kpartial.

TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a (partial) orthogonal factorization of MAT is computed. Otherwise, a simple QR factorization with column pivoting of MAT is computed.

On exit, the scalars factors of the elementary reflectors defining the orthogonal matrix Z in the (partial) orthogonal factorization of MAT.

See Further Details.

The size of TAU must verify:

• size( TAU ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).

RNG_ALG (INPUT, OPTIONAL) integer(i4b)

On entry, a scalar integer to select the random (uniform) number generator used to build the random gaussian matrix in the randomized partial QR algorithm.

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 current random number generator and its current state are not changed. Note further, that, on exit, the current random number generator is not reset to its previous value before the call to PARTIAL_RQR_CMP2.

See the documentation of subroutine RANDOM_SEED_ in module Random for further information.

BLK_SIZE (INPUT, OPTIONAL) integer(i4b)

On entry, the block size used in the randomized partial QR factorization.

BLK_SIZE must be greater or equal to one and less than min(m,n) and must be set to a much smaller value than min(m,n) usually, depending also on the architecture of the computer.

By default, BLK_SIZE is set to min( BLKSZ_QR, min(m,n) ), where parameter BLKSZ_QR is the default block size for QR related algorithms specified in module Select_Parameters.

NOVER (INPUT, OPTIONAL) integer(i4b)

The oversampling size used in the randomized partial QR algorithm.

NOVER must be positive or null and verifies the relationship:

NOVER + BLK_SIZE <= size(MAT,1)

and is adjusted if necessary to verify this relationship in all cases.

See Further Details and the cited references for the meaning and usefulness of the oversampling size in the partial randomized QR algorithm.

By default, the oversampling size is set to 10.

Further Details

The matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(krank), where krank = size( BETA ) <= min( m , n ).

Each H(i) has the form

H(i) = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).

On exit of PARTIAL_RQR_CMP2, the orthonormal matrix Q (or its first n columns) can be computed explicitly by a call to subroutine ORTHO_GEN_QR with arguments MAT and BETA.

The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

On exit, if the optional argument TAU is absent:

• The elements above the diagonal of the array MAT(:krank,:krank) contain the corresponding elements of the triangular matrix R11.
• The elements of the diagonal of R11 are stored in the array DIAGR.
• The submatrix R12 is stored in MAT(:krank,krank+1:n).
• The submatrix R22 is stored in MAT(krank+1:m,krank+1:n).
• krank is stored in the real argument KRANK.

If TAU is present, a partial or complete orthogonal factorization of MAT is computed. The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R (e.g., in R12), is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-krank) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.

The Z n-by-n orthogonal matrix is given by

Z = Z(1) * Z(2) * … * Z(krank)

On exit, if the optional argument TAU is present:

• The scalar tau defining T(k) is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,krank+1:n). The other elements of u(k) are not stored.
• The elements above the diagonal of the array section MAT(:krank,:krank) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array DIAGR.
• The submatrix T22 (=R22) is stored in MAT(krank+1:m,krank+1:n).
• krank is stored in the real argument KRANK.

In both cases, norm(MAT(krank+1:m,krank+1:n)) gives the error of the associated matrix approximation in the Frobenius norm, on exit.

If it is possible that MAT may not be of full rank (i.e., certain columns of MAT are linear combinations of other columns), then the eventual linearly dependent columns in the partial QR decomposition of MAT, which is sought, can be determined by using TOL=relative precision of the elements in MAT. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.

The computations are parallelized if OPENMP is used. Note also that PARTIAL_RQR_CMP2 uses a randomized “BLAS3” algorithm described in the references (3), (4), (5) and (6), which has about the same efficiency of a “BLAS3” QR algorithm without column pivoting and is thus highly efficient on large matrices.

For further details, see:

1. Lawson, C.L., and Hanson, R.J., 1974:

   Solving least square problems. Prentice-Hall.

2. Golub, G.H., and Van Loan, C.F., 1996:

   Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

3. Duersch, J.A., and Gu, M., 2017:

   Randomized QR with column pivoting. SIAM J. Sci. Comput., Volume 39, Issue 4, C263-C291.

4. Martinsson, P.G., Quintana-Orti, G., Heavner, N., and Van de Geijn, R., 2017:

   Householder QR factorization with randomization for column pivoting (HQRRP). SIAM J. Sci. Comput., Volume 39, Issue 2, C96-C115.

5. Xiao, J., Gu, M., and Langou, J., 2017:

   Fast parallel randomized QR with column pivoting algorithms for reliable low-rank matrix approximations. IEEE 24th International Conference on High Performance Computing (HiPC), IEEE, 2017, 233-242.

6. Duersch, J.A., and Gu, M., 2020:

   Randomized projection for rank-revealing matrix factorizations and low-rank approximations. SIAM Review, Volume 62, Issue 3, 661-682.

subroutine partial_rtqr_cmp ( mat, diagr, beta, ip, krank, tol, tau, rng_alg, niter, nover )

Purpose

PARTIAL_RTQR_CMP computes a randomized partial and truncated QR factorization with column pivoting, or an orthogonal factorization, of a m-by-n matrix MAT:

MAT * P = Q * R

, where P is a n-by-n permutation matrix, R is an upper triangular or trapezoidal kpartial-by-n matrix and Q is a m-by-kpartial matrix with orthogonal columns.

The randomized QR factorization is only partial, e.g., the subroutine ends when the numbers of columns of Q is equal to a predefined value equals to kpartial = size( DIAGR ) = size( BETA ) <= min(m,n).

This leads implicitly to the following partition of R:

[ R11 R12 ]

where R11 is a kpartial-by-kpartial triangular matrix and R12 is a full kpartial-by-(n-kpartial) matrix.

Then, if the optional scalar argument TOL is present and:

• is in ]0,1[, the rank of R11 is determined by finding the submatrix of R11 which is defined as the largest leading submatrix whose estimated condition number, in the 1-norm, is less than 1/TOL.
• is equal to 0, the rank of R11, krank, is determined by finding the largest submatrix of R11 such that abs(R11[j,j])>0.

In both cases, the order of this submatrix, krank, is the effective rank of R11 (and MAT if krank is less than kpartial or if krank=kpartial=min(m,n)).

If TOL is absent or outside [0,1[, the rank of R11 is not checked and is assumed to be equal to kpartial.

If krank is less than kpartial, then MAT is not of full rank (i.e., certain columns of MAT(:m,:kpartial) are linear combinations of other columns of MAT(:m,:kpartial)) and krank is also an estimate of the rank of MAT.

This leads to a redefinition of the partition of Q = [ Q1 Q2 ], where Q1 and Q2 are now m-by-krank and m-by-(kpartial-krank) orthonormal matrices, and a corresponding redefinition of the associated partition of R, where R11 is now a krank-by-krank triangular matrix and R12 is a full krank-by-(n-krank) matrix.

In a final step, if TAU is present, R12 is annihilated by orthogonal transformations from the right, arriving at the partial orthogonal factorization:

MAT * P = Q * T * Z

, where P is a n-by-n permutation matrix, Q is a m-by-krank orthogonal matrix, Z is a n-by-n orthogonal matrix and T is a krank-by-n matrix and has the form:

[ T11 T12 ]

Here T12 is all zero and T11 is a krank-by-krank upper triangular matrix.

On exit, P is stored compactly in the vector argument IP and if:

• TAU is absent, PARTIAL_RTQR_CMP computes Q and submatrices R11 and R12. Submatrices R11, and R12 are stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA.
• TAU is present, PARTIAL_RTQR_CMP computes Q, Z and submatrix T11. Submatrix T11 is stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA. Z is stored compactly in factored form in the array arguments MAT and TAU.

See Further Details for more information on how the partial QR or orthogonal decomposition is stored in MAT on exit.

PARTIAL_RTQR_CMP performs the same task as subroutines PARTIAL_RQR_CMP and PARTIAL_RQR_CMP2 subroutines in module Random, but is significantly faster when krank is relatively small and MAT is a very large matrix. This is due to the fact that PARTIAL_RTQR_CMP computes Q (in factored form), R11 and R12, but not R22 (where R22 is the bottom left (m-krank)-by-(n-krank) submatrix in the QR factorization of MAT) as PARTIAL_RQR_CMP and PARTIAL_RQR_CMP2 subroutines. Furthermore, only an estimate of R12 is computed by PARTIAL_RTQR_CMP, using a randomization algorithm described in the reference (7), while the computation of R12 is exact in PARTIAL_RQR_CMP and PARTIAL_RQR_CMP2 subroutines.

In other words, PARTIAL_RTQR_CMP avoids trailing updates of MAT (e.g., R22) during the randomized QR factorization, which reduces significantly the CPU time compared to PARTIAL_RQR_CMP and PARTIAL_RQR_CMP2 subroutines, if krank is small and MAT is a very large matrix, as these subroutines perform these trailing updates.

Note that PARTIAL_RTQR_CMP also does not recompute or update the compression matrix at each iteration as PARTIAL_RQR_CMP and PARTIAL_RQR_CMP2 subroutines. See references (3), (4), (5), (6) and (7) for further information.

In summary, PARTIAL_RTQR_CMP subroutine is usually faster than PARTIAL_RQR_CMP and PARTIAL_RQR_CMP2 subroutines but is less accurate, especially for matrices with a slow decay of their singular values. See reference (7) for details.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real m-by-n matrix to be decomposed.

On exit, MAT has been overwritten by details of its partial and truncated QR or orthogonal factorization.

See Further Details.

DIAGR (OUTPUT) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R11 or T11.

See Further Details.

The size of DIAGR must verify:

• size( DIAGR ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).

BETA (OUTPUT) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q.

See Further Details.

The size of BETA must verify:

• size( BETA ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).

IP (OUTPUT) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

See Further Details.

The size of IP must be equal to size( MAT, 2 ) = n.

KRANK (OUTPUT) integer(i4b)

On exit, KRANK contains the effective rank of R11, i.e., krank, which is the order of this submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of MAT and is also the rank of MAT if krank is less than kpartial = size( BETA ).

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, if TOL is present then:

• if TOL is in ]0,1[, the calculations to determine the condition number of R11 are performed. Then, TOL is used to determine the effective pseudo-rank of R11, which is defined as the order of the largest leading triangular submatrix in the partial QR factorization with column pivoting of MAT, whose estimated condition number in the 1-norm is less than 1/TOL. On exit, the reciprocal of the condition number is returned in TOL.
• if TOL=0 is specified, the calculations to determine the condition number of R11 are not performed and crude tests on R(j,j) are done to determine the numerical pseudo-rank of R11. On exit, TOL is not changed.

If TOL is not specified or is outside [0,1[, the calculations to determine the rank of R11 are not performed and this rank is assumed to be equal to kpartial.

TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a (partial) complete orthogonal factorization of MAT is computed. Otherwise, a simple QR factorization with column pivoting of MAT is computed.

On exit, the scalars factors of the elementary reflectors defining the orthogonal matrix Z in the partial orthogonal factorization of MAT.

See Further Details.

The size of TAU must verify:

• size( TAU ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).

RNG_ALG (INPUT, OPTIONAL) integer(i4b)

On entry, a scalar integer to select the random (uniform) number generator used to build the random gaussian matrix in the randomized partial QR algorithm.

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 current random number generator and its current state are not changed. Note further, that, on exit, the current random number generator is not reset to its previous value before the call to PARTIAL_RTQR_CMP.

See the documentation of subroutine RANDOM_SEED_ in module Random for further information.

NITER (INPUT, OPTIONAL) integer(i4b)

The number of randomized subspace iterations performed in the subroutine for improving the accuracy of the compression matrix before computing its partial QR factorization with column pivoting.

NITER must be positive or null.

By default, 0 randomized subspace iterations are performed.

NOVER (INPUT, OPTIONAL) integer(i4b)

The oversampling size used in the randomized partial QR algorithm.

NOVER must be positive or null and verifies the relationship:

NOVER + kpartial <= size(MAT,1)

and is adjusted if necessary to verify this relationship in all cases.

See Further Details and the cited references for the meaning and usefulness of the oversampling size in the randomized partial and truncated QR algorithm.

By default, the oversampling size is set to max( kpartial/2_i4b, 10 ).

Further Details

The matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(krank), where krank <= min( m , n ).

Each H(i) has the form

H(i) = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).

On exit of PARTIAL_RTQR_CMP, the orthonormal matrix Q (or its first krank columns) can be computed explicitly by a call to subroutine ORTHO_GEN_QR with arguments MAT(:m,:krank) and BETA(:krank).

The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

On exit, if the optional argument TAU is absent:

• The elements above the diagonal of the array MAT(:krank,:krank) contain the corresponding elements of the triangular matrix R11.
• The elements of the diagonal of R11 are stored in the array DIAGR.
• The approximation of the submatrix R12 is stored in MAT(:krank,krank+1:n).
• krank is stored in the real argument KRANK.

If TAU is present, a partial complete orthogonal factorization of MAT is computed. The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R (e.g., in R12), is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-krank) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.

The Z n-by-n orthogonal matrix is given by

Z = Z(1) * Z(2) * … * Z(krank)

On exit, if the optional argument TAU is present:

• The scalar tau defining T(k) is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,krank+1:n). The other elements of u(k) are not stored.
• The elements above the diagonal of the array section MAT(:krank,:krank) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array DIAGR.
• krank is stored in the real argument KRANK.

If it is possible that MAT may not be of full rank (i.e., certain columns of MAT are linear combinations of other columns), then the eventual linearly dependent columns in the partial QR decomposition of MAT, which is sought, can be determined by using TOL=relative precision of the elements in MAT. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.

The computations are parallelized if OPENMP is used. Note also that PARTIAL_RTQR_CMP uses a randomized “BLAS3” algorithm described in the reference (7), which has about the same efficiency of a “BLAS3” QR algorithm without column pivoting and is thus highly efficient on large matrices.

For further details, see:

1. Lawson, C.L., and Hanson, R.J., 1974:

   Solving least square problems. Prentice-Hall.

2. Golub, G.H., and Van Loan, C.F., 1996:

   Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

3. Duersch, J.A., and Gu, M., 2017:

   Randomized QR with column pivoting. SIAM J. Sci. Comput., Volume 39, Issue 4, C263-C291.

4. Martinsson, P.G., Quintana-Orti, G., Heavner, N., and Van de Geijn, R., 2017:

   Householder QR factorization with randomization for column pivoting (HQRRP). SIAM J. Sci. Comput., Volume 39, Issue 2, C96-C115.

5. Xiao, J., Gu, M., and Langou, J., 2017:

   Fast parallel randomized QR with column pivoting algorithms for reliable low-rank matrix approximations. IEEE 24th International Conference on High Performance Computing (HiPC), IEEE, 2017, 233-242.

6. Duersch, J.A., and Gu, M., 2020:

   Randomized projection for rank-revealing matrix factorizations and low-rank approximations. SIAM Review, Volume 62, Issue 3, 661-682.

7. Mary, T., Yamazaki, I., Kurzak, J., Luszczek, P., Tomov, S., and Dongarra, J., 2015:

   Performance of Random Sampling for Computing Low-rank Approximations of a Dense Matrix on GPUs. International Conference for High Performance Computing, Networking, Storage and Analysis (SC’15).

subroutine partial_rqr_cmp_fixed_precision ( mat, relerr, diagr, beta, ip, krank, tau, rng_alg, blk_size, nover )

Purpose

PARTIAL_RQR_CMP_FIXED_PRECISION computes a randomized partial QR factorization with column pivoting or orthogonal factorization of a m-by-n matrix MAT:

MAT * P = Q * R

, where P is a n-by-n permutation matrix, R is a krank-by-n (upper trapezoidal) matrix and Q is a m-by-krank matrix with orthonormal columns. This leads to the following matrix approximation of MAT of rank krank:

MAT = Q * (R * P’)

krank is the target rank of the matrix approximation, which is sought, and this partial factorization must have an approximation error which fulfills:

|| MAT - Q * ( R * P’ ) ||_F <= ||MAT||_F * relerr

|| ||_F is the Frobenius norm and relerr is a prescribed accuracy tolerance for the relative error of the computed matrix approximation, specified in the input argument RELERR.

PARTIAL_RQR_CMP_FIXED_PRECISION searches incrementally the best (e.g., of smallest rank) Q * (R * P’) approximation, which fulfills the prescribed accuracy tolerance for the relative error. More precisely, the rank of the matrix approximation is increased progressively of BLK_SIZE by BLK_SIZE until the prescribed accuracy tolerance is satisfied.

In other words, the rank, krank, of the matrix approximation is not known in advance and is determined in the subroutine. krank is stored in the argument KRANK and the relative error of the computed matrix approximation is output in argument RELERR on exit.

The computed matrix approximation leads implicitly to the following partition of R:

[ R1 R2 ]

where R1 is a krank-by-krank triangular matrix and R2 is a full krank-by-(n-krank) matrix.

In a final step, if TAU is present, R2 is annihilated by orthogonal transformations from the right, arriving at the partial orthogonal factorization:

MAT * P = Q * T1 * Z

, where P is a n-by-n permutation matrix, Q is a m-by-krank matrix with orthonormal columns, Z is a krank-by-n matrix with orthonormal rows and T1 is a krank-by-krank upper triangular matrix.

Note, however, that this final step does not change the matrix approximation and its relative error, only the output format of this matrix approximation, which is now composed of four factors instead of three.

On exit, P is stored compactly in the vector argument IP, krank is stored in the scalar argument KRANK and if:

• TAU is absent, PARTIAL_RQR_CMP_FIXED_PRECISION computes Q and submatrices R1 and R2. Submatrices R1 and R2 are stored in the array arguments MAT and DIAGR(:KRANK). Q is stored compactly in factored form in the array arguments MAT and BETA(:KRANK).
• TAU is present, PARTIAL_RQR_CMP_FIXED_PRECISION computes Q, Z and submatrice T1. Submatrice T1 is stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA(:KRANK). Z is stored compactly in factored form in the array arguments MAT and TAU(:KRANK).

In all cases, the relative error of the computed matrix approximation is output in argument RELERR.

See Further Details for more information.

PARTIAL_RQR_CMP_FIXED_PRECISION performs the same task as subroutine PARTIAL_QR_CMP_FIXED_PRECISION in module QR_Procedures, but is much faster on large matrices because of the use of a randomized and blocked “BLAS3” algorithm instead of a standard “BLAS2” algorithm. Another difference is that, in PARTIAL_RQR_CMP_FIXED_PRECISION, the rank of the matrix approximation is increased progressively of BLK_SIZE by BLK_SIZE until the prescribed tolerance for the relative error is satisfied while in PARTIAL_QR_CMP_FIXED_PRECISION subroutine, the rank of the matrix approximation is increased one by one until the prescribed tolerance for the relative error is satisfied. In other words, the rank of the matrix approximation found by PARTIAL_RQR_CMP_FIXED_PRECISION is always larger than the one found by PARTIAL_QR_CMP_FIXED_PRECISION and is a multiple of BLK_SIZE.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real m-by-n matrix to be decomposed.

On exit, MAT has been overwritten by details of its partial QR or orthogonal factorization.

See Further Details.

RELERR (INPUT/OUTPUT) real(stnd)

On entry, the requested accuracy tolerance for the relative error of the computed partial matrix approximation.

The preset value for RELERR must be greater than 4*epsilon( RELERR ) and less than one.

On exit, RELERR contains the relative error of the computed partial matrix approximation:

• RELERR = || MAT - Q * ( R * P’ ) ||_F / ||MAT||_F

DIAGR (OUTPUT) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R1 or T1 are stored in the array section DIAGR(:KRANK). Other elements of DIAGR are set to zero on exit.

See Further Details.

The size of DIAGR must verify:

• size( DIAGR ) = min( size(MAT,1) , size(MAT,2) ).

BETA (OUTPUT) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q are stored in the array section BETA(:KRANK). Other elements of BETA are set to zero on exit.

See Further Details.

The size of BETA must verify:

• size( BETA ) = min( size(MAT,1) , size(MAT,2) ).

IP (OUTPUT) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

See Further Details.

The size of IP must be equal to size( MAT, 2 ) = n.

KRANK (OUTPUT) integer(i4b)

On exit, KRANK contains the rank of R1, i.e., krank, which is the order of this submatrix R1. This is the same as the order of the submatrix T1 in the “partial” complete orthogonal factorization of MAT and is also the rank of the computed matrix approximation.

In all cases, norm(MAT(KRANK+1:m,KRANK+1:n)) gives the error of the associated matrix approximation in the Frobenius norm, on exit.

TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a partial complete orthogonal factorization of MAT is computed. Otherwise, a simple QR factorization with column pivoting of MAT is computed.

On exit, the scalars factors of the elementary reflectors defining the orthogonal matrix Z in the partial complete orthogonal factorization of MAT are stored in the array section TAU(:KRANK). Other elements of TAU are set to zero on exit.

See Further Details.

The size of TAU must verify:

• size( TAU ) = min( size(MAT,1) , size(MAT,2) ).

RNG_ALG (INPUT, OPTIONAL) integer(i4b)

On entry, a scalar integer to select the random (uniform) number generator used to build the random gaussian matrix in the randomized partial QR algorithm.

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 current random number generator and its current state are not changed. Note further, that, on exit, the current random number generator is not reset to its previous value before the call to PARTIAL_RQR_CMP_FIXED_PRECISION.

See the documentation of subroutine RANDOM_SEED_ in module Random for further information.

BLK_SIZE (INPUT, OPTIONAL) integer(i4b)

On entry, the block size used in the randomized partial QR factorization.

BLK_SIZE must be greater or equal to one and less than min(m,n) and must be set to a much smaller value than min(m,n) usually, depending also on the architecture of the computer.

By default, BLK_SIZE is set to min( BLKSZ_QR, min(m,n) ), where parameter BLKSZ_QR is the default block size for QR related algorithms specified in module Select_Parameters.

NOVER (INPUT, OPTIONAL) integer(i4b)

The oversampling size used in the randomized partial QR algorithm.

NOVER must be positive or null and verifies the relationship:

NOVER + BLK_SIZE <= size(MAT,1)

and is adjusted if necessary to verify this relationship in all cases.

See Further Details and the cited references for the meaning and usefulness of the oversampling size in the partial randomized QR algorithm.

By default, the oversampling size is set to 10.

Further Details

The matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(krank), where krank <= min( m , n ).

Each H(i) has the form

H(i) = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).

On exit of PARTIAL_RQR_CMP_FIXED_PRECISION, the matrix Q can be computed explicitly by a call to subroutine ORTHO_GEN_QR with arguments MAT and BETA(:KRANK).

The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

On exit, if the optional argument TAU is absent:

• The elements above the diagonal of the array MAT(:krank,:krank) contain the corresponding elements of the triangular matrix R1.
• The elements of the diagonal of R1 are stored in the array DIAGR.
• The submatrix R2 is stored in MAT(:krank,krank+1:n).
• krank is stored in the real argument KRANK.

If TAU is present, a “partial” complete orthogonal factorization of MAT is computed. The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R (e.g., in R2), is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-krank) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R2.

The Z n-by-n orthogonal matrix is given by

Z = Z(1) * Z(2) * … * Z(krank)

On exit, if the optional argument TAU is present:

• The scalar tau defining T(k) is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,krank+1:n). The other elements of u(k) are not stored.
• The elements above the diagonal of the array section MAT(:krank,:krank) contain the corresponding elements of the triangular matrix T1. The elements of the diagonal of T1 are stored in the array DIAGR.
• krank is stored in the real argument KRANK.

In both cases, norm(MAT(KRANK+1:m,KRANK+1:n)) gives the error of the associated partial matrix approximation in the Frobenius norm, and argument RELERR stores the relative error in the Frobenius norm of the matrix approximation on exit.

The computations are parallelized if OPENMP is used. Note also that PARTIAL_RQR_CMP_FIXED_PRECISION uses a randomized “BLAS3” algorithm described in the references (3), (4), (5) and (6), which has about the same efficiency of a “BLAS3” QR algorithm without column pivoting and is thus highly efficient on large matrices.

For further details, see:

1. Lawson, C.L., and Hanson, R.J., 1974:

   Solving least square problems. Prentice-Hall.

2. Golub, G.H., and Van Loan, C.F., 1996:

   Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

3. Duersch, J.A., and Gu, M., 2017:

   Randomized QR with column pivoting. SIAM J. Sci. Comput., Volume 39, Issue 4, C263-C291.

4. Martinsson, P.G., Quintana-Orti, G., Heavner, N., and Van de Geijn, R., 2017:

   Householder QR factorization with randomization for column pivoting (HQRRP). SIAM J. Sci. Comput., Volume 39, Issue 2, C96-C115.

5. Xiao, J., Gu, M., and Langou, J., 2017:

   Fast parallel randomized QR with column pivoting algorithms for reliable low-rank matrix approximations. IEEE 24th International Conference on High Performance Computing (HiPC), IEEE, 2017, 233-242.

6. Duersch, J.A., and Gu, M., 2020:

   Randomized projection for rank-revealing matrix factorizations and low-rank approximations. SIAM Review, Volume 62, Issue 3, 661-682.

subroutine rqb_cmp ( mat, q, b, niter, rng_alg, ortho, comp_qr, ip, tol, tau )

Purpose

RQB_CMP computes a partial QB, QR (eventually with column pivoting) or complete orthogonal factorization of a full m-by-n real matrix MAT using randomized power or subspace iterations.

nqb is the target rank of the partial QB, QR or complete orthogonal decomposition, which is sought, and is equal to the number of columns of the output real matrix argument Q, i.e., nqb = size( Q, 2 ).

The routine first computes a partial QB factorization of MAT with the help of a randomized algorithm:

MAT = Q * B

, where Q is a m-by-nqb orthonormal matrix, B is a nqb-by-n matrix and the product Q*B is a good approximation of MAT according to the spectral or Frobenius norm.

In a second step:

• if the optional logical argument COMP_QR is used with the value true and the optional array arguments IP and TAU are absent, a QR factorization of B is computed to obtain an approximate QR factorization of MAT:

  MAT = Q * B = Q * ( O * R ) = (Q * O ) * R

  , where O is an nqb-by-nqb orthogonal matrix and R is an nqb-by-n upper trapezoidal matrix.

• if the optional array argument IP is present, a QR factorization of B with column pivoting is performed to obtain an approximate QR factorization with column pivoting of MAT :

  MAT * P = Q * ( O * R ) = (Q * O ) * R

  , where P is an n-by-n permutation matrix, O is an nqb-by-nqb orthogonal matrix and R is an nqb-by-n upper trapezoidal matrix.

• if the optional array argument TAU is present, a complete orthogonal factorization of B is performed to obtain an approximate complete orthogonal factorization of MAT:

  MAT = Q * ( O * T * Z ) = (Q * O ) * T * Z

  , where O is an nqb-by-nqb orthogonal matrix, Z is a n-by-n orthogonal matrix and T is a nqb-by-n matrix, which has the form:

  [ T1 0 ]

  , where T1 is a nqb-by-nqb upper triangular matrix.

• if, finally, both the optional array arguments IP and TAU are present, a complete orthogonal factorization of B with column pivoting is performed to obtain an approximate complete orthogonal factorization with column pivoting of MAT :

  MAT * P = Q * ( O * T * Z ) = (Q * O ) * T * Z

  , where P is a n-by-n permutation matrix, O is an nqb-by-nqb orthogonal matrix, Z is a n-by-n orthogonal matrix and T a nqb-by-n matrix, which has the form:

  [ T1 0 ]

  , where T1 is a nqb-by-nqb upper triangular matrix.

If a partial QR or complete orthogonal factorization is computed, we have the following partition of R:

R = [ R1 R2 ]

where R1 is a nqb-by-nqb upper triangular matrix and R2 is a full nqb-by-(n-nqb) matrix. Then, if:

• the optional scalar argument TOL is present and is in ]0,1[, the rank of R1 is determined by finding the submatrix of R1 which is defined as the largest leading submatrix whose estimated condition number, in the 1-norm, is less than 1/TOL. The order of this submatrix, krank, is the effective rank of R1.
• the optional scalar argument TOL is present and is equal to 0, the numerical rank of R1, krank, is determined.
• the optional scalar argument TOL is absent or outside [0,1[, the numerical rank of R1, krank, is determined by finding the largest leading submatrix of R1 such that abs(R1[j,j])>0.

In all cases, if krank < nqb, this indicates that column pivoting must be used or, if column pivoting has been already specified, that the rank of MAT is probably less than nqb and nqb = size( Q, 2 ) has been set to a too large value. In such cases, the subroutine will exit with an error message.

Arguments

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-n matrix MAT.

MAT is not modified by the routine.

Q (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed m-by-nqb orthonormal matrix of the partial QB, QR or complete orthogonal factorization of MAT.

See Further Details.

The shape of Q must verify:

• size( Q, 1 ) = m = size( MAT, 1 ),
• size( Q, 2 ) = nqb <= min( size(MAT,1) , size(MAT,2) ) .

B (OUTPUT) real(stnd), dimension(:,:)

On exit:

• the computed B matrix of the partial QB factorization of MAT if the optional arguments COMP_QR, IP and TAU are absent or if COMP_QR is used with the value false;
• the upper trapezoidal matrix R of the partial QR factorization if the optional logical argument COMP_QR is used with the value true or if the optional array argument IP is present;
• the computed matrices T1 and Z of the partial complete orthogonal factorization of MAT if the optional argument TAU is present.

See Further Details.

The shape of B must verify:

• size( B, 1 ) = size( Q, 2 ) = nqb ,
• size( B, 2 ) = size( MAT, 2 ) = n .

NITER (INPUT, OPTIONAL) integer(i4b)

The number of randomized power or subspace iterations performed in the first phase of the randomized algorithm for computing the preliminary QB factorization.

NITER must be positive or null.

By default, 5 randomized power or subspace iterations are performed.

RNG_ALG (INPUT, OPTIONAL) integer(i4b)

On entry, a scalar integer to select the random (uniform) number generator used to build the random gaussian test matrix in the randomized partial QB or QR algorithm.

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 current random number generator and its current state are not changed. Note further, that, on exit, the current random number generator is not reset to its previous value before the call to RQB_CMP.

See the documentation of subroutine RANDOM_SEED_ in module Random for further information.

ORTHO (INPUT, OPTIONAL) logical(lgl)

On entry, if:

• ORTHO=true, orthonormalization is carried out between each step of the power iterations, to avoid loss of accuracy due to rounding errors. This means that subspace iterations are used instead of power iterations in the QB phase of the algorithm,
• ORTHO=false, orthonormalization is not performed.

The default is to use orthonormalization, e.g., ORTHO=true.

COMP_QR (INPUT, OPTIONAL) logical(lgl)

The optional logical argument COMP_QR determines if a partial QB or QR factorization is computed.

On entry, if:

• COMP_QR=true, a partial QR factorization is computed;
• COMP_QR=false, a partial QB factorization is computed.

The default is to compute a partial QB factorization, e.g., COMP_QR=false.

IP (OUTPUT, OPTIONAL) integer(i4b), dimension(:)

On entry, if IP is present a partial QR factorization with column pivoting of MAT is performed instead of a QB factorization (e.g., this implies COMP_QR=true).

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

See Further Details.

The size of IP must verify: size( MAT, 2 ) = n.

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

If COMP_QR=true and TOL is present and is in [0,1[, then:

• the calculations to determine the condition number of of R1 in the 1-norm are performed. Then, TOL is used to determine the effective rank of R1, which is defined as the order of the largest leading triangular submatrix of R1, whose estimated condition number is less than 1/TOL.
• if TOL=0 is specified the numerical rank of R1 is determined.
• on exit, the reciprocal of the condition number is returned in TOL.

If COMP_QR=true, but TOL is not specified or is outside [0,1[ :

• the calculations to determine the condition number of R1 are not performed and crude tests on R1(j,j) are done to determine the rank of R1. If TOL is present, it is not changed.

If each element of MAT is correct to, say, 5 digits then TOL=0.00001 should be used on entry.

TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a complete orthogonal factorization of MAT is computed instead of a QB factorization (e.g., this implies COMP_QR=true).

On exit, the scalars factors of the elementary reflectors defining Z.

See Further Details.

The size of TAU must verify: size( TAU ) = nqb = size( Q, 2 ).

Further Details

For a good introduction to randomized linear algebra, see the references (1) and (2).

The randomized power or subspace iteration was proposed in (3; see Algorithm 4.4) to compute an orthonormal matrix whose range approximates the range of MAT. An approximate partial QB, QR or complete orthogonal factorization can then be computed using the aforementioned orthonormal matrix, see the references (1) and (3) for details.

The orthonormal m-by-nqb matrix Q of the partial QB, QR or complete orthogonal factorization of MAT is computed explicitly and stored on exit in the real array argument Q.

The nqb-by-n matrix B of the partial QB factorization or the upper trapezoidal matrix R of the partial QR or orthogonal factorization of MAT is stored on exit in the real array argument B.

If the integer array argument IP is present, a (partial) QR factorization of MAT with column pivoting is computed, as described above, and on exit the permutation matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

Next, if the real array argument TAU is present, a (partial) complete orthogonal factorization of MAT is computed from the partial QR factorization of MAT. The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R of the (partial) QR factorization of MAT, is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-nqb) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.

The Z n-by-n orthogonal matrix is finally given by the product

Z = Z(1) * Z(2) * … * Z(nqb)

On exit:

• the orthogonal matrix Z is stored in factored form. The scalar tau, which defines Z(k) is returned in the kth element of TAU and the vector u(k) in the kth row of B, such that the elements of z(k) are in B(k,nqb+1:n).
• the upper triangular nqb-by-nqb matrix T1, which defines the T factor of the (partial) complete orthogonal factorization of MAT (see above) is stored in the upper triangle of the real array argument B.

Finally, if the optional arguments IP and TAU are both present, a (partial) complete orthogonal factorization with column pivoting of MAT is computed and on exit this factorization is stored with the same conventions as described above.

For further details on randomized linear algebra, computing low-rank matrix approximations using randomized power or subspace iterations, see:

1. Martinsson, P.G., 2019:

   Randomized methods for matrix computations. arXiv.1607.01649

2. Erichson, N.B., Voronin, S., Brunton, S.L., and Kutz, J.N., 2019:

   Randomized matrix decompositions using R. arXiv.1608.02148

3. Halko, N., Martinsson, P.G., and Tropp, J.A., 2011:

   Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev., 53, 217-288.

4. Gu, M., 2015:

   Subspace iteration randomization and singular value problems. SIAM J. Sci. Comput., 37, A1139-A1173.

subroutine rqb_cmp_fixed_precision ( mat, relerr, q, b, failure_relerr, niter, rng_alg, blk_size, maxiter_qb, ortho, reortho, niter_qb, comp_qr, ip, tol, tau )

Purpose

RQB_CMP_FIXED_PRECISION computes a partial QB, QR (eventually with column pivoting) or complete orthogonal factorization of a full m-by-n real matrix MAT using randomized power or subspace iterations.

nqb is the target rank of the partial QB, QR or complete orthogonal factorization, which is sought, and this partial factorization must have an approximation error which fulfills:

|| MAT - Q*B ||_F <= ||MAT||_F * relerr

, where Q is a m-by-nqb orthonormal matrix, B is a nqb-by-n matrix and the matrix product Q*B is the computed matrix approximation. || ||_F is the Frobenius norm and relerr is a prescribed accuracy tolerance for the relative error of the computed matrix approximation, specified in the input argument RELERR.

In other words, the rank, nqb, of the matrix approximation is not known in advance and is determined in the subroutine. This explains why the output array arguments Q and B, which contain the computed factors of the matrix approximation, must be declared in the calling program as pointers.

On exit, nqb is equal to the numbers of columns of the output array pointer argument Q, which is also equal to the numbers of rows of the output array pointer argument B. In other words, nqb = size( Q, 2 ) = size( B, 1 ) and the relative error, in the Frobenius norm, of the computed matrix approximation Q * B is output in argument RELERR.

RQB_CMP_FIXED_PRECISION searches incrementally the best (e.g., of smallest rank) Q * B approximation, which fulfills the prescribed accuracy tolerance for the relative error. More precisely, the rank of the matrix approximation is increased progressively of BLK_SIZE by BLK_SIZE until the prescribed accuracy tolerance is satisfied and than adjusted precisely to obtain the Q * B matrix approximation of smallest rank, which satisfies the prescribed tolerance.

Note that the product of the two integer arguments BLK_SIZE and MAXITER_QB (see below for their precise meaning), BLK_SIZE*MAXITER_QB, determines the maximum allowable rank of the matrix approximation, which is sought. In other words, the subroutine will stop the search for the best (e.g., smallest) matrix approximation, which fulfills the requested tolerance, if the rank of this matrix approximation exceeds BLK_SIZE*MAXITER_QB. In that case, the subroutine will return the current matrix approximation (with a rank equal to BLK_SIZE*MAXITER_QB).

In all cases the relative error of the computed matrix approximation is output in argument RELERR.

If, finally, the optional logical argument FAILURE_RELERR is used, it will be set to true if the computed matrix approximation does not fulfill the requested relative error specified on entry in the argument RELERR and to false otherwise.

In a second step:

• if the optional logical argument COMP_QR is used with the value true and the optional array arguments IP and TAU are absent, a QR factorization of B is computed to obtain an approximate QR factorization of MAT:

  MAT = Q * B = Q * ( O * R ) = (Q * O ) * R

  , where O is an nqb-by-nqb orthogonal matrix and R is an nqb-by-n upper trapezoidal matrix.

• if the optional array argument IP is present, a QR factorization of B with column pivoting is performed to obtain an approximate QR factorization with column pivoting of MAT :

  MAT * P = Q * ( O * R ) = (Q * O ) * R

  , where P is an n-by-n permutation matrix, O is an nqb-by-nqb orthogonal matrix and R is an nqb-by-n upper trapezoidal matrix.

• if the optional array argument TAU is present, a complete orthogonal factorization of B is performed to obtain an approximate complete orthogonal factorization of MAT:

  MAT = Q * ( O * T * Z ) = (Q * O ) * T * Z

  , where O is an nqb-by-nqb orthogonal matrix, Z is a n-by-n orthogonal matrix and T is a nqb-by-n matrix, which has the form:

  [ T1 0 ]

  , where T1 is a nqb-by-nqb upper triangular matrix.

• if, finally, both the optional array arguments IP and TAU are present, a complete orthogonal factorization of B with column pivoting is performed to obtain an approximate complete orthogonal factorization with column pivoting of MAT :

  MAT * P = Q * ( O * T * Z ) = (Q * O ) * T * Z

  , where P is a n-by-n permutation matrix, O is an nqb-by-nqb orthogonal matrix, Z is a n-by-n orthogonal matrix and T a nqb-by-n matrix, which has the form:

  [ T1 0 ]

  , where T1 is a nqb-by-nqb upper triangular matrix.

If a partial QR or complete orthogonal factorization is computed, we have the following partition of R:

R = [ R1 R2 ]

where R1 is a nqb-by-nqb upper triangular matrix and R2 is a full nqb-by-(n-nqb) matrix. Then, if:

• the optional scalar argument TOL is present and is in ]0,1[, the rank of R1 is determined by finding the submatrix of R1 which is defined as the largest leading submatrix whose estimated condition number, in the 1-norm, is less than 1/TOL. The order of this submatrix, krank, is the effective rank of R1.
• the optional scalar argument TOL is present and is equal to 0, the numerical rank of R1, krank, is determined.
• the optional scalar argument TOL is absent or outside [0,1[, the numerical rank of R1, krank, is determined by finding the largest leading submatrix of R1 such that abs(R1[j,j])>0.

In all cases, if krank < nqb, this indicates that column pivoting must be used or, if column pivoting has been already specified, that the rank of MAT is probably less than nqb and nqb = size( Q, 2 ) has been set to a too large value. In such cases, the subroutine will exit with an error message.

Arguments

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-n matrix MAT.

MAT is not modified by the routine.

RELERR (INPUT/OUTPUT) real(stnd)

On entry, the requested accuracy tolerance for the relative error of the computed matrix approximation.

The preset value for RELERR must be greater than 4*epsilon( RELERR ), less than one and verifies:

• RELERR >= 2 * sqrt( epsilon( RELERR )/RELERR )

and is forced to be greater than 2*sqrt( epsilon( RELERR )/RELERR ) if this is not the case to avoid loss of accuracy in the algorithm. See reference (6) for more details.

On exit, RELERR contains the relative error of the computed matrix approximation:

• RELERR = ||MAT-Q*B||_F / ||MAT||_F

Q (OUTPUT) real(stnd), dimension(:,:), pointer

On exit, the computed m-by-nqb orthonormal matrix of the partial QB, QR or complete orthogonal factorization of MAT.

The statut of the pointer Q must not be undefined on entry. If, on entry, the pointer Q is already allocated, it will be first deallocated and then reallocated with the correct shape.

On exit, the shape of the pointer Q will verify:

• size( Q, 1 ) = m = size( MAT, 1 ),
• size( Q, 2 ) = nqb <= min( size(MAT,1) , size(MAT,2) ) .

B (OUTPUT) real(stnd), dimension(:,:), pointer

On exit:

• the computed B matrix of the partial QB factorization of MAT if the optional arguments COMP_QR, IP and TAU are absent or if COMP_QR is used with the value false;
• the upper trapezoidal matrix R of the partial QR factorization if the optional logical argument COMP_QR is used with the value true or if the optional array argument IP is present;
• the computed matrices T1 and Z of the partial complete orthogonal factorization of MAT if the optional argument TAU is present.

The statut of the pointer B must not be undefined on entry. If, on entry, the pointer B is already allocated, it will be first deallocated and then reallocated with the correct shape.

On exit, the shape of the pointer B will verify:

• size( B, 1 ) = size( Q, 2 ) = nqb ,
• size( B, 2 ) = size( MAT, 2 ) = n .

FAILURE_RELERR (OUTPUT, OPTIONAL) logical(lgl)

On exit, if the optional logical argument FAILURE_RELERR is present, it is set as follows:

• FAILURE_RELERR = false : indicates successful exit and the computed matrix approximation fulfills the requested relative error specified on entry in the argument RELERR,
• FAILURE_RELERR = true : indicates that the computed matrix approximation has a relative error larger than the requested relative error. This means that the requested accuracy tolerance for the relative error is too small (i.e., RELERR < 2 * sqrt( epsilon( RELERR )/RELERR ) or that the input parameters BLK_SIZE and/or MAXITER_QB have a too small value, given the distribution of the singular values of MAT, and must be increased to fullfill the preset accuracy tolerance for the relative error of the matrix approximation.

NITER (INPUT, OPTIONAL) integer(i4b)

The number of randomized power or subspace iterations performed in the first phase of the randomized algorithm for computing the preliminary QB factorization.

NITER must be positive or null.

By default, 1 randomized power or subspace iteration is performed.

RNG_ALG (INPUT, OPTIONAL) integer(i4b)

On entry, a scalar integer to select the random (uniform) number generator used to build the random gaussian test matrix in the randomized partial QB or QR algorithm.

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 current random number generator and its current state are not changed. Note further, that, on exit, the current random number generator is not reset to its previous value before the call to RQB_CMP_FIXED_PRECISION.

See the documentation of subroutine RANDOM_SEED_ in module Random for further information.

BLK_SIZE (INPUT, OPTIONAL) integer(i4b)

On entry, the block size used in the randomized QB factorization, which is used in the first phase of the subroutine.

BLK_SIZE must be greater or equal to one and less than min(m,n) and must be set to a much smaller value than min(m,n) usually, depending also on the architecture of the computer.

By default, BLK_SIZE is set to min( 10, min(m,n) ).

MAXITER_QB (INPUT, OPTIONAL) integer(i4b)

MAXITER_QB controls the maximum number of allowed iterations in the randomized QB algorithm, which is used in the first phase of the subroutine.

MAXITER_QB must be set greater or equal to one and less than int( min(m,n)/BLK_SIZE ).

By default, MAXITER_QB is set to max( 1, int( min(m,n)/(4*BLK_SIZE) ) ).

ORTHO (INPUT, OPTIONAL) logical(lgl)

On entry, if:

• ORTHO=true, orthonormalization is carried out between each step of the power iterations, to avoid loss of accuracy due to rounding errors. This means that subspace iterations are used instead of power iterations in the QB phase of the algorithm,
• ORTHO=false, orthonormalization is not performed.

The default is to use orthonormalization, e.g., ORTHO=true.

REORTHO (INPUT, OPTIONAL) logical(lgl)

On entry, if:

• REORTHO=true, a reorthogonalization step is performed to avoid the loss of orthogonality in the Gram-Schmidt procedure, which is used in the randomized QB factorization;
• REORTHO=false, a reorthogonalization step is not performed in the Gram-Schmidt procedure.

The default is to use a reorthogonalization step, e.g., REORTHO=true.

NITER_QB (INPUT, OPTIONAL) integer(i4b)

The number of subspace iterations performed in the last phase of the QB algorithm for improving the initial QB factorization of MAT.

NITER_QB must be greater or equal to 0.

By default, 2 final subspace iterations are performed.

COMP_QR (INPUT, OPTIONAL) logical(lgl)

The optional logical argument COMP_QR determines if a partial QB or QR factorization is computed.

On entry, if:

• COMP_QR=true, a partial QR factorization is computed;
• COMP_QR=false, a partial QB factorization is computed.

The default is to compute a partial QB factorization, e.g., COMP_QR=false.

IP (OUTPUT, OPTIONAL) integer(i4b), dimension(:)

On entry, if IP is present a partial QR factorization with column pivoting of MAT is performed instead of a QB factorization (e.g., this implies COMP_QR=true).

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

See Further Details.

The size of IP must verify: size( MAT, 2 ) = n.

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

If COMP_QR=true and TOL is present and is in [0,1[, then:

• the calculations to determine the condition number of of R1 in the 1-norm are performed. Then, TOL is used to determine the effective rank of R1, which is defined as the order of the largest leading triangular submatrix of R1, whose estimated condition number is less than 1/TOL.
• if TOL=0 is specified the numerical rank of R1 is determined.
• on exit, the reciprocal of the condition number is returned in TOL.

If COMP_QR=true, but TOL is not specified or is outside [0,1[ :

• the calculations to determine the condition number of R1 are not performed and crude tests on R1(j,j) are done to determine the rank of R1. If TOL is present, it is not changed.

If each element of MAT is correct to, say, 5 digits then TOL=0.00001 should be used on entry.

TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:), pointer

On entry, if TAU is present, a complete orthogonal factorization of MAT is computed instead of a QB factorization (e.g., this implies COMP_QR=true).

On exit, the scalars factors of the elementary reflectors defining Z.

See Further Details.

The statut of the pointer TAU must not be undefined on entry. If, on entry, the pointer TAU is already allocated, it will be first deallocated and then reallocated with the correct size.

On exit, the size of the pointer TAU will verify:

• size( TAU ) = nqb = size( Q, 2 ) = size( B, 1 ).

Further Details

For a good introduction to randomized linear algebra , see the references (1) and (2).

The randomized subspace iteration was proposed in (3; see Algorithm 4.4) to compute an orthonormal matrix whose range approximates the range of MAT. An approximate partial QB, QR or complete orthogonal factorization can then be computed using the aforementioned orthonormal matrix, see the references (1) and (3) for details.

Usually, the problem of low-rank matrix approximation falls into two categories:

• the fixed-rank problem, where the rank parameter nqb is given;

• the fixed-precision problem, where we seek a partial matrix factorization, Q * B, of rank as small as possible such that

  || MAT - Q*B ||_F <= eps

  , where eps is a given accuracy tolerance.

RQB_CMP_FIXED_PRECISION is dedicated to solve the fixed-precision problem. The fixed-rank problem can be solved by subroutine RQB_CMP.

RQB_CMP_FIXED_PRECISION uses an improved version of the “randQB_FP” algorithm described in the reference (6) to solve the fixed-precision problem.

The orthonormal m-by-nqb matrix Q of the partial QB, QR or complete orthogonal factorization of MAT is computed explicitly and stored on exit in the real array pointer Q.

The nqb-by-n matrix B of the partial QB factorization or the upper trapezoidal matrix R of the partial QR or orthogonal factorization of MAT is stored on exit in the real array pointer B.

If the integer array argument IP is present, a (partial) QR factorization of MAT with column pivoting is computed, as described above, and on exit the permutation matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

Next, if the real array pointer TAU is present, a (partial) complete orthogonal factorization of MAT is computed from the partial QR factorization of MAT. The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R of the (partial) QR factorization of MAT, is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-nqb) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.

The Z n-by-n orthogonal matrix is finally given by the product

Z = Z(1) * Z(2) * … * Z(nqb)

On exit:

• the orthogonal matrix Z is stored in factored form. The scalar tau, which defines Z(k) is returned in the kth element of TAU and the vector u(k) in the kth row of B, such that the elements of z(k) are in B(k,nqb+1:n).
• the upper triangular nqb-by-nqb matrix T1, which defines the T factor of the (partial) complete orthogonal factorization of MAT (see above) is stored in the upper triangle of the real array argument B.

Finally, if the optional arguments IP and TAU are both present, a (partial) complete orthogonal factorization with column pivoting of MAT is computed and on exit this factorization is stored with the same conventions as described above.

For further details, on randomized linear algebra, computing low-rank matrix approximations using randomized power or subspace iterations or solving the fixed-precision problem, see:

1. Martinsson, P.G., 2019:

   Randomized methods for matrix computations. arXiv.1607.01649

2. Erichson, N.B., Voronin, S., Brunton, S.L., and Kutz, J.N., 2019:

   Randomized matrix decompositions using R. arXiv.1608.02148

3. Halko, N., Martinsson, P.G., and Tropp, J.A., 2011:

   Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev., 53, 217-288.

4. GU, M., 2015:

   Subspace iteration randomization and singular value problems. SIAM J. Sci. Comput., 37, A1139-A1173.

5. Martinsson, P.-G., and Voronin, S., 2016:

   A randomized blocked algorithm for efficiently computing rank-revealing factorizations of matrices. SIAM J. Sci. Comput., 38:5, S485-S507.

6. Yu, W., Gu, Y., and Li, Y., 2018:

   Efficient randomized algorithms for the fixed-precision low-rank matrix approximation. SIAM J. Mat. Ana. Appl., 39:3, 1339-1359.

subroutine id_cmp ( mat, ip, t, c, v, diagr, beta, rnorm, tol, random_qr, rng_alg, blk_size, nover )

Purpose

ID_CMP computes a (partial) column Interpolative Decomposition (ID) of a m-by-n real matrix MAT.

A column ID factorization of rank krank approximates MAT as:

MAT = C * V

where C is an m-by-krank matrix, which consists of a subset of krank columns of MAT and V is a krank-by-n matrix, which contains a krank-by-krank identity matrix as a submatrix.

Such column ID factorization can be computed with the help of a (randomized) (partial) QR factorization with column pivoting of MAT (see the references (1), (2) and (4)-(7) for details), which is defined as

MAT * P = Q * R

, where P is a n-by-n permutation matrix, R is an upper triangular or trapezoidal (i.e., if n>m) matrix and Q is a m-by-m orthogonal matrix if the QR factorization is complete.

For computing a column ID decomposition of rank, krank, a partial QR factorization with column pivoting of MAT of rank, krank, is sufficient, e.g., the QR decomposition can be stopped when the numbers of computed columns of Q is equal to krank.

This leads implicitly to the following partition of Q:

Q = [ Q1 Q2 ]

where Q1 is an m-by-krank orthonormal matrix and Q2 is m-by-(m-krank) orthonormal matrix orthogonal to Q1, and to the following corresponding partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R11 is a krank-by-krank triangular matrix, R21 is zero by construction, R12 is a full krank-by-(n-krank) matrix and R22 is a full (m-krank)-by-(n-krank) matrix. This leads to the following approximation of MAT:

MAT = Q1 * [ R11 R12 ] * P’

Here, we see that Q1 * R11 equals the first krank columns of A * P, and so we can define the matrix C in the partial ID factorization of MAT as

C = Q1 * R11 = MAT * P(:,:krank)

, then the dominant term Q1 * [ R11 R12 ] in the partial QR decomposition of MAT * P can be written as

Q1 * [ R11 R12 ] = Q1 * R11 * [ I T ] = C * [ I T ]

where I is the identity matrix of order krank and T a krank-by-(n-krank) matrix and a solution to the matrix equation:

R11 * T = R12 ,e.g., T = inv(R11) * R12

Thus, we can obtain a column ID factorization of MAT from its partial QR factorization with column pivoting as

MAT = (Q1 * R11) * [ I T ] * P’ = C * V

where V = [ I T ] * P’ is a krank-by-n matrix.

Finally, observe that the error matrix associated with the column ID decomposition of MAT is given by

MAT - C * V = [ 0 Q2*R22] * P’

and that the Frobenius norm of this error matrix can be computed efficiently as

|| MAT - C * V ||_F = || Q2 * R22 ||_F = || R22 ||_F

krank is the target rank of the column ID, which is sought, and is set to the number of rows of the output array argument T.

If the optional logical argument RANDOM_QR is used with the value true, a fast randomized partial QR factorization with column pivoting is used in the first phase of the ID algorithm.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-n matrix MAT.

On exit, MAT is overwritten by details of its partial QR factorization with column pivoting.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures or of PARTIAL_RQR_CMP subroutine in module Random for further details on how the partial QR decomposition is stored compactly in MAT (and arguments IP, DIAGR and BETA) on exit.

Note that, on exit, the Frobenius norm of the error associated with the computed column ID is equal to norm( mat(krank+1:m,krank+1:n) ), which is the same as the Frobenius norm of the error associated with the partial QR decomposition with column pivoting (of rank krank) of MAT.

IP (OUTPUT) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

The matrix C in the (column) ID of MAT corresponds to the subset of the columns of MAT with the indices IP(:krank).

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP subroutine in module Random for further details.

The size of IP must be equal to size( MAT, 2 ) = n

T (OUTPUT) real(stnd), dimension(:,:)

On exit, the krank-by-(n-krank) submatrix T = inv(R11) * R12 in the column ID factorization of MAT.

The shape of T must verify:

• size( T, 1 ) = krank <= min( size(MAT,1) , size(MAT,2) ) = min(m,n)
• size( T, 2 ) = size( MAT, 2 ) - size( T, 1 ) = n - krank

C (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, the m-by-krank matrix C in the ID factorization of MAT. Note that C is computed as C = Q * R11 as MAT is overwritten by its partial QR factorization with column pivoting before C can be estimated.

The shape of C must verify:

• size( C, 1 ) = size( MAT, 1 ) = m
• size( C, 2 ) = size( T, 1 ) = krank

V (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, the krank-by-n matrix V in the ID factorization of MAT.

The shape of V must verify:

• size( V, 1 ) = size( T, 1 ) = krank
• size( V, 2 ) = size( MAT, 2 ) = n

RNORM (OUTPUT, OPTIONAL) real(stnd)

On exit, the Frobenius norm of the error matrix associated with the column ID computed as || R22 ||_F .

DIAGR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R11 in the partial QR factorization with column pivoting of MAT.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP subroutine in module Random for further details.

The size of DIAGR must verify:

• size( DIAGR ) = krank <= min( size(MAT,1) , size(MAT,2) )

BETA (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q in the partial QR factorization with column pivoting of MAT.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP subroutine in module Random for further details.

The size of BETA must verify:

• size( BETA ) = krank <= min( size(MAT,1) , size(MAT,2) ).

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, if TOL is present then:

• if TOL is in ]0,1[, calculations to determine the condition number of submatrix R11 in the partial QR factorization of MAT are performed. TOL is used to determine the effective pseudo-rank of R11, which is defined as the order of the largest leading triangular submatrix in the partial QR factorization with column pivoting of MAT, whose estimated condition number in the 1-norm is less than 1/TOL. On exit, the reciprocal of the condition number is returned in TOL.
• if TOL=0 is specified, the calculations to determine the condition number of R11 are not performed and crude tests on R(j,j) are done to determine the numerical pseudo-rank of R11. On exit, TOL is not changed.

The TOL argument is useful to check if the matrix R11 in the QR decomposition of MAT, which must be inverted to compute T (in the column ID of MAT), is sufficiently well conditioned to obtain a stable and robust column ID of MAT.

If the computed pseudo-rank of R11 is less than krank = size( T, 1 ), the subroutine will exit with an error message as a stable solution to the matrix equation R11 * T = R12 cannot be computed.

On the other hand, if TOL is not specified or is outside [0,1[, the calculations to determine the rank of R11 are not performed and this rank is assumed to be equal to krank = size( T, 1 ).

If it is possible that MAT may not be of full rank then the stability of the column ID decomposition of MAT, which is sought, can be checked by using TOL=relative precision of the elements in MAT. If each element of MAT is correct to, say, 5 digits then TOL=0.00001 should be used.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP subroutine in module Random for further details.

RANDOM_QR (INPUT, OPTIONAL) logical(lgl)

On entry, if RANDOM_QR is used with the value true, a fast randomized partial QR factorization with column pivoting is used in the first phase of the ID algorithm.

By default, RANDOM_QR = false, i.e., a standard (partial) QR factorization with column pivoting is used in the first phase of the ID algorithm.

RNG_ALG (INPUT, OPTIONAL) integer(i4b)

On entry, a scalar integer to select the random (uniform) number generator used to build the random gaussian matrix in the randomized (partial) QR phase of the ID algorithm if RANDOM_QR = true.

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 current random number generator and its current state are not changed. Note further, that, on exit, the current random number generator is not reset to its previous value before the call to ID_CMP subroutine.

See the documentation of subroutine RANDOM_SEED_ in module Random for further information.

This optional argument has no effect if logical argument RANDOM_QR is set to false or is absent in the call to ID_CMP subroutine.

BLK_SIZE (INPUT, OPTIONAL) integer(i4b)

On entry, the block size used in the randomized partial QR phase of the ID algorithm if RANDOM_QR = true.

BLK_SIZE must be greater or equal to one and less than min(m,n) and must be set to a much smaller value than min(m,n) usually, depending also on the architecture of the computer.

By default, BLK_SIZE is set to min( BLKSZ_QR, min(m,n) ), where parameter BLKSZ_QR is the default block size for QR related algorithms specified in module Select_Parameters.

See description of subroutine PARTIAL_RQR_CMP in module Random for the meaning of the the block size in the randomized (partial) QR algorithm.

This optional argument has no effect if logical argument RANDOM_QR is set to false or is absent in the call to ID_CMP.

NOVER (INPUT, OPTIONAL) integer(i4b)

The oversampling size used in the randomized partial QR phase of the ID algorithm if RANDOM_QR = true.

NOVER must be positive or null and verifies the relationship:

NOVER + BLK_SIZE <= size(MAT,1)

and is adjusted if necessary to verify this relationship in all cases.

See description of subroutine PARTIAL_RQR_CMP in module Random for the meaning and usefulness of the oversampling size in the randomized (partial) QR algorithm.

By default, the oversampling size is set to 10.

This optional argument has no effect if logical argument RANDOM_QR is set to false or is absent in the call to ID_CMP.

Further Details

For further details on the ID decomposition and computing a column ID from a (randomized) partial QR factorization with column pivoting, see:

1. Martinsson, P.G., 2019:

   Randomized methods for matrix computations. arXiv.1607.01649

2. Erichson, N.B., Voronin, S., Brunton, S.L., and Kutz, J.N., 2019:

   Randomized matrix decompositions using R. arXiv.1608.02148

3. Voronin, S., Martinsson, P.G., 2015:

   Rsvdpack: Subroutines for computing partial singular value decompositions via randomized sampling on single core, multi core, and gpu architectures. arXiv.1502.05366

4. Voronin, S., Martinsson, P.G., 2017:

   Efficient algorithms for cur and interpolative matrix decompositions Adv Comput Math, Volume 43, 495-516.

5. Xiao, J., Gu, M., and Langou, J., 2017:

   Fast parallel randomized QR with column pivoting algorithms for reliable low-rank matrix approximations. IEEE 24th International Conference on High Performance Computing (HiPC), IEEE, 2017, 233-242.

6. Stewart, G.W., 1999:

   Four algorithms for the the efficient computation of truncated pivoted qr approximations to a sparse matrix. Numerische Mathematik, Volume 83, 313-323.

7. Berry, M.W., Pulatova, S.A., and Stewart, G.W., 2005:

   Algorithm 844: Computing sparse reduced-rank approximations to sparse matrices. ACM Transactions on Mathematical Software, Volume 31, No. 2, 252-269.

subroutine ts_id_cmp ( mat, ip_row, ip_col, w, v, skelmat, diagr, beta, rnorm, tol, random_qr, rng_alg, blk_size, nover )

Purpose

TS_ID_CMP computes a (partial) two-sided Interpolative Decomposition (tsID) of a full m-by-n real matrix MAT.

A tsID factorization of rank krank approximates MAT as the matrix product

MAT = W * MAT_skel * V

where W is a m-by-krank matrix, V is a krank-by-n matrix and MAT_skel consists of a squared krank-by-krank submatrix of MAT, which defines the so-called skeleton of MAT.

The tsID factorization can be computed with the help of (randomized) (partial) QR factorizations with column pivoting of MAT and of a matrix derived from its partial QR decomposition, more precisely with a column ID of MAT and a row ID of a subset of the columns of MAT.

The first step is thus to compute a partial column ID decomposition of MAT as

MAT = C * V

where C is a m-by-krank subset of the columns of MAT and V is a krank-by-n matrix. See description of subroutine ID_CMP for more details about the ID decomposition and how such decomposition can be computed from a QR decomposition with column pivoting of MAT.

In a second step, a complete column ID decomposition of C’ (e.g., a row ID decomposition of C) is computed as

C’ = MAT_skel’ * W’

where MAT_skel is a squared krank-by-krank matrix, which is a submatrix of MAT, and W is a m-by-krank matrix. This gives the desired tsID decomposition of MAT as

MAT = W * MAT_skel * V

See the references (1) and (4) and also description of ID_CMP subroutine in module Random for more details on the ID and tsID decompositions.

krank is the target rank of the tsID decomposition, which is sought, and is equal to the number of rows or columns of the output array argument SKELMAT, which stores the skeleton of MAT, e.g., the matrix MAT_skel.

If the optional logical argument RANDOM_QR is used with the value true, fast randomized (partial) QR factorizations with column pivoting are used for computing the column and row ID decompositions in the two phases of the tsID algorithm.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-n matrix MAT.

On exit, MAT is overwritten by details of its partial QR factorization with column pivoting.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures or of PARTIAL_RQR_CMP subroutine in module Random for further details on how the partial QR decomposition is stored compactly in MAT (and arguments IP, DIAGR and BETA) on exit.

Note that, the Frobenius norm of the error associated with the computed tsID is the same as the one associated with this partial QR factorization with column pivoting of MAT or its column ID decomposition and is given by norm( mat(krank+1:m,krank+1:n) ) on exit. See description of ID_CMP subroutine for further details.

IP_ROW (OUTPUT) integer(i4b), dimension(:)

On exit, IP_ROW stores the permutation matrix N in the partial LQ decomposition with row pivoting of C, which is computed to obtain the row ID decomposition of C.

On exit, if IP_ROW(j)=k, then the j-th row of N*C was the k-th row of C, where N = I(IP_ROW,:) and I is the identity matrix of order m.

The matrix MAT_skel in the tsID decomposition of MAT is the submatrix defined as the intersection of the rows IP_ROW(:krank) and the columns IP_COL(:krank) of MAT.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP and ID_CMP subroutines in module Random for further details.

The size of IP_ROW must be equal to size( MAT, 1 ) = m.

IP_COL (OUTPUT) integer(i4b), dimension(:)

On exit, IP_COL stores the permutation matrix P in the partial QR decomposition with column pivoting of MAT, which is computed to obtain the column ID decomposition of MAT.

If IP_COL(j)=k, then the j-th column of MAT*P was the k-th column of MAT, where P = I(:,IP_COL) and I is the identity matrix of order n.

The matrix C in the (column) ID of MAT corresponds to the subset of the columns of MAT with the indices IP_COL(:krank). Furthermore, the matrix MAT_skel in the tsID decomposition of MAT is the submatrix defined as the intersection of the rows IP_ROW(:krank) and the columns IP_COL(:krank) of MAT.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP and ID_CMP subroutines in module Random for further details.

The size of IP_COL must be equal to size( MAT, 2 ) = n.

W (OUTPUT) real(stnd), dimension(:,:)

On exit, the m-by-krank matrix W in the tsID factorization of MAT.

The shape of W must verify:

• size( W, 1 ) = size( MAT, 1 ) = m
• size( W, 2 ) = size( SKELMAT, 1 ) = krank

V (OUTPUT) real(stnd), dimension(:,:)

On exit, the krank-by-n matrix V in the tsID factorization of MAT.

The shape of V must verify:

• size( V, 1 ) = size( SKELMAT, 1 ) = krank
• size( V, 2 ) = size( MAT, 2 ) = n

SKELMAT (OUTPUT) real(stnd), dimension(:,:)

On exit, the squared krank-by-krank matrix MAT_skel (e.g., the skeleton of MAT) in the tsID factorization of MAT.

MAT_skel is the submatrix of MAT defines by the rows IP_ROW(:krank) and the columns IP_COL(:krank) of MAT on entry. However, in the routine, MAT_skel is recomputed from the column ID of MAT and row ID of C to save space as MAT is overwritten by details of its partial QR factorization when computing the column ID of MAT.

The shape of SKELMAT must verify:

• size( SKELMAT, 1 ) = size( SKELMAT, 2 ) = krank <= min(m,n)

RNORM (OUTPUT, OPTIONAL) real(stnd)

On exit, the Frobenius norm of the error matrix associated with the tsID, which is defined as || MAT - W * MAT_skel * V ||_F . Note that this error matrix is the same as the one associated with the partial QR decomposition with column pivoting of MAT or its column ID decomposition.

DIAGR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R11 in the partial QR factorization with column pivoting of MAT.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP subroutine in module Random for further details.

The size of DIAGR must verify:

• size( DIAGR ) = krank <= min( size(MAT,1) , size(MAT,2) ) = min(m,n)

BETA (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q in the partial QR factorization with column pivoting of MAT.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP subroutine in module Random for further details.

The size of BETA must verify:

• size( BETA ) = krank <= min( size(MAT,1) , size(MAT,2) ) = min(m,n)

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, if TOL is present then:

• if TOL is in ]0,1[, calculations to determine the condition number of submatrix R11 in the partial QR factorization of MAT are performed. TOL is used to determine the effective pseudo-rank of R11, which is defined as the order of the largest leading triangular submatrix in the partial QR factorization with column pivoting of MAT, whose estimated condition number in the 1-norm is less than 1/TOL. On exit, the reciprocal of the condition number is returned in TOL.
• if TOL=0 is specified, the calculations to determine the condition number of R11 are not performed and crude tests on R(j,j) are done to determine the numerical pseudo-rank of R11. On exit, TOL is not changed.

The TOL argument is useful to check if the matrix R11 in the QR decomposition of MAT, which must be inverted to compute T (in the column ID of MAT), is sufficiently well conditioned to obtain a stable and robust column ID of MAT and, consequently, a robust tsID decomposition.

If the computed pseudo-rank of R11 is less than krank = size( SKELMAT, 1 ), the subroutine will exit with an error message.

On the other hand, if TOL is not specified or is outside [0,1[, the calculations to determine the rank of R11 are not performed and this rank is assumed to be equal to krank = size( SKELMAT, 1 ).

If it is possible that MAT may not be of full rank then the stability of the tsID decomposition of MAT, which is sought, can be checked by using TOL=relative precision of the elements in MAT. If each element of MAT is correct to, say, 5 digits then TOL=0.00001 should be used.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP subroutine in module Random for further details.

RANDOM_QR (INPUT, OPTIONAL) logical(lgl)

On entry, if RANDOM_QR is used with the value true, fast randomized (partial) QR factorizations are used in the two phases of the tsID algorithm.

By default, RANDOM_QR = false, i.e., a standard (partial) QR factorization with column pivoting is used in the two phases of the tsID algorithm.

RNG_ALG (INPUT, OPTIONAL) integer(i4b)

On entry, a scalar integer to select the random (uniform) number generator used to build the random gaussian matrix in the two randomized (partial) QR phases of the tsID algorithm if RANDOM_QR = true.

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 current random number generator and its current state are not changed. Note further, that, on exit, the current random number generator is not reset to its previous value before the call to TS_ID_CMP subroutine.

See the documentation of subroutine RANDOM_SEED_ in module Random for further information.

This optional argument has no effect if logical argument RANDOM_QR is set to false or is absent in the call to TS_ID_CMP subroutine.

BLK_SIZE (INPUT, OPTIONAL) integer(i4b)

On entry, the block size used in the randomized partial QR phases of the tsID algorithm if RANDOM_QR = true.

BLK_SIZE must be greater or equal to one and less than min(m,n) and must be set to a much smaller value than min(m,n) usually, depending also on the architecture of the computer.

By default, BLK_SIZE is set to min( BLKSZ_QR, min(m,n) ), where parameter BLKSZ_QR is the default block size for QR related algorithms specified in module Select_Parameters.

See description of subroutine PARTIAL_RQR_CMP in module Random for the meaning of the the block size in the randomized (partial) QR algorithm.

This optional argument has no effect if logical argument RANDOM_QR is set to false or is absent in the call to TS_ID_CMP subroutine.

NOVER (INPUT, OPTIONAL) integer(i4b)

The oversampling size used in the randomized partial QR phases of the tsID algorithm if RANDOM_QR = true.

NOVER must be positive or null and verifies the relationship:

NOVER + BLK_SIZE <= size(MAT,1)

and is adjusted if necessary to verify this relationship in all cases.

See description of subroutine PARTIAL_RQR_CMP in module Random for the meaning and usefulness of the oversampling size in the randomized (partial) QR algorithm.

By default, the oversampling size is set to 10.

This optional argument has no effect if logical argument RANDOM_QR is set to false or is absent in the call to TS_ID_CMP subroutine.

Further Details

For further details on the tsID and computing the tsID decomposition from the (randomized) partial QR factorizations with column pivoting of a matrix, see:

1. Martinsson, P.G., 2019:

   Randomized methods for matrix computations. arXiv.1607.01649

2. Erichson, N.B., Voronin, S., Brunton, S.L., and Kutz, J.N., 2019:

   Randomized matrix decompositions using R. arXiv.1608.02148

3. Voronin, S., Martinsson, P.G., 2015:

   Rsvdpack: Subroutines for computing partial singular value decompositions via randomized sampling on single core, multi core, and gpu architectures. arXiv.1502.05366

4. Voronin, S., Martinsson, P.G., 2017:

   Efficient algorithms for cur and interpolative matrix decompositions Adv Comput Math, volume 43, 495-516.

5. Stewart, G.W., 1999:

   Four algorithms for the the efficient computation of truncated pivoted qr approximations to a sparse matrix. Numerische Mathematik, Volume 83, 313-323.

6. Berry, M.W., Pulatova, S.A., and Stewart, G.W., 2005:

   Algorithm 844: Computing sparse reduced-rank approximations to sparse matrices. ACM Transactions on Mathematical Software, Volume 31, No. 2, 252-269.

subroutine cur_cmp ( mat, ip_row, ip_col, u, c, r, rnorm_row, rnorm_col, tol, random_qr, rng_alg, blk_size, nover )

Purpose

CUR_CMP computes a (partial) CUR decomposition of a m-by-n real matrix MAT. A CUR decomposition provides a reduced-rank approximate decomposition of a full m-by-n real matrix MAT of the form:

MAT = C * U * R

where C and R are m-by-krank and krank-by-n matrices, which are, respectively, subsets of the columns and rows of MAT, and U is a krank-by-krank matrix, which is estimated to make the matrix product C * U * R a good approximation of MAT according to the Frobenius norm. The CUR factorization is an important tool for handling large-scale data sets, offering several advantages over the Singular Value Decomposition (SVD): the columns and rows that comprise C and R are representative of the data and they are sparse if MAT is sparse. See references (1) and (2) for a discussion.

Computing an approximate CUR decomposition is generally a three steps process.

The C and R submatrices in the CUR factorization can be first estimated with the help of (randomized) partial QR factorizations with column pivoting of MAT and MAT’, respectively.

The initial step is thus to compute a (partial) QR factorization with column pivoting of MAT, which is defined as

MAT * P = Q * T

, where P is an n-by-n permutation matrix, T is an upper triangular or trapezoidal (i.e., if n>m) matrix and Q is a m-by-m orthogonal matrix if the QR factorization is complete.

For computing a CUR decomposition of rank, krank, a partial QR factorization with column pivoting of MAT of rank, krank, is sufficient, e.g., the QR decomposition can be stopped when the numbers of columns of Q is equal to krank .

This leads implicitly to the following partition of Q:

Q = [ Q1 Q2 ]

where Q1 is an m-by-krank orthonormal matrix and Q2 is m-by-(m-krank) orthonormal matrix orthogonal to Q1, and to the following corresponding partition of T:

[ T11 T12 ]

[ T21 T22 ]

where T11 is a krank-by-krank triangular matrix, T21 is zero by construction, T12 is a full krank-by-(n-krank) matrix and T22 is a full (m-krank)-by-(n-krank) matrix. This leads to the following (partial QR) approximation of MAT:

MAT = Q1 * [ T11 T12 ] * P’

Here, we see that Q1 * T11 equals the first krank columns of A * P, and we can define the matrix C in the CUR factorization of MAT as

C = Q1 * T11 = MAT * P(:,:krank)

This choice leads to a subset of the columns of MAT, which will minimize the Frobenius norm, RNORM_COL, of the error associated with the reduced-rank QR approximation or column ID of MAT and, also, of the resulting CUR approximation of MAT as we will illustrate below.

Furthermore, remember that C is also the selected subset of the columns of MAT, which defines the left factor in the column ID of MAT of rank krank:

MAT = (Q1 * T11) * V = C * V

where V = [ I Z ] * P’ is a krank-by-n matrix, I is a krank-by-krank identity matrix and Z is a krank-by-(n-krank) matrix, which is the solution to the matrix equation:

T11 * Z = T12 ,e.g., Z = inv(T11) * T12

see description of ID_CMP subroutine for details.

In a second step, the same partial QR algorithm is applied to MAT’ to get a representative subset of the rows of MAT, R, such that

R = L11 * K1 = N(:krank,:) * MAT

where R is a krank-by-n matrix, which is a subset of the rows of MAT, L11 is a krank-by-krank lower triangular matrix, K1 is a krank-by-n matrix with orthonormal rows and N is a m-by-m permutation matrix. Again, this choice leads to a subset of the rows of MAT, R, which will minimize the Frobenius norm, RNORM_ROW, of the error associated with the reduced-rank LQ approximation of MAT (and also of the resulting CUR approximation of MAT).

In a final step, we then seek a krank-by-krank matrix U such that

|| MAT - C * U * R ||_F = min

Once C and U are fixed, it can be shown that the minimum is attained for a matrix U, which is defined as

U = pseudo-inv(C) * MAT * pseudo-inv(R) = inv(T11) * Q1’ * MAT * K1’ * inv( L11 )

where pseudo-inv(C) is the pseudo-inverse of C, which is equal to inv(T11) * Q1’ if T11 is of full rank and with similar results for the matrix pseudo-inv(R). See reference (5) for details.

However, such direct computation of U involves the explicit inversion of two triangular matrices, which can be severely affected by ill-conditioning of the matrices T11 and L11. Using the column ID of MAT, we have

Q1’ * MAT = T11 * V

, we deduce that

U = inv(T11) * Q1’ * MAT * K1’ * inv( L11 ) = V * K1’ * inv( L11 ) = V * pseudo-inv(R)

and we observe that we can also determine U as the solution of the least squares problem:

U * R = V

Such a least squares problem must have an accurate solution as both the rows of R and V should span roughly the same space, namely, the space spanned by the krank leading right singular vectors of MAT (see reference (4) for discussion). If it is the case, the Frobenius norm of the error associated with the resulting CUR decomposition will be about of the same order of the error associated with the column ID of MAT, only slightly larger.

Using this approach, the matrix U can thus be estimated by solving the linear least squares problem:

U * R = V

and this computation never involves the explicit inversion of the matrices T11 and L11, but rather solving two linear least squares problems associated with these two triangular matrices.

Moreover, we have:

|| MAT - C * U * R ||_F <= sqrt( RNORM_COL**2 + RNORM_ROW**2 )

See references (5) and (6) for further details.

Note that for matrices whose singular values experience a fast decay, the accuracy of CUR factorization, as computed by CUR_CMP, can deteriorate due to ill-conditioning and the need to solve linear (least squares) problems associated with the triangular matrices T11 and L11 in order to estimate the matrix U. See reference (4) for a more detailed discussion of this problem.

The condition numbers (in the 1-norm) and the ranks of the matrices T11 and L11 are estimated if the optional argument TOL is present in order to check the robustness of the computed CUR decomposition.

krank is the target rank of the partial CUR decomposition, which is sought, and is equal to the number of rows or columns of the output array argument U, which stores the factor U in the CUR of MAT.

If the optional logical argument RANDOM_QR is used with the value true, fast randomized (partial) QR factorizations with column pivoting are used for computing the partial column QR and row LQ decompositions of MAT in the first two phases of the CUR algorithm.

Arguments

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-n matrix MAT.

IP_ROW (OUTPUT) integer(i4b), dimension(:)

On exit, IP_ROW stores the permutation matrix N in the partial QR decomposition with column pivoting of MAT’, which is computed to obtain the row LQ decomposition of MAT.

On exit, if IP_ROW(j)=k, then the j-th row of N*MAT was the k-th row of MAT, where N = I(IP_ROW,:) and I is the identity matrix of order m.

The matrix R in the CUR of MAT corresponds to the subset of the rows of MAT with the indices IP_ROW(:krank).

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP subroutine in module Random for further details.

The size of IP_ROW must be equal to size( MAT, 1 ) = m.

IP_COL (OUTPUT) integer(i4b), dimension(:)

On exit, IP_COL stores the permutation matrix P in the partial QR decomposition with column pivoting of MAT, which is computed to obtain the column QR decomposition of MAT.

If IP_COL(j)=k, then the j-th column of MAT*P was the k-th column of MAT, where P = I(:,IP_COL) and I is the identity matrix of order n.

The matrix C in the CUR of MAT corresponds to the subset of the columns of MAT with the indices IP_COL(:krank).

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP subroutine in module Random for further details.

The size of IP_COL must be equal to size( MAT, 2 ) = n.

U (OUTPUT) real(stnd), dimension(:,:)

On exit, the krank-by-krank matrix U in the approximate CUR factorization of MAT.

The shape of U must verify:

• size( U, 1 ) = size( U, 2 ) = krank <= min( size(MAT,1) , size(MAT,2) ) = min(m,n)

C (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, the m-by-krank matrix C in the approximate CUR factorization of MAT.

C consists of the subset of the columns of MAT defined by the indices IP_COL(:krank).

The shape of C must verify:

• size( C, 1 ) = size( MAT, 1 ) = m
• size( C, 2 ) = size( U, 2 ) = krank

R (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, the krank-by-n matrix R in the approximate CUR factorization of MAT.

R consists of the subset of the rows of MAT defined by the indices IP_ROW(:krank).

The shape of R must verify:

• size( R, 1 ) = size( U, 1 ) = krank
• size( R, 2 ) = size( MAT, 2 ) = n

RNORM_ROW (OUTPUT, OPTIONAL) real(stnd)

On exit, the Frobenius norm of the error matrix associated with the row partial LQ decomposition of MAT (e.g., the QR decomposition of MAT’), computed as || L22 ||_F .

RNORM_COL (OUTPUT, OPTIONAL) real(stnd)

On exit, the Frobenius norm of the error matrix associated with the column partial QR decomposition of MAT, computed as || T22 ||_F .

TOL (INPUT, OPTIONAL) real(stnd)

If TOL is present and is in [0,1[, then calculations to determine the condition number of submatrices T11 and L11 in the partial QR factorizations of MAT and MAT’ are performed.

TOL is used to determine the effective pseudo-ranks of T11 and L11, which are defined as the order of the largest leading triangular submatrices in the partial QR factorizations with column pivoting of MAT and MAT’ whose estimated condition numbers in the 1-norm are less than 1/TOL. If TOL=0 is specified, the calculations to determine the condition numbers of T11 and L11 are not performed and crude tests on T(j,j) and L(j,j) are done to determine the numerical pseudo-ranks of T11 and L11.

The TOL argument is useful to check if the matrices T11 and L11 are sufficiently well conditioned to obtain a stable and robust CUR decomposition.

If the computed pseudo-ranks of T11 or L11 are less than krank = size( U, 1 ), the subroutine will exit with an error message.

On the other hand, if TOL is not specified or is outside [0,1[, the calculations to determine the ranks of T11 and L11 are not performed and these ranks are assumed to be equal to krank = size( U, 1 ).

If it is possible that MAT, T11 and L11 may not be of full rank, then the stability of the CUR decomposition of MAT, which is sought, can be checked by using TOL=relative precision of the elements in MAT. If each element of MAT is correct to, say, 5 digits then TOL=0.00001 should be used.

See description of PARTIAL_QR_CMP subroutine in module QR_Procedures and PARTIAL_RQR_CMP subroutine in module Random for further details.

RANDOM_QR (INPUT, OPTIONAL) logical(lgl)

On entry, if RANDOM_QR is used with the value true, fast randomized (partial) QR factorizations are used in the first two phases of the CUR algorithm.

By default, RANDOM_QR = false, i.e., a standard (partial) QR factorization with column pivoting is used in the first two phases of the CUR algorithm.

RNG_ALG (INPUT, OPTIONAL) integer(i4b)

On entry, a scalar integer to select the random (uniform) number generator used to build the random gaussian matrix in the two randomized (partial) QR phases of the CUR algorithm if RANDOM_QR = true.

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 current random number generator and its current state are not changed. Note further, that, on exit, the current random number generator is not reset to its previous value before the call to CUR_CMP subroutine.

See the documentation of subroutine RANDOM_SEED_ in module Random for further information.

This optional argument has no effect if logical argument RANDOM_QR is set to false or is absent in the call to CUR_CMP subroutine.

BLK_SIZE (INPUT, OPTIONAL) integer(i4b)

On entry, the block size used in the randomized partial QR phases of the CUR algorithm if RANDOM_QR = true.

BLK_SIZE must be greater or equal to one and less than min(m,n) and must be set to a much smaller value than min(m,n) usually, depending also on the architecture of the computer.

By default, BLK_SIZE is set to min( BLKSZ_QR, min(m,n) ), where parameter BLKSZ_QR is the default block size for QR related algorithms specified in module Select_Parameters.

See description of subroutine PARTIAL_RQR_CMP in module Random for the meaning of the the block size in the randomized (partial) QR algorithm.

This optional argument has no effect if logical argument RANDOM_QR is set to false or is absent in the call to CUR_CMP subroutine.

NOVER (INPUT, OPTIONAL) integer(i4b)

The oversampling size used in the randomized partial QR phases of the CUR algorithm if RANDOM_QR = true.

NOVER must be positive or null and verifies the relationship:

NOVER + BLK_SIZE <= size(MAT,1)

and is adjusted if necessary to verify this relationship in all cases.

See description of subroutine PARTIAL_RQR_CMP in module Random for the meaning and usefulness of the oversampling size in the randomized (partial) QR algorithm.

By default, the oversampling size is set to 10.

This optional argument has no effect if logical argument RANDOM_QR is set to false or is absent in the call to CUR_CMP subroutine.

Further Details

For further details on the CUR and computing the CUR decomposition from (randomized) partial QR factorizations with column pivoting of a matrix and its transpose, see:

1. Mahoney, M.W., and Drineas, P., 2009:

   CUR matrix decompositions for improved data analysis. PNAS, Volume 106, No. 3, 697-702.

2. Martinsson, P.G., 2019:

   Randomized methods for matrix computations. arXiv.1607.01649

3. Voronin, S., Martinsson, P.G., 2015:

   Rsvdpack: Subroutines for computing partial singular value decompositions via randomized sampling on single core, multi core, and gpu architectures. arXiv.1502.05366

4. Voronin, S., Martinsson, P.G., 2017:

   Efficient algorithms for cur and interpolative matrix decompositions Adv Comput Math, Volume 43, 495-516.

5. Stewart, G.W., 1999:

   Four algorithms for the the efficient computation of truncated pivoted qr approximations to a sparse matrix. Numerische Mathematik, Volume 83, 313-323.

6. Berry, M.W., Pulatova, S.A., and Stewart, G.W., 2005:

   Algorithm 844: Computing sparse reduced-rank approximations to sparse matrices. ACM Transactions on Mathematical Software, Volume 31, No. 2, 252-269.

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:

1. 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:

1. 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:

1. 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:

1. 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:

1. Noreen, E.W., 1989:

   Computer-intensive methods for testing hypotheses: an introduction. Wiley and Sons, New York, USA, ISBN:978-0-471-61136-3