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 ...
- 18.5k
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 ...
- 4,794
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 ...
- 53.7k
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 ...
- 53.7k
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{...
- 3,063
11
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 ...
- 2,305
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 ...
- 8,592
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 ...
- 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 ...
- 53.7k
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 ...
- 2,229
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=...
- 254
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 ...
- 3,108
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 ...
- 5,954
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 ...
- 8,592
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 ...
- 10.6k
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 ...
- 53.7k
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.
...
- 1,380
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. ...
- 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 ...
- 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 ...
- 1,193
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 ...
- 791
6
votes
Accepted
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 ...
- 5,954
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 ...
- 8,592
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 ...
- 1,147
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:
...
- 12.1k
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$), ...
- 4,794
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 ...
- 9,813
6
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 ...
6
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 ...
- 561
Only top scored, non community-wiki answers of a minimum length are eligible