diff --git a/Makefile.am b/Makefile.am index 9cd354b62017949467a3083f6fb7e0ccda1db3ce..5552c44b775c92d8c83003c98eca7b826c54ab56 100644 --- a/Makefile.am +++ b/Makefile.am @@ -84,6 +84,7 @@ endif EXTRA_libelpa@SUFFIX@_private_la_DEPENDENCIES = \ src/elpa1/elpa_reduce_add_vectors.F90 \ src/elpa1/elpa_transpose_vectors.F90 \ + src/elpa_api_math_template.F90 \ src/elpa_impl_template.F90 \ src/elpa1/elpa1_compute_template.F90 \ src/elpa2/elpa2_compute_real_template.F90 \ @@ -657,6 +658,7 @@ EXTRA_DIST = \ manual_cpp \ nvcc_wrap \ src/GPU/cuUtils_template.cu \ + src/elpa_api_math_template.F90 \ src/elpa_impl_template.F90 \ src/elpa1/elpa1_compute_template.F90 \ src/elpa1/elpa1_merge_systems_real_template.F90 \ diff --git a/src/elpa_api.F90 b/src/elpa_api.F90 index 52377232111e4f09ee7b39b5f1ed210068bf07d6..e4ea3747f28823aca9528dea4db701a75699a831 100644 --- a/src/elpa_api.F90 +++ b/src/elpa_api.F90 @@ -506,992 +506,39 @@ module elpa_api ! Actual math routines - !> \brief abstract definition of interface to solve double real eigenvalue problem - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a double real matrix a: defines the problem to solve - !> \param ev double real: on output stores the eigenvalues - !> \param q double real matrix q: on output stores the eigenvalues - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_eigenvectors_d_i(self, a, ev, q, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - real(kind=c_double) :: a(self%local_nrows, *), q(self%local_nrows,*) -#else - real(kind=c_double) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols) -#endif - real(kind=c_double) :: ev(self%na) - -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to solve single real eigenvalue problem - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a single real matrix a: defines the problem to solve - !> \param ev single real: on output stores the eigenvalues - !> \param q single real matrix q: on output stores the eigenvalues - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_eigenvectors_f_i(self, a, ev, q, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - real(kind=c_float) :: a(self%local_nrows, *), q(self%local_nrows, *) -#else - real(kind=c_float) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols) -#endif - real(kind=c_float) :: ev(self%na) -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to solve double complex eigenvalue problem - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a double complex matrix a: defines the problem to solve - !> \param ev double real: on output stores the eigenvalues - !> \param q double complex matrix q: on output stores the eigenvalues - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_eigenvectors_dc_i(self, a, ev, q, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self - -#ifdef USE_ASSUMED_SIZE - complex(kind=c_double_complex) :: a(self%local_nrows, *), q(self%local_nrows, *) -#else - complex(kind=c_double_complex) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols) -#endif - real(kind=c_double) :: ev(self%na) -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to solve single complex eigenvalue problem - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a single complex matrix a: defines the problem to solve - !> \param ev single real: on output stores the eigenvalues - !> \param q single complex matrix q: on output stores the eigenvalues - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_eigenvectors_fc_i(self, a, ev, q, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - complex(kind=c_float_complex) :: a(self%local_nrows, *), q(self%local_nrows, *) -#else - complex(kind=c_float_complex) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols) -#endif - real(kind=c_float) :: ev(self%na) -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to solve double real eigenvalue problem - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a double real matrix a: defines the problem to solve - !> \param ev double real: on output stores the eigenvalues - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_eigenvalues_d_i(self, a, ev, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - real(kind=c_double) :: a(self%local_nrows, *) -#else - real(kind=c_double) :: a(self%local_nrows, self%local_ncols) -#endif - real(kind=c_double) :: ev(self%na) -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to solve single real eigenvalue problem - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a single real matrix a: defines the problem to solve - !> \param ev single real: on output stores the eigenvalues - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_eigenvalues_f_i(self, a, ev, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - real(kind=c_float) :: a(self%local_nrows, *) -#else - real(kind=c_float) :: a(self%local_nrows, self%local_ncols) -#endif - real(kind=c_float) :: ev(self%na) -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to solve double complex eigenvalue problem - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a double complex matrix a: defines the problem to solve - !> \param ev double real: on output stores the eigenvalues - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_eigenvalues_dc_i(self, a, ev, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self - -#ifdef USE_ASSUMED_SIZE - complex(kind=c_double_complex) :: a(self%local_nrows, *) -#else - complex(kind=c_double_complex) :: a(self%local_nrows, self%local_ncols) -#endif - real(kind=c_double) :: ev(self%na) -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to solve single complex eigenvalue problem - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a single complex matrix a: defines the problem to solve - !> \param ev single real: on output stores the eigenvalues - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_eigenvalues_fc_i(self, a, ev, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - complex(kind=c_float_complex) :: a(self%local_nrows, *) -#else - complex(kind=c_float_complex) :: a(self%local_nrows, self%local_ncols) -#endif - real(kind=c_float) :: ev(self%na) -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - !> \brief abstract definition of interface to solve double real generalized eigenvalue problem - !> - !> The dimensions of the matrix a and b (locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a double real matrix a: defines the problem to solve - !> \param b double real matrix b: defines the problem to solve - !> \param ev double real: on output stores the eigenvalues - !> \param q double real matrix q: on output stores the eigenvalues - !> \param is_already_decomposed logical, input: is it repeated call with the same b (decomposed in the fist call)? - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_generalized_eigenvectors_d_i(self, a, b, ev, q, is_already_decomposed, error) - use iso_c_binding - use elpa_constants - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - real(kind=c_double) :: a(self%local_nrows, *), b(self%local_nrows, *), q(self%local_nrows, *) -#else - real(kind=c_double) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols), & - q(self%local_nrows, self%local_ncols) -#endif - real(kind=c_double) :: ev(self%na) - logical :: is_already_decomposed - integer, optional :: error - end subroutine - end interface - - !> \brief abstract definition of interface to solve single real generalized eigenvalue problem - !> - !> The dimensions of the matrix a and b(locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a single real matrix a: defines the problem to solve - !> \param b single real matrix b: defines the problem to solve - !> \param ev single real: on output stores the eigenvalues - !> \param q single real matrix q: on output stores the eigenvalues - !> \param is_already_decomposed logical, input: is it repeated call with the same b (decomposed in the fist call)? - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_generalized_eigenvectors_f_i(self, a, b, ev, q, is_already_decomposed, error) - use iso_c_binding - use elpa_constants - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - real(kind=c_float) :: a(self%local_nrows, *), b(self%local_nrows, *), q(self%local_nrows, *) -#else - real(kind=c_float) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols), & - q(self%local_nrows, self%local_ncols) -#endif - real(kind=c_float) :: ev(self%na) - logical :: is_already_decomposed - - integer, optional :: error - end subroutine - end interface - - !> \brief abstract definition of interface to solve double complex generalized eigenvalue problem - !> - !> The dimensions of the matrix a and b(locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a double complex matrix a: defines the problem to solve - !> \param b double complex matrix b: defines the problem to solve - !> \param ev double real: on output stores the eigenvalues - !> \param q double complex matrix q: on output stores the eigenvalues - !> \param is_already_decomposed logical, input: is it repeated call with the same b (decomposed in the fist call)? - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_generalized_eigenvectors_dc_i(self, a, b, ev, q, is_already_decomposed, error) - use iso_c_binding - use elpa_constants - import elpa_t - implicit none - class(elpa_t) :: self - -#ifdef USE_ASSUMED_SIZE - complex(kind=c_double_complex) :: a(self%local_nrows, *), b(self%local_nrows, *), q(self%local_nrows, *) -#else - complex(kind=c_double_complex) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols), & - q(self%local_nrows, self%local_ncols) -#endif - real(kind=c_double) :: ev(self%na) - logical :: is_already_decomposed - - integer, optional :: error - end subroutine - end interface - - !> \brief abstract definition of interface to solve single complex generalized eigenvalue problem - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution - !> blocksize, the number of eigenvectors - !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> It is possible to change the behaviour of the method by setting tunable parameters with the - !> class method "set" - !> Parameters - !> \details - !> \param self class(elpa_t), the ELPA object - !> \param a single complex matrix a: defines the problem to solve - !> \param b single complex matrix b: defines the problem to solve - !> \param ev single real: on output stores the eigenvalues - !> \param q single complex matrix q: on output stores the eigenvalues - !> \param is_already_decomposed logical, input: is it repeated call with the same b (decomposed in the fist call)? - !> \result error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_generalized_eigenvectors_fc_i(self, a, b, ev, q, is_already_decomposed, error) - use iso_c_binding - use elpa_constants - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - complex(kind=c_float_complex) :: a(self%local_nrows, *), b(self%local_nrows, *), q(self%local_nrows, *) -#else - complex(kind=c_float_complex) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols), & - q(self%local_nrows, self%local_ncols) -#endif - real(kind=c_float) :: ev(self%na) - logical :: is_already_decomposed - - integer, optional :: error - end subroutine - end interface - - - !> \brief abstract definition of interface to compute C : = A**T * B for double real matrices - !> where A is a square matrix (self%a,self%na) which is optionally upper or lower triangular - !> B is a (self%na,ncb) matrix - !> C is a (self%na,ncb) matrix where optionally only the upper or lower - !> triangle may be computed - !> - !> the MPI commicators are already known to the type. Thus the class method "setup" must be called - !> BEFORE this method is used - !> \details - !> - !> \param self class(elpa_t), the ELPA object - !> \param uplo_a 'U' if A is upper triangular - !> 'L' if A is lower triangular - !> anything else if A is a full matrix - !> Please note: This pertains to the original A (as set in the calling program) - !> whereas the transpose of A is used for calculations - !> If uplo_a is 'U' or 'L', the other triangle is not used at all, - !> i.e. it may contain arbitrary numbers - !> \param uplo_c 'U' if only the upper diagonal part of C is needed - !> 'L' if only the upper diagonal part of C is needed - !> anything else if the full matrix C is needed - !> Please note: Even when uplo_c is 'U' or 'L', the other triangle may be - !> written to a certain extent, i.e. one shouldn't rely on the content there! - !> \param ncb Number of columns of global matrices B and C - !> \param a matrix a - !> \param self%local_nrows number of rows of local (sub) matrix a, set with method set("local_nrows,value") - !> \param self%local_ncols number of columns of local (sub) matrix a, set with method set("local_ncols,value") - !> \param b matrix b - !> \param nrows_b number of rows of local (sub) matrix b - !> \param ncols_b number of columns of local (sub) matrix b - !> \param nblk blocksize of cyclic distribution, must be the same in both directions! - !> \param c matrix c - !> \param nrows_c number of rows of local (sub) matrix c - !> \param ncols_c number of columns of local (sub) matrix c - !> \param error optional argument, error code which can be queried with elpa_strerr - abstract interface - subroutine elpa_hermitian_multiply_d_i (self,uplo_a, uplo_c, ncb, a, b, nrows_b, ncols_b, & - c, nrows_c, ncols_c, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self - character*1 :: uplo_a, uplo_c - integer(kind=c_int), intent(in) :: nrows_b, ncols_b, nrows_c, ncols_c, ncb -#ifdef USE_ASSUMED_SIZE - real(kind=c_double) :: a(self%local_nrows,*), b(nrows_b,*), c(nrows_c,*) -#else - real(kind=c_double) :: a(self%local_nrows,self%local_ncols), b(nrows_b,ncols_b), c(nrows_c,ncols_c) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to compute C : = A**T * B - !> where A is a square matrix (self%na,self%na) which is optionally upper or lower triangular - !> B is a (self%na,ncb) matrix - !> C is a (self%na,ncb) matrix where optionally only the upper or lower - !> triangle may be computed - !> - !> the MPI commicators are already known to the type. Thus the class method "setup" must be called - !> BEFORE this method is used - !> \details - !> - !> \param self class(elpa_t), the ELPA object - !> \param uplo_a 'U' if A is upper triangular - !> 'L' if A is lower triangular - !> anything else if A is a full matrix - !> Please note: This pertains to the original A (as set in the calling program) - !> whereas the transpose of A is used for calculations - !> If uplo_a is 'U' or 'L', the other triangle is not used at all, - !> i.e. it may contain arbitrary numbers - !> \param uplo_c 'U' if only the upper diagonal part of C is needed - !> 'L' if only the upper diagonal part of C is needed - !> anything else if the full matrix C is needed - !> Please note: Even when uplo_c is 'U' or 'L', the other triangle may be - !> written to a certain extent, i.e. one shouldn't rely on the content there! - !> \param ncb Number of columns of global matrices B and C - !> \param a matrix a - !> \param self%local_nrows number of rows of local (sub) matrix a, set with method set("local_nrows",value) - !> \param self%local_ncols number of columns of local (sub) matrix a, set with method set("local_nrows",value) - !> \param b matrix b - !> \param nrows_b number of rows of local (sub) matrix b - !> \param ncols_b number of columns of local (sub) matrix b - !> \param nblk blocksize of cyclic distribution, must be the same in both directions! - !> \param c matrix c - !> \param nrows_c number of rows of local (sub) matrix c - !> \param ncols_c number of columns of local (sub) matrix c - !> \param error optional argument, error code which can be queried with elpa_strerr - abstract interface - subroutine elpa_hermitian_multiply_f_i (self,uplo_a, uplo_c, ncb, a, b, nrows_b, ncols_b, & - c, nrows_c, ncols_c, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self - character*1 :: uplo_a, uplo_c - integer(kind=c_int), intent(in) :: nrows_b, ncols_b, nrows_c, ncols_c, ncb -#ifdef USE_ASSUMED_SIZE - real(kind=c_float) :: a(self%local_nrows,*), b(nrows_b,*), c(nrows_c,*) -#else - real(kind=c_float) :: a(self%local_nrows,self%local_ncols), b(nrows_b,ncols_b), c(nrows_c,ncols_c) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to compute C : = A**H * B - !> where A is a square matrix (self%na,self%a) which is optionally upper or lower triangular - !> B is a (self%na,ncb) matrix - !> C is a (self%na,ncb) matrix where optionally only the upper or lower - !> triangle may be computed - !> - !> the MPI commicators are already known to the type. Thus the class method "setup" must be called - !> BEFORE this method is used - !> \details - !> - !> \param self class(elpa_t), the ELPA object - !> \param uplo_a 'U' if A is upper triangular - !> 'L' if A is lower triangular - !> anything else if A is a full matrix - !> Please note: This pertains to the original A (as set in the calling program) - !> whereas the transpose of A is used for calculations - !> If uplo_a is 'U' or 'L', the other triangle is not used at all, - !> i.e. it may contain arbitrary numbers - !> \param uplo_c 'U' if only the upper diagonal part of C is needed - !> 'L' if only the upper diagonal part of C is needed - !> anything else if the full matrix C is needed - !> Please note: Even when uplo_c is 'U' or 'L', the other triangle may be - !> written to a certain extent, i.e. one shouldn't rely on the content there! - !> \param ncb Number of columns of global matrices B and C - !> \param a matrix a - !> \param self%local_nrows number of rows of local (sub) matrix a, set with the method set("local_nrows",value) - !> \param self%local_ncols number of columns of local (sub) matrix a, set with the method set("local_ncols",value) - !> \param b matrix b - !> \param nrows_b number of rows of local (sub) matrix b - !> \param ncols_b number of columns of local (sub) matrix b - !> \param nblk blocksize of cyclic distribution, must be the same in both directions! - !> \param c matrix c - !> \param nrows_c number of rows of local (sub) matrix c - !> \param ncols_c number of columns of local (sub) matrix c - !> \param error optional argument, error code which can be queried with elpa_strerr - abstract interface - subroutine elpa_hermitian_multiply_dc_i (self,uplo_a, uplo_c, ncb, a, b, nrows_b, ncols_b, & - c, nrows_c, ncols_c, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self - character*1 :: uplo_a, uplo_c - integer(kind=c_int), intent(in) :: nrows_b, ncols_b, nrows_c, ncols_c, ncb -#ifdef USE_ASSUMED_SIZE - complex(kind=c_double_complex) :: a(self%local_nrows,*), b(nrows_b,*), c(nrows_c,*) -#else - complex(kind=c_double_complex) :: a(self%local_nrows,self%local_ncols), b(nrows_b,ncols_b), c(nrows_c,ncols_c) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to compute C : = A**H * B - !> where A is a square matrix (self%na,self%na) which is optionally upper or lower triangular - !> B is a (self%na,ncb) matrix - !> C is a (self%na,ncb) matrix where optionally only the upper or lower - !> triangle may be computed - !> - !> the MPI commicators are already known to the type. Thus the class method "setup" must be called - !> BEFORE this method is used - !> \details - !> - !> \param self class(elpa_t), the ELPA object - !> \param uplo_a 'U' if A is upper triangular - !> 'L' if A is lower triangular - !> anything else if A is a full matrix - !> Please note: This pertains to the original A (as set in the calling program) - !> whereas the transpose of A is used for calculations - !> If uplo_a is 'U' or 'L', the other triangle is not used at all, - !> i.e. it may contain arbitrary numbers - !> \param uplo_c 'U' if only the upper diagonal part of C is needed - !> 'L' if only the upper diagonal part of C is needed - !> anything else if the full matrix C is needed - !> Please note: Even when uplo_c is 'U' or 'L', the other triangle may be - !> written to a certain extent, i.e. one shouldn't rely on the content there! - !> \param ncb Number of columns of global matrices B and C - !> \param a matrix a - !> \param self%local_nrows number of rows of local (sub) matrix a, set with class method set("local_nrows",value) - !> \param self%local_ncols number of columns of local (sub) matrix a, set with class method set("local_ncols",value) - !> \param b matrix b - !> \param nrows_b number of rows of local (sub) matrix b - !> \param ncols_b number of columns of local (sub) matrix b - !> \param nblk blocksize of cyclic distribution, must be the same in both directions! - !> \param c matrix c - !> \param nrows_c number of rows of local (sub) matrix c - !> \param ncols_c number of columns of local (sub) matrix c - !> \param error optional argument, error code which can be queried with elpa_strerr - abstract interface - subroutine elpa_hermitian_multiply_fc_i (self, uplo_a, uplo_c, ncb, a, b, nrows_b, ncols_b, & - c, nrows_c, ncols_c, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self - character*1 :: uplo_a, uplo_c - integer(kind=c_int), intent(in) :: nrows_b, ncols_b, nrows_c, ncols_c, ncb -#ifdef USE_ASSUMED_SIZE - complex(kind=c_float_complex) :: a(self%local_nrows,*), b(nrows_b,*), c(nrows_c,*) -#else - complex(kind=c_float_complex) :: a(self%local_nrows,self%local_ncols), b(nrows_b,ncols_b), c(nrows_c,ncols_c) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to do a cholesky decomposition of a double real matrix - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution - !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> Parameters - !> \param self class(elpa_t), the ELPA object - !> \param a double real matrix: the matrix to be decomposed - !> \param error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_cholesky_d_i (self, a, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - real(kind=c_double) :: a(self%local_nrows,*) -#else - real(kind=c_double) :: a(self%local_nrows,self%local_ncols) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to do a cholesky decomposition of a single real matrix - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution - !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> Parameters - !> \param self class(elpa_t), the ELPA object - !> \param a single real matrix: the matrix to be decomposed - !> \param error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_cholesky_f_i(self, a, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - real(kind=c_float) :: a(self%local_nrows,*) -#else - real(kind=c_float) :: a(self%local_nrows,self%local_ncols) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to do a cholesky decomposition of a double complex matrix - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution - !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> Parameters - !> \param self class(elpa_t), the ELPA object - !> \param a double complex matrix: the matrix to be decomposed - !> \param error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_cholesky_dc_i (self, a, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - complex(kind=c_double_complex) :: a(self%local_nrows,*) -#else - complex(kind=c_double_complex) :: a(self%local_nrows,self%local_ncols) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to do a cholesky decomposition of a single complex matrix - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution - !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> Parameters - !> \param self class(elpa_t), the ELPA object - !> \param a single complex matrix: the matrix to be decomposed - !> \param error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_cholesky_fc_i (self, a, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - complex(kind=c_float_complex) :: a(self%local_nrows,*) -#else - complex(kind=c_float_complex) :: a(self%local_nrows,self%local_ncols) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to invert a triangular double real matrix - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution - !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> Parameters - !> \param self class(elpa_t), the ELPA object - !> \param a double real matrix: the matrix to be inverted - !> \param error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_invert_trm_d_i (self, a, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - real(kind=c_double) :: a(self%local_nrows,*) -#else - real(kind=c_double) :: a(self%local_nrows,self%local_ncols) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to invert a triangular single real matrix - !> Parameters - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution - !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> \param self class(elpa_t), the ELPA object - !> \param a single real matrix: the matrix to be inverted - !> \param error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_invert_trm_f_i (self, a, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - real(kind=c_float) :: a(self%local_nrows,*) -#else - real(kind=c_float) :: a(self%local_nrows,self%local_ncols) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to invert a triangular double complex matrix - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution - !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> Parameters - !> \param self class(elpa_t), the ELPA object - !> \param a double complex matrix: the matrix to be inverted - !> \param error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_invert_trm_dc_i (self, a, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - complex(kind=c_double_complex) :: a(self%local_nrows,*) -#else - complex(kind=c_double_complex) :: a(self%local_nrows,self%local_ncols) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to invert a triangular single complex matrix - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution - !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> Parameters - !> \param self class(elpa_t), the ELPA object - !> \param a single complex matrix: the matrix to be inverted - !> \param error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_invert_trm_fc_i (self, a, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self -#ifdef USE_ASSUMED_SIZE - complex(kind=c_float_complex) :: a(self%local_nrows,*) -#else - complex(kind=c_float_complex) :: a(self%local_nrows,self%local_ncols) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to solve the eigenvalue problem for a double-precision real valued tridiangular matrix - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution - !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> Parameters - !> \param self class(elpa_t), the ELPA object - !> \param d double real 1d array: the diagonal elements of a matrix defined in setup, on output the eigenvalues - !> in ascending order - !> \param e double real 1d array: the subdiagonal elements of a matrix defined in setup - !> \param q double real matrix: on output contains the eigenvectors - !> \param error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_solve_tridiagonal_d_i (self, d, e, q, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self - real(kind=c_double) :: d(self%na), e(self%na) -#ifdef USE_ASSUMED_SIZE - real(kind=c_double) :: q(self%local_nrows,*) -#else - real(kind=c_double) :: q(self%local_nrows,self%local_ncols) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - - - !> \brief abstract definition of interface to solve the eigenvalue problem for a single-precision real valued tridiangular matrix - !> - !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution - !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE - !> with the class method "setup" - !> - !> Parameters - !> \param self class(elpa_t), the ELPA object - !> \param d single real 1d array: the diagonal elements of a matrix defined in setup, on output the eigenvalues - !> in ascending order - !> \param e single real 1d array: the subdiagonal elements of a matrix defined in setup - !> \param q single real matrix: on output contains the eigenvectors - !> \param error integer, optional : error code, which can be queried with elpa_strerr - abstract interface - subroutine elpa_solve_tridiagonal_f_i (self, d, e, q, error) - use iso_c_binding - import elpa_t - implicit none - class(elpa_t) :: self - real(kind=c_float) :: d(self%na), e(self%na) -#ifdef USE_ASSUMED_SIZE - real(kind=c_float) :: q(self%local_nrows,*) -#else - real(kind=c_float) :: q(self%local_nrows,self%local_ncols) -#endif -#ifdef USE_FORTRAN2008 - integer, optional :: error -#else - integer :: error -#endif - end subroutine - end interface - +#define REALCASE 1 +#define DOUBLE_PRECISION 1 +#include "general/precision_macros.h" +#include "elpa_api_math_template.F90" +#undef REALCASE +#undef DOUBLE_PRECISION + +#ifdef WANT_SINGLE_PRECISION_REAL +#define REALCASE 1 +#define SINGLE_PRECISION 1 +#include "general/precision_macros.h" +#include "elpa_api_math_template.F90" +#undef REALCASE +#undef SINGLE_PRECISION +#endif /* WANT_SINGLE_PRECISION_REAL */ + +#define COMPLEXCASE 1 +#define DOUBLE_PRECISION 1 +#include "general/precision_macros.h" +#include "elpa_api_math_template.F90" +#undef DOUBLE_PRECISION +#undef COMPLEXCASE + +#ifdef WANT_SINGLE_PRECISION_COMPLEX +#define COMPLEXCASE 1 +#define SINGLE_PRECISION +#include "general/precision_macros.h" +#include "elpa_api_math_template.F90" +#undef COMPLEXCASE +#undef SINGLE_PRECISION +#endif /* WANT_SINGLE_PRECISION_COMPLEX */ + +! end of math routines !> \brief abstract definition of interface to destroy an ELPA object !> Parameters diff --git a/src/elpa_api_math_template.F90 b/src/elpa_api_math_template.F90 new file mode 100644 index 0000000000000000000000000000000000000000..058c6316aed8679e5340f932cc5c85900cc545a9 --- /dev/null +++ b/src/elpa_api_math_template.F90 @@ -0,0 +1,374 @@ + + !> \brief abstract definition of interface to solve double real eigenvalue problem + !> + !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution + !> blocksize, the number of eigenvectors + !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE + !> with the class method "setup" + !> + !> It is possible to change the behaviour of the method by setting tunable parameters with the + !> class method "set" + !> Parameters + !> \details + !> \param self class(elpa_t), the ELPA object +#if ELPA_IMPL_SUFFIX == d + !> \param a double real matrix a: defines the problem to solve + !> \param ev double real: on output stores the eigenvalues + !> \param q double real matrix q: on output stores the eigenvalues +#endif +#if ELPA_IMPL_SUFFIX == f + !> \param a single real matrix a: defines the problem to solve + !> \param ev single real: on output stores the eigenvalues + !> \param q single real matrix q: on output stores the eigenvalues +#endif +#if ELPA_IMPL_SUFFIX == dc + !> \param a double complex matrix a: defines the problem to solve + !> \param ev double real: on output stores the eigenvalues + !> \param q double complex matrix q: on output stores the eigenvalues +#endif +#if ELPA_IMPL_SUFFIX ==fc + !> \param a single complex matrix a: defines the problem to solve + !> \param ev single real: on output stores the eigenvalues + !> \param q single complex matrix q: on output stores the eigenvalues +#endif + !> \result error integer, optional : error code, which can be queried with elpa_strerr + abstract interface + subroutine elpa_eigenvectors_& + &ELPA_IMPL_SUFFIX& + &_i(self, a, ev, q, error) + use iso_c_binding + import elpa_t + implicit none + class(elpa_t) :: self + +#ifdef USE_ASSUMED_SIZE + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows, *), q(self%local_nrows,*) +#else + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows, self%local_ncols), q(self%local_nrows, self%local_ncols) +#endif + real(kind=C_REAL_DATATYPE) :: ev(self%na) + +#ifdef USE_FORTRAN2008 + integer, optional :: error +#else + integer :: error +#endif + end subroutine + end interface + + + + !> \brief abstract definition of interface to solve a eigenvalue problem + !> + !> The dimensions of the matrix a (locally ditributed and global), the block-cyclic distribution + !> blocksize, the number of eigenvectors + !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE + !> with the class method "setup" + !> + !> It is possible to change the behaviour of the method by setting tunable parameters with the + !> class method "set" + !> Parameters + !> \details + !> \param self class(elpa_t), the ELPA object +#if ELPA_IMPL_SUFFIX == d + !> \param a double real matrix a: defines the problem to solve + !> \param ev double real: on output stores the eigenvalues +#endif +#if ELPA_IMPL_SUFFIX == f + !> \param a single real matrix a: defines the problem to solve + !> \param ev single real: on output stores the eigenvalues +#endif +#if ELPA_IMPL_SUFFIX == dc + !> \param a double complex matrix a: defines the problem to solve + !> \param ev double real: on output stores the eigenvalues +#endif +#if ELPA_IMPL_SUFFIX == fc + !> \param a single complex matrix a: defines the problem to solve + !> \param ev single real: on output stores the eigenvalues +#endif + + !> \result error integer, optional : error code, which can be queried with elpa_strerr + abstract interface + subroutine elpa_eigenvalues_& + &ELPA_IMPL_SUFFIX& + &_i(self, a, ev, error) + use iso_c_binding + import elpa_t + implicit none + class(elpa_t) :: self +#ifdef USE_ASSUMED_SIZE + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows, *) +#else + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows, self%local_ncols) +#endif + real(kind=C_REAL_DATATYPE) :: ev(self%na) + +#ifdef USE_FORTRAN2008 + integer, optional :: error +#else + integer :: error +#endif + end subroutine + end interface + + + + !> \brief abstract definition of interface to solve a generalized eigenvalue problem + !> + !> The dimensions of the matrix a and b (locally ditributed and global), the block-cyclic distribution + !> blocksize, the number of eigenvectors + !> to be computed and the MPI communicators are already known to the object and MUST be set BEFORE + !> with the class method "setup" + !> + !> It is possible to change the behaviour of the method by setting tunable parameters with the + !> class method "set" + !> Parameters + !> \details + !> \param self class(elpa_t), the ELPA object +#if ELPA_IMPL_SUFFIX == d + !> \param a double real matrix a: defines the problem to solve + !> \param b double real matrix b: defines the problem to solve + !> \param ev double real: on output stores the eigenvalues + !> \param q double real matrix q: on output stores the eigenvalues +#endif +#if ELPA_IMPL_SUFFIX == f + !> \param a single real matrix a: defines the problem to solve + !> \param b single real matrix b: defines the problem to solve + !> \param ev single real: on output stores the eigenvalues + !> \param q single real matrix q: on output stores the eigenvalues +#endif +#if ELPA_IMPL_SUFFIX == dc + !> \param a double complex matrix a: defines the problem to solve + !> \param b double complex matrix b: defines the problem to solve + !> \param ev double real: on output stores the eigenvalues + !> \param q double complex matrix q: on output stores the eigenvalues +#endif +#if ELPA_IMPL_SUFFIX == fc + !> \param a single complex matrix a: defines the problem to solve + !> \param b single complex matrix b: defines the problem to solve + !> \param ev single real: on output stores the eigenvalues + !> \param q single complex matrix q: on output stores the eigenvalues +#endif + + !> \param is_already_decomposed logical, input: is it repeated call with the same b (decomposed in the fist call)? + !> \result error integer, optional : error code, which can be queried with elpa_strerr + abstract interface + subroutine elpa_generalized_eigenvectors_& + &ELPA_IMPL_SUFFIX& + &_i(self, a, b, ev, q, is_already_decomposed, error) + use iso_c_binding + use elpa_constants + import elpa_t + implicit none + class(elpa_t) :: self +#ifdef USE_ASSUMED_SIZE + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows, *), b(self%local_nrows, *), q(self%local_nrows, *) +#else + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows, self%local_ncols), b(self%local_nrows, self%local_ncols), & + q(self%local_nrows, self%local_ncols) +#endif + real(kind=C_REAL_DATATYPE) :: ev(self%na) + + logical :: is_already_decomposed + integer, optional :: error + end subroutine + end interface + + + + !> \brief abstract definition of interface to compute C : = A**T * B + !> where A is a square matrix (self%a,self%na) which is optionally upper or lower triangular + !> B is a (self%na,ncb) matrix + !> C is a (self%na,ncb) matrix where optionally only the upper or lower + !> triangle may be computed + !> + !> the MPI commicators are already known to the type. Thus the class method "setup" must be called + !> BEFORE this method is used + !> \details + !> + !> \param self class(elpa_t), the ELPA object + !> \param uplo_a 'U' if A is upper triangular + !> 'L' if A is lower triangular + !> anything else if A is a full matrix + !> Please note: This pertains to the original A (as set in the calling program) + !> whereas the transpose of A is used for calculations + !> If uplo_a is 'U' or 'L', the other triangle is not used at all, + !> i.e. it may contain arbitrary numbers + !> \param uplo_c 'U' if only the upper diagonal part of C is needed + !> 'L' if only the upper diagonal part of C is needed + !> anything else if the full matrix C is needed + !> Please note: Even when uplo_c is 'U' or 'L', the other triangle may be + !> written to a certain extent, i.e. one shouldn't rely on the content there! + !> \param ncb Number of columns of global matrices B and C + !> \param a matrix a + !> \param self%local_nrows number of rows of local (sub) matrix a, set with method set("local_nrows,value") + !> \param self%local_ncols number of columns of local (sub) matrix a, set with method set("local_ncols,value") + !> \param b matrix b + !> \param nrows_b number of rows of local (sub) matrix b + !> \param ncols_b number of columns of local (sub) matrix b + !> \param nblk blocksize of cyclic distribution, must be the same in both directions! + !> \param c matrix c + !> \param nrows_c number of rows of local (sub) matrix c + !> \param ncols_c number of columns of local (sub) matrix c + !> \param error optional argument, error code which can be queried with elpa_strerr + abstract interface + subroutine elpa_hermitian_multiply_& + &ELPA_IMPL_SUFFIX& + &_i (self,uplo_a, uplo_c, ncb, a, b, nrows_b, ncols_b, & + c, nrows_c, ncols_c, error) + use iso_c_binding + import elpa_t + implicit none + class(elpa_t) :: self + character*1 :: uplo_a, uplo_c + integer(kind=c_int), intent(in) :: nrows_b, ncols_b, nrows_c, ncols_c, ncb +#ifdef USE_ASSUMED_SIZE + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows,*), b(nrows_b,*), c(nrows_c,*) +#else + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows,self%local_ncols), b(nrows_b,ncols_b), c(nrows_c,ncols_c) +#endif + +#ifdef USE_FORTRAN2008 + integer, optional :: error +#else + integer :: error +#endif + end subroutine + end interface + + + !> \brief abstract definition of interface to do a cholesky decomposition of a matrix + !> + !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution + !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE + !> with the class method "setup" + !> + !> Parameters + !> \param self class(elpa_t), the ELPA object +#if ELPA_IMPL_SUFFIX == d + !> \param a double real matrix: the matrix to be decomposed +#endif +#if ELPA_IMPL_SUFFIX == f + !> \param a single real matrix: the matrix to be decomposed +#endif +#if ELPA_IMPL_SUFFIX == dc + !> \param a double complex matrix: the matrix to be decomposed +#endif +#if ELPA_IMPL_SUFFIX == fc + !> \param a single complex matrix: the matrix to be decomposed +#endif + !> \param error integer, optional : error code, which can be queried with elpa_strerr + abstract interface + subroutine elpa_cholesky_& + &ELPA_IMPL_SUFFIX& + &_i (self, a, error) + use iso_c_binding + import elpa_t + implicit none + class(elpa_t) :: self +#ifdef USE_ASSUMED_SIZE + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows,*) +#else + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows,self%local_ncols) +#endif + +#ifdef USE_FORTRAN2008 + integer, optional :: error +#else + integer :: error +#endif + end subroutine + end interface + + + + !> \brief abstract definition of interface to invert a triangular matrix + !> + !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution + !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE + !> with the class method "setup" + !> + !> Parameters + !> \param self class(elpa_t), the ELPA object +#if ELPA_IMPL_SUFFIX == d + !> \param a double real matrix: the matrix to be inverted +#endif +#if ELPA_IMPL_SUFFIX == f + !> \param a single real matrix: the matrix to be inverted +#endif +#if ELPA_IMPL_SUFFIX == dc + !> \param a double complex matrix: the matrix to be inverted +#endif +#if ELPA_IMPL_SUFFIX == fc + !> \param a single complex matrix: the matrix to be inverted +#endif + + !> \param error integer, optional : error code, which can be queried with elpa_strerr + abstract interface + subroutine elpa_invert_trm_& + &ELPA_IMPL_SUFFIX& + &_i (self, a, error) + use iso_c_binding + import elpa_t + implicit none + class(elpa_t) :: self +#ifdef USE_ASSUMED_SIZE + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows,*) +#else + MATH_DATATYPE(kind=C_DATATYPE_KIND) :: a(self%local_nrows,self%local_ncols) +#endif + +#ifdef USE_FORTRAN2008 + integer, optional :: error +#else + integer :: error +#endif + end subroutine + end interface + + + + !> \brief abstract definition of interface to solve the eigenvalue problem for a valued tridiangular matrix + !> + !> The dimensions of the matrix a (locally ditributed and global), the block-cylic-distribution + !> block size, and the MPI communicators are already known to the object and MUST be set BEFORE + !> with the class method "setup" + !> + !> Parameters + !> \param self class(elpa_t), the ELPA object +#if ELPA_IMPL_SUFFIX == d + !> \param d double real 1d array: the diagonal elements of a matrix defined in setup, on output the eigenvalues + !> in ascending order + !> \param e double real 1d array: the subdiagonal elements of a matrix defined in setup + !> \param q double real matrix: on output contains the eigenvectors +#endif +#if ELPA_IMPL_SUFFIX == f + !> \param d single real 1d array: the diagonal elements of a matrix defined in setup, on output the eigenvalues + !> in ascending order + !> \param e single real 1d array: the subdiagonal elements of a matrix defined in setup + !> \param q single real matrix: on output contains the eigenvectors +#endif + !> \param error integer, optional : error code, which can be queried with elpa_strerr + abstract interface + subroutine elpa_solve_tridiagonal_& + &ELPA_IMPL_SUFFIX& + &_i (self, d, e, q, error) + use iso_c_binding + import elpa_t + implicit none + class(elpa_t) :: self + real(kind=C_REAL_DATATYPE) :: d(self%na), e(self%na) +#ifdef USE_ASSUMED_SIZE + real(kind=C_REAL_DATATYPE) :: q(self%local_nrows,*) +#else + real(kind=C_REAL_DATATYPE) :: q(self%local_nrows,self%local_ncols) +#endif + +#ifdef USE_FORTRAN2008 + integer, optional :: error +#else + integer :: error +#endif + end subroutine + end interface + diff --git a/src/general/precision_macros.h b/src/general/precision_macros.h index f05061a657202d64f88495fe6415e9b7eb42bb33..72c9a643b8e48dd1198c2b5efd9078fe8bf2dbd9 100644 --- a/src/general/precision_macros.h +++ b/src/general/precision_macros.h @@ -9,6 +9,7 @@ #undef SPECIAL_COMPLEX_DATATYPE #undef PRECISION_STR #undef REAL_DATATYPE +#undef C_REAL_DATATYPE #undef PRECISION_TRTRI #undef PRECISION_POTRF @@ -63,6 +64,7 @@ #define PRECISION_SUFFIX "_double" #define ELPA_IMPL_SUFFIX d #define REAL_DATATYPE rk8 +#define C_REAL_DATATYPE c_double #define BLAS_CHAR D #define BLAS_CHAR_AND_SY_OR_HE DSY #define SPECIAL_COMPLEX_DATATYPE ck8 @@ -116,6 +118,7 @@ #define PRECISION_SUFFIX "_single" #define ELPA_IMPL_SUFFIX f #define REAL_DATATYPE rk4 +#define C_REAL_DATATYPE c_float #define BLAS_CHAR S #define BLAS_CHAR_AND_SY_OR_HE SSY #define SPECIAL_COMPLEX_DATATYPE ck4 @@ -174,6 +177,7 @@ #undef COMPLEX_DATATYPE /* in the complex case also sometime real valued variables are needed */ #undef REAL_DATATYPE +#undef C_REAL_DATATYPE #undef PRECISION_TRTRI #undef PRECISION_POTRF @@ -238,6 +242,7 @@ #define BLAS_CHAR Z #define BLAS_CHAR_AND_SY_OR_HE ZHE #define REAL_DATATYPE RK8 +#define C_REAL_DATATYPE c_double #define PRECISION_TRTRI ZTRTRI #define PRECISION_POTRF ZPOTRF @@ -295,6 +300,7 @@ #define BLAS_CHAR C #define BLAS_CHAR_AND_SY_OR_HE CHE #define REAL_DATATYPE RK4 +#define C_REAL_DATATYPE c_float #define PRECISION_TRTRI CTRTRI #define PRECISION_POTRF CPOTRF