People who code: we want your input. Take the Survey
11

Such an effect happens because of how the data of the int** a is stored in memory (as per C/C++). This question on StackOverflow has answers with some more details (in particular, a difference between int** and int[][] that many users noted in the comments), and how it looks like an array of arrays - it's just laid out contiguously in memory. It's worth ...


7

You have a size mismatch issue: A is a count x 2*count matrix, and you are trying to solve Ax=B with B a 2*count x 1 vector. Moreover, if you compile without -DNDEBUG, you should get a nice assertion telling you this is wrong. To resize a matrix or vector: A.resize(count, 2*count); B.resize(count);


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 ...


5

I can't give you a good answer for 1, but I can give you decent answers for 2 and 3. I'm not terribly familiar with the Boost interface to UMFPACK. In the C interface, normally, you call UMFPACK routines that will allocate memory. If the memory cannot be allocated because there is not enough free memory available, UMFPACK will return a null pointer. You can ...


4

You could try using UPC++, which sets up a globally accessible address space distributed across your nodes. A more standard approach would be to learn how to use MPI.


4

Let me start with the following: Even very good and experienced programmers have a very hard time estimating whether a particular piece of code is performance critical or not. This has given rise to the adage "Premature optimization is the root of all evil", which can be translated as "Unnecessarily optimizing code leads to obscure code that will be ...


3

That might also have been a trick question. Let's say you want to solve the normal equations for $Ax=b$, i.e., $(A^T A) x = A^T b$. Let's assume for a moment that the questioner meant that $A$ is actually already stored in memory, so we know that that much memory is already available. Let's also assume that $A$ is tall and narrow (more specifically, has ...


3

You probably want a factorization of the form $\mathbf A = \mathbf L \mathbf D \mathbf L^T$, it can certainly be applied to a complex symmetric $\mathbf A$. LAPACK implements this factorization within [zsytrf] and provides a corresponding backsolution routine within [zsytrs]. There are sparse-direct versions of this algorithm too. PARDISO, TAUCS, and ...


3

Noting correctly that this is an array of pointers to arrays of ints: In the first version, the code is accessing element 0 of each successive array of ints, then element 1, etc. That means it is changing which pointer (to an array of ints) is being dereferenced on each inner iteration. In the second version, the code is accessing each element of the 0'th ...


3

