blas3pp.cc

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

#ifdef HAVE_CONFIG_H
# include <config.h>
#endif

#include "lafnames.h"
#include LA_EXCEPTION_H
#include LA_VECTOR_DOUBLE_H
#include LA_SYMM_MAT_DOUBLE_H
#include LA_UNIT_UPPER_TRIANG_MAT_DOUBLE_H
#include LA_UPPER_TRIANG_MAT_DOUBLE_H
#include LA_UNIT_LOWER_TRIANG_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_TRIDIAG_MAT_DOUBLE_H
#include LA_BAND_MAT_DOUBLE_H
#include LA_SYMM_TRIDIAG_MAT_DOUBLE_H

#ifdef LA_COMPLEX_SUPPORT
#include "f2c.h"
#include "lapackc.h"
#include LA_VECTOR_COMPLEX_H
#include LA_GEN_MAT_COMPLEX_H
#endif

#include "blas3pp.h"
#include "blas3.h"
#include "blas1.h"
#include "blas1pp.h" // for Blas_Norm2
#include "laindex.h"


void Blas_Mat_Mat_Mult(const LaGenMatDouble &A, 
		       const LaGenMatDouble &B, LaGenMatDouble &C, 
		       bool transpose_A, bool transpose_B,
		       double alpha, double beta)
{
   char transa = transpose_A ? 'T' : 'N';
   char transb = transpose_B ? 'T' : 'N';
   // m, n, k have to be determined according to op(A), op(B)!
   integer m = transpose_A ? A.size(1) : A.size(0);
   integer k = transpose_A ? A.size(0) : A.size(1);
   integer n = transpose_B ? B.size(0) : B.size(1);
   integer lda = A.gdim(0), ldb = B.gdim(0), ldc = C.gdim(0);
   if (alpha != 0.0)
   {
      assert(m == C.size(0));
      assert(n == C.size(1));
      assert(k == (transpose_B ? B.size(1) : B.size(0)) );
   }

   F77NAME(dgemm)(&transa, &transb, &m, &n, &k,
		  &alpha, &A(0,0), &lda, &B(0,0), &ldb,
		  &beta, &C(0,0), &ldc);
}

void Blas_Mat_Mat_Mult(const LaGenMatDouble &A, 
            const LaGenMatDouble &B, LaGenMatDouble &C, 
            double alpha , double beta )
{
   Blas_Mat_Mat_Mult(A, B, C, false, false, alpha, beta);
}

void Blas_Mat_Trans_Mat_Mult(const LaGenMatDouble &A, 
            const LaGenMatDouble &B, LaGenMatDouble &C, 
            double alpha, double beta )
{
   Blas_Mat_Mat_Mult(A, B, C, true, false, alpha, beta);
}

void Blas_Mat_Mat_Trans_Mult(const LaGenMatDouble &A, 
            const LaGenMatDouble &B, LaGenMatDouble &C, 
            double alpha,  double beta )
{
   Blas_Mat_Mat_Mult(A, B, C, false, true, alpha, beta);
}



double my_Dot_Prod(const LaGenMatDouble &dx, const LaGenMatDouble &dy)
{
  assert(dx.size(0)*dx.size(1)==dy.size(0)*dy.size(1));
  integer n = dx.size(0)*dx.size(1);
  integer incx = dx.inc(0)*dx.inc(1);
  integer incy = dy.inc(0)*dy.inc(1);

  return F77NAME(ddot)(&n, &dx(0,0), &incx, &dy(0,0), &incy);
}

void Blas_Mat_Mat_Mult(const LaGenMatDouble &A, 
		       const LaGenMatDouble &B, LaVectorDouble &C)
{
  // calculate only the diagonal of A times B
  int msize = std::min(A.size(0), B.size(1));
  assert(A.size(1) == B.size(0));
  assert(C.size() >= msize);
  for (int i=0; i < msize; i++)
    C(i) = my_Dot_Prod( A.row(i), B.col(i) );
}

