Matrix Algorithms

Algorithms for processing matrices in sparse and dense formats. More...

Functions
template<typename SourceType , typename DestinationType >
void	cusp::convert (const SourceType &src, DestinationType &dst)
	Convert between matrix formats.
template<typename SourceType , typename DestinationType >
void	cusp::copy (const SourceType &src, DestinationType &dst)
	Copy one array or matrix to another.
template<typename MatrixType1 , typename MatrixType2 , typename MatrixType3 , typename BinaryFunction >
void	cusp::elementwise (const MatrixType1 &A, const MatrixType2 &B, MatrixType3 &C, BinaryFunction op)
	Perform transform operation on two matrices.
template<typename MatrixType1 , typename MatrixType2 , typename MatrixType3 >
void	cusp::add (const MatrixType1 &A, const MatrixType2 &B, MatrixType3 &C)
	Compute the sum of two matrices.
template<typename MatrixType1 , typename MatrixType2 , typename MatrixType3 >
void	cusp::subtract (const MatrixType1 &A, const MatrixType2 &B, MatrixType3 &C)
	Compute the difference of two matrices.
template<typename OffsetArray , typename IndexArray >
void	cusp::offsets_to_indices (const OffsetArray &offsets, IndexArray &indices)
	Expand CSR row offsets to COO row indices.
template<typename IndexArray , typename OffsetArray >
void	cusp::indices_to_offsets (const IndexArray &indices, OffsetArray &offsets)
	Compress COO row indices to CSR row offsets.
template<typename MatrixType , typename ArrayType >
void	cusp::extract_diagonal (const MatrixType &A, ArrayType &output)
	Extract the main diagonal of a matrix.
template<typename ArrayType1 , typename ArrayType2 >
size_t	cusp::count_diagonals (const size_t num_rows, const size_t num_cols, const ArrayType1 &row_indices, const ArrayType2 &column_indices)
	Count the number of occupied diagonals in the input matrix.
template<typename ArrayType >
size_t	cusp::compute_max_entries_per_row (const ArrayType &row_offsets)
	Compute the maximum row length of a matrix.
template<typename ArrayType >
size_t	cusp::compute_optimal_entries_per_row (const ArrayType &row_offsets, float relative_speed=3.0f, size_t breakeven_threshold=4096)
	Compute the optimal number of entries per row of HYB matrix.
template<typename LinearOperator , typename MatrixOrVector1 , typename MatrixOrVector2 >
void	cusp::multiply (const LinearOperator &A, const MatrixOrVector1 &B, MatrixOrVector2 &C)
	Implements matrix-matrix and matrix-vector multiplication.
template<typename LinearOperator , typename MatrixOrVector1 , typename MatrixOrVector2 , typename UnaryFunction , typename BinaryFunction1 , typename BinaryFunction2 >
void	cusp::multiply (const LinearOperator &A, const MatrixOrVector1 &B, MatrixOrVector2 &C, UnaryFunction initialize, BinaryFunction1 combine, BinaryFunction2 reduce)
	Implements matrix-vector multiplication with custom combine and reduce functionality.
template<typename LinearOperator , typename MatrixOrVector1 , typename MatrixOrVector2 , typename UnaryFunction , typename BinaryFunction1 , typename BinaryFunction2 >
void	cusp::generalized_spgemm (const LinearOperator &A, const MatrixOrVector1 &B, MatrixOrVector2 &C, UnaryFunction initialize, BinaryFunction1 combine, BinaryFunction2 reduce)
	Implements generalized matrix-matrix multiplication.
template<typename LinearOperator , typename Vector1 , typename Vector2 , typename Vector3 , typename BinaryFunction1 , typename BinaryFunction2 >
void	cusp::generalized_spmv (const LinearOperator &A, const Vector1 &x, const Vector2 &y, Vector3 &z, BinaryFunction1 combine, BinaryFunction2 reduce)
	Implements generalized matrix-vector multiplication.
template<typename ArrayType >
void	cusp::counting_sort (ArrayType &keys, typename ArrayType::value_type min, typename ArrayType::value_type max)
	Use counting sort to order an array.
template<typename ArrayType1 , typename ArrayType2 >
void	cusp::counting_sort_by_key (ArrayType1 &keys, ArrayType2 &vals, typename ArrayType1::value_type min, typename ArrayType1::value_type max)
	Use counting sort to order an array and permute an array of values.
template<typename ArrayType1 , typename ArrayType2 , typename ArrayType3 >
void	cusp::sort_by_row (ArrayType1 &row_indices, ArrayType2 &column_indices, ArrayType3 &values, typename ArrayType1::value_type min_row=0, typename ArrayType1::value_type max_row=0)
	Sort matrix indices by row.
