dbcsr_filter_anytype Subroutine

   SUBROUTINE dbcsr_filter_anytype(matrix, eps, method, &
                                   use_absolute, filter_diag)
      !! filter a dbcsr matrix

      TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix
         !! the matrix
      TYPE(dbcsr_scalar_type), INTENT(IN)                :: eps
         !! the threshold
      INTEGER, INTENT(IN), OPTIONAL                      :: method
         !! how the matrix is filtered
      LOGICAL, INTENT(in), OPTIONAL                      :: use_absolute, filter_diag
         !! NYI

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

      COMPLEX(KIND=real_4), DIMENSION(:), POINTER        :: data_c
      COMPLEX(KIND=real_8), DIMENSION(:), POINTER        :: data_z
      INTEGER                                            :: blk, blk_nze, col, col_size, handle, &
                                                            my_method, row, row_size, data_type
      LOGICAL                                            :: gt0, my_filter_diag, tr
      REAL(KIND=real_4)                                  :: nrm_s
      REAL(KIND=real_4), DIMENSION(:), POINTER           :: data_s
      REAL(KIND=real_8)                                  :: my_absolute, nrm_d
      REAL(KIND=real_8), DIMENSION(:), POINTER           :: data_d
      TYPE(dbcsr_iterator)                               :: iter

      REAL(KIND=real_8), EXTERNAL                        :: DZNRM2
#if defined (__ACCELERATE)
      REAL(KIND=real_8), EXTERNAL                        :: SCNRM2
#else
      REAL(KIND=real_4), EXTERNAL                        :: SCNRM2
#endif

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

      CALL timeset(routineN, handle)
      my_method = dbcsr_filter_frobenius
      IF (PRESENT(method)) my_method = method
      my_absolute = 1.0_dp
      IF (PRESENT(use_absolute)) my_absolute = dbcsr_maxabs(matrix)
      my_filter_diag = .TRUE.
      IF (PRESENT(filter_diag)) my_filter_diag = filter_diag

      SELECT CASE (eps%data_type)
      CASE (dbcsr_type_real_4)
         gt0 = eps%r_sp .GT. 0.0_real_4
      CASE (dbcsr_type_real_8)
         gt0 = eps%r_dp .GT. 0.0_real_8
      CASE (dbcsr_type_complex_4)
         gt0 = ABS(eps%c_sp) .GT. 0.0_real_4
      CASE (dbcsr_type_complex_8)
         gt0 = ABS(eps%c_dp) .GT. 0.0_real_8
      CASE default
         gt0 = .FALSE.
      END SELECT

      IF (gt0) THEN
         data_type = dbcsr_get_data_type(matrix)

!$OMP        PARALLEL DEFAULT(NONE) PRIVATE(iter,row,col,data_s,data_d,data_c,data_z,tr, &
!$OMP                                       blk,row_size,col_size,blk_nze,nrm_d,nrm_s) &
!$OMP                               SHARED(my_method,my_absolute,eps,matrix,data_type)

         CALL dbcsr_iterator_start(iter, matrix, contiguous_pointers=.TRUE.)
         DO WHILE (dbcsr_iterator_blocks_left(iter))
            SELECT CASE (data_type)

            CASE (dbcsr_type_real_4)
               CALL dbcsr_iterator_next_block(iter, row, col, data_s, tr, blk, &
                                              row_size, col_size)
               blk_nze = row_size*col_size
               IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks
               SELECT CASE (my_method)
               CASE (dbcsr_filter_frobenius)
                  !
                  ! Frobenius based
                  nrm_s = norm2(data_s)
                  IF (nrm_s .LT. my_absolute*eps%r_sp) &
                     CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk)
               CASE DEFAULT
                  DBCSR_ABORT("Only Frobenius based filtering")
               END SELECT

            CASE (dbcsr_type_real_8)
               CALL dbcsr_iterator_next_block(iter, row, col, data_d, tr, blk, &
                                              row_size, col_size)
               blk_nze = row_size*col_size
               IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks
               SELECT CASE (my_method)
               CASE (dbcsr_filter_frobenius)
                  !
                  ! Frobenius based
                  nrm_d = norm2(data_d)
                  IF (nrm_d .LT. my_absolute*eps%r_dp) &
                     CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk)
               CASE DEFAULT
                  DBCSR_ABORT("Only Frobenius based filtering")
               END SELECT

            CASE (dbcsr_type_complex_4)
               CALL dbcsr_iterator_next_block(iter, row, col, data_c, tr, blk, &
                                              row_size, col_size)
               blk_nze = row_size*col_size
               IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks
               SELECT CASE (my_method)
               CASE (dbcsr_filter_frobenius)
                  !
                  ! Frobenius based
                  nrm_d = SCNRM2(SIZE(data_c), data_c(1), 1)
                  IF (nrm_d .LT. my_absolute*eps%r_dp) &
                     CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk)
               CASE DEFAULT
                  DBCSR_ABORT("Only Frobenius based filtering")
               END SELECT

            CASE (dbcsr_type_complex_8)
               CALL dbcsr_iterator_next_block(iter, row, col, data_z, tr, blk, &
                                              row_size, col_size)
               blk_nze = row_size*col_size
               IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks
               SELECT CASE (my_method)
               CASE (dbcsr_filter_frobenius)
                  !
                  ! Frobenius based
                  nrm_d = DZNRM2(SIZE(data_z), data_z(1), 1)
                  IF (nrm_d .LT. my_absolute*eps%r_dp) &
                     CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk)
               CASE DEFAULT
                  DBCSR_ABORT("Only Frobenius based filtering")
               END SELECT

            CASE DEFAULT
               DBCSR_ABORT("Wrong data type")
            END SELECT
         END DO
         CALL dbcsr_iterator_stop(iter)
         CALL dbcsr_finalize(matrix, reshuffle=.TRUE.)
!$OMP        END PARALLEL

         CALL dbcsr_index_compact(matrix)

      END IF
      CALL timestop(handle)
   END SUBROUTINE dbcsr_filter_anytype