Ffpack_interface

---

title: Ffpack interface permalink: wiki/Ffpack_interface/ layout: wiki

FFPACK source

svn://linalg.org/fflas-ffpack/

FFPACK documentation

ffpack.h
fflas.h

Types to potentially include in M2, or at least at the engine level:

FFPACK Field types:
"modular-balanced.h":
ModularBalanced
optional: ModularBalanced optional: "modular-positive.h":
Modular
Modular "modular-int.h":
Modular

A matrix is just an array of Field::Element's (i.e. Field::Element*), stored row-wise (row-major matrix format)

General function properties

if nothing is said about the input matrix: the input matrix is modified by the called routine!
The most FFLAS routines don't utilize BLAS, but Field (Modular<>) operations.
The FFLAS routine ftrsm (triangular solve) is build on BLAS, so only routines relying on ftrsm will reach the maximal possible performance.

Functions to interface to:

FFLAS functions (wrappers for BLAS functions)

see fflas.h   for online available documentation; for recent documentation checkout the latest FFPACK trunk!

MatCopy allocates new space C and makes a copy of matrix A into the new allocated space. Leading dimension of matrix C is #cols(A) !

Field::Element* MatCopy (
const Field& F,
#rows(A),
#cols(A),
const Field::Element * matrixA,
const size_t leadingDimension(matrixA) 
)

---

level 1 BLAS :
FFLAS::fscal * Computes x = alpha* x.

static void FFLAS::fscal (const Field    F,
        size of vector X,
         const typename Field::Element alpha,
         typename Field::Element* X (dense vector of size N), const size_t incX (stride of X), 
 ) 

FFLAS::faxpy * Computes y = alpha* x+y.
static void FFLAS::fscal (const Field F,
size of vector X and Y,
const typename Field::Element alpha,
const typename Field::Element* X (dense vector of size N), const size_t incX (stride of X),
typename Field::Element* Y (dense vector of size N), const size_t incY (stride of Y),
) FFLAS::fdot * Computes x^T * y.
static typename Field::Element FFLAS::fdot (const Field F,
size of vector X and Y,
const typename Field::Element* X (dense vector of size N), const size_t incX (stride of X),
const typename Field::Element* Y (dense vector of size N), const size_t incY (stride of Y),
) optional: fswap, fcopy

---

level 2 BLAS :
FFLAS::fadd * Computes C = A + B.
static void FFLAS::fadd (const Field F,
#rows,
#cols,
const matrix A, const size_t lda,
const matrix B, const size_t ldb,
[inout] matrix C, const size_t ldc
) TODO: is it mandatory to initialize matrix C?

FFLAS::fsub * Computes C = A - B; same signature as fadd
FFLAS::faddin * Computes C = B + C;
static void FFLAS::fadd (const Field F,
#rows,
#cols,
const matrix B, const size_t ldb,
[inout] matrix C, const size_t ldc
)

FFLAS::fsubin * Computes C = B - C; same signature as faddin
FFLAS::fgemv computes Y = alpha * op(A)*X + beta * Y where op(A)=A or op(A)=transpose(A)
fgemv (const Field& F,
const FFLAS_TRANSPOSE opA (FflasNoTrans or FflasTrans),
const size_t M,
const size_t N,
const typename Field::Element alpha,
const typename Field::Element* A, const size_t leadingDimension(A),
const typename Field::Element* X (dense vector of size N), const size_t incX (stride of X),
const typename Field::Element beta,
[in][out] typename Field::Element* Y (dense vector of size M), const size_t incY (stride of Y)
)

FFLAS::fger computes A = alpha*X . Y^T + A
fger (const Field& F,
const size_t M,
const size_t N,
const typename Field::Element alpha,
const typename Field::Element* X (dense vector of size M), const size_t incX (stride of X),
const typename Field::Element* Y (dense vector of size N), const size_t incY (stride of Y)
[in][out] typename Field::Element* A (size MxN), const size_t leadingDimension(A),
)

