Commit fe6c4a87 authored by Benjamin Huber's avatar Benjamin Huber

Merge branch 'v3' of hemioGit:xerus/xerus into v3

parents fcc8e1d4 09d0ae8c
Pipeline #629 passed with stages
in 8 minutes and 4 seconds
......@@ -75,7 +75,7 @@ static misc::UnitTest tensor_solve_sparse("Tensor", "solve_sparse", [](){
std::mt19937_64 &rnd = xerus::misc::randomEngine;
std::normal_distribution<double> dist(0.0, 1.0);
const size_t N = 100;
std::uniform_int_distribution<size_t> eDist(1,N*N-1);
std::uniform_int_distribution<size_t> eDist(1, N*N-1);
Index i,j,k;
......@@ -100,19 +100,18 @@ static misc::UnitTest tensor_solve_sparse("Tensor", "solve_sparse", [](){
internal::CholmodSparse idt(id.get_unsanitized_sparse_data(), N, N, false);
Tensor id2({N,N});
id2.get_unsanitized_sparse_data() = idt.to_map();
MTEST(frob_norm(id-id2) < 1e-15, frob_norm(id-id2));
MTEST(frob_norm(id-id2) < 1e-12, frob_norm(id-id2));
Tensor fid(id);
fid.use_dense_representation();
Tensor fx;
// id.use_sparse_representation();
TEST(id.is_sparse());
fx(i) = r(j) / fid(j,i);
x(i) = r(j) / id(j,i);
MTEST(frob_norm(id(j,i)*x(i) - r(j))/frob_norm(x)<1e-13, frob_norm(id(j,i)*x(i) - r(j))/frob_norm(x));
MTEST(frob_norm(fid(j,i)*fx(i) - r(j))/frob_norm(x)<1e-13, frob_norm(fid(j,i)*fx(i) - r(j))/frob_norm(x));
MTEST(frob_norm(fx-x)/frob_norm(x)<1e-13, frob_norm(fx-x)/frob_norm(x));
MTEST(frob_norm(id(j,i)*x(i) - r(j))/frob_norm(x)<1e-12, frob_norm(id(j,i)*x(i) - r(j))/frob_norm(x));
MTEST(frob_norm(fid(j,i)*fx(i) - r(j))/frob_norm(x)<1e-12, frob_norm(fid(j,i)*fx(i) - r(j))/frob_norm(x));
MTEST(frob_norm(fx-x)/frob_norm(x)<1e-12, frob_norm(fx-x)/frob_norm(x));
});
static misc::UnitTest tensor_solve_trans("Tensor", "solve_transposed", [](){
......@@ -136,7 +135,7 @@ static misc::UnitTest tensor_solve_trans("Tensor", "solve_transposed", [](){
Tensor x1, x2;
x1(i) = r(j) / A(i,j);
x2(i) = r(j) / At(j,i);
MTEST(frob_norm(x1-x2) < 1e-14, "s " << frob_norm(x1-x2));
MTEST(frob_norm(x1-x2) < 1e-12, "s " << frob_norm(x1-x2));
A.use_dense_representation();
At.use_dense_representation();
......@@ -146,10 +145,10 @@ static misc::UnitTest tensor_solve_trans("Tensor", "solve_transposed", [](){
x4(i) = r(j) / At(j,i);
residual(i) = A(i,j) * (x3(j) - x4(j));
MTEST(frob_norm(x3-x4) < 1e-14, "d " << frob_norm(x3-x4) << " residual: " << frob_norm(residual));
MTEST(frob_norm(x3-x4) < 1e-12, "d " << frob_norm(x3-x4) << " residual: " << frob_norm(residual));
residual(i) = A(i,j) * (x1(j) - x3(j));
MTEST(frob_norm(x1-x3)/frob_norm(x1) < 5e-14, "sd " << frob_norm(x1-x3)/frob_norm(x1) << " residual: " << frob_norm(residual));
MTEST(frob_norm(x1-x3)/frob_norm(x1) < 1e-12, "sd " << frob_norm(x1-x3)/frob_norm(x1) << " residual: " << frob_norm(residual));
});
......@@ -172,7 +171,6 @@ static misc::UnitTest tensor_solve_matrix("Tensor", "solve_matrix", [](){
for(size_t xi = 0; xi < degN; ++xi) { nDims.push_back(n2Dist(rnd)); }
for(size_t xi = 0; xi < degP; ++xi) { pDims.push_back(n2Dist(rnd)); }
XERUS_LOG(bla, "Set " << mDims << "x" << nDims << " * " << pDims);
auto A = Tensor::random(mDims | nDims);
A *= realDist(rnd);
......@@ -194,7 +192,6 @@ static misc::UnitTest tensor_solve_matrix("Tensor", "solve_matrix", [](){
residual(i^degM, k^degP) = A(i^degM, j^degN)*X(j^degN, k^degP) - B(i^degM, k^degP);
MTEST(frob_norm(residual) < 1e-10, frob_norm(residual));
LOG(asd, degN << ' '<< degP);
X(j^degN, k^degP) = B(i^degM, k^degP) / A(i^degM, j^degN); //solve_least_squares(X, A, B, degP);
......
......@@ -121,7 +121,7 @@ static misc::UnitTest alg_adf_rnd("Algorithm", "adf_random_low_rank", [](){
}
TEST(test);
ADFVariant ourADF(500, 1e-6, 1e-6);
ADFVariant ourADF(500, 1e-6, 0.999);
TTTensor X = TTTensor::ones(std::vector<size_t>(D, N));
PerformanceData perfData([&](const TTTensor& _x) {return frob_norm(_x - trueSolution)/frob_norm(trueSolution);}, true, false);
......
......@@ -565,7 +565,7 @@ namespace xerus {
// If we follow a rank increasing strategie, increase the ransk until we reach the targetResidual, the maxRanks or the maxIterations.
while(residualNorm > targetResidualNorm && x.ranks() != maxRanks && (maxIterations == 0 || iteration < maxIterations)) {
LOG(xRanKResBefore, measurments.test(x));
// LOG(xRanKResBefore, measurments.test(x));
// Increase the ranks
x.move_core(0, true);
const auto rndTensor = TTTensor::random(x.dimensions, std::vector<size_t>(x.degree()-1, 1));
......
......@@ -540,7 +540,7 @@ namespace xerus {
std::unique_ptr<int[]> pivot(new int[_n]);
misc::set_zero(pivot.get(), _n);
// std::unique_ptr<double[]> signulars(new double[std::min(_n, _m)]);
std::unique_ptr<double[]> signulars(new double[std::min(_n, _m)]);
int rank;
......@@ -553,20 +553,7 @@ namespace xerus {
misc::set_zero(bOrX+_m*_p, (_n-_m)*_p); // Lapacke is unhappy if the array contains NANs...
}
IF_CHECK( int lapackAnswer = ) LAPACKE_dgelsy(
LAPACK_ROW_MAJOR,
static_cast<int>(_m), // Left dimension of A
static_cast<int>(_n), // Right dimension of A
static_cast<int>(_p), // Number of rhss
_A, // Matrix A
static_cast<int>(_n), // LDA
bOrX, // On input b, on output x
static_cast<int>(_p), // LDB
pivot.get(), // Pivot, entries must be zero to allow pivoting
xerus::EPSILON, // Used to determine the accuracy of the Lapacke call. Basically all singular values smaller than RCOND*s[0] are ignored. (s[0] is the largest signular value)
&rank); // Outputs the rank of A
// IF_CHECK( int lapackAnswer = ) LAPACKE_dgelsd(
// IF_CHECK( int lapackAnswer = ) LAPACKE_dgelsy(
// LAPACK_ROW_MAJOR,
// static_cast<int>(_m), // Left dimension of A
// static_cast<int>(_n), // Right dimension of A
......@@ -574,11 +561,24 @@ namespace xerus {
// _A, // Matrix A
// static_cast<int>(_n), // LDA
// bOrX, // On input b, on output x
// static_cast<int>(_p), // LDB
// signulars.get(), // Pivot, entries must be zero to allow pivoting
// static_cast<int>(_p), // LDB
// pivot.get(), // Pivot, entries must be zero to allow pivoting
// xerus::EPSILON, // Used to determine the accuracy of the Lapacke call. Basically all singular values smaller than RCOND*s[0] are ignored. (s[0] is the largest signular value)
// &rank); // Outputs the rank of A
IF_CHECK( int lapackAnswer = ) LAPACKE_dgelsd(
LAPACK_ROW_MAJOR,
static_cast<int>(_m), // Left dimension of A
static_cast<int>(_n), // Right dimension of A
static_cast<int>(_p), // Number of rhss
_A, // Matrix A
static_cast<int>(_n), // LDA
bOrX, // On input b, on output x
static_cast<int>(_p), // LDB
signulars.get(), // Pivot, entries must be zero to allow pivoting
xerus::EPSILON, // Used to determine the accuracy of the Lapacke call. Basically all singular values smaller than RCOND*s[0] are ignored. (s[0] is the largest signular value)
&rank); // Outputs the rank of A
CHECK(lapackAnswer == 0, error, "Unable to solves min ||Ax - b||_2 for x. Lapacke says: " << lapackAnswer << " sizes are " << _m << " x " << _n << " * " << _p);
if(_m >= _n) { // I.e. bOrX is _b
......
......@@ -251,7 +251,6 @@ namespace xerus {
_out.assign_indices();
// Assign the dimensions
_base.assign_index_dimensions();
_out.assign_index_dimensions();
......
......@@ -22,156 +22,59 @@
* @brief Implementation of the indexed tensor / operator.
*/
#include <xerus/misc/internal.h>
#include <xerus/misc/containerSupport.h>
#include <xerus/indexedTensor_tensor_factorisations.h>
#include <xerus/index.h>
#include <xerus/indexedTensor.h>
#include <xerus/indexedTensorMoveable.h>
#include <xerus/tensor.h>
#include <xerus/blasLapackWrapper.h>
#include <xerus/cholmod_wrapper.h>
namespace xerus {
void solve(internal::IndexedTensorWritable<Tensor>&& _x, internal::IndexedTensorReadOnly<Tensor>&& _a, internal::IndexedTensorReadOnly<Tensor>&& _b) {
// x takes the dimensions of A -- also ensures that every index of x is contained in A
_x.tensorObject->reset(_a.get_evaluated_dimensions(_x.indices), Tensor::Initialisation::None);
_a.assign_indices();
internal::IndexedTensorMoveable<Tensor> operator/ (internal::IndexedTensorReadOnly<Tensor>&& _b, internal::IndexedTensorReadOnly<Tensor>&& _A) {
_A.assign_indices();
_b.assign_indices();
// IF_CHECK( _x.check_indices(false); )
// for(size_t i = 0; i < bAssIndices.numIndices; ++i) {
// REQUIRE(!bAssIndices.indexOpen[i] || contains(AAssIndices.indices, bAssIndices.indices[i]), "Every open index of b must be contained in A.");
// REQUIRE(!contains(_x.indices, bAssIndices.indices[i]), "x and b must not have indices in common.");
// }
// If possible we don't want to reorder A
size_t extraDims = 0;
std::vector<Index> orderA;
std::vector<Index> orderB;
std::vector<Index> orderX;
std::vector<size_t> dimensionsA;
std::vector<size_t> dimensionsB;
std::vector<size_t> dimensionsX;
size_t dimensionsCount = 0;
for(const Index& idx : _a.indices) {
// If possible we don't want to reorder A, so first divide A into those shared with b and x.
for(const Index& idx : _A.indices) {
if(misc::contains(_b.indices, idx)) {
orderA.push_back(idx);
orderB.push_back(idx);
for(size_t i = 0; i < idx.span; ++i) {
dimensionsA.push_back(_a.tensorObjectReadOnly->dimensions[dimensionsCount]);
dimensionsB.push_back(_a.tensorObjectReadOnly->dimensions[dimensionsCount]);
dimensionsCount++;
}
} else {
orderX.push_back(idx);
for(size_t i = 0; i < idx.span; ++i) {
dimensionsX.push_back(_a.tensorObjectReadOnly->dimensions[dimensionsCount++]);
}
}
}
orderA.insert(orderA.end(), orderX.begin(), orderX.end());
dimensionsA.insert(dimensionsA.end(), dimensionsX.begin(), dimensionsX.end());
//Save slot for eventual tmpX and tmpB
std::unique_ptr<internal::IndexedTensor<Tensor>> saveSlotX;
internal::IndexedTensorWritable<Tensor>* usedX;
std::unique_ptr<internal::IndexedTensor<Tensor>> saveSlotB;
internal::IndexedTensorReadOnly<Tensor>* usedB;
// We need tmp objects for A and b, because Lapacke wants to destroys the input
// and we need to cast b if it is sparse as we cannot handle that yet
if (!_a.tensorObjectReadOnly->is_sparse()) {
saveSlotB.reset(new internal::IndexedTensor<Tensor>(
new Tensor(std::move(dimensionsB), Tensor::Representation::Dense, Tensor::Initialisation::None), orderB, true
));
evaluate(std::move(*saveSlotB), std::move(_b));
saveSlotB->tensorObject->use_dense_representation();
usedB = saveSlotB.get();
} else {
usedB = &_b;
}
bool sparseResult = (_a.tensorObjectReadOnly->is_sparse() && usedB->tensorObjectReadOnly->is_sparse());
if (orderX != _x.indices || (_x.tensorObject->is_sparse() != sparseResult)) {
if (sparseResult) {
saveSlotX.reset(new internal::IndexedTensor<Tensor>(
new Tensor(std::move(dimensionsX), Tensor::Representation::Sparse, Tensor::Initialisation::None), orderX, true));
} else {
saveSlotX.reset(new internal::IndexedTensor<Tensor>(
new Tensor(std::move(dimensionsX), Tensor::Representation::Dense, Tensor::Initialisation::None), orderX, true));
}
usedX = saveSlotX.get();
} else {
usedX = &_x;
}
// Assume A is an MxN matrix
const size_t M = usedB->tensorObjectReadOnly->size;
const size_t N = usedX->tensorObjectReadOnly->size;
// This far the dimensions of A and B coincide
orderB = orderA;
internal::IndexedTensor<Tensor> tmpA(new Tensor(std::move(dimensionsA), _a.tensorObjectReadOnly->representation, Tensor::Initialisation::None), orderA, true);
evaluate(std::move(tmpA), std::move(_a));
// Add the second part of indices in the order obtained for X
orderA.insert(orderA.end(), orderX.begin(), orderX.end());
if (tmpA.tensorObjectReadOnly->is_sparse()) {
if (usedB->tensorObjectReadOnly->is_sparse()) {
internal::CholmodSparse::solve_sparse_rhs(
usedX->tensorObject->get_unsanitized_sparse_data(), N,
tmpA.tensorObjectReadOnly->get_unsanitized_sparse_data(), false,
usedB->tensorObjectReadOnly->get_unsanitized_sparse_data(), M);
} else {
internal::CholmodSparse::solve_dense_rhs(
usedX->tensorObject->get_unsanitized_dense_data(), N,
tmpA.tensorObjectReadOnly->get_unsanitized_sparse_data(), false,
usedB->tensorObjectReadOnly->get_unsanitized_dense_data(), M);
// Now complete indices of b and x with those not shared with A ( in order of b as we don't want to reorder b if possible).
for(const Index& idx : _b.indices) {
if(!misc::contains(_A.indices, idx)) {
orderB.push_back(idx);
orderX.push_back(idx);
for(size_t i = 0; i < idx.span; ++i) {
extraDims++;
}
}
// Propagate the constant factor
usedX->tensorObject->factor = usedB->tensorObjectReadOnly->factor / tmpA.tensorObjectReadOnly->factor;
} else {
// dense A
// We need tmp objects for A and b, because Lapacke wants to destroys the input
tmpA.tensorObject->ensure_own_data();
blasWrapper::solve_least_squares_destructive(
usedX->tensorObject->override_dense_data(),
tmpA.tensorObject->get_unsanitized_dense_data(), M, N,
saveSlotB->tensorObject->get_unsanitized_dense_data(), 1);
// Propagate the constant factor
usedX->tensorObject->factor = usedB->tensorObjectReadOnly->factor / tmpA.tensorObjectReadOnly->factor;
}
if(saveSlotX) { evaluate(std::move(_x), std::move(*usedX)); }
}
internal::IndexedTensorMoveable<Tensor> operator/ (internal::IndexedTensorReadOnly<Tensor>&& _b, internal::IndexedTensorReadOnly<Tensor>&& _A) {
_A.assign_indices();
_b.assign_indices();
// If indices coincide no reordering occours (only shared data pointer is set).
Tensor reorderedA, reorderedB;
reorderedA(orderA) = std::move(_A);
reorderedB(orderB) = std::move(_b);
std::vector<Index> indicesX;
std::vector<size_t> dimensionsX;
internal::IndexedTensorMoveable<Tensor> tmpX(new Tensor(), std::move(orderX));
size_t dimensionsCount = 0;
for(const Index& idx : _A.indices) {
if(!misc::contains(_b.indices, idx)) {
indicesX.push_back(idx);
for(size_t i = 0; i < idx.span; ++i) {
dimensionsX.push_back(_A.tensorObjectReadOnly->dimensions[dimensionsCount++]);
}
} else {
dimensionsCount += idx.span;
}
}
internal::IndexedTensorMoveable<Tensor> tmpX(new Tensor(std::move(dimensionsX), Tensor::Representation::Dense, Tensor::Initialisation::None), std::move(indicesX));
solve_least_squares(*tmpX.tensorObject, reorderedA, reorderedB, extraDims);
solve(std::move(tmpX), std::move(_A), std::move(_b));
return tmpX;
}
} // namespace xerus
......@@ -1574,31 +1574,36 @@ namespace xerus {
const size_t p = misc::product(_B.dimensions, degM, degM+_extraDegree);
if (_A.is_sparse()) {
LOG(fatal, "Sparse not yet impl");
// if (usedB->tensorObjectReadOnly->is_sparse()) {
// internal::CholmodSparse::solve_sparse_rhs(
// usedX->tensorObject->get_unsanitized_sparse_data(), N,
// tmpA.tensorObjectReadOnly->get_unsanitized_sparse_data(), false,
// usedB->tensorObjectReadOnly->get_unsanitized_sparse_data(), M);
// } else {
// internal::CholmodSparse::solve_dense_rhs(
// usedX->tensorObject->get_unsanitized_dense_data(), N,
// tmpA.tensorObjectReadOnly->get_unsanitized_sparse_data(), false,
// usedB->tensorObjectReadOnly->get_unsanitized_dense_data(), M);
// }
//
// // Propagate the constant factor
// usedX->tensorObject->factor = usedB->tensorObjectReadOnly->factor / tmpA.tensorObjectReadOnly->factor;
REQUIRE( p == 1, "Matrix least squares not supported in sparse");
if (_B.is_sparse()) {
internal::CholmodSparse::solve_sparse_rhs(
_X.override_sparse_data(),
n,
_A.get_unsanitized_sparse_data(),
false,
_B.get_unsanitized_sparse_data(),
m);
} else {
internal::CholmodSparse::solve_dense_rhs(
_X.override_dense_data(),
n,
_A.get_unsanitized_sparse_data(),
false,
_B.get_unsanitized_dense_data(),
m);
}
} else { // Dense A
REQUIRE(_B.is_dense(), "Not yet implemented");
blasWrapper::solve_least_squares(
_X.override_dense_data(),
_A.get_unsanitized_dense_data(), m, n,
_B.get_unsanitized_dense_data(), p);
// Propagate the constant factor
_X.factor = _B.factor / _A.factor;
}
// Propagate the constant factor
_X.factor = _B.factor / _A.factor;
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment