dbcsr_function_of_elements Subroutine

   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