SCPPNT: SCPPNT Namespace Reference


Classes
class	Matrix
	Simple concrete dense matrix class for numerical computation. More...
class	Index1D
	Represents a consecutive index range (a 1-offset lower and upper index). More...
class	Region2D_iterator
	Foward iterator over all elements of a matrix region. More...
class	Region2D
	A rectangular subset of a matrix. More...
class	const_Region2D
	Constant version of Region2D. More...
class	Exception
	General exception class from which classes indicating more specific types are errors are derived. More...
class	LogicError
	Errors in program logic that should be caught during program debugging. More...
class	InvalidArgument
	A function argument is invalid. More...
class	RuntimeError
	General run-time errors. Classes for specific types of run-time errors are derived from this class. More...
class	BoundsError
	Attempt to access element outside vector or matrix bounds. More...
class	BadDimension
	Bad dimension of vector or matrix argument of function. More...
class	slice_pointer_base
	Acts like T* but allows consecutive elements to be stored non-consecutively in memory. More...
class	Standardize
	Standardize a set of values so they sum to a particular value. More...
class	Sum
	Compute the sum of a sequence of values. More...
class	Mean
	Compute the mean of a sequence of values. More...
class	Transpose_View
	Allows use of the transpose of a matrix without creating a new matrix. More...
class	LowerTriangularView
	Lower triangular view of a matrix. More...
class	UnitLowerTriangularView
	Lower triangular view of a matrix with a unit diagonal. More...
class	UpperTriangularView
	Upper triangular view of a matrix. More...
class	UnitUpperTriangularView
	Upper triangular view of a matrix with a unit diagonal. More...
class	Vector
	Simple concrete vector class for numerical computation. More...
Typedefs
typedef SCPPNT_SUBSCRIPT_TYPE	Subscript
	Subscript index type.
Functions
template<class SPDMatrix, class SymmMatrix>
int	Cholesky_upper_factorization (SPDMatrix &A, SymmMatrix &L)
	Compute Cholesky factorization.
template<class T>
std::ostream &	operator<< (std::ostream &s, const Matrix< T > &A)
	Write a matrix to an output stream.
template<class T>
std::istream &	operator>> (std::istream &s, Matrix< T > &A)
	Read elements of a matrix to an input stream.
template<class T, class M>
Matrix< T >	operator+ (const Matrix< T > &A, const M &B)
	Add two matrices and return new matrix giving sum.
template<class T, class M>
Matrix< T >	operator- (const Matrix< T > &A, const M &B)
	Subtract matrix B from matrix A and return matrix giving difference.
template<class T, class M>
Matrix< T >	operator * (const Matrix< T > &A, const M &B)
	Multiplication operator returning a new matrix containing the matrix product A*B.
template<class T>
Vector< T >	operator * (const Matrix< T > &A, const Vector< T > &x)
	Return a new vector containing the product of the matrix A and the vector x.
Index1D	operator+ (const Index1D &D, Subscript i)
	Increment index range by i (increase both lower and upper limits by i).
Index1D	operator+ (Subscript i, const Index1D &D)
	Increment index range by i (increase both lower and upper limits by i).
Index1D	operator- (Index1D &D, Subscript i)
	Decrement index range by i (decrease both lower and upper limits by i).
Index1D	operator- (Subscript i, Index1D &D)
	Decrement index range by i (decrease both lower and upper limits by i).
template<class MaTRiX, class VecToRSubscript>
int	LU_factor (MaTRiX &A, VecToRSubscript &indx)
	Right-looking LU factorization algorithm (unblocked).
template<class MaTRiX, class VecToR, class VecToRSubscripts>
int	LU_solve (const MaTRiX &A, const VecToRSubscripts &indx, VecToR &b)
	Solve system of linear equations using LU factorization computed from LU_factor.
template<class M, class FUNC>
void	apply_rows (M &matrix, FUNC &f)
	Apply a function object to each row of a matrix.
template<class M, class FUNC>
void	apply_columns (M &matrix, FUNC &f)
	Apply a function object to each column of a matrix.
template<class M, class V, class FUNC>
V	over_columns (M &matrix, FUNC &f)
template<class M, class V, class FUNC>
V	over_rows (M &matrix, FUNC &f)
template<class MAT>
MAT	transpose (const MAT &A)
	Return a new matrix containing the transpose of A.
template<class MA, class MB>
void	matcopy (MA &A, const MB &B)
	Copy all elements of matrix B to matrix A (B must be the same size as A).
template<class MA, class MB, class MC>
void	matadd (const MA &A, const MB &B, MC &C)
	Add matrices A and B and place sum in C.
template<class MA, class MB>
void	matadd (MA &A, const MB &B)
	Add matrices A and B and replace A with result.
template<class MA, class MB, class MC>
void	matsub (const MA &A, const MB &B, MC &C)
	Substract B from A and place difference in C.
template<class MA, class MB>
void	matsub (MA &A, const MB &B)
	Subtract B from A and replace A with result.
template<class MR, class MA, class MB>
MR	mult_element (const MA &A, const MB &B)
	Return new matrix containing the products of the corresponding elements of A and B (element-by-element product).
template<class MR, class MA, class MB>
MR	matmult (const MA &A, const MB &B)
	Return a new matrix containing the matrix product A*B.
template<class MC, class MA, class MB>
void	matmult (MC &C, const MA &A, const MB &B)
	Calculate matrix product A*B and put result in C.
template<class MA, class MB>
void	matmult_assign (MA &A, const MB &B)
	Matrix multiplication of AB where A is overwritten to store result (requires AB and A to be same size).
