4.5. Linear System Solvers

These functions and classes solve linear systems and also perform singular value decomposition.

4.5.1. Cholesky Decomposition Solver

This section describes the Cholesky decomposition processing object provided by VSIPL++.

4.5.1.1. Class template chold<>

The class chold uses Cholesky decomposition to solve linear systems.

template <typename T,
          return_mechanism_type ReturnMechanism = by_value>
class chold;

Template parameters

T

The value type used for the decomposition object. May be single- or double-precision floating-point, and either real or complex.

ReturnMechanism

The return mechanism type indicates whether the Cholesky decomposition object's solve() function returns results by value (by_value) or by reference (by_reference) into matrices provided by the caller. See section 2.7.7 for details.

4.5.1.2. Constructor

chold(mat_uplo uplo, length_type len);

Description: Constructs a chold object that will decompose len by len positive definite matrices. The parameter uplo controls whether an upper LU or lower LU decomposition will be performed (see section 2.7.6). Note also that chold objects may also be copied (constructed) from other chold objects.

Requirements:  The parameter len must be greater than zero.

4.5.1.3. Accessor functions

mat_uplo uplo() const;
length_type length() const;

Description: Report the attributes of this chold object. The length() function returns the length, equal to the number of rows as well as the number of columns in the decomposed matrix. The uplo() function indicates whether the lower half or upper half of a decomposed matrix is referenced.

4.5.1.4. Solve Systems

template <typename Block>
bool 
decompose(Matrix<T, Block> A);

Description: Performs Cholesky decomposition of A. When it can be used, Cholesky decomposition is twice as efficient as normal LU decomposition.

Requirements: The matrix A must be the same size as specified in the constructor. Note that the contents may be overwritten, therefore A should not be modified until all solve() calls have been performed. The matrix A must be both symmetric and positive definite for Cholesky decomposition to work.

Result: False is returned if the decomposition fails.

4.5.1.5. Solve Systems (by_value)

This function is available only if the chold class template is parameterized with ReturnMechanism=by_value.

template <mat_op_type tr,
          typename    Block>
const_Matrix<T, unspecified>
solve(const_Matrix<T, Block> B);

Description: Solves a linear system.

The parameter tr controls what type of operation is performed on A when solving the system (see section 2.7.3).

If tr == mat_trans and T is not a specialization of complex, then ATX = B is solved.

If tr == mat_herm and T is a specialization of complex, then AHX = B is solved.

Otherwise, AX = B is solved.

Requirements: The number of rows of B must be equal to the value returned by length(). A successful call to decompose() must have occurred.

Result: Returns the solution to the linear system. The returned matrix’s block type may be a different type from Block.

4.5.1.6. Solve Systems (by_reference)

This function is available only if the chold class template is parameterized with ReturnMechanism=by_reference.

template <mat_op_type tr,
          typename    Block0,
          typename    Block1>
bool
solve(
  const_Matrix<T, Block0> B,
  Matrix<T, Block1>       X);

Description: Solves a linear system.

The parameter tr controls what type of operation is performed on A when solving the system (see section 2.7.3).

If tr == mat_trans and T is not a specialization of complex, then ATX = B is solved.

If tr == mat_herm and T is a specialization of complex, then AHX = B is solved.

Otherwise, AX = B is solved.

Requirements: The number of rows of both B and X must be equal to the value returned by length(). A successful call to decompose() must have occurred.

Result: Stores the solution in X. True is returned if the computation succeeds.

4.5.2. Covariance Solver

This section describes the covariance linear system solver function provided by VSIPL++.

4.5.2.1. Solve Systems (return by value)

template <typename T,
          typename Block0,
          typename Block1>
Matrix<T, unspecified>
covsol(
  Matrix<T, Block0>       A,
  const_Matrix<T, Block1> B);

Description: Solves the covariance linear system ATAX = B for X if type T is real and AHAX = B for X if type T is complex.

Requirements: The matrix A is of size M by N (where M >= N) and is of rank N. The matrix B is of size N by K. The type T may be single- or double-precision floating-point, and either real or complex. Temporary workspace may be allocated, which may result in nondeterministic execution time. As an alternative, use the QR routines.

