LAPACK

Interface to LAPACK routines. More...

Functions
template<typename Array2d , typename Array1d >
void	cusp::lapack::getrf (Array2d &A, Array1d &piv)
	Compute LU factorization of matrix. More...
template<typename Array2d >
void	cusp::lapack::potrf (Array2d &A, char uplo= 'U')
	Computes the Cholesky factorization of a symmetric (Hermitian) positive-definite matrix. More...
template<typename Array2d , typename Array1d >
void	cusp::lapack::sytrf (Array2d &A, Array1d &piv, char uplo= 'U')
	Computes the Bunch-Kaufman factorization of a symmetric matrix. More...
template<typename Array2d , typename Array1d >
void	cusp::lapack::getrs (const Array2d &A, const Array1d &piv, Array2d &B, char trans= 'N')
	Solves a system of linear equations with an LU-factored square matrix, with multiple right-hand sides. More...
template<typename Array2d >
void	cusp::lapack::potrs (const Array2d &A, Array2d &B, char uplo= 'U')
	Solves a system of linear equations with a Cholesky-factored symmetric (Hermitian) positive-definite matrix. More...
template<typename Array2d , typename Array1d >
void	cusp::lapack::sytrs (const Array2d &A, const Array1d &piv, Array2d &B, char uplo= 'U')
	Solves a system of linear equations with a UDU- or LDL-factored symmetric matrix. More...
template<typename Array2d >
void	cusp::lapack::trtrs (const Array2d &A, Array2d &B, char uplo= 'U', char trans= 'N', char diag= 'N')
	Solves a system of linear equations with a triangular matrix, with multiple right-hand sides. More...
template<typename Array2d >
void	cusp::lapack::trtri (Array2d &A, char uplo= 'U', char diag= 'N')
	This routine computes the inverse of a triangular matrix. More...
template<typename Array2d , typename Array1d >
void	cusp::lapack::syev (const Array2d &A, Array1d &eigvals, Array2d &eigvecs, char uplo= 'U')
	Computes all eigenvalues and eigenvectors of a real symmetric matrix. More...
template<typename Array1d1 , typename Array1d2 , typename Array1d3 , typename Array2d >
void	cusp::lapack::stev (const Array1d1 &alphas, const Array1d2 &betas, Array1d3 &eigvals, Array2d &eigvecs, char job= 'V')
	Computes all eigenvalues and, optionally, eigenvectors of a real symmetric tridiagonal matrix. More...
template<typename Array2d1 , typename Array2d2 , typename Array1d , typename Array2d3 >
void	cusp::lapack::sygv (const Array2d1 &A, const Array2d2 &B, Array1d &eigvals, Array2d3 &eigvecs)
	Computes all eigenvalues and, optionally, eigenvectors of a real generalized symmetric definite eigenproblem. More...
template<typename Array2d , typename Array1d >
void	cusp::lapack::gesv (const Array2d &A, Array2d &B, Array1d &piv)
	Computes the solution of a system of linear equations with a square matrix A and multiple right-hand sides. More...

Function Documentation

template<typename Array2d , typename Array1d >

void cusp::lapack::gesv	(	const Array2d &	A,
		Array2d &	B,
		Array1d &	piv
	)

Computes the solution of a system of linear equations with a square matrix A and multiple right-hand sides.

Template Parameters

Array2d	Type of the input matrix
Array1d	Type of the input pivots array

Parameters

A	Contains the upper or lower triangle of a symmetric matrix
B	Contains the upper or lower triangle of a symmetric positive definite matrix
piv	Array containing pivots

Overview

This routine solves the system of linear equations A*X = B, where A is an n-by-n matrix for X, the columns of matrix B are individual right-hand sides, and the columns of X are the corresponding solutions.

An LU decomposition with partial pivoting and row interchanges is used to factor A as A = P*L*U, where P is a permutation matrix, L is unit lower triangular, and U is upper triangular. The factored form of A is then used to solve the system of equations A*X = B.

Example

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;
  cusp::array1d<float,cusp::host_memory> piv;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // solve multiple RHS vectors
  cusp::lapack::gesv(A, B, piv);

  // print the contents of B
  cusp::print(B);
}

