root/trunk/multiarray/dflat/Blas.d

/*==========================================================================
 * Blas.d
 *    Written in the D Programming Language (http://www.digitalmars.com/d)
 */
/***************************************************************************
 * Wrappers of BLAS operations taking dFlat types.
 *
 * <TODO: Description>
 *
 * Authors:  William V. Baxter III, OLM Digital, Inc.
 * Date: 18 Feb 2008
 * Copyright: (C) 2008  William Baxter, OLM Digital, Inc.
 *            Based in part on FLENS, Copyright (c) 2007, Michael Lehn
 *            All rights reserved.
 * License:
 *
 *   Redistribution and use in source and binary forms, with or without
 *   modification, are permitted provided that the following conditions
 *   are met:
 *
 *   1) Redistributions of source code must retain the above copyright
 *      notice, this list of conditions and the following disclaimer.
 *   2) Redistributions in binary form must reproduce the above copyright
 *      notice, this list of conditions and the following disclaimer in
 *      the documentation and/or other materials provided with the
 *      distribution.
 *   3) Neither the name of the DFLAT development group nor the names of
 *      its contributors may be used to endorse or promote products derived
 *      from this software without specific prior written permission.
 *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 *   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 *   OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 *   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 *   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 *   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 *   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
//===========================================================================

module dflat.Blas;

import dflat.Storage;
import dflat.DenseVector;
import dflat.GeneralMatrix;
import dflat.Range;
import dflat.IdMatrix;

import dfb = dflat.dflat_blas;

alias dflat.Storage.Transpose Transpose;

private template ReT(T) {
    alias typeof(T.init.re) ReT;
}
/** Some private templates to unify access to Matrices and DenseVectors 
 *  as RHSs 
 */
private int _NumCols(MT)(ref MT mat)
{
    static if(is(typeof(mat.numCols))) {
        return mat.numCols();
    }
    else static if (is(MT MTc : DenseVector!(MTc))) {
        return 1;
    }
    else static if (is(MT MTc: MTc[])) {
        return 1;
    }
    else {
        static assert(false, "No NumCols for type: "~typeof(mat).stringof);
    }
}
/// Ditto
private int _NumRows(MT)(ref MT mat)
{
    static if(is(typeof(mat.numRows))) {
        return mat.numRows();
    }
    else static if (is(MT MTc : DenseVector!(MTc))) {
        return mat.length;
    }
    else static if (is(MT : MT[])) {
        return mat.length;
    }
    else {
        static assert(false, "No NumRows for type: "~typeof(mat).stringof);
    }
}
/// Ditto
private int _RowStride(MT)(ref MT mat)
{
    static if(is(typeof(mat.strideRow))) {
        return mat.strideRow;
    }
    else static if (is(MT MTc : DenseVector!(MTc))) {
        return mat.stride();
    }
    else static if (is(MT MTc: MTc[])) {
        return 1;
    }
    else {
        static assert(false, "No RowStride for type: "~typeof(mat).stringof);
    }
}
/// Ditto
private int _LeadingDim(MT)(ref MT mat)
{
    static if(is(typeof(mat.leadingDimension))) {
        return mat.leadingDimension;
    }
    else static if (is(MT MTc : DenseVector!(MTc))) {
        return mat.length;
    }
    else static if (is(MT MTc : MTc[])) {
        return mat.length;
    }
    else {
        static assert(false, "No leading dimension for type: "~typeof(mat).stringof);
    }
}


//----------------------------------------------------------------------------
// Level 1 routines for DenseVector

X.ElementType
dot_dv(X,Y)(/*const*/ ref DenseVector!(X) x, /*const*/ ref DenseVector!(Y) y) 
{
    alias X.ElementType T;
    static if(is(T==float)||is(T==double)) {
        // x.T * y
        return dfb.dot(x.length, x.ptr, x.stride, y.ptr, y.stride);
    }
    else {
        // x.H * y
        return dfb.dotc(x.length, x.ptr, x.stride, y.ptr, y.stride);
    }
}

void copy_dv(X,Y)(/*const*/ ref DenseVector!(X) x, ref DenseVector!(Y) y)
{
//    assert(&x!=&y);
    if (y.length!=x.length) {
        y.engine().resize(x.length, x.beginIndex());
    }
    dfb.copyT(x.length, x.ptr, x.stride(), y.ptr, y.stride());
}

void scal_dv(X)(DenseVector!(X).ElementType alpha, ref DenseVector!(X) x)
{
    dfb.scalT(x.length, alpha, x.ptr, x.stride());
}

DenseVector!(I).ElementType  asum_dv(I)(/*const*/ ref DenseVector!(I) x)
{
    return dfb.asum(x.length, x.ptr, x.stride());
}