template<class MAT, class IT1, class IT2>
void	MatTimesVec (const MAT &A, IT1 begin, IT2 result)
	Multiply a matrix times a vector using iterators to access the vector elements.
template<class M, class V>
V	matrix_times_vector (const M &A, const V &x)
	Multiply a matrix times a vector and return the result in a new vector.
template<class MaTRiX, class Vector>
int	QR_factor (MaTRiX &A, Vector &C, Vector &D)
	Classical QR factorization, based on Stewart[1973].
template<class MaTRiX, class Vector>
int	R_solve (const MaTRiX &A, Vector &D, Vector &b)
	Modified form of upper triangular solve, except that the main diagonal of R (upper portion of A) is in D.
template<class MaTRiX, class Vector>
int	QR_solve (const MaTRiX &A, const Vector &c, Vector &d, Vector &b)
	Solve Rx = Q'b, where A, c, and d are output from QR_factor.
template<class Array2D>
std::ostream &	operator<< (std::ostream &s, const Region2D< Array2D > &A)
	Write a matrix to an output stream.
template<class Array2D, class M>
M	operator+ (const Region2D< Array2D > &A, const M &B)
	Add two matrices and return new matrix giving sum.
template<class Array2D, class M>
M	operator- (const Region2D< Array2D > &A, const M &B)
	Subtract matrix B from matrix A and return matrix giving difference.
template<class Array2D, class M>
M	operator * (const Region2D< Array2D > &A, const M &B)
	Multiplication operator returning a new matrix containing the matrix product A*B.
template<class T, class Array2D>
Vector< T >	operator * (const Region2D< Array2D > &A, const Vector< T > &x)
	Return a new vector containing the product of the matrix A and the vector x.
template<class Array2D>
std::ostream &	operator<< (std::ostream &s, const const_Region2D< Array2D > &A)
	Write a matrix to an output stream.
template<class Array2D, class M>
M	operator+ (const const_Region2D< Array2D > &A, const M &B)
	Add two matrices and return new matrix giving sum.
template<class Array2D, class M>
M	operator- (const const_Region2D< Array2D > &A, const M &B)
	Subtract matrix B from matrix A and return matrix giving difference.
template<class Array2D, class M>
M	operator * (const const_Region2D< Array2D > &A, const M &B)
	Multiplication operator returning a new matrix containing the matrix product A*B.
template<class T, class Array2D>
Vector< T >	operator * (const const_Region2D< Array2D > &A, const Vector< T > &x)
	Return a new vector containing the product of the matrix A and the vector x.
template<class T>
T	abs (T x)
	Returns absolute value of argument.
template<class T>
T	min (T a, T b)
	Returns minimum of two scalars.
template<class T>
T	max (T a, T b)
	Returns maximum of two scalars.
template<class T>
T	sign (T a)
	Returns 1 if argument is positive, otherwise returns -1.
template<class Array2D>
std::ostream &	operator<< (std::ostream &s, const Transpose_View< Array2D > &A)
	Write a matrix to an output stream.
template<class Array2D, class M>
M	operator+ (const Transpose_View< Array2D > &A, const M &B)
	Add two matrices and return new matrix giving sum.
template<class Array2D, class M>
M	operator- (const Transpose_View< Array2D > &A, const M &B)
	Subtract matrix B from matrix A and return matrix giving difference.
template<class Array2D, class M>
M	operator * (const Transpose_View< Array2D > &A, const M &B)
	Multiplication operator returning a new matrix containing the matrix product A*B.
template<class Array2D>
Vector< typename Array2D::value_type >	operator * (const Transpose_View< Array2D > &A, const Vector< typename Array2D::value_type > &x)
	Return a new vector containing the product of the matrix A and the vector x.
template<class MaTRiX, class VecToR>
VecToR	matmult (LowerTriangularView< MaTRiX > &A, VecToR &x)
	Multiply lower triangular view times a vector.
template<class MaTRiX, class VecToR>
VecToR	operator * (LowerTriangularView< MaTRiX > &A, VecToR &x)
	Multiplication operator for product of lower triangular view and a vector.
template<class MaTRiX>
LowerTriangularView< MaTRiX >	Lower_triangular_view (MaTRiX &A)
	Adaptor to create a lower trangular view of a matrix.
template<class MaTRiX>
UnitLowerTriangularView< MaTRiX >	Unit_lower_triangular_view (MaTRiX &A)
	Adaptor to create a unit lower trangular view of a matrix.
template<class MaTRiX, class VecToR>
VecToR	matmult (UnitLowerTriangularView< MaTRiX > &A, VecToR &x)
	Multiply a lower triangular view times a vector.
template<class MaTRiX, class VecToR>
VecToR	operator * (UnitLowerTriangularView< MaTRiX > &A, VecToR &x)
	Multiplication operator giving the product of a lower triangular view and a vector,.
template<class MaTRiX>
std::ostream &	operator<< (std::ostream &s, const LowerTriangularView< MaTRiX > &A)
	Assign values of a lower triangular view from a stream.
template<class MaTRiX>
std::ostream &	operator<< (std::ostream &s, const UnitLowerTriangularView< MaTRiX > &A)
	Write values of a lower triangular view to a stream.
template<class MaTRiX, class VecToR>
VecToR	matmult (UpperTriangularView< MaTRiX > &A, VecToR &x)
	Multiply an upper triangular view times a vector.
template<class MaTRiX, class VecToR>
VecToR	operator * (UpperTriangularView< MaTRiX > &A, VecToR &x)
	Multiplication operator giving project of an upper triangular view and a vector.
