18 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 ...
user avatar
14 votes

20% performance penalty for a nice software design

It's a question on what you spend your time on. For most of us, we spend 3/4 of the time programming and 1/4 of the time waiting for results. (Your numbers may vary, but I think the number is not ...
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 ...
user avatar
  • 4,316
13 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 ...
user avatar
12 votes
Accepted

Is the Thomas algorithm the fastest way to solve a symmetric diagonally dominant sparse tridiagonal linear system

I believe comparing an iterative method (multigrid) to a direct/exact method (Thomas) in terms of exact operation count isn't really meaningful. IIRC, Thomas operation count is $8N$ for any ...
user avatar
  • 2,203
12 votes
Accepted

How to efficiently implement Dirichlet boundary conditions in global sparse finite element stiffnes matrices

In deal.II (http://www.dealii.org -- disclaimer: I'm one of the principal authors of that library), we do eliminate whole rows and columns, and it is not too expensive overall. The trick is to use the ...
user avatar
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{...
user avatar
  • 3,013
11 votes

Is the Thomas algorithm the fastest way to solve a symmetric diagonally dominant sparse tridiagonal linear system

The short answer is that the Thomas algorithm will be faster than any iterative scheme for almost all cases. The exception would perhaps be applying a single iteration of a very simple iterative ...
user avatar
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 ...
user avatar
  • 8,382
10 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 ...
user avatar
  • 341
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 ...
user avatar
9 votes
Accepted

20% performance penalty for a nice software design

Very few scientific software developers understand good principles of design, so I apologize if this answer is a bit long-winded. From a software engineering perspective, the goal of the scientific ...
user avatar
  • 6,841
9 votes

Solving Lx = b for big sparse Laplacian matrices

Daniel Shapero's answer is excellent, but I felt I should add the following: Properly preconditioned iterative methods will almost certainly win here, unless your systems have a very, very special ...
user avatar
  • 372
9 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 ...
user avatar
  • 2,185
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=...
user avatar
  • 254
8 votes
Accepted

Eigenvalues of a sparse banded nonsymmetric matrix from an elliptic operator

The simple answer is that you would use inverse iteration (subspace or with deflation). This is basically the power method (repeatedly multiplying the matrix by a vector and normalizing, singling out ...
user avatar
  • 4,410
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 ...
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 ...
user avatar
  • 5,744
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 ...
user avatar
  • 8,382
8 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 ...
user avatar
  • 2,199
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 ...
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 ...
user avatar
7 votes
Accepted

Solving Lx = b for big sparse Laplacian matrices

Which one was faster when you tried using Eigen? You may also consider Trilinos if you're using C++. By "generic Laplacian matrix", do you mean a finite difference / finite element discretization of ...
user avatar
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. ...
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. ...
user avatar
  • 188
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 ...
user avatar
  • 1,522
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 ...
user avatar
  • 1,133
6 votes
Accepted

What is the cost of factorization for one-dimensional sparse problems?

It may help to define $N$, the number of discretization points along a 1D edge, and relate it to $n$, the number of unknowns in the system. In 2D on a square grid of points, $n = O(N^2)$. Nested ...
user avatar
  • 2,961
6 votes
Accepted

Eigen - store sparse matrix as binary

I've rolled my own. Here is a MCVE: ...
user avatar
6 votes

How to compute the rank of a large sparse matrix in MATLAB

There are two things you can do that may significantly reduce the computation time and memory used. First, you aren't using the Q matrix, so don't ask MATLAB to compute it. Since it is dense, that ...
user avatar
  • 5,744

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