Result: The solution X is returned and is of size N by K. Note that A may be overwritten. The returned matrix’s block type may be a different type from Block0 or Block1.

4.5.2.2. Solve Systems (return by reference)

template <typename T,
          typename Block0,
          typename Block1,
          typename Block2>
Matrix<T, unspecified>
covsol(
  Matrix<T, Block0>       A,
  const_Matrix<T, Block1> B,
  Matrix<T, Block2>       X);

Description: Solves the covariance linear system ATAX = B for X if type T is real and AHAX = B for X if type T is complex.

Requirements: The matrix A is of size M by N (where M >= N) and is of rank N. The matrix B is of size N by K. The matrix X is also of size N by K. The type T may be single- or double-precision floating-point, and either real or complex. Temporary workspace may be allocated, which may result in nondeterministic execution time. As an alternative, use the QR routines.

Result: The solution is placed in X. Note that A may be overwritten.

4.5.3. Linear Least Squares Solver

This section describes the linear least squares solver function provided by VSIPL++.

4.5.3.1. Solve Systems (return by value)

template <typename T,
          typename Block0,
          typename Block1>
Matrix<T, unspecified>
llsqsol(
  Matrix<T, Block0>       A,
  const_Matrix<T, Block1> B);

Description: Solves the linear least squares problem minX ||AX − B||2 for X.

Requirements: The matrix A is of size M by N (where M >= N) and is of rank N. The matrix B is of size M by K. The type T may be single- or double-precision floating-point, and either real or complex. Temporary workspace may be allocated, which may result in nondeterministic execution time. As an alternative, use the QR routines.

Result: Returns the solution X which is of size N by K. Note that A may be overwritten. The returned matrix’s block type may be a different type from Block0 or Block1.

4.5.3.2. Solve Systems (return by reference)

template <typename T,
          typename Block0,
          typename Block1,
          typename Block2>
Matrix<T, Block2>
llsqsol(
  Matrix<T, Block0>       A,
  const_Matrix<T, Block1> B,
  Matrix<T, Block2>       X);

Description: Solves the linear least squares problem minX ||AX − B||2 for X.

Requirements: The matrix A is of size M by N (where M >= N) and is of rank N. The matrix B is of size M by K. The matrix X is of size N by K. The type T may be single- or double-precision floating-point, and either real or complex. Temporary workspace may be allocated, which may result in nondeterministic execution time. As an alternative, use the QR routines.

Result: Stores the solution in X and returns it. Note that A may be overwritten.

4.5.4. LU Decomposition Solver

This section describes the LU (lower and upper triangular) decomposition processing object provided by VSIPL++.

4.5.4.1. Class template lud<>

The class lud performs LU decomposition to solve linear systems.

template <typename T,
          return_mechanism_type ReturnMechanism = by_value>
class lud;

Template parameters

T

The value type used for the decomposition object. May be real or complex, single- or double-precision floating-point types.

ReturnMechanism

The return mechanism type indicates whether the LU decomposition object's solve() function returns results by value (by_value) or by reference (by_reference) into matrices provided by the caller. See section 2.7.7 for details.

4.5.4.2. Constructor

lud(length_type len);

Description: Constructs an lud object that will decompose len by len matrices. Note also that lud objects may also be copied (constructed) from other lud objects.

Requirements:  The parameter len must be greater than zero.

4.5.4.3. Accessor functions

length_type length() const;

Description: Report the length attribute of this lud object, equal to the number of rows as well as the number of columns in the decomposed matrix.

4.5.4.4. Solve Systems

template <typename Block>
bool 
decompose(Matrix<T, Block> A);

Description: Performs LU decomposition of A.

Requirements: The matrix A must be the same size as specified in the constructor. Note that the contents may be overwritten, therefore A should not be modified until all solve() calls have been performed.

Result: False is returned if the decomposition fails.

4.5.4.5. Solve Systems (by_value)

This function is available only if the lud class template is parameterized with ReturnMechanism=by_value.

template <mat_op_type tr,
          typename    Block>