int amax_dv(I)(/*const*/ ref DenseVector!(I) x)
{
    return dfb.amax(x.length, x.ptr, x.stride()) + x.beginIndex();
}

int amin_dv(T)(/*const*/ ref DenseVector!(I) x)
{
    return dfb.amin(x.length, x.ptr, x.stride()) + x.beginIndex();
}


void axpy_dv(X,Y)(DenseVector!(X).ElementType alpha,
                  /*const*/ref DenseVector!(X) x,
                  ref DenseVector!(Y) y)
{
    assert(y.length==x.length);

    dfb.axpy(x.length, alpha, x.ptr, x.stride(), y.ptr, y.stride());
}

ReT!(X.ElementType)
nrm2_dv(X)(/*const*/ref DenseVector!(X) x)
{
    return dfb.nrm2(x.length, x.ptr, x.stride());
}

//----------------------------------------
// Level 1 for ge matrices
void scal_ge(X)(GeMatrix!(X).ElementType alpha, ref GeMatrix!(X) x)
{
    for (int i=x.beginRow(); i<x.endRow(); ++i) {
        scal_dv!(X.VectorView)( alpha, x[i,_] );
    }
}

void copy_ge(X,Y)(/*const*/ ref GeMatrix!(X) x, ref GeMatrix!(Y) y)
{
    if ((y.numRows()!=x.numRows()) || (y.numCols()!=x.numCols())) {
        y.resize(x.numRows(), x.numCols(),
                 x.beginRow(), x.beginCol());
    }

    for (int i=x.beginRow(), I=y.beginRow(); i<x.endRow(); ++i, ++I) {
        y[I,_] = x[i,_];
    }
}

void copyTrans_ge(X,Y)(/*const*/ ref GeMatrix!(X) x, ref GeMatrix!(Y) y)
{
    if ((y.numRows()!=x.numCols()) || (y.numCols()!=x.numRows())) {
        y.resize(x.numCols(), x.numRows(),
                 x.beginCol(), x.beginRow());
    }

    for (int i=x.beginRow(), I=y.beginCol(); i<x.endRow(); ++i, ++I) {
        y[_,I] = x[i,_];
    }
}

void copyConjugateTrans_ge(X,Y)(/*const*/ ref GeMatrix!(X) x, ref GeMatrix!(Y) y)
{
    if ((y.numRows()!=x.numCols()) || (y.numCols()!=x.numRows())) {
        y.resize(x.numCols(), x.numRows(),
                 x.beginCol(), x.beginRow());
    }
    // there is no blas function to calculate the conjugate transpose ...
    for (int i=x.beginRow(), I=y.beginCol(); i<x.endRow(); ++i, ++I) {
        for(int jj = x.beginCol(), JJ = y.beginRow(); jj < x.endCol(); jj++, JJ++) {
            y[JJ,I] = conjugate(x[i,jj]);
        }
    }
}

void copy_ge(Y)(/*const*/ ref IdMatrix A, ref GeMatrix!(Y) B)
{
    if (A.dim()!=B.numRows()) {
        B.resize(A.dim(),A.dim(),1,1);
    } else {
        B = 0.;
    }
    B.diag(0) = 1.;
}


void axpy_ge(X,Y)(GeMatrix!(X).ElementType alpha,
               /*const*/ ref GeMatrix!(X) x,
               ref GeMatrix!(Y) y)
{
    assert(y.numRows()==x.numRows());
    assert(y.numCols()==x.numCols());

    alias GeMatrix!(X).VectorView.Storage AX;
    alias GeMatrix!(Y).VectorView.Storage AY;


    int I = y.beginRow();
    for (int i=x.beginRow(); i<x.endRow(); ++i, ++I) {
        axpy(alpha, x[i,_], y[I,_]);
        //y[I,_] += alpha*x[i,_];
    }
}

//----------------------------------------
// Level 1 for sy matrices


void scal_sy(X)(SyMatrix!(X).ElementType alpha, ref SyMatrix!(X) x)
{
    if (x.upLo()==StorageUpLo.Upper) {
        int offset = x.beginCol() - x.beginRow();
        for (int i=x.beginRow(); i<x.endRow(); ++i) {
            for (int j=i+offset; j<x.endCol(); ++j) {
                x[i,j] *= alpha;
            }
        }
    } else {
        int offset = x.endCol() - x.endRow();
        for (int i=x.beginRow(); i<x.endRow(); ++i) {
            for (int j=x.beginCol; j<i+offset; ++j) {
                x[i,j] *= alpha;
            }
        }
    }
}


//---------------------------------------------------------------------------
// Generic versions
//   D is incapable of picking the best match from multiple templates.
//   You have to have one big template that inspects all the 
//   cases one by one in a big static if - else if - else block.

