Hot answers tagged

11

The way I have been measuring whether the eco-system is ready is how things are going with the transition for the homebrew package manager. They have been carefully documenting the progress of getting things running on apple silicon via a github issue. Beyond @Federico Poloni's point the biggest problem is that GCC itself is not yet working and is the ...


9

In broad terms, algorithms that run faster on the GPU are ones where you are doing the same type of instruction on many different data points. An easy example to illustrate this is with matrix multiplication. Suppose we are doing the matrix computation $A \times B = C$ A simple CPU algorithm might look something like //starting with C = 0 for (...


9

You should consider giving Julia a try. Let me explain what's going on in the design space right now that would be of interest to you. Full disclosure I am the lead developer of JuliaDiffEq. JuliaDiffEq and DifferentialEquations.jl has a large feature set dedicated to efficiently integrating computationally-difficult differential equations. It has a simple ...


9

From this, it looks like there is no functional native Fortran compiler yet. If that is really the case, things look bleak. Almost anything that uses linear algebra includes some Fortran code (Lapack), and it has to run fast.


8

The benchmarks on the Julia website 1 2 include R and Matlab as competitors. Note that these are benchmarks focusing on testing the pure speed of the language, not the quality of the underlying linear algebra or FFT libraries. The speed for operations that are outsourced to these libraries (such as a large matrix multiplication) can vary a lot depending on ...


7

I don't know of any libraries that are still in use today but the keyword in literature searches you will be looking for is "out of memory solver" or "out of core solver" -- linear solvers (and LU decompositions) that work on matrices stored on disk (or tape, at the time) were quite popular in the 1960s, 70s and 80s when memory was expensive and small. That ...


6

GPUs are sensitive beasts. Although Nvidia's beefiest card can theoretically execute any of the operations you listed 100x faster than the fastest CPU, about a million things can get in the way of that speedup. Every part of the relevant algorithm, and of the program which runs it, has to be extensively tweaked and optimized in order to get anywhere near ...


5

Which storage scheme is optimal depends very strongly on the access pattern. You already said that you're accessing them one particle at a time, but are you going to access the particles sequentially in order, or randomly? Another conventional name for this issue is struct-of-arrays vs array-of-structs. The key thing here is to know how the memory hierarchy ...


5

The problem posed is a multiobjective optimization problem, and the usual notion of optimality for these types of problems is Pareto optimality. Scalarization (as proposed in the comments by ChristianClason, TheNobleSunfish, Paul, and DougLipinski) is one way to solve the problem. This approach leverages the large body of theory and algorithms for single ...


5

Your matrix does not have an inverse, which is why dgesv returns an error.


5

The various Fortran standards allow a lot of compiler dependent behaviour in terms of function binary interfaces when being called with "complicated" data types such as Fortran90 style arrays and complex numbers. This means calling code compiled with one compiler from another one is not guaranteed to do what you expect, and can lead to grabbing the wrong bit ...


5

Iain Sandoe has been working on porting both gcc and gFortran to this architecture. Based on this, François-Xavier Coudert has created an experimental gFortran release for the M1. I provide a few Fortran benchmarks here. Presumably this was used for the miniforge Python release that provides a M1 native numpy. One also suspects this is what leman used to ...


4

For all of the applications you mentioned, GPUs should be more capable (from a hardware perspective) than CPUs for sufficiently large matrices. I don't know anything about R's implementation, but I've used cuBLAS and Magma with great success for inversions around $n = 2^{10}$ and multiplication/correlation for rectangular matrices with $n,m \approx 2^{10}, ...


3

I know time may be limited for you, but if you have ~1-3 weeks to really learn CUDA, I highly recommend CUDA by Example: An Introduction to General-Purpose GPU Programming (Amazon.com). It does a fantastic job of explaining the general concepts of GPU programming, and does a great job of getting one up to speed with NVIDIA's GPU language, CUDA. CUDA is ...


3

If time is of grave concern, I would highly suggest looking at Intels Xeon phi coprocessor. Not only are they nearly or as fast, only require openmp to code for, but Intels customer service on the Intel developer forums is fantastic. I don't know if you can use R, but standards languages such as c, c++ , and fortran can be used. You could also use Intels mkl ...


3