const_Matrix<T, unspecified>
solve(const_Matrix<T, Block> B);

Description: Solves a linear system. The parameter tr controls what type of operation is performed on A when solving the system (see section 2.7.3).

If tr == mat_trans and T is not a specialization of complex, then ATX = B is solved.

If tr == mat_herm and T is a specialization of complex, then AHX = B is solved.

Otherwise, AX = B is solved.

Requirements: The number of rows of B must be equal to the value returned by length(). A successful call to decompose() must have occurred.

Result: Returns the solution to the linear system.

4.5.4.6. Solve Systems (by_reference)

This function is available only if the lud class template is parameterized with ReturnMechanism=by_reference.

template <mat_op_type tr,
          typename    Block0,
          typename    Block1>
bool
solve(
  const_Matrix<T, Block0> B,
  Matrix<T, Block1>       X);

Description: Solves a linear system. The parameter tr controls what type of operation is performed on A when solving the system (see section 2.7.3).

If tr == mat_trans and T is not a specialization of complex, then ATX = B is solved.

If tr == mat_herm and T is a specialization of complex, then AHX = B is solved.

Otherwise, AX = B is solved.

Requirements: The number of rows of both B and X must be equal to the value returned by length(). A successful call to decompose() must have occurred.

Result: Stores the solution in X. True is returned if the computation succeeds.

Exceptions: If the backends enabled do not support the requested LU decomposition, a vsip::impl::unimplemented exception will be thrown. This is a deviation from the VSIPL++ spec.

4.5.5. QR Decomposition Solver

This section describes the QR decomposition processing object provided by VSIPL++.

4.5.5.1. Class template qrd<>

The class qrd performs QR decomposition and solves linear systems.

template <typename T,
          return_mechanism_type ReturnMechanism = by_value>
class qrd;

Template parameters

T

The value type used for the decomposition object. May be real or complex, single- or double-precision floating-point types.

ReturnMechanism

The return mechanism type indicates whether the QR solver and product functions return results by value (by_value) or by reference (by_reference) into matrices provided by the caller. See section 2.7.7 for details.

4.5.5.2. Constructor

qrd(length_type rows, length_type columns, storage_type st);

Description: Constructs a qrd object. The parameters rows and columns refer to the size of the Q matrix. The parameter st controls how much of the Q matrix is stored (see section 2.7.8). Note also that qrd objects may also be copied (constructed) from other qrd objects.

Requirements:  The number of rows must be greater than or equal to the number of columns.

4.5.5.3. Accessor functions

length_type rows() const;
length_type columns() const;
storage_type qstorage() const;

Description: Report the various attributes of this qrd object. The number of rows is returned by rows() and the number of columns by columns(). The function st() returns how the decomposition matrix Q is stored by this object, if at all.

4.5.5.4. Solve Systems

template <typename Block>
bool 
decompose(Matrix<T, Block> A);

Description: Performs a QR decomposition of the matrix A into matrices Q and R.

Requirements: The matrix A must be the same size as specified in the constructor. Note that the contents may be overwritten therefore A should not be modified prior to calling any other function.

Result: False is returned if the decomposition fails because A does not have full column rank. Note: If T is a specialization of complex, Q is unitary. Otherwise, Q is orthogonal. R is an upper triangular matrix. If A has full rank, then R is a nonsingular matrix. No column interchanges are performed.

4.5.5.5. Solve Systems (by_value)

This function is available only if the qrd class template is parameterized with ReturnMechanism=by_value

template <mat_op_type       tr,
          product_side_type ps,
          typename          Block>
const_Matrix<T, unspecified>
prodq(const_Matrix<T, Block> C);

Description: Calculates the product of Q and C. The parameter tr controls what type of operation is performed on C before the product is computed (see section 2.7.3) and ps determines what side of the product Q is placed on (see section 2.7.5). The actual product and its number of rows and columns (shown in the table below) depends on the values of tr, ps, and qstorage() and whether T is not or is a specialization of complex.

For qstorage() == qrd_saveq1,

 ps == mat_lsideps == mat_rside
