// -*-C++-*-
// Copyright (C) 2004
// Christian Stimming <stimming@tuhh.de>
// 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)
//-------------------------------------
//@{
/** 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. */
LaGenMatDouble operator*(const LaGenMatDouble &A,
const LaGenMatDouble &B)
{
char t = 'N';
integer m = A.size(0), k = A.size(1), n = B.size(1);
integer lda = A.gdim(0), ldb = B.gdim(0);
double alpha = 1.0, beta = 0.0;
LaGenMatDouble C(m,n);
integer ldc = A.gdim(0);
//C = 0.0; -- beta is zero, doesn't need to be set
F77NAME(dgemm)(&t, &t, &m, &n, &k, &alpha, &A(0,0), &lda, &B(0,0), &ldb,
&beta, &C(0,0), &ldc);
return C.shallow_assign();
}
#ifdef _LA_UNIT_LOWER_TRIANG_MAT_DOUBLE_H_
LaGenMatDouble operator*(const LaUnitLowerTriangMatDouble &A,
const LaGenMatDouble &B)
{
char side = 'L', uplo = 'L', transa = 'N', diag = 'U';
double alpha = 1.0;
integer m = B.size(0), n = B.size(1),
lda = A.gdim(0), ldb = B.gdim(0);
LaGenMatDouble C(B);
F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
&A(0,0), &lda, &C(0,0), &ldb);
return C;
}
#endif
#ifdef _LA_UNIT_UPPER_TRIANG_MAT_DOUBLE_H_
LaGenMatDouble operator*(const LaUnitUpperTriangMatDouble &A,
const LaGenMatDouble &B)
{
char side = 'L', uplo = 'U', transa = 'N', diag = 'U';
double alpha = 1.0;
integer m = B.size(0), n = B.size(1),
lda = A.gdim(0), ldb = B.gdim(0);
LaGenMatDouble C(B);
F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
&A(0,0), &lda, &C(0,0), &ldb);
return C;
}
#endif
#ifdef _LA_LOWER_TRIANG_MAT_DOUBLE_H_
LaGenMatDouble operator*(const LaLowerTriangMatDouble &A,
const LaGenMatDouble &B)
{
char side = 'L', uplo = 'L', transa = 'N', diag = 'N';
double alpha = 1.0;
integer m = B.size(0), n = B.size(1),
lda = A.gdim(0), ldb = B.gdim(0);
LaGenMatDouble C(B);
F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
&A(0,0), &lda, &C(0,0), &ldb);
return C;
}
#endif
#ifdef _LA_UPPER_TRIANG_MAT_DOUBLE_H_
LaGenMatDouble operator*(const LaUpperTriangMatDouble &A,
const LaGenMatDouble &B)
{
char side = 'L', uplo = 'U', transa = 'N', diag = 'N';
double alpha = 1.0;
integer m = B.size(0), n = B.size(1),
lda = A.gdim(0), ldb = B.gdim(0);
LaGenMatDouble C(B);
F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
&A(0,0), &lda, &C(0,0), &ldb);
return C;
}
#endif
#ifdef _LA_SYMM_MAT_DOUBLE_H_
LaGenMatDouble operator*(const LaSymmMatDouble &A,
const LaGenMatDouble &B)
{
char side = 'L', uplo = 'L';
double alpha = 1.0, beta = 0.0;
LaGenMatDouble C(B.size(1),A.size(1));
integer m = C.size(0), n = C.size(1), lda = A.gdim(0),
ldb = B.gdim(0), ldc = C.gdim(0);
F77NAME(dsymm)(&side, &uplo, &m, &n, &alpha, &A(0,0), &lda,
&B(0,0), &ldb, &beta, &C(0,0), &ldc);
return C;
}
#endif
#ifdef _LA_SYMM_TRIDIAG_MAT_DOUBLE_H_
LaVectorDouble operator*(const LaSymmTridiagMatDouble& A,
const LaVectorDouble& X)
{
integer M = A.size();
integer N = X.size();
LaVectorDouble R(M);
R(0) = ((A.diag(0)(0) * X(0)) + (A.diag(-1)(0) * X(1)));
for (integer i = 1; i < M-2; i++)
{
R(i) = ((A.diag(-1)(i-1) * X(i-1)) +
(A.diag(0)(i) * X(i)) +
(A.diag(-1)(i) * X(i+1)));
}
R(M-1) = ((A.diag(0)(M-1) * X(N-1)) + (A.diag(-1)(M-2) * X(N-2)));
return R;
}
#endif
#ifdef _LA_TRIDIAG_MAT_DOUBLE_H_
LaVectorDouble operator*(const LaTridiagMatDouble& A,
const LaVectorDouble& X)
{
integer M = A.size();
integer N = X.size();
LaVectorDouble R(M);
R(0) = ((A.diag(0)(0) * X(0)) + (A.diag(-1)(0) * X(1)));
for (integer i = 1; i < M-2; i++)
{
R(i) = ((A.diag(-1)(i-1) * X(i-1)) +
(A.diag(0)(i) * X(i)) +
(A.diag(1)(i) * X(i+1)));
}
R(M-1) = ((A.diag(0)(M-1) * X(N-1)) + (A.diag(1)(M-2) * X(N-2)));
return R;
}
#endif
/** 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. */
LaGenMatDouble operator-(const LaGenMatDouble &A,
const LaGenMatDouble &B)
{
#ifndef HPPA
const char fname[] = "operator+(A,B)";
#else
char *fname = NULL;
#endif
integer M = A.size(0);
integer N = A.size(1);
if (M != B.size(0) || N != B.size(1))
{
throw(LaException(fname, "matrices non-conformant."));
}
LaGenMatDouble C(M,N);
// slow mode
// we'll hook the BLAS in later
for (integer i=0; i<M; i++)
for(integer j=0; j<N; j++)
C(i,j) = A(i,j) - B(i,j);
return C;
}
/** 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. */
LaGenMatDouble operator+(const LaGenMatDouble &A,
const LaGenMatDouble &B)
{
#ifndef HPPA
const char fname[] = "operator+(A,B)";
#else
char *fname = NULL;
#endif
integer M = A.size(0);
integer N = A.size(1);
if (M != B.size(0) || N != B.size(1))
{
throw(LaException(fname, "matrices non-conformant."));
}
LaGenMatDouble C(M,N);
// slow mode
// we'll hook the BLAS in later
for (integer i=0; i<M; i++)
for(integer j=0; j<N; j++)
C(i,j) = A(i,j) + B(i,j);
return C;
}
# ifdef LA_COMPLEX_SUPPORT
LaGenMatComplex operator+(const LaGenMatComplex &A,
const LaGenMatComplex &B)
{
#ifndef HPPA
const char fname[] = "operator+(A,B)";
#else
char *fname = NULL;
#endif
integer M = A.size(0);
integer N = A.size(1);
if (M != B.size(0) || N != B.size(1))
{
throw(LaException(fname, "matrices non-conformant."));
}
LaGenMatComplex C(M,N);
// slow mode
// we'll hook the BLAS in later
for (integer i=0; i<M; i++)
for(integer j=0; j<N; j++)
C(i,j) = LaComplex(A(i,j)) + LaComplex(B(i,j));
return C;
}
/** 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. */
LaGenMatComplex operator-(const LaGenMatComplex &A,
const LaGenMatComplex &B)
{
#ifndef HPPA
const char fname[] = "operator+(A,B)";
#else
char *fname = NULL;
#endif
integer M = A.size(0);
integer N = A.size(1);
if (M != B.size(0) || N != B.size(1))
{
throw(LaException(fname, "matrices non-conformant."));
}
LaGenMatComplex C(M,N);
// slow mode
// we'll hook the BLAS in later
for (integer i=0; i<M; i++)
for(integer j=0; j<N; j++)
C(i,j) = LaComplex(A(i,j)) - LaComplex(B(i,j));
return C;
}
//@}
# endif // LA_COMPLEX_SUPPORT
#endif // LA_NO_DEPRECATED
syntax highlighted by Code2HTML, v. 0.9.1