59 votes
Accepted

A geometric multigrid solver for Mathematica?

Background Details about multigrid solvers can be found in this pretty neat script by Volker John. That's basically the source from which I drew the information to implement the V-cycle solver below. ...
Henrik Schumacher's user avatar
37 votes
Accepted

Complex eigenvalues from a sparse Hermitian matrix

Very good observation. Indeed, this issue is really frustrating. To single out the issue: It seems that Arnoldi's method is to blame: ...
Henrik Schumacher's user avatar
22 votes
Accepted

Efficiently defining a SparseArray function

The way the code is written, you can neither exploit packed arrays nor any vectorization. There are two major reasons: Using Rule prevents using packed arrays ...
Henrik Schumacher's user avatar
21 votes
Accepted

Implementation of MATLAB's "numgrid" function

The source code of some of the MATLAB functions is available to its users. In this example, if you type edit numgrid, its file will be open, and to my surprise, it ...
Ben Izd's user avatar
  • 9,239
20 votes

A geometric multigrid solver for Mathematica?

3D Example The problem with direct solvers is that starting in 3 dimensions, their performance for dealing with matrices stemming from PDEs drops rapidly. This is why I wanted to show at least one 3-...
Henrik Schumacher's user avatar
19 votes
Accepted

Solving "Resistance between two nodes on a grid" problem in Mathematica

In addition to Carl Woll's post: Computing the pseudoinverse of a the graph Laplacian matrix (a.k.a. the KirchhoffMatrix) is very expensive and in general leads to ...
Henrik Schumacher's user avatar
18 votes
Accepted

Quickly creating a sparse array

You can use the first documented usage for SparseArray: So what you want to do is collect all of the rules {i,j}->val before ...
Jason B.'s user avatar
  • 68.4k
16 votes
Accepted

Most efficient way to build block / banded SparseArray

This is definitely faster. It took me some time to figure out the combinatorics and there might still be some potential for improvement within the compiled functions. ...
Henrik Schumacher's user avatar
16 votes

Solving "Resistance between two nodes on a grid" problem in Mathematica

Based on rcampion2012's answer to Efficient Implementation of Resistance Distance for graphs?, you could use: ...
Carl Woll's user avatar
  • 131k
15 votes

Optimising subdiagonal shift matrix generation time

...
kglr's user avatar
  • 395k
14 votes
Accepted

How to convert a sparse matrix to a list of sparse arrays?

List @@ s will do the trick: ...
WReach's user avatar
  • 68.9k
14 votes
Accepted

Why is IdentityMatrix not defined as a SparseArray by default?

The fastest way to get the identity matrix as a sparse array is simply this: IdentityMatrix[10000, SparseArray]; // AbsoluteTiming Here are just some thoughts on ...
Jens's user avatar
  • 97.3k
14 votes

Using dense matrix is faster than sparse matrix in LinearSolve

You cannot benefit from SparseArray if your matrix is not sparse! The data structure storing a SparseArray needs a certain ...
Henrik Schumacher's user avatar
14 votes
Accepted

Set every element to zero in matrix unless it's 1 or 1/2

Be careful with assigning a "form" (e.g. MatrixForm) to a variable. Another observation is that a is a ...
gwr's user avatar
  • 13.5k
14 votes
Accepted

Reusing PARDISO symbolic factorization

No, this is not built-in and I have been waiting quite long for such a feature. PardisoLink With the help of Szabolcs's package "LTemplate`" I had written such an interface for the Intel MKL ...
Henrik Schumacher's user avatar
13 votes
Accepted

What are valid settings for the TreatRepeatedEntries sparse array option?

The complete list is (V11.2) First or symbolic 0 Total or symbolic 1 List or symbolic 2 ...
user21's user avatar
  • 39.8k
13 votes

Most efficient way to build block / banded SparseArray

Not faster, but maybe easier to implement ...
Henrik Schumacher's user avatar
12 votes

Lowest Magnitude Eigenvalues of Large Sparse Matrices

Use the Arnoldi method with shift-inversion: Eigenvalues[A, 3, Method -> {"Arnoldi", "Criteria" -> "Magnitude", "Shift" -> 0}] gives you the three ...
Roman's user avatar
  • 47.5k
12 votes
Accepted

Efficiently populate a Sparsearray for a set of rules for a constrained basis

There are a couple of things that can be improved: Binary search can simply be replaced by a lookup table that is very cheap to build. Moreover, IntegerPart[i/len] ...
Henrik Schumacher's user avatar
12 votes

Implementation of MATLAB's "numgrid" function

After searching for a while I found the source code of numgrid here. So, here's my trial with FunctionCompile: ...
xzczd's user avatar
  • 66.2k
12 votes

Implementation of MATLAB's "numgrid" function

...
chyanog's user avatar
  • 15.6k
12 votes
Accepted

Fast Sparse Tensor Addition

Mathematica uses a modified CSR format for SparseArrays of dimension two or higher. But the CSR format was actually only developed for matrices (dimension equal to ...
Henrik Schumacher's user avatar
11 votes
Accepted

Constructing a large tridiagonal matrix with alternating signs

Here is a way using Band: ...
Jens's user avatar
  • 97.3k
11 votes
Accepted

Wrong eigenvalues from a sparse matrix

It appears that the proximity of the eigenvalues causes Eigensystem with the default parameters to be inaccurate. This can be fixed by increasing the basis size to 30 ...
Omrie's user avatar
  • 361
11 votes

Set every element to zero in matrix unless it's 1 or 1/2

A = KroneckerProduct[PotentialTilde[1], PotentialTilde[1]]; B = A (1 - Unitize[A - 1/2] Unitize[A - 1]); B // MatrixForm It is even much more efficient to use ...
Henrik Schumacher's user avatar
11 votes

Quickly creating a sparse array

Let's assume we are given the following data: ...
Henrik Schumacher's user avatar
11 votes

Optimising subdiagonal shift matrix generation time

tminus[n_Integer /; n >= 2] := SparseArray[Band[{2, 1}] -> 1, {n, n}]
Roman's user avatar
  • 47.5k
11 votes
Accepted

Optimising subdiagonal shift matrix generation time

Since your matrix is large and has very few nonzero elements it is much more efficient to create and store it as a SparseArray. In addition many matrix operations ...
Ray Shadow's user avatar
  • 7,826
11 votes

Generate coefficient array from general formula of linear equation system

Partial ND It's possible to again extract all this weight data and stuff in the higher dimensional cases using effectively the same procedure as for the 1D: ...
b3m2a1's user avatar
  • 46.9k
11 votes
Accepted

Speeding up tensor contractions and multiplication

You seem to be coming from Matlab as you try to transpose a vector, a concept that is not that useful in Mathematica. We will see in a second why that is. Dense tensor example ...
Henrik Schumacher's user avatar

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