tr == mat_ntransQC, rows(), sCQ, s, columns()
tr == mat_trans, TQTC, columns(), sCQT, s, rows()
tr == mat_herm, complex<T>QHC, columns(), sCQH, s, rows()

where s is an arbitrary positive value.

For qstorage() == qrd_saveq,

 ps == mat_lsideps == mat_rside
tr == mat_ntransQC, rows(), sCQ, s, rows()
tr == mat_trans, TQTC, rows(), sCQT, s, rows()
tr == mat_herm, complex<T>QHC, rows(), sCQH, s, rows()

Requirements: A successful call to decompose() must have occurred for this object with qstorage() equaling either qrd_saveq1 or qrd_saveq. Otherwise, the behavior is undefined. The number of rows and columns (shown in the table below) of C depend on the values of tr, ps, and qstorage().

For qstorage() == qrd_saveq1,

 ps == mat_lsideps == mat_rside
tr == mat_ntranscolumns(), ss, rows()
tr == mat_transrows(), ss, columns()
tr == mat_hermrows(), ss, columns()

where s is the same variable as above.

For qstorage() == qrd_saveq,

 ps == mat_lsideps == mat_rside
tr == mat_ntransrows(), ss, rows()
tr == mat_transrows(), ss, rows()
tr == mat_hermrows(), ss, rows()

Result: Returns the product of Q and C. The returned matrix’s block type may be a different type from Block.

template <mat_op_type tr,
          typename    Block>
const_Matrix<T, unspecified>
rsol(
  const_Matrix<T, Block> B,
  T const                 alpha);

Description: Solves a linear system for X. The parameter tr controls what type of operation is performed on R before the system is solved (see section 2.7.3).

If tr == mat_trans and T is not a specialization of complex, then RTX = alpha B is solved.

If tr == mat_herm and T is a specialization of complex, then RHX = alpha B is solved.

Otherwise, RX = alpha B is solved.

Requirements: The number of rows in B must be equal to the value returned by columns(). A successful call to decompose() must have occurred.

Result: Returns a constant matrix X containing the solution.

template <typename Block>
const_Matrix<T, unspecified>
covsol(const_Matrix<T, Block> B);

Description: Solves a covariance linear system for X.

If T is not a specialization of complex, then ATAX = B is solved, where A is the matrix given to the most recent call to decompose().

If T is a specialization of complex, then AHAX = B is solved.

Requirements: The number of rows in B must be equal to the value returned by columns(). Note also that X and B are element-conformant

Result: Returns a matrix X containing the solution.

template <typename Block>
Matrix<T, unspecified>
lsqsol(const_Matrix<T, Block> B)

Description: Solves the linear least squares problem minX ||AX − B||2 for X, where A is the matrix given to the most recent call to decompose().

Requirements: The number of rows in B must be equal to the value returned by rows(). The number of rows in X will equal the value returned by columns().

Result: Returns a constant matrix X containing the solution.

4.5.5.6. Solve Systems (by_reference)

This function is available only if the qrd class template is parameterized with ReturnMechanism=by_reference

template <mat_op_type       tr,
          product_side_type ps,
          typename          Block0,
          typename          Block1>
bool
prodq(
  const_Matrix<T, Block0> C,
  Matrix<T, Block1>       X);

Description: Calculates the product of Q and C. The parameter tr controls what type of operation is performed on C before the product is computed (see section 2.7.3) and ps determines what side of the product Q is placed on (see section 2.7.5). The actual product depends on the values of tr, and whether T is not or is a specialization of complex:

 ps == mat_lsideps == mat_rside
tr == mat_ntransQCCQ
tr == mat_trans, TQTCCQT
tr == mat_herm, complex<T>QHCCQH

Requirements: A successful call to decompose() must have occurred for this object with qstorage() equaling either qrd_saveq1 or qrd_saveq. Otherwise, the behavior is undefined. The number of rows and columns (shown in the table below) of C depend on the values of tr, ps, and qstorage().

For qstorage() == qrd_saveq1,

 ps == mat_lsideps == mat_rside
tr == mat_ntranscolumns(), ss, rows()
tr == mat_transrows(), ss, columns()
tr == mat_hermrows(), ss, columns()

