Skip to main content
20 votes

Rule of thumb for sparse vs dense matrix storage

For what it is worth, for random sparse matrices of size 10,000 by 10,000 vs. dense matrices of the same size, on my Xeon workstation using MATLAB and Intel MKL as the BLAS, the sparse matrix-vector ...
Brian Borchers's user avatar
16 votes
Accepted

Rule of thumb for sparse vs dense matrix storage

All matrix operations are memory bound (and not compute bound) on today's processors. So basically, you have to ask which format stores fewer bytes. This is easy to compute: For a full matrix, you ...
Wolfgang Bangerth's user avatar
13 votes
Accepted

Why does sparse linear algebra have a low arithmetic intensity?

BLAS1-operations, BLAS2-operations, and sparse-operations share the same curse of low arithmetic intensity, that they perform $O(1)$ flops for each memory read (contrast this to a BLAS3-operation like ...
rchilton1980's user avatar
  • 4,862
12 votes
Accepted

Compute all eigenvalues of a very big and very sparse adjacency matrix

You can use the shift-invert spectral transform [1] and compute the spectrum band by band. The technique is also explained in my article [2]. Besides the implementation in [1], an implementation is ...
BrunoLevy's user avatar
  • 2,315
12 votes
Accepted

Solving linear system of the form $ABx=b$

Defining the auxiliary variable $y=Bx$ yields the following algebraically equivalent expanded system, $$\underbrace{\begin{bmatrix} 0 & A \\ B & -I \end{bmatrix}}_{K} \underbrace{\begin{...
Nick Alger's user avatar
  • 3,143
11 votes

Which C++ linear algebra library is probably the fastest on solving huge sparse [square matrix] linear system?

Eigen 3 is a nice C++ template library some of whose routines are parallelized. c.f. Eigen documentation The parallelization is OMP only, so if you intend to parallelise using MPI (and OMP) it is ...
Nox's user avatar
  • 351
10 votes
Accepted

How can a CG solver solve a non positive definite sparse matrix

I highly recommend the following read: J.R. Shewchuk, "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain" In short, if the matrix is non-positive definite, there is no ...
Anton Menshov's user avatar
  • 8,672
10 votes
Accepted

Is there an iterative solver for dense matrices with possible zero diagonal entries?

Iterative Krylov-subspace solvers generally only require matrix-vector products and don't care whether or where there are zeros in the matrix. In your case, unless you have other information about the ...
Wolfgang Bangerth's user avatar
9 votes
Accepted

C standard for computational science

In theory, as the original authors, you're free to pick and name a standard, then expect others to follow it. In practise, if you're supporting an HPC system, then your choice is likely to be ...
origimbo's user avatar
  • 2,249
9 votes
Accepted

Efficiently computing $e^{tX}$ for many different values of $t$

An anti-Hermitian matrix is diagonalizable, with orthogonal eigenvectors (ref). Hence you can write $X = PDP^{-1}$, where $D$ is a diagonal matrix. Therefore the exponential can be calculated as $e^X=...
rpm2718's user avatar
  • 254
9 votes

How large is large for direct solvers?

The answer here mostly depends on how good the preconditioner you have for the iterative solver. If you don't have a good preconditioner, direct methods tend to be the best until you run out of ram. ...
Oscar Smith's user avatar
8 votes
Accepted

Solving Ax = b with sparse A and sparse b

When looking at the solution of your system, you will find that almost all entries of $x$ are nonzero although the right-hand side is "sparse". Hence, whatever algorithm you use, it'll have to visit ...
Nico Schlömer's user avatar
8 votes
Accepted

Eigen - Max and minimum eigenvalues of a sparse matrix

There are two relatively convenient options for calculating selected (e.g. a few largest or smallest) eigenvalues using Eigen. The first is Spectra, a header-only C++ library based on Eigen that uses ...
Bill Greene's user avatar
  • 6,064
8 votes
Accepted

Integer operations vs floating point operations

There is a nice discussion on StackOverflow regarding floating point vs integer operations. In short, the performance of the operations depends a lot on processor architecture how the data is stored ...
Anton Menshov's user avatar
  • 8,672
8 votes
Accepted

How to compute all the eigenvalues of a large sparse matrix using matlab?

"Get more RAM" may be one of your best options. :) Prices are reasonably low right now, and it's one of the best upgrades you can gift your computer anyway. 10k x 10k is borderline but still doable ...
Federico Poloni's user avatar
8 votes
Accepted

Fast way to build stiffness directly as CSC matrix

COO is an unsuitable matrix format except for particular purposes (e.g., if there is a substantial number of rows that have no entries at all, possibly with the exception of the diagonal). The way ...
Wolfgang Bangerth's user avatar
7 votes
Accepted

Applying the result of Cuthill-McKee in SciPy

The reverse Cuthill-McKee algorithm produces a reordering that applies to both the rows and columns. This is because it works by considering matrices as graphs of (undirected) connected nodes. ...
Nick C.'s user avatar
  • 188
7 votes

Solve for $C$ such that $C^{T}AC$ is banded of given width

Yes. The block Lanczos algorithm http://www.netlib.org/utk/people/JackDongarra/etemplates/node250.html produces a block triangular matrix where you control the block size, hence the bandwidth. ...
Carl Christian's user avatar
7 votes
Accepted

bit-packing and compression of data structures in scientific computing

Reducing Memory for Sparse Matrices One method (that they mention in the first paper you linked, but is worth emphasizing) is the Block Compressed Sparse Row (BCSR) storage format. If your problem ...
Tyler Olsen's user avatar
  • 1,522
7 votes

Rule of thumb for sparse vs dense matrix storage

Even if a matrix is very sparse, its matrix product with itself can be dense. Take for example a diagonal matrix and fill its first row and column with nonzero entries; its product with itself will be ...
Henrik Schumacher's user avatar
7 votes

C standard for computational science

You should definitely jump to C99, or newer(!). The C99 standard introduced the restrict keyword. Loosely speaking, with this keyword you can inform the compiler ...
wim's user avatar
  • 571
7 votes
Accepted

Sparse matrix inversion

For a matrix that small, you're probably not going to do better than using dense methods. I wrote up a quick test in C++ for an 18x18 matrix with your sparse structure and randomly generated values ...
LedHead's user avatar
  • 1,253
7 votes
Accepted

What is the format of saving sparse matrix in MATLAB?

Matlab internally uses compressed sparse column (CSC) format for sparse matrices. The design and implementation of Matlab's sparse matrices are described in this document. As a consequence of using ...
Will P.'s user avatar
  • 821
7 votes

How large is large for direct solvers?

Here is a "real live" test that I did: A few weeks ago I did a simulation of a large grid of 1Ohm resistors. Both with an iterative relaxation solver. ( The code is gitlab: grid-of-1ohm-...
mond's user avatar
  • 71
6 votes
Accepted

Fastest way to solve a sparse unsymmetric system many times

Generally speaking, for many right-hand side (RHS) problems, a direct solver is a more feasible solution for several reasons: Major computations are performed during the factorization step (which is ...
Anton Menshov's user avatar
  • 8,672
6 votes
Accepted

Fast c++ library to solve very big sparse systems

I second the idea of using Eigen, which is pretty efficient, but also very simple to include. If you need a lot more performance, you could try to use PETSc or Trilinos. They are very powerful ...
BlaB's user avatar
  • 1,157
6 votes
Accepted

Why is the speed of the parts of the LU-decomposition so different?

First, don't forget to also time the LU decomposition in a loop! Otherwise it's not really a fair comparison. If I do that, I get the following timings: ...
Christian Clason's user avatar
6 votes
Accepted

Correct use of scipy's sparse.linalg.spilu

If $\mathbf L$ and $\mathbf U$ give an approximate factorization of $\mathbf A$, you wouldn't want to use $\mathbf P = \mathbf L\cdot \mathbf U$ as a preconditioner (that's approximately $\mathbf A$), ...
rchilton1980's user avatar
  • 4,862
6 votes
Accepted

Ways to solve $Ax=b$ for a sparse (banded) $A$ with updates

If the only non-zero entries of $A_{ij}$ have $j$ in $\{i - 1, i, i + 1\}$, then $A$ is a banded matrix with bandwidth 1. More generally, you can talk about matrices of bandwidth $k$ where $k$ is any ...
Daniel Shapero's user avatar
6 votes
Accepted

(FEM) Nodes reordering for sparse matrix storing techniques

You should use a reordering. Although it's true that storing a sparse matrix requires the same amount of memory whether or not you reorder it using RCM, reordering it should lead to faster ...
rchilton1980's user avatar
  • 4,862

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