optional interface: FFLAS::ftrsv triangular solve; computes X = (op(A)^-1) * X
ftrsv (const Field& F,
const FFLAS_FFLAS_UPLO Uplo( is A upper (FflasUpper) or lower (FflasLower) triangular?),
const FFLAS_TRANSPOSE opA (FflasNoTrans or FflasTrans),
const FFLAS_FFLAS_DIAG Diag ( A has explicit general diagonal (FflasNonUnit) or implicit unitary diagonal ( FflasUnit ),
const size_t N ( rows (op(A)),
const typename Field::Element* A (size NxN), const size_t leadingDimension(A),
[in][out] typename Field::Element* X (dense vector of size N), const size_t incX (stride of X),
)

---

level 3 BLAS :
FFLAS::fgemm computes C = alpha.op(A)*op(B) + beta.C where op(A)=A or op(A)=transpose(A)
fgemm (const Field& F,
const FFLAS_TRANSPOSE ta (FflasNoTrans or FflasTrans),
const FFLAS_TRANSPOSE tb (FflasNoTrans or FflasTrans),
const size_t m (size op(A) is m x k ),
const size_t n (size op(B) is k x n ),
const size_t k (size op(A) is m x k ),
const typename Field::Element alpha,
const typename Field::Element* A, const size_t leadingDimension(A),
const typename Field::Element* B, const size_t leadingDimension(B),
const typename Field::Element beta,
typename Field::Element* C, const size_t leadingDimension(C)
(optional) const size_t w ( recursive levels of Winograd's algorithm are used),
)

FFLAS::fsquare computes C = alpha.op(A)*op(A) + beta.C where op(A)=A or op(A)=transpose(A)
fsquare (const Field& F,
const FFLAS_TRANSPOSE ta (FflasNoTrans or FflasTrans),
const size_t n,
const typename Field::Element alpha,
const typename Field::Element* A, const size_t leadingDimension(A),
const typename Field::Element beta,
typename Field::Element* C, const size_t leadingDimension(C)
)

optional interface to ftrmm: Triangular Matrix Multiply
optional interface to ftrsm: Triangular System Solve

---

Linear Algebra Functions

see ffpack.h for online available documentation; for recent 'ffpack.h' checkout the FFPACK trunk

size_t rank =  FFPACK::Rank( 
const Field &F, 
#rows(A),  
#colsA, 
matrix A, 
leading dim of A
)

Field::Element det =  FFPACK::Det (
const Field &F, 
#rows(A),  
#cols(A),
matrix A,
leading dim of A
)

bool FFPACK::IsSingular(

const Field& F,
#rows(A),  
#cols(A),
matrix A,
leading dim of A
)

rank(A) =  FFPACK::RowRankProfile (
const Field& F, 
#rows(A),  
#cols(A),
matrix A,
leading dim of A,
[out] size_t* rkprofile ( size(rkprofile) = rank(A) )
)

rank(A) =  FFPACK::ColumnRankProfile (
const Field& F, 
#rows(A),  
#cols(A),
matrix A,
leading dim of A,
[out] size_t* rkprofile ( size(rkprofile) = rank(A) )
)

compute a permutation of A: ( P A, A P, invP (A) or A inv(P) depending on parameters 'side' and 'trans').

void  FFPACK::applyP( 
const Field& F,
const FFLAS_SIDE Side: FflasLeft for row permutation and  FflasRight for a column permutation
const FFLAS_TRANSPOSE Trans==FflasTrans for the inverse permutation P,
#rows(A),  
const int ibeg (TODO), 
const int iend (TODO),
[in][out] square matrix A,
leading dim of A,
const size_t * P ( dim(P)= #rows(A) )
)

dim( Nullspace(A) ) = FFPACK::NullSpaceBasis(
const Field &F, 
side: FFLAS::FflasRight or FFLAS::FflasLeft,
#rows(A), 
#cols(A), 
matrix A, 
leading dim of A ,
[out] NullSpace, 
[out] leadingDimension(NullSpace), 
[out] dim of nullspace 
)

type (NullSpace): pointer to Field::Element; NullSpace is allocated inside of NullSpaceBasis().

 dim( Nullspace(A) ) = FFPACK::Invert2(
const Field &F, 
#rows(A), 
matrix A to invert, leading dim (A), 
pointer to the pre-allocated storage for the inverse AInv, 
leading dim (AInv),
[out] nullspacedim 
)

'Invert2' is faster than the 'Invert' function.

Solve Ax=b
pointer to x = FFPACK::Solve (

const Field &F, 
rows (A), 
matrix A, leading dim (A),
preallocated pointer for the solution X , 
incx (stride), 
pointer to the right hand side matrix B,  
incb (stride)
);

the preallocated memory for the solution x must be initialized with zeroes.
FFPACK::fgesv solves AX=B (side=FFPACK::FflasLeft) or XA=B (side=FFPACK::FflasRight)

rank(A) = fgesv`` ``(  [in] const Field &F
[in] side, 
[in] #rows(A),[in] #cols(A), 
[in] NRHS: if side==FFPACK::FflasLeft NRHS=#cols(X)(=cols(B)); if side==FFPACK::FflasRight NRHS=#rows(X)(=rows(B)),
[in] matrix A, [in]  leading dim (A), 
[in][out] preallocated solution X matrix ,
[in] leading dim (X), 
[in] matrix B, [in] leading dim (B), 
[in][out] &info )

info: is a reference to a locally defined int. in case of success info is =0 and >0 if the system is inconsistent.
TODO:

 FFPACK::ColumnEchelonForm, FFPACK::RowEchelonForm,FFPACK::ReducedColumnEchelonForm,FFPACK::ReducedRowEchelonForm  : Documentation is for me hard to understand.
 FFPACK::CharPoly
 FFPACK::MinPoly
 FFPACK::LUdivine  Compute the LQUP factorization of the given matrix

Optional interface to

ftrtri : compute the inverse of a triangular matrix.
 ftrtrm : compute U*L of the Square Matrix A=[L\U].

Terms/abbreviations:

ld leading dimension of a matrix:
in case of a row-major matrix format (as in FFPACK) it is the number of elements in one row of the entire matrix (=#cols!)
in case of a column-major matrix format it is the number of elements in one column of the entire matrix (=#rows!)
Information about the leading dimension is mandatory due to the fact that given matrix parameters could be submatrices of a bigger matrix !