where s is an arbitrary positive value.

For qstorage() == qrd_saveq,

 ps == mat_lsideps == mat_rside
tr == mat_ntransrows(), ss, rows()
tr == mat_transrows(), ss, rows()
tr == mat_hermrows(), ss, rows()

The number of rows and columns of X (shown in the table below) depend on the values of tr, ps, and qstorage().

For qstorage() == qrd_saveq1,

 ps == mat_lsideps == mat_rside
tr == mat_ntransrows(), ss, columns()
tr == mat_transcolumns(), ss, rows()
tr == mat_hermcolumns(), ss, rows()

where s is the same variable as above.

For qstorage() == qrd_saveq,

 ps == mat_lsideps == mat_rside
tr == mat_ntransrows(), ss, rows()
tr == mat_transrows(), ss, rows()
tr == mat_hermrows(), ss, rows()

Result: Calculates the product of Q and C stores it in X.

template <mat_op_type tr,
          typename    Block0,
          typename    Block1>
bool
rsol(
  const_Matrix<T, Block0> B,
  T const                 alpha,
  Matrix<T, Block1>       X);

Description: Solves a linear system for X. The parameter tr controls what type of operation is performed on R before the system is solved (see section 2.7.3).

If tr == mat_trans and T is not a specialization of complex, then RTX = alpha B is solved.

If tr == mat_herm and T is a specialization of complex, then RHX = alpha B is solved.

Otherwise, RX = alpha * B is solved.

Requirements: The number of rows in B must be equal to the value returned by columns(). A successful call to decompose() must have occurred.

Result: Stores the solution in X.

template <typename Block0,
          typename Block1>
bool
covsol(
  const_Matrix<T, Block0> B,
  Matrix<T, Block1>       X);

Description: Solves a covariance linear system for X. If T is not a specialization of complex, then ATAX = B is solved, where A is the matrix given to the most recent call to decompose(). If T is a specialization of complex, then AHAX = B is solved.

Requirements: The number of rows in B must be equal to the value returned by columns(). Note also that X is modifiable and element-conformant with B.

Result: The solution is stored in X

template <typename Block0,
          typename Block1>
bool
lsqsol(
  const_Matrix<T, Block0> B,
  Matrix<T, Block1>       X)   

Description: Solves the linear least squares problem minX ||AX − B||2 for X, where A is the matrix given to the most recent call to decompose().

Requirements: The number of rows in B must be equal to the value returned by rows(). The number of rows in X must equal the value returned by columns().

Result: Stores the solution in the matrix X.

4.5.6. Singular Value Decomposition

This section describes the singular value decomposition processing object provided by VSIPL++.

4.5.6.1. Class template svd<>

The class svd performs singular value decomposition to decompose a matrix into orthogonal or unitary matrixes and singular values.

template <typename T,
          return_mechanism_type ReturnMechanism = by_value>
class svd;

Template parameters

T

The value type used for the decomposition object. May be real or complex, single- or double-precision floating-point types.

ReturnMechanism

The return mechanism type indicates whether the SV decomposition and product functions return results by value (by_value) or by reference (by_reference) into matrices provided by the caller. See section 2.7.7 for details.

4.5.6.2. Constructor

svd(
  length_type  rows, 
  length_type  columns, 
  storage_type ustorage, 
  storage_type vstorage);

Description: Constructs an svd object. The parameters rows and columns refer to the size of the matrix to be decomposed. The parameters ustorage and vstorage control how much of the U and V matrices are stored, respectively (see section 2.7.8). Note also that svd objects may also be copied (constructed) from other svd objects.

Requirements:  The number of rows and columns must both be positive.

4.5.6.3. Accessor functions

length_type rows() const;
length_type columns() const;
storage_type ustorage() const;
storage_type vstorage() const;

Description: Report the various attributes of this svd object. The number of rows is returned by rows() and the number of columns by columns(). The function ustorage() returns how the decomposition matrix U is stored by this object, if at all. The function vstorage() returns how the decomposition matrix VT or VH is stored by this object, if at all.

4.5.6.4. Solve Systems (by_value)