template<class MaTRiX>
UpperTriangularView< MaTRiX >	Upper_triangular_view (MaTRiX &A)
	Adapter to create upper triangular view of matrix.
template<class MaTRiX, class VecToR>
VecToR	matmult (UnitUpperTriangularView< MaTRiX > &A, VecToR &x)
	Multiply a unit upper triangular view times a vector.
template<class MaTRiX, class VecToR>
VecToR	operator * (UnitUpperTriangularView< MaTRiX > &A, VecToR &x)
	Multiplication operator giving the project of a unit upper triangular view and a vector.
template<class MaTRiX>
std::ostream &	operator<< (std::ostream &s, UpperTriangularView< MaTRiX > &A)
	Write the contents of an upper triangular view to a stream.
template<class MaTRiX>
std::ostream &	operator<< (std::ostream &s, UnitUpperTriangularView< MaTRiX > &A)
	Write the contents of a unit upper triangular view to a stream.
template<class MaTriX, class VecToR>
VecToR	Lower_triangular_solve (const MaTriX &A, const VecToR &b)
	Solve linear system Ax=b using lower triangular portion of A including diagonal.
template<class MaTriX, class VecToR>
VecToR	Unit_lower_triangular_solve (const MaTriX &A, const VecToR &b)
	Solve linear system Ax=b using lower triangular portion of A assuming unit diagonal.
template<class MaTriX, class VecToR>
VecToR	linear_solve (const LowerTriangularView< MaTriX > &A, const VecToR &b)
	Solve linear system using lower triangular view A.
template<class MaTriX, class VecToR>
VecToR	linear_solve (const UnitLowerTriangularView< MaTriX > &A, const VecToR &b)
	Solve linear system using unit lower triangular view A.
template<class MaTriX, class VecToR>
VecToR	Upper_triangular_solve (const MaTriX &A, const VecToR &b)
	Solve linear system Ax=b using upper triangular portion of A including diagonal.
template<class MaTriX, class VecToR>
VecToR	Unit_upper_triangular_solve (const MaTriX &A, const VecToR &b)
	Solve linear system Ax=b using upper triangular portion of A assuming unit diagonal.
template<class MaTriX, class VecToR>
VecToR	linear_solve (const UpperTriangularView< MaTriX > &A, const VecToR &b)
	Solve linear system using upper triangular view A.
template<class MaTriX, class VecToR>
VecToR	linear_solve (const UnitUpperTriangularView< MaTriX > &A, const VecToR &b)
	Solve linear system using unit upper triangular view A.
template<class T>
std::ostream &	operator<< (std::ostream &s, const Vector< T > &A)
	Write a vector to an output stream.
template<class T>
std::istream &	operator>> (std::istream &s, Vector< T > &A)
	Read elements of a vector from an input stream.
template<class T>
Vector< T >	operator+ (const Vector< T > &lhs, const Vector< T > &rhs)
	Add a vector to a vector and return a new vector containing the sum.
template<class T>
Vector< T >	operator+ (const Vector< T > &lhs, const T value)
	Add add a scalar to each element of a vector and return a new vector containing the sum.
template<class T>
Vector< T >	operator- (const Vector< T > &lhs, const Vector< T > &rhs)
	Subtract the vector rhs from the vector lhs and return a new vector containing the difference.
template<class T>
Vector< T >	operator- (const Vector< T > &lhs, const T value)
	Subtract scalar from each element of a vector a return a new vector containing the difference.
template<class T>
Vector< T >	operator * (const Vector< T > &lhs, const Vector< T > &rhs)
	Multiply corresponding elements of two vectors and return a new vector containing the product.
template<class T>
Vector< T >	operator * (const Vector< T > &lhs, const T value)
	Multiply each element of a vector by a scalar and return a new vector containing the product.
template<class IT1, class IT2>
std::iterator_traits< IT1 > ::value_type	dot_prod (Subscript N, IT1 i1, IT2 i2)
	Dot product of values pointed to by two iterators.
template<class T>
T	dot_prod (const Vector< T > &A, const Vector< T > &B)
	Dot product of two vectors.

Detailed Description

Namespace containing Simple C++ Numerical Toolkit (SCPPNT) classes and functions.

Typedef Documentation

typedef SCPPNT_SUBSCRIPT_TYPE SCPPNT::Subscript

Subscript index type.

This definition describes the default SCPPNT data type used for indexing into SCPPNT matrices and vectors. The data type should signed and be wide enough to index into large arrays. It defaults to an "int", but can be overriden at compile time redefining SCPPNT_SUBSCRIPT_TYPE, e.g.

g++ -DSCPPNT_SUBSCRIPT_TYPE='int_fast64_t' ...

Definition at line 58 of file subscript.h.

Function Documentation

template<class M, class FUNC>

void SCPPNT::apply_columns	(	M &	matrix,
		FUNC &	f
	)			`[inline]`

Apply a function object to each column of a matrix.

Applies the function object f to each column of a matrix. f can modify the elements of a column.

Here is an example of using apply_columns:

include "scppnt/cmat.h" include "scppnt/rowcolfunc.h" include <iterator>

using namespace SCPPNT;

Function object that standardizes a sequence of elements to sum to 1. template <class it>=""> class Standardize { public: void operator ()(Subscript num_elements, IT iterator) { typename std::iterator_traits<IT>::value_type sum = 0; Subscript n = num_elements; IT it = iterator; for (; n--; ++it) sum += *it; for (it = iterator; num_elements--; ++it) *it /= sum; } };

Matrix<double> m(2, 2, "1.0 2.0 3.0 4.0");

