songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
cgal/Kernel_d/test/Kernel_d/Linear_algebra-test.cpp

416 lines
12 KiB
C++
Raw Blame History

#include <CGAL/Linear_algebraHd.h>
#include <CGAL/Linear_algebraCd.h>
#include <CGAL/random_selection.h>
#include <CGAL/test_macros.h>
#include <CGAL/double.h>
#include <cstdlib>
#define VEC_DIM 10
#define MAT_DIM 10
#include <CGAL/Exact_integer.h>
#include <CGAL/Exact_rational.h>
typedef CGAL::Exact_integer RT;
typedef CGAL::Exact_rational FT;
int main(int argc, char* argv[])
{
CGAL_KD_SETDTHREAD(151);
CGAL::IO::set_pretty_mode ( std::cerr );
CGAL_TEST_START;
{
typedef RT NT;
typedef CGAL::Linear_algebraHd<RT> LA;
typedef LA::Matrix Matrix;
typedef LA::Vector Vector;
bool IOTEST = true;
{
int vec_dim;
if (argc == 2) vec_dim = std::atoi(argv[1]);
else vec_dim = VEC_DIM;
/* some construction and access ops */
Vector v0(vec_dim), v1(vec_dim), v2(vec_dim);
int F[] = { 1,2,3,4,5 };
Vector v11(F,F+2), v12(F,F+3), v13(F,F+4), v14(F,F+5),
v15(v13.begin(),v13.end()),
v16(vec_dim,NT(1));
CGAL_TEST(v13==v15){}
CGAL_TEST(v0==v1){}
for (int i = 0; i < vec_dim; i++) {
v1[i] = i;
v2[i] = vec_dim-i;
}
CGAL_TEST(v0!=v1){}
Vector v3(v1);
CGAL_TEST(v3==v1){}
v1 += NT(2)*v2;
v1 -= v2;
v1 = -v1;
CGAL_TEST((v1*v1 == NT(vec_dim*vec_dim*vec_dim))){}
/* squared length of (v1 + v2) = v1*v1
should be equal to dim*dim*dim */
NT res1 = (v1*v1 - v2*v2);
NT res2 = ((v1-v2)*(v1+v2));
CGAL_TEST(res1==res2){}
/* we test the third binomial formula */
Vector v21(v13),v22(v13);
v21 *= 13;
CGAL_TEST(v21 == NT(13)*v22){}
if (IOTEST) CGAL_IO_TEST(v1,v2,CGAL::IO::ASCII);
}
{
int mat_dim;
if (argc == 2) mat_dim = std::atoi(argv[1]);
else mat_dim = MAT_DIM;
/* some construction and access ops */
Matrix::Identity ID;
Matrix A(mat_dim,mat_dim), B(mat_dim),
I(mat_dim,ID), One(mat_dim,mat_dim,NT(2));
CGAL_TEST(A==B){}
CGAL_TEST(One*I==One){}
std::vector<Vector> F(mat_dim);
int i,j;
for (i = 0; i < mat_dim; i++) {
Vector v(mat_dim);
for (j = 0; j < mat_dim; j++) {
A(i,j) = i;
B(i,j) = j;
v[j] = j;
}
F[i] = v;
}
#ifndef CGAL_SIMPLE_INTERFACE
Matrix C(F), D(A), E(F.begin(),F.end());
#else
Matrix C(F), D(A), E(F);
#endif
CGAL_TEST(C==F && E==F){}
/* A = (0,0,0,0...)
(1,1,1,1...)
(2,2,2,2...)
...
B = (0,1,2,3...)
(0,1,2,3...)
(0,1,2,3...)
...
C = A = D = E;
*/
Matrix::iterator it;
for (it = A.begin(), i = 0; it != A.end(); ++i,++it) {
CGAL_TEST(A( i/mat_dim, i%mat_dim ) == *it){}
}
CGAL_TEST(A==C){}
CGAL_TEST(A!=B){}
CGAL_TEST(A==D){}
CGAL_TEST(A.row_dimension()==mat_dim && A.column_dimension()==mat_dim){}
CGAL_TEST(A.row(1)*B.column(1)==NT(mat_dim)){}
/* some basic arithmetic testing */
C = A; C += A;
C -= NT(3)*A;
C = -C;
CGAL_TEST(C==A){}
/* row sum test: */
Vector ones(mat_dim), row_sum_vec(mat_dim);
for (i = 0; i < mat_dim; i++) {
ones[i] = 1;
row_sum_vec[i] = i*mat_dim;
}
CGAL_TEST(A*ones==row_sum_vec){}
C = A+B;
/* matrix operations + ,* and |transpose|, |rank| */
CGAL_TEST(C==LA::transpose(C)){}
C = C-B;
CGAL_TEST(C==A){}
C = I+A;
CGAL_TEST(LA::rank(C)==mat_dim){}
/* matrix operations 2* , |determinant| */
C = NT(2) * I;
C = C * NT(2);
Matrix L,U;
Vector c;
std::vector<int> q;
NT det = LA::determinant(C, L, U, q, c); // det must be 2^{2*mat_dim}
NT pot = 1; for (i=0; i<mat_dim; ++i) pot *= NT(4);
CGAL_TEST(det == pot){}
CGAL_TEST(LA::verify_determinant(C, det, L, U, q, c)){}
CGAL_TEST(det == LA::determinant(C)){}
CGAL_TEST(CGAL_NTS sign(det) == LA::sign_of_determinant(C)){}
if (IOTEST) CGAL_IO_TEST(A,C,CGAL::IO::ASCII);
/* a random linear solver task: */
Vector b(mat_dim),x(mat_dim),e;
NT denom;
for (i = 0; i < mat_dim; i++) {
for (j = 0; j < mat_dim; j++)
C(i,j) = CGAL::get_default_random().get_int(-mat_dim,mat_dim);
b[i] = CGAL::get_default_random().get_int(-mat_dim,mat_dim);
}
for (double f = 1.0; f > 0.0; f-=0.1) {
Matrix E = C;
for (i = 0; i< mat_dim; ++i)
for (j = 0; j < mat_dim; ++j)
if ( CGAL::get_default_random().get_double() > f ) E(i,j)=0;
if (LA::linear_solver(E, b, x, denom, A, e)) {
CGAL_TEST(E*x == b*denom){}
LA::linear_solver(E,b,x,denom,e);
CGAL_TEST(E*x == b*denom){}
LA::linear_solver(E,b,x,denom);
CGAL_TEST(E*x == b*denom){}
CGAL_TEST(LA::is_solvable(E,b)){}
}
else {
Vector null(mat_dim);
CGAL_TEST(LA::transpose(E)*e==null && b*e!=NT(0)){}
CGAL_TEST(!LA::is_solvable(E,b)){}
}
Matrix SV;
if (LA::homogeneous_linear_solver(E,x)) {
CGAL_TEST(E*x==Vector(mat_dim)){}
int r = LA::homogeneous_linear_solver(E,SV);
CGAL_TEST(r==LA::rank(SV)){}
}
// inverse
if (LA::inverse(E,D,denom,c)) {
CGAL_TEST(E*D==denom*Matrix(mat_dim,Matrix::Identity())){}
CGAL_TEST(D==LA::inverse(E,denom)){}
} else {
CGAL_TEST(LA::transpose(E)*c==Vector(mat_dim)){}
}
LA::independent_columns(E,q);
}
if (mat_dim > 1) { // ueberbestimmt:
Matrix N(mat_dim,mat_dim-1);
for (i = 0; i < mat_dim; i++) {
for (j = 0; j < mat_dim-1; j++)
N(i,j) = CGAL::get_default_random().get_int(-mat_dim,mat_dim);
b[i] = CGAL::get_default_random().get_int(-mat_dim,mat_dim);
}
if (LA::linear_solver(N,b,x,denom,e)) {
CGAL_TEST(N*x == b*denom){}
} else {
Vector null(mat_dim-1);
CGAL_TEST(LA::transpose(N)*e==null && b*e!=NT(0)){}
}
}
}
}
{
typedef FT NT;
typedef CGAL::Linear_algebraCd<FT> LA;
typedef LA::Matrix Matrix;
typedef LA::Vector Vector;
bool IOTEST = false;
{
int vec_dim;
if (argc == 2) vec_dim = std::atoi(argv[1]);
else vec_dim = VEC_DIM;
/* some construction and access ops */
Vector v0(vec_dim), v1(vec_dim), v2(vec_dim);
int F[] = { 1,2,3,4,5 };
Vector v11(F,F+2), v12(F,F+3), v13(F,F+4), v14(F,F+5),
v15(v13.begin(),v13.end()),
v16(vec_dim,NT(1));
CGAL_TEST(v13==v15){}
CGAL_TEST(v0==v1){}
for (int i = 0; i < vec_dim; i++) {
v1[i] = i;
v2[i] = vec_dim-i;
}
CGAL_TEST(v0!=v1){}
Vector v3(v1);
CGAL_TEST(v3==v1){}
v1 += NT(2)*v2;
v1 -= v2;
v1 = -v1;
CGAL_TEST((v1*v1 == NT(vec_dim*vec_dim*vec_dim))){}
/* squared length of (v1 + v2) = v1*v1
should be equal to dim*dim*dim */
NT res1 = (v1*v1 - v2*v2);
NT res2 = ((v1-v2)*(v1+v2));
CGAL_TEST(res1==res2){}
/* we test the third binomial formula */
Vector v21(v13),v22(v13);
v21 *= 13;
CGAL_TEST(v21 == NT(13)*v22){}
v21 /= 13;
CGAL_TEST(v21 == v22){}
if (IOTEST) CGAL_IO_TEST(v1,v2,CGAL::IO::ASCII);
}
{
int mat_dim;
if (argc == 2) mat_dim = std::atoi(argv[1]);
else mat_dim = MAT_DIM;
/* some construction and access ops */
Matrix::Identity ID;
Matrix A(mat_dim,mat_dim), B(mat_dim),
I(mat_dim,ID), One(mat_dim,mat_dim,NT(2));
CGAL_TEST(A==B){}
CGAL_TEST(One*I==One){}
std::vector<Vector> F(mat_dim);
int i,j;
for (i = 0; i < mat_dim; i++) {
Vector v(mat_dim);
for (j = 0; j < mat_dim; j++) {
A(i,j) = i;
B(i,j) = j;
v[j] = j;
}
F[i] = v;
}
#ifndef CGAL_SIMPLE_INTERFACE
Matrix C(F), D(A), E(F.begin(),F.end());
#else
Matrix C(F), D(A), E(F);
#endif
CGAL_TEST(C==F && E==F){}
/* A = (0,0,0,0...)
(1,1,1,1...)
(2,2,2,2...)
...
B = (0,1,2,3...)
(0,1,2,3...)
(0,1,2,3...)
...
C = A = D = E;
*/
Matrix::iterator it;
for (it = A.begin(), i = 0; it != A.end(); ++i,++it) {
CGAL_TEST(A( i/mat_dim, i%mat_dim ) == *it){}
}
CGAL_TEST(A==C){}
CGAL_TEST(A!=B){}
CGAL_TEST(A==D){}
CGAL_TEST(A.row_dimension()==mat_dim && A.column_dimension()==mat_dim){}
CGAL_TEST(A.row(1)*B.column(1)==NT(mat_dim)){}
/* some basic arithmetic testing */
C = A; C += A;
C -= NT(3)*A;
C = -C;
CGAL_TEST(C==A){}
/* row sum test: */
Vector ones(mat_dim), row_sum_vec(mat_dim);
for (i = 0; i < mat_dim; i++) {
ones[i] = 1;
row_sum_vec[i] = i*mat_dim;
}
CGAL_TEST(A*ones==row_sum_vec){}
C = A+B;
/* matrix operations + ,* and |transpose|, |rank| */
CGAL_TEST(C==LA::transpose(C)){}
C = C-B;
CGAL_TEST(C==A){}
C = I+A;
CGAL_TEST(LA::rank(C)==mat_dim){}
/* matrix operations 2* , |determinant| */
C = NT(2) * I;
C = C * NT(2);
Matrix L,U;
Vector c;
std::vector<int> q;
NT det = LA::determinant(C, L, U, q, c); // det must be 2^{2*mat_dim}
NT pot = 1; for (i=0; i<mat_dim; ++i) pot *= NT(4);
CGAL_TEST(det == pot){}
CGAL_TEST(LA::verify_determinant(C, det, L, U, q, c)){}
CGAL_TEST(det == LA::determinant(C)){}
CGAL_TEST(CGAL_NTS sign(det) == LA::sign_of_determinant(C)){}
if (IOTEST) CGAL_IO_TEST(A,C,CGAL::IO::ASCII);
// add binary test later: if (IOTEST) CGAL_IO_TEST(A,C,CGAL::IO::BINARY);
/* a random linear solver task: */
Vector b(mat_dim),x(mat_dim),e;
NT denom;
for (i = 0; i < mat_dim; i++) {
for (j = 0; j < mat_dim; j++)
C(i,j) = CGAL::get_default_random().get_int(-mat_dim,mat_dim);
b[i] = CGAL::get_default_random().get_int(-mat_dim,mat_dim);
}
for (double f = 1.0; f > 0.0; f-=0.1) {
Matrix E = C;
for (i = 0; i< mat_dim; ++i)
for (j = 0; j < mat_dim; ++j)
if ( CGAL::get_default_random().get_double() > f ) E(i,j)=0;
if (LA::linear_solver(E, b, x, denom, A, e)) {
CGAL_TEST(E*x == b*denom){}
LA::linear_solver(E,b,x,denom,e);
CGAL_TEST(E*x == b*denom){}
LA::linear_solver(E,b,x,denom);
CGAL_TEST(E*x == b*denom){}
CGAL_TEST(LA::is_solvable(E,b)){}
}
else {
Vector null(mat_dim);
CGAL_TEST(LA::transpose(E)*e==null && b*e!=NT(0)){}
CGAL_TEST(!LA::is_solvable(E,b)){}
}
Matrix SV;
if (LA::homogeneous_linear_solver(E,x)) {
CGAL_TEST(E*x==Vector(mat_dim)){}
int r = LA::homogeneous_linear_solver(E,SV);
CGAL_TEST(r==LA::rank(SV)){}
}
// inverse
if (LA::inverse(E,D,denom,c)) {
CGAL_TEST(E*D==denom*Matrix(mat_dim,Matrix::Identity())){}
CGAL_TEST(D==LA::inverse(E,denom)){}
} else {
CGAL_TEST(LA::transpose(E)*c==Vector(mat_dim)){}
}
LA::independent_columns(E,q);
}
if (mat_dim > 1) { // ueberbestimmt:
Matrix N(mat_dim,mat_dim-1);
for (i = 0; i < mat_dim; i++) {
for (j = 0; j < mat_dim-1; j++)
N(i,j) = CGAL::get_default_random().get_int(-mat_dim,mat_dim);
b[i] = CGAL::get_default_random().get_int(-mat_dim,mat_dim);
}
if (LA::linear_solver(N,b,x,denom,e)) {
CGAL_TEST(N*x == b*denom){}
} else {
Vector null(mat_dim-1);
CGAL_TEST(LA::transpose(N)*e==null && b*e!=NT(0)){}
}
}
}
}
CGAL_TEST_END;
}
Reference in New Issue View Git Blame Copy Permalink