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.
136
questions
2
votes
0
answers
47
views
MATLAB `ordschur` producing different results from Lapack's `dtrsen` in Simulink
My question is very similar to the one previously posted here. I'm using the same MATLAB code for embedding Lapack functions into my code, which I also compiled with Mingw64. The main difference is ...
- 21
0
votes
0
answers
20
views
Use of scipy.linalg.cython_lapack.dgbsv in cython script
I was working into accellerating a solver for baded matrix using cython.
The base case is the following pure python code :
...
- 35
3
votes
1
answer
111
views
2x2 complex symmetric eigendecomposition - LAPACK subroutine CLAESY
Asking here because I searched the LAPACK user forums and found nothing.
I have a problem that requires the computation of the eigendecomposition $A=A^T=Q \Lambda Q^T$ for the 2x2 complex symmetric ...
- 1,031
4
votes
1
answer
249
views
Matrix Diagonalization and Computational Requirements
I have some questions about diagonalizing matrices. My interest lies in computing all eigenvalues of a given matrix. To avoid wasting time and improve my research efficiency, I want to understand the ...
- 143
0
votes
1
answer
77
views
Weird runtime behavior of `scipy.linalg.solve_triangular` and `trtrs`
I want to understand the time complexity of scipy.linalg.solve_triangular, which calls trtrs from LAPACK under the hood, so I ...
- 181
2
votes
1
answer
138
views
What is the correct way of computing $LL^T$ in packed format with blas/lapack
I have a triangular matrix stored in packed format (ie $L$). I need to compute $LL^T$ (not the decomposition, just the multiplication). What would be the preferred way of computing this with blas/...
- 23
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
186
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
2
votes
1
answer
205
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 ...
- 433
4
votes
2
answers
271
views
Reordering eigenvalues in Schur factorization - MATLAB ordschur 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 <...
- 41
1
vote
1
answer
212
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
5
votes
1
answer
189
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....
- 380
2
votes
1
answer
236
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
139
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
93
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
142
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
93
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
319
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
246
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
645
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,735
1
vote
1
answer
485
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} + ...
- 637
4
votes
1
answer
207
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 ...
- 12k
4
votes
0
answers
77
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
747
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
342
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
298
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 &...
- 215
2
votes
2
answers
910
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
1
answer
156
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
143
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 ...
- 211
9
votes
3
answers
888
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 ...
- 211
1
vote
1
answer
242
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
230
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
434
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
194
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 ...
- 12k
1
vote
1
answer
148
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
687
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,068
1
vote
1
answer
816
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
86
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
337
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
273
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
99
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 ...
- 735
3
votes
1
answer
497
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
291
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,450
4
votes
1
answer
250
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
189
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
162
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