For the sake of completeness of the answer, you can use ls() to list all the assigned variables and data frames. This is similar to who and whos in Matlab. If you want to get the list data.frames only you can use something like this: ls()[grepl('data.frame', sapply(ls(), function(x) class(get(x))))] You should probably ask questions like this on R forums.


3

As AdamO said, you have an issue with non-positive values. Purely for optimization purposes, exponentiating to make sure things are non-negative and adding some salt to avoid zero generally does the trick. So... something like: myfun<-function(par, x){ par <- exp(par) + 10^-10 f<- sum(x)*length(x)+sum(log(gamma(par))*x)+1 return(-f) } Then none ...


3

I suggest to read the Quadpack book (Quadpack, Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. ISBN 978-3-540-12553-2) and Pedro Gonnet's PhD thesis "Adaptive Quadrature Re-revisited (available as pdf here). Pedro is a contributor to ...


3

I did not yet succeed to use the sde.sim() function of the SDE package, however, I succeeded to solve the system (with and without noise) using the suggestions of Chris and the diffeqr package in R. library(plotly) # Lotka-Volterra Model (SDE without noise) f <- function(u,p,t) { du1 = p[1]*u[1]-p[2]*u[1]*u[2] du2 = p[3]*u[1]*u[2]-p[4]*u[2] return(...


3

The parameter can be any type, so here I pass in a time-dependent function for p and use it in the differential equation: # Packages library(tidyverse) library(diffeqr) library(JuliaCall) diffeq_setup() # Drift function f <- function(u,p,t){ du1 = p(t) return(c(du1)) } # Diffusion function g <- function(u,p,t){ du1 = 0 # note that there is ...


2

The AIC function need an 'lm' or 'glm' object (linear models). See functions lm and glm So just do : AIC(lm(logCPK~dataPOW))


2

Actually, there are two considerable packages in R to handle neural networks with easiness. Here they are: nnet and neuralnet. Install them via install.packages('nnet') install.packages('neuralnet') in R. To get help and see examples, see ?neuralnet::neuralnet ?nnet::nnet You can look at a neural network as a function f(x) where x is a vector of inputs. ...


2

Why do it yourself when experts have already done it for you? See polynom and polynomF for specific implementations in R.


2

R, Matlab, and Armadillo (Mat class) matrices use column-major order so you should create 3xN matrices so that each set of 3 values is stored consecutively in memory. A bit more info: You should store your matrices such that elements which are accessed together are stored sequentially in memory. The memory order of matrix storage is language, data structure,...


2

There is a tutorial on GPU computing in R at r-tutor.com. It has various examples you can look at and primarily uses the RPUD package which is open source and also makes use of the non-free RPUDPLUS. Additionally this website has a discussion of a few different packages that aid in GPU computing in R. The packages mentioned are gputools HiPLARM rpud magma ...


2

Without knowing anything about inverse interpolation, this problem should be a straight-forward numerical task given the information about $g(x)$. Firstly, $g(x)$ is said to be monotonously increasing. Therefore, given $y$, find two nodes s.t. $g(x_i)\leq y\leq g(x_{i+1})$ [this is simple since all is ordered]. From there, use any root-finder which should ...


2

We decided to go with C++, because it is free is fast has a freely available and widely supported numerical integration library (boost/odeint) is easy to integrate with our larger simulation driver/GUI is widely known makes standalone simulation easily available to biological scientists (because they can use RStudio as an interactive simulation environment ...


2

It seems hard to write a general enough FD library thas has wide applicability, since FD methods are not as easy to write for general domains, unlike FEM which uses unstructured grids, for which there is a standard approach. I know of two libraries which might be useful for you Overture: An Object-Oriented Toolkit for Solving Partial Differential Equations ...


1

If you want to know more about efficient numerical procedures to compute Lagrange polynomials, I recommend the following reference. J.-P. Berrut and L. N. Trefethen, “Barycentric Lagrange Interpolation,” SIAM Rev., vol. 46, no. 3, pp. 501–517, 2004. The main point of the "barycentric" Lagrange formula is to pre-compute barycentric weights (in $O(n^2$) ...


1

In this line: EuclideanDist = 1 / sqrt(sum(EuclideanDist^2 )) the sqrt should not be present: the integrand is one over the squared norm, not one over norm. In this line: arcLength = integrate(FUN, x,y)$value this will give you a negative (incorrect) value for arc length whenever $x>y$, which confuses the rest of the code. The correct evaluation ...


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