template<typename ArrayType1 , typename ArrayType2 , typename ArrayType3 >
void	cusp::sort_by_row_and_column (ArrayType1 &row_indices, ArrayType2 &column_indices, ArrayType3 &values, typename ArrayType1::value_type min_row=0, typename ArrayType1::value_type max_row=0, typename ArrayType2::value_type min_col=0, typename ArrayType2::value_type max_col=0)
	Sort matrix indices by row and column.
template<typename MatrixType1 , typename MatrixType2 >
void	cusp::transpose (const MatrixType1 &A, MatrixType2 &At)
	Transpose a matrix.
template<typename MatrixType >
bool	cusp::is_valid_matrix (const MatrixType &matrix)
	Validate format of a given matrix.
template<typename MatrixType , typename OutputStream >
bool	cusp::is_valid_matrix (const MatrixType &matrix, OutputStream &ostream)
	Validate format of a given matrix.
template<typename MatrixType >
void	cusp::assert_is_valid_matrix (const MatrixType &matrix)
	Validate format of a given matrix and exit if invalid.
template<typename Array1 , typename Array2 >
void	cusp::assert_same_dimensions (const Array1 &array1, const Array2 &array2)
	Assert dimensions of two arrays are equal. Throw invalid input exception if dimensions do not match.
template<typename Array1 , typename Array2 , typename Array3 >
void	cusp::assert_same_dimensions (const Array1 &array1, const Array2 &array2, const Array3 &array3)
	Assert dimensions of two arrays are equal. Throw invalid input exception if dimensions do not match.
template<typename Array1 , typename Array2 , typename Array3 , typename Array4 >
void	cusp::assert_same_dimensions (const Array1 &array1, const Array2 &array2, const Array3 &array3, const Array4 &array4)
	Assert dimensions of two arrays are equal. Throw invalid input exception if dimensions do not match.

Detailed Description

Algorithms for processing matrices in sparse and dense formats.

format

Function Documentation

add()

template<typename MatrixType1 , typename MatrixType2 , typename MatrixType3 >

void cusp::add	(	const MatrixType1 &	A,
		const MatrixType2 &	B,
		MatrixType3 &	C
	)

Compute the sum of two matrices.

Template Parameters

MatrixType1	Type of first matrix
MatrixType2	Type of second matrix
MatrixType3	Type of output matrix

Parameters

A	first input matrix
B	second input matrix
C	output matrix

Example

The following code snippet demonstrates how to use add.

#include <cusp/elementwise.h>
#include <cusp/array2d.h>
#include <cusp/print.h>
 
int main(void)
{
    // initialize first 2x3 matrix
    cusp::array2d<float, cusp::host_memory> A(2,3);
    A(0,0) = 10;  A(0,1) = 20;  A(0,2) = 30;
    A(1,0) = 40;  A(1,1) = 50;  A(1,2) = 60;
 
    // print A
    cusp::print(A);
 
    // initialize second 2x3 matrix
    cusp::array2d<float, cusp::host_memory> B(2,3);
    B(0,0) = 60;  B(0,1) = 50;  B(0,2) = 40;
    B(1,0) = 30;  B(1,1) = 20;  B(1,2) = 10;
 
    // print B
    cusp::print(B);
 
    // compute the sum
    cusp::array2d<float, cusp::host_memory> C;
    cusp::add(A, B, C);
 
    // print C
    cusp::print(C);
 
    return 0;
}

assert_is_valid_matrix()

template<typename MatrixType >

void cusp::assert_is_valid_matrix

(

const MatrixType &

matrix

)

Validate format of a given matrix and exit if invalid.

Template Parameters

MatrixType

matrix container

Parameters

matrix

A matrix container (e.g. csr_matrix or coo_matrix)

assert_same_dimensions() [1/3]

template<typename Array1 , typename Array2 >

void cusp::assert_same_dimensions	(	const Array1 &	array1,
		const Array2 &	array2
	)

Assert dimensions of two arrays are equal. Throw invalid input exception if dimensions do not match.

Template Parameters

Array1	First array container type
Array2	Second array container type

Parameters

array1	First array used in comparison
array2	Second array used in comparison

assert_same_dimensions() [2/3]

template<typename Array1 , typename Array2 , typename Array3 >

void cusp::assert_same_dimensions	(	const Array1 &	array1,
		const Array2 &	array2,
		const Array3 &	array3
	)

Assert dimensions of two arrays are equal. Throw invalid input exception if dimensions do not match.

Template Parameters

Array1	First array container type
Array2	Second array container type
Array3	Third array container type

Parameters

array1	First array used in comparison
array2	Second array used in comparison
array3	Third array used in comparison

assert_same_dimensions() [3/3]

template<typename Array1 , typename Array2 , typename Array3 , typename Array4 >

