Commit a991bebc authored by COULAUD Olivier's avatar COULAUD Olivier

Add CMAKE_CXX_STANDARD, fix some issues detected by sonarQube

parent 26e3813d
......@@ -13,7 +13,6 @@ set(FUSE_LIST " MPI;BLAS;FFT;STARPU;CUDA;OPENCL;OMP4;SSE;AVX;AVX2;MIC;MPI2")
# Project Declaration
#===========================================================================
project(SCALFMM C CXX )
# check if compiling into source directories
string(COMPARE EQUAL "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_BINARY_DIR}" insource)
if(insource)
......@@ -28,7 +27,7 @@ INCLUDE(CMakeDependentOption)
# Add to check CPU info
include(GetCpuInfos)
GetCpuInfos()
#
#===========================================================================
# Version Number
#===========================================================================
......@@ -39,6 +38,10 @@ set(SCALFMM_MINOR_VERSION 0)
set(SCALFMM_PATCH_VERSION rc0)
set(SCALFMM_VERSION "${SCALFMM_MAJOR_VERSION}.${SCALFMM_MINOR_VERSION}.${SCALFMM_PATCH_VERSION}" )
SET(CXX_STANDARD_REQUIRED ON)
SET(CMAKE_CXX_STANDARD 14)
set( MORSE_DISTRIB_DIR "" CACHE PATH "Directory of MORSE distribution")
if (MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/CMakeModules/morse/")
......@@ -140,7 +143,7 @@ if (MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/CMakeModules/morse/
message(STATUS "CXX ${CMAKE_CXX_COMPILER_ID}" )
# Set scalfmm to default libraries
set(SCALFMM_LIBRARIES "")
set(SCALFMM_CXX_FLAGS "-std=c++14 -fpic -Wall")
set(SCALFMM_CXX_FLAGS "-fpic -Wall")
set(SCALFMM_CXX_FLAGS "${SCALFMM_CXX_FLAGS} -I${CMAKE_CURRENT_SOURCE_DIR}/Contribs")
#
#
......
......@@ -18,7 +18,7 @@
template <class FReal, int ORDER>
unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
FReal* &U, FReal* &C, FReal* &B);
FReal* &U, FReal* &C, FReal* &B);
/**
......@@ -53,144 +53,150 @@ unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
template <class FReal, int ORDER, class MatrixKernelClass>
class FChebM2LHandler : FNoCopyable
{
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes,
ninteractions = 316}; // 7^3 - 3^3 (max num cells in far-field)
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes,
ninteractions = 316}; // 7^3 - 3^3 (max num cells in far-field)
const MatrixKernelClass *const MatrixKernel;
const MatrixKernelClass *const MatrixKernel;
FReal *U, *C, *B;
const FReal epsilon; //<! accuracy which determines trucation of SVD
unsigned int rank; //<! truncation rank, satisfies @p epsilon
FReal *U, *C, *B;
const FReal epsilon; //<! accuracy which determines trucation of SVD
unsigned int rank; //<! truncation rank, satisfies @p epsilon
static const std::string getFileName(FReal epsilon)
{
const char precision_type = (typeid(FReal)==typeid(double) ? 'd' : 'f');
std::stringstream stream;
stream << "m2l_k"<< MatrixKernelClass::getID() << "_" << precision_type
<< "_o" << order << "_e" << epsilon << ".bin";
return stream.str();
}
static const std::string getFileName(FReal epsilon)
{
const char precision_type = (typeid(FReal)==typeid(double) ? 'd' : 'f');
std::stringstream stream;
stream << "m2l_k"<< MatrixKernelClass::getID() << "_" << precision_type
<< "_o" << order << "_e" << epsilon << ".bin";
return stream.str();
}
public:
FChebM2LHandler(const MatrixKernelClass *const inMatrixKernel, const FReal _epsilon)
: MatrixKernel(inMatrixKernel), U(nullptr), C(nullptr), B(nullptr), epsilon(_epsilon), rank(0)
{}
~FChebM2LHandler()
{
if (U != nullptr) delete [] U;
if (B != nullptr) delete [] B;
if (C != nullptr) delete [] C;
}
/**
* Computes, compresses and sets the matrices \f$Y, C_t, B\f$
*/
void ComputeAndCompressAndSet()
{
// measure time
FTic time; time.tic();
// check if aready set
if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
rank = ComputeAndCompress(MatrixKernel, epsilon, U, C, B);
unsigned long sizeM2L = 343*rank*rank*sizeof(FReal);
// write info
std::cout << "Compressed and set M2L operators (" << long(sizeM2L) << " B) in "
<< time.tacAndElapsed() << "sec." << std::endl;
}
/**
* Computes, compresses, writes to binary file, reads it and sets the matrices \f$Y, C_t, B\f$
*/
void ComputeAndCompressAndStoreInBinaryFileAndReadFromFileAndSet()
{
FChebM2LHandler<FReal, ORDER,MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(epsilon);
this->ReadFromBinaryFileAndSet();
}
/**
* Computes and compressed all \f$K_t\f$.
*
* @param[in] epsilon accuracy
* @param[out] U matrix of size \f$\ell^3\times r\f$
* @param[out] C matrix of size \f$r\times 316 r\f$ storing \f$[C_1,\dots,C_{316}]\f$
* @param[out] B matrix of size \f$\ell^3\times r\f$
*/
static unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel, const FReal epsilon, FReal* &U, FReal* &C, FReal* &B);
/**
* Computes, compresses and stores the matrices \f$Y, C_t, B\f$ in a binary
* file
*/
static void ComputeAndCompressAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel, const FReal epsilon);
/**
* Reads the matrices \f$Y, C_t, B\f$ from the respective binary file
*/
void ReadFromBinaryFileAndSet();
/**
* @return rank of the SVD compressed M2L operators
*/
unsigned int getRank() const {return rank;}
FChebM2LHandler(const MatrixKernelClass *const inMatrixKernel, const FReal _epsilon)
: MatrixKernel(inMatrixKernel), U(nullptr), C(nullptr), B(nullptr), epsilon(_epsilon), rank(0)
{}
~FChebM2LHandler()
{
if (U != nullptr)
{delete [] U;}
if (B != nullptr)
{ delete [] B;}
if (C != nullptr)
{delete [] C; }
}
/**
* Computes, compresses and sets the matrices \f$Y, C_t, B\f$
*/
void ComputeAndCompressAndSet()
{
// measure time
FTic time;
time.tic();
// check if aready set
if (U||C||B) {
throw std::runtime_error("Compressed M2L operator already set"); }
rank = ComputeAndCompress(MatrixKernel, epsilon, U, C, B);
unsigned long sizeM2L = 343*rank*rank*sizeof(FReal);
// write info
std::cout << "Compressed and set M2L operators (" << long(sizeM2L) << " B) in "
<< time.tacAndElapsed() << "sec." << std::endl;
}
/**
* Computes, compresses, writes to binary file, reads it and sets the matrices \f$Y, C_t, B\f$
*/
void ComputeAndCompressAndStoreInBinaryFileAndReadFromFileAndSet()
{
FChebM2LHandler<FReal, ORDER,MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(epsilon);
this->ReadFromBinaryFileAndSet();
}
/**
* Computes and compressed all \f$K_t\f$.
*
* @param[in] epsilon accuracy
* @param[out] U matrix of size \f$\ell^3\times r\f$
* @param[out] C matrix of size \f$r\times 316 r\f$ storing \f$[C_1,\dots,C_{316}]\f$
* @param[out] B matrix of size \f$\ell^3\times r\f$
*/
static unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel, const FReal epsilon, FReal* &U, FReal* &C, FReal* &B);
/**
* Computes, compresses and stores the matrices \f$Y, C_t, B\f$ in a binary
* file
*/
static void ComputeAndCompressAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel, const FReal epsilon);
/**
* Reads the matrices \f$Y, C_t, B\f$ from the respective binary file
*/
void ReadFromBinaryFileAndSet();
/**
* @return rank of the SVD compressed M2L operators
*/
unsigned int getRank() const {return rank;}
/**
* Expands potentials \f$x+=UX\f$ of a target cell. This operation can be
* seen as part of the L2L operation.
*
* @param[in] X compressed local expansion of size \f$r\f$
* @param[out] x local expansion of size \f$\ell^3\f$
*/
* Expands potentials \f$x+=UX\f$ of a target cell. This operation can be
* seen as part of the L2L operation.
*
* @param[in] X compressed local expansion of size \f$r\f$
* @param[out] x local expansion of size \f$\ell^3\f$
*/
void applyU(const FReal *const X, FReal *const x) const
{
FBlas::gemva(nnodes, rank, 1., U, const_cast<FReal*>(X), x);
constexpr FReal ONE = 1.0 ;
FBlas::gemva(nnodes, rank, ONE, U, const_cast<FReal*>(X), x);
}
/**
* Compressed M2L operation \f$X+=C_tY\f$, where \f$Y\f$ is the compressed
* multipole expansion and \f$X\f$ is the compressed local expansion, both
* of size \f$r\f$. The index \f$t\f$ denotes the transfer vector of the
* target cell to the source cell.
*
* @param[in] transfer transfer vector
* @param[in] Y compressed multipole expansion
* @param[out] X compressed local expansion
* @param[in] CellWidth needed for the scaling of the compressed M2L operators which are based on a homogeneous matrix kernel computed for the reference cell width \f$w=2\f$, ie in \f$[-1,1]^3\f$.
*/
* Compressed M2L operation \f$X+=C_tY\f$, where \f$Y\f$ is the compressed
* multipole expansion and \f$X\f$ is the compressed local expansion, both
* of size \f$r\f$. The index \f$t\f$ denotes the transfer vector of the
* target cell to the source cell.
*
* @param[in] transfer transfer vector
* @param[in] Y compressed multipole expansion
* @param[out] X compressed local expansion
* @param[in] CellWidth needed for the scaling of the compressed M2L operators which are based on a homogeneous matrix kernel computed for the reference cell width \f$w=2\f$, ie in \f$[-1,1]^3\f$.
*/
void applyC(const int transfer[3], FReal CellWidth,
const FReal *const Y, FReal *const X) const
const FReal *const Y, FReal *const X) const
{
const unsigned int idx
= static_cast<int>((transfer[0]+3)*7*7 + (transfer[1]+3)*7 )+ (transfer[2]+3);
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
const unsigned int idx
= static_cast<int>((transfer[0]+3)*7*7 + (transfer[1]+3)*7 )+ (transfer[2]+3);
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
}
void applyC(const unsigned int idx, FReal CellWidth,
const FReal *const Y, FReal *const X) const
const FReal *const Y, FReal *const X) const
{
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
}
void applyC(FReal CellWidth,
const FReal *const Y, FReal *const X) const
const FReal *const Y, FReal *const X) const
{
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
FBlas::gemva(rank, rank * 343, scale, C, const_cast<FReal*>(Y), X);
}
/**
* Compresses densities \f$Y=B^\top y\f$ of a source cell. This operation
* can be seen as part of the M2M operation.
*
* @param[in] y multipole expansion of size \f$\ell^3\f$
* @param[out] Y compressed multipole expansion of size \f$r\f$
*/
* Compresses densities \f$Y=B^\top y\f$ of a source cell. This operation
* can be seen as part of the M2M operation.
*
* @param[in] y multipole expansion of size \f$\ell^3\f$
* @param[out] Y compressed multipole expansion of size \f$r\f$
*/
void applyB(FReal *const y, FReal *const Y) const
{
FBlas::gemtv(nnodes, rank, 1., B, y, Y);
......@@ -216,96 +222,108 @@ public:
template <class FReal, int ORDER, class MatrixKernelClass>
unsigned int
FChebM2LHandler<FReal, ORDER, MatrixKernelClass>::ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
const FReal epsilon,
FReal* &U,
FReal* &C,
FReal* &B)
const FReal epsilon,
FReal* &U,
FReal* &C,
FReal* &B)
{
// allocate memory and store compressed M2L operators
if (U||C||B) throw std::runtime_error("Compressed M2L operators are already set");
// interpolation points of source (Y) and target (X) cell
FPoint<FReal> X[nnodes], Y[nnodes];
// set roots of target cell (X)
FChebTensor<FReal, order>::setRoots(FPoint<FReal>(0.,0.,0.), FReal(2.), X);
// allocate memory and compute 316 m2l operators
FReal *_U, *_C, *_B;
_U = _B = nullptr;
_C = new FReal [nnodes*nnodes * ninteractions];
unsigned int counter = 0;
for (int i=-3; i<=3; ++i) {
for (int j=-3; j<=3; ++j) {
for (int k=-3; k<=3; ++k) {
if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
// set roots of source cell (Y)
const FPoint<FReal> cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
FChebTensor<FReal, order>::setRoots(cy, FReal(2.), Y);
// evaluate m2l operator
for (unsigned int n=0; n<nnodes; ++n)
for (unsigned int m=0; m<nnodes; ++m)
_C[counter*nnodes*nnodes + n*nnodes + m]
= MatrixKernel->evaluate(X[m], Y[n]);
// increment interaction counter
counter++;
}
}
}
}
if (counter != ninteractions)
throw std::runtime_error("Number of interactions must correspond to 316");
//////////////////////////////////////////////////////////
FReal weights[nnodes];
FChebTensor<FReal, order>::setRootOfWeights(weights);
for (unsigned int i=0; i<316; ++i)
for (unsigned int n=0; n<nnodes; ++n) {
FBlas::scal(nnodes, weights[n], _C+i*nnodes*nnodes + n, nnodes); // scale rows
FBlas::scal(nnodes, weights[n], _C+i*nnodes*nnodes + n * nnodes); // scale cols
}
//////////////////////////////////////////////////////////
// svd compression of M2L
const unsigned int rank = Compress<FReal, ORDER>(epsilon, ninteractions, _U, _C, _B);
if (!(rank>0)) throw std::runtime_error("Low rank must be larger then 0!");
// store U
U = new FReal [nnodes * rank];
FBlas::copy(rank*nnodes, _U, U);
delete [] _U;
// store B
B = new FReal [nnodes * rank];
FBlas::copy(rank*nnodes, _B, B);
delete [] _B;
// store C
counter = 0;
C = new FReal [343 * rank*rank];
for (int i=-3; i<=3; ++i)
for (int j=-3; j<=3; ++j)
for (int k=-3; k<=3; ++k) {
const unsigned int idx = (i+3)*7*7 + (j+3)*7 + (k+3);
if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
FBlas::copy(rank*rank, _C + counter*rank*rank, C + idx*rank*rank);
counter++;
} else FBlas::setzero(rank*rank, C + idx*rank*rank);
}
if (counter != ninteractions)
throw std::runtime_error("Number of interactions must correspond to 316");
delete [] _C;
//////////////////////////////////////////////////////////
for (unsigned int n=0; n<nnodes; ++n) {
FBlas::scal(rank, FReal(1.) / weights[n], U+n, nnodes); // scale rows
FBlas::scal(rank, FReal(1.) / weights[n], B+n, nnodes); // scale rows
}
//////////////////////////////////////////////////////////
// return low rank
return rank;
// allocate memory and store compressed M2L operators
if (U||C||B){
throw std::runtime_error("Compressed M2L operators are already set");
}
// interpolation points of source (Y) and target (X) cell
FPoint<FReal> X[nnodes], Y[nnodes];
// set roots of target cell (X)
FChebTensor<FReal, order>::setRoots(FPoint<FReal>(0.,0.,0.), FReal(2.), X);
// allocate memory and compute 316 m2l operators
FReal *_U, *_C, *_B;
_U = _B = nullptr;
_C = new FReal [nnodes*nnodes * ninteractions];
//
unsigned int counter = 0;
for (int i=-3; i<=3; ++i) {
for (int j=-3; j<=3; ++j) {
for (int k=-3; k<=3; ++k) {
if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
// set roots of source cell (Y)
const FPoint<FReal> cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
FChebTensor<FReal, order>::setRoots(cy, FReal(2.), Y);
// evaluate m2l operator
for (unsigned int n=0; n<nnodes; ++n){
for (unsigned int m=0; m<nnodes; ++m){
_C[counter*nnodes*nnodes + n*nnodes + m]
= MatrixKernel->evaluate(X[m], Y[n]);
}
}
// increment interaction counter
++counter;
}
}
}
}
if (counter != ninteractions){
delete [] _C ;
throw std::runtime_error("Number of interactions must correspond to 316");
}
//////////////////////////////////////////////////////////
FReal weights[nnodes];
FChebTensor<FReal, order>::setRootOfWeights(weights);
for (unsigned int i=0; i<316; ++i)
for (unsigned int n=0; n<nnodes; ++n) {
FBlas::scal(nnodes, weights[n], _C+i*nnodes*nnodes + n, nnodes); // scale rows
FBlas::scal(nnodes, weights[n], _C+i*nnodes*nnodes + n * nnodes); // scale cols
}
//////////////////////////////////////////////////////////
// svd compression of M2L
const unsigned int rank = Compress<FReal, ORDER>(epsilon, ninteractions, _U, _C, _B);
if (!(rank>0)) {
throw std::runtime_error("Low rank must be larger then 0!");
}
// store U
U = new FReal [nnodes * rank];
FBlas::copy(rank*nnodes, _U, U);
delete [] _U;
// store B
B = new FReal [nnodes * rank];
FBlas::copy(rank*nnodes, _B, B);
delete [] _B;
// store C
counter = 0;
C = new FReal [343 * rank*rank];
for (int i=-3; i<=3; ++i)
for (int j=-3; j<=3; ++j)
for (int k=-3; k<=3; ++k) {
const unsigned int idx = (i+3)*7*7 + (j+3)*7 + (k+3);
if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
FBlas::copy(rank*rank, _C + counter*rank*rank, C + idx*rank*rank);
++counter;
}
else {
FBlas::setzero(rank*rank, C + idx*rank*rank);
}
}
if (counter != ninteractions){
throw std::runtime_error("Number of interactions must correspond to 316");
}
delete [] _C;
//////////////////////////////////////////////////////////
for (unsigned int n=0; n<nnodes; ++n) {
FBlas::scal(rank, FReal(1.) / weights[n], U+n, nnodes); // scale rows
FBlas::scal(rank, FReal(1.) / weights[n], B+n, nnodes); // scale rows
}
//////////////////////////////////////////////////////////
// return low rank
return rank;
}
......@@ -317,39 +335,42 @@ template <class FReal, int ORDER, class MatrixKernelClass>
void
FChebM2LHandler<FReal, ORDER, MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel, const FReal epsilon)
{
// measure time
FTic time; time.tic();
// start computing process
FReal *U, *C, *B;
U = C = B = nullptr;
const unsigned int rank = ComputeAndCompress(MatrixKernel, epsilon, U, C, B);
// store into binary file
const std::string filename(getFileName(epsilon));
std::ofstream stream(filename.c_str(),
std::ios::out | std::ios::binary | std::ios::trunc);
if (stream.good()) {
stream.seekp(0);
// 1) write number of interpolation points (int)
int _nnodes = nnodes;
stream.write(reinterpret_cast<char*>(&_nnodes), sizeof(int));
// 2) write low rank (int)
int _rank = rank;
stream.write(reinterpret_cast<char*>(&_rank), sizeof(int));
// 3) write U (rank*nnodes * FReal)
stream.write(reinterpret_cast<char*>(U), sizeof(FReal)*rank*nnodes);
// 4) write B (rank*nnodes * FReal)
stream.write(reinterpret_cast<char*>(B), sizeof(FReal)*rank*nnodes);
// 5) write 343 C (343 * rank*rank * FReal)
stream.write(reinterpret_cast<char*>(C), sizeof(FReal)*rank*rank*343);
} else throw std::runtime_error("File could not be opened to write");
stream.close();
// free memory
if (U != nullptr) delete [] U;
if (B != nullptr) delete [] B;
if (C != nullptr) delete [] C;
// write info
std::cout << "Compressed M2L operators ("<< rank << ") stored in binary file " << filename
<< " in " << time.tacAndElapsed() << "sec." << std::endl;
// measure time
FTic time; time.tic();
// start computing process
FReal *U, *C, *B;
U = C = B = nullptr;
const unsigned int rank = ComputeAndCompress(MatrixKernel, epsilon, U, C, B);
// store into binary file
const std::string filename(getFileName(epsilon));
std::ofstream stream(filename.c_str(),
std::ios::out | std::ios::binary | std::ios::trunc);
if (stream.good()) {
stream.seekp(0);
// 1) write number of interpolation points (int)
int _nnodes = nnodes;
stream.write(reinterpret_cast<char*>(&_nnodes), sizeof(int));
// 2) write low rank (int)
int _rank = rank;
stream.write(reinterpret_cast<char*>(&_rank), sizeof(int));
// 3) write U (rank*nnodes * FReal)
stream.write(reinterpret_cast<char*>(U), sizeof(FReal)*rank*nnodes);
// 4) write B (rank*nnodes * FReal)
stream.write(reinterpret_cast<char*>(B), sizeof(FReal)*rank*nnodes);
// 5) write 343 C (343 * rank*rank * FReal)
stream.write(reinterpret_cast<char*>(C), sizeof(FReal)*rank*rank*343);
} else throw std::runtime_error("File could not be opened to write");
stream.close();
// free memory
if (U != nullptr) {
delete [] U;}
if (B != nullptr) {
delete [] B;}
if (C != nullptr) {
delete [] C; }
// write info
std::cout << "Compressed M2L operators ("<< rank << ") stored in binary file " << filename
<< " in " << time.tacAndElapsed() << "sec." << std::endl;
}
......@@ -357,42 +378,50 @@ template <class FReal, int ORDER, class MatrixKernelClass>
void
FChebM2LHandler<FReal, ORDER, MatrixKernelClass>::ReadFromBinaryFileAndSet()
{
// measure time
FTic time; time.tic();
// start reading process
if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
const std::string filename(getFileName(epsilon));
std::ifstream stream(filename.c_str(),
std::ios::in | std::ios::binary | std::ios::ate);
const std::ifstream::pos_type size = stream.tellg();
if (size<=0) {
std::cout << "Info: The requested binary file " << filename
<< " does not yet exist. Compute it now ... " << std::endl;
this->ComputeAndCompressAndStoreInBinaryFileAndReadFromFileAndSet();
return;
}
if (stream.good()) {
stream.seekg(0);
// 1) read number of interpolation points (int)
int npts;
stream.read(reinterpret_cast<char*>(&npts), sizeof(int));
if (npts!=nnodes) throw std::runtime_error("nnodes and npts do not correspond");
// 2) read low rank (int)
stream.read(reinterpret_cast<char*>(&rank), sizeof(int));
// 3) write U (rank*nnodes * FReal)
U = new FReal [rank*nnodes];
stream.read(reinterpret_cast<char*>(U), sizeof(FReal)*rank*nnodes);
// 4) write B (rank*nnodes * FReal)
B = new FReal [rank*nnodes];
stream.read(reinterpret_cast<char*>(B), sizeof(FReal)*rank*nnodes);