Standardize each column of m to sum to 1 Standardize<Matrix<double>::column_iterator> f; apply_columns(m, f);

Definition at line 37 of file matalg.h.

00038   {
00039     typename M::columns_iterator icolumns = matrix.begin_columns();
00040 
00041     Subscript nrows = matrix.num_rows();
00042 
00043     for (Subscript i = matrix.num_columns(); i--; ++icolumns)
00044     {
00045       f(nrows, *icolumns);
00046     }
00047   }

template<class M, class FUNC>

void SCPPNT::apply_rows	(	M &	matrix,
		FUNC &	f
	)			`[inline]`

Apply a function object to each row of a matrix.

Applies the function object f to each row of a matrix. f can modify the elements of a row.

Here is an example of using apply_rows:

include "scppnt/cmat.h" include "scppnt/rowcolfunc.h" include <iterator>

using namespace SCPPNT;

Matrix<double> m(2, 2, "1.0 2.0 3.0 4.0");

Standardize each row of m to sum to 1 apply_rows(m, f); Standardize<Matrix<double>::row_iterator> f;

Definition at line 24 of file matalg.h.

00025   {
00026     typename M::rows_iterator irows = matrix.begin_rows();
00027 
00028     Subscript ncolumns = matrix.num_columns();
00029 
00030     for (Subscript i = matrix.num_rows(); i--; ++irows)
00031     {
00032       f(ncolumns, *irows);
00033     }
00034   }

template<class SPDMatrix, class SymmMatrix>

int SCPPNT::Cholesky_upper_factorization	(	SPDMatrix &	A,
		SymmMatrix &	L
	)			`[inline]`

Compute Cholesky factorization.

Modified version of Cholesky_upper_factorization from the Template Numerical Toolkit (TNT) using matrix iterators rather than matrix subscripting.

Only upper part of A is used. Cholesky factor is returned in lower part of L. Returns 0 if successful, 1 otherwise.

Definition at line 75 of file cholesky.h.

00077   {
00078     Subscript M = A.dim(1);
00079     Subscript N = A.dim(2);
00080 
00081     // assert(M == N);                 // make sure A is square
00082     if (M != N)
00083       throw BoundsError("Cholesky_upper_factorization");
00084 
00085     // readjust size of L, if necessary
00086 
00087     if (M != L.dim(1) || N != L.dim(2))
00088       L = SymmMatrix(N, N);
00089 
00090     Subscript i, j, k;
00091 
00092     typename SPDMatrix::element_type dot;
00093 
00094     typename SymmMatrix::diag_iterator diagL = L.begin_diagonal(1, 1);
00095     typename SPDMatrix::diag_iterator diagA = A.begin_diagonal(1, 1);
00096 
00097     for (j=1; j<=N; j++, ++diagL, ++diagA) // form column j of L
00098     {
00099       dot= 0;
00100 
00101       typename SymmMatrix::row_iterator irL = L.begin_row(j);
00102 
00103       for (i=1; i<j; i++) // for k= 1 TO j-1
00104       {
00105         dot += *irL * *irL; //dot = dot +  L(j,i)*L(j,i);
00106         ++irL;
00107       }
00108 
00109       *diagL = *diagA - dot; // L(j,j) = A(j,j) - dot;
00110 
00111       typename SymmMatrix::column_iterator icL = L.begin_column(j) + j;
00112       typename SPDMatrix::row_iterator irA = A.begin_row(j) + j;
00113 
00114       irL = L.begin_row(j);
00115 
00116       for (i=j+1; i<=N; i++, ++icL, ++irA)
00117       {
00118         dot = 0;
00119         typename SymmMatrix::row_iterator iL = L.begin_row(i);
00120         typename SymmMatrix::row_iterator jL =irL;
00121         for (k=1; k<j; k++, ++iL, ++jL)
00122           dot += *iL * *jL; // dot = dot +  L(i,k)*L(j,k);
00123         *icL = *irA - dot; // L(i,j) = A(j,i) - dot;
00124       }
00125 
00126       if (*diagL <= 0.0)
00127         return 1; // if (L(j,j) <= 0.0) return 1;
00128 
00129       *diagL = std::sqrt(*diagL); // L(j,j) = sqrt( L(j,j) );
00130 
00131       icL = L.begin_column(j) + j;
00132 
00133       for (i=j+1; i<=N; i++, ++icL)
00134         *icL /= *diagL; // L(i,j) = L(i,j) / L(j,j);
00135 
00136     }
00137     return 0;
00138   }

template<class IT1, class IT2>

std::iterator_traits<IT1>::value_type SCPPNT::dot_prod	(	Subscript	N,
		IT1	i1,
		IT2	i2
	)			`[inline]`

Dot product of values pointed to by two iterators.

Returns sum of products of corresponding elements pointed to by the iterators i1 and i2.

Parameters:

	N	Number of elements in dot product
	i1	Iterator to first sequence of elements to use for product.
	i2	Iterator to second sequence of elements to use for product.

Definition at line 773 of file vec.h.

Referenced by dot_prod().