void cusp::assert_same_dimensions	(	const Array1 &	array1,
		const Array2 &	array2,
		const Array3 &	array3,
		const Array4 &	array4
	)

Assert dimensions of two arrays are equal. Throw invalid input exception if dimensions do not match.

Template Parameters

Array1	First array container type
Array2	Second array container type
Array3	Third array container type
Array4	Fourth array container type

Parameters

array1	First array used in comparison
array2	Second array used in comparison
array3	Third array used in comparison
array4	Fourth array used in comparison

compute_max_entries_per_row()

template<typename ArrayType >

size_t cusp::compute_max_entries_per_row

(

const ArrayType &

row_offsets

)

Compute the maximum row length of a matrix.

Template Parameters

ArrayType

Type of input row offsets

Parameters

row_offsets

row offsets of input matrix

Returns

maximum row length

Example

#include <cusp/csr_matrix.h>
#include <cusp/gallery/poisson.h>
 
#include <cusp/format_utils.h>
 
#include <iostream>
 
int main()
{
  // initialize 5x5 poisson matrix
  cusp::csr_matrix<int,float,cusp::host_memory> A;
  cusp::gallery::poisson5pt(A, 5, 5);
 
  // count the number of occupied diagonals of A
  std::cout << compute_max_entries_per_row(A.row_offsets) << std::endl;
}

compute_optimal_entries_per_row()

template<typename ArrayType >

size_t cusp::compute_optimal_entries_per_row	(	const ArrayType &	row_offsets,
		float	relative_speed = 3.0f,
		size_t	breakeven_threshold = 4096
	)

Compute the optimal number of entries per row of HYB matrix.

Template Parameters

ArrayType

Type of input row offsets

Parameters

row_offsets	row offsets of input matrix
relative_speed	Estimated performance difference between ELL and COO row storage schemes
breakeven_threshold	Threshold value separating ELL and COO row classification

Returns

optimal number of columns to store in ELL matrix

Example

#include <cusp/csr_matrix.h>
#include <cusp/gallery/poisson.h>
 
#include <cusp/format_utils.h>
 
#include <iostream>
 
int main()
{
  // initialize 5x5 poisson matrix
  cusp::csr_matrix<int,float,cusp::host_memory> A;
  cusp::gallery::poisson5pt(A, 5, 5);
 
  // count the number of occupied diagonals of A
  std::cout << compute_optimal_entries_per_row(A.row_offsets) << std::endl;
}

convert()

template<typename SourceType , typename DestinationType >

void cusp::convert	(	const SourceType &	src,
		DestinationType &	dst
	)

Convert between matrix formats.

Template Parameters

SourceType	Type of the input matrix to convert
DestinationType	Type of the output matrix to create

Parameters

src	Input matrix to convert
dst	Output matrix created by converting src to the specified format

Note

DestinationType will be resized as necessary

Example

#include <cusp/array2d.h>
#include <cusp/coo_matrix.h>
#include <cusp/print.h>
 
#include <cusp/gallery/poisson.h>
 
// include cusp convert header file
#include <cusp/convert.h>
 
int main()
{
  // create an empty sparse matrix structure
  cusp::coo_matrix<int,float,cusp::host_memory> A;
 
  // create 2D Poisson problem
  cusp::gallery::poisson5pt(A, 4, 4);
 
  // create an empty array2d structure
  cusp::array2d<float,cusp::host_memory> B;
 
  // convert coo_matrix to array2d
  cusp::convert(A, B);
 
  // print the contents of B
  cusp::print(B);
}

See also

copy

copy()

template<typename SourceType , typename DestinationType >

void cusp::copy	(	const SourceType &	src,
		DestinationType &	dst
	)

Copy one array or matrix to another.

Template Parameters

SourceType	Type of the input matrix to copy
DestinationType	Type of the output matrix

Parameters

src	Input matrix to copy
dst	Output matrix created by copying src to dst

Note

SourceType and DestinationType must have the same format

DestinationType will be resized as necessary

Example

#include <cusp/array1d.h>
#include <cusp/print.h>
 
// include cusp copy header file
#include <cusp/copy.h>
 
int main()
{
  // Allocate a array of size 10
  cusp::array1d<int,cusp::host_memory> a(10);
 
  // Create a random array with 10 entries
  cusp::random_array<int> rand(10);
 
  // Copy random values from rand into a
  cusp::copy(rand, a);
 
  // print the contents of a
  cusp::print(a);
}

See also

convert

count_diagonals()

template<typename ArrayType1 , typename ArrayType2 >

size_t cusp::count_diagonals	(	const size_t	num_rows,
		const size_t	num_cols,
		const ArrayType1 &	row_indices,
		const ArrayType2 &	column_indices
	)

Count the number of occupied diagonals in the input matrix.

