# include # include # include # include # include # include "test_eigen.h" /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the smaller of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, two integers to be compared. Output, int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_uniform_ab ( int a, int b, int *seed ) /******************************************************************************/ /* Purpose: I4_UNIFORM_AB returns a scaled pseudorandom I4 between A and B. Discussion: The pseudorandom number should be uniformly distributed between A and B. Licensing: This code is distributed under the MIT license. Modified: 24 May 2012 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Parameters: Input, int A, B, the limits of the interval. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, int I4_UNIFORM_AB, a number between A and B. */ { int c; const int i4_huge = 2147483647; int k; float r; int value; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_UNIFORM_AB - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } /* Guaranteee A <= B. */ if ( b < a ) { c = a; a = b; b = c; } k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( float ) ( *seed ) * 4.656612875E-10; /* Scale R to lie between A-0.5 and B+0.5. */ r = ( 1.0 - r ) * ( ( float ) ( a ) - 0.5 ) + r * ( ( float ) ( b ) + 0.5 ); /* Round R to the nearest integer. */ value = round ( r ); /* Guarantee that A <= VALUE <= B. */ if ( value < a ) { value = a; } if ( b < value ) { value = b; } return value; } /******************************************************************************/ double r8_normal_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_NORMAL_01 samples the standard normal probability distribution. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. The Box-Muller method is used, which is efficient, but generates two values at a time. Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input/output, int *SEED, a seed for the random number generator. Output, double R8_NORMAL_01, a normally distributed random value. */ { double pi = 3.141592653589793; double r1; double r2; static int used = -1; double x; static double y = 0.0; if ( used == -1 ) { used = 0; } /* If we've used an even number of values so far, generate two more, return one, and save one. */ if ( ( used % 2 )== 0 ) { for ( ; ; ) { r1 = r8_uniform_01 ( seed ); if ( r1 != 0.0 ) { break; } } r2 = r8_uniform_01 ( seed ); x = sqrt ( -2.0 * log ( r1 ) ) * cos ( 2.0 * pi * r2 ); y = sqrt ( -2.0 * log ( r1 ) ) * sin ( 2.0 * pi * r2 ); } else { x = y; } used = used + 1; return x; } /******************************************************************************/ double r8_sign ( double x ) /******************************************************************************/ /* Purpose: R8_SIGN returns the sign of an R8. Licensing: This code is distributed under the MIT license. Modified: 08 May 2006 Author: John Burkardt Parameters: Input, double X, the number whose sign is desired. Output, double R8_SIGN, the sign of X. */ { double value; if ( x < 0.0 ) { value = - 1.0; } else { value = + 1.0; } return value; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01 returns a unit pseudorandom R8. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r8_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R8_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the MIT license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. P A Lewis, A S Goodman, J M Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } r = ( ( double ) ( *seed ) ) * 4.656612875E-10; return r; } /******************************************************************************/ void r8bin_print ( int bin_num, int bin[], double bin_limit[], char *title ) /******************************************************************************/ /* Purpose: R8BIN_PRINT prints the bins of a real vector. Licensing: This code is distributed under the MIT license. Modified: 22 February 2012 Author: John Burkardt Parameters: Input, int BIN_NUM, the number of bins. Input, int BIN[BIN_NUM+2]. BIN(0) counts entries of X less than BIN_LIMIT(0). BIN(BIN_NUM+1) counts entries greater than or equal to BIN_LIMIT(BIN_NUM). For 1 <= I <= BIN_NUM, BIN(I) counts the entries X(J) such that BIN_LIMIT(I-1) <= X(J) < BIN_LIMIT(I). where H is the bin spacing. Input, double BIN_LIMIT[BIN_NUM+1], the "limits" of the bins. BIN(I) counts the number of entries X(J) such that BIN_LIMIT(I-1) <= X(J) < BIN_LIMIT(I). Input, char *TITLE, a title. */ { int i; printf ( "\n" ); printf ( "%s\n", title ); printf ( "\n" ); printf ( " Index Lower Limit Count Upper Limit\n" ); printf ( "\n" ); printf ( " %6d %6d %14g\n", 0, bin[0], bin_limit[0] ); for ( i = 1; i <= bin_num; i++ ) { printf ( " %6d %14g %6d %14g\n", i, bin_limit[i-1], bin[i], bin_limit[i] ); } printf ( " %6d %14g %6d\n", bin_num + 1, bin_limit[bin_num], bin[bin_num+1] ); return; } /******************************************************************************/ void r8mat_copy ( int m, int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8MAT_COPY copies one R8MAT to another. Discussion: An R8MAT is a doubly dimensioned array of R8's, which may be stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 26 July 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A1[M*N], the matrix to be copied. Output, double A2[M*N], the copy of A1. */ { int i; int j; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { a2[i+j*m] = a1[i+j*m]; } } return; } /******************************************************************************/ double *r8mat_house_axh_new ( int n, double a[], double v[] ) /******************************************************************************/ /* Purpose: R8MAT_HOUSE_AXH_NEW computes A*H where H is a compact Householder matrix. Discussion: An R8MAT is a doubly dimensioned array of double precision values, which may be stored as a vector in column-major order. The Householder matrix H(V) is defined by H(V) = I - 2 * v * v' / ( v' * v ) This routine is not particularly efficient. Licensing: This code is distributed under the MIT license. Modified: 05 October 2010 Author: John Burkardt Parameters: Input, int N, the order of A. Input, double A[N*N], the matrix to be postmultiplied. Input, double V[N], a vector defining a Householder matrix. Output, double R8MAT_HOUSE_AXH[N*N], the product A*H. */ { double *ah; int i; int j; int k; double v_normsq; v_normsq = 0.0; for ( i = 0; i < n; i++ ) { v_normsq = v_normsq + v[i] * v[i]; } /* Compute A*H' = A*H */ ah = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { for ( j = 0; j < n; j++ ) { ah[i+j*n] = a[i+j*n]; for ( k = 0; k < n; k++ ) { ah[i+j*n] = ah[i+j*n] - 2.0 * a[i+k*n] * v[k] * v[j] / v_normsq; } } } return ah; } /******************************************************************************/ void r8mat_identity ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8MAT_IDENTITY sets an R8MAT to the identity matrix. Discussion: An R8MAT is a doubly dimensioned array of R8's, which may be stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 06 September 2005 Author: John Burkardt Parameters: Input, int N, the order of A. Output, double A[N*N], the N by N identity matrix. */ { int i; int j; int k; k = 0; for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { if ( i == j ) { a[k] = 1.0; } else { a[k] = 0.0; } k = k + 1; } } return; } /******************************************************************************/ double *r8mat_mm_new ( int n1, int n2, int n3, double a[], double b[] ) /******************************************************************************/ /* Purpose: R8MAT_MM_NEW multiplies two matrices. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. For this routine, the result is returned as the function value. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Parameters: Input, int N1, N2, N3, the order of the matrices. Input, double A[N1*N2], double B[N2*N3], the matrices to multiply. Output, double R8MAT_MM[N1*N3], the product matrix C = A * B. */ { double *c; int i; int j; int k; c = ( double * ) malloc ( n1 * n3 * sizeof ( double ) ); for ( i = 0; i < n1; i ++ ) { for ( j = 0; j < n3; j++ ) { c[i+j*n1] = 0.0; for ( k = 0; k < n2; k++ ) { c[i+j*n1] = c[i+j*n1] + a[i+k*n1] * b[k+j*n2]; } } } return c; } /******************************************************************************/ void r8mat_orth_uniform ( int n, int *seed, double q[] ) /******************************************************************************/ /* Purpose: R8MAT_ORTH_UNIFORM returns a random orthogonal matrix. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. The inverse of A is equal to A'. A * A' = A' * A = I. Columns and rows of A have unit Euclidean norm. Distinct pairs of columns of A are orthogonal. Distinct pairs of rows of A are orthogonal. The L2 vector norm of A*x = the L2 vector norm of x for any vector x. The L2 matrix norm of A*B = the L2 matrix norm of B for any matrix B. The determinant of A is +1 or -1. All the eigenvalues of A have modulus 1. All singular values of A are 1. All entries of A are between -1 and 1. Discussion: Thanks to Eugene Petrov, B I Stepanov Institute of Physics, National Academy of Sciences of Belarus, for convincingly pointing out the severe deficiencies of an earlier version of this routine. Essentially, the computation involves saving the Q factor of the QR factorization of a matrix whose entries are normally distributed. However, it is only necessary to generate this matrix a column at a time, since it can be shown that when it comes time to annihilate the subdiagonal elements of column K, these (transformed) elements of column K are still normally distributed random values. Hence, there is no need to generate them at the beginning of the process and transform them K-1 times. For computational efficiency, the individual Householder transformations could be saved, as recommended in the reference, instead of being accumulated into an explicit matrix format. Licensing: This code is distributed under the MIT license. Modified: 20 January 2020 Author: John Burkardt Reference: Pete Stewart, Efficient Generation of Random Orthogonal Matrices With an Application to Condition Estimators, SIAM Journal on Numerical Analysis, Volume 17, Number 3, June 1980, pages 403-409. Parameters: Input, int N, the order of A. Input/output, int *SEED, a seed for the random number generator. Output, double Q[N*N], the orthogonal matrix. */ { double *a_col; int i; int j; int k; double *q2; double *v; /* Start with Q = the identity matrix. */ r8mat_identity ( n, q ); /* Now behave as though we were computing the QR factorization of some other random matrix A. Generate the N elements of the first column, compute the Householder matrix H1 that annihilates the subdiagonal elements, and set Q := Q * H1' = Q * H. On the second step, generate the lower N-1 elements of the second column, compute the Householder matrix H2 that annihilates them, and set Q := Q * H2' = Q * H2 = H1 * H2. On the N-1 step, generate the lower 2 elements of column N-1, compute the Householder matrix HN-1 that annihilates them, and and set Q := Q * H(N-1)' = Q * H(N-1) = H1 * H2 * ... * H(N-1). This is our random orthogonal matrix. */ a_col = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 1; j < n; j++ ) { /* Set the vector that represents the J-th column to be annihilated. */ for ( i = 1; i < j; i++ ) { a_col[i-1] = 0.0; } for ( i = j; i <= n; i++ ) { a_col[i-1] = r8_normal_01 ( seed ); } /* Compute the vector V that defines a Householder transformation matrix H(V) that annihilates the subdiagonal elements of A. */ v = r8vec_house_column ( n, a_col, j ); /* Postmultiply the matrix Q by H'(V) = H(V). */ q2 = r8mat_house_axh_new ( n, q, v ); free ( v ); /* Randomly insert reflections. */ i = i4_uniform_ab ( 0, n - 1, seed ); for ( k = 0; k < n; k++ ) { q2[i+k*n] = - q2[i+k*n]; } r8mat_copy ( n, n, q2, q ); free ( q2 ); } free ( a_col ); return; } /******************************************************************************/ void r8mat_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT prints an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8's, which may be stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*M] Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, double A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { r8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT_SOME prints some of an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8's, which may be stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, double A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2hi; int i2lo; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); if ( m <= 0 || n <= 0 ) { fprintf ( stdout, "\n" ); fprintf ( stdout, " (None)\n" ); return; } /* Print the columns of the matrix, in strips of 5. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; j2hi = i4_min ( j2hi, n ); j2hi = i4_min ( j2hi, jhi ); fprintf ( stdout, "\n" ); /* For each column J in the current range... Write the header. */ fprintf ( stdout, " Col: "); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %7d ", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ i2lo = i4_max ( ilo, 1 ); i2hi = i4_min ( ihi, m ); for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) 5 entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %14f", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ void r8nsymm_gen ( int n, double lambda_mean, double lambda_dev, int *seed, double a[], double ql[], double qr[], double lambda[] ) /******************************************************************************/ /* Purpose: R8NSYMM_GEN generates a nonsymmetric matrix with a certain eigenstructure. Discussion: An R8NSYMM is a real nonsymmetric matrix stored using full storage, and using R8 arithmetic. Licensing: This code is distributed under the MIT license. Modified: 09 March 2018 Author: John Burkardt Parameters: Input, int N, the order of the matrix. Input, double LAMBDA_MEAN, the mean value of the normal distribution from which the eigenvalues will be chosen. Input, double LAMBDA_DEV, the standard deviation of the normal distribution from which the eigenvalues will be chosen. Input/output, int *SEED, a seed for the random number generator. Output, double A[N*N], the test matrix. Output, double QL[N*N], the left eigenvector matrix. Output, double QR[N*N], the right eigenvector matrix. Output, double LAMBDA[N], the eigenvalue vector. */ { int i; int j; int k; /* Choose the eigenvalues LAMBDA. */ r8vec_normal_ab ( n, lambda_mean, lambda_dev, seed, lambda ); /* Get random orthogonal matrices QL and QR. */ r8mat_orth_uniform ( n, seed, ql ); r8mat_orth_uniform ( n, seed, qr ); /* Set A = QL * Lambda*I * QR'. */ for ( i = 0; i < n; i++ ) { for ( j = 0; j < n; j++ ) { a[i+j*n] = 0.0; for ( k = 0; k < n; k++ ) { a[i+j*n] = a[i+j*n] + ql[i+k*n] * lambda[k] * qr[j+k*n]; } } } return; } /******************************************************************************/ void r8symm_gen ( int n, double lambda_mean, double lambda_dev, int *seed, double a[], double q[], double lambda[] ) /******************************************************************************/ /* Purpose: R8SYMM_GEN generates a symmetric matrix with a certain eigenstructure. Discussion: An R8SYMM is a real symmetric matrix stored using full storage, and using R8 arithmetic. Licensing: This code is distributed under the MIT license. Modified: 22 February 2012 Author: John Burkardt Parameters: Input, int N, the order of the matrix. Input, double LAMBDA_MEAN, the mean value of the normal distribution from which the eigenvalues will be chosen. Input, double LAMBDA_DEV, the standard deviation of the normal distribution from which the eigenvalues will be chosen. Input/output, int *SEED, a seed for the random number generator. Output, double A[N*N], the test matrix. Output, double Q[N*N], the eigenvector matrix. Output, double LAMBDA[N], the eigenvalue vector. */ { int i; int j; int k; /* Choose the eigenvalues LAMBDA. */ r8vec_normal_ab ( n, lambda_mean, lambda_dev, seed, lambda ); /* Get a random orthogonal matrix Q. */ r8mat_orth_uniform ( n, seed, q ); /* Set A = Q * Lambda*I * Q'. */ for ( i = 0; i < n; i++ ) { for ( j = 0; j < n; j++ ) { a[i+j*n] = 0.0; for ( k = 0; k < n; k++ ) { a[i+j*n] = a[i+j*n] + q[i+k*n] * lambda[k] * q[j+k*n]; } } } return; } /******************************************************************************/ void r8vec_bin ( int n, double x[], int bin_num, double bin_min, double bin_max, int bin[], double bin_limit[] ) /******************************************************************************/ /* Purpose: R8VEC_BIN computes bins based on a given R8VEC. Discussion: The user specifies minimum and maximum bin values, BIN_MIN and BIN_MAX, and the number of bins, BIN_NUM. This determines a "bin width": H = ( BIN_MAX - BIN_MIN ) / BIN_NUM so that bin I will count all entries X(J) such that BIN_LIMIT(I-1) <= X(J) < BIN_LIMIT(I). The array X does NOT have to be sorted. Licensing: This code is distributed under the MIT license. Modified: 22 February 2012 Author: John Burkardt Parameters: Input, int N, the number of entries of X. Input, double X[N], an (unsorted) array to be binned. Input, int BIN_NUM, the number of bins. Two extra bins, #0 and #BIN_NUM+1, count extreme values. Input, double BIN_MIN, BIN_MAX, define the range and size of the bins. BIN_MIN and BIN_MAX must be distinct. Normally, BIN_MIN < BIN_MAX, and the documentation will assume this, but proper results will be computed if BIN_MIN > BIN_MAX. Output, int BIN[BIN_NUM+2]. BIN(0) counts entries of X less than BIN_MIN. BIN(BIN_NUM+1) counts entries greater than or equal to BIN_MAX. For 1 <= I <= BIN_NUM, BIN(I) counts the entries X(J) such that BIN_LIMIT(I-1) <= X(J) < BIN_LIMIT(I). where H is the bin spacing. Output, double BIN_LIMIT[BIN_NUM+1], the "limits" of the bins. BIN(I) counts the number of entries X(J) such that BIN_LIMIT(I-1) <= X(J) < BIN_LIMIT(I). */ { int i; int j; double t; if ( bin_max == bin_min ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_BIN - Fatal error!\n" ); fprintf ( stderr, " BIN_MIN = BIN_MAX = %g\n", bin_max ); exit ( 1 ); } for ( i = 0; i <= bin_num + 1; i++ ) { bin[i] = 0; } for ( i = 0; i < n; i++ ) { t = ( x[i] - bin_min ) / ( bin_max - bin_min ); if ( t < 0.0 ) { j = 0; } else if ( 1.0 <= t ) { j = bin_num + 1; } else { j = 1 + ( int ) ( ( double ) ( bin_num ) * t ); } bin[j] = bin[j] + 1; } /* Compute the bin limits. */ for ( i = 0; i <= bin_num; i++ ) { bin_limit[i] = ( ( double ) ( bin_num - i ) * bin_min + ( double ) ( i ) * bin_max ) / ( double ) ( bin_num ); } return; } /******************************************************************************/ void r8vec_copy ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8VEC_COPY copies an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 03 July 2005 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], the vector to be copied. Input, double A2[N], the copy of A1. */ { int i; for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return; } /******************************************************************************/ double *r8vec_house_column ( int n, double a[], int k ) /******************************************************************************/ /* Purpose: R8VEC_HOUSE_COLUMN defines a Householder premultiplier that "packs" a column. Discussion: An R8VEC is a vector of R8's. The routine returns a vector V that defines a Householder premultiplier matrix H(V) that zeros out the subdiagonal entries of column K of the matrix A. H(V) = I - 2 * v * v' Licensing: This code is distributed under the MIT license. Modified: 21 August 2010 Author: John Burkardt Parameters: Input, int N, the order of the matrix A. Input, double A[N], column K of the matrix A. Input, int K, the column of the matrix to be modified. Output, double R8VEC_HOUSE_COLUMN[N], a vector of unit L2 norm which defines an orthogonal Householder premultiplier matrix H with the property that the K-th column of H*A is zero below the diagonal. */ { int i; double s; double *v; v = r8vec_zero_new ( n ); if ( k < 1 || n <= k ) { return v; } s = r8vec_norm_l2 ( n+1-k, a+k-1 ); if ( s == 0.0 ) { return v; } v[k-1] = a[k-1] + fabs ( s ) * r8_sign ( a[k-1] ); r8vec_copy ( n-k, a+k, v+k ); s = r8vec_norm_l2 ( n-k+1, v+k-1 ); for ( i = k-1; i < n; i++ ) { v[i] = v[i] / s; } return v; } /******************************************************************************/ double r8vec_max ( int n, double r8vec[] ) /******************************************************************************/ /* Purpose: R8VEC_MAX returns the value of the maximum element in a R8VEC. Licensing: This code is distributed under the MIT license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double R8VEC[N], a pointer to the first entry of the array. Output, double R8VEC_MAX, the value of the maximum element. This is set to 0.0 if N <= 0. */ { int i; double value; if ( n <= 0 ) { value = 0.0; return value; } value = r8vec[0]; for ( i = 1; i < n; i++ ) { if ( value < r8vec[i] ) { value = r8vec[i]; } } return value; } /******************************************************************************/ double r8vec_min ( int n, double r8vec[] ) /******************************************************************************/ /* Purpose: R8VEC_MIN returns the value of the minimum element in a R8VEC. Licensing: This code is distributed under the MIT license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double R8VEC[N], the array to be checked. Output, double R8VEC_MIN, the value of the minimum element. */ { int i; double value; value = r8vec[0]; for ( i = 1; i < n; i++ ) { if ( r8vec[i] < value ) { value = r8vec[i]; } } return value; } /******************************************************************************/ double r8vec_norm_l2 ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_NORM_L2 returns the L2 norm of an R8VEC. Discussion: The vector L2 norm is defined as: R8VEC_NORM_L2 = sqrt ( sum ( 1 <= I <= N ) A(I)^2 ). Licensing: This code is distributed under the MIT license. Modified: 01 March 2003 Author: John Burkardt Parameters: Input, int N, the number of entries in A. Input, double A[N], the vector whose L2 norm is desired. Output, double R8VEC_NORM_L2, the L2 norm of A. */ { int i; double v; v = 0.0; for ( i = 0; i < n; i++ ) { v = v + a[i] * a[i]; } v = sqrt ( v ); return v; } /******************************************************************************/ void r8vec_normal_ab ( int n, double b, double c, int *seed, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_NORMAL_AB returns a scaled pseudonormal R8VEC. Discussion: The scaled normal probability distribution function (PDF) has mean A and standard deviation B. This routine can generate a vector of values on one call. It has the feature that it should provide the same results in the same order no matter how we break up the task. Before calling this routine, the user may call RANDOM_SEED in order to set the seed of the random number generator. The Box-Muller method is used, which is efficient, but generates an even number of values each time. On any call to this routine, an even number of new values are generated. Depending on the situation, one value may be left over. In that case, it is saved for the next call. Licensing: This code is distributed under the MIT license. Modified: 02 February 2005 Author: John Burkardt Parameters: Input, int N, the number of values desired. If N is negative, then the code will flush its internal memory; in particular, if there is a saved value to be used on the next call, it is instead discarded. This is useful if the user has reset the random number seed, for instance. Input, double B, C, the mean and standard deviation. Input/output, int *SEED, a seed for the random number generator. Output, double X[N], a sample of the standard normal PDF. Local parameters: Local, int MADE, records the number of values that have been computed. On input with negative N, this value overwrites the return value of N, so the user can get an accounting of how much work has been done. Local, double R[N+1], is used to store some uniform random values. Its dimension is N+1, but really it is only needed to be the smallest even number greater than or equal to N. Local, int SAVED, is 0 or 1 depending on whether there is a single saved value left over from the previous call. Local, int X_LO, X_HI, records the range of entries of X that we need to compute. This starts off as 1:N, but is adjusted if we have a saved value that can be immediately stored in X(1), and so on. Local, double Y, the value saved from the previous call, if SAVED is 1. */ { # define R8_PI 3.141592653589793 int i; int m; static int made = 0; double *r; static int saved = 0; int x_hi; int x_lo; static double y = 0.0; /* I'd like to allow the user to reset the internal data. But this won't work properly if we have a saved value Y. I'm making a crock option that allows the user to signal explicitly that any internal memory should be flushed, by passing in a negative value for N. */ if ( n < 0 ) { made = 0; saved = 0; y = 0.0; return; } else if ( n == 0 ) { return; } /* Record the range of X we need to fill in. */ x_lo = 1; x_hi = n; /* Use up the old value, if we have it. */ if ( saved == 1 ) { x[0] = y; saved = 0; x_lo = 2; } /* Maybe we don't need any more values. */ if ( x_hi - x_lo + 1 == 0 ) { } /* If we need just one new value, do that here to avoid null arrays. */ else if ( x_hi - x_lo + 1 == 1 ) { r = r8vec_uniform_01_new ( 2, seed ); x[x_hi-1] = sqrt ( - 2.0 * log ( r[0] ) ) * cos ( 2.0 * R8_PI * r[1] ); y = sqrt ( - 2.0 * log ( r[0] ) ) * sin ( 2.0 * R8_PI * r[1] ); saved = 1; made = made + 2; free ( r ); } /* If we require an even number of values, that's easy. */ else if ( ( x_hi - x_lo + 1 ) % 2 == 0 ) { m = ( x_hi - x_lo + 1 ) / 2; r = r8vec_uniform_01_new ( 2*m, seed ); for ( i = 0; i <= 2*m-2; i = i + 2 ) { x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); x[x_lo+i ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); } made = made + x_hi - x_lo + 1; free ( r ); } /* If we require an odd number of values, we generate an even number, and handle the last pair specially, storing one in X(N), and saving the other for later. */ else { x_hi = x_hi - 1; m = ( x_hi - x_lo + 1 ) / 2 + 1; r = r8vec_uniform_01_new ( 2*m, seed ); for ( i = 0; i <= 2*m-4; i = i + 2 ) { x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); x[x_lo+i ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); } i = 2 * m - 2; x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); y = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); saved = 1; made = made + x_hi - x_lo + 2; free ( r ); } for ( i = 0; i < n; i++ ) { x[i] = b + c * x[i]; } return; # undef R8_PI } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_PRINT prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14f\n", i, a[i] ); } return; } /******************************************************************************/ double *r8vec_uniform_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 19 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Parameters: Input, int N, the number of entries in the vector. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; int i4_huge = 2147483647; int k; double *r; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_UNIFORM_01_NEW - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i] = ( double ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ double *r8vec_zero_new ( int n ) /******************************************************************************/ /* Purpose: R8VEC_ZERO_NEW creates and zeroes an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 25 March 2009 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, double R8VEC_ZERO_NEW[N], a vector of zeroes. */ { double *a; int i; a = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { a[i] = 0.0; } return a; } /******************************************************************************/ void r8vec2_print ( int n, double a1[], double a2[], char *title ) /******************************************************************************/ /* Purpose: R8VEC2_PRINT prints an R8VEC2. Discussion: An R8VEC2 is a dataset consisting of N pairs of real values, stored as two separate vectors A1 and A2. Licensing: This code is distributed under the MIT license. Modified: 26 March 2009 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A1[N], double A2[N], the vectors to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %4d: %14f %14f\n", i, a1[i], a2[i] ); } return; } /******************************************************************************/ void timestamp ( void ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }