# Tag Info

13

If you have constants that will not change before runs, declare them in a header file: //constants.hpp #ifndef PROJECT_NAME_constants_hpp #define PROJECT_NAME_constants_hpp namespace constants { constexpr double G = 6.67408e-11; constexpr double M_EARTH = 5.972e24; constexpr double GM_EARTH = G*M_EARTH; } #endif //main.cpp using namespace ...

12

I was able to reproduce the behavior reported in the question, and traced the observed inaccuracies to the following line: return y*sin(pi<Real>()*x)/pi<Real>(); The explicit multiplication with a floating-point approximation of π introduces a small error into the argument to sin, which comprises the representational error in the constant and ...

12

Another alternative that may be in line with your train of thought is to use a namespace (or nested namespaces) to properly group constants. An example might be: namespace constants { namespace earth { constexpr double G = 6.67408e-11; constexpr double Mass_Earth = 5.972e24; constexpr double GM = G*Mass_Earth; }// constant properties ...

7

The 14th order methods due to Feagin can be found in DifferentialEquations.jl. An example using them with 128-bit floating point arithmetic is as follows: using OrdinaryDiffEq, DoubleFloats function lorenz(du,u,p,t) du[1] = 10.0(u[2]-u[1]) du[2] = u[1]*(28.0-u[3]) - u[2] du[3] = u[1]*u[2] - (8/3)*u[3] end u0 = [1.0;0.0;0.0] tspan = (0.0,100.0) prob = ...

5

I have looked at your simple code example, and my suspicion is that what you observe in loss of speed is due to the C-heritage requirement that right-hand side expressions be evaluated using intermediate double precision operations, even when the source variables are single precision and the output location (left-hand side destination) is single precision. ...

4

Numpy has a file format that is pretty simple, which makes it perfectly compatible with basically every other high level language. (https://www.numpy.org/devdocs/reference/generated/numpy.lib.format.html) It looks like the format is much lighter than boost or hdf5. The docs say it should be easy enough to write a parser yourself, if necessary, which I would ...

4

I don't recommend it, but not because of performance issues. It will be a little less performant than a traditional matrix, which are usually allocated as a big chunk of contiguous data that is indexed using a single pointer dereference and integer arithmetic . The reason for the performance hit is mostly caching differences, but once your matrix size gets ...

4

There are multiple possible explanations: with floats and doubles the different number of computations happen due to, say, number of iterations/function evaluations (or something else, as pointed out by @Richard in the comments) there is some type conversions/implementations (say templated code) that are non-optimal with float types as opposed to doubles. I ...

3

I can only give you my personal impressions and opinions. My impression is that many people find the class somewhat redundant, because functionally it doesn't really give you much that the vector class doesn't. The algorithms it implements will have the same asymptotic complexity, although optimizations can in principle be made, because it restricts the type ...

3

I think you need to know that when you discretize a volume by using tetrahedral meshes, you will get just an approximation of your perfect surfaces because of triangulation. You started with STEP file format, which is NOT a mesh format, but it's a CAD format. The difference between CAD formats and mesh formats is that CAD formats like STEP, IGES, SAT, etc. ...

3

The goal is to compute an $x\in\mathbb{R}^{n}$ to satisfy the strict inequalities $$f(z_{i})^{T}x>0\qquad\forall z_{i}\in\{1,\ldots,m\}.\tag{1}$$ If we write $\epsilon>0$ as the smallest margin of feasibility, that is $$\epsilon=\min_{i}\{f(z_{i})^{T}x\},$$ then (1) is equivalent to $$f(z_{i})^{T}x\ge\epsilon>0\qquad\forall z_{i}\in\{1,\ldots,... 3 So, you are comparing a generally slower Matlab implementation of algorithm A to a generally faster C++ implementation of algorithm B, and still getting the advantage for A. I would say, congratulations, you certainly have a stronger point now, since the "competing" algorithm is given an advantage. It's worth no note though, that implementation in C++ is ... 2 Using the recurrence step and the normalized incomplete gamma function definition I simplified the formula to:$$\theta e^{k \log \left(\frac{o}{\theta }\right)-\frac{o}{\theta }- \ln{\Gamma (\kappa)}} +\theta \kappa Q \left(\kappa,\frac{o}{\theta }\right)-o Q \left(\kappa,\frac{o}{\theta }\right)+o+s$$In this equation Q stands for the normalized ... 2 The library Boost.Numeric.uBlas recently added a tensor extension which is shipped with Boost version 1.70. Please have a look at https://github.com/boostorg/ublas. It provides standard matrix and tensor operations with runtime-variable order (number of dimensions), dimensions for the first- and last-order storage formats (column- and row-major). You can ... 2 Cereal is a straight-forward and easy-to-use serialization library for C++ data structures. It knows all of the STL data structures and adapting custom classes for use with it is easy. Example code: #include <cereal/types/unordered_map.hpp> #include <cereal/types/memory.hpp> #include <cereal/archives/binary.hpp> #include <fstream> ... 2 You can use boost's serialization facilities, which can be straightforwardly extended to support custom data structures. If you use Eigen you can adapt something like this in order to suit your needs. 2 Before you start down this path it's important to determine whether there's enough data parallelism in your current code to make using a GPU worthwhile. I'd encourage you to start by describing your application and algorithms in more detail in the question. Depending on the application, there may be computational tasks for which library routines are ... 2 It is recommended to think about parallelization first and then discuss the implementation. Think about what the code does, what data dependencies exist, and what operations can be carried out in parallel. Then there are several C++ frameworks (alpaka, kokkos, or the ArrayFire library mentioned in another answer) that help you to introduce a layer of ... 2 OpenCL is runnable on multicore cpu, intel hd graphics and even DSP cards. It was pretty much the standard for cross platform gpu computing until compute shaders were introduced. There are various libraries that have OpenCL as a backend such as viennaCL or ArrayFire. Some of these libraries can use other backends for gpu computation such as CUDA, which runs ... 2 Are you sure it is supposed to be sign function ? According to this https://projects.coin-or.org/ADOL-C/browser/stable/2.1/ADOL-C/doc/adolc-manual.pdf?format=raw in section 1.8, condassign(a,b,c,d) is equal to a = (b > 0) ? c : d So the function posted is actually giving Heaviside function and not the sign function. It implements this function$$ f(x) =...

2

Here are two considerations that might help you figure out whether a solution exists: First, if there is a vector $X$ so that $F(z_i) \cdot X > 0$, then all vectors $Y=\alpha X$ for any choice of $\alpha>0$ are also solutions because $F(z_i) \cdot Y = \alpha \underbrace{(F(z_i)\cdot X)}_{\ge 0} > 0$. In other words, if a solution exists, then there ...

2

The solution I found for PETSc 3.7 was to use MatGetSubMatrices. Warning : this method is not present in the last API. I discovered that BAIJ is actually not intended for block matrices but rather to keep rows by group when distributing the matrix over processors. Thus I used AIJ format. (MPIAIJ to get a parallel matrix). First have the rights rows and ...

1

For anyone interested in this problem, I found the following solution: #include "itensor/all.h" using namespace itensor; int main() { int N = 100; // // Initialize the site degrees of freedom. // auto sites = SpinHalf(N,{"ConserveQNs=",false}); //make a chain of N spin 1/2's //Transverse field Real h = 4.0; // //...

1

I would second the opinion expressed in the context. For that small problem and limited usage, you don't need a library. Generation of a structured grid and retrieving points can be coded up in at most several screens of code. However, I would point out for you ViennaGrid library, which can be used and actually provides STL iterators. In addition, the ...

1

Definitely use C++ and Eigen. Reasons: Sounds like it fits best for your case. Eigen beats LAPACK using optimization flags (-O3) and a good compiler (GCC, Clang). At least for most tests I recently performed (dense linear algebra) Once you debug your application use the -ffast-math compiler flag for a huge speedup, but test your output to see if it remains ...

1

One way that I do is to use singleton. When you start your program you initiate your singleton and fill it with the constant data (probably from a properties file that you have for the run). You get this in every class that you need the values and just use it.

1

OP's comment on OP: yes,numbers are intergers For arbitrarily large integer calculations, MAPLE is an option worth trying. it has very efficient implementations for a number of such calculations.

1

If you have a structured way of enumerating the graphs you wish to calculate the determinants of, perhaps you could find low-rank updates which transfer you from one graph to another. If so, then you could use the matrix determinant lemma to cheaply calculate the determinant of the subsequent graph to be enumerated using your knowledge of the current graph'...

1

I also write a benchmark. For matrix of small size (<100*100), the performance is similar for vector< vector< double>> and wrapped 1D vector. For matrix of large size (~1000*1000), wrapped 1D vector is better. The Eigen matrix behaves worse. It is surprise to me that the Eigen is the worst. #include <iostream> #include <iomanip> #...

Only top voted, non community-wiki answers of a minimum length are eligible