Template Parameters

ArrayType1	Type of input row indices
ArrayType2	Type of input column indices

Parameters

num_rows	Number of rows.
num_cols	Number of columns.
row_indices	row indices of input matrix
column_indices	column indices of input matrix

Returns

number of occupied diagonals

Example

#include <cusp/coo_matrix.h>
#include <cusp/gallery/poisson.h>
 
#include <cusp/format_utils.h>
 
#include <iostream>
 
int main()
{
  // initialize 5x5 poisson matrix
  cusp::coo_matrix<int,float,cusp::host_memory> A;
  cusp::gallery::poisson5pt(A, 5, 5);
 
  // count the number of occupied diagonals of A
  std::cout << count_diagonals(A.num_rows,
                               A.num_cols,
                               A.row_indices,
                               A.column_indices) << std::endl;
}

counting_sort()

template<typename ArrayType >

void cusp::counting_sort	(	ArrayType &	keys,
		typename ArrayType::value_type	min,
		typename ArrayType::value_type	max
	)

Use counting sort to order an array.

Template Parameters

ArrayType

Type of input array

Parameters

keys	input of keys to sort
min	minimum key in input array
max	maximum key in input array

Example

The following code snippet demonstrates how to use counting_sort.

#include <cusp/array1d.h>
#include <cusp/print.h>
#include <cusp/sort.h>
 
#include <thrust/transform.h>
 
struct mod
{
    int operator()(int i)
    {
        // transform i to range [0,5)
        return (i >= 0 ? i : -i) % 5;
    }
};
 
int main(void)
{
    // create random array with 10 elements
    cusp::random_array<int> rand(10);
    // initialize v to random array
    cusp::array1d<int,cusp::host_memory> v(rand);
    // transform entries to interval [0,5)
    thrust::transform(v.begin(), v.end(), v.begin(), mod());
    // print array
    cusp::print(v);
    // sort
    cusp::counting_sort(v, 0, 5);
    // print the sorted array
    cusp::print(v);
 
    return 0;
}

counting_sort_by_key()

template<typename ArrayType1 , typename ArrayType2 >

void cusp::counting_sort_by_key	(	ArrayType1 &	keys,
		ArrayType2 &	vals,
		typename ArrayType1::value_type	min,
		typename ArrayType1::value_type	max
	)

Use counting sort to order an array and permute an array of values.

Template Parameters

ArrayType1	Type of keys array
ArrayType2	Type of values array

Parameters

keys	input array of keys to sort
vals	input array of values to permute
min	minimum key in keys array
max	maximum key in keys array

Example

The following code snippet demonstrates how to use counting_sort_by_key.

#include <cusp/array1d.h>
#include <cusp/print.h>
#include <cusp/sort.h>
 
#include <thrust/transform.h>
 
struct mod
{
    int operator()(int i)
    {
        // transform i to range [0,5)
        return (i >= 0 ? i : -i) % 5;
    }
};
 
int main(void)
{
    // create random array with 10 elements
    cusp::random_array<int> rand1(10,0);
    cusp::random_array<int> rand2(10,3);
    // initialize v to random array
    cusp::array1d<int,cusp::host_memory> v1(rand1);
    cusp::array1d<int,cusp::host_memory> v2(rand2);
    // transform entries to interval [0,5)
    thrust::transform(v1.begin(), v1.end(), v1.begin(), mod());
    // print array
    cusp::print(v1);
    cusp::print(v2);
    // sort
    cusp::counting_sort_by_key(v1, v2, 0, 5);
    // print the sorted array
    cusp::print(v1);
    cusp::print(v2);
 
    return 0;
}

elementwise()

template<typename MatrixType1 , typename MatrixType2 , typename MatrixType3 , typename BinaryFunction >

void cusp::elementwise	(	const MatrixType1 &	A,
		const MatrixType2 &	B,
		MatrixType3 &	C,
		BinaryFunction	op
	)

Perform transform operation on two matrices.

Template Parameters

MatrixType1	Type of first matrix
MatrixType2	Type of second matrix
MatrixType3	Type of output matrix
BinaryFunction	Type of binary transform to apply

Parameters

A	first input matrix
B	second input matrix
C	output matrix
op	binary transform to apply

Example

The following code snippet demonstrates how to use elementwise.

#include <cusp/elementwise.h>
#include <cusp/array2d.h>
#include <cusp/print.h>
 