void Blas_Mat_Trans_Mat_Mult(const LaGenMatDouble &A, 
			     const LaGenMatDouble &B, LaVectorDouble &C)
{
  // calculate only the diagonal of A times B
  int msize = std::min(A.size(0), B.size(1));
  assert(A.size(0) == B.size(0));
  assert(C.size() >= msize);
  for (int i=0; i < msize; i++)
    C(i) = my_Dot_Prod( A.col(i), B.col(i) );
}

void Blas_Mat_Mat_Trans_Mult(const LaGenMatDouble &A, 
			     const LaGenMatDouble &B, LaVectorDouble &C)
{
  // calculate only the diagonal of A times B
  int msize = std::min(A.size(0), B.size(1));
  assert(A.size(1) == B.size(1));
  assert(C.size() >= msize);
  for (int i=0; i < msize; i++)
    C(i) = my_Dot_Prod( A.row(i), B.row(i) );
}






#ifdef _LA_UNIT_LOWER_TRIANG_MAT_DOUBLE_H_

void Blas_Mat_Mat_Solve(LaUnitLowerTriangMatDouble &A, 
            LaGenMatDouble &B, double alpha)
{
        char side = 'L', uplo = 'L', transa = 'N', diag = 'U';
        integer m = B.size(0), n = B.size(1), 
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrsm)(&side, &uplo, &transa, &diag, &m, &n, &alpha, 
                &A(0,0), &lda, &B(0,0), &ldb);
}

#endif

#ifdef _LA_UNIT_UPPER_TRIANG_MAT_DOUBLE_H_
void Blas_Mat_Mat_Mult(LaUnitUpperTriangMatDouble &A,
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'U', transa = 'N', diag = 'U';
        integer m = B.size(0), n = B.size(1),
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
                &A(0,0), &lda, &B(0,0), &ldb);
}

void Blas_Mat_Mat_Solve(LaUnitUpperTriangMatDouble &A, 
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'U', transa = 'N', diag = 'U';
        integer m = B.size(0), n = B.size(1), 
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrsm)(&side, &uplo, &transa, &diag, &m, &n, &alpha, 
                &A(0,0), &lda, &B(0,0), &ldb);
}

#endif

#ifdef _LA_LOWER_TRIANG_MAT_DOUBLE_H_
void Blas_Mat_Mat_Mult(LaLowerTriangMatDouble &A,
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'L', transa = 'N', diag = 'N';
        integer m = B.size(0), n = B.size(1),
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
                &A(0,0), &lda, &B(0,0), &ldb);
}

void Blas_Mat_Mat_Solve(LaLowerTriangMatDouble &A, 
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'L', transa = 'N', diag = 'N';
        integer m = B.size(0), n = B.size(1), 
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrsm)(&side, &uplo, &transa, &diag, &m, &n, &alpha, 
                &A(0,0), &lda, &B(0,0), &ldb);
}
#endif


#ifdef _LA_UPPER_TRIANG_MAT_DOUBLE_H_
void Blas_Mat_Mat_Mult(LaUpperTriangMatDouble &A,
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'U', transa = 'N', diag = 'N';
        integer m = B.size(0), n = B.size(1),
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
                &A(0,0), &lda, &B(0,0), &ldb);
}

void Blas_Mat_Mat_Solve(LaUpperTriangMatDouble &A, 
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'U', transa = 'N', diag = 'N';
        integer m = B.size(0), n = B.size(1), 
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrsm)(&side, &uplo, &transa, &diag, &m, &n, &alpha, 
                &A(0,0), &lda, &B(0,0), &ldb);
}
#endif


#ifdef _LA_UNIT_LOWER_TRIANG_MAT_DOUBLE_H_
void Blas_Mat_Trans_Mat_Solve(LaUnitLowerTriangMatDouble &A,
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'L', transa = 'T', diag = 'U';
        integer m = B.size(0), n = B.size(1),
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrsm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
                &A(0,0), &lda, &B(0,0), &ldb);
}
#endif

