Commit a7c9cb84 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Cleanup and fix minor issues in tests and classes related to interpolation...

Cleanup and fix minor issues in tests and classes related to interpolation based FMM + Update gaussian correlation kernel to include derivatives.
parent 9d132c52
......@@ -440,7 +440,7 @@ public:
FReal CellWidth = RootCellWidth / FReal(2.); // at level 1
CellWidth /= FReal(2.); // at level 2
for (unsigned int l=2; l<TreeHeight; ++l) {
precompute<ORDER>(MatrixKernel, CellWidth, Epsilon, K[l], LowRank[l]);
precompute<FReal,ORDER>(MatrixKernel, CellWidth, Epsilon, K[l], LowRank[l]);
CellWidth /= FReal(2.); // at level l+1
}
}
......
......@@ -50,9 +50,9 @@ enum CORRELATION_SUPPORT_EXTENSION{INFINITE,FINITE};
// hole effect correlation (nonstrictly decreasing)
// Whittle model? (slow decay => important constraint on size of grid / length)
template <class FReal>
struct AbstractCorrelationKernel : FNoCopyable
struct FAbstractCorrelationKernel : FNoCopyable
{
virtual ~AbstractCorrelationKernel(){}
virtual ~FAbstractCorrelationKernel(){}
virtual FReal evaluate(const FReal*, const FReal*) const = 0;
};
......@@ -80,7 +80,7 @@ struct AbstractCorrelationKernel : FNoCopyable
/// Generic Gaussian correlation function
/// Special case of Matern function with $\nu \rightarrow \infty$
template<class FReal>
struct CK_Gauss : AbstractCorrelationKernel<FReal>
struct FInterpMatrixKernelGauss : FAbstractCorrelationKernel<FReal>
{
static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
static const CORRELATION_SUPPORT_EXTENSION Extension = INFINITE;
......@@ -91,24 +91,30 @@ struct CK_Gauss : AbstractCorrelationKernel<FReal>
static const unsigned int NLHS = 1; //< dim of loc exp
FReal lengthScale_;
CK_Gauss(const FReal lengthScale = FReal(1.))
FInterpMatrixKernelGauss(const FReal lengthScale = FReal(1.))
: lengthScale_(lengthScale)
{
std::cout<< "Gaussian " << 3 <<"D"
<< ", i.e. r(x)=exp(-0.5*(x_i/l*x_i/l)) with l="
<< lengthScale_ << "." ;
}
{}
// copy ctor
CK_Gauss(const CK_Gauss& other)
FInterpMatrixKernelGauss(const FInterpMatrixKernelGauss& other)
: lengthScale_(other.lengthScale_)
{}
// ID accessor
static const char* getID() { return "GAUSS"; }
static void printInfo() { std::cout << "K(x,y)=exp(-0.5*(r_i/l*r_i/l)) with r=|x-y|" << std::endl; }
// returns position in reduced storage
int getPosition(const unsigned int) const
{return 0;}
// returns coefficient of mutual interaction
// 1 for symmetric kernels
// -1 for antisymmetric kernels
// somethings else if other property of symmetry
FReal getMutualCoefficient() const{ return FReal(1.); }
/*
* r(x)=exp(-(|x|/l)^2)
*/
......@@ -158,10 +164,9 @@ struct CK_Gauss : AbstractCorrelationKernel<FReal>
ValueClass block[1], ValueClass blockDerivative[3]) const
{
block[0]=this->evaluate(x1,y1,z1,x2,y2,z2);
// derivative not needed
blockDerivative[0] = FMath::Zero<ValueClass>();
blockDerivative[1] = FMath::Zero<ValueClass>();
blockDerivative[2] = FMath::Zero<ValueClass>();
blockDerivative[0] = -block[0]*(x1-x2)/(lengthScale_*lengthScale_);
blockDerivative[1] = -block[0]*(y1-y2)/(lengthScale_*lengthScale_);
blockDerivative[2] = -block[0]*(z1-z2)/(lengthScale_*lengthScale_);
}
/*
......@@ -179,10 +184,6 @@ struct CK_Gauss : AbstractCorrelationKernel<FReal>
return FReal(1.);
}
// returns position in reduced storage
int getPosition(const unsigned int) const
{return 0;}
FReal evaluate(const FPoint<FReal>& p1, const FPoint<FReal>& p2) const{
return evaluate<FReal>(p1.getX(), p1.getY(), p1.getZ(), p2.getX(), p2.getY(), p2.getZ());
}
......
......@@ -31,6 +31,7 @@
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel_Covariance.hpp"
#include "Kernels/Chebyshev/FChebKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
......@@ -63,10 +64,10 @@ int main(int argc, char* argv[])
FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::NbThreads);
typedef double FReal;
const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 4);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
FTic time;
#ifdef _OPENMP
......@@ -77,9 +78,10 @@ int main(int argc, char* argv[])
#endif
// interaction kernel evaluator
//typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
//typedef FInterpMatrixKernelGauss<FReal> MatrixKernelClass;
//const MatrixKernelClass MatrixKernel(.5);
// init particles position and physical value
struct TestParticle{
......@@ -118,11 +120,11 @@ int main(int argc, char* argv[])
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential,&MatrixKernel);
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential,&MatrixKernel);
}
}
}
......@@ -134,13 +136,17 @@ int main(int argc, char* argv[])
////////////////////////////////////////////////////////////////////
{ // begin Chebyshev kernel
const unsigned int ORDER = 7;
// accuracy
const unsigned int ORDER = 7 ;
// typedefs
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
typedef FChebCell<FReal,ORDER> CellClass;
typedef FOctree<FReal,CellClass,ContainerClass,LeafClass> OctreeClass;
//typedef FChebKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FChebKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
//typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
#ifdef _OPENMP
typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#else
......@@ -170,7 +176,7 @@ int main(int argc, char* argv[])
{ // -----------------------------------------------------
std::cout << "\nChebyshev FMM (ORDER="<< ORDER << "... " << std::endl;
time.tic();
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel,FReal(1.e-15));
FmmClass algorithm(&tree, &kernels);
algorithm.execute();
time.tac();
......@@ -211,6 +217,9 @@ int main(int argc, char* argv[])
} // end Chebyshev kernel
std::cout << "Memory used at the end " << FMemStats::controler.getCurrentAllocated() << " Bytes (" << FMemStats::controler.getCurrentAllocatedMB() << "MB)\n";
std::cout << "Max memory used " << FMemStats::controler.getMaxAllocated() << " Bytes (" << FMemStats::controler.getMaxAllocatedMB() << "MB)\n";
std::cout << "Total memory used " << FMemStats::controler.getTotalAllocated() << " Bytes (" << FMemStats::controler.getTotalAllocatedMB() << "MB)\n";
return 0;
}
......@@ -38,7 +38,7 @@
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel_Covariance.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Components/FSimpleLeaf.hpp"
......@@ -68,11 +68,12 @@ int main(int argc, char* argv[])
FHelpDescribeAndExit(argc, argv,
"Test Uniform kernel and compare it with the direct computation.",
FParameterDefinitions::OctreeHeight,FParameterDefinitions::NbThreads,
FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile);
FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
FParameterDefinitions::SeparationCriterion);
typedef double FReal;
const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 3);
const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 4);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
......@@ -87,9 +88,9 @@ int main(int argc, char* argv[])
FTic time;
// interaction kernel evaluator
// typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
typedef FInterpMatrixKernelAPLUSRR<FReal> MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
typedef FInterpMatrixKernelGauss<FReal> MatrixKernelClass;
const FReal kernelParameter=FReal(0.5);
const MatrixKernelClass MatrixKernel(kernelParameter);
// init particles position and physical value
struct TestParticle{
......@@ -117,23 +118,31 @@ int main(int argc, char* argv[])
particles[idxPart].forces[2] = 0.0;
}
////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
{ // begin direct computation
std::cout << "\nDirect computation ... " << std::endl;
std::cout << "\nDirect computation (including self-interactions) ... " << std::endl;
time.tic();
{
for(FSize idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(FSize idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
FP2P::NonMutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&MatrixKernel);
for(FSize idxOther = idxTarget+1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential,
&MatrixKernel);
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential,
&MatrixKernel);
}
}
}
......@@ -148,7 +157,7 @@ int main(int argc, char* argv[])
{ // begin Lagrange kernel
// accuracy
const unsigned int ORDER = 5 ;
const unsigned int ORDER = 7 ;
// typedefs
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
......@@ -156,14 +165,16 @@ int main(int argc, char* argv[])
typedef FUnifCell<FReal,ORDER> CellClass;
typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FUnifKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
#ifdef _OPENMP
typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#else
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#endif
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
// Separation criterion for the leaf clusters
const int LeafLevelSeparationCriterion = FParameters::getValue(argc, argv, FParameterDefinitions::SeparationCriterion.options, 1);
const int LeafLevelSeparationCriterion = -1;//FParameters::getValue(argc, argv, FParameterDefinitions::SeparationCriterion.options, 1);
{ // -----------------------------------------------------
......@@ -183,7 +194,7 @@ int main(int argc, char* argv[])
} // -----------------------------------------------------
{ // -----------------------------------------------------
std::cout << "\nLagrange/Uniform grid FMM (ORDER="<< ORDER << ") ... " << std::endl;
std::cout << "\nUniform grid FMM (ORDER="<< ORDER << ") ... " << std::endl;
time.tic();
KernelClass* kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel, LeafLevelSeparationCriterion);
FmmClass algorithm(&tree, kernels, LeafLevelSeparationCriterion);
......@@ -198,8 +209,6 @@ int main(int argc, char* argv[])
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz;
FReal checkPotential[20000];
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
......@@ -212,8 +221,6 @@ int main(int argc, char* argv[])
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize indexPartOrig = indexes[idxPart];
//PB: store potential in nbParticles array
checkPotential[indexPartOrig]=potentials[idxPart];
potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
......
......@@ -39,6 +39,7 @@
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel_Covariance.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Components/FSimpleLeaf.hpp"
......@@ -69,7 +70,7 @@ int main(int argc, char* argv[])
typedef double FReal;
const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 3);
const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 4);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
......@@ -86,6 +87,8 @@ int main(int argc, char* argv[])
// interaction kernel evaluator
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
//typedef FInterpMatrixKernelGauss<FReal> MatrixKernelClass;
//const MatrixKernelClass MatrixKernel(.5);
// init particles position and physical value
struct TestParticle{
......@@ -124,12 +127,12 @@ int main(int argc, char* argv[])
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential,
&MatrixKernel);
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential,
&MatrixKernel);
}
}
}
......@@ -144,7 +147,7 @@ int main(int argc, char* argv[])
{ // begin Lagrange kernel
// accuracy
const unsigned int ORDER = 5 ;
const unsigned int ORDER = 7 ;
// typedefs
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
......@@ -152,9 +155,11 @@ int main(int argc, char* argv[])
typedef FUnifCell<FReal,ORDER> CellClass;
typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FUnifKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
#ifdef _OPENMP
typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#else
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#endif
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
......@@ -176,7 +181,7 @@ int main(int argc, char* argv[])
} // -----------------------------------------------------
{ // -----------------------------------------------------
std::cout << "\nLagrange/Uniform grid FMM (ORDER="<< ORDER << ") ... " << std::endl;
std::cout << "\nUniform grid FMM (ORDER="<< ORDER << ") ... " << std::endl;
time.tic();
KernelClass* kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
FmmClass algorithm(&tree, kernels);
......@@ -191,8 +196,6 @@ int main(int argc, char* argv[])
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz;
FReal checkPotential[20000];
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
......@@ -205,8 +208,6 @@ int main(int argc, char* argv[])
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize indexPartOrig = indexes[idxPart];
//PB: store potential in nbParticles array
checkPotential[indexPartOrig]=potentials[idxPart];
potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
......
......@@ -117,7 +117,7 @@ int main(int argc, char ** argv){
////////////////////////////////////////////////////////////////////
// approximative computation
const unsigned int ORDER = 4;
const unsigned int ORDER = 2;
const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
typedef FUnifInterpolator<FReal,ORDER,MatrixKernelClass> InterpolatorClass;
InterpolatorClass S;
......@@ -185,25 +185,35 @@ int main(int argc, char ** argv){
TensorType::setNodeIds(node_ids);
unsigned int node_ids_pairs[rc][2];
TensorType::setNodeIdsPairs(node_ids_pairs);
unsigned int permC[rc];
TensorType::setStoragePermutation(permC);
// std::cout << "2*order-1=" << 2*ORDER-1 << std::endl;
unsigned int ido=0;
for(unsigned int l=0; l<2*ORDER-1; ++l)
for(unsigned int m=0; m<2*ORDER-1; ++m)
for(unsigned int l=0; l<2*ORDER-1; ++l){
for(unsigned int m=0; m<2*ORDER-1; ++m){
for(unsigned int n=0; n<2*ORDER-1; ++n){
std::cout<< permC[ido] << ", ";
C[ido]
= MatrixKernel.evaluate(rootsX[node_ids_pairs[ido][0]],
rootsY[node_ids_pairs[ido][1]]);
C[permC[ido]] = MatrixKernel.evaluate(rootsX[node_ids_pairs[ido][0]], rootsY[node_ids_pairs[ido][1]]);
ido++;
}
// // Display C (gathers every values of K that need to be stored,
// // corresponds to the first line of the padded matrix (reverse order?))
// std::cout<<"C="<<std::endl;
// for (unsigned int n=0; n<rankpad; ++n)
// std::cout<< C[ido] << ", ";
// std::cout<<std::endl;
}
}
// Display C (gathers every values of K that need to be stored,
// corresponds to the first line of the padded matrix (reverse order?))
std::cout<<"C="<<std::endl;
ido=0;
for(unsigned int l=0; l<2*ORDER-1; ++l){
for(unsigned int m=0; m<2*ORDER-1; ++m){
for(unsigned int n=0; n<2*ORDER-1; ++n){
std::cout<< C[ido] << ", ";
ido++;
}
}
}
std::cout<<std::endl;
//////////////////////////////////////////////////////////////////////////////////////////////////////
// K is a block Toeplitz matrix
......@@ -217,14 +227,14 @@ int main(int argc, char ** argv){
// In order to actually embed K into a circulant matrix C one just
// needs to insert (ORDER-1) extra lines/columns (to each block).
// std::cout<< "K=" <<std::endl;
// for (unsigned int i=0; i<nnodes; ++i){
// for (unsigned int j=0; j<nnodes; ++j){
// std::cout<< K[i*nnodes+j]<<", ";
// }
// std::cout<<std::endl;
// }
// std::cout<<std::endl;
std::cout<< "K=" <<std::endl;
for (unsigned int i=0; i<nnodes; ++i){
for (unsigned int j=0; j<ORDER*ORDER*ORDER/*nnodes*/; ++j){
std::cout<< K[i*nnodes+j]<<", ";
}
std::cout<<std::endl;
}
std::cout<<std::endl;
if(ORDER<4){// display some extra results for low orders
......
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