14

There has been considerable recent interest in numerical linear algebra using mixed precision with some combination of 16, 32, 64, and 128 bit floating point arithmetic. For example, a low precision factorization of a matrix can be used used to precondition a higher precision iterative solution. Since the factorization takes $O(N^{3})$ operations and the ...


4

SYRK is not really relevant here I think; it is just something else that happens to have the same name "rank k update". In your case, you need to know how to update a QR factorization by inserting rows; a good reference is Golub, Van Loan, section 6.5.3: Appending or Deleting a Row. Many computational environments have it already implemented for you, see e....


3

As the commenters have already said, for your specific problem, QR factorization might be a better approach not because of speed but because of numerical stability. It's still a good question to ask in general. One of the advantages you cite is that $LDL^*$ can be used for indefinite matrices, which is definitely a point in its favor. The linear algebra ...


2

If you will be adding lots of rows, then you will want to use the tall skinny QR algorithm (TSQR) of Demmel et al, 2008, https://arxiv.org/abs/0806.2159 This algorithm can be combined with the level 3 BLAS QR algorithm of Elmroth and Gustavson, 2000 to efficiently update the factorization with each new block of rows. It is best to save the new rows into a ...


2

You generally don't gain much for sparse matrix-matrix and sparse matrix-vector products using things such as SSE/AVX/... if the matrices are large. That's because these instructions offer the ability to do some floating point operations in parallel -- but for large sparse matrices, you are limited by the time it takes to get data from memory onto the ...


2

Yes, those Givens rotations do not seem correctly implemented to me. Since you are doing this as a learning project: learn from the masters http://www.netlib.org/lapack/explore-html/de/da4/group__double__blas__level1_ga54d516b6e0497df179c9831f2554b1f9.html. Also make sure that the givens rotations you are applying from the left are the same as the ones you ...


1

The example you cited appears to be generating random Householder vectors and multiplying them out using backwards accumulation. Another simple thing to do would be to generate a random matrix $\mathbf A$, then compute its $\mathbf A=\mathbf Q \mathbf R$ decomposition and discard the $\mathbf R$ factor. The two LAPACK functions that you need are [geqrf] (to ...


1

Not a full answer, but I will just mention two subtle details: OP is applying LDLT to matrices that are positive semidefinite in exact arithmetic; hence one would expect that, barring catastrophic cancellation errors, LDLT always uses 1x1 pivots rather than 2x2 pivots. Hence the benchmark results for those matrices may differ from generic ones (based on ...


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