This function is available only if the svd class template is parameterized with ReturnMechanism=by_value

template <typename Block>
const_Vector<T, unspecified>
decompose(Matrix<T, Block> A);

Description: Performs a singular value decomposition of the matrix A containing M rows and N columns.

If T is not a specialization of complex, A = U SVT, where square orthogonal matrix U has the same number of rows as A, S is a matrix with the same shape as A and all zero values except its first p diagonal elements are real, nonincreasing, nonnegative values, and V is a square orthogonal matrix with the same number of columns as A. The number of diagonal elements are given by p = min(M, N).

If T is a specialization of complex, A = U SVH, where U, S, and V, are similar to those described above except U and V are unitary, not orthogonal, matrices.

Requirements: The matrix A must be the same size as specified in the constructor. Note that the contents may be overwritten therefore A should not be modified prior to calling any other function.

Result: Returns a vector of length p containing the sinugular values of A in non-increasing order. Note that memory may be allocated. The returned vector’s block type may be a different type from Block.

template <mat_op_type       tr,
          product_side_type ps,
          typename          Block>
const_Matrix<T, unspecified>
produ(const_Matrix<T, Block> C) const;

Description: Calculates the product of U and C. The parameter tr controls what type of operation is performed on C before the product is computed (see section 2.7.3) and ps determines what side of the product U is placed on (see section 2.7.5). The actual product and its number of rows and columns (shown in the table below) depends on the values of tr, ps, and ustorage() and whether T is not or is a specialization of complex.

For ustorage() == qrd_uvpart,

 ps == mat_lsideps == mat_rside
tr == mat_ntransUC, columns(), sCU, s, p
tr == mat_trans, TUTC, p, sCUT, s, columns()
tr == mat_herm, complex<T>UHC, p, sCUH, s, columns()

where s is an arbitrary positive value and p = min(M, N).

For ustorage() == svd_uvfull,

 ps == mat_lsideps == mat_rside
tr == mat_ntransUC, columns(), sCU, s, columns()
tr == mat_trans, TUTC, columns(), sCUT, s, columns()
tr == mat_herm, complex<T>UHC, columns(), sCUH, s, columns()

Requirements: A successful call to decompose() must have occurred for this object with ustorage() equaling either svd_uvpart or svd_uvfull. Otherwise, the behavior is undefined. The number of rows and columns (shown in the table below) of C depend on the values of tr, ps, and ustorage().

For ustorage() == svd_uvpart,

 ps == mat_lsideps == mat_rside
tr == mat_ntransp, ss, columns()
tr == mat_transcolumns(), ss, p
tr == mat_hermcolumns(), ss, p

where s is an arbitrary positive value and p = min(M, N).

For ustorage() == svd_uvfull,

 ps == mat_lsideps == mat_rside
tr == mat_ntranscolumns(), ss, columns()
tr == mat_transcolumns(), ss, columns()
tr == mat_hermcolumns(), ss, columns()

Result: Returns the product of U and C. The returned matrix’s block type may be a different type from Block.

template <mat_op_type       tr,
          product_side_type ps,
          typename          Block>
const_Matrix<T, unspecified>
prodv(const_Matrix<T, Block> C) const;

Description: Calculates the product of V and C. The parameter tr controls what type of operation is performed on C before the product is computed (see section 2.7.3) and ps determines what side of the product V is placed on (see section 2.7.5). The actual product and its number of rows and columns (shown in the table below) depends on the values of tr, ps, and vstorage() and whether T is not or is a specialization of complex.

For vstorage() == qrd_uvpart,

 ps == mat_lsideps == mat_rside
tr == mat_ntransVC, columns(), sCV, s, p
tr == mat_trans, TVTC, p, sCVT, s, columns()
tr == mat_herm, complex<T>VHC, p, sCVH, s, columns()

where s is an arbitrary positive value and p = min(M, N).

For vstorage() == svd_uvfull,

 ps == mat_lsideps == mat_rside
tr == mat_ntransVC, columns(), sCV, s, columns()
tr == mat_trans, TVTC, columns(), sCVT, s, columns()
tr == mat_herm, complex<T>VHC, columns(), sCVH, s, columns()