X.ElementType dot(X,Y)(/*const*/ ref X x, /*const*/ ref Y y) 
{
    alias X.ElementType XT;
    alias Y.ElementType YT;
    static if (is(X Xc: DenseVector!(Xc)) && is(Y Yc: DenseVector!(Yc))) { 
        return dot_dv!(Xc,Yc)(x,y);
    } else {
        static assert("Bad arg types for template dot()");
    }
}

ReT!(X.ElementType) nrm2(X)(/*const*/ref X x)
{
    static if (is(X Xc: DenseVector!(Xc))) { 
        return nrm2_dv!(Xc)(x);
    } else {
        static assert("Bad arg types for template dot()");
    }
}

void copy(X,Y)(/*const*/ ref X x, ref Y y)
{
    static if(is(X Xc : GeMatrix!(Xc)) && is(Y Yc : GeMatrix!(Yc))) {
        copy_ge!(Xc,Yc)(x,y);
    } else static if (is(X Xc: DenseVector!(Xc)) && is(Y Yc: DenseVector!(Yc))) { 
        copy_dv!(Xc,Yc)(x,y);
    } else {
        static assert("Bad arg types for template axpy()");
    }
}

void copyTrans(X,Y)(/*const*/ ref X x, ref Y y)
{
    static if(is(X Xc : GeMatrix!(Xc)) && is(Y Yc : GeMatrix!(Yc))) {
        copyTrans_ge!(Xc,Yc)(x,y);
    } else {
        static assert("Bad arg types for template copyTrans()");
    }
}

void copyConjugateTrans(X,Y)(/*const*/ ref X x, ref Y y)
{
    static if(is(X Xc : GeMatrix!(Xc)) && is(Y Yc : GeMatrix!(Yc))) {
        copyConjugateTrans_ge!(Xc,Yc)(x,y);
    } else {
        static assert("Bad arg types for template copyConjugateTrans()");
    }
}


void axpy(ALPHA,X,Y)(ALPHA alpha, ref X x, ref Y y) 
{
    static if(is(X Xc : GeMatrix!(Xc)) && is(Y Yc : GeMatrix!(Yc))) {
        axpy_ge!(Xc,Yc)(alpha,x,y);
    } else static if (is(X Xc: DenseVector!(Xc)) && is(Y Yc: DenseVector!(Yc))) { 
        axpy_dv!(Xc,Yc)(alpha,x,y);
    } else {
        static assert("Bad arg types for template axpy()");
    }
}

void scal(ALPHA,X)(ALPHA alpha, ref X x)
{
    static if(is(X Xc : GeMatrix!(Xc))) {
        scal_ge!(Xc)(alpha,x);
    }
    else static if(is(X Xc : SyMatrix!(Xc))) {
        scal_sy!(Xc)(alpha,x);
    }
    else static if (is(X Xc: DenseVector!(Xc))) { 
        scal_dv!(Xc)(alpha,x);
    } else {
        static assert("Bad arg types for template scal()");
    }
}


//- Level 2 --------------------------------------------------------------------

// gemv
/** matrix vector multiply
        y = alpha*A*x + beta*y
   OR   y = alpha*A.T*x + beta*y
   OR   y = alpha*A.H*x + beta*y,  with A an mxn matrix
*/
void mv(ALPHA,MA,VX,BETA,VY)(
    Transpose trans,
    ALPHA alpha, /*const*/ ref GeMatrix!(MA) A, /*const*/ ref VX x,
    BETA beta, ref VY y)
{
    assert(y.ptr !is x.ptr);
    assert(x.length==((trans==Transpose.NoTrans) 
                      ? A.numCols()
                      : A.numRows()));

    int yLength = (trans==Transpose.NoTrans) 
        ? A.numRows()
        : A.numCols();

    assert((beta==0) || (y.length==yLength));

    if (y.length!=yLength) {
        y.length = yLength;
    }

    dfb.gemv(MA.order,
             trans, A.numRows(), A.numCols(),
             alpha,
             A.ptr, _LeadingDim(A),
             x.ptr, _RowStride(x),
             beta,
             y.ptr, _RowStride(y));
}


// symv
void mv(ALPHA, MA, VX, BETA, VY)(
    ALPHA alpha, /*const*/ ref SyMatrix!(MA) A, /*const*/ ref VX x,
    BETA beta, ref VY y)
{
    assert(x.length==A.dim());
    assert((beta==0) || (y.length==A.dim));

    if (y.length!=A.dim) {
        y.resize(A.dim);
    }

    dfb.symv(MA.order,
             A.upLo, A.dim,
             alpha,
             A.ptr, _LeadingDim(A),
             x.ptr, _RowStride(x),
             beta,
             y.ptr, _RowStride(y));
}


//- Level 3 --------------------------------------------------------------------

// gemm
/// matrix matrix multiply
///     C := alpha * transa(A) * transb(B) + beta * C
void mm(ALPHA,MA,MB,BETA,MC)(
    Transpose transA, Transpose transB,
    ALPHA alpha, /*const*/ ref GeMatrix!(MA) A, /*const*/ ref GeMatrix!(MB) B,
    BETA beta, ref GeMatrix!(MC) C)
{
    assert(MA.order==MB.order);
    assert(MA.order==MC.order);

    // M: op(A) - M x K
    // N: op(B) - K x N
    version(NDEBUG) {
    } else {
        int K_A = (transA==Transpose.NoTrans) ? A.numCols() : A.numRows();
        int K_B = (transB==Transpose.NoTrans) ? B.numRows() : B.numCols();
        assert(K_A==K_B);
    }

    int m = (transA==Transpose.NoTrans) ? A.numRows() : A.numCols();
    int n = (transB==Transpose.NoTrans) ? B.numCols() : B.numRows();

    assert((beta==0) || (C.numRows()==m));
    assert((beta==0) || (C.numCols()==n));

    if ((C.numRows()!=m) || (C.numCols()!=n)) {
        C.resize(m,n);
    }

    dfb.gemm(MA.order,
             transA, transB,
             C.numRows(),
             C.numCols(),
             (transA==Transpose.NoTrans) ? A.numCols() : A.numRows(),   // K
             alpha,
             A.ptr, _LeadingDim(A),
             B.ptr, _LeadingDim(B),
             beta,
             C.ptr, _LeadingDim(C));
}


// symm
void mm(ALPHA, MA, MB, BETA, MC, dum=void)(
    BlasSide side,
    ALPHA alpha, /*const*/ ref SyMatrix!(MA) A, /*const*/ ref GeMatrix!(MB) B,
    BETA beta, ref GeMatrix!(MC) C)
{
    assert(MA.order==MB.order);
    assert(MA.order==MC.order);
    debug{ 
        if (side==BlasSide.Left) {
            assert(A.dim==B.numRows);
        } else {
            assert(B.numCols==A.dim);
        }
    }

    int m = (side==BlasSide.Left) ? A.dim() : B.numRows();
    int n = (side==BlasSide.Left) ? B.numCols() : A.dim();

    assert((beta==0) || (C.numRows()==m));
    assert((beta==0) || (C.numCols()==n));

    if ((C.numRows()!=m) || (C.numCols()!=n)) {
        C.resize(m,n);
    }

    symm(MA.order,
         side, A.upLo,
         C.numRows, C.numCols,
         alpha,
         A.ptr, A.leadingDimension,
         B.ptr, B.leadingDimension,
         beta, C.ptr, C.leadingDimension);
}

// trmm
void mm(ALPHA, MA, MB)(
    BlasSide side,
    Transpose transA, ALPHA alpha, /*const*/ ref TrMatrix!(MA) A,
    ref GeMatrix!(MB) B)
{
    assert(MA.order==MB.order);
    debug{
        if (side==BlasSide.Left) {
            assert(A.dim==B.numRows);
        } else {
            assert(B.numCols==A.dim);
        }
    }

    trmm(MA.order,
         side, A.upLo, transA,
         A.unitDiag,
         B.numRows, B.numCols,
         alpha,
         A.ptr, A.leadingDimension,
         B.ptr, B.leadingDimension);
    
}

// trsv
/** solving triangular matrix problems
        x := A.inv * x
    OR  x := A.inv.T * x
    OR  x := A.inv.H * x
 */
void trsv(E, S)(
    Transpose transpose, /*const*/ ref TrMatrix!(S) A, ref DenseVector!(E) xb)
{
    assert(xb.length==A.dim());

    trsv(S.order,                    // ORDER
         A.upLo,                     // UPLO
         transA,                     // TRANSA
         diag,                       // DIAG
         A.dim,                      // N
         A.ptr,                      // A
         A.leadingDimension,         // LDA
         xb.ptr,                     // X
         xb.stride);                 // INCX
}

// trsm
/// solving triangular matrix with multiple right hand sides
///     B := alpha * transa(A.inv) * B
/// OR  B := alpha * B * transa(A.inv)
void trsm(ALPHA, MA, MB)(
    BlasSide side,
    Transpose transA, ALPHA alpha, /*const*/ ref TrMatrix!(MA) A,
    ref GeMatrix!(MB) B)
{
    assert(MA.order==MB.order);
    debug {
        if (side==BlasSide.Left) {
            assert(A.dim==B.numRows);
        } else {
            assert(B.numCols==A.dim);
        }
    }

    trsm(MA.order,
         side, A.upLo, transA,
         A.unitDiag,
         B.numRows, B.numCols,
         alpha,
         A.ptr, A.leadingDimension,
         B.ptr, B.leadingDimension);

}



//--- Emacs setup ---
// Local Variables:
// c-basic-offset: 4
// indent-tabs-mode: nil
// End: