Hot answers tagged

24

I'm going to disagree with some of the other answers and say that I believe that figuring out how to use LAPACK is important in the field of scientific computing. However, there is a large learning curve to using LAPACK. This is because it is written at a very low level. The disadvantage of that is that it seems very cryptic, and not pleasant to the senses. ...


10

This can be interpreted as summing over an index of a tensor when the vector $x$ is reshaped into a box of numbers instead of a list. In particular, if $X$ is the $d\text{-by-}d$ folded version of $x$, then the operation you are doing is, \begin{align} Dx &= \mathrm{vec}\left((I \otimes \mathbf{1})\mathrm{vec}(X)\right) \\ &= \mathrm{vec}(\mathbf{1}^...


9

There's no reason to append a row of 1's. You should just perform a rank-revealing QR factorization (like with routine SGEQP3) on $A^T$, and the last column of $Q$ should be in the nullspace. This has the added advantage that the relative magnitude of the last element on the diagonal of $R$ gives you some idea of how singular the solution is. Even better ...


8

The problem, of course, is that computing the true rank (e.g., via a QR decomposition) is not really any cheaper than computing a low-rank representation of the matrix. The best you can probably do is to use a randomized algorithm to find low-rank approximations. These can, at least in theory, be significantly faster than working on the entire matrix ...


8

There is a neat trick I have recently learned from this paper. You start doing rank-revealing QR, and stop after the first $k$ Householder reflections, when you have a matrix of the form $$ \begin{bmatrix} R_1 & R_{12}\\ 0 & R_{22} \end{bmatrix}, $$ with $R_1$ triangular of size $k\times k$, and $R_{22}$ typically not triangular (since we stopped ...


8

When one says an algorithm is of order $O(n)$, that may mean that the complexity is given by: $c + b*n$. With every new element you add you increase in runtime (effectively). What mathematically minded people often forget is that these statements do not include how large the constants are. That of course carries over to $O(n²)$ and such. I can not answer ...


8

The matrix B (M in the documentation) needs to positive definite according to the documentation: "If sigma is None, M is positive definite", this is in addition to the first requirement "M must represent a real, symmetric matrix if A is real" which your B follows. The eigenvalues of your current matrix B are -1, 1 and 6. So matrix B is ...


7

Took me quite a while to figure this out and as usual it becomes obvious after you find the culprit. After checking the problematic cases reported in David S. Watkins. A case where balancing is harmful. Electron. Trans. Numer. Anal, 23:1–4, 2006 and also the discussion here (both being cited in arXiv:1401.5766v1), it turns out that matlab uses the ...


7

The cost of the multiplication is almost insignificant compared to the cost of loading the data from memory (and writing it back). If you're worried about performance, you should be thinking about data locality. Perform more flops with each of your $x_i$ and $z_i$ values (if that's possible) when you load them from memory and before proceeding on. Setting ...


7

I usually resist telling people what I think they should do rather than answering their question but in this case I'm going to make an exception. Lapack is written in FORTRAN and the API is very FORTRAN-like. There is a C API to Lapack that makes the interface slightly less painful but it will never be a pleasant experience to use Lapack from C++. ...


7

Here's another answer in the same vein as the above. You should look into the Armadillo C++ linear algebra library. Pros: The function syntax is high-level (similar to that of MATLAB). So no DGESV mumbo-jumbo, just X = solve( A, B ) (although there is a reason behind those oddly-looking LAPACK function names...). Implements various matrix decompositions ...


7

It's the blocked variant of Householder-QR that is driving this design. If you look in Golub and Van Loan's book (Ch 5.2 or so) they talk about how k-iterations of the algorithm can be blocked together by accumulating the individual reflectors into a rank-k reflector of the form $\mathbf I + \mathbf W \mathbf Y^{\mathrm T}$, where both $\mathbf W$ and $\...


7

I think you are overestimating the overhead of computing L. There are zero extra operations needed; the only additional cost is writing to RAM some numbers that you have already computed anyway. The algorithms commonly used (in Lapack, for instance) to compute U also compute L along the way, and you'd save 0 flops by omitting it. For instance, if you think ...


6

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 libraries to store and solve sparse systems, they allow for a large number of iterative or direct solvers and are compatible with MPI for added performance. However, ...


6

Seems that you have a duplicate eigenvalue. Thus, you have two eigenpairs $(\lambda_1, x_1)$ and $(\lambda_2, x_2)$ where $\lambda_1 = \lambda_2$. Denote $\lambda = \lambda_1 = \lambda_2$. Let $\alpha$ and $\beta$ be arbitrary complex numbers. Then $$A (\alpha x_1 + \beta x_2) = \alpha A x_1 + \beta A x_2 = \alpha \lambda x_1 + \beta \lambda x_2 = \lambda (\...


6

LAPACK has been on the cutting edge for just about three decades, and probably still is for its niche. However, given given recent developments in libraries for the simpler BLAS-type matrix operations that LAPACK traditionally builds upon, it is perhaps conceivable that we could see the emergence of serious competitors to the traditional FORTRAN-based LAPACK ...


6

The trick is trying to find out why that matrix has real eigenvalues in the first place. Usually it is because a suitable set of conjugations turns it into a symmetric matrix, and then you can reduce to a symmetric computation. Multiplying and dividing by $(AC)^{-1}$ you can rewrite $$ D_1 = C^{-1}BAC (C^{-1}BAC+I)^{-1}, $$ so your computation is equivalent ...


5

Yes, you want to call the BLAS routine DGEMM. The place to start for how to call it from C is to look at the documentation for DGEMM, which you can find online. Then you want to understand how to call FORTRAN routines from C (DGEMM, like all the standard BLAS routines has a FORTRAN calling convention). For example, this document https://computing.llnl.gov/...


5

The LAPACK routine zgesv first computes the LU factorization, and then solves the system making use of the factorization. It is a simple driver for calling the two routines zgetrf (compute the LU factorization) and zgetrs (solve the system). The LU factorization is computed using partial pivoting and row interchanges. Then, the system is solved using simple ...


5

tl;dr Yes. But your question doesn't make it clear that you understand what LAPACK is about. LAPACK is both a software as well as an interface. That is, the operations that LAPACK defines are standard enough that they can be replaced by other software, such as ATLAS, MKL, and so on. Another way of looking at this is that any software that does linear algebra ...


5

The other answers already tell you what went wrong, but I will add a terminology note: the term for what is happening is that the pencil $A - \lambda B$ is a singular matrix pencil, i.e., $\det (A - \lambda B)$ is identically equal to zero. So there are no generalized eigenvalues (or, at least, they cannot be defined as usual as the roots of the generalized ...


4

First of all, yes, these are all based on an initial tridiagonalization (often quoted to be $\frac{4}{3}n^3$ flops). DSYEV is just an easier to use version of DSYEVX, so let's ignore it for now. The basic breakdown is such: DSYEVX: Tridiagonal implicitly shifted QR (bulge chasing) DSYEVD: Divide and Conquer (stitching by rank 1 modifications) DSYEVR: Use ...


4

The MATLAB syntax qz(A,B,'real') is consistent with schur(A,'real'), so we might as well ask why the default is complex in the Schur form. Two reasons spring to mind. Backward compatibility. Probably there was a time when only the complex Schur form was implemented in Matlab (possibly from the pre-LAPACK times), and the default is retaining that behavior, ...


4

When you call Lapack's zgtsv, it doesn't just solve a tridiagonal system $Ax=b$. What it does first is perform an LU factorization (zgttrf) $A = LU$, where $L,U$ are lower- and upper-tridiagonal matrices, and only then proceeds to solve $LUx=b$. When you give it the lower, main, and upper diagonals of the matrix $A$, those diagonals are overwritten by the ...


4

Note: I haven't run your code. Perhaps this is a problem with the row-major/column-major conventions at play here: http://docs.nvidia.com/cuda/cublas/index.html#data-layout. You seem to use row-major matrices, but BLAS uses the column-major convention. When you pass row-major $A$ and $B$ to dgemm, it implicitly interprets them as $A^t$ and $B^t$, giving ...


4

In my experience the answer to this question is not clear-cut; it is dependent on the span of the eigenvalues and relative matrix structure itself. That said, your current approach evokes an implictly shifted QR solver that is essentially the standard for this exact type of problem. With some experimentation, however, you may squeak out a modest performance ...


4

Supposing you already have an LU factorization, you can save a half of a forward substitution step. In the system $Lx=b$, you would have $$ \begin{pmatrix} l_{11}&0&&&&\cdots&0\\ \vdots&&&&&&0\\ l_{k-1,1}&\cdots&l_{k-1,k-1}&0&\cdots&\cdots&0\\ l_{k,1}&\cdots&l_{k,k-1}&l_{k,k}&...


4

Another approach, which might be of interest to you is randomized sampling. This is of particular interest if you can quickly compute matrix-vector products $x\rightarrow Ax$ and $x\rightarrow A^* x$. The core idea is to form a small sampling matrix $S = A\Omega$, where $\Omega$ is a Gaussian random matrix. If the sampling matrix is large enough, $S$ will ...


4

You write, that you are computing the eigenvalues of a symmetric matrix. Does the matrix have real entries? In this case all eigenvalues are real, and you can use a symmetric eigenvalue solver, which returns only real entries. Hence, sorting them should not be a problem. When your matrix has complex entries, you have to track the eigenvalues. I am assuming ...


4

Well, other than the usual "don't invert your matrices unless you need the inverse itself" you can still use the banded routines ?gbtrf and then use ?gbtrs with the right hand side being an identity matrix.


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