Hot answers tagged

24 votes
Accepted

How to start using LAPACK in c++?

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 ...
user avatar
  • 1,133
10 votes
Accepted

Fast vector - "diagonal" matrix multiplication

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$,...
user avatar
  • 3,013
8 votes

Rapidly determining whether or not a dense matrix is of low rank

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 ...
user avatar
8 votes

Rapidly determining whether or not a dense matrix is of low rank

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{...
user avatar
8 votes

Is LAPACK behind the cutting edge of dense linear algebra?

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 ...
user avatar
  • 1,976
8 votes
Accepted

Python scipy eigh(Arpack) giving wrong eigenvalues for generalized eigenvalue problem

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 &...
user avatar
7 votes
Accepted

Matrix Balancing Algorithm

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 ...
user avatar
  • 383
7 votes

What's the fastest implementation of elementwise vector multiplication in Fortran?

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 ...
user avatar
7 votes

How to start using LAPACK in c++?

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 ...
user avatar
  • 5,744
7 votes

How to start using LAPACK in c++?

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 <...
user avatar
  • 2,106
7 votes

Matlab, Mathematica & LAPACK returning 3 different eigenvectors

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$...
user avatar
  • 1,876
7 votes

What is the reason that LAPACK uses $\tau$ in QR decomposition (instead of normalizing the reflection vector)?

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 ...
user avatar
  • 4,316
7 votes
Accepted

Do most statistical packages and libraries in high-level programming languages rely on LAPACK for their matrix inversion operations?

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 ...
user avatar
  • 3,111
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 ...
user avatar
  • 1,127
6 votes
Accepted

Functions from Scipy, Blas, or Lapack that compute only upper triangular matrix

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 ...
user avatar
6 votes

Is LAPACK behind the cutting edge of dense linear algebra?

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 ...
user avatar
  • 161
6 votes
Accepted

Faster eigenvector routine for non-symmetric matrices with real eigensystem?

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 ...
user avatar
5 votes

BLAS, LAPACK or ATLAS for Matrix Multiplication in C

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 ...
user avatar
  • 5,744
5 votes
Accepted

How many operations are needed for LAPACK's zgesv to solve a linear system?

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 ...
user avatar
  • 611
5 votes

A misunderstanding or a bug in LAPACK's solver for generalized eigenvalue problems?

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 - \...
user avatar
4 votes
Accepted

Schur(QZ) Decomposition Differences

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 ...
user avatar
4 votes
Accepted

Flop counts for LAPACK symmetric eigenvalue routines DSYEV, DSYEVD, DSYEVX and DSYEVR

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 ...
user avatar
  • 4,410
4 votes
Accepted

Finite Difference Beam Propagation Method problem

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 ...
user avatar
  • 11.4k
4 votes
Accepted

What is wrong with this matrix multiplication?

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-...
user avatar
  • 11.4k
4 votes

Compute all eigenvectors and eigenvalues of small symmetric matrices

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 ...
user avatar
4 votes
Accepted

Efficient computation of AX=B where B has special structure (block-diagonal)

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\...
user avatar
  • 11.4k
4 votes

Rapidly determining whether or not a dense matrix is of low rank

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$. ...
user avatar
  • 81
4 votes
Accepted

LAPACK sorting eigenvalues differently each time

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 ...
user avatar
4 votes

Are there any packaged routines (in lapack or elsewhere) for inverting a banded matrix?

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 ...
user avatar
  • 383
4 votes

LAPACK non-convergent eigenvalue algorithm for complex but not double matrix

Your matrix is not diagonalizable, in the Jordan decomposition of it there is a block for the eigenvalue $0$ of the form $$\begin{pmatrix}0&0&0\\0&0&1\\0&0&0\end{pmatrix},$$ ...
user avatar
  • 11.4k

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