template<typename Array2d , typename Array1d >

void cusp::lapack::getrf	(	Array2d &	A,
		Array1d &	piv
	)

Compute LU factorization of matrix.

Template Parameters

Array2d	Type of the input matrix to factor
Array1d	Type of pivot array

Parameters

A	Input matrix to factor
piv	Array containing pivots

Overview

This routine computes the LU factorization of a general m-by-n matrix A as A = P*L*U, where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n) and U is upper triangular (upper trapezoidal if m < n). The routine uses partial pivoting, with row interchanges.

Example

#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;
  cusp::array1d<float,cusp::host_memory> piv;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // compute LU factorization of A
  cusp::lapack::getrf(A, piv);

  // print the contents of A
  cusp::print(A);
  // print the contents of piv
  cusp::print(piv);
}

template<typename Array2d , typename Array1d >

void cusp::lapack::getrs	(	const Array2d &	A,
		const Array1d &	piv,
		Array2d &	B,
		char	trans = 'N'
	)

Solves a system of linear equations with an LU-factored square matrix, with multiple right-hand sides.

Template Parameters

Array2d	Type of the input matrices
Array1d	Type of pivot array

Parameters

A	LU factored input matrix
piv	Array containing pivots
B	matrix containing multiple right-hand side vectors
trans	If 'N', then A*X = B is solved for X.

Overview

This routine solves for X the following systems of linear equations: A*X = B if trans='N', AT*X = B if trans='T', AH*X = B if trans='C' (for complex matrices only).

Note

Before calling this routine, you must call getrf to compute the LU factorization of A.

Example

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;
  cusp::array1d<float,cusp::host_memory> piv;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // create initial RHS of 2 vectors and initialize
  cusp::array2d<float,cusp::host_memory> B(A.num_rows, 2);
  B.values = cusp::random_array<float>(B.values.size());

  // compute LU factorization of A
  cusp::lapack::getrf(A, piv);
  // solve multiple RHS vectors
  cusp::lapack::getrs(A, piv, B);

  // print the contents of B
  cusp::print(B);
}

template<typename Array2d >

void cusp::lapack::potrf	(	Array2d &	A,
		char	uplo = 'U'
	)

Computes the Cholesky factorization of a symmetric (Hermitian) positive-definite matrix.

Template Parameters

Array2d

Type of the input matrix to factor

Parameters

A	Input matrix to factor
uplo	Indicates whether A is upper or lower triangular

Overview

This routine forms the Cholesky factorization of a symmetric positive-definite or, for complex data, Hermitian positive-definite matrix A: A = UT*U for real data, A = UH*U for complex data if uplo='U' A = L*LT for real data, A = L*LH for complex data if uplo='L' where L is a lower triangular matrix and U is upper triangular.

Example

#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // compute Cholesky factorization of A
  cusp::lapack::potrf(A);

  // print the contents of A
  cusp::print(A);
}

template<typename Array2d >

void cusp::lapack::potrs	(	const Array2d &	A,
		Array2d &	B,
		char	uplo = 'U'
	)

Solves a system of linear equations with a Cholesky-factored symmetric (Hermitian) positive-definite matrix.

Template Parameters

Array2d

Type of the input matrices

Parameters

A	Cholesky factored input matrix
B	matrix containing multiple right-hand side vectors
uplo	Indicates whether A is upper or lower triangular

Overview

This routine solves for X the system of linear equations A*X = B with a symmetric positive-definite or, for complex data, Hermitian positive-definite matrix A, given the Cholesky factorization of A:

A = UT*U for real data, A = UH*U for complex dataif UHplo='U' A = L*LT for real data, A = L*LH for complex dataif uplo='L'

where L is a lower triangular matrix and U is upper triangular. The system is solved with multiple right-hand sides stored in the columns of the matrix B.

Note

Before calling this routine, you must call potrf to compute the Cholesky factorization of A.

