Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions src/core/CrossCorrelation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -106,11 +106,11 @@ void CrossCorrelation::readCCCBin() {
int K = this->order + 1;
this->Left = MatrixXd::Zero(K * K, 2 * K);
this->Right = MatrixXd::Zero(K * K, 2 * K);
double dL[2 * K];
double dR[2 * K];
std::vector<double> dL(2 * K);
std::vector<double> dR(2 * K);
for (int i = 0; i < K * K; i++) {
L_fis.read((char *)dL, sizeof(double) * 2 * K);
R_fis.read((char *)dR, sizeof(double) * 2 * K);
L_fis.read(reinterpret_cast<char *>(dL.data()), sizeof(double) * 2 * K);
R_fis.read(reinterpret_cast<char *>(dR.data()), sizeof(double) * 2 * K);
for (int j = 0; j < 2 * K; j++) {
if (std::abs(dL[j]) < MachinePrec) dL[j] = 0.0;
if (std::abs(dR[j]) < MachinePrec) dR[j] = 0.0;
Expand Down
8 changes: 4 additions & 4 deletions src/core/MWFilter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -206,14 +206,14 @@ void MWFilter::generateBlocks() {

int K = this->order + 1;

double dH[K];
double dG[K];
std::vector<double> dH(K);
std::vector<double> dG(K);
/* read H0 and G0 from disk */
this->G0 = Eigen::MatrixXd::Zero(K, K);
this->H0 = Eigen::MatrixXd::Zero(K, K);
for (int i = 0; i < K; i++) {
H_fis.read((char *)dH, sizeof(double) * K);
G_fis.read((char *)dG, sizeof(double) * K);
H_fis.read(reinterpret_cast<char *>(dH.data()), sizeof(double) * K);
G_fis.read(reinterpret_cast<char *>(dG.data()), sizeof(double) * K);
for (int j = 0; j < K; j++) {
this->G0(i, j) = dG[j]; // G0
this->H0(i, j) = dH[j]; // H0
Expand Down
1 change: 0 additions & 1 deletion src/functions/Gaussian.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -194,7 +194,6 @@ template <int D> GaussExp<D> Gaussian<D>::periodify(const std::array<double, D>
auto needed_cells_vec = std::vector<int>();
for (auto i = 0; i < D; i++) {
auto upper_bound = pos[i] + x_std;
auto lower_bound = pos[i] - x_std;
// number of cells upp and down relative to the center of the Gaussian
needed_cells_vec.push_back(std::ceil(upper_bound / period[i]));
}
Expand Down
5 changes: 2 additions & 3 deletions src/treebuilders/ABGVCalculator.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ ABGVCalculator::ABGVCalculator(const ScalingBasis &basis, double a, double b)

void ABGVCalculator::calcValueVectors(const ScalingBasis &basis) {
int kp1 = basis.getQuadratureOrder();
double sqrtCoef[kp1];
std::vector<double> sqrtCoef(kp1);
for (int i = 0; i < kp1; i++) { sqrtCoef[i] = std::sqrt(2.0 * i + 1.0); }

switch (basis.getScalingType()) {
Expand All @@ -74,7 +74,7 @@ void ABGVCalculator::calcValueVectors(const ScalingBasis &basis) {

void ABGVCalculator::calcKMatrix(const ScalingBasis &basis) {
int kp1 = basis.getQuadratureOrder();
double sqrtCoef[kp1];
std::vector<double> sqrtCoef(kp1);
for (int i = 0; i < kp1; i++) { sqrtCoef[i] = std::sqrt(2.0 * i + 1.0); }
getQuadratureCache(qCache);
const VectorXd &roots = qCache.getRoots(kp1);
Expand Down Expand Up @@ -105,7 +105,6 @@ void ABGVCalculator::calcNode(MWNode<2> &node) {
node.zeroCoefs();

const auto &idx = node.getNodeIndex();
int l = idx[1] - idx[0];
int np1 = idx.getScale() + 1;
int kp1 = node.getKp1();
int kp1_d = node.getKp1_d();
Expand Down
5 changes: 2 additions & 3 deletions src/treebuilders/ConvolutionCalculator.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,6 @@ template <int D, typename T> MWNodeVector<D, T> *ConvolutionCalculator<D, T>::ma
auto *band = new MWNodeVector<D, T>;

int o_depth = gNode.getScale() - this->oper->getOperatorRoot();
int g_depth = gNode.getDepth();
int width = this->oper->getMaxBandWidth(o_depth);

bool periodic = gNode.getMWTree().isPeriodic();
Expand Down Expand Up @@ -228,8 +227,8 @@ template <int D, typename T> void ConvolutionCalculator<D, T>::calcNode(MWNode<D

int o_depth = gNode.getScale() - this->oper->getOperatorRoot();
if (manipulateOperator and this->oper->getOperatorRoot() < 0) o_depth = gNode.getDepth();
T tmpCoefs[gNode.getNCoefs()];
OperatorState<D, T> os(gNode, tmpCoefs);
std::vector<T> tmpCoefs(gNode.getNCoefs());
OperatorState<D, T> os(gNode, tmpCoefs.data());
this->operStat.incrementGNodeCounters(gNode);

// Get all nodes in f within the bandwith of O in g
Expand Down
14 changes: 4 additions & 10 deletions src/treebuilders/DerivativeCalculator.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -90,8 +90,8 @@ template <int D, typename T> void DerivativeCalculator<D, T>::calcNode(MWNode<D,
// if (this->oper->getMaxBandWidth() > 1) MSG_ABORT("Only implemented for zero bw");
outNode.zeroCoefs();
int nComp = (1 << D);
T tmpCoefs[outNode.getNCoefs()];
OperatorState<D, T> os(outNode, tmpCoefs);
std::vector<T> tmpCoefs(outNode.getNCoefs());
OperatorState<D, T> os(outNode, tmpCoefs.data());

os.setFNode(inpNode);
os.setFIndex(inpNode.nodeIndex);
Expand All @@ -116,8 +116,8 @@ template <int D, typename T> void DerivativeCalculator<D, T>::calcNode(MWNode<D,
gNode.zeroCoefs();

int nComp = (1 << D);
T tmpCoefs[gNode.getNCoefs()];
OperatorState<D, T> os(gNode, tmpCoefs);
std::vector<T> tmpCoefs(gNode.getNCoefs());
OperatorState<D, T> os(gNode, tmpCoefs.data());
this->operStat.incrementGNodeCounters(gNode);

// Get all nodes in f within the bandwith of O in g
Expand Down Expand Up @@ -182,11 +182,8 @@ template <int D, typename T> void DerivativeCalculator<D, T>::applyOperator_bw0(
// cout<<" applyOperator "<<endl;
MWNode<D, T> &gNode = *os.gNode;
MWNode<D, T> &fNode = *os.fNode;
const NodeIndex<D> &fIdx = *os.fIdx;
const NodeIndex<D> &gIdx = gNode.getNodeIndex();
int depth = gNode.getDepth();

double oNorm = 1.0;
double **oData = os.getOperData();

for (int d = 0; d < D; d++) {
Expand Down Expand Up @@ -218,7 +215,6 @@ template <int D, typename T> void DerivativeCalculator<D, T>::applyOperator(Oper
const NodeIndex<D> &gIdx = gNode.getNodeIndex();
int depth = gNode.getDepth();

double oNorm = 1.0;
double **oData = os.getOperData();

for (int d = 0; d < D; d++) {
Expand All @@ -236,8 +232,6 @@ template <int D, typename T> void DerivativeCalculator<D, T>::applyOperator(Oper

const OperatorNode &oNode = oTree.getNode(depth, oTransl);
int oIdx = os.getOperIndex(d);
double ocn = oNode.getComponentNorm(oIdx);
oNorm *= ocn;
if (this->applyDir == d) {
oData[d] = const_cast<double *>(oNode.getCoefs()) + oIdx * os.kp1_2;
} else {
Expand Down
1 change: 0 additions & 1 deletion src/treebuilders/grid.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,6 @@ template <int D> void build_grid(FunctionTree<D> &out, const GaussExp<D> &inp, i
builder.build(out, calculator, adaptor, maxIter);
}
} else {
auto period = out.getMRA().getWorldBox().getScalingFactors();
for (auto i = 0; i < inp.size(); i++) {
auto *gauss = inp.getFunc(i).copy();
build_grid(out, *gauss, maxIter);
Expand Down
8 changes: 4 additions & 4 deletions src/treebuilders/multiply.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -334,8 +334,8 @@ template <int D, typename T> double node_norm_dot(FunctionTree<D, T> &bra, Funct

double result = 0.0;
int ncoef = bra.getKp1_d() * bra.getTDim();
T valA[ncoef];
T valB[ncoef];
std::vector<T> valA(ncoef);
std::vector<T> valB(ncoef);
int nNodes = bra.getNEndNodes();

for (int n = 0; n < nNodes; n++) {
Expand All @@ -345,8 +345,8 @@ template <int D, typename T> double node_norm_dot(FunctionTree<D, T> &bra, Funct
// convert to interpolating coef, take abs, convert back
FunctionNode<D, T> *mwNode = static_cast<FunctionNode<D, T> *>(ket.findNode(idx));
if (mwNode == nullptr) MSG_ABORT("Trees must have same grid");
node.getAbsCoefs(valA);
mwNode->getAbsCoefs(valB);
node.getAbsCoefs(valA.data());
mwNode->getAbsCoefs(valB.data());
for (int i = 0; i < ncoef; i++) result += std::norm(valA[i] * valB[i]);
} else {
// approximate by product of node norms
Expand Down
4 changes: 2 additions & 2 deletions src/trees/FunctionNode.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -130,7 +130,7 @@ template <int D, typename T> T FunctionNode<D, T>::integrateInterpolating() cons
int qOrder = this->getKp1();
getQuadratureCache(qc);
const VectorXd &weights = qc.getWeights(qOrder);
double sqWeights[qOrder];
std::vector<double> sqWeights(qOrder);
for (int i = 0; i < qOrder; i++) sqWeights[i] = std::sqrt(weights[i]);

int kp1_p[D];
Expand Down Expand Up @@ -170,7 +170,7 @@ template <int D, typename T> T FunctionNode<D, T>::integrateValues() const {
this->getCoefs(coefs);
int ncoefs = coefs.size();
int ncoefChild = ncoefs / (1 << D);
T cc[ncoefChild];
std::vector<T> cc(ncoefChild);
// factorize out the children
for (int i = 0; i < ncoefChild; i++) cc[i] = coefs[i];
for (int j = 1; j < (1 << D); j++)
Expand Down
24 changes: 4 additions & 20 deletions src/trees/FunctionTree.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -290,9 +290,6 @@ template <int D, typename T> void FunctionTree<D, T>::loadTreeTXT(const std::str
* @param[in] file: File name
*/
template <int D, typename T> void FunctionTree<D, T>::saveTreeTXT(const std::string &fname) {
int nRoots = this->getRootBox().size();
MWNode<D, T> **roots = this->getRootBox().getNodes();

std::ofstream out(fname);
out << std::setprecision(14);
out << D << std::endl;
Expand All @@ -306,7 +303,7 @@ template <int D, typename T> void FunctionTree<D, T>::saveTreeTXT(const std::str
int ncoefs = 1;
for (int d = 0; d < D; d++) ncoefs *= kp1;
int Tdim = std::pow(2, D);
T values[ncoefs * Tdim];
std::vector<T> values(ncoefs * Tdim);

int nout = this->endNodeTable.size();
out << Tdim * nout << std::endl; // could output only scaling coeff?
Expand Down Expand Up @@ -334,7 +331,7 @@ template <int D, typename T> void FunctionTree<D, T>::saveTreeTXT(const std::str
MWNode<D, T> *node = &(this->getNode(idx, false));
T *coefs = node->getCoefs();
for (int i = 0; i < ncoefs * Tdim; i++) values[i] = coefs[i];
node->attachCoefs(values);
node->attachCoefs(values.data());
int n = idx.getScale();
node->mwTransform(Reconstruction);
node->cvTransform(Forward);
Expand Down Expand Up @@ -443,7 +440,6 @@ template <> double FunctionTree<3, double>::integrateEndNodes(RepresentableFunct
// traverse tree, and treat end nodes only
std::vector<FunctionNode<3> *> stack; // node from this
for (int i = 0; i < this->getRootBox().size(); i++) stack.push_back(&(this->getRootFuncNode(i)));
int basis = getMRA().getScalingBasis().getScalingType();
double result = 0.0;
while (stack.size() > 0) {
FunctionNode<3> *Node = stack.back();
Expand All @@ -456,9 +452,9 @@ template <> double FunctionTree<3, double>::integrateEndNodes(RepresentableFunct
double *coefs = Node->getCoefs(); // save position of coeff, but do not use them!
// The data in fmat is not organized so that two consecutive points are stored after each other in memory, so needs to copy before mwtransform, cannot use memory adress directly.
int nc = fmat.cols();
double cc[nc];
std::vector<double> cc(nc);
for (int i = 0; i < nc; i++) cc[i] = fmat(0, i);
Node->attachCoefs(cc);
Node->attachCoefs(cc.data());
result += Node->integrateValues();
Node->attachCoefs(coefs); // put back original coeff
}
Expand Down Expand Up @@ -880,8 +876,6 @@ void FunctionTree<D, T>::makeCoeffVector(std::vector<T *> &coefs,
indices.clear();
parent_indices.clear();
max_index = 0;
int sizecoeff = (1 << refTree.getDim()) * refTree.getKp1_d();
int sizecoeffW = ((1 << refTree.getDim()) - 1) * refTree.getKp1_d();
std::vector<MWNode<D, double> *> refstack; // nodes from refTree
std::vector<MWNode<D, T> *> thisstack; // nodes from this Tree
for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) {
Expand Down Expand Up @@ -1096,7 +1090,6 @@ template <> int FunctionTree<3, double>::saveNodesAndRmCoeff() {
for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { stack.push_back(this->getRootBox().getNodes()[rIdx]); }
while (stack.size() > stack_p) {
MWNode<3, double> *Node = stack[stack_p++];
int id = 0;
NodesCoeff->put_data(Node->getNodeIndex(), sizecoeff, Node->getCoefs());
for (int i = 0; i < Node->getNChildren(); i++) { stack.push_back(Node->children[i]); }
}
Expand All @@ -1119,7 +1112,6 @@ template <> int FunctionTree<3, ComplexDouble>::saveNodesAndRmCoeff() {
for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { stack.push_back(this->getRootBox().getNodes()[rIdx]); }
while (stack.size() > stack_p) {
MWNode<3, ComplexDouble> *Node = stack[stack_p++];
int id = 0;
NodesCoeff->put_data(Node->getNodeIndex(), sizecoeff, Node->getCoefs());
for (int i = 0; i < Node->getNChildren(); i++) { stack.push_back(Node->children[i]); }
}
Expand Down Expand Up @@ -1195,7 +1187,6 @@ template <int D, typename T> FunctionTree<D, double> *FunctionTree<D, T>::Imag()

template <> void FunctionTree<3, double>::CopyTreeToComplex(FunctionTree<3, ComplexDouble> *&outTree) {
delete outTree;
double ref = 0.0;
outTree = new FunctionTree<3, ComplexDouble>(this->getMRA());
std::vector<MWNode<3, double> *> instack; // node from this
std::vector<MWNode<3, ComplexDouble> *> outstack; // node from outTree
Expand All @@ -1204,7 +1195,6 @@ template <> void FunctionTree<3, double>::CopyTreeToComplex(FunctionTree<3, Comp
instack.push_back(this->getRootBox().getNodes()[rIdx]);
outstack.push_back(outTree->getRootBox().getNodes()[rIdx]);
}
int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk());
int ncoefs = this->getNodeAllocator().getNCoefs();
while (instack.size() > 0) {
// inNode and outNode are the same node in space, but on different trees
Expand Down Expand Up @@ -1235,7 +1225,6 @@ template <> void FunctionTree<3, double>::CopyTreeToComplex(FunctionTree<3, Comp

template <> void FunctionTree<2, double>::CopyTreeToComplex(FunctionTree<2, ComplexDouble> *&outTree) {
delete outTree;
double ref = 0.0;
outTree = new FunctionTree<2, ComplexDouble>(this->getMRA());
std::vector<MWNode<2, double> *> instack; // node from this
std::vector<MWNode<2, ComplexDouble> *> outstack; // node from outTree
Expand All @@ -1244,7 +1233,6 @@ template <> void FunctionTree<2, double>::CopyTreeToComplex(FunctionTree<2, Comp
instack.push_back(this->getRootBox().getNodes()[rIdx]);
outstack.push_back(outTree->getRootBox().getNodes()[rIdx]);
}
int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk());
int ncoefs = this->getNodeAllocator().getNCoefs();
while (instack.size() > 0) {
// inNode and outNode are the same node in space, but on different trees
Expand Down Expand Up @@ -1275,7 +1263,6 @@ template <> void FunctionTree<2, double>::CopyTreeToComplex(FunctionTree<2, Comp

template <> void FunctionTree<1, double>::CopyTreeToComplex(FunctionTree<1, ComplexDouble> *&outTree) {
delete outTree;
double ref = 0.0;
outTree = new FunctionTree<1, ComplexDouble>(this->getMRA());
std::vector<MWNode<1, double> *> instack; // node from this
std::vector<MWNode<1, ComplexDouble> *> outstack; // node from outTree
Expand All @@ -1284,7 +1271,6 @@ template <> void FunctionTree<1, double>::CopyTreeToComplex(FunctionTree<1, Comp
instack.push_back(this->getRootBox().getNodes()[rIdx]);
outstack.push_back(outTree->getRootBox().getNodes()[rIdx]);
}
int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk());
int ncoefs = this->getNodeAllocator().getNCoefs();
while (instack.size() > 0) {
// inNode and outNode are the same node in space, but on different trees
Expand Down Expand Up @@ -1316,7 +1302,6 @@ template <> void FunctionTree<1, double>::CopyTreeToComplex(FunctionTree<1, Comp
// for testing
template <> void FunctionTree<3, double>::CopyTreeToReal(FunctionTree<3, double> *&outTree) {
delete outTree;
double ref = 0.0;
// FunctionTree<3, double>* inTree = this;
outTree = new FunctionTree<3, double>(this->getMRA());
std::vector<MWNode<3, double> *> instack; // node from this
Expand All @@ -1326,7 +1311,6 @@ template <> void FunctionTree<3, double>::CopyTreeToReal(FunctionTree<3, double>
instack.push_back(this->getRootBox().getNodes()[rIdx]);
outstack.push_back(outTree->getRootBox().getNodes()[rIdx]);
}
int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk());
int ncoefs = this->getNodeAllocator().getNCoefs();
while (instack.size() > 0) {
// inNode and outNode are the same node in space, but on different trees
Expand Down
15 changes: 15 additions & 0 deletions src/trees/HilbertPath.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -114,6 +114,21 @@ const int HilbertPath<3>::hTable[12][8] = {
{4,5,7,6,3,2,0,1}
};

template <int D>
short int HilbertPath<D>::getChildPath(int hIdx) const {
return pTable[this->path][hIdx];
}

template <int D>
int HilbertPath<D>::getZIndex(int hIdx) const {
return zTable[this->path][hIdx];
}

template <int D>
int HilbertPath<D>::getHIndex(int zIdx) const {
return hTable[this->path][zIdx];
}

// clang-format on

template class HilbertPath<1>;
Expand Down
8 changes: 4 additions & 4 deletions src/trees/HilbertPath.h
Original file line number Diff line number Diff line change
Expand Up @@ -42,10 +42,10 @@ template <int D> class HilbertPath final {
}

short int getPath() const { return this->path; }
short int getChildPath(int hIdx) const { return this->pTable[this->path][hIdx]; }
short int getChildPath(int hIdx) const;

int getZIndex(int hIdx) const { return this->zTable[this->path][hIdx]; }
int getHIndex(int zIdx) const { return this->hTable[this->path][zIdx]; }
int getZIndex(int hIdx) const;
int getHIndex(int zIdx) const;

private:
short int path{0};
Expand All @@ -54,4 +54,4 @@ template <int D> class HilbertPath final {
static const int hTable[][8];
};

} // namespace mrcpp
} // namespace mrcpp
Loading
Loading