int main(void)
{
    // initialize first 2x3 matrix
    cusp::array2d<float, cusp::host_memory> A(2,3);
    A(0,0) = 10;  A(0,1) = 20;  A(0,2) = 30;
    A(1,0) = 40;  A(1,1) = 50;  A(1,2) = 60;
 
    // print A
    cusp::print(A);
 
    // initialize second 2x3 matrix
    cusp::array2d<float, cusp::host_memory> B(2,3);
    B(0,0) = 60;  B(0,1) = 50;  B(0,2) = 40;
    B(1,0) = 30;  B(1,1) = 20;  B(1,2) = 10;
 
    // print B
    cusp::print(B);
 
    // compute the sum
    cusp::array2d<float, cusp::host_memory> C;
    cusp::elementwise(A, B, C, ::cuda::std::plus<int>());
 
    // print C
    cusp::print(C);
 
    return 0;
}

extract_diagonal()

template<typename MatrixType , typename ArrayType >

void cusp::extract_diagonal	(	const MatrixType &	A,
		ArrayType &	output
	)

Extract the main diagonal of a matrix.

Template Parameters

MatrixType	Type of input matrix
ArrayType	Type of input diagonal array

Parameters

A	The input matrix
output	On return contains the main diagonal of A with zeros inserted for missing entries.

Example

// include cusp array1d header file
#include <cusp/array1d.h>
#include <cusp/coo_matrix.h>
#include <cusp/print.h>
#include <cusp/gallery/poisson.h>
 
#include <cusp/format_utils.h>
 
int main()
{
  // initialize 5x5 poisson matrix
  cusp::coo_matrix<int,float,cusp::host_memory> A;
  cusp::gallery::poisson5pt(A, 5, 5);
 
  // allocate array to hold diagonal entries
  cusp::array1d<float, cusp::host_memory> diagonal(A.num_rows);
 
  // extract diagonal of A
  cusp::extract_diagonal(A, diagonal);
 
  // print diagonal entries
  cusp::print(diagonal);
}

generalized_spgemm()

template<typename LinearOperator , typename MatrixOrVector1 , typename MatrixOrVector2 , typename UnaryFunction , typename BinaryFunction1 , typename BinaryFunction2 >

void cusp::generalized_spgemm	(	const LinearOperator &	A,
		const MatrixOrVector1 &	B,
		MatrixOrVector2 &	C,
		UnaryFunction	initialize,
		BinaryFunction1	combine,
		BinaryFunction2	reduce
	)

Implements generalized matrix-matrix multiplication.

Overview

generalized multiply can be used with dense and sparse matrices, and user-defined linear_operator objects. This function compute nonzeros only for the entries present in the output matrix. Entries that are not specified in the output matrix are disregarded (annihilated). This specification significantly reduces the computational work corresponding to computing the sparsity of the output and performing unnecessary (combine, reduce) operations.

Template Parameters

LinearOperator	Type of matrix
MatrixOrVector1	Type of second matrix
MatrixOrVector2	Type of output matrix
UnaryFunction	Type of unary function to initialize RHS
BinaryFunction1	Type of binary function to combine entries
BinaryFunction2	Type of binary function to reduce entries

Parameters

A	first input matrix
B	second input matrix
C	output matrix
initialize	unary function used to initialize the output
combine	binary function used to combine entries
reduce	binary function used to reduce entries

Example

The following code snippet demonstrates how to use generalized_spgemm to compute a matrix-matrix product.

#include <cusp/coo_matrix.h>
#include <cusp/functional.h>
#include <cusp/multiply.h>
#include <cusp/print.h>
 
#include <cusp/gallery/poisson.h>
 
int main(void)
{
    // define multiply functors
    ::cuda::std::identity   identity;
    cusp::constant_functor<float> zero;
    ::cuda::std::multiplies<float> combine;
    ::cuda::std::plus<float>       reduce;
 
    // initialize matrix
    cusp::coo_matrix<int,float,cusp::host_memory> A;
    cusp::gallery::poisson5pt(A, 3, 3);
 
    // allocate output matrices and initialize output nonzeros
    cusp::coo_matrix<int, float, cusp::host_memory> B(A);
    cusp::coo_matrix<int, float, cusp::host_memory> C(A);
 
    // compute B = A * A
    cusp::generalized_spgemm(A, A, B, zero, combine, reduce);
    // compute C += A * A
    cusp::generalized_spgemm(A, A, C, identity, combine, reduce);
 
    // print output matrices
    cusp::print(B);
    cusp::print(C);
 
    return 0;
}

generalized_spmv()

template<typename LinearOperator , typename Vector1 , typename Vector2 , typename Vector3 , typename BinaryFunction1 , typename BinaryFunction2 >

void cusp::generalized_spmv	(	const LinearOperator &	A,
		const Vector1 &	x,
		const Vector2 &	y,
		Vector3 &	z,
		BinaryFunction1	combine,
		BinaryFunction2	reduce
	)

Implements generalized matrix-vector multiplication.

Overview

generalized multiply can be used with dense and sparse matrices, and user-defined linear_operator objects.

Template Parameters