Example

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // create initial RHS of 2 vectors and initialize
  cusp::array2d<float,cusp::host_memory> B(A.num_rows, 2);
  B.values = cusp::random_array<float>(B.values.size());

  // compute Bunch-Kaufman factorization of A
  cusp::lapack::potrf(A);

  // solve multiple RHS vectors
  cusp::lapack::potrs(A, B);

  // print the contents of B
  cusp::print(B);
}

template<typename Array1d1 , typename Array1d2 , typename Array1d3 , typename Array2d >

void cusp::lapack::stev	(	const Array1d1 &	alphas,
		const Array1d2 &	betas,
		Array1d3 &	eigvals,
		Array2d &	eigvecs,
		char	job = 'V'
	)

Computes all eigenvalues and, optionally, eigenvectors of a real symmetric tridiagonal matrix.

Template Parameters

Array1d1	Type of the input array
Array1d2	Type of the input array
Array1d3	Type of the input array
Array2d	Type of the input array

Parameters

alphas	Main diagonal entries
betas	First sub-diagonal entries
eigvals	On return contain eigenvalues of matrix
eigvecs	On return contain eigenvectors of the matrix
job	If 'V', then eigenvalues and eigenvectors are computed

Overview

This routine computes all eigenvalues and eigenvectors of a real symmetric matrix.

Example

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // solve multiple RHS vectors
  cusp::lapack::syev(A, eigvals, eigvecs);

  // print the contents of B
  cusp::print(eigvals);
}

template<typename Array2d , typename Array1d >

void cusp::lapack::syev	(	const Array2d &	A,
		Array1d &	eigvals,
		Array2d &	eigvecs,
		char	uplo = 'U'
	)

Computes all eigenvalues and eigenvectors of a real symmetric matrix.

Template Parameters

Array2d	Type of the input matrices
Array1d	Type of the input array

Parameters

A	Symmetric input matrix
eigvals	On return contain eigenvalues of matrix
eigvecs	On return contain eigenvectors of the matrix
uplo	Indicates whether upper or lower portion of symmetric matrix is stored.

Overview

This routine computes all eigenvalues and eigenvectors of a real symmetric matrix.

Example

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // solve multiple RHS vectors
  cusp::lapack::syev(A, eigvals, eigvecs);

  // print the contents of B
  cusp::print(eigvals);
}

template<typename Array2d1 , typename Array2d2 , typename Array1d , typename Array2d3 >

void cusp::lapack::sygv	(	const Array2d1 &	A,
		const Array2d2 &	B,
		Array1d &	eigvals,
		Array2d3 &	eigvecs
	)

Computes all eigenvalues and, optionally, eigenvectors of a real generalized symmetric definite eigenproblem.

Template Parameters

Array2d1	Type of the input array
Array2d2	Type of the input array
Array1d	Type of the input array
Array2d3	Type of the input array

Parameters

A	Contains the upper or lower triangle of a symmetric matrix
B	Contains the upper or lower triangle of a symmetric positive definite matrix
eigvals	On return contain eigenvalues of matrix
eigvecs	On return contain eigenvectors of the matrix

Overview

This routine computes all the eigenvalues, and optionally, the eigenvectors of a real generalized symmetric-definite eigenproblem, of the form

A*x = λ*B*x, A*B*x = λ*x, or B*A*x = λ*x.

Here A and B are assumed to be symmetric and B is also positive definite.

Example

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // compute eigvals and eigvecs
  cusp::lapack::sygv(A, A, eigvals, eigvecs);

  // print the eigenvalues
  cusp::print(eigvals);
}

template<typename Array2d , typename Array1d >

void cusp::lapack::sytrf	(	Array2d &	A,
		Array1d &	piv,
		char	uplo = 'U'
	)

Computes the Bunch-Kaufman factorization of a symmetric matrix.

Template Parameters

Array2d	Type of the input matrix to factor
Array1d	Type of pivot array

Parameters

A	Input matrix to factor
piv	Array containing pivots
uplo	Indicates whether A is upper or lower triangular

Overview

