# dbcsr_function_of_elements Subroutine

## public subroutine dbcsr_function_of_elements(matrix_a, func, a0, a1, a2)

Computes various functions (defined by func) of matrix elements

Note

sign(A,B) returns the value of A with the sign of B dbcsr_func_inverse: 1/(a1x+a0) fails if the inversion produces infinite numbers dbcsr_func_inverse_special: 1/(x+sign(a0,x)) safe inverse: if a0>0 then the denominator is never zero dbcsr_func_tanh: tanh(a1x+a0) dbcsr_func_dtanh: d(tanh(a1x+a0)) / dx dbcsr_func_ddtanh: d2(tanh(a1x+a0)) / dx2 dbcsr_func_artanh: artanh(a1x+a0)=ln[(1+(a1x+a0))/(1-(a1x+a0))]/2 fails if |a1x+a0| >= 1 dbcsr_func_sread_from_zero: if |x|<|a0| then x=sign(a0,x) dbcsr_func_truncate: if |x|>|a0| then x=sign(a0,x) dbcsr_func_sin: sin(a1x+a0) dbcsr_func_cos: cos(a1x+a0) dbcsr_func_dsin: d(sin(a1x+a0)) / dx = a1cos(a1x+a0) dbcsr_func_ddsin: d2(sin(a1x+a0)) / dx2 = -a1a1sin(a1x+a0) dbcsr_func_asin: asin(a1x+a0) fails if |a1*x+a0| > 1

### Arguments

Type IntentOptional Attributes Name
type(dbcsr_type), intent(inout) :: matrix_a

DBCSR matrix

integer, intent(in) :: func
real(kind=dp), intent(in), optional :: a0
real(kind=dp), intent(in), optional :: a1
real(kind=dp), intent(in), optional :: a2

## Source Code

   SUBROUTINE dbcsr_function_of_elements(matrix_a, func, a0, a1, a2)
!! Computes various functions (defined by func) of matrix elements
!! @note  sign(A,B) returns the value of A with the sign of B
!! dbcsr_func_inverse:   1/(a1*x+a0)
!! fails if the inversion produces infinite numbers
!! dbcsr_func_inverse_special: 1/(x+sign(a0,x))
!! safe inverse: if a0>0 then the denominator is never zero
!! dbcsr_func_tanh:    tanh(a1*x+a0)
!! dbcsr_func_dtanh:   d(tanh(a1*x+a0)) / dx
!! dbcsr_func_ddtanh:  d2(tanh(a1*x+a0)) / dx2
!! dbcsr_func_artanh:  artanh(a1*x+a0)=ln[(1+(a1*x+a0))/(1-(a1*x+a0))]/2
!! fails if |a1*x+a0| >= 1
!! dbcsr_func_sread_from_zero:  if |x|<|a0| then x=sign(a0,x)
!! dbcsr_func_truncate:  if |x|>|a0| then x=sign(a0,x)
!! dbcsr_func_sin:     sin(a1*x+a0)
!! dbcsr_func_cos:     cos(a1*x+a0)
!! dbcsr_func_dsin:    d(sin(a1*x+a0)) / dx = a1*cos(a1*x+a0)
!! dbcsr_func_ddsin:   d2(sin(a1*x+a0)) / dx2 = -a1*a1*sin(a1*x+a0)
!! dbcsr_func_asin:    asin(a1*x+a0)
!! fails if |a1*x+a0| > 1

TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_a
!! DBCSR matrix
INTEGER, INTENT(IN)                                :: func
REAL(kind=dp), INTENT(IN), OPTIONAL                :: a0, a1, a2

CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_function_of_elements'

INTEGER                                            :: blk, col, col_size, data_type, handle, &
ielem, nze, row, row_size
LOGICAL                                            :: tr_a
REAL(kind=dp)                                      :: p0, p1, p2
TYPE(dbcsr_data_obj)                               :: a_data
TYPE(dbcsr_iterator)                               :: iter

!   ---------------------------------------------------------------------------

CALL timeset(routineN, handle)

IF (PRESENT(a0)) THEN
p0 = a0
ELSE
p0 = 0.0_dp
END IF
IF (PRESENT(a1)) THEN
p1 = a1
ELSE
p1 = 1.0_dp
END IF
IF (PRESENT(a2)) THEN
p2 = a2
ELSE
p2 = 0.0_dp
END IF

