7

I'm the primary author of many Julia libraries geared toward "architecture-specific optimizations", including LoopVectorization.jl and Octavian.jl. For BLAS-like operations, one of the most important optimizations LoopVectorization.jl does is "register tiling". While CPUs may have a huge number of actual registers (used for register ...


6

For a single (maybe a few) $b$'s, I have found Conjugate Gradient can beat an $LL^T$ Cholesky factorization . To do this, I use MKL's Inspector-Executor sparse matrix-vector product mkl_sparse_d_mv with the following: a good re-ordering (COLAMD and METIS worked for me). This speeds up each matrix-vector product. To avoid a permutation operation each time, I ...


6

For problems this small, sparse direct solvers are faster than most iterative solvers even if you include the cost of factorization. As a result, I don't believe that you will be able to find a preconditioner for an iterative solver that can work in less than twice the cost of the forward-backward solve. You either have to pay the memory cost of a sparse ...


6

Domain decomposition was developed in the late 1990s and early 2000s because it allowed the re-use of sequential PDE solvers: You only have to write a wrapper around it that sends the computed solution to other processors, receives other processors' solutions, and uses these as boundary values for the next iteration. This works reasonably well for the small ...


5

If we assume that $D$ is nonsingular, then there is a relatively straightforward (and efficient) solution based on an $LU$ decomposition. If we write $$ \pmatrix{D & B \\ B^T & A} = \pmatrix{ L_{11} & \\ L_{21} & L_{22}} \pmatrix{U_{11} & U_{12} \\ & U_{22}} = \pmatrix{L_{11} U_{11} & L_{11} U_{12} \\ L_{21} U_{11} & L_{21} U_{...


4

Most people use (algebraic as well as geometric) multigrid as preconditioners these days. It's an empirical observation that that leads to faster convergence in terms of iterations, given that in a typical CG iteration, applying a multigrid preconditioner takes up the vast majority of the effort compared to all other operations. In other words, you get the ...


4

I'd recommend a BLAS1/2-style implementation of LU decomposition with partial pivoting, along with forward/backward triangular solves. If you only have single digit numbers of unknowns, it should be fine to run dense/O(N^3) algorithms even on structured systems (ie banded/sparse/?). All these algorithms can be written in terms of "axpy" operations, ...


3

This landmark paper by George proves that a nested dissection ordering of a regular, four-node element, finite element mesh produces minimum fill-in. Although it is straightforward to produce such an ordering by inspection, the graph algorithms that attempt to do this for a general sparse structure only approximate this ordering. Assuming you are obtaining ...


3

Abdullah, thank you for the plug for our materials. We have repackaged these as a Massive Open Online Course (MOOC) on edX titled "LAFF-On Programming for High Performance". It is free for auditors. For info, see http://ulaff.net. The premier open source implementation of the BLAS is the BLAS-like Library Instantiation Software (BLIS). There is ...


3

This is a response/answer to the discussion in the comments to the question. I have been looking for the sources too, and I couldn't find them either. I am wondering if I am misremembering. Let me maybe discuss it in my own words. First of all, we know that a complex matrix $A$ is positive definite, e.g. $$Re[x^{\ast}Ax]>0 \qquad \forall x\in\mathbb{C}^n,$...


3

One such case is if the sparse matrix is banded. For example, tridiagonal linear systems can be solved in linear time using Thomas' algorithm. For small bandwidths, you can find an algorithm of linear time cost. Note that as the bandwidth grows, the hidden coefficient grows too. The literature on the topic is active and there are many characterizations as ...


3

I encountered a similar problem in the past and I could find no simple solution either. One of the terms is a Kronecker product, another is a rank-1 modification, but the rest makes the problem more difficult. I don't think there is a closed-form solution; you could try using an iterative method, dropping some terms to get a preconditioner. But if someone ...


2

Most of this was already discussed in the comments, but I would like to elaborate and put a detailed answer. There are no elementary characteristics (definiteness, symmetry, bandwidth) which can tell you whether the underlying (mixed or not) FEM/FVM is stable to solve the continuous problem. You can not tell anything about that just by looking at those ...


2

We can do some transformations to your problem to show that it's easily solvable via a linear program: Given a matrix $M$ with non-negative real entries and a vector $v$ you wish to solve the problem: $$ \begin{align} \min_v \quad & \lVert Mv \rVert_\infty \\ s.t. \quad & v_i\ge0 \\ & \sum_i v_i = 1 \end{align} $$ Now, note that $\lVert Mv \...


2

As correctly pointed out in the comments, HSL_MA57 is just an interface extended wrapper version of MA57. It is a Fortran 95 encapsulation of its original (Fortran 77) and offers a more feature-rich interface than the original MA57 (in addition to the Fortran 95 syntax). This info can be found in the user documentation of HSL_MA57. Note also, that MA57 is ...


2

If you don't care too much about which $b$ you're working with (say just for linear algebra), then you can use the Method of Manufactured Solutions (MMS) in Linear Algebra much like we differential equations geeks would for our problems to start with a known exact solution, $x^*$, derive a $b$, and then apply your solver to see how decent an $x$ you can ...


2

TL;DR For an image of size $m\times n$ you can solve this problem in $O(nm(\log(n) + \log(m)))$. In fact, there is nothing to "solve", the solution can be written down analytically. The differencing operator is linear space invariant, i.e. a convolution, and you are asking how to do a "de-convolution". This can be done efficiently in the ...


2

You can use Krylov iterative solvers with preconditioners. Since you mentioned that your problem of interest is linear elasticity, you end up with symmetric positive definite matrices, assuming that you use standard techniques. For this case, you can use the Conjugate Gradient method with the ILU preconditioner. These solvers are available in many languages. ...


2

This question is as ill-posed as the matrices you may want to invert :-) Even the most complicated iterative methods, say GMRES to name an example, are not terribly difficult to implement and require only a couple of hundred lines of code. That's probably only a factor of 4 or 5 worse than the arguably simplest method, Richardson defect correction. But it is ...


2

I would just compute the Cholesky factorization and then solve in batches using it. This will get technical, though: you will need to call Lapack functions by hand, I am afraid (*potrf and *potrs), since Python doesn't help you here, so to use the exact same algorithm you may want to check how it is done in the source of linalg.solve and dposv.f (good luck ...


1

They only deal with the 2D case, but I can point you to the work of some of my colleagues on https://doi.org/10.1137/17M1157155 . Using the off-diagonal low-rank structure in a recursive fashion, they can reach quasi-linear cost for the 2D case.


1

There has been progress both in terms of the speed of computers (basically driven by Moore's law) and in the algorithms used to solve LP's and especially MILP's. Overall, improvements in algorithms have had at least as large an effect as improvements in hardware. How this will work out for your particular model is a different question. Especially for MILP'...


1

If each row is a tensor product $v_i \otimes v_i$, then any vector $x$ that is a column-stacked version of a skew-symmetric matrix is in the kernel of $A$. For any $x$, you have $$(v_i \otimes v_i)x = v_i X v_i^T,$$ where $x = \mathrm{vec}(X)$ (see: https://en.wikipedia.org/wiki/Kronecker_product). If $X = -X^T$ then $$v_i X v_i^T = -v_i X^T v_i^T = -v_i X ...


1

There are many ways to solve systems with sparse matrices, so there is no way to answer this definitively and exhaust all possibilities. I'll add Krylov methods as one answer though. Many Krylov methods achieve fast results under the right conditions. The conditions for Krylov methods to perform well have been asked here before. Every Krylov solver is a ...


1

For a linear PDE, like the Laplace equation, when you discretize it you should get a linear system. Since you're 1D, the Thomas algorithm should be able to solve the system, and it's executed by running over the system once; Thomas algorithm is a direct solver, not an iterative one. If I understand your question correctly, you're asking what happens if you ...


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