18

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 multiply was faster for densities of 15% or less. At 67% (as proposed by another answer), the dense matrix-vector multiplication was about three time faster.


14

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 completely without merits.) So, if you have a design that allowed you to program twice as fast (3/4 of a time unit instead of 1.5 time units), then you can can take ...


13

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 gemm, which performs $O(N^3)$ flops over $O(N^2)$ reads and only becomes more and more arithmetic-intensive/compute-bound for large $N$). However, sparse ...


12

The LU factors of a sparse matrix are at least somewhat sparse. The $Q$ matrix in QR can also somewhat preserve sparsity, and is typically used when the matrix is very long and skinny. The SVD of a sparse matrix will almost always have fully dense $U$ and $V$ factors, so it destroys any reason to perform the computations treating the matrix sparsely.


12

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 tridiagonal system. The only time I can imagine multigrid conceivably beating that is for a trivial case of having a linear solution, and even then the cost of evaluating ...


12

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{bmatrix} x \\ y \end{bmatrix}}_{u} = \underbrace{\begin{bmatrix} b \\ 0 \end{bmatrix}}_{f},$$ which you could solve with GMRES or another nonsymmetric Krylov method. ...


12

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 store 8 bytes (one double) per entry For a sparse matrix, you store 12 bytes per entry (one double for the value, and one integer for the column index of the ...


11

For eigenvalues, simply take $k$ largest or smallest eigenvalues of $T$. They are good approximations of $A$, provided that the number of Lanczos iterations is large compared to $k$. Things are a little trickier if we want eigenvectors as well. The simplest way is to multiply each eigenvector $\mathbf{u}_i$ of $T$ by $V$ to the left, where $V$ is, as you ...


11

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 scheme such as Gauss-Seidel, but this is highly unlikely to give an acceptable solution. Also, this is ignoring parallel processing concerns. Multigrid is an ...


11

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 fact that the sparsity pattern is typically symmetric, so you know which rows you need to look into when eliminating a whole column. The better approach, in ...


10

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 matrix (e.g., symmetry), you could for example use GMRES. What you probably had in mind is the question of preconditioning, and that you can't use things such ...


9

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 software developer is to design a solution that satisfies a set of constraints that are often conflicting. Here are some typical examples of these constraints, ...


9

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 available in C++ in my Graphite software [3] (update Jan 17: now everything is ported to geogram/graphite version 3), that I used to compute the eigenfunctions ...


9

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 probably not suitable for your purpose. The nice feature of Eigen is that you can swap in a high performance BLAS library (like MKL or OpenBLAS) for some routines ...


9

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=Pe^DP^{-1}$, and $e^{tX} = Pe^{tD}P^{-1}$. If $d_1$, $d_2$, etc.... are the diagonal elements of $D$, then for each value of $t_i \in t$, $e^{t_iD}$ is just a ...


8

The matrix product $B = AA^T$ is generally faster to apply as $A (A^T x)$ even when the product matrix $B$ is already available. Only for peculiar graphs with extremely low expansion factor would the product matrix be faster. PDE graphs in two or more dimensions have higher expansion factors. For example, if $A$ is a 9-point differencing operation on a 2D ...


8

Yes, the most common approach is to rebuild. Data structures that are modifiable in-place tend to be less efficient once set up, and reallocation is actually quite cheap compared to reassembly (e.g., due to nonlinearity) so it's really a fine solution. Outside of relatively rare niches with very easy solves, attempts to use dynamic data structures in the ...


8

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 the eigenvector corresponding to the eigenvalue of largest magnitude) applied to $A^{-1}$. Since you desire the $k$ closest to the origin, you need to use some ...


8

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 each and every entry at least once, so one wouldn't expect that you can save a lot of time using the sparsity of $b$. Right-hand sides where you can save a lot ...


8

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 algorithms similar to ARPACK (implicitly-restarted Arnoldi) to calculate a few eigensolutions. Since it is header-only, you simply download and include the ...


8

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 in memory and in which order it is accessed if (and which) SSE/AVX/AVX2/etc instructions are used (and how efficiently) This probably provides some insight ...


8

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 restricted among the compiler standards that the given system you're using for testing has for its tool stack (you are testing your software, aren't you?). This might ...


8

"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 on modern computers: that matrix takes $10^4 \times 10^4 \times 8$ bytes, that is, 760 MiB. On my laptop that code runs without problems. Another option is ...


8

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 all large-scale codes I know of build the system matrix is through a three-step process: In the first step, you loop over all cells and figure out which entries ...


7

The simplest/fastest way to solve ill-conditioned problems is to increase precision of computations (by brute force). Another (yet not always possible) way is to re-formulate your problem. You might need to use quadruple precision (34 decimal digits). Even though 20 digits will be lost in a course (because of condition number) you will still get 14 correct ...


7

You may want to watch lecture 34 here: http://www.math.tamu.edu/~bangerth/videos.html


7

Your matrix is not "banded" in the sense of "banded" direct solvers, so don't bother with those. Multigrid is absolutely the best way to solve these "pressure Poisson" problems. There are lots of libraries and it's not difficult to implement a simple multigrid algorithm for structured grids. CG and direct solvers are fundamentally non-scalable, so the ...


7

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. Certainly, one can argue that a block tridiagonal matrix is not a "proper" banded matrix as there regular patches of certain zeros within the band. If you want ...


7

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. According to the function's documentation in SciPy, the output array is the permuted row/column indices, so you can simply do the following perm = ...


7

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 creates dense $n \times n$ blocks (common in e.g. FEM with multiple DoFs per node), you modify the typical CSR (CSC) storage scheme to store only a single column (...


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