base.pure

/*  Copyright (c) 2012 by Enrico Avventi, <eavventi yahoo it>

    This file is part of pure-linalg.

    pure-linalg is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    pure-linalg is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with pure-linalg.  If not, see <http://www.gnu.org/licenses/>.*/

using namespace __cblas, __linalg;

public conj;


/**********************************************************
****            preliminary definitions                ****
**********************************************************/


/* transpose and adjoint operators */
postfix (^) ^' ^*;

dim (A::matrix^') |
dim (A::matrix^*) = m,n when n,m = dim A; end;

// indexing
A::matrix^'!(k,l) = A!(l,k);
A::matrix^*!(k,l) = conj(A!(l,k));

A::matrix^'!k = A!(k mod n,k div n) when n,m = dim A; end;
A::matrix^*!k = conj(A!(k mod n,k div n)) when n,m = dim A; end;

// slicing
A::matrix^'!!(ks,ls) = (A!!(ls,ks))^';
A::matrix^*!!(ks,ls) = (A!!(ls,ks))^*;

A::matrix^'!!ks = (transpose A)!!ks;
A::matrix^*!!ks = (transpose $ conj A)!!ks;

// these operations do not require data movement so can resolved automatically
(A::matrix^')^' = A;
v::matrix^' = rowvector v if colvectorp v;
v::matrix^' = colvector v if rowvectorp v;

(A::matrix^*)^* = A;
A::dmatrix^* = A^';
A::imatrix^* = A^'; // integer matrices are not supported by blas

// packing (will force the evaluation)
pack (A::matrix^') = transpose A;
pack (A::matrix^*) = transpose $ conj A; // FIX ME: this needs to be done in one go


/* identity matrix */

dim (eye n::int) = (n,n);

// indexing
(eye n)!(k,l) = throw out_of_bounds if k<0 || k>=n || l<0 || l>=n;
(eye n)!(k,l) = double (k==l);

(eye n)!k = throw out_of_bounds if k<0 || k>=n*n;
(eye n)!k = double (k div n == k mod n);

// slicing (will force the evaluation)
(eye n)!!(ks,ls) = {double (i==j)|i=ks; j=ls};

// transpose
(eye n)^' |
(eye n)^* = eye n;

// packing (will force the evaluation)
pack (eye n) = {double (i==j)|i=1..n; j=1..n};


/* sub-identity matrix 

subeye (n,m) k represent a n by m matrix whose 
k-th subdiagonal is composed by ones.
*/

subeye n::int k = subeye (n,n) k; 

dim (subeye (n,m) k) = (n,m);

// indexing
(subeye (n,m) j)!(k,l) = throw out_of_bounds if k<0 || k>=n || l<0 || l>=m;
(subeye (n,m) j)!(k,l) = double (k==l-j);

(subeye (n,m) j)!k = throw out_of_bounds if k<0 || k>=n*m;
(subeye (n,m) j)!k = double (k div n == (l mod n)-j);

// slicing (will force the evaluation)
(subeye (n,m) k)!!(ks,ls) = {double (i==j-k)|i=ks; j=ls};

// transpose
(subeye (n,m) k)^' |
(subeye (n,m) k)^* = subeye (m,n) (-k);

// packing (will force the evaluation)
pack (subeye (n,m) k) = {double (i==j-k)|i=1..n; j=1..m};


/* constant matrix

cstmat (n,m) alpha represent a n by m matrix
whose elements are all equat to alpha.
*/

// shortcut for a matrix of ones
ones _dim = cstmat _dim 1.0;

cstmat n::int alpha = cstmat (n,n) alpha;

dim (cstmat (n,m) k) = (n,m);

// indexing
(cstmat (n,m) alpha)!(k,l) = throw out_of_bounds if k<0 || k>=n || l<0 || l>=m;
(cstmat (n,m) alpha)!(k,l) = alpha;

(cstmat (n,m) alpha)!k = throw out_of_bounds if k<0 || k>=n*m;
(cstmat (n,m) alpha)!k = alpha;

// slicing (will force the evaluation)
(cstmat (n,m) alpha::double)!!(ks,ls) = {alpha|i=1..#ks; j=1..#ls};

// transpose
(cstmat (n,m) alpha)^' = cstmat (m,n) alpha;
(cstmat (n,m) alpha)^* = cstmat (m,n) conj k;

// packing (will force the evaluation)
pack (cstmat (n,m) alpha::double) = {alpha|i=1..n; j=1..m};


/* zero matrix */

/*zeros n::int = zeros (n,n);

dim (zeros (n,m)) = (n,m);

// indexing
(zeros (n,m))!(k,l) = throw out_of_bounds if k<0 || k>=n || l<0 || l>=m;
(zeros (n,m))!(k,l) = 0.0;

(zeros (n,m))!k = throw out_of_bounds if k<0 || k>=n*m;
(zeros (n,m))!k = 0.0;

// slicing (will force the evaluation)
(zeros (n,m))!!(ks,ls) = dmatrix (#ks,#ls);

// transpose
(zeros (n,m))^' |
(zeros (n,m))^* = zeros (m,n);

// packing (will force the evaluation)
pack (zeros (n,m)) = dmatrix (n,m);*/


/* definition of extended and lazy matrix types. 

extended matrices are the union of numerical matrices and 
lazy matrices */

lmatrixp A::dmatrix = false;

lmatrixp (A::dmatrix^*) |
lmatrixp (A::dmatrix^') = true;

lmatrixp (eye n::int) = true;

lmatrixp (cstmat (n::int,m::int) alpha::double) = true;

/*lmatrixp (zeros (n::int,m::int)) = true;*/

type lmatrix = lmatrixp;

xmatrixp A::lmatrix = true;
xmatrixp A::dmatrix = true;

type xmatrix = xmatrixp;

/**********************************************************
****                  blas wrappers                    ****
**********************************************************/


namespace __cblas;
using "lib:libblas";

// getting various enums values for controlling blas routines
using "lib:cblas_enums";

// row or colum major enums
private extern void* row_major();
private extern void* col_major();

const CblasRowMajor = row_major;
const CblasColMajor = col_major;

// transpose and adjoint enums
private extern void* no_trans();
private extern void* trans();
private extern void* conj_trans();

const CblasNoTrans = no_trans;
const CblasTrans = trans;
const CblasConjTrans = conj_trans;

/* daxpy */
extern void cblas_daxpy(int n,double alpha,void* x,int incx,void* y,int incy);

/* dgemv */
extern void cblas_dgemv(void* order, void* transA, int M, int N, double alpha,
    double* A, int lda, double* x, int incx, double beta, double* y, int incy);

/* dgemm */
extern void cblas_dgemm(void* order, void* transA, void* transB, int M, int N, 
    int K, double alpha, double* A, int lda, double* B, int ldb, double beta, 
    double* C, int ldc);
namespace;



/**********************************************************
****          inplace function definitions             ****
**********************************************************/

/* vector-vector inplace multiplication and summation

TODO
*/


/* matrix-vector inplace multiplication and summation

    alpha*M*x + y -> y

where M is an xmatrix and alpha is a scalar.
*/

inplace_mulsum alpha::double A::dmatrix x::dmatrix y::dmatrix = ()
    when
        mA, nA = dim A;
        _ = cblas_dgemv CblasRowMajor CblasNoTrans mA nA alpha A nA x 1 1.0 y 1;
    end
    if dim x!1 == 1 && dim A!1 == dim x!0 && dim y!1 == 1 && dim y!0 == dim A!0;

inplace_mulsum alpha::double (A::dmatrix^') x::dmatrix y::dmatrix = ()
    when
        mA, nA = dim A;
        _ = cblas_dgemv CblasRowMajor CblasTrans mA nA alpha A nA x 1 1.0 y 1;
    end 
    if dim x!1 == 1 && dim A!0 == dim x!0 && dim y!1 == 1 && dim y!0 == dim A!1;


/* vector-matrix inplace multiplication and summation

    alpha*x*M + y -> y

where M is an xmatrix and alpha is a skalar.
*/

inplace_mulsum alpha::double x::dmatrix A::dmatrix y::dmatrix = ()
    when
        mA, nA = dim A;
        _ = cblas_dgemv CblasRowMajor CblasTrans mA nA alpha A nA x 1 1.0 y 1;
    end
    if dim x!0 == 1 && dim A!0 == dim x!1 && dim y!0 == 1 && dim y!1 == dim A!1;

inplace_mulsum alpha::double x::dmatrix (A::dmatrix^') y::dmatrix = ()
    when
        mA, nA = dim A;
        _ = cblas_dgemv CblasRowMajor CblasNoTrans mA nA alpha A nA x 1 1.0 y 1;
    end
    if dim x!0 == 1 && dim A!1 == dim x!1 && dim y!0 == 1 && dim y!1 == dim A!0;


/* matrix-matrix inplace multiplication and summation

    alpha*H*F + C -> C

where H and F are xmatrix and alpha is a scalar.
*/

inplace_mulsum alpha::double A::dmatrix B::dmatrix C::dmatrix = ()
    when
        mA, nA = dim A;
        mB, nB = dim B;
        _ = cblas_dgemm CblasRowMajor CblasNoTrans CblasNoTrans mA nB nA alpha A nA B nB 1.0 C nB;
    end
    if dim A!1 == dim B!0 && dim A!0 == dim C!0 && dim B!1 == dim C!1;

inplace_mulsum alpha::double (A::dmatrix^') B::dmatrix C::dmatrix = ()
    when
        mA, nA = dim A;
        mB, nB = dim B;
        _ = cblas_dgemm CblasRowMajor CblasTrans CblasNoTrans nA nB mA alpha A nA B nB 1.0 C nB;
    end 
    if dim A!0 == dim B!0 && dim A!1 == dim C!0 && dim B!1 == dim C!1;

inplace_mulsum alpha::double A::dmatrix (B::dmatrix^') C::dmatrix = ()
    when
        mA, nA = dim A;
        mB, nB = dim B;
        _ = cblas_dgemm CblasRowMajor CblasNoTrans CblasTrans mA mB nA alpha A nA B nB 1.0 C mB;
    end
    if dim A!1 == dim B!1 && dim A!0 == dim C!0 && dim B!0 == dim C!1;

inplace_mulsum alpha::double (A::dmatrix^') (B::dmatrix^') C::dmatrix = ()
    when
        mA, nA = dim A;
        mB, nB = dim B;
        _ = cblas_dgemm CblasRowMajor CblasTrans CblasTrans nA mB mA alpha A nA B nB 1.0 C mB;
    end
    if dim A!0 == dim B!1 && dim A!1 == dim C!0 && dim B!0 == dim C!1;


/* multiplication with identity resolves into a inplace summation */

inplace_mulsum alpha::double A::xmatrix (eye n::int) C::dmatrix = inplace_sum alpha A C
    if dim A!1 == n; 

inplace_mulsum alpha::double (eye n::int) A::xmatrix C::dmatrix = inplace_sum alpha A C
    if dim A!0 == n;


/* multiplication with a constant matrix */

// matrix-vector

inplace_mulsum alpha::double A::xmatrix (cstmat (n,m) beta) C::dmatrix = 
    inplace_mulsum alpha A B C
    when
        B = pack (cstmat (n,m) beta);
    end
    if dim A!1 == n && m == 1 && dim C!0 == dim A!0 && dim C!1 == 1; 

inplace_mulsum alpha::double (cstmat (n,m) beta) A::xmatrix C::dmatrix = 
    inplace_sum alpha B A C
    when
        B = pack (cstmat (n,m) beta);
    end
    if dim A!0 == m && n == 1 && dim C!0 == 1 && dim C!1 == dim A!1;

/* matrix-matrix, resolve using the following

alpha*cstmat (n,m) b = (cstmat (n,1) alpha) * (cstmat (1.m) b) 
*/

/*inplace_mulsum alpha::double A::xmatrix (cstmat (n::int,m::int) b::double) C::dmatrix = 
    inplace_mulsum 1.0 AB D C
    when
        B = pack (cstmat (n,1) alpha);
        AB = dmatrix (dim A!0,1);
        _ = inplace_mulsum 1.0 A B AB;
        D = pack (cstmat (1,m) b);
    end
    if dim A!1 == n && dim C!0 == dim A!0 && dim C!1 == m; 

inplace_mulsum alpha::double (cstmat (n::int,m::int) b::double) A::xmatrix C::dmatrix = 
    inplace_sum 1.0 D BA C
    when
        B = pack (cstmat (1,m) alpha);
        BA = dmatrix (1,dim A!1);
        _ = inplace_mulsum 1.0 B A BA;
        D = pack (cstmat (n,1) b);
    end
    if dim A!0 == m && dim C!0 == n && dim C!1 == dim A!1;*/


/* multiplication with zero matrix */

/*inplace_mulsum alpha::double (zeros (n::int,m::int)) A::xmatrix C::dmatrix = ()
    if dim A!0 == m && dim C!0 == n && dim C!1 == dim A!1;

inplace_mulsum alpha::double A::xmatrix (zeros (n::int,m::int)) C::dmatrix = ()
    if dim A!1 == n && dim C!0 == dim A!0 && dim C!1 == m;*/

// catch all error
inplace_mulsum alpha A::xmatrix B::xmatrix C::dmatrix = throw dimensions_mismatch;


/* inplace summation */

inplace_sum alpha::double A::dmatrix B::dmatrix = ()
    when
        _ = cblas_daxpy (#A) alpha A 1 B 1;
    end
    if dim A == dim B;

// maybe a version without temporary is needed
inplace_sum alpha::double (A::dmatrix^') B::dmatrix = ()
    when
        C = transpose A;
        _ = cblas_daxpy (#B) alpha C 1 B 1;
    end
    if dim (A^') == dim B;

inplace_sum alpha::double (eye n::int) B::dmatrix = ()
    when 
        // really need to have a better way to create constant matrices
        C = {1.0| i=1..n}; 
        _ = cblas_daxpy (#C) alpha C 1 B (n+1);
    end
    if dim B == (n,n);

// last resort, just pack the lmatrix
inplace_sum alpha::double A::lmatrix B::dmatrix = ()
    when 
         C = pack A; 
        _ = cblas_daxpy (#B) alpha C 1 B 1;
    end
    if dim A == dim B;

// catch all error
inplace_sum alpha A::xmatrix B::dmatrix = throw dimensions_mismatch;