From what I can see, on a normal machine, you would run out of memory if you are using a Dense matrix. The implementation you indicated in the question uses a dense matrix ( storing all the matrix entries ). Generally in most simulations, each element computes fluxes with its adjacent elements and the resulting matrix structure is sparse ( the number of non-...


3

I took a look at the sizes of the matrices you listed. I'm not sure I would characterize their size as "modest". ohne2 has only around 181K equations but nearly 7M non-zeros. Hamrle3 has around 1.5M equations which is a significant number for an unsymmetric, direct solver. If you are running a 32-bit executable on a 64-bit OS, you likely have a full 4GB ...


3

ADI is not a very good parallel algorithm. You should seriously consider formulating the problem in 3D and solving with multigrid. You could get a start with src/ksp/ksp/examples/tutorials/ex45.c. If you insist on using ADI, you should seriously consider allocating the three matrices separately. It's more memory, but then you won't have to reassemble on ...


3

Geoff mentioned the limitations of sparse direct solvers. To solve a 3D problem of size $N = n^3$, an optimal direct solver will need $N^{4/3}$ memory and $N^2$ time. See this question for more details. But it sounds like you have a structured grid diffusion problem. If the diffusion coefficient is smooth (or constant), you can solve the problem using ...


3

Presumably you are using numpy arrays for your data. The key is to avoid making a copy when you reshape them. For that, you want to ensure that you are just generating new "views" of your arrays. Take a look at the help page for numpy.reshape.


3

In essence, you are asking whether you can enumerate the integer lattice sites within your domain from $1$ to $N$ in such a way that accessing the east/west/north/south neighbors of a location $n$ requires accessing positions $n_e,n_w,n_n,n_s$ so that the distance between index $n$ and these four indices is minimal. This problem is equivalent to shuffling ...


2

For direct solution in 3D, you should probably be using some flavor of nested dissection (ND) or minimum degree (MD). These attack the storage requirements of A=LL' factorization directly, not the bandwidth (which has only an indirect effect on fill-in). On it's own, bandwidth reduction is just not strong enough to make 3D direct solve tractable. Good ND ...


2

It's impossible to tell without knowing the code you are using. But fortunately, segmentation faults are easy to debug: basically, a segfault means that you are accessing memory you should not access, and the operating system stops your program at the point where this is happening. This means, that if you run your program under a debugger, then you will see ...


2

Warning: this answer is only going to give a brief overview, for the real details, the one source that won't be wrong is the source code. The core matrix AIJ format is basically the same as the one known as compressed sparse row (CSR) or Yale format. This stores a sparse matrix as a list (normally concretely implemented as an array), $A$, of the nonzero ...


2

Unfortunately, GPUs will be of no help to you in this particular situation. Your problem is in the memory limitation; thus, you just do not have enough RAM resources to allocate/factorize/solve the system. GPUs usually have a much smaller amount of memory available on them, so they are not used to help with memory-limited problems. Suggestions to try: It ...


2

You can probably speed things up a little bit by storing the array in 4x4 square subarrays so that each of them fit in cache line (64 bytes = 4x4 32-bit integers). This changes the probability distribution of number of memory accesses from [0,0,1] (mean number of accesses = 3) (for 1,2,or 3 accesses) to [1/4,1/2,1/4] (mean number of accesses = 2). Maybe ...


2

Really the only way to run simulations on a multi-node cluster is to use MPI. And you don't distribute data because that would mean it would start out in one place. With MPI everything is distributed all the time, including being created distributed. It also means you don't create one global array: in your case you have 120 processes that each can store a ...


1

Among approximate techniques: Gradient descent should do a decent job, given the limitations. Randomized SVD is an effective technique. It is quick to implement, and there are ready-to-use error bounds (see e.g. https://doi.org/10.1137/090771806 for a review of that area). Among techniques that produce an exact solution: There is research on out-of-core ...


1

So called “matrix free” methods, relying primarily on the ability to perform multiplication of the matrix by a vector, lend themselves nicely to iterative techniques such as GMRES. The matrix itself might be on disk, but portions retrieved selectively to compute the necessary matrix-vector products.


1

The GPU consists of several streaming multiprocessors. Each SM has 64-96kB of shared memory that can be accessed by up to 1024-2048 threads. This shared memory allows these threads to communicate. To communicate between SMs you must write to and read from the GPU's global memory, which is 4-32GB in size. But it is best to think of your problem as consisting ...


1

A lot of this will depend on the details of your do_big_calculation function. In general you want to avoid pushing data to disk for performance reasons. Disk I/O speed is significantly slower than memory speed. There are some strategies that might help avoid creating that huge matrix in the first place. If the output value at each location depends only ...


1

It seems that you are going from one extreme to the other: you probably want to generate all $N$ particles at once without the for-loop; however, you don't want to generate all the num_steps at once since you only need two last ones. So, I think you are looking for something like: import numpy as np import numpy.random as npr N = 10000000 dim = 3 ...


1

You seem to have given the 1D equations for the discretizations, even though the problem is in 2D. Regardless, the explicit method requires the least memory since you don't even have to form a matrix to compute the solution at the next time step. If you have the solution vector at time $t_i$, you can apply simple stencil operations and directly obtain the ...


1

In my opinion there is not a memory type better than other. Simply they are different things. Normal ram are used only by the cpu, and the gpu ram is used onnly by the gpu. This is quite clear when you code direct in CUDA, i.e. if you want that the gpu use some data in normal ram you must move them to the gpu memory (host -> device). And if you want use ...


1

The scipy.integrate.odeint function does not take a vector as input. It takes a function f, an initial state y0 and a vector of points in time t (and possibly a bunch of other optional arguments to control the integration strategy). So you will need to provide it with a function that calculates the right-hand-side of your ODE system, which will boil down to ...


1

Beyond hwloc there are a few tools that can report on a HPC cluster's memory environment and which can be used to set a variety of NUMA configurations. I would recommend LIKWID as one such tool as it avoids a code based approach allowing you for instance to pin a process to a core. This approach of tooling to address machine specific memory configuration ...


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