LinearOperator	Type of first matrix
Vector1	Type of second input vector
Vector2	Type of third input vector
Vector3	Type of output vector

Parameters

A	input matrix
x	input vector
y	input vector
z	output vector
combine	binary function used to combine entries
reduce	binary function used to reduce entries

Example

The following code snippet demonstrates how to use multiply to compute a matrix-vector product.

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/multiply.h>
#include <cusp/print.h>
 
int main(void)
{
    // define multiply functors
    ::cuda::std::multiplies<float>          combine;
    ::cuda::std::plus<float>                reduce;
 
    // initialize matrix
    cusp::array2d<float, cusp::host_memory> A(2,2);
    A(0,0) = 10;  A(0,1) = 20;
    A(1,0) = 40;  A(1,1) = 50;
 
    // initialize input vector
    cusp::array1d<float, cusp::host_memory> x(2);
    x[0] = 1;
    x[1] = 2;
 
    // initial RHS filled with 2's
    cusp::constant_array<float> y(2, 2);
 
    // allocate output vector
    cusp::array1d<float, cusp::host_memory> z(2);
 
    // compute z = y + (A * x)
    cusp::generalized_spmv(A, x, y, z, combine, reduce);
 
    // print z
    cusp::print(z);
 
    return 0;
}

indices_to_offsets()

template<typename IndexArray , typename OffsetArray >

void cusp::indices_to_offsets	(	const IndexArray &	indices,
		OffsetArray &	offsets
	)

Compress COO row indices to CSR row offsets.

Template Parameters

IndexArray	Type of input row indices
OffsetType	Type of input row offsets

Parameters

indices	The input row indices
offsets	The output row offsets

Example

// include cusp array1d header file
#include <cusp/array1d.h>
#include <cusp/format_utils.h>
#include <cusp/print.h>
 
int main()
{
  // Allocate a array of size 5
  cusp::array1d<int,cusp::host_memory> a(5);
 
  // Set the elements to [0, 0, 1, 1, 1]
  a[0] = 0; a[1] = 0; a[2] = 1; a[3] = 1; a[4] = 1;
 
  // Allocate a seceond array equal to the last entry in
  // a (number of rows) + 2 (zero entry + last entry)
  cusp::array1d<int,cusp::host_memory> b(a.back() + 2);
 
  // compute the expanded entries
  cusp::indices_to_offsets(a, b);
 
  // print [0, 2, 5]
  cusp::print(b);
}

is_valid_matrix() [1/2]

template<typename MatrixType >

bool cusp::is_valid_matrix

(

const MatrixType &

matrix

)

Validate format of a given matrix.

Template Parameters

MatrixType

matrix container

Parameters

matrix

A matrix container (e.g. csr_matrix or coo_matrix)

Returns

true if format is valid otherwise false

is_valid_matrix() [2/2]

template<typename MatrixType , typename OutputStream >

bool cusp::is_valid_matrix	(	const MatrixType &	matrix,
		OutputStream &	ostream
	)

Validate format of a given matrix.

Template Parameters

MatrixType	matrix container
OutputStream	stream type

Parameters

matrix	A matrix container (e.g. csr_matrix or coo_matrix)
ostream	Stream to which the validate stream should print

Returns

true if format is valid otherwise false

multiply() [1/2]

template<typename LinearOperator , typename MatrixOrVector1 , typename MatrixOrVector2 >

void cusp::multiply	(	const LinearOperator &	A,
		const MatrixOrVector1 &	B,
		MatrixOrVector2 &	C
	)

Implements matrix-matrix and matrix-vector multiplication.

Overview

multiply can be used with dense matrices, sparse matrices, and user-defined linear_operator objects.

Template Parameters

LinearOperator	Type of first matrix
MatrixOrVector1	Type of second matrix or vector
MatrixOrVector2	Type of output matrix or vector

Parameters

A	input matrix
B	input matrix or vector
C	output matrix or vector

Example

The following code snippet demonstrates how to use multiply to compute a matrix-vector product.

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/multiply.h>
#include <cusp/print.h>
 
int main(void)
{
    // initialize matrix
    cusp::array2d<float, cusp::host_memory> A(2,2);
    A(0,0) = 10;  A(0,1) = 20;
    A(1,0) = 40;  A(1,1) = 50;
 
    // initialize input vector
    cusp::array1d<float, cusp::host_memory> x(2);
    x[0] = 1;
    x[1] = 2;
 
    // allocate output vector
    cusp::array1d<float, cusp::host_memory> y(2);
 
    // compute y = A * x
    cusp::multiply(A, x, y);
 
    // print y
    cusp::print(y);
 
    return 0;
}

multiply() [2/2]

template<typename LinearOperator , typename MatrixOrVector1 , typename MatrixOrVector2 , typename UnaryFunction , typename BinaryFunction1 , typename BinaryFunction2 >

