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