From 969d02bd210b860b1cf18bb1c8519244595fe642 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 23 Jan 2009 17:27:39 +0000 Subject: [PATCH] We can do the writing into Trilinos vectors better. Create a reference to daxpy in Lapack templates (not used anywhere right now, but can be nice to have). git-svn-id: https://svn.dealii.org/trunk@18282 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/config.h.in | 6 ++ deal.II/configure | 99 ++++++++++++++++++- deal.II/configure.in | 7 +- .../blastemplates/lapack_templates.h.in | 12 ++- deal.II/lac/include/lac/lapack_templates.h | 51 ++++++++-- .../lac/include/lac/trilinos_vector_base.h | 24 ++--- deal.II/lac/source/trilinos_vector_base.cc | 13 ++- 7 files changed, 176 insertions(+), 36 deletions(-) diff --git a/deal.II/base/include/base/config.h.in b/deal.II/base/include/base/config.h.in index eca10363c9..97b4bb871d 100644 --- a/deal.II/base/include/base/config.h.in +++ b/deal.II/base/include/base/config.h.in @@ -238,6 +238,9 @@ /* Define if the compiler provides __builtin_expect */ #undef HAVE_BUILTIN_EXPECT +/* Define to 1 if you have the `daxpy_' function. */ +#undef HAVE_DAXPY_ + /* Define to 1 if you have the `dgeevx_' function. */ #undef HAVE_DGEEVX_ @@ -342,6 +345,9 @@ /* Define to 1 if you have the header file. */ #undef HAVE_SACADO_HPP +/* Define to 1 if you have the `saxpy_' function. */ +#undef HAVE_SAXPY_ + /* Define to 1 if you have the `sgeevx_' function. */ #undef HAVE_SGEEVX_ diff --git a/deal.II/configure b/deal.II/configure index 8686662851..1aca20dff9 100755 --- a/deal.II/configure +++ b/deal.II/configure @@ -16003,7 +16003,7 @@ fi -for ac_func in dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_ +for ac_func in daxpy_ saxpy_ dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ do as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh` { echo "$as_me:$LINENO: checking for $ac_func" >&5 @@ -16102,9 +16102,104 @@ done +for ac_func in dgesvd_ sgesvd_ dgetrf_ sgetrf_ dgetri_ sgetri_ +do +as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh` +{ echo "$as_me:$LINENO: checking for $ac_func" >&5 +echo $ECHO_N "checking for $ac_func... $ECHO_C" >&6; } +if { as_var=$as_ac_var; eval "test \"\${$as_var+set}\" = set"; }; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ +/* Define $ac_func to an innocuous variant, in case declares $ac_func. + For example, HP-UX 11i declares gettimeofday. */ +#define $ac_func innocuous_$ac_func + +/* System header to define __stub macros and hopefully few prototypes, + which can conflict with char $ac_func (); below. + Prefer to if __STDC__ is defined, since + exists even on freestanding compilers. */ + +#ifdef __STDC__ +# include +#else +# include +#endif + +#undef $ac_func + +/* Override any GCC internal prototype to avoid an error. + Use char because int might match the return type of a GCC + builtin and then its argument prototype would still apply. */ +#ifdef __cplusplus +extern "C" +#endif +char $ac_func (); +/* The GNU C library defines this for functions which it implements + to always fail with ENOSYS. Some functions are actually named + something starting with __ and the normal name is an alias. */ +#if defined __stub_$ac_func || defined __stub___$ac_func +choke me +#endif + +int +main () +{ +return $ac_func (); + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext conftest$ac_exeext +if { (ac_try="$ac_link" +case "(($ac_try" in + *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;; + *) ac_try_echo=$ac_try;; +esac +eval "echo \"\$as_me:$LINENO: $ac_try_echo\"") >&5 + (eval "$ac_link") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_cxx_werror_flag" || + test ! -s conftest.err + } && test -s conftest$ac_exeext && + $as_test_x conftest$ac_exeext; then + eval "$as_ac_var=yes" +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + eval "$as_ac_var=no" +fi + +rm -f core conftest.err conftest.$ac_objext conftest_ipa8_conftest.oo \ + conftest$ac_exeext conftest.$ac_ext +fi +ac_res=`eval echo '${'$as_ac_var'}'` + { echo "$as_me:$LINENO: result: $ac_res" >&5 +echo "${ECHO_T}$ac_res" >&6; } +if test `eval echo '${'$as_ac_var'}'` = yes; then + cat >>confdefs.h <<_ACEOF +#define `echo "HAVE_$ac_func" | $as_tr_cpp` 1 +_ACEOF + +fi +done + + + -for ac_func in dgetrf_ sgetrf_ dgetri_ sgetri_ dgetrs_ sgetrs_ dstev_ sstev_ +for ac_func in dgetrs_ sgetrs_ dstev_ sstev_ do as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh` { echo "$as_me:$LINENO: checking for $ac_func" >&5 diff --git a/deal.II/configure.in b/deal.II/configure.in index 4a3a6430f5..59135b6ef5 100644 --- a/deal.II/configure.in +++ b/deal.II/configure.in @@ -9,7 +9,7 @@ dnl is stored. dnl dnl dnl Copyright: The deal.II authors -dnl 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 +dnl 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 dnl @@ -550,8 +550,9 @@ if test "x$with_blas" != "x" -a "x$with_blas" != "xno" ; then AC_CHECK_FUNC(daxpy_,,[AC_MSG_ERROR([BLAS library not installed correctly($with_blas)])]) fi -AC_CHECK_FUNCS([dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_]) -AC_CHECK_FUNCS([dgetrf_ sgetrf_ dgetri_ sgetri_ dgetrs_ sgetrs_ dstev_ sstev_]) +AC_CHECK_FUNCS([daxpy_ saxpy_ dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_]) +AC_CHECK_FUNCS([dgesvd_ sgesvd_ dgetrf_ sgetrf_ dgetri_ sgetri_]) +AC_CHECK_FUNCS([dgetrs_ sgetrs_ dstev_ sstev_]) dnl ------------------------------------------------------------- dnl set include paths of several libraries diff --git a/deal.II/contrib/blastemplates/lapack_templates.h.in b/deal.II/contrib/blastemplates/lapack_templates.h.in index c795262b03..83571349a3 100644 --- a/deal.II/contrib/blastemplates/lapack_templates.h.in +++ b/deal.II/contrib/blastemplates/lapack_templates.h.in @@ -1,3 +1,7 @@ +// vector update of the form y += alpha*x with a scalar, x,y vectors +void daxpy_ (const int* n, const double* alpha, const double* x, + const int* incx, double* y, const int* incy); + // General Matrix // Matrix vector product void dgemv_ (const char* trans, const int* m, const int* n, @@ -7,10 +11,10 @@ void dgemv_ (const char* trans, const int* m, const int* n, // Matrix matrix product void dgemm_ (const char* transa, const char* transb, - const int* m, const int* n, const int* k, - const double* alpha, const double* A, const int* lda, - const double* B, const int* ldb, - const double* beta, double* C, const int* ldc); + const int* m, const int* n, const int* k, + const double* alpha, const double* A, const int* lda, + const double* B, const int* ldb, + const double* beta, double* C, const int* ldc); // Compute LU factorization void dgetrf_ (const int* m, const int* n, double* A, diff --git a/deal.II/lac/include/lac/lapack_templates.h b/deal.II/lac/include/lac/lapack_templates.h index 71c8d40234..148ab3353c 100644 --- a/deal.II/lac/include/lac/lapack_templates.h +++ b/deal.II/lac/include/lac/lapack_templates.h @@ -23,6 +23,11 @@ using namespace dealii; extern "C" { +// vector update of the form y += alpha*x with a scalar, x,y vectors +void daxpy_ (const int* n, const double* alpha, const double* x, + const int* incx, double* y, const int* incy); +void saxpy_ (const int* n, const float* alpha, const float* x, + const int* incx, float* y, const int* incy); // General Matrix // Matrix vector product void dgemv_ (const char* trans, const int* m, const int* n, @@ -35,15 +40,15 @@ void sgemv_ (const char* trans, const int* m, const int* n, const float* b, float* y, const int* incy); // Matrix matrix product void dgemm_ (const char* transa, const char* transb, - const int* m, const int* n, const int* k, - const double* alpha, const double* A, const int* lda, - const double* B, const int* ldb, - const double* beta, double* C, const int* ldc); + const int* m, const int* n, const int* k, + const double* alpha, const double* A, const int* lda, + const double* B, const int* ldb, + const double* beta, double* C, const int* ldc); void sgemm_ (const char* transa, const char* transb, - const int* m, const int* n, const int* k, - const float* alpha, const float* A, const int* lda, - const float* B, const int* ldb, - const float* beta, float* C, const int* ldc); + const int* m, const int* n, const int* k, + const float* alpha, const float* A, const int* lda, + const float* B, const int* ldb, + const float* beta, float* C, const int* ldc); // Compute LU factorization void dgetrf_ (const int* m, const int* n, double* A, const int* lda, int* ipiv, int* info); @@ -128,6 +133,36 @@ void sstev_ (const char* jobz, const int* n, +#ifdef HAVE_DAXPY_ +inline void +axpy (const int* n, const double* alpha, const double* x, const int* incx, double* y, const int* incy) +{ + daxpy_ (n,alpha,x,incx,y,incy); +} +#else +inline void +axpy (const int*, const double*, const double*, const int*, double*, const int*) +{ + Assert (false, LAPACKSupport::ExcMissing("daxpy")); +} +#endif + + +#ifdef HAVE_SAXPY_ +inline void +axpy (const int* n, const float* alpha, const float* x, const int* incx, float* y, const int* incy) +{ + saxpy_ (n,alpha,x,incx,y,incy); +} +#else +inline void +axpy (const int*, const float*, const float*, const int*, float*, const int*) +{ + Assert (false, LAPACKSupport::ExcMissing("saxpy")); +} +#endif + + #ifdef HAVE_DGEMV_ inline void gemv (const char* trans, const int* m, const int* n, const double* alpha, const double* A, const int* lda, const double* x, const int* incx, const double* b, double* y, const int* incy) diff --git a/deal.II/lac/include/lac/trilinos_vector_base.h b/deal.II/lac/include/lac/trilinos_vector_base.h index 6bc7d2087e..268f6e2cc8 100644 --- a/deal.II/lac/include/lac/trilinos_vector_base.h +++ b/deal.II/lac/include/lac/trilinos_vector_base.h @@ -1104,22 +1104,20 @@ namespace TrilinosWrappers if (last_action != Insert) last_action = Insert; - int ierr; for (unsigned int i=0; iMap().LID(indices[i]); if (local_row == -1) { - ierr = vector->ReplaceGlobalValues (1, - (const int*)(&row), - &values[i]); + const int ierr = vector->ReplaceGlobalValues (1, + (const int*)(&row), + &values[i]); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); compressed = false; } else - ierr = vector->ReplaceMyValue (local_row, 0, values[i]); - - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + (*vector)[0][local_row] = values[i]; } } @@ -1163,22 +1161,20 @@ namespace TrilinosWrappers if (last_action != Add) last_action = Add; - int ierr; for (unsigned int i=0; iMap().LID(indices[i]); if (local_row == -1) { - ierr = vector->SumIntoGlobalValues (1, - (const int*)(&row), - &values[i]); + const int ierr = vector->SumIntoGlobalValues (1, + (const int*)(&row), + &values[i]); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); compressed = false; } else - ierr = vector->SumIntoMyValue (local_row, 0, values[i]); - - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + (*vector)[0][local_row] += values[i]; } } diff --git a/deal.II/lac/source/trilinos_vector_base.cc b/deal.II/lac/source/trilinos_vector_base.cc index 6147ec10ed..ba7b126271 100644 --- a/deal.II/lac/source/trilinos_vector_base.cc +++ b/deal.II/lac/source/trilinos_vector_base.cc @@ -39,12 +39,15 @@ namespace TrilinosWrappers // these bounds by ourselves, so // we can use []. Note that we // can only get local values. - AssertThrow (vector.in_local_range(index), - ExcAccessToNonLocalElement (index, - vector.vector->Map().MinMyGID(), - vector.vector->Map().MaxMyGID())); - return (*(vector.vector))[0][index]; + const int local_index = vector.vector->Map().LID(index); + Assert (local_index >= 0, + ExcAccessToNonLocalElement (index, + vector.vector->Map().MinMyGID(), + vector.vector->Map().MaxMyGID())); + + + return (*(vector.vector))[0][local_index]; } } -- 2.39.5