dbcsr_fill_wm_from_matrix Subroutine

   SUBROUTINE dbcsr_fill_wm_from_matrix(wm, matrix, size_used, limits)
      !! Fills index and data of the work matrix from the
      !! previously-finalized one.
      !!
      !! limits
      !! The limits is a 4-tuple
      !! (lower_row, higher_row, lower_column, higher_column).

      TYPE(dbcsr_work_type), DIMENSION(:), INTENT(INOUT) :: wm
         !! the work matrix to fill
      TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix
         !! DBCSR matrix
      INTEGER, INTENT(IN)                                :: size_used
      INTEGER, DIMENSION(4), INTENT(IN), OPTIONAL        :: limits
         !! only fills blocks within this range

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

      INTEGER                                            :: blk, blk_p, col, handle, ithread, &
                                                            nthreads, nwms, nze, row, wblk_p, &
                                                            which_wm, wm_first, wm_last
      LOGICAL                                            :: careful, limit, mt, tr
      LOGICAL, SAVE                                      :: mutable
      TYPE(dbcsr_data_obj)                               :: data_block
      TYPE(dbcsr_iterator)                               :: iter

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

      CALL timeset(routineN, handle)
      nwms = SIZE(matrix%wms)
      mt = .FALSE.
!$    IF (nwms .GT. 1) mt = omp_get_num_threads() .GT. 1
      ithread = 0; nthreads = 1
!$    ithread = omp_get_thread_num()
!$    nthreads = omp_get_num_threads()
      limit = PRESENT(limits)
      careful = size_used + size_used/8 &
                .LT. dbcsr_data_get_size_referenced(matrix%data_area)
      CALL dbcsr_data_init(data_block)
      CALL dbcsr_data_new(data_block, dbcsr_data_get_type(matrix%data_area))
      IF (mt) THEN
         wm_first = ithread + 1
         wm_last = ithread + 1
      ELSE
         wm_first = 1
         wm_last = nwms
      END IF
      ! Prepares the work matrices to accept the main data.
!$OMP     MASTER
      mutable = .FALSE.
      DO which_wm = 1, nwms
         mutable = mutable .OR. dbcsr_wm_use_mutable(wm(which_wm))
      END DO
!$OMP     END MASTER
!$OMP     BARRIER
      DO which_wm = wm_first, wm_last
         IF (dbcsr_wm_use_mutable(wm(which_wm))) &
            DBCSR_ABORT("Adding main matrix into mutable not supported.")
         IF (mutable) THEN
            IF (.NOT. dbcsr_mutable_instantiated(wm(which_wm)%mutable)) THEN
               CALL dbcsr_mutable_new(wm(which_wm)%mutable, matrix%data_type)
            END IF
         ELSE
            ! We don't know how much data we'll get so we have to be generous.
            CALL dbcsr_data_ensure_size(wm(which_wm)%data_area, &
                                        size_used/nwms)
            CALL dbcsr_data_set_size_referenced(wm(which_wm)%data_area, 0)
         END IF
      END DO
      ! Now copy the data
      CALL dbcsr_iterator_start(iter, matrix, shared=mt, &
                                contiguous_pointers=.FALSE., read_only=.TRUE.)
      DO WHILE (dbcsr_iterator_blocks_left(iter))
         CALL dbcsr_iterator_next_block(iter, row, col, data_block, &
                                        transposed=tr, block_number=blk)
         IF (limit) THEN
            IF (.NOT. within_limits(row, col, limits)) CYCLE
         END IF
         blk_p = matrix%blk_p(blk)
         which_wm = ithread + 1
         wblk_p = SIGN(wm(which_wm)%datasize + 1, blk_p)
         nze = dbcsr_data_get_size(data_block)
         IF (mt .OR. limit .OR. careful .OR. mutable) THEN
            ! The data gets copied block by block so the block pointers
            ! are ordered accordingly.
            IF (.NOT. mutable) THEN
               CALL add_work_coordinate(wm(which_wm), row, col, wblk_p)
               CALL dbcsr_data_ensure_size(wm(which_wm)%data_area, &
                                           ABS(wblk_p) + nze - 1, factor=default_resize_factor)
               CALL dbcsr_data_set_size_referenced(wm(which_wm)%data_area, &
                                                   ABS(wblk_p) + nze - 1)
               CALL dbcsr_data_set(wm(which_wm)%data_area, &
                                   lb=ABS(wblk_p), &
                                   data_size=nze, &
                                   src=data_block, source_lb=1)
            END IF
         ELSE
            ! The data gets copied all at once so the block pointers
            ! should remain the same as they were.
            CALL add_work_coordinate(wm(which_wm), row, col, blk_p)
         END IF
         IF (.NOT. mutable) &
            wm(which_wm)%datasize = wm(which_wm)%datasize + nze
      END DO
      CALL dbcsr_iterator_stop(iter)
      CALL dbcsr_data_clear_pointer(data_block)
      CALL dbcsr_data_release(data_block)
      ! Copy all blocks at once
      IF (.NOT. mt .AND. .NOT. limit .AND. .NOT. careful .AND. .NOT. mutable) THEN
         DO which_wm = 1, nwms
            CALL dbcsr_data_ensure_size(wm(which_wm)%data_area, &
                                        dbcsr_data_get_size_referenced(matrix%data_area))
            CALL dbcsr_data_copyall(wm(which_wm)%data_area, matrix%data_area)
            wm(which_wm)%datasize = size_used
         END DO
      END IF
      CALL timestop(handle)
   END SUBROUTINE dbcsr_fill_wm_from_matrix