00775   {
00776 
00777     typename std::iterator_traits<IT1>::value_type sum = 0;
00778     while(N--)
00779     {
00780       sum += *i1 * *i2;
00781       ++i1;
00782       ++i2;
00783     }
00784 
00785     return sum;

Here is the caller graph for this function:

template<class MaTRiX, class VecToRSubscript>

int SCPPNT::LU_factor	(	MaTRiX &	A,
		VecToRSubscript &	indx
	)			`[inline]`

Right-looking LU factorization algorithm (unblocked).

Factors matrix A into lower and upper triangular matrices (L and U respectively) in solving the linear equation Ax=b. LU_solve can be used to solve a system of linear equating using factorization computed.

Return value:

int (0 if successful, 1 otherwise)

Args:

Parameters:

	A	(input/output) Matrix(1:n, 1:n) In input, matrix to be factored. On output, overwritten with lower and upper triangular factors.
	indx	(output) Vector(1:n) Pivot vector. Describes how the rows of A were reordered to increase numerical stability.

Definition at line 106 of file lu.h.

00107   {
00108     // currently for 1-offset
00109     // vectors and matrices
00110     //if (A.lbound() != 1) throw InvalidArgument("Matrix not 1-offset", "SCPPNT::LU_factor");
00111     //if (indx.lbound() != 1) throw InvalidArgument("Subscript vector not 1-offset", "SCPPNT::LU_factor");
00112 
00113     typedef typename MaTRiX::column_iterator column_iterator;
00114     typedef typename MaTRiX::row_iterator row_iterator;
00115 
00116     Subscript M = A.num_rows();
00117     Subscript N = A.num_columns();
00118 
00119     if (M == 0 || N==0)
00120       return 0;
00121     if (indx.dim() != M)
00122       indx.newsize(M);
00123 
00124     Subscript i, j, k;
00125     Subscript jp;
00126 
00127     typename MaTRiX::element_type t;
00128 
00129     Subscript minMN = (M < N ? M : N); // min(M,N);
00130 
00131     typename MaTRiX::diag_iterator adiag = A.begin_diagonal(1, 1);
00132     typename MaTRiX::columns_iterator icolumns = A.begin_columns();
00133     for (j=1; j<= minMN; j++, ++adiag, ++icolumns)
00134     {
00135       // find pivot in column j and  test for singularity.
00136 
00137       jp = j;
00138       t = std::fabs(*adiag); //t = fabs(A(j,j));
00139 
00140       column_iterator ic = *icolumns + j;
00141       for (i=j+1; i<=M; i++, ++ic)
00142       {
00143         double ab = std::fabs(*ic);
00144         if (ab > t) // if ( fabs(A(i,j)) > t)
00145         {
00146           jp = i;
00147           t = ab; //t = fabs(A(i,j));
00148         }
00149       }
00150 
00151       indx(j) = jp;
00152 
00153       // jp now has the index of maximum element 
00154       // of column j, below the diagonal
00155 
00156       if ((*icolumns)[jp-1] == 0) // if ( A(jp,j) == 0 )
00157         return 1; // factorization failed because of zero pivot
00158 
00159       if (jp != j) // swap rows j and jp
00160       {
00161         row_iterator irowj = A.begin_row(j);
00162         row_iterator irowjp = A.begin_row(jp);
00163 
00164         for (k=1; k<=N; k++, ++irowj, ++irowjp)
00165         {
00166           t = *irowj; //  t = A(j,k);
00167           *irowj = *irowjp; //  A(j,k) = A(jp,k);
00168           *irowjp = t; // A(jp,k) =t;
00169         }
00170       }
00171 
00172       if (j<M) // compute elements j+1:M of jth column
00173       {
00174         // note A(j,j), was A(jp,p) previously which was
00175         // guarranteed not to be zero (Label #1)
00176         //
00177         typename MaTRiX::element_type recp = 1.0 / *adiag; // 1.0 / A(j,j)
00178 
00179         ic = *icolumns + j;
00180         for (k=j+1; k<=M; k++, ++ic)
00181           *ic *= recp; //A(k,j) *= recp;
00182 
00183       }
00184 
00185       if (j < minMN)
00186       {
00187         // rank-1 update to trailing submatrix:   E = E - x*y;
00188         //
00189         // E is the region A(j+1:M, j+1:N)
00190         // x is the column vector A(j+1:M,j)
00191         // y is row vector A(j,j+1:N)
00192 
00193         Subscript ii, jj;
00194 
00195         column_iterator ix = *icolumns + j;
00196         for (ii=j+1; ii<=M; ii++, ++ix) // for (ii=j+1; ii<=M; ii++)
00197         {
00198           row_iterator iE = A.begin_row(ii) + j;
00199           row_iterator iy = A.begin_row(j) + j;
00200           for (jj=j+1; jj<=N; jj++, ++iE, ++iy)
00201             // for (jj=j+1; jj<=N; jj++)
00202             *iE -= *ix * *iy; // A(ii,jj) -= A(ii,j)*A(j,jj);
00203         }
00204       }
00205     }
00206     return 0;
00207   }

template<class MaTRiX, class VecToR, class VecToRSubscripts>

int SCPPNT::LU_solve	(	const MaTRiX &	A,
		const VecToRSubscripts &	indx,
		VecToR &	b
	)			`[inline]`

Solve system of linear equations using LU factorization computed from LU_factor.

Typical usage for solving system of linear equations Ax = b:

Matrix<double> A; Vector<Subscript> ipiv; Vector<double> b;

LU_Factor(A,ipiv); LU_Solve(A,ipiv,b);

Now b has the solution x. Note that both A and b are overwritten.

Definition at line 224 of file lu.h.

00226   {
00227     // currently for 1-offset
00228     // vectors and matrices
00229     //if (A.lbound() != 1) throw InvalidArgument("Matrix not 1-offset", "SCPPNT::LU_solve");
00230     //if (indx.lbound() != 1) throw InvalidArgument("Subscript vector not 1-offset", "SCPPNT::LU_solve");
00231     //if (b.lbound() != 1) throw InvalidArgument("Vector not 1-offset", "SCPPNT::LU_solve");
00232 
00233     typedef typename MaTRiX::const_row_iterator row_iterator;
00234     typedef typename VecToR::iterator vector_iterator;
00235 
00236     Subscript i, ii=0, ip, j;
00237     Subscript n = b.dim();
00238     typename MaTRiX::element_type sum;
00239 
00240     for (i=1; i<=n; i++)
00241     {
00242       ip=indx(i);
00243       sum=b(ip);
00244       b(ip)=b(i);
00245       if (ii)
00246       {
00247         row_iterator Aj = A.begin_row(i) + ii-1;
00248         vector_iterator bj = b.begin() + ii - 1;
00249         for (j=ii; j<=i-1; j++, ++Aj, ++bj)
00250           sum -= *Aj * *bj; //sum -= A(i,j)*b(j);
00251       }
00252       else if (sum)
00253         ii=i;
00254       b(i)=sum;
00255     }
00256 
00257     typename MaTRiX::const_diag_iterator adiag = A.begin_diagonal(1, 1) + n - 1;
00258     for (i=n; i>=1; --i, --adiag)
00259     {
00260       sum=b(i);
00261       row_iterator Aj = A.begin_row(i) + i;
00262       vector_iterator bj = b.begin() + i;
00263       for (j=i+1; j<=n; j++, ++Aj, ++bj)
00264         sum -= *Aj * *bj; //sum -= A(i,j)*b(j);
00265       b(i)=sum / *adiag; //  b(i)=sum/A(i,i)
00266     }
00267     return 0;
00268   }

template<class MAT, class IT1, class IT2>

void SCPPNT::MatTimesVec	(	const MAT &	A,
		IT1	begin,
		IT2	result
	)			`[inline]`

Multiply a matrix times a vector using iterators to access the vector elements.

Parameters:

	A	Matrix to multiply times vector.
	begin	Iterator over elements of vector to be multiplied by matrix
	result	Iterator over vector to hold result of multiplication.

Definition at line 325 of file matop.h.

Referenced by matrix_times_vector().

00326   {
00327 
00328     typedef typename MAT::value_type value_type;
00329 
00330     Subscript N = A.num_columns();
00331     Subscript M = A.num_rows();
00332 
00333     typename MAT::const_rows_iterator irows = A.begin_rows();
00334     for (Subscript i=M; i--; ++irows, ++result)
00335     {
00336       value_type sum = 0;
00337       typename MAT::const_row_iterator rowi = *irows;
00338       IT1 iv = begin;
00339       for (Subscript j=N; j--; ++iv, ++rowi)
00340       {
00341         sum += *rowi * *iv;
00342       }
00343 
00344       *result = sum;
00345     }
00346   }

Here is the caller graph for this function:

template<class T>

std::ostream& SCPPNT::operator<<	(	std::ostream &	s,
		const Vector< T > &	A
	)			`[inline]`

Write a vector to an output stream.

Writes the number of elements in the vector, followed by a new line. The elements of the vector are then written on consecutive lines.

Definition at line 673 of file vec.h.

00675   {
00676     Subscript N=A.dim();
00677 
00678     s << N << std::endl;
00679 
00680     for (Subscript i=0; i<N; i++)
00681     s << A[i] << " " << std::endl;
00682     s << std::endl;
00683 
00684     return s;

template<class Array2D>

std::ostream& SCPPNT::operator<<	(	std::ostream &	s,
		const Transpose_View< Array2D > &	A
	)			`[inline]`

Write a matrix to an output stream.

Writes the number of rows and columns, followed by a new line. The following lines contain the elements in the rows of the matrix.

Definition at line 275 of file transv.h.

References SCPPNT::Transpose_View< Array2D >::lbound(), SCPPNT::Transpose_View< Array2D >::num_columns(), and SCPPNT::Transpose_View< Array2D >::num_rows().

00277   {
00278     Subscript M=A.num_rows();
00279     Subscript N=A.num_columns();
00280 
00281     Subscript start = A.lbound();
00282     Subscript Mend = M + A.lbound() - 1;
00283     Subscript Nend = N + A.lbound() - 1;
00284 
00285     s << M << "  " << N << std::endl;
00286     for (Subscript i=start; i<=Mend; i++)
00287     {
00288       for (Subscript j=start; j<=Nend; j++)
00289       {
00290         s << A(i, j) << " ";
00291       }
00292       s << std::endl;
00293     }
00294 
00295     return s;
00296   }

Here is the call graph for this function:

template<class Array2D>

std::ostream& SCPPNT::operator<<	(	std::ostream &	s,
		const const_Region2D< Array2D > &	A
	)			`[inline]`

Write a matrix to an output stream.

Writes the number of rows and columns, followed by a new line. The following lines contain the elements in the rows of the matrix.

Definition at line 1481 of file region2d.h.

01482   {
01483     Subscript M=A.num_rows();
01484     Subscript N=A.num_columns();
01485 
01486     Subscript start = A.lbound();
01487     Subscript Mend = M + A.lbound() - 1;
01488     Subscript Nend = N + A.lbound() - 1;
01489 
01490     s << M << "  " << N << std::endl;
01491     for (Subscript i=start; i<=Mend; i++)
01492     {
01493       for (Subscript j=start; j<=Nend; j++)
01494       {
01495         s << A(i,j) << " ";
01496       }
01497       s << std::endl;
01498     }
01499 
01500     return s;
01501   }

template<class Array2D>

std::ostream& SCPPNT::operator<<	(	std::ostream &	s,
		const Region2D< Array2D > &	A
	)			`[inline]`

Write a matrix to an output stream.

Writes the number of rows and columnumns, followed by a new line. The following lines contain the elements in the rows of the matrix.

Definition at line 958 of file region2d.h.

00959   {
00960     Subscript M=A.num_rows();
00961     Subscript N=A.num_columns();
00962 
00963     Subscript start = A.lbound();
00964     Subscript Mend = M + A.lbound() - 1;
00965     Subscript Nend = N + A.lbound() - 1;
00966 
00967     s << M << "  " << N << std::endl;
00968     for (Subscript i=start; i<=Mend; i++)
00969     {
00970       for (Subscript j=start; j<=Nend; j++)
00971       {
00972         s << A(i,j) << " ";
00973       }
00974       s << std::endl;
00975     }
00976     return s;
00977   }

template<class T>

std::ostream& SCPPNT::operator<<	(	std::ostream &	s,
		const Matrix< T > &	A
	)			`[inline]`

Write a matrix to an output stream.

Writes the number of rows and columns, followed by a new line. The following lines contain the elements in the rows of the matrix.

Definition at line 1204 of file cmat.h.

01205   {
01206     Subscript M = A.num_rows();
01207     Subscript N = A.num_columns();
01208 
01209     s << M << " " << N << "\n";
01210 
01211     for (Subscript i=0; i<M; i++)
01212     {
01213       for (Subscript j=0; j<N; j++)
01214       {
01215         s << A[i][j] << " ";
01216       }
01217       s << "\n";
01218     }
01219     return s;
01220   }

template<class T>

std::istream& SCPPNT::operator>>	(	std::istream &	s,
		Vector< T > &	A
	)			`[inline]`

Read elements of a vector from an input stream.

Reads the number of elements in the vector followed by the elements of the vector. All values read should be separated by white space.

Definition at line 694 of file vec.h.

References SCPPNT::Vector< T >::newsize().

00696   {
00697 
00698     Subscript N;
00699 
00700     s >> N;
00701 
00702     A.newsize(N);
00703 
00704     for (Subscript i=0; i<N; i++)
00705     s >> A[i];
00706 
00707     return s;

Here is the call graph for this function:

template<class T>

std::istream& SCPPNT::operator>>	(	std::istream &	s,
		Matrix< T > &	A
	)			`[inline]`

Read elements of a matrix to an input stream.

Reads the number of rows and columns, followed by the elements in each row of the matrix.

Definition at line 1229 of file cmat.h.

References SCPPNT::Matrix< T >::newsize().

01230   {
01231     Subscript M, N;
01232 
01233     s >> M >> N;
01234 
01235     A.newsize(M,N);
01236 
01237     for (Subscript i=0; i<M; i++)
01238     for (Subscript j=0; j<N; j++)
01239     {
01240       s >> A[i][j];
01241     }
01242     return s;
01243   }

Here is the call graph for this function:

template<class M, class V, class FUNC>

V SCPPNT::over_columns	(	M &	matrix,
		FUNC &	f
	)			`[inline]`

Apply a function object to each row of a matrix to produce a scalar, and return a vector containing these scalars.

Apply a function to each row of a matrix that returns a single value (the function is applied over columns of each row). These values are put into a vector that is returned (the number of elements in the vector is the number of rows in the matrix).

Here is an example of using over_columns:

include "scppnt/cmat.h" include "scppnt/vec.h" include "scppnt/rowcolfunc.h" include <iterator>

using namespace SCPPNT;

Function object that returns the sum of each element of an array. template <class it>=""> class Sum { public: typename std::iterator_traits<IT>::value_type operator ()(Subscript num_elements, IT iterator) { typename std::iterator_traits<IT>::value_type sum = 0; for(; num_element--; ++iterator) sum += *iterator; } };

Matrix<double> m(2, 2, "1.0 2.0 3.0 4.0");

Get sum of each column of m Sum<Matrix<double>::row_iterator> f; Vector<double> sums = over_columns<Matrix<double>,Vector<double>,Sum>(m, f);

Definition at line 55 of file matalg.h.

00056   {
00057     typename M::rows_iterator irows = matrix.begin_rows();
00058     V vec(matrix.num_rows());
00059 
00060     Subscript ncolumns = matrix.num_columns();
00061 
00062     typename V::iterator iv = vec.begin();
00063     for (Subscript i = matrix.num_rows(); i--; ++irows, ++iv)
00064     {
00065       *iv = f(ncolumns, *irows);
00066     }
00067 
00068     return vec;
00069   }

template<class M, class V, class FUNC>

V SCPPNT::over_rows	(	M &	matrix,
		FUNC &	f
	)			`[inline]`

Apply a function object to each column of a matrix to produce a scalar, and return a vector containing these scalars.

Apply a function to each column of a matrix that returns a single value (the function is applied over rows of each column). These values are put into a vector that is returned (the number of elements in the vector is the number of columns in the matrix).

Here is an example of using over_rows:

include "scppnt/cmat.h" include "scppnt/vec.h" include "scppnt/rowcolfunc.h" include <iterator>

using namespace SCPPNT;

Matrix<double> m(2, 2, "1.0 2.0 3.0 4.0");

Get sum of each row of m Sum<Matrix<double>::column_iterator> f; Vector<double> sums = over_columns<Matrix<double>,Vector<double>,Sum>(m, f);

Definition at line 77 of file matalg.h.

00078   {
00079     typename M::columns_iterator icolumns = matrix.begin_columns();
00080     V vec(matrix.num_columns());
00081 
00082     Subscript nrows = matrix.num_rows();
00083 
00084     typename V::iterator iv = vec.begin();
00085     for (Subscript i = matrix.num_columns(); i--; ++icolumns, ++iv)
00086     {
00087       *iv = f(nrows, *icolumns);
00088     }
00089 
00090     return vec;
00091   }

template<class MaTRiX, class Vector>

int SCPPNT::QR_factor	(	MaTRiX &	A,
		Vector &	C,
		Vector &	D
	)			`[inline]`

Classical QR factorization, based on Stewart[1973].

This algorithm computes the factorization of a matrix A into a product of an orthognal matrix (Q) and an upper triangular matrix (R), such that QR = A.

Returns:

0 if successful, 1 if A is detected to be singular

Parameters:

Parameters:

	A	(in/out): On input, A is square, Matrix(1:N, 1:N), that represents the matrix to be factored. On output, Q and R is encoded in the same Matrix(1:N,1:N) in the following manner: R is contained in the upper triangular section of A, except that R's main diagonal is in D. The lower triangular section of A represents Q, where each column j is the vector Qj = I - uj*uj'/pi_j.
	C	(output): vector of Pi[j]
	D	(output): main diagonal of R, i.e. D(i) is R(i,i)

Definition at line 93 of file qr.h.

References SCPPNT::Vector< T >::begin(), SCPPNT::Vector< T >::lbound(), SCPPNT::Vector< T >::newsize(), sign(), and SCPPNT::Vector< T >::size().

00094   {
00095     //assert(A.lbound() == 1);        // ensure these are all
00096     if (A.lbound() != 1)
00097       throw BoundsError("SCPPNT::QR_factor");
00098 
00099     //assert(C.lbound() == 1);        // 1-based arrays and vectors
00100     if (C.lbound() != 1)
00101       throw BoundsError("SCPPNT::QR_factor");
00102 
00103     //assert(D.lbound() == 1);
00104     if (D.lbound() != 1)
00105       throw BoundsError("SCPPNT::R_factor");
00106 
00107     Subscript M = A.num_rows();
00108     Subscript N = A.num_columns();
00109 
00110     // assert(M == N);                 // make sure A is square
00111     if (M != N)
00112       throw BoundsError("QR_factor");
00113 
00114     Subscript i, j, k;
00115     typename MaTRiX::element_type eta, sigma, sum;
00116 
00117     // adjust the shape of C and D, if needed...
00118 
00119     if (N != C.size())
00120       C.newsize(N);
00121     if (N != D.size())
00122       D.newsize(N);
00123 
00124     typename MaTRiX::diag_iterator diagA = A.begin_diagonal(1, 1);
00125     typename Vector::iterator iC = C.begin();
00126     typename Vector::iterator iD = D.begin();
00127     typename MaTRiX::columns_iterator columnsA = A.begin_columns();
00128     typename Vector::iterator dn = D.begin() + N - 1;
00129     typename MaTRiX::row_iterator ann = A.begin_row(N) + N - 1;
00130 
00131     for (k=1; k<N; k++, ++columnsA, ++iC, ++iD)
00132     {
00133       // eta = max |M(i,k)|,  for k <= i <= n
00134       //
00135       eta = 0;
00136       typename MaTRiX::column_iterator columnA = *columnsA + k - 1;
00137 
00138       for (i=k; i<=N; i++, ++columnA)
00139       {
00140         double absA = std::fabs(*columnA); // double absA = std::fabs(A(i,k));
00141         eta = (absA > eta ? absA : eta );
00142       }
00143 
00144       if (eta == 0) // matrix is singular
00145       {
00146         // cerr << "QR: k=" << k << "\n";
00147         return 1;
00148       }
00149 
00150       // form Qk and premiltiply M by it
00151       //
00152       columnA = *columnsA + k - 1;
00153       for (i=k; i<=N; i++, ++columnA)
00154         *columnA /= eta; //A(i,k)  = A(i,k) / eta;
00155 
00156       sum = 0;
00157       columnA = *columnsA + k - 1;
00158       for (i=k; i<=N; i++, ++columnA)
00159         sum += *columnA * *columnA; //sum = sum + A(i,k)*A(i,k);
00160 
00161       sigma = sign(*diagA) * std::sqrt(sum); // sigma = sign(A(k,k)) *  sqrt(sum);
00162 
00163       *diagA += sigma; // A(k,k) = A(k,k) + sigma;
00164       *iC = sigma * *diagA; // C(k) = sigma * A(k,k);
00165       ++diagA;
00166       *iD = -eta * sigma; //D(k) = -eta * sigma;
00167 
00168       typename MaTRiX::columns_iterator columnsj = columnsA+1;
00169       for (j=k+1; j<=N; j++, ++columnsj)
00170       {
00171         sum = 0;
00172         columnA = *columnsA + k - 1;
00173         typename MaTRiX::column_iterator columnj = *columnsj + k - 1;
00174 
00175         for (i=k; i<=N; i++, ++columnj, ++columnA)
00176           sum += *columnA * *columnj; // sum = sum + A(i,k)*A(i,j);
00177 
00178         sum /= *iC; // sum = sum / C(k);
00179 
00180         columnA = *columnsA + k - 1;
00181         columnj = *columnsj + k - 1;
00182         for (i=k; i<=N; i++, ++columnA, ++columnj)
00183           *columnj -= sum * *columnA; // A(i,j) = A(i,j) - sum * A(i,k);
00184       }
00185 
00186       *dn = *ann; // D(N) = A(N,N);
00187     }
00188     return 0;
00189   }

Here is the call graph for this function:

SCPPNT Namespace Reference

Classes

Typedefs

Functions

Detailed Description

Typedef Documentation

Function Documentation