Requirements: A successful call to decompose() must have occurred for this object with vstorage() equaling either svd_uvpart or svd_uvfull. Otherwise, the behavior is undefined. The number of rows and columns (shown in the table below) of C depend on the values of tr, ps, and vstorage().

For vstorage() == svd_uvpart,

 ps == mat_lsideps == mat_rside
tr == mat_ntransp, ss, columns()
tr == mat_transcolumns(), ss, p
tr == mat_hermcolumns(), ss, p

where s is an arbitrary positive value and p = min(M, N).

For vstorage() == svd_uvfull,

 ps == mat_lsideps == mat_rside
tr == mat_ntranscolumns(), ss, columns()
tr == mat_transcolumns(), ss, columns()
tr == mat_hermcolumns(), ss, columns()

Result: Returns the product of V and C. The returned matrix’s block type may be a different type from Block.

const_Matrix<T, unspecified>
u(
  index_type low,
  index_type high) const;

Description: Returns consecutive columns in the matrix U from a singular value decomposition.

Requirements: A successful call to decompose() must have occurred with ustorage() equaling either svd_uvpart or svd_uvfull. Otherwise, the behavior is undefined.

Result: Returns the constant matrix U containing columns low, low+1, ..., high, inclusive.

const_Matrix<T, unspecified>
v(
  index_type low,
  index_type high) const;

Description: Returns consecutive columns in the matrix V from a singular value decomposition.

Requirements: A successful call to decompose() must have occurred with vstorage() equaling either svd_uvpart or svd_uvfull. Otherwise, the behavior is undefined.

Result: Returns the constant matrix V containing columns low, low+1, ..., high, inclusive.

4.5.6.5. Solve Systems (by_reference)

This function is available only if the svd class template is parameterized with ReturnMechanism=by_reference

template <typename Block0,
          typename Block1>
bool
decompose(
  Matrix<T, Block0> A,
  Vector<T, Block1> x);

Description: Performs a singular value decomposition of the matrix A containing M rows and N columns.

If T is not a specialization of complex, A = U SVT, where square orthogonal matrix U has the same number of rows as A, S is a matrix with the same shape as A and all zero values except its first p diagonal elements are real, nonincreasing, nonnegative values, and V is a square orthogonal matrix with the same number of columns as A. The number of diagonal elements are given by p = min(M, N).

If T is a specialization of complex, A = U SVH, where U, S, and V, are similar to those described above except U and V are unitary, not orthogonal, matrices.

Requirements: The matrix A may be overwritten. It must be the same size as specified in the constructor. The vector x must be of length p.

Result: Returns true if the decomposition succeeds. The vector x is filled with the sinugular values of A in non-increasing order. Note that memory may be allocated.

template <mat_op_type       tr,
          product_side_type ps,
          typename          Block0,
          typename          Block1>
bool
produ(
  const_Matrix<T, Block0> C,
  Matrix<T, Block1>       X) const;

Description: Calculates the product of U and C. The parameter tr controls what type of operation is performed on C before the product is computed (see section 2.7.3) and ps determines what side of the product U is placed on (see section 2.7.5). The actual product (shown in the table below) depends on the values of tr, ps, and whether T is not or is a specialization of complex.

 ps == mat_lsideps == mat_rside
tr == mat_ntransUCCU
tr == mat_trans, TUTCCUT
tr == mat_herm, complex<T>UHCCUH

Requirements: A successful call to decompose() must have occurred for this object with ustorage() equaling either svd_uvpart or svd_uvfull. Otherwise, the behavior is undefined. The number of rows and columns (shown in the table below) of C depend on the values of tr, ps, and ustorage().

For ustorage() == svd_uvpart,

 ps == mat_lsideps == mat_rside
tr == mat_ntransp, ss, rows()
tr == mat_transrows(), ss, p
tr == mat_hermrows(), ss, p

where s is an arbitrary positive value and p = min(M, N).

For ustorage() == svd_uvfull,

 ps == mat_lsideps == mat_rside
tr == mat_ntransrows(), ss, rows()
tr == mat_transrows(), ss, rows()
tr == mat_hermrows(), ss, rows()

The number of rows and columns of X (shown in the table below) depend on the values of tr, ps, and qstorage().

For ustorage() == svd_uvpart,

 ps == mat_lsideps == mat_rside
tr == mat_ntransrows(), ss, p
tr == mat_transp, ss, rows()
tr == mat_hermp, ss, rows()

where s is the same variable as above.

For ustorage() == svd_uvfull,

 ps == mat_lsideps == mat_rside
tr == mat_ntransrows(), ss, rows()
tr == mat_transrows(), ss, rows()
tr == mat_hermrows(), ss, rows()

Result: Returns true if the product succeeds. The product of U and C is stored in X.

template <mat_op_type       tr,
          product_side_type ps,
          typename          Block0,
          typename          Block1>
bool
prodv(
  const_Matrix<T, Block0> C,
  Matrix<T, Block1>       X) const;

Description: Calculates the product of V and C. The parameter tr controls what type of operation is performed on C before the product is computed (see section 2.7.3) and ps determines what side of the product V is placed on (see section 2.7.5). The actual product (shown in the table below) depends on the values of tr, ps, and whether T is not or is a specialization of complex.

For vstorage() == qrd_uvpart,

 ps == mat_lsideps == mat_rside
tr == mat_ntransVCCV
tr == mat_trans, TVTCCVT
tr == mat_herm, complex<T>VHCCVH
template <typename Block>
bool
u(
  index_type       low,
  index_type       high
  Matrix<T, Block> X) const;

Description: Returns consecutive columns in the matrix U from a singular value decomposition.

Requirements: A successful call to decompose() must have occurred with ustorage() equaling either svd_uvpart or svd_uvfull.

Result: Returns true if the matrix is stored. Stores the constant matrix U containing columns low, low+1, ..., high, inclusive in X.

template <typename Block>
bool
v(
  index_type       low,
  index_type       high
  Matrix<T, Block> X) const;

Description: Returns consecutive columns in the matrix V from a singular value decomposition.

Requirements: A successful call to decompose() must have occurred with vstorage() equaling either svd_uvpart or svd_uvfull.

Result: Returns true if the matrix is stored. Stores the constant matrix V containing columns low, low+1, ..., high, inclusive in X.

4.5.7. Toeplitz Solver

This section describes the Toeplitz linear system solver function provided by VSIPL++.

4.5.7.1. Solve Systems (return by value)

template <typename T,
          typename Block0,
          typename Block1,
          typename Block2>
const_Vector<T, unspecified>
toepsol(
  const_Vector<T, Block0> a,
  const_Vector<T, Block1> b,
  Vector<T, Block2>       w);

Description: Solve a real symmetric or complex Hermitian positive definite Toeplitz linear system Ax = b, where a specifies the Toeplitz matrix A. The Toeplitz linear system is real symmetric if type T is real and is Hermitian if type T is complex.

Requirements: The N by N Toeplitz matrix formed from a must have full rank and be positive definite. The sizes of a, b and w are equal to N. The type T may be single- or double-precision floating-point, and either real or complex. The vector w is used as a temporary workspace.

Result: The solution x is returned. The returned matrix’s block type may be a different type from Block0, Block1 or Block2.

4.5.7.2. Solve Systems (return by reference)

template <typename T,
          typename Block0,
          typename Block1,
          typename Block2,
          typename Block3>
const_Vector<T, Block3>
toepsol(
  const_Vector<T, Block0> a,
  const_Vector<T, Block1> b,
  Vector<T, Block2>       w,
  Vector<T, Block3>       x);

Description: Solve a real symmetric or complex Hermitian positive definite Toeplitz linear system Ax = b, where a specifies the Toeplitz matrix A. The Toeplitz linear system is real symmetric if type T is real and is Hermitian if type T is complex.

Requirements: The N by N Toeplitz matrix formed from a must have full rank and be positive definite. The sizes of a, b, w and x are equal to N. The type T may be single- or double-precision floating-point, and either real or complex. The vector w is used as a temporary workspace.

Result: The solution is stored in x and returned.