This routine computes the factorization of a symmetric or hermitian matrix using the Bunch-Kaufman diagonal pivoting method. The form of the factorization is: A = P*U*D*UT*PT, if uplo='U' A = P*L*D*LT*PT, if uplo='L' where A is the input matrix, P is a permutation matrix, U and L are upper and lower triangular matrices with unit diagonal, and D is a symmetric block-diagonal matrix with 1-by-1 and 2-by-2 diagonal blocks. U and L have 2-by-2 unit diagonal blocks corresponding to the 2-by-2 blocks of D.

Example

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;
  cusp::array1d<float,cusp::host_memory> piv;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // compute Bunch-Kaufman factorization of A
  cusp::lapack::sytrf(A, piv);

  // print the contents of A
  cusp::print(A);
  // print the contents of piv
  cusp::print(piv);
}

template<typename Array2d , typename Array1d >

void cusp::lapack::sytrs	(	const Array2d &	A,
		const Array1d &	piv,
		Array2d &	B,
		char	uplo = 'U'
	)

Solves a system of linear equations with a UDU- or LDL-factored symmetric matrix.

Template Parameters

Array2d	Type of the input matrices
Array1d	Type of pivot array

Parameters

A	Input matrix to factor
piv	Array containing pivots
B	matrix containing multiple right-hand side vectors
uplo	Indicates whether A is upper or lower triangular

Overview

The routine solves for X the system of linear equations A*X = B with a symmetric matrix A, given the Bunch-Kaufman factorization of A:

A = P*U*D*UT*PT, if uplo='U' A = P*L*D*LT*PT, if uplo='L'

where P is a permutation matrix, U and L are upper and lower triangular matrices with unit diagonal, and D is a symmetric block-diagonal matrix. The system is solved with multiple right-hand sides stored in the columns of the matrix B.

Note

You must supply the factor U (or L) and the array of pivots returned by the factorization routine sytrf.

Example

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;
  cusp::array1d<float,cusp::host_memory> piv;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // create initial RHS of 2 vectors and initialize
  cusp::array2d<float,cusp::host_memory> B(A.num_rows, 2);
  B.values = cusp::random_array<float>(B.values.size());

  // compute Bunch-Kaufman factorization of A
  cusp::lapack::sytrf(A, piv);

  // solve multiple RHS vectors
  cusp::lapack::sytrs(A, piv, B);

  // print the contents of B
  cusp::print(B);
}

template<typename Array2d >

void cusp::lapack::trtri	(	Array2d &	A,
		char	uplo = 'U',
		char	diag = 'N'
	)

This routine computes the inverse of a triangular matrix.

Template Parameters

Array2d

Type of the input matrix

Parameters

A	Triangular input matrix
uplo	Indicates whether A is upper or lower triangular
diag	If 'N', then A is not a unit triangular matrix.

Example

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // print the contents of A
  cusp::print(A);
}

template<typename Array2d >

void cusp::lapack::trtrs	(	const Array2d &	A,
		Array2d &	B,
		char	uplo = 'U',
		char	trans = 'N',
		char	diag = 'N'
	)

Solves a system of linear equations with a triangular matrix, with multiple right-hand sides.

Template Parameters

Array2d

Type of the input matrices

Parameters

A	Triangular input matrix
B	matrix containing multiple right-hand side vectors
uplo	Indicates whether A is upper or lower triangular
trans	If 'N', then A*X = B is solved for X.
diag	If 'N', then A is not a unit triangular matrix.

Overview

This routine solves for X the following systems of linear equations with a triangular matrix A, with multiple right-hand sides stored in B: A*X = B if trans='N', AT*X = B if trans='T', AH*X = B if trans='C' (for complex matrices only).

Example

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>

// include cusp lapack header file
#include <cusp/lapack/lapack.h>

int main()
{
  // create an empty dense matrix structure
  cusp::array2d<float,cusp::host_memory> A;

  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);

  // create initial RHS of 2 vectors and initialize
  cusp::array2d<float,cusp::host_memory> B(A.num_rows, 2);
  B.values = cusp::random_array<float>(B.values.size());

  // solve multiple RHS vectors
  cusp::lapack::trtrs(A, B);

  // print the contents of B
  cusp::print(B);
}