anchor.diff:         & RedistributeTripletLists, ShiftTripletList, SortTripletList
     USE TripletModule, ONLY : Triplet_t
     USE MPI
  +#ifdef _OPENMP
  +  USE omp_lib, ONLY : omp_get_num_threads, omp_get_thread_num
  +#endif
     IMPLICIT NONE
     PRIVATE
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
positive.code:        CALL ComputePartners(jdata, jdata%final_swap, swap1, swap0)
         jdata%phase_array(stage_counter) = 4
      END IF
  
      !! Cleanup
      DEALLOCATE(swap0)
  
negative.code:     !! Temporary
      INTEGER :: counter
      INTEGER :: ind
      ! TYPE(DenseMatrix_t) :: TempMat
      ! TYPE(DenseMatrix_t) :: AMat
      INTEGER, DIMENSION(2) :: row_sizes, col_sizes
      TYPE(DenseMatrix_t), DIMENSION(jdata%num_threads) :: TempArray
      TYPE(DenseMatrix_t), DIMENSION(jdata%num_threads) :: AArray
      INTEGER :: thread_id
  
      CALL BlockingDenseMatrixGather(TargetV, receive_list, jdata%communicator)
  
      !$OMP PARALLEL PRIVATE(thread_id, ind, row_sizes, col_sizes)
      !$OMP DO
      DO counter = 1, jdata%num_processes
  #ifdef _OPENMP
         thread_id = omp_get_thread_num() + 1
  #else
         thread_id = 1
  #endif
         ind = (counter-1)*2 + 1
         row_sizes(1) = ABlocks(1,ind)%rows
         row_sizes(2) = ABlocks(2,ind+1)%rows
         col_sizes(1) = ABlocks(1,ind)%columns
         col_sizes(2) = ABlocks(2,ind+1)%columns
  
         CALL ComposeDenseMatrix(ABlocks(1:2,ind:ind+1), 2, 2, AArray(thread_id))
         CALL MultiplyDense(AArray(thread_id), receive_list(counter), &
              & TempArray(thread_id))
         CALL SplitDenseMatrix(TempArray(thread_id), 2, 2, &
              & ABlocks(1:2,ind:ind+1), block_size_row_in=row_sizes, &
              & block_size_column_in=col_sizes)
      END DO
      !$OMP END DO
      !$OMP END PARALLEL
      ! DO counter = 1, jdata%num_processes
      !    ind = (counter-1)*2 + 1
      !    row_sizes(1) = ABlocks(1,ind)%rows
      !    row_sizes(2) = ABlocks(2,ind+1)%rows
      !    col_sizes(1) = ABlocks(1,ind)%columns
      !    col_sizes(2) = ABlocks(2,ind+1)%columns
      !
      !    CALL ComposeDenseMatrix(ABlocks(1:2,ind:ind+1), 2, 2, AMat)
      !    CALL MultiplyDense(AMat, receive_list(counter), TempMat)
      !    CALL SplitDenseMatrix(TempMat, 2, 2, ABlocks(1:2,ind:ind+1), &
      !         & block_size_row_in=row_sizes, block_size_column_in=col_sizes)
      ! END DO
  
      DO counter = 1, jdata%num_threads
         CALL DestructDenseMatrix(TempArray(counter))
         CALL DestructDenseMatrix(AArray(counter))
      END DO
  
      DO counter = 1, jdata%num_processes
         CALL DestructDenseMatrix(receive_list(counter))
      END DO
  
      ! CALL DestructDenseMatrix(TempMat)
  
    END SUBROUTINE ApplyToColumns
  END MODULE JacobiEigenSolverModule
  