Skip to main content
Share Your Experience: Take the 2024 Developer Survey

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.

Filter by
Sorted by
Tagged with
2 votes
1 answer
159 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 ...
yosuga's user avatar
  • 123
0 votes
1 answer
67 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 ...
nalzok's user avatar
  • 181
2 votes
1 answer
129 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/...
atwood's user avatar
  • 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 ...
tisPrimeTime's user avatar
2 votes
1 answer
173 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 ...
Bananach's user avatar
  • 799
1 vote
0 answers
163 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 ...
Natasha's user avatar
  • 421
4 votes
2 answers
243 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 <...
two_Thomas's user avatar
1 vote
1 answer
152 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 ...
alpx's user avatar
  • 113
5 votes
1 answer
160 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....
eepperly16's user avatar
2 votes
1 answer
222 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. ...
Filippo Caleca's user avatar
3 votes
0 answers
130 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{...
Bananach's user avatar
  • 799
3 votes
1 answer
89 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 ...
Bananach's user avatar
  • 799
1 vote
0 answers
131 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. ...
Urwald's user avatar
  • 111
0 votes
0 answers
90 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/...
Debasish Das's user avatar
4 votes
2 answers
295 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 ...
Dimitar Slavchev's user avatar
2 votes
1 answer
234 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 ...
Matthias Beaupère's user avatar
2 votes
1 answer
597 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 ...
Yaroslav Bulatov's user avatar
1 vote
1 answer
472 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} + ...
IPribec's user avatar
  • 617
4 votes
1 answer
185 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 ...
Federico Poloni's user avatar
4 votes
0 answers
76 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}) = \...
nwknoblauch's user avatar
2 votes
2 answers
704 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. ...
prananna's user avatar
2 votes
1 answer
315 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 ...
Roopayan Ghosh's user avatar
8 votes
2 answers
283 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 &...
user3677630's user avatar
2 votes
2 answers
858 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 ...
Himanshu Chaudhary's user avatar
3 votes
1 answer
151 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 ...
Daniel Kaplan's user avatar
2 votes
0 answers
139 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 ...
Grayscale's user avatar
  • 211
9 votes
3 answers
791 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 ...
Grayscale's user avatar
  • 211
1 vote
1 answer
233 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 ...
Akira's user avatar
  • 207
0 votes
1 answer
215 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 ...
Riri's user avatar
  • 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. ...
user avatar
5 votes
2 answers
430 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,...
StatsScared's user avatar
1 vote
1 answer
188 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 ...
Federico Poloni's user avatar
1 vote
1 answer
146 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, ...
Bananach's user avatar
  • 799
4 votes
2 answers
654 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: ...
vibe's user avatar
  • 1,058
1 vote
1 answer
778 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 ...
jecht300's user avatar
1 vote
0 answers
81 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 ...
tisPrimeTime's user avatar
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 ...
user30200's user avatar
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 ...
geza's user avatar
  • 211
4 votes
1 answer
327 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 ...
fcdimitr's user avatar
  • 141
1 vote
1 answer
267 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 ...
Blablablu's user avatar
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 ...
Endulum's user avatar
  • 735
3 votes
1 answer
482 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 ...
Pedro H. N. Vieira's user avatar
10 votes
2 answers
286 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 ...
Carl Christian's user avatar
4 votes
1 answer
237 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 ...
brubeck's user avatar
  • 63
3 votes
1 answer
180 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 &...
Salsiccio's user avatar
  • 131
1 vote
1 answer
158 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 ...
Blablablu's user avatar
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 ...
MAh2014's user avatar
  • 131
1 vote
1 answer
462 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 ...
Alice Schwarze's user avatar
1 vote
0 answers
291 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.
Puneet Sharma's user avatar
0 votes
1 answer
69 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 ...
Artur's user avatar
  • 3