void cusp::multiply	(	const LinearOperator &	A,
		const MatrixOrVector1 &	B,
		MatrixOrVector2 &	C,
		UnaryFunction	initialize,
		BinaryFunction1	combine,
		BinaryFunction2	reduce
	)

Implements matrix-vector multiplication with custom combine and reduce functionality.

Overview

multiply can be used with dense matrices, sparse matrices, and user-defined linear_operator objects.

Template Parameters

LinearOperator	Type of matrix
MatrixOrVector1	Type of second vector
MatrixOrVector2	Type of output vector
UnaryFunction	Type of unary function to initialize RHS
BinaryFunction1	Type of binary function to combine entries
BinaryFunction2	Type of binary function to reduce entries

Parameters

A	input matrix
B	input vector
C	output vector
initialize	unary function used to initialize the output
combine	binary function used to combine entries
reduce	binary function used to reduce entries

Example

The following code snippet demonstrates how to use multiply to compute a matrix-vector product.

#include <cusp/array1d.h>
#include <cusp/array2d.h>
#include <cusp/functional.h>
#include <cusp/multiply.h>
#include <cusp/print.h>
 
int main(void)
{
    // define multiply functors
    cusp::constant_functor<float> initialize;
    ::cuda::std::multiplies<float> combine;
    ::cuda::std::plus<float>       reduce;
 
    // initialize matrix
    cusp::array2d<float, cusp::host_memory> A(2,2);
    A(0,0) = 10;  A(0,1) = 20;
    A(1,0) = 40;  A(1,1) = 50;
 
    // initialize input vector
    cusp::array1d<float, cusp::host_memory> x(2);
    x[0] = 1;
    x[1] = 2;
 
    // allocate output vector
    cusp::array1d<float, cusp::host_memory> y(2);
 
    // compute y = A * x
    cusp::multiply(A, x, y, initialize, combine, reduce);
 
    // print y
    cusp::print(y);
 
    return 0;
}

offsets_to_indices()

template<typename OffsetArray , typename IndexArray >

void cusp::offsets_to_indices	(	const OffsetArray &	offsets,
		IndexArray &	indices
	)

Expand CSR row offsets to COO row indices.

Template Parameters

OffsetType	Type of input row offsets
IndexArray	Type of input row indices

Parameters

offsets	The input row offsets
indices	The output row indices

Example

// include cusp array1d header file
#include <cusp/array1d.h>
#include <cusp/format_utils.h>
 
int main()
{
  // Allocate a array of size 3
  cusp::array1d<int,cusp::host_memory> a(3);
 
  // Set the elements to [0, 2, 5]
  a[0] = 0; a[1] = 2; a[2] = 5;
 
  // Allocate a seceond array equal to the total number
  // of expanded entries from vector a (last entry in the vector)
  cusp::array1d<int,cusp::host_memory> b(a.back());
 
  // compute the expanded entries
  cusp::offsets_to_indices(a, b);
 
  // print [0, 0, 1, 1, 1]
  cusp::print(b);
}

sort_by_row()

template<typename ArrayType1 , typename ArrayType2 , typename ArrayType3 >

void cusp::sort_by_row	(	ArrayType1 &	row_indices,
		ArrayType2 &	column_indices,
		ArrayType3 &	values,
		typename ArrayType1::value_type	min_row = 0,
		typename ArrayType1::value_type	max_row = 0
	)

Sort matrix indices by row.

Template Parameters

ArrayType1	Type of input matrix row indices
ArrayType2	Type of input matrix column indices
ArrayType3	Type of input matrix values

Parameters

row_indices	input matrix row indices
column_indices	input matrix column indices
values	input matrix values
min_row	minimum row index
max_row	maximum row index

Example

The following code snippet demonstrates how to use sort_by_row.

#include <cusp/coo_matrix.h>
#include <cusp/print.h>
#include <cusp/sort.h>
 
int main(void)
{
    // allocate storage for (4,3) matrix with 6 nonzeros
    cusp::coo_matrix<int,float,cusp::host_memory> A(4,3,6);
 
    // initialize matrix entries on host
    A.row_indices[0] = 3; A.column_indices[0] = 0; A.values[0] = 10;
    A.row_indices[1] = 3; A.column_indices[1] = 2; A.values[1] = 20;
    A.row_indices[2] = 2; A.column_indices[2] = 0; A.values[2] = 30;
    A.row_indices[3] = 2; A.column_indices[3] = 2; A.values[3] = 40;
    A.row_indices[4] = 0; A.column_indices[4] = 1; A.values[4] = 50;
    A.row_indices[5] = 0; A.column_indices[5] = 2; A.values[5] = 60;
 
    // sort A by row
    cusp::sort_by_row(A.row_indices, A.column_indices, A.values);
 
    // print the sorted matrix
    cusp::print(A);
 
    return 0;
}

sort_by_row_and_column()

template<typename ArrayType1 , typename ArrayType2 , typename ArrayType3 >

void cusp::sort_by_row_and_column	(	ArrayType1 &	row_indices,
		ArrayType2 &	column_indices,
		ArrayType3 &	values,
		typename ArrayType1::value_type	min_row = 0,
		typename ArrayType1::value_type	max_row = 0,
		typename ArrayType2::value_type	min_col = 0,
		typename ArrayType2::value_type	max_col = 0
	)

Sort matrix indices by row and column.

Template Parameters

ArrayType1	Type of input matrix row indices
ArrayType2	Type of input matrix column indices
ArrayType3	Type of input matrix values

Parameters

row_indices	input matrix row indices
column_indices	input matrix column indices
values	input matrix values
min_row	minimum row index
max_row	maximum row index
min_col	minimum column index
max_col	maximum column index

Example

The following code snippet demonstrates how to use sort_by_row_and_column.

#include <cusp/coo_matrix.h>
#include <cusp/print.h>
#include <cusp/sort.h>
 
int main(void)
{
    // allocate storage for (4,3) matrix with 6 nonzeros
    cusp::coo_matrix<int,float,cusp::host_memory> A(4,3,6);
 
    // initialize matrix entries on host
    A.row_indices[0] = 3; A.column_indices[0] = 2; A.values[0] = 10;
    A.row_indices[1] = 3; A.column_indices[1] = 0; A.values[1] = 20;
    A.row_indices[2] = 2; A.column_indices[2] = 0; A.values[2] = 30;
    A.row_indices[3] = 2; A.column_indices[3] = 2; A.values[3] = 40;
    A.row_indices[4] = 0; A.column_indices[4] = 2; A.values[4] = 50;
    A.row_indices[5] = 0; A.column_indices[5] = 1; A.values[5] = 60;
 
    // sort A by row
    cusp::sort_by_row_and_column(A.row_indices, A.column_indices, A.values);
 
    // print the sorted matrix
    cusp::print(A);
 
    return 0;
}

subtract()

template<typename MatrixType1 , typename MatrixType2 , typename MatrixType3 >

void cusp::subtract	(	const MatrixType1 &	A,
		const MatrixType2 &	B,
		MatrixType3 &	C
	)

Compute the difference of two matrices.

Template Parameters

MatrixType1	Type of first matrix
MatrixType2	Type of second matrix
MatrixType3	Type of output matrix

Parameters

A	first input matrix
B	second input matrix
C	output matrix

Example

The following code snippet demonstrates how to use subtract.

#include <cusp/elementwise.h>
#include <cusp/array2d.h>
#include <cusp/print.h>
 
int main(void)
{
    // initialize first 2x3 matrix
    cusp::array2d<float, cusp::host_memory> A(2,3);
    A(0,0) = 10;  A(0,1) = 20;  A(0,2) = 30;
    A(1,0) = 40;  A(1,1) = 50;  A(1,2) = 60;
 
    // print A
    cusp::print(A);
 
    // initialize second 2x3 matrix
    cusp::array2d<float, cusp::host_memory> B(2,3);
    B(0,0) = 60;  B(0,1) = 50;  B(0,2) = 40;
    B(1,0) = 30;  B(1,1) = 20;  B(1,2) = 10;
 
    // print B
    cusp::print(B);
 
    // compute the subtract
    cusp::array2d<float, cusp::host_memory> C;
    cusp::subtract(A, B, C);
 
    // print C
    cusp::print(C);
 
    return 0;
}

transpose()

template<typename MatrixType1 , typename MatrixType2 >

void cusp::transpose	(	const MatrixType1 &	A,
		MatrixType2 &	At
	)

Transpose a matrix.

Template Parameters

MatrixType1	Type of input matrix to transpose
MatrixType2	Type of output matrix

Parameters

A	input matrix
At	output matrix (transpose of A)

Example

The following code snippet demonstrates how to use transpose.

#include <cusp/transpose.h>
#include <cusp/array2d.h>
#include <cusp/print.h>
 
int main(void)
{
    // initialize a 2x3 matrix
    cusp::array2d<float, cusp::host_memory> A(2,3);
    A(0,0) = 10;  A(0,1) = 20;  A(0,2) = 30;
    A(1,0) = 40;  A(1,1) = 50;  A(1,2) = 60;
 
    // print A
    cusp::print(A);
 
    // compute the transpose
    cusp::array2d<float, cusp::host_memory> At;
    cusp::transpose(A, At);
 
    // print A^T
    cusp::print(At);
 
    return 0;
}