Computes various functions (defined by func) of matrix elements
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
| Type | Intent | Optional | 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 |
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)
CASE (dbcsr_func_spread_from_zero)
! 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