blaspp.cc

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

/** 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