make_index_triangular Subroutine

   SUBROUTINE make_index_triangular(new_row_p, new_col_i, new_blk_p, &
                                    old_row_p, old_col_i, old_blk_p, matrix)
      !! Makes a CP2K triangular index given the index arrays
      !!
      !! Description of canonical ordering
      !! A non-(anti)symmetric matrix is left as is. Otherwise, the row and column
      !! are stored in the position prescribed by the distribution.
      !! @note
      !! This routine uses hard-coded logic as to what constitutes a
      !! canonical ordering

      INTEGER, DIMENSION(:), INTENT(OUT)                 :: new_row_p, new_col_i, new_blk_p
      INTEGER, DIMENSION(:), INTENT(IN)                  :: old_row_p, old_col_i, old_blk_p
      TYPE(dbcsr_type), INTENT(IN)                       :: matrix

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

      INTEGER                                            :: blk, col, nblks, row, stored_col, &
                                                            stored_row
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: row_i
      LOGICAL                                            :: tr

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

      nblks = SIZE(old_blk_p)
      ALLOCATE (row_i(nblks))
      IF (debug_mod) THEN
         WRITE (*, *) "old row_p", old_row_p
         WRITE (*, *) "old col_i", old_col_i
         WRITE (*, *) "old blk_p", old_blk_p
      END IF
      DO row = 1, SIZE(old_row_p) - 1
         DO blk = old_row_p(row) + 1, old_row_p(row + 1)
            col = old_col_i(blk)
            stored_row = row
            stored_col = col
            tr = .FALSE.
            CALL get_stored_canonical(matrix, stored_row, stored_col, tr)
            IF (stored_row .GT. stored_col) THEN
               CALL swap(stored_row, stored_col)
               tr = .NOT. tr
            END IF
            IF (debug_mod) &
               WRITE (*, '(A,2(1X,I5),A,2(1X,I5),";",I7,1X,L1)') &
               routineN//" X->", row, col, "->", &
               stored_row, stored_col, blk, tr
            row_i(blk) = stored_row
            new_col_i(blk) = stored_col
            IF (.NOT. tr) THEN
               new_blk_p(blk) = old_blk_p(blk)
            ELSE
               new_blk_p(blk) = -old_blk_p(blk)
            END IF
         END DO
      END DO
      CALL dbcsr_sort_indices(nblks, row_i, new_col_i, blk_p=new_blk_p)
      ! Re-create the index
      CALL dbcsr_make_dbcsr_index(new_row_p, row_i, SIZE(new_row_p) - 1, nblks)
      IF (debug_mod) THEN
         WRITE (*, *) "new row_p", new_row_p
         WRITE (*, *) "new row_i", row_i
         WRITE (*, *) "new col_i", new_col_i
         WRITE (*, *) "new blk_p", new_blk_p
      END IF
   END SUBROUTINE make_index_triangular