ivan@140: /**************************************************************************
ivan@140:  *
ivan@140:  * File name: myblas.h
ivan@140:  *
ivan@140:  * Ron Rubinstein
ivan@140:  * Computer Science Department
ivan@140:  * Technion, Haifa 32000 Israel
ivan@140:  * ronrubin@cs
ivan@140:  *
ivan@140:  * Version: 1.1
ivan@140:  * Last updated: 17.8.2009
ivan@140:  *
ivan@140:  * A collection of basic linear algebra functions, in the spirit of the
ivan@140:  * BLAS/LAPACK libraries.
ivan@140:  *
ivan@140:  *************************************************************************/
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: #ifndef __MY_BLAS_H__
ivan@140: #define __MY_BLAS_H__
ivan@140: 
ivan@140: 
ivan@140: #include "mex.h"
ivan@140: #include <math.h>
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Squared value.
ivan@140:  **************************************************************************/
ivan@140: #define SQR(X) ((X)*(X))
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Matrix-vector multiplication. 
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   y := alpha*A*x
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   A - matrix of size n X m
ivan@140:  *   x - vector of length m
ivan@140:  *   y - output vector of length n
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m - dimensions of A
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of y.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Matrix-transpose-vector multiplication. 
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   y := alpha*A'*x
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   A - matrix of size n X m
ivan@140:  *   x - vector of length n
ivan@140:  *   y - output vector of length m
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m - dimensions of A
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of y.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Sparse-matrix-vector multiplication. 
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   y := alpha*A*x
ivan@140:  *
ivan@140:  * where A is a sparse matrix.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   pr,ir,jc - sparse representation of the matrix A, of size n x m
ivan@140:  *   x - vector of length m
ivan@140:  *   y - output vector of length n
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m - dimensions of A
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of y.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Sparse-matrix-transpose-vector multiplication. 
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   y := alpha*A'*x
ivan@140:  *
ivan@140:  * where A is a sparse matrix.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   pr,ir,jc - sparse representation of the matrix A, of size n x m
ivan@140:  *   x - vector of length m
ivan@140:  *   y - output vector of length n
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m - dimensions of A
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of y.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Matrix-sparse-vector multiplication. 
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   y := alpha*A*x
ivan@140:  *
ivan@140:  * where A is a matrix and x is a sparse vector.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   A - matrix of size n X m
ivan@140:  *   pr,ir,jc - sparse representation of the vector x, of length m
ivan@140:  *   y - output vector of length n
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m - dimensions of A
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of y.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Matrix-transpose-sparse-vector multiplication. 
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   y := alpha*A'*x
ivan@140:  *
ivan@140:  * where A is a matrix and x is a sparse vector.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   A - matrix of size n X m
ivan@140:  *   pr,ir,jc - sparse representation of the vector x, of length n
ivan@140:  *   y - output vector of length m
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m - dimensions of A
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of y.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Sparse-matrix-sparse-vector multiplication. 
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   y := alpha*A*x
ivan@140:  *
ivan@140:  * where A is a sparse matrix and x is a sparse vector.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   pr,ir,jc - sparse representation of the matrix A, of size n x m
ivan@140:  *   prx,irx,jcx - sparse representation of the vector x (of length m)
ivan@140:  *   y - output vector of length n
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m - dimensions of A
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of y.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Sparse-matrix-transpose-sparse-vector multiplication. 
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   y := alpha*A'*x
ivan@140:  *
ivan@140:  * where A is a sparse matrix and x is a sparse vector.
ivan@140:  *
ivan@140:  * Importnant note: this function is provided for completeness, but is NOT efficient.
ivan@140:  * If possible, convert x to non-sparse representation and use matT_vec_sp instead.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   pr,ir,jc - sparse representation of the matrix A, of size n x m
ivan@140:  *   prx,irx,jcx - sparse representation of the vector x (of length n)
ivan@140:  *   y - output vector of length n
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m - dimensions of A
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of y.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Matrix-matrix multiplication. 
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   X := alpha*A*B
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   A - matrix of size n X m
ivan@140:  *   B - matrix of size m X k
ivan@140:  *   X - output matrix of size n X k
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m, k - dimensions of A, B
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of X.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Matrix-transpose-matrix multiplication. 
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   X := alpha*A*B
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   A - matrix of size n X m
ivan@140:  *   B - matrix of size m X k
ivan@140:  *   X - output matrix of size n X k
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m, k - dimensions of A, B
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of X.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Tensor-matrix multiplication. 
ivan@140:  *
ivan@140:  * This function accepts a 3-D tensor A of size n X m X k
ivan@140:  * and a 2-D matrix B of size l X k.
ivan@140:  * The function computes the 3-D tensor X of size n X m X l, where
ivan@140:  *
ivan@140:  *   X(i,j,:) = B*A(i,j,:)
ivan@140:  *
ivan@140:  * for all i,j.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   A - tensor of size n X m X k
ivan@140:  *   B - matrix of size l X k
ivan@140:  *   X - output tensor of size n X m X l
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m, k, l - dimensions of A, B
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of X.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Tensor-matrix-transpose multiplication. 
ivan@140:  *
ivan@140:  * This function accepts a 3-D tensor A of size n X m X k
ivan@140:  * and a 2-D matrix B of size k X l.
ivan@140:  * The function computes the 3-D tensor X of size n X m X l, where
ivan@140:  *
ivan@140:  *   X(i,j,:) = B'*A(i,j,:)
ivan@140:  *
ivan@140:  * for all i,j.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   A - tensor of size n X m X k
ivan@140:  *   B - matrix of size k X l
ivan@140:  *   X - output tensor of size n X m X l
ivan@140:  *   alpha - real constant
ivan@140:  *   n, m, k, l - dimensions of A, B
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of X.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Vector-vector sum.
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   y := alpha*x + y
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   x - vector of length n
ivan@140:  *   y - output vector of length n
ivan@140:  *   alpha - real constant
ivan@140:  *   n - length of x,y
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of y.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void vec_sum(double alpha, double x[], double y[], mwSize n);
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Vector-vector scalar multiply.
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   y := alpha* x.*y
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   x - vector of length n
ivan@140:  *   y - output vector of length n
ivan@140:  *   alpha - real constant
ivan@140:  *   n - length of x,y
ivan@140:  *
ivan@140:  * Note: This function re-writes the contents of y.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: 
ivan@140: 
ivan@140: void vec_smult(double alpha, double x[], double y[], mwSize n);
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Triangular back substitution.
ivan@140:  *
ivan@140:  * Solve the set of linear equations
ivan@140:  *
ivan@140:  *   T*x = b
ivan@140:  *
ivan@140:  * where T is lower or upper triangular.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   ul - 'U' for upper triangular, 'L' for lower triangular
ivan@140:  *   A  - matrix of size n x m containing T
ivan@140:  *   b  - vector of length k
ivan@140:  *   x  - output vector of length k
ivan@140:  *   n  - size of first dimension of A
ivan@140:  *   k  - the size of the equation set, k<=n,m
ivan@140:  *
ivan@140:  * Note:
ivan@140:  *   The matrix A can be of any size n X m, as long as n,m >= k. 
ivan@140:  *   Only the lower/upper triangle of the submatrix A(1:k,1:k) defines the
ivan@140:  *   matrix T (depending on the parameter ul).
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Solve a set of equations using a Cholesky decomposition.
ivan@140:  *
ivan@140:  * Solve the set of linear equations
ivan@140:  *
ivan@140:  *   M*x = b
ivan@140:  *
ivan@140:  * where M is positive definite with a known Cholesky decomposition:
ivan@140:  * either M=L*L' (L lower triangular) or M=U'*U (U upper triangular).
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   ul - 'U' for upper triangular, 'L' for lower triangular decomposition
ivan@140:  *   A  - matrix of size n x m with the Cholesky decomposition of M
ivan@140:  *   b  - vector of length k
ivan@140:  *   x  - output vector of length k
ivan@140:  *   n  - size of first dimension of A
ivan@140:  *   k  - the size of the equation set, k<=n,m
ivan@140:  *
ivan@140:  * Note:
ivan@140:  *   The matrix A can be of any size n X m, as long as n,m >= k. 
ivan@140:  *   Only the lower/upper triangle of the submatrix A(1:k,1:k) is used as
ivan@140:  *   the Cholesky decomposition of M (depending on the parameter ul).
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Maximum absolute value.
ivan@140:  *
ivan@140:  * Returns the index of the coefficient with maximal absolute value in a vector.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   x - vector of length n
ivan@140:  *   n - length of x
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: mwIndex maxabs(double x[], mwSize n);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Maximum vector element.
ivan@140:  *
ivan@140:  * Returns the index of the maximal coefficient in a vector.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   x - vector of length n
ivan@140:  *   n - length of x
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: mwIndex maxpos(double x[], mwSize n);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Vector-vector dot product.
ivan@140:  *
ivan@140:  * Computes an operation of the form:
ivan@140:  *
ivan@140:  *   c = a'*b
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   a, b - vectors of length n
ivan@140:  *   n - length of a,b
ivan@140:  *
ivan@140:  * Returns: The dot product c.
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: double dotprod(double a[], double b[], mwSize n);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Indexed vector assignment.
ivan@140:  *
ivan@140:  * Perform a permutation assignment of the form
ivan@140:  *
ivan@140:  *   y = x(ind)
ivan@140:  *
ivan@140:  * where ind is an array of indices to x.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   y - output vector of length k
ivan@140:  *   x - input vector of arbitrary length
ivan@140:  *   ind - array of indices into x (indices begin at 0)
ivan@140:  *   k - length of the array ind
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void vec_assign(double y[], double x[], mwIndex ind[], mwSize k);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Matrix transpose.
ivan@140:  *
ivan@140:  * Computes Y := X'
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   X - input matrix of size n X m
ivan@140:  *   Y - output matrix of size m X n
ivan@140:  *   n, m - dimensions of X
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void transpose(double X[], double Y[], mwSize n, mwSize m);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Print a matrix.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   A - matrix of size n X m
ivan@140:  *   n, m - dimensions of A
ivan@140:  *   matname - name of matrix to display
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void printmat(double A[], int n, int m, char* matname);
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /**************************************************************************
ivan@140:  * Print a sparse matrix.
ivan@140:  *
ivan@140:  * Parameters:
ivan@140:  *   A - sparse matrix of type double
ivan@140:  *   matname - name of matrix to display
ivan@140:  *
ivan@140:  **************************************************************************/
ivan@140: void printspmat(mxArray *A, char* matname);
ivan@140: 
ivan@140: 
ivan@140: #endif
ivan@140: