blaspp.cc

LAPACK++ (Linear Algebra PACKage in C++) is a software library for numerical linear algebra that sol

// -*-C++-*- 

// Copyright (C) 2004 
// Christian Stimming <[e-mail]>

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2, or (at
// your option) any later version.

// This library 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 Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public License along
// with this library; see the file COPYING.  If not, write to the Free
// Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
// USA.

//      LAPACK++ (V. 1.1)

// requires
//

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include "lafnames.h"
#include LA_EXCEPTION_H
#include "blas1pp.h"
#include "blas2pp.h"
#include "blas3pp.h"
#include <cmath>
#include LA_VECTOR_DOUBLE_H

#include LA_BAND_MAT_DOUBLE_H
#include LA_LOWER_TRIANG_MAT_DOUBLE_H
#include LA_SPD_MAT_DOUBLE_H
#include LA_SYMM_BAND_MAT_DOUBLE_H 
#include LA_SYMM_MAT_DOUBLE_H
#include LA_SYMM_TRIDIAG_MAT_DOUBLE_H
#include LA_TRIDIAG_MAT_DOUBLE_H
#include LA_UNIT_LOWER_TRIANG_MAT_DOUBLE_H
#include LA_UNIT_UPPER_TRIANG_MAT_DOUBLE_H
#include LA_UPPER_TRIANG_MAT_DOUBLE_H

#include "blaspp.h"
#include "blas3.h"


// Only enable this when LA_NO_DEPRECATED is not defined
#ifndef LA_NO_DEPRECATED

//-------------------------------------
// Vector/Vector operators
//-------------------------------------

LaVectorDouble operator*(const LaVectorDouble &x, double a)
{
    int N = x.size();
    LaVectorDouble t(N);

    for (int i=0; i<N; i++)
    {
        t(i) = a * x(i);
    }

    return t;
}
    

/** DEPRECATED. Use the Blas functions from blas1pp.h, blas2pp.h and
 * blas3pp.h instead because they are much faster. These operators can
 * already be disabled when you #define LA_NO_DEPRECATED. */
double operator*(const LaVectorDouble &dx, 
                                const LaVectorDouble &dy)
{
    assert(dx.size()==dy.size());
    integer incx = dx.inc(), incy = dy.inc(), n = dx.size();

    return F77NAME(ddot)(&n, &dx(0), &incx, &dy(0), &incy);
}
                      
/** DEPRECATED. Use the Blas functions from blas1pp.h, blas2pp.h and
 * blas3pp.h instead because they are much faster. These operators can
 * already be disabled when you #define LA_NO_DEPRECATED. */
LaVectorDouble operator+(const LaVectorDouble &dx, 
                                const LaVectorDouble &dy)
{
    assert(dx.size()==dy.size());
    integer incx = dx.inc(), incy = dx.inc(), n = dx.size();
    double da = 1.0;

    LaVectorDouble tmp((int) n);
    tmp = dy;

    F77NAME(daxpy)(&n, &da, &dx(0), &incx, &tmp(0), &incy);
    return tmp;
}

/** DEPRECATED. Use the Blas functions from blas1pp.h, blas2pp.h and
 * blas3pp.h instead because they are much faster. These operators can
 * already be disabled when you #define LA_NO_DEPRECATED. */
LaVectorDouble operator-(const LaVectorDouble &dx, 
                                const LaVectorDouble &dy)
{
    assert(dx.size()==dy.size());
    integer incx = dx.inc(), incy = dy.inc(), n = dx.size();
    double da = -1.0;

    LaVectorDouble tmp(n);
    tmp = dx;

    F77NAME(daxpy)(&n, &da, &dy(0), &incx, &tmp(0), &incy);
    return tmp;
}
//@}

//-------------------------------------
/// @name Matrix/Vector operators (deprecated) 
//-------------------------------------
//@{

/** DEPRECATED. Use the Blas functions from blas1pp.h, blas2pp.h and
 * blas3pp.h instead because they are much faster. These operators can
 * already be disabled when you #define LA_NO_DEPRECATED. */
LaVectorDouble operator*(const LaGenMatDouble &A, 
                                const LaVectorDouble &dx)
{
    char trans = 'N';
    double alpha = 1.0, beta = 0.0;
    integer M = A.size(0), N = A.size(1), lda = A.gdim(0);

    LaVectorDouble dy(M);
    integer incx = dx.inc();
    integer incy = dy.inc();

    // dy = 0.0; -- unneeded since beta is zero

    F77NAME(dgemv)(&trans, &M, &N, &alpha, &A(0,0), &lda, &dx(0), &incx, 
        &beta, &dy(0), &incy); 
    return dy; 
        
}

#ifdef _LA_BAND_MAT_DOUBLE_H_
LaVectorDouble operator*(const LaBandMatDouble &A, 
                                const LaVectorDouble &dx)
{
    char trans = 'N';
    double alpha = 1.0, beta = 0.0;
    integer M = A.size(0), N = A.size(1), lda = A.gdim(0),
        kl = A.subdiags(), ku = A.superdiags(); 

    LaVectorDouble dy(N);
    integer incx = dx.inc(), incy = dy.inc();

    F77NAME(dgbmv)(&trans, &M, &N, &kl, &ku, &alpha, &A(0,0), &lda,
                   &dx(0), &incx, &beta, &dy(0), &incy);
    return dy;
}
#endif

#ifdef _LA_SYMM_MAT_DOUBLE_H_
LaVectorDouble operator*(const LaSymmMatDouble &A, 
                                const LaVectorDouble &dx)
{
    char uplo = 'L';
    double alpha = 1.0, beta = 0.0;
    integer N = A.size(1), lda = A.gdim(0);

    LaVectorDouble dy(N);
    integer incx = dx.inc(), incy = dy.inc();

    F77NAME(dsymv)(&uplo, &N, &alpha, &A(0,0), &lda, &dx(0), &incx,
                   &beta, &dy(0), &incy);
    return dy;
}
#endif

#ifdef _LA_SYMM_BAND_MAT_DOUBLE_H_ 
LaVectorDouble operator*(const LaSymmBandMatDouble &A, 
        const LaVectorDouble &dx) 
{
    char uplo = 'L';
    double alpha = 1.0, beta = 0.0;
    integer N = A.size(1), lda = A.gdim(0), k = A.subdiags();

    LaVectorDouble dy(N);
    integer incx = dx.inc(), incy = dy.inc();

    F77NAME(dsbmv)(&uplo, &N, &k, &alpha, &A(0,0), &lda, &dx(0), &incx,
                   &beta, &dy(0), &incy);
    return dy;
}

#endif


#ifdef _LA_SPD_MAT_DOUBLE_H_
LaVectorDouble operator*(const LaSpdMatDouble &AP, 
                                const LaVectorDouble &dx)
{
    char uplo = 'L';
    double alpha = 1.0, beta = 0.0;
    integer N = AP.size(1), incx = dx.inc(); 
    integer lda = AP.gdim(0);

    LaVectorDouble dy(N);
    integer incy = dy.inc();

    F77NAME(dsymv)(&uplo, &N, &alpha, &AP(0,0), &lda, &dx(0), &incx, &beta,
                    &dy(0), &incy);
    return dy;
}
#endif

#ifdef _LA_LOWER_TRIANG_MAT_DOUBLE_H_
LaVectorDouble operator*(const LaLowerTriangMatDouble &A, 
                                const LaVectorDouble &dx)
{
    char uplo = 'L', trans = 'N', diag = 'N';
    integer N = A.size(1), lda = A.gdim(0),
        incx = dx.inc();

    LaVectorDouble dy(dx);

    F77NAME(dtrmv)(&uplo, &trans, &diag, &N, &A(0,0), &lda, &dy(0), &incx);

    return dy;
}
#endif

#ifdef _LA_UPPER_TRIANG_MAT_DOUBLE_H_
LaVectorDouble operator*(const LaUpperTriangMatDouble &A, 
                                const LaVectorDouble &dx)
{
    char uplo = 'U', trans = 'N', diag = 'N';
    integer N = A.size(1), lda = A.gdim(0),
        incx = dx.inc();

    LaVectorDouble dy(dx);

    F77NAME(dtrmv)(&uplo, &trans, &diag, &N, &A(0,0), &lda, &dy(0), &incx);

    return dy;
}
#endif

#ifdef _LA_UNIT_LOWER_TRIANG_MAT_DOUBLE_H_
LaVectorDouble operator*(const LaUnitLowerTriangMatDouble &A,
                                const LaVectorDouble &dx)
{
    char uplo = 'L', trans = 'N', diag = 'U';
    integer N = A.size(1), lda = A.gdim(0),
        incx = dx.inc();

    LaVectorDouble dy(dx);

    F77NAME(dtrmv)(&uplo, &trans, &diag, &N, &A(0,0), &lda, &dy(0), &incx);

    return dy;
}

#endif

#ifdef _LA_UNIT_UPPER_TRIANG_MAT_DOUBLE_H_
LaVectorDouble operator*(const LaUnitUpperTriangMatDouble &A,
                                const LaVectorDouble &dx)
{
    char uplo = 'U', trans = 'N', diag = 'U';
    integer N = A.size(1), lda = A.gdim(0),
        incx = dx.inc();

    LaVectorDouble dy(dx);

    F77NAME(dtrmv)(&uplo, &trans, &diag, &N, &A(0,0), &lda, &dy(0), &incx);

    return dy;
}
#endif
//@}

//-------------------------------------
/// @name Matrix/Matrix operators (deprecated) 
//-------------------------------------
//@{