#ifdef _LA_UNIT_UPPER_TRIANG_MAT_DOUBLE_H_
void Blas_Mat_Trans_Mat_Solve(LaUnitUpperTriangMatDouble &A,
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'U', transa = 'T', diag = 'U';
        integer m = B.size(0), n = B.size(1),
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrsm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
                &A(0,0), &lda, &B(0,0), &ldb);
}
#endif

#ifdef LA_LOWER_TRIANG_MAT_DOUBLE_H
void Blas_Mat_Mat_Mult(LaUnitLowerTriangMatDouble &A,
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'L', transa = 'N', diag = 'U';
        integer m = B.size(0), n = B.size(1),
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
                &A(0,0), &lda, &B(0,0), &ldb);
}

void Blas_Mat_Trans_Mat_Solve(LaLowerTriangMatDouble &A,
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'L', transa = 'T', diag = 'N';
        integer m = B.size(0), n = B.size(1),
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrsm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
                &A(0,0), &lda, &B(0,0), &ldb);
}
#endif


#ifdef _LA_UPPER_TRIANG_MAT_DOUBLE_H 
void Blas_Mat_Trans_Mat_Solve(LaUpperTriangMatDouble &A,
            LaGenMatDouble &B, double alpha )
{
        char side = 'L', uplo = 'U', transa = 'T', diag = 'N';
        integer m = B.size(0), n = B.size(1),
                lda = A.gdim(0), ldb = B.gdim(0);

  F77NAME(dtrsm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
                &A(0,0), &lda, &B(0,0), &ldb);
}

#endif

#ifdef _LA_SYMM_MAT_DOUBLE_H_
void Blas_Mat_Mat_Mult(LaSymmMatDouble &A, LaGenMatDouble &B, 
		       LaGenMatDouble &C, 
		       double alpha , double beta,
		       bool b_left_side )
{
  char side = b_left_side ? 'L' : 'R';
  
  if (side=='L')
  {
    assert ( B.size(1)==C.size(0) && A.size(0)==B.size(0) && A.size(0)==C.size(1) );
  }
  else
  {
    assert ( B.size(0)==C.size(1) && A.size(0)==B.size(1) && A.size(0)==C.size(0) );
  }
  
  char uplo = 'L';
  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);
}

void Blas_R1_Update(LaSymmMatDouble &C, LaGenMatDouble &A,
		    double alpha , double beta, bool right_transposed )
{
  char trans = right_transposed ? 'N' : 'T';

  char uplo = 'L';
  integer n = C.size(0), k,
          lda = A.gdim(0), ldc = C.gdim(0);

  if (trans=='N')
  {
    k = A.size(1);
    assert ( A.size(0)==n );
  }
  else
  {
    k = A.size(0);
    assert ( A.size(1)==n );
  }

  F77NAME(dsyrk)(&uplo, &trans, &n, &k, &alpha, &A(0,0), &lda,
                &beta, &C(0,0), &ldc);
}

void Blas_R2_Update(LaSymmMatDouble &C, LaGenMatDouble &A,
            LaGenMatDouble &B, double alpha , double beta, bool right_transposed )
{
  char trans = right_transposed ? 'N' : 'T';
  char uplo = 'L';
  integer n = C.size(0), k, lda = A.gdim(0),
          ldb = B.gdim(0), ldc = C.gdim(0);

  if (trans=='N')
  {
    k = A.size(1);
    assert( A.size(0)==n && B.size(0)==n && B.size(1)==k);
  }
  else
  {
    k = A.size(0);
    assert( A.size(1)==n && B.size(1)==n && B.size(0)==k);
  }

  F77NAME(dsyr2k)(&uplo, &trans, &n, &k, &alpha, &A(0,0), &lda,
                &B(0,0), &ldb, &beta, &C(0,0), &ldc);
}

#endif


#ifdef LA_COMPLEX_SUPPORT