data_type = dbcsr_get_data_type(matrix_a)
CALL dbcsr_data_init(a_data)
CALL dbcsr_data_new(a_data, data_type)
CALL dbcsr_iterator_start(iter, matrix_a)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, row, col, a_data, tr_a, blk, &
row_size=row_size, col_size=col_size)
nze = row_size*col_size
SELECT CASE (data_type)
!CASE (dbcsr_type_real_4)
!   a_data%d%r_sp(1:nze) = 1.0_real_4/a_data%d%r_sp(1:nze)
!   IF(MAXVAL(ABS(a_data%d%r_sp)).GE.HUGE(0.0_real_4))&
!      DBCSR_ABORT("Division by zero")
CASE (dbcsr_type_real_8)
SELECT CASE (func)
! if |x|<|a0| then x=|a0|*sign(x)
DO ielem = 1, nze
IF (ABS(a_data%d%r_dp(ielem)) .LT. ABS(p0)) THEN
a_data%d%r_dp(ielem) = SIGN(p0, a_data%d%r_dp(ielem))
END IF
END DO
CASE (dbcsr_func_truncate)
! if |x|>|a0| then x=|a0|*sign(x)
DO ielem = 1, nze
IF (ABS(a_data%d%r_dp(ielem)) .GT. ABS(p0)) THEN
a_data%d%r_dp(ielem) = SIGN(p0, a_data%d%r_dp(ielem))
END IF
END DO
CASE (dbcsr_func_inverse_special)
!IF (MINVAL(ABS(a_data%d%r_dp)).le.ABS(p2)) THEN
!   ! there is at least one near-zero element,
!   ! invert element-by-element
!   DO ielem=1,nze
!     IF (a_data%d%r_dp(ielem).le.ABS(p2)) THEN
!        a_data%d%r_dp(ielem) = 0.0_real_8
!     ELSE
!        a_data%d%r_dp(ielem) = &
!           1.0_real_8/(p1*a_data%d%r_dp(ielem)+p0)
!     ENDIF
!   ENDDO
!ELSE
!   a_data%d%r_dp(1:nze) = 1.0_real_8/(p1*a_data%d%r_dp(1:nze)+p0)
!ENDIF
a_data%d%r_dp(1:nze) = 1.0_real_8/(a_data%d%r_dp(1:nze) + SIGN(p0, a_data%d%r_dp(1:nze)))
CASE (dbcsr_func_inverse)
a_data%d%r_dp(1:nze) = 1.0_real_8/(p1*a_data%d%r_dp(1:nze) + p0)
IF (MAXVAL(ABS(a_data%d%r_dp)) .GE. HUGE(0.0_real_8)) &
DBCSR_ABORT("Division by zero")
CASE (dbcsr_func_tanh)
a_data%d%r_dp(1:nze) = TANH(p1*a_data%d%r_dp(1:nze) + p0)
CASE (dbcsr_func_dtanh)
a_data%d%r_dp(1:nze) = TANH(p1*a_data%d%r_dp(1:nze) + p0)
a_data%d%r_dp(1:nze) = a_data%d%r_dp(1:nze)**2
a_data%d%r_dp(1:nze) = p1*(1.0_real_8 - a_data%d%r_dp(1:nze))
CASE (dbcsr_func_ddtanh)
a_data%d%r_dp(1:nze) = TANH(p1*a_data%d%r_dp(1:nze) + p0)
a_data%d%r_dp(1:nze) = a_data%d%r_dp(1:nze)**3 - a_data%d%r_dp(1:nze)
a_data%d%r_dp(1:nze) = 2.0_real_8*(p1**2)*a_data%d%r_dp(1:nze)
CASE (dbcsr_func_artanh)
a_data%d%r_dp(1:nze) = p1*a_data%d%r_dp(1:nze) + p0
IF (MAXVAL(ABS(a_data%d%r_dp)) .GE. 1.0_real_8) &
DBCSR_ABORT("ARTANH is undefined for |x|>=1")
a_data%d%r_dp(1:nze) = (1.0_real_8 + a_data%d%r_dp(1:nze)) &
/(1.0_real_8 - a_data%d%r_dp(1:nze))
a_data%d%r_dp(1:nze) = 0.5_real_8*LOG(a_data%d%r_dp(1:nze))
CASE (dbcsr_func_sin)
a_data%d%r_dp(1:nze) = SIN(p1*a_data%d%r_dp(1:nze) + p0)
CASE (dbcsr_func_cos)
a_data%d%r_dp(1:nze) = COS(p1*a_data%d%r_dp(1:nze) + p0)
CASE (dbcsr_func_dsin)
a_data%d%r_dp(1:nze) = p1*COS(p1*a_data%d%r_dp(1:nze) + p0)
CASE (dbcsr_func_ddsin)
a_data%d%r_dp(1:nze) = -p1*p1*SIN(p1*a_data%d%r_dp(1:nze) + p0)
CASE (dbcsr_func_asin)
a_data%d%r_dp(1:nze) = p1*a_data%d%r_dp(1:nze) + p0
IF (MAXVAL(ABS(a_data%d%r_dp)) .GT. 1.0_real_8) &
DBCSR_ABORT("ASIN is undefined for |x|>1")
a_data%d%r_dp(1:nze) = ASIN(a_data%d%r_dp(1:nze))
CASE DEFAULT
DBCSR_ABORT("Unknown function of matrix elements")
END SELECT
!CASE (dbcsr_type_complex_4)
!CASE (dbcsr_type_complex_8)
CASE DEFAULT
DBCSR_ABORT("Operation is implemented only for dp real values")
END SELECT
END DO
CALL dbcsr_iterator_stop(iter)
CALL dbcsr_data_clear_pointer(a_data)
CALL dbcsr_data_release(a_data)
CALL timestop(handle)

END SUBROUTINE dbcsr_function_of_elements