运行：
./a.out 7
编译：
g++ ./hello_svd_jacobi.cpp
#include <iostream>
#include <cmath>
#include <cassert>
#include <iomanip>
#include <cstdlib>
#include <chrono>
#include <random>
#include <algorithm>
#include <vector>
#include <array>
using std::cout;
using std::cerr;
using std::endl;
using std::setprecision;
using std::vector;
using std::array;
#include<cassert>
namespace matrix_alloc_jpd {
/// @brief Allocate a 2-dimensional array. (Uses row-major order.)
template<typename Entry>
void Alloc2D(size_t nrows, //!< size of the array (number of rows)
size_t ncols, //!< size of the array (number of columns)
Entry ***paaX //!< pointer to a 2D C-style array
);
/// @brief Deallocate arrays that were created using Alloc2D().
template<typename Entry>
void Dealloc2D(Entry ***paaX //!< pointer to a 2D C-style array
);
// ---- IMPLEMENTATION ----
template<typename Entry>
void Alloc2D(size_t nrows, // size of the array (number of rows)
size_t ncols, // size of the array (number of columns)
Entry ***paaX) // pointer to a 2D C-style array
{
assert(paaX);
*paaX = new Entry* [nrows]; //conventional 2D C array (pointer-to-pointer)
(*paaX)[0] = new Entry [nrows * ncols]; // 1D C array (contiguous memory)
for(size_t iy=0; iy<nrows; iy++)
(*paaX)[iy] = (*paaX)[0] + iy*ncols;
// The caller can access the contents using (*paaX)[i][j]
}
template<typename Entry>
void Dealloc2D(Entry ***paaX) // pointer to a 2D C-style array
{
if (paaX && *paaX) {
delete [] (*paaX)[0];
delete [] (*paaX);
*paaX = nullptr;
}
}
} // namespace matrix_alloc_jpd
//
#include <algorithm>
#include <cmath>
//#include <cassert>
namespace jacobi_pd {
using namespace matrix_alloc_jpd;
/// @class Jacobi
/// @brief Calculate the eigenvalues and eigevectors of a symmetric matrix
/// using the Jacobi eigenvalue algorithm.
/// @note The "Vector" and "Matrix" type arguments can be any
/// C or C++ object that support indexing, including pointers or vectors.
template<typename Scalar,
typename Vector,
typename Matrix,
typename ConstMatrix=Matrix>
class Jacobi
{
int n; //!< the size of the matrix
Scalar **M; //!< local copy of the matrix being analyzed
// Precomputed cosine, sine, and tangent of the most recent rotation angle:
Scalar c; //!< = cos(θ)
Scalar s; //!< = sin(θ)
Scalar t; //!< = tan(θ), (note |t|<=1)
int *max_idx_row; //!< for row i, the index j of the maximum element where j>i
public:
/// @brief Specify the size of the matrices you want to diagonalize later.
/// @param n the size (ie. number of rows) of the square matrix
void SetSize(int n);
Jacobi(int n = 0) {
Init();
SetSize(n);
}
~Jacobi() {
Dealloc();
}
// @typedef choose the criteria for sorting eigenvalues and eigenvectors
typedef enum eSortCriteria {
DO_NOT_SORT,
SORT_DECREASING_EVALS,
SORT_INCREASING_EVALS,
SORT_DECREASING_ABS_EVALS,
SORT_INCREASING_ABS_EVALS
} SortCriteria;
/// @brief Calculate all the eigenvalues and eigevectors of a symmetric matrix
/// using the Jacobi eigenvalue algorithm:
/// https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
/// @returns The number of Jacobi iterations attempted, which should be > 0.
/// If the return value is not strictly > 0 then convergence failed.
/// @note To reduce the computation time further, set calc_evecs=false.
int
Diagonalize(ConstMatrix mat, //!< the matrix you wish to diagonalize (size n)
Vector eval, //!< store the eigenvalues here
Matrix evec, //!< store the eigenvectors here (in rows)
SortCriteria sort_criteria=SORT_DECREASING_EVALS,//!<sort results?
bool calc_evecs=true, //!< calculate the eigenvectors?
int max_num_sweeps = 50); //!< limit the number of iterations
private:
/// @brief Calculate the components of a rotation matrix which performs a
/// rotation in the i,j plane by an angle (θ) causing M[i][j]=0.
/// The resulting parameters will be stored in this->c, this->s, and
/// this->t (which store cos(θ), sin(θ), and tan(θ), respectively).
void CalcRot(Scalar const *const *M, //!< matrix
int i, //!< row index
int j); //!< column index
/// @brief Apply the (previously calculated) rotation matrix to matrix M
/// by multiplying it on both sides (a "similarity transform").
/// (To save time, only the elements in the upper-right-triangular
/// region of the matrix are updated. It is assumed that i < j.)
void ApplyRot(Scalar **M, //!< matrix
int i, //!< row index
int j); //!< column index
/// @brief Multiply matrix E on the left by the (previously calculated)
/// rotation matrix.
void ApplyRotLeft(Matrix E, //!< matrix
int i, //!< row index
int j); //!< column index
///@brief Find the off-diagonal index in row i whose absolute value is largest
int MaxEntryRow(Scalar const *const *M, int i) const;
/// @brief Find the indices (i_max, j_max) marking the location of the
/// entry in the matrix with the largest absolute value. This
/// uses the max_idx_row[] array to find the answer in O(n) time.
/// @returns This function does not return a value. However after it is
/// invoked, the location of the largest matrix element will be
/// stored in the i_max and j_max arguments.
void MaxEntry(Scalar const *const *M, int& i_max, int& j_max) const;
// @brief Sort the rows in M according to the numbers in v (also sorted)
void SortRows(Vector v, //!< vector containing the keys used for sorting
Matrix M, //!< matrix whose rows will be sorted according to v
int n, //!< size of the vector and matrix
SortCriteria s=SORT_DECREASING_EVALS //!< sort decreasing order?
) const;
// memory management:
void Alloc(int N);
void Init();
void Dealloc();
public:
// memory management: copy and move constructor, swap, and assignment operator
Jacobi(const Jacobi<Scalar, Vector, Matrix, ConstMatrix>& source);
Jacobi(Jacobi<Scalar, Vector, Matrix, ConstMatrix>&& other);
void swap(Jacobi<Scalar, Vector, Matrix, ConstMatrix> &other) noexcept;
Jacobi<Scalar, Vector, Matrix, ConstMatrix>& operator = (Jacobi<Scalar, Vector, Matrix, ConstMatrix> source);
}; // class Jacobi
// -------------- IMPLEMENTATION --------------
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
int Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Diagonalize(ConstMatrix mat, // the matrix you wish to diagonalize (size n)
Vector eval, // store the eigenvalues here
Matrix evec, // store the eigenvectors here (in rows)
SortCriteria sort_criteria, // sort results?
bool calc_evec, // calculate the eigenvectors?
int max_num_sweeps) // limit the number of iterations ("sweeps")
{
// -- Initialization --
for (int i = 0; i < n; i++)
for (int j = i; j < n; j++) //copy mat[][] into M[][]
M[i][j] = mat[i][j]; //(M[][] is a local copy we can modify)
if (calc_evec)
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
evec[i][j] = (i==j) ? 1.0 : 0.0; //Set evec equal to the identity matrix
for (int i = 0; i < n-1; i++) //Initialize the "max_idx_row[]" array
max_idx_row[i] = MaxEntryRow(M, i); //(which is needed by MaxEntry())
// -- Iteration --
int n_iters;
int max_num_iters = max_num_sweeps*n*(n-1)/2; //"sweep" = n*(n-1)/2 iters
for (n_iters=1; n_iters <= max_num_iters; n_iters++) {
int i,j;
MaxEntry(M, i, j); // Find the maximum entry in the matrix. Store in i,j
// If M[i][j] is small compared to M[i][i] and M[j][j], set it to 0.
if ((M[i][i] + M[i][j] == M[i][i]) && (M[j][j] + M[i][j] == M[j][j])) {
M[i][j] = 0.0;
max_idx_row[i] = MaxEntryRow(M,i); //update max_idx_row[i]
}
if (M[i][j] == 0.0)
break;
// Otherwise, apply a rotation to make M[i][j] = 0
CalcRot(M, i, j); // Calculate the parameters of the rotation matrix.
ApplyRot(M, i, j); // Apply this rotation to the M matrix.
if (calc_evec) // Optional: If the caller wants the eigenvectors, then
ApplyRotLeft(evec,i,j); // apply the rotation to the eigenvector matrix
} //for (int n_iters=0; n_iters < max_num_iters; n_iters++)
// -- Post-processing --
for (int i = 0; i < n; i++)
eval[i] = M[i][i];
// Optional: Sort results by eigenvalue.
SortRows(eval, evec, n, sort_criteria);
if ((n_iters > max_num_iters) && (n>1)) // If we exceeded max_num_iters,
return 0; // indicate an error occured.
//assert(n_iters > 0);
return n_iters;
}
/// @brief Calculate the components of a rotation matrix which performs a
/// rotation in the i,j plane by an angle (θ) that (when multiplied on
/// both sides) will zero the ij'th element of M, so that afterwards
/// M[i][j] = 0. The results will be stored in c, s, and t
/// (which store cos(θ), sin(θ), and tan(θ), respectively).
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
CalcRot(Scalar const *const *M, // matrix
int i, // row index
int j) // column index
{
t = 1.0; // = tan(θ)
Scalar M_jj_ii = (M[j][j] - M[i][i]);
if (M_jj_ii != 0.0) {
// kappa = (M[j][j] - M[i][i]) / (2*M[i][j])
Scalar kappa = M_jj_ii;
t = 0.0;
Scalar M_ij = M[i][j];
if (M_ij != 0.0) {
kappa /= (2.0*M_ij);
// t satisfies: t^2 + 2*t*kappa - 1 = 0
// (choose the root which has the smaller absolute value)
t = 1.0 / (std::sqrt(1 + kappa*kappa) + std::abs(kappa));
if (kappa < 0.0)
t = -t;
}
}
//assert(std::abs(t) <= 1.0);
c = 1.0 / std::sqrt(1 + t*t);
s = c*t;
}
/// brief Perform a similarity transformation by multiplying matrix M on both
/// sides by a rotation matrix (and its transpose) to eliminate M[i][j].
/// details This rotation matrix performs a rotation in the i,j plane by
/// angle θ. This function assumes that c=cos(θ). s=sin(θ), t=tan(θ)
/// have been calculated in advance (using the CalcRot() function).
/// It also assumes that i<j. The max_idx_row[] array is also updated.
/// To save time, since the matrix is symmetric, the elements
/// below the diagonal (ie. M[u][v] where u>v) are not computed.
/// verbatim
/// M' = R^T * M * R
/// where R the rotation in the i,j plane and ^T denotes the transpose.
/// i j
/// _ _
/// | 1 |
/// | . |
/// | . |
/// | 1 |
/// | c ... s |
/// | . . . |
/// R = | . 1 . |
/// | . . . |
/// | -s ... c |
/// | 1 |
/// | . |
/// | . |
/// |_ 1 _|
/// endverbatim
///
/// Let M' denote the matrix M after multiplication by R^T and R.
/// The components of M' are:
///
/// verbatim
/// M'_uv = Σ_w Σ_z R_wu * M_wz * R_zv
/// endverbatim
///
/// Note that a the rotation at location i,j will modify all of the matrix
/// elements containing at least one index which is either i or j
/// such as: M[w][i], M[i][w], M[w][j], M[j][w].
/// Check and see whether these modified matrix elements exceed the
/// corresponding values in max_idx_row[] array for that row.
/// If so, then update max_idx_row for that row.
/// This is somewhat complicated by the fact that we must only consider
/// matrix elements in the upper-right triangle strictly above the diagonal.
/// (ie. matrix elements whose second index is > the first index).
/// The modified elements we must consider are marked with an "X" below:
///
/// @verbatim
/// i j
/// _ _
/// | . X X |
/// | . X X |
/// | . X X |
/// | . X X |
/// | X X X X X 0 X X X X | i
/// | . X |
/// | . X |
/// M = | . X |
/// | . X |
/// | X X X X X | j
/// | . |
/// | . |
/// | . |
/// |_ . _|
/// @endverbatim
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
ApplyRot(Scalar **M, // matrix
int i, // row index
int j) // column index
{
// Recall that:
// c = cos(θ)
// s = sin(θ)
// t = tan(θ) (which should be <= 1.0)
// Compute the diagonal elements of M which have changed:
M[i][i] -= t * M[i][j];
M[j][j] += t * M[i][j];
// Note: This is algebraically equivalent to:
// M[i][i] = c*c*M[i][i] + s*s*M[j][j] - 2*s*c*M[i][j]
// M[j][j] = s*s*M[i][i] + c*c*M[j][j] + 2*s*c*M[i][j]
//Update the off-diagonal elements of M which will change (above the diagonal)
//assert(i < j);
M[i][j] = 0.0;
//compute M[w][i] and M[i][w] for all w!=i,considering above-diagonal elements
for (int w=0; w < i; w++) { // 0 <= w < i < j < n
M[i][w] = M[w][i]; //backup the previous value. store below diagonal (i>w)
M[w][i] = c*M[w][i] - s*M[w][j]; //M[w][i], M[w][j] from previous iteration
if (i == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(M, w);
else if (std::abs(M[w][i])>std::abs(M[w][max_idx_row[w]])) max_idx_row[w]=i;
//assert(max_idx_row[w] == MaxEntryRow(M, w));
}
for (int w=i+1; w < j; w++) { // 0 <= i < w < j < n
M[w][i] = M[i][w]; //backup the previous value. store below diagonal (w>i)
M[i][w] = c*M[i][w] - s*M[w][j]; //M[i][w], M[w][j] from previous iteration
}
for (int w=j+1; w < n; w++) { // 0 <= i < j+1 <= w < n
M[w][i] = M[i][w]; //backup the previous value. store below diagonal (w>i)
M[i][w] = c*M[i][w] - s*M[j][w]; //M[i][w], M[j][w] from previous iteration
}
// now that we're done modifying row i, we can update max_idx_row[i]
max_idx_row[i] = MaxEntryRow(M, i);
//compute M[w][j] and M[j][w] for all w!=j,considering above-diagonal elements
for (int w=0; w < i; w++) { // 0 <= w < i < j < n
M[w][j] = s*M[i][w] + c*M[w][j]; //M[i][w], M[w][j] from previous iteration
if (j == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(M, w);
else if (std::abs(M[w][j])>std::abs(M[w][max_idx_row[w]])) max_idx_row[w]=j;
//assert(max_idx_row[w] == MaxEntryRow(M, w));
}
for (int w=i+1; w < j; w++) { // 0 <= i+1 <= w < j < n
M[w][j] = s*M[w][i] + c*M[w][j]; //M[w][i], M[w][j] from previous iteration
if (j == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(M, w);
else if (std::abs(M[w][j])>std::abs(M[w][max_idx_row[w]])) max_idx_row[w]=j;
//assert(max_idx_row[w] == MaxEntryRow(M, w));
}
for (int w=j+1; w < n; w++) { // 0 <= i < j < w < n
M[j][w] = s*M[w][i] + c*M[j][w]; //M[w][i], M[j][w] from previous iteration
}
// now that we're done modifying row j, we can update max_idx_row[j]
max_idx_row[j] = MaxEntryRow(M, j);
} //Jacobi::ApplyRot()
///@brief Multiply matrix M on the LEFT side by a transposed rotation matrix R^T
/// This matrix performs a rotation in the i,j plane by angle θ (where
/// the arguments "s" and "c" refer to cos(θ) and sin(θ), respectively).
/// @verbatim
/// E'_uv = Σ_w R_wu * E_wv
/// @endverbatim
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
ApplyRotLeft(Matrix E, // matrix
int i, // row index
int j) // column index
{
// Recall that c = cos(θ) and s = sin(θ)
for (int v = 0; v < n; v++) {
Scalar Eiv = E[i][v]; //backup E[i][v]
E[i][v] = c*E[i][v] - s*E[j][v];
E[j][v] = s*Eiv + c*E[j][v];
}
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
int Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
MaxEntryRow(Scalar const *const *M, int i) const {
int j_max = i+1;
for(int j = i+2; j < n; j++)
if (std::abs(M[i][j]) > std::abs(M[i][j_max]))
j_max = j;
return j_max;
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
MaxEntry(Scalar const *const *M, int& i_max, int& j_max) const {
// find the maximum entry in the matrix M in O(n) time
i_max = 0;
j_max = max_idx_row[i_max];
Scalar max_entry = std::abs(M[i_max][j_max]);
int nm1 = n-1;
for (int i=1; i < nm1; i++) {
int j = max_idx_row[i];
if (std::abs(M[i][j]) > max_entry) {
max_entry = std::abs(M[i][j]);
i_max = i;
j_max = j;
}
}
//#ifndef NDEBUG
make sure that the maximum element really is stored at i_max, j_max
(comment out the next loop before using. it slows down the code)
//for (int i = 0; i < nm1; i++)
// for (int j = i+1; j < n; j++)
// assert(std::abs(M[i][j]) <= max_entry);
//#endif
}
//Sort the rows of a matrix "evec" by the numbers contained in "eval"
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
SortRows(Vector eval, Matrix evec, int n, SortCriteria sort_criteria) const
{
for (int i = 0; i < n-1; i++) {
int i_max = i;
for (int j = i+1; j < n; j++) {
// find the "maximum" element in the array starting at position i+1
switch (sort_criteria) {
case SORT_DECREASING_EVALS:
if (eval[j] > eval[i_max])
i_max = j;
break;
case SORT_INCREASING_EVALS:
if (eval[j] < eval[i_max])
i_max = j;
break;
case SORT_DECREASING_ABS_EVALS:
if (std::abs(eval[j]) > std::abs(eval[i_max]))
i_max = j;
break;
case SORT_INCREASING_ABS_EVALS:
if (std::abs(eval[j]) < std::abs(eval[i_max]))
i_max = j;
break;
default:
break;
}
}
std::swap(eval[i], eval[i_max]); // sort "eval"
for (int k = 0; k < n; k++)
std::swap(evec[i][k], evec[i_max][k]); // sort "evec"
}
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Init() {
n = 0;
M = nullptr;
max_idx_row = nullptr;
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
SetSize(int n) {
Dealloc();
Alloc(n);
}
// memory management:
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Alloc(int N) {
n = N;
if (n > 0) {
max_idx_row = new int[n];
Alloc2D(n, n, &M);
}
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Dealloc() {
if (max_idx_row) {
delete [] max_idx_row;
max_idx_row = nullptr;
}
Dealloc2D(&M);
Init();
}
// memory management: copy and move constructor, swap, and assignment operator:
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Jacobi(const Jacobi<Scalar, Vector, Matrix, ConstMatrix>& source)
{
Init();
SetSize(source.n);
//assert(n == source.n);
// The following lines aren't really necessary, because the contents
// of source.M and source.max_idx_row are not needed (since they are
// overwritten every time Jacobi::Diagonalize() is invoked).
std::copy(source.max_idx_row,
source.max_idx_row + n,
max_idx_row);
for (int i = 0; i < n; i++)
std::copy(source.M[i],
source.M[i] + n,
M[i]);
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
swap(Jacobi<Scalar, Vector, Matrix, ConstMatrix> &other) noexcept {
std::swap(n, other.n);
std::swap(max_idx_row, other.max_idx_row);
std::swap(M, other.M);
}
// Move constructor (C++11)
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Jacobi(Jacobi<Scalar, Vector, Matrix, ConstMatrix>&& other)
: Jacobi() // delegate to the default constructor (which invokes Init())
{
this->swap(other);
}
// Using the "copy-swap" idiom for the assignment operator
// The assignment operator below implements both copy and move assignments.
// It accepts its argument by value, which invokes the appropriate
// (copy or move) constructor.
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
Jacobi<Scalar, Vector, Matrix, ConstMatrix>&
Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
operator = (Jacobi<Scalar, Vector, Matrix, ConstMatrix> source) {
this->swap(source);
return *this;
}
} // namespace jacobi
// THIS FILE USED TO BE EASY TO READ until I added "#if defined" statements.
// (They were added to test for many different kinds of array formats.)
using namespace matrix_alloc_jpd;
using namespace jacobi_pd;
// This code works with various types of C++ matrices (for example,
// double **, vector<vector<double>> array<array<double,5>,5>).
// I use "#if defined" statements to test different matrix types.
// For some of these (eg. array<array<double,5>,5>), the size of the matrix
// must be known at compile time. I specify that size now.
#if defined USE_ARRAY_OF_ARRAYS
const int NF=5; //(the array size must be known at compile time)
#elif defined USE_C_FIXED_SIZE_ARRAYS
const int NF=5; //(the array size must be known at compile time)
#endif
// @brief Are two numbers "similar"?
template<typename Scalar>
inline static bool Similar(Scalar a, Scalar b,
Scalar eps=1.0e-06,
Scalar ratio=1.0e-06,
Scalar ratio_denom=1.0)
{
return ((std::abs(a-b)<=std::abs(eps))
||
(std::abs(ratio_denom)*std::abs(a-b)
<=
std::abs(ratio)*0.5*(std::abs(a)+std::abs(b))));
}
/// @brief Are two vectors (containing n numbers) similar?
template<typename Scalar, typename Vector>
inline static bool SimilarVec(Vector a, Vector b, int n,
Scalar eps=1.0e-06,
Scalar ratio=1.0e-06,
Scalar ratio_denom=1.0)
{
for (int i = 0; i < n; i++)
if (not Similar(a[i], b[i], eps, ratio, ratio_denom))
return false;
return true;
}
/// @brief Are two vectors (or their reflections) similar?
template<typename Scalar, typename Vector>
inline static bool SimilarVecUnsigned(Vector a, Vector b, int n,
Scalar eps=1.0e-06,
Scalar ratio=1.0e-06,
Scalar ratio_denom=1.0)
{
if (SimilarVec(a, b, n, eps))
return true;
else {
for (int i = 0; i < n; i++)
if (not Similar(a[i], -b[i], eps, ratio, ratio_denom))
return false;
return true;
}
}
/// @brief Multiply two matrices A and B, store the result in C. (C = AB).
template<typename Matrix, typename ConstMatrix>
void mmult(ConstMatrix A, //<! input array
ConstMatrix B, //<! input array
Matrix C, //<! store result here
int m, //<! number of rows of A
int n=0, //<! optional: number of columns of B (=m by default)
int K=0 //<! optional: number of columns of A = num rows of B (=m by default)
)
{
if (n == 0) n = m; // if not specified, then assume the matrices are square
if (K == 0) K = m; // if not specified, then assume the matrices are square
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
C[i][j] = 0.0;
// perform matrix multiplication
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
for (int k = 0; k < K; k++)
C[i][j] += A[i][k] * B[k][j];
}
/// @brief
///Sort the rows of a matrix "evec" by the numbers contained in "eval"
///(This is a simple O(n^2) sorting method, but O(n^2) is a lower bound anyway.)
///This is the same as the Jacobi::SortRows(), but that function is private.
template<typename Scalar, typename Vector, typename Matrix>
void
SortRows(Vector eval,
Matrix evec,
int n,
bool sort_decreasing=true,
bool sort_abs=false)
{
for (int i = 0; i < n-1; i++) {
int i_max = i;
for (int j = i+1; j < n; j++) {
if (sort_decreasing) {
if (sort_abs) { //sort by absolute value?
if (std::abs(eval[j]) > std::abs(eval[i_max]))
i_max = j;
}
else if (eval[j] > eval[i_max])
i_max = j;
}
else {
if (sort_abs) { //sort by absolute value?
if (std::abs(eval[j]) < std::abs(eval[i_max]))
i_max = j;
}
else if (eval[j] < eval[i_max])
i_max = j;
}
}
std::swap(eval[i], eval[i_max]); // sort "eval"
for (int k = 0; k < n; k++)
std::swap(evec[i][k], evec[i_max][k]); // sort "evec"
}
}
/// @brief Generate a random orthonormal n x n matrix
template<typename Scalar, typename Matrix>
void GenRandOrth(Matrix R,
int n,
std::default_random_engine &rand_generator)
{
std::normal_distribution<Scalar> gaussian_distribution(0,1);
std::vector<Scalar> v(n);
for (int i = 0; i < n; i++) {
// Generate a vector, "v", in a random direction subject to the constraint
// that it is orthogonal to the first i-1 rows-vectors of the R matrix.
Scalar rsq = 0.0;
while (rsq == 0.0) {
// Generate a vector in a random direction
// (This works because we are using a normal (Gaussian) distribution)
for (int j = 0; j < n; j++)
v[j] = gaussian_distribution(rand_generator);
//Now subtract from v, the projection of v onto the first i-1 rows of R.
//This will produce a vector which is orthogonal to these i-1 row-vectors.
//(They are already normalized and orthogonal to each other.)
for (int k = 0; k < i; k++) {
Scalar v_dot_Rk = 0.0;
for (int j = 0; j < n; j++)
v_dot_Rk += v[j] * R[k][j];
for (int j = 0; j < n; j++)
v[j] -= v_dot_Rk * R[k][j];
}
// check if it is linearly independent of the other vectors and non-zero
rsq = 0.0;
for (int j = 0; j < n; j++)
rsq += v[j]*v[j];
}
// Now normalize the vector
Scalar r_inv = 1.0 / std::sqrt(rsq);
for (int j = 0; j < n; j++)
v[j] *= r_inv;
// Now copy this vector to the i'th row of R
for (int j = 0; j < n; j++)
R[i][j] = v[j];
} //for (int i = 0; i < n; i++)
} //void GenRandOrth()
/// @brief Generate a random symmetric n x n matrix, M.
/// This function generates random numbers for the eigenvalues ("evals_known")
/// as well as the eigenvectors ("evecs_known"), and uses them to generate M.
/// The "eval_magnitude_range" argument specifies the the base-10 logarithm
/// of the range of eigenvalues desired. The "n_degeneracy" argument specifies
/// the number of repeated eigenvalues desired (if any).
/// @returns This function does not return a value. However after it is
/// invoked, the M matrix will be filled with random numbers.
/// Additionally, the "evals" and "evecs" arguments will contain
/// the eigenvalues and eigenvectors (one eigenvector per row)
/// of the matrix. Later, they can be compared with the eigenvalues
/// and eigenvectors calculated by Jacobi::Diagonalize()
template <typename Scalar, typename Vector, typename Matrix>
void GenRandSymm(Matrix M, //<! store the matrix here
int n, //<! matrix size
Vector evals, //<! store the eigenvalues of here
Matrix evecs, //<! store the eigenvectors here
std::default_random_engine &rand_generator,//<! makes random numbers
Scalar min_eval_size=0.1, //<! minimum possible eigenvalue size
Scalar max_eval_size=10.0,//<! maximum possible eigenvalue size
int n_degeneracy=1//<!number of repeated eigevalues(1disables)
)
{
assert(n_degeneracy <= n);
std::uniform_real_distribution<Scalar> random_real01;
std::normal_distribution<Scalar> gaussian_distribution(0, max_eval_size);
bool use_log_uniform_distribution = false;
if (min_eval_size > 0.0)
use_log_uniform_distribution = true;
#if defined USE_VECTOR_OF_VECTORS
vector<vector<Scalar> > D(n, vector<Scalar>(n));
vector<vector<Scalar> > tmp(n, vector<Scalar>(n));
#elif defined USE_ARRAY_OF_ARRAYS
array<array<Scalar, NF>, NF> D;
array<array<Scalar, NF>, NF> tmp;
#elif defined USE_C_FIXED_SIZE_ARRAYS
Scalar D[NF][NF], tmp[NF][NF];
#else
#define USE_C_POINTER_TO_POINTERS
Scalar **D, **tmp;
Alloc2D(n, n, &D);
Alloc2D(n, n, &tmp);
#endif
// Randomly generate the eigenvalues
for (int i = 0; i < n; i++) {
if (use_log_uniform_distribution) {
// Use a "log-uniform distribution" (a.k.a. "reciprocal distribution")
// (This is a way to specify numbers with a precise range of magnitudes.)
assert((min_eval_size > 0.0) && (max_eval_size > 0.0));
Scalar log_min = std::log(std::abs(min_eval_size));
Scalar log_max = std::log(std::abs(max_eval_size));
Scalar log_eval = (log_min + random_real01(rand_generator)*(log_max-log_min));
evals[i] = std::exp(log_eval);
// also consider both positive and negative eigenvalues:
if (random_real01(rand_generator) < 0.5)
evals[i] = -evals[i];
}
else {
evals[i] = gaussian_distribution(rand_generator);
}
}
// Does the user want us to force some of the eigenvalues to be the same?
if (n_degeneracy > 1) {
int *permutation = new int[n]; //a random permutation from 0...n-1
for (int i = 0; i < n; i++)
permutation[i] = i;
std::shuffle(permutation, permutation+n, rand_generator);
for (int i = 1; i < n_degeneracy; i++) //set the first n_degeneracy to same
evals[permutation[i]] = evals[permutation[0]];
delete [] permutation;
}
// D is a diagonal matrix whose diagonal elements are the eigenvalues
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
D[i][j] = ((i == j) ? evals[i] : 0.0);
// Now randomly generate the (transpose of) the "evecs" matrix
GenRandOrth<Scalar, Matrix>(evecs, n, rand_generator); //(will transpose it later)
// Construct the test matrix, M, where M = Rt * D * R
// Original code:
//mmult(evecs, D, tmp, n); // <--> tmp = Rt * D
// Unfortunately, C++ guesses the types incorrectly. Must manually specify:
// #ifdefs making the code ugly again:
#if defined USE_VECTOR_OF_VECTORS
mmult<vector<vector<Scalar> >&, const vector<vector<Scalar> >&>
#elif defined USE_ARRAY_OF_ARRAYS
mmult<array<array<Scalar,NF>,NF>&, const array<array<Scalar,NF>,NF>&>
#elif defined USE_C_FIXED_SIZE_ARRAYS
mmult<Scalar (*)[NF], Scalar (*)[NF]>
#else
mmult<Scalar**, Scalar const *const *>
#endif
(evecs, D, tmp, n);
for (int i = 0; i < n-1; i++)
for (int j = i+1; j < n; j++)
std::swap(evecs[i][j], evecs[j][i]); //transpose "evecs"
// Original code:
//mmult(tmp, evecs, M, n);
// Unfortunately, C++ guesses the types incorrectly. Must manually specify:
// #ifdefs making the code ugly again:
#if defined USE_VECTOR_OF_VECTORS
mmult<vector<vector<Scalar> >&, const vector<vector<Scalar> >&>
#elif defined USE_ARRAY_OF_ARRAYS
mmult<array<array<Scalar,NF>,NF>&, const array<array<Scalar,NF>,NF>&>
#elif defined USE_C_FIXED_SIZE_ARRAYS
mmult<Scalar (*)[NF], Scalar (*)[NF]>
#else
mmult<Scalar**, Scalar const *const *>
#endif
(tmp, evecs, M, n);
//at this point M = Rt*D*R (where "R"="evecs")
#if defined USE_C_POINTER_TO_POINTERS
Dealloc2D(&D);
Dealloc2D(&tmp);
#endif
} // GenRandSymm()
template <typename Scalar>
void TestJacobi(int n, //<! matrix size
int n_matrices=100, //<! number of matrices to test
Scalar min_eval_size=0.1, //<! minimum possible eigenvalue sizw
Scalar max_eval_size=10.0, //<! maximum possible eigenvalue size
int n_tests_per_matrix=1, //<! repeat test for benchmarking?
int n_degeneracy=1, //<! repeated eigenvalues?
unsigned seed=0, //<! random seed (if 0 then use the clock)
Scalar eps=1.0e-06
)
{
bool test_code_coverage = false;
if (n_tests_per_matrix < 1) {
cout << "-- Testing code-coverage --" << endl;
test_code_coverage = true;
n_tests_per_matrix = 1;
}
cout << endl << "-- Diagonalization test (real symmetric) --" << endl;
// construct a random generator engine using a time-based seed:
if (seed == 0) // if the caller did not specify a seed, use the system clock
seed = std::chrono::system_clock::now().time_since_epoch().count();
std::default_random_engine rand_generator(seed);
// Create an instance of the Jacobi diagonalizer, and allocate the matrix
// we will test it on, as well as the arrays that will store the resulting
// eigenvalues and eigenvectors.
// The way we do this depends on what version of the code we are using.
// This is controlled by "#if defined" statements.
#if defined USE_VECTOR_OF_VECTORS
Jacobi<Scalar,
vector<Scalar>&,
vector<vector<Scalar> >&,
const vector<vector<Scalar> >& > ecalc(n);
// allocate the matrix, eigenvalues, eigenvectors
vector<vector<Scalar> > M(n, vector<Scalar>(n));
vector<vector<Scalar> > evecs(n, vector<Scalar>(n));
vector<vector<Scalar> > evecs_known(n, vector<Scalar>(n));
vector<Scalar> evals(n);
vector<Scalar> evals_known(n);
vector<Scalar> test_evec(n);
#elif defined USE_ARRAY_OF_ARRAYS
n = NF;
cout << "Testing std::array (fixed size).\n"
"(Ignoring first argument, and setting matrix size to " << n << ")" << endl;
Jacobi<Scalar,
array<Scalar, NF>&,
array<array<Scalar, NF>, NF>&,
const array<array<Scalar, NF>, NF>&> ecalc(n);
// allocate the matrix, eigenvalues, eigenvectors
array<array<Scalar, NF>, NF> M;
array<array<Scalar, NF>, NF> evecs;
array<array<Scalar, NF>, NF> evecs_known;
array<Scalar, NF> evals;
array<Scalar, NF> evals_known;
array<Scalar, NF> test_evec;
#elif defined USE_C_FIXED_SIZE_ARRAYS
n = NF;
cout << "Testing C fixed size arrays.\n"
"(Ignoring first argument, and setting matrix size to " << n << ")" << endl;
Jacobi<Scalar,
Scalar*,
Scalar (*)[NF],
Scalar const (*)[NF]> ecalc(n);
// allocate the matrix, eigenvalues, eigenvectors
Scalar M[NF][NF];
Scalar evecs[NF][NF];
Scalar evecs_known[NF][NF];
Scalar evals[NF];
Scalar evals_known[NF];
Scalar test_evec[NF];
#else
#define USE_C_POINTER_TO_POINTERS
// Note: Normally, you would just use this to instantiate Jacobi:
// Jacobi<Scalar, Scalar*, Scalar**, Scalar const*const*> ecalc(n);
// -------------------------
// ..but since Jacobi manages its own memory using new and delete, I also want
// to test that the copy constructors, copy operators, and destructors work.
// The following lines do this:
Jacobi<Scalar, Scalar*, Scalar**, Scalar const*const*> ecalc_test_mem1;
ecalc_test_mem1.SetSize(n);
Jacobi<Scalar, Scalar*, Scalar**, Scalar const*const*> ecalc_test_mem2(2);
// test the = operator
ecalc_test_mem2 = ecalc_test_mem1;
// test the copy constructor
Jacobi<Scalar, Scalar*, Scalar**, Scalar const*const*> ecalc(ecalc_test_mem2);
// allocate the matrix, eigenvalues, eigenvectors
Scalar **M, **evecs, **evecs_known;
Alloc2D(n, n, &M);
Alloc2D(n, n, &evecs);
Alloc2D(n, n, &evecs_known);
Scalar *evals = new Scalar[n];
Scalar *evals_known = new Scalar[n];
Scalar *test_evec = new Scalar[n];
#endif
// --------------------------------------------------------------------
// Now, generate random matrices and test Jacobi::Diagonalize() on them.
// --------------------------------------------------------------------
for(int imat = 0; imat < n_matrices; imat++) {
// Create a randomly generated symmetric matrix.
//This function generates random numbers for the eigenvalues ("evals_known")
//as well as the eigenvectors ("evecs_known"), and uses them to generate M.
#if defined USE_VECTOR_OF_VECTORS
GenRandSymm<Scalar, vector<Scalar>&, vector<vector<Scalar> >&>
#elif defined USE_ARRAY_OF_ARRAYS
GenRandSymm<Scalar, array<Scalar,NF>&, array<array<Scalar,NF>,NF>&>
#elif defined USE_C_FIXED_SIZE_ARRAYS
GenRandSymm<Scalar, Scalar*, Scalar (*)[NF]>
#else
GenRandSymm<Scalar, Scalar*, Scalar**>
#endif
(M,
n,
evals_known,
evecs_known,
rand_generator,
min_eval_size,
max_eval_size,
n_degeneracy);
// Sort the matrix evals and eigenvector rows:
// Original code:
//SortRows<Scalar>(evals_known, evecs_known, n);
// Unfortunately, C++ guesses the types incorrectly. Must use #ifdefs again:
#if defined USE_VECTOR_OF_VECTORS
SortRows<Scalar, vector<Scalar>&, vector<vector<Scalar> >&>
#elif defined USE_ARRAY_OF_ARRAYS
SortRows<Scalar, array<Scalar,NF>&, array<array<Scalar,NF>,NF>&>
#elif defined USE_C_FIXED_SIZE_ARRAYS
SortRows<Scalar, Scalar*, Scalar (*)[NF]>
#else
SortRows<Scalar, Scalar*, Scalar**>
#endif
(evals_known, evecs_known, n);
if (n_matrices == 1) {
cout << "Matrix:\n";
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
cout << M[i][j] << " ";
cout << "\n";
}
cout << "Eigenvalues (after sorting):\n";
for (int i = 0; i < n; i++)
cout << evals_known[i] << " ";
cout << "\n";
cout << "Eigenvectors (rows) which are known in advance:\n";
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
cout << evecs_known[i][j] << " ";
cout << "\n";
}
cout << " (The eigenvectors calculated by Jacobi::Diagonalize() should match these.)\n";
}
for (int i_test = 0; i_test < n_tests_per_matrix; i_test++) {
if (test_code_coverage) {
// test SORT_INCREASING_ABS_EVALS:
#if defined USE_VECTOR_OF_VECTORS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
vector<Scalar>&,
vector<vector<Scalar> >&,
const vector<vector<Scalar> >& >::SORT_INCREASING_ABS_EVALS);
#elif defined USE_ARRAY_OF_ARRAYS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
array<Scalar,NF>&,
array<array<Scalar,NF>,NF>&,
const array<array<Scalar,NF>,NF>&>::SORT_INCREASING_ABS_EVALS);
#elif defined USE_C_FIXED_SIZE_ARRAYS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
Scalar*,
Scalar (*)[NF],
Scalar const (*)[NF]>::SORT_INCREASING_ABS_EVALS);
#else
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
Scalar*,
Scalar**,
Scalar const*const*>::SORT_INCREASING_ABS_EVALS);
// Also make sure the code considers a scenario where convergence fails:
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
Scalar*,
Scalar**,
Scalar const*const*>::SORT_INCREASING_ABS_EVALS,
true,
0); //<-- set the maximum allowed iterations to 0
#endif
for (int i = 1; i < n; i++)
assert(std::abs(evals[i-1])<=std::abs(evals[i]));
// test SORT_DECREASING_ABS_EVALS:
#if defined USE_VECTOR_OF_VECTORS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
vector<Scalar>&,
vector<vector<Scalar> >&,
const vector<vector<Scalar> >& >::SORT_DECREASING_ABS_EVALS);
#elif defined USE_ARRAY_OF_ARRAYS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
array<Scalar,NF>&,
array<array<Scalar,NF>,NF>&,
const array<array<Scalar,NF>,NF>&>::SORT_DECREASING_ABS_EVALS);
#elif defined USE_C_FIXED_SIZE_ARRAYS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
Scalar*,
Scalar (*)[NF],
Scalar const (*)[NF]>::SORT_DECREASING_ABS_EVALS);
#else
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
Scalar*,
Scalar**,
Scalar const*const*>::SORT_DECREASING_ABS_EVALS);
#endif
for (int i = 1; i < n; i++)
assert(std::abs(evals[i-1])>=std::abs(evals[i]));
// test SORT_INCREASING_EVALS:
#if defined USE_VECTOR_OF_VECTORS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
vector<Scalar>&,
vector<vector<Scalar> >&,
const vector<vector<Scalar> >& >::SORT_INCREASING_EVALS);
#elif defined USE_ARRAY_OF_ARRAYS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
array<Scalar,NF>&,
array<array<Scalar,NF>,NF>&,
const array<array<Scalar,NF>,NF>&>::SORT_INCREASING_EVALS);
#elif defined USE_C_FIXED_SIZE_ARRAYS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
Scalar*,
Scalar (*)[NF],
Scalar const (*)[NF]>::SORT_INCREASING_EVALS);
#else
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
Scalar*,
Scalar**,
Scalar const*const*>::SORT_INCREASING_EVALS);
#endif
for (int i = 1; i < n; i++)
assert(evals[i-1] <= evals[i]);
// test DO_NOT_SORT
#if defined USE_VECTOR_OF_VECTORS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
vector<Scalar>&,
vector<vector<Scalar> >&,
const vector<vector<Scalar> >& >::DO_NOT_SORT);
#elif defined USE_ARRAY_OF_ARRAYS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
array<Scalar,NF>&,
array<array<Scalar,NF>,NF>&,
const array<array<Scalar,NF>,NF>&>::DO_NOT_SORT);
#elif defined USE_C_FIXED_SIZE_ARRAYS
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
Scalar*,
Scalar (*)[NF],
Scalar const (*)[NF]>::DO_NOT_SORT);
#else
ecalc.Diagonalize(M,
evals,
evecs,
Jacobi<Scalar,
Scalar*,
Scalar**,
Scalar const*const*>::DO_NOT_SORT);
#endif
} //if (test_code_coverage)
// Now (finally) calculate the eigenvalues and eigenvectors
int n_sweeps = ecalc.Diagonalize(M, evals, evecs);
if ((n_matrices == 1) && (i_test == 0)) {
cout <<"Jacobi::Diagonalize() ran for "<<n_sweeps<<" iters (sweeps).\n";
cout << "Eigenvalues calculated by Jacobi::Diagonalize()\n";
for (int i = 0; i < n; i++)
cout << evals[i] << " ";
cout << "\n";
cout << "Eigenvectors (rows) calculated by Jacobi::Diagonalize()\n";
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
cout << evecs[i][j] << " ";
cout << "\n";
}
}
assert(SimilarVec(evals, evals_known, n, eps*max_eval_size, eps));
//Check that each eigenvector satisfies Mv = λv
// <--> Σ_b M[a][b]*evecs[i][b] = evals[i]*evecs[i][b] (for all a)
for (int i = 0; i < n; i++) {
for (int a = 0; a < n; a++) {
test_evec[a] = 0.0;
for (int b = 0; b < n; b++)
test_evec[a] += M[a][b] * evecs[i][b];
assert(Similar(test_evec[a],
evals[i] * evecs[i][a],
eps, // tolerance (absolute difference)
eps*max_eval_size, // tolerance ratio (numerator)
evals_known[i] // tolerance ration (denominator)
));
}
}
} //for (int i_test = 0; i_test < n_tests_per_matrix; i++)
} //for(int imat = 0; imat < n_matrices; imat++) {
#if defined USE_C_POINTER_TO_POINTERS
Dealloc2D(&M);
Dealloc2D(&evecs);
Dealloc2D(&evecs_known);
delete [] evals;
delete [] evals_known;
delete [] test_evec;
#endif
} //TestJacobi()
int main(int argc, char **argv) {
int n_size = 2;
int n_matr = 1;
double emin = 0.0;
double emax = 1.0;
int n_tests = 1;
int n_degeneracy = 1;
unsigned seed = 0;
if (argc <= 1) {
cerr <<
"Error: This program requires at least 1 argument.\n"
"\n"
"Description: Run Jacobi::Diagonalize() on randomly generated matrices.\n"
"\n"
"Arguments: n_size [n_matr emin emax n_degeneracy n_tests seed eps]\n"
" n_size = the size of the matrices\n"
" (NOTE: The remaining arguments are optional.)\n"
" n_matr = the number of randomly generated matrices to test\n"
" emin = the smallest possible eigenvalue magnitude (eg. 1e-05)\n"
" emax = the largest possible eigenvalue magnitude (>0 eg. 1e+05)\n"
" (NOTE: If emin=0, a normal distribution is used centered at 0.\n"
" Otherwise a log-uniform distribution is used from emin to emax.)\n"
" n_degeneracy = the number of repeated eigenvalues (1 disables, default)\n"
" n_tests = the number of times the eigenvalues and eigenvectors\n"
" are calculated for EACH matrix. By default this is 1.\n"
" (Increase this to at least 20 if you plan to use this\n"
" program for benchmarking (speed testing), because the time\n"
" needed for generating a random matrix is not negligible.)\n"
" (IF THIS NUMBER IS 0, it will test CODE-COVERAGE instead.)\n"
" seed = the seed used by the random number \"rand_generator\".\n"
" (If this number is 0, which is the default, the system\n"
" clock is used to choose a random seed.)\n"
" eps = the tolerance. The difference between eigenvalues and their\n"
" true value, cannot exceed this (multiplied by the eigenvalue\n"
" of maximum magnitude). Similarly, the difference between\n"
" the eigenvectors after multiplication by the matrix and by\n"
" and after multiplication by the eigenvalue, cannot exceed\n"
" eps*maximum_eigenvalue/eigenvalue. The default value is\n"
" 1.0e-06 (which works well for double precision numbers).\n"
<< endl;
return 1;
}
n_size = std::stoi(argv[1]);
if (argc > 2)
n_matr = std::stoi(argv[2]);
if (argc > 3)
emin = std::stof(argv[3]);
if (argc > 4)
emax = std::stof(argv[4]);
if (argc > 5)
n_degeneracy = std::stoi(argv[5]);
if (argc > 6)
n_tests = std::stoi(argv[6]);
if (argc > 7)
seed = std::stoi(argv[7]);
double eps = 1.0e-06;
if (argc > 8)
eps = std::stof(argv[8]);
TestJacobi(n_size, n_matr, emin, emax, n_tests, n_degeneracy, seed, eps);
cout << "test passed\n" << endl;
return EXIT_SUCCESS;
}
1525
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
