Questions tagged [lapack]
LAPACK (Linear Algebra PACKage) is a commonly used library of subroutines for numerical linear algebra tasks, including solutions of linear sets of equations, linear least squares, eigenvalue problems, and singular value decomposition. LAPACK routines may be used with fortran, C and relatives and a variety of other languages.
131
questions
9
votes
1
answer
1k
views
Why is it that SVD routines on nearly square matrices run significantly faster than if the matrix was highly non square?
In Python / Matlab, if you run a routine for SVD on a significantly non-square matrix, X, such as X.shape = (2,15000) you will ...
- 203
2
votes
1
answer
161
views
Fast weighted vector inner product x*A*y with BLAS/LAPACK
Is there a way to compute the weighted vector inner product xAy with vectors x and y and Matrix A using BLAS/LAPACK while avoiding additional allocations or overwriting the inputs?
I'm happy with ...
- 799
1
vote
0
answers
132
views
Using Sundials CVODE in MATLAB
I'm currently using ode15s to solve a set of stiff differential equations.
I am trying to use the MATLAB profiler to understand the section of the ode solver code which calls BLAS routines.
Since the ...
- 421
2
votes
2
answers
187
views
Reordering eigenvalues in Schur factorization - MATLAB orschur and LAPACK dtrsen not producing the same results
Disclaimer: I previously posted this on SO, but though it would be more relevant for scicomp. The original post has been deleted.
I have been trying to recreate the functionality provided by MATLABs <...
- 21
1
vote
1
answer
112
views
Conditioning and Stability of generalized eigenvalue problem
The (generalized) eigenvalue problems with a multiple eigenvalue are the ill-posed ones.
I have two questions that should be simple for experts:
(1) Is the eigenvalue problem much more sensitive to ...
- 113
0
votes
0
answers
56
views
Rank-one updates for symmetric matrix eigen-system
Are there existing implementations for rank-one updating of symmetric matrices eigensystems?
This is the mathematical statement of the problem. Let $S=QDQ^T$
$$S + vv^T = QDQ^T +vv^T = Q_{new}D_{new}...
4
votes
1
answer
136
views
Eigenvalues of diagonal plus rank-one
I need to compute an eigendecomposition of an $n\times n$ matrix
$$
D + c vv^\top = Q\Lambda Q^\top \tag{1}
$$
in MATLAB, where $D$ is a real diagonal matrix, $c$ is a scalar, and $v$ is a real vector....
- 350
2
votes
1
answer
199
views
Parallelize pseudo inverse of a matrix using Lapacke
I am currently using the protocol described in
https://stackoverflow.com/questions/55599950/computation-of-pseidoinverse-with-svd-in-c-using-blas-and-lapacke to compute the pseudo inverse of a matrix.
...
3
votes
0
answers
106
views
Compute orthogonal complement using BLAS / LAPACK
Is there a fast method to compute an orthogonal complement of an arbitrary matrix $U\in\mathbb{R}^{m \times n}$ in BLAS / LAPACK?
Specifically, I want any matrix $V\in \mathbb{R}^{m \times (m - \text{...
- 799
3
votes
1
answer
87
views
Reuse linear mapping that provides the solution to least squares problem using LAPACK
LAPACK.gglse allows me to solve
min x^T Q x
s.t. A x = y
(in my present use case, $Q$ is symmetric positive definite)
without having to think about the numerical ...
- 799
1
vote
0
answers
112
views
Does cblas_dgemm mutate my input matrices?
I have written a matrix class Matrix<T> for which I have implemented a wrapper function for cblas_dgemm.
...
- 111
0
votes
0
answers
86
views
PLASMA usage (Linear algebra routines that supports multithreading)
I have been looking for linear algebra libraries that support multithreading. I have found PLASMA which looks promising. It is from the same group that developed LAPACK.
http://icl.cs.utk.edu/...
- 45
4
votes
2
answers
265
views
Is there a simple way to add a sparse matrix to an LU decomposition of a dense matrix?
I am solving a parabolic equation in the form:
$$
\left( {M \over\tau_j} + A \right) u^{j+1} = M f^j + {u^j \over \tau_j},
$$
where $A$ and $f$ are a dense stiffness matrix and the right hand side of ...
- 183
2
votes
1
answer
215
views
Efficiently compute a projection matrix from Householders reflectors
Let $A \in \mathbb{R}^{m \times n}$ where $m \geq n$.
Let $B$ and $\tau$ be the result of applying LAPACK's dgeqrfp method (R on the upper right triangle, and the ...
2
votes
1
answer
537
views
Faster eigenvector routine for non-symmetric matrices with real eigensystem?
I have non-symmetric real-valued matrices with real-valued eigensystems. How to compute eigenvectors efficiently?
Using scipy.linalg.eig (which calls ...
- 2,273
1
vote
1
answer
451
views
Factorization of cubic spline interpolation matrix
In cubic spline interpolation, we use the set of knots and function values $(x_i,y_i),i=1,...,n$ to construct a (tridiagonal) system of equations for the unknowns $\sigma_i$:
$$
h_{i-1}\sigma_{i-1} + ...
- 544
4
votes
1
answer
174
views
Trace of inverse from LU decomposition
Given an LU decomposition of $A\in \mathbb{R}^{n\times n}$, is there a way to compute $\operatorname{trace}(A^{-1})$ with lower complexity than that of the inversion ($O(n^3)$ in practice)?
This ...
- 11.2k
4
votes
0
answers
74
views
Efficient computation of marginalized multivariate normal likelihood
In general,if we know that the marginal Gaussian distribution for some variable $\textbf{x}$ and a conditional Gaussian distribution for some $\textbf{y}|\textbf{x}$ of the forms:
$$p(\textbf{x}) = \...
- 41
2
votes
2
answers
637
views
Does LAPACK offer routines for Krylov sub-space based solvers and nonlinear solvers?
I have skimmed through the LAPACK user guide, but I could not find if LAPACK offers routines for Krylov Subspace based methods (such as CG or BiCGSTAB etc) and Newton method based nonlinear solvers. ...
- 33
2
votes
1
answer
293
views
Diagonalization of Hermitian matrices vs Unitary matrices
What are the general algorithms used for diagonalization of large Hermitian matrices and Unitary matrices? ($>5000 \times 5000$)
LAPACK seems to diagonalize Hermitian matrices almost 20 times as ...
- 133
8
votes
2
answers
256
views
A misunderstanding or a bug in LAPACK's solver for generalized eigenvalue problems?
In my application, I have two general real matrices $A$,$B$ defined as follows,
$$
A=\begin{bmatrix}
-s I_3 & A_0 & 0 & 0 \\
A_0^T & -s I_3 & 0 & 0 \\
0 &...
- 205
2
votes
2
answers
777
views
Python scipy eigh(Arpack) giving wrong eigenvalues for generalized eigenvalue problem
I am trying to solve a generalized eigenvalue problem using Arpack, right now the code is using LAPACK but that's too slow, we only need a few eigenvalues and the matrices are sparse so using Arpack ...
3
votes
0
answers
103
views
Numerical calculation of the Berry connection
I'm doing some numerical calculations involving Hermitian matrices, and derivatives of the eigenvectors.
Essentially, I have an n x n, Hermitian matrix H(x), which is dependent on some continuous ...
2
votes
0
answers
126
views
Why is LAPACK (seemingly) suboptimal for packed and banded eigenvalue problems?
Based on this LAPACK routines list, it looks like there is no relatively robust representation (RRR) driver routine for either packed or banded symmetric eigenvalue problems. According to the relevant ...
- 201
8
votes
3
answers
647
views
Is LAPACK behind the cutting edge of dense linear algebra?
I have been digging into some numerical linear algebra lately, and reading in particular about how LAPACK solves symmetric eigenvalue problems. I noticed that the ...
- 201
1
vote
1
answer
202
views
Functions from Scipy, Blas, or Lapack that compute only upper triangular matrix
My goal is to transform a matrix into upper triangular form in Python. I know the function scipy.linalg.lu will do LU decomposition and get both upper and lower ...
- 207
0
votes
1
answer
191
views
Library to solve dense linear system with GMRES
I have a fortran 90 code and I want to solve a dense linear system with GMRES. I would prefer the restarted GMRES with preconditioning. Is there some library that you know of that I could use? Now I ...
- 43
5
votes
1
answer
1k
views
How to know which LAPACK's function is used by Scipy's eig function?
As far as I understood, scipy.linalg.eig use wrappers from scipy.lapack to compute the eigenvalues and eigenvectors of a matrix. ...
user33403
5
votes
2
answers
413
views
Do most statistical packages and libraries in high-level programming languages rely on LAPACK for their matrix inversion operations?
Possible an open-ended question, but I am wondering if most statistical packages and libraries, for instance, Stata, R, Python's NumPy and MATLAB rely on LAPACK algorithms to perform matrix operations,...
- 153
1
vote
1
answer
177
views
Lapack symmetric update $B^{-1}AB^{-T}$
Does Lapack have a routine that, given symmetric $A=A^T$ and $B$, computes the symmetric matrix $B^{-1}AB^{-T}$ (while preserving symmetry exactly)?
It would be enough to have this routine for ...
- 11.2k
1
vote
1
answer
142
views
Using LAPACK to compute $B^{-1}AB^{-T}$ for thin $B$
How can I use BLAS/LAPACK to compute
$$
B^{-1}AB^{-T}
$$
where $A\in\mathbb{R}^{n,n}$, $B\in\mathbb{R}^{m,n}$ is full rank matrix with $m>n$, and $B^{-1}y:=\arg \min_{x} \|Bx-y\|_{2}$.
In theory, ...
- 799
4
votes
2
answers
621
views
Inverting a matrix from LU decomposition
The LAPACK routines xGETRI compute the inverse of a matrix $A = PLU$ in its LU decomposed form by first computing $U^{-1}$, and then solving the system:
$$
(A^{-1} P) L = U^{-1}
$$
My question is: ...
- 1,048
1
vote
1
answer
704
views
How to use LAPACK function (DGELSY) in Fortran
I am trying to use Least Squares Minimization to solve a the matrix problem: b = A*x for x. The system is overdetermined, and A is a dense matrix.
In the LAPACK library, I believe the routine DGELSY ...
- 13
1
vote
0
answers
76
views
Stability of SVD, Eigendecompositions, and pseudoinverse procedures in modern LAPACK routines
I have proposed an optimisation algorithm which I claim has improved upon the previous algorithm in a number of ways.
One of these claims is that my proposed solution requires no explicit SVD and ...
- 203
2
votes
1
answer
3k
views
LAPACK equivalent on c++ , which is the best one? [duplicate]
I am following a course of computational material physics. The professor uses fortran to code and uses lapack to solve eigenvalue problems.
So far I just know c++. There is an equivalent library that ...
- 21
11
votes
3
answers
1k
views
What is the reason that LAPACK uses $\tau$ in QR decomposition (instead of normalizing the reflection vector)?
LAPACK's QR routine stores Q as Householder reflectors. It scales the reflection vector $v$ with $1/v_1$, so the first element of the result becomes $1$, so it doesn't have to be stored. And it stores ...
- 211
4
votes
1
answer
314
views
Fast matrix multiplication with matrix elements computed on-the-fly (without forming the matrix)
Is there any library or routine for high-performance matrix-matrix product, where the matrix elements are computed on-the-fly using a given function of $i$ and $j$?
More specifically, in the problem ...
- 141
1
vote
1
answer
253
views
Wrong result of 'ddot' from BLAS
I am having trouble with a C/C++ program that uses the BLAS routine ddot.
I am running Linux and so far LAPACK routines worked without any problems.
I get a wrong ...
- 25
2
votes
1
answer
98
views
Inverting small matrices: canned factorization versus explicit formula
I am interested in solving a large number of small linear systems of equations, $Ax=b$, with $A$ either $2\times2$ or $3\times3$. Assuming none of these systems are actually singular, is there ...
- 725
3
votes
1
answer
453
views
How many operations are needed for LAPACK's zgesv to solve a linear system?
I have a linear system of complex numbers. I am using LAPACK' zgesv (actually I am using intel MKL LAPACKE, but I am assuming the algorithm is the same). No assumption can be made about the system.
I ...
10
votes
2
answers
273
views
Benchmark problems for eigenvalue reordering algorithms sought
Every real matrix $A$ can be reduce to real Schur form $T = U^T A U$ using an orthogonal similiary transform $U$. Here the matrix $T$ is quasi-triangular form with 1 by 1 or 2 by 2 blocks on the main ...
- 1,391
4
votes
1
answer
226
views
Compute bilinear form with LAPACK
I need to compute a bilinear form for a set of left and right vectors
$$ w_k = \sum_{i,j} V_{ik}^*A_{ij}U_{jk},$$
where $A_{ij}\in\mathbb{R}$ and $U_{jk}, V_{ik} \in \mathbb{C}$ (Assume that all the ...
- 63
3
votes
1
answer
176
views
Diagonalizing a block tridiagonal Toepliz Hermitean matrix
I have to diagonalize, within a Fortran-written code, a block tridiagonal Toeplitz Hermitian matrix, e.g.
$$
\left[
\begin{array}{ccccc}
\ddots & \hat{A} & & & \\
\hat{A}^\dagger &...
- 131
1
vote
1
answer
147
views
Problem with 'dsysv' from LAPACK
I am having trouble with a C program that uses the function dsysv from LAPACK. Everything compiles and works without any errors, my only problem is, that the ...
- 25
0
votes
2
answers
2k
views
Does armadillo library slow down the execution of matrix operations?
I've converted a MATLAB code to C++ to speed it up, using the Armadillo library to handle matrix operations in C++, but surprisingly it is 10 times slower than the MATLAB code!
So I test the ...
- 131
1
vote
1
answer
452
views
C++: Efficient library for dense linear algebra operations (determinant & principal minors)
I usually work with Python, but my basic knowledge of c++ allows me to switch when I need to increase the speed of my code.
Currently, I have a python script that (among other things) computes the ...
- 113
1
vote
0
answers
267
views
Kronecker product of matrices
I have to use Kronecker product of a matrix with a unit matrix is there any routine in ScaLapack or Lapack which can do so efficiently.
0
votes
1
answer
66
views
LAPACK zlapmt_ freezes code [closed]
I'm using LAPACK's zggev_ routine to solve some generalized eigenvalue problem.
While it produces the correct results, I want the eigenvalues and according eigenvectors sorted by absolute value.
For ...
- 3
0
votes
1
answer
433
views
How to do a Generalized Complex Schur (or QZ) Decomposition with Eigen C++? [closed]
I would like to do a Generalized Schur (or QZ) decomposition for a pair of complex matrices $A$ and $B$.
I found the following class:
...
- 101
1
vote
0
answers
126
views
Eigenvalue problem (LAPACK)
I am working on a project in numerical analysis which I have to program in C (using Lapack and Blas). Matrix is given which is tridiagonal and "almost" symmetric (one element is to be changed to make ...
- 111