anchor.diff:      CALL MPI_Allgatherv(MPI_IN_PLACE, diags_per_proc(my_row+1), MPINTREAL, &
            & diag, diags_per_proc, diag_displ, MPINTREAL, column_comm, grid_error)
   
  +    !! Root Lookups
  +    ALLOCATE(row_root_lookup(AMat%logical_matrix_dimension))
  +    row_root_lookup = 0
  +    ALLOCATE(col_root_lookup(AMat%logical_matrix_dimension))
  +    col_root_lookup = 0
  +    rows_per_proc(my_row+1) = AMat%local_rows
  +    cols_per_proc(my_column+1) = AMat%local_columns
  +    CALL MPI_Allgather(MPI_IN_PLACE, 1, MPI_INTEGER, rows_per_proc, 1, &
  +         & MPI_INTEGER, column_comm, grid_error)
  +    CALL MPI_Allgather(MPI_IN_PLACE, 1, MPI_INTEGER, cols_per_proc, 1, &
  +         & MPI_INTEGER, row_comm, grid_error)
  +    displ_rows_per_proc(1) = 0
  +    DO II = 2, num_process_rows
  +       displ_rows_per_proc(II) = displ_rows_per_proc(II-1) + rows_per_proc(II-1)
  +    END DO
  +    displ_cols_per_proc(1) = 0
  +    DO II = 2, num_process_columns
  +       displ_cols_per_proc(II) = displ_cols_per_proc(II-1) + cols_per_proc(II-1)
  +    END DO
  +    row_root_lookup(AMat%start_row:AMat%end_row - 1) = column_rank
  +    col_root_lookup(AMat%start_column:AMat%end_column - 1) = row_rank
  +    CALL MPI_Allgatherv(MPI_IN_PLACE, AMat%local_rows, MPI_INTEGER, &
  +         & row_root_lookup, rows_per_proc, displ_rows_per_proc, MPI_INTEGER, &
  +         & column_comm, grid_error)
  +    CALL MPI_Allgatherv(MPI_IN_PLACE, AMat%local_columns, MPI_INTEGER, &
  +         & col_root_lookup, cols_per_proc, displ_cols_per_proc, MPI_INTEGER, &
  +         & row_comm, grid_error)
  +     WRITE(*,*) "::", row_root_lookup
  +     WRITE(*,*) ",,", num_process_rows, num_process_columns
  +
       !! Allocate space for L
       ALLOCATE(values_per_column_l(sparse_a%columns))
       ALLOCATE(index_l(sparse_a%rows, sparse_a%columns))
  
positive.code:     ALLOCATE(recv_index(sparse_a%rows))
      ALLOCATE(recv_values(sparse_a%rows))
      ALLOCATE(dot_values(sparse_a%columns))
      ALLOCATE(diag_correction_index(sparse_a%columns))
      ALLOCATE(diag_correction_values(sparse_a%columns))
      ALLOCATE(local_pivots(sparse_a%columns))
      !! Main Loop
      DO JJ = 1, rank_in
         !! Pick a pivot vector
         local_JJ  = JJ - AMat%start_row + 1
         CALL GetPivot(AMat, JJ, pivot_vector, diag, pi_j, pivot_value)
         !! l[pi[j],j] = sqrt(d[pi[j]])
         IF (pi_j .GE. AMat%start_column .AND. pi_j .LT. AMat%end_column) THEN
            root = row_rank
            local_pi_j = pi_j - AMat%start_column + 1
            insert_value = SQRT(diag(local_pi_j))
            inverse_factor = 1.0/insert_value
  
negative.code:     !! Which pivots to treat locally
      INTEGER, DIMENSION(:), ALLOCATABLE :: local_pivots
      INTEGER :: num_local_pivots
      !! Root Lookups
      INTEGER, DIMENSION(:), ALLOCATABLE :: row_root_lookup
      INTEGER, DIMENSION(:), ALLOCATABLE :: col_root_lookup
      INTEGER, DIMENSION(num_process_rows) :: rows_per_proc
      INTEGER, DIMENSION(num_process_columns) :: cols_per_proc
      INTEGER, DIMENSION(num_process_rows) :: displ_rows_per_proc
      INTEGER, DIMENSION(num_process_columns) :: displ_cols_per_proc
      !! Temporary Variables
      REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: a_buf
      INTEGER :: II, JJ, local_JJ, local_II
      INTEGER :: local_pi_i, local_pi_j
      REAL(NTREAL) :: Aval, insert_value, inverse_factor
      REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: dot_values
      INTEGER :: row_root, col_root
      INTEGER :: fill_counter
      INTEGER, DIMENSION(num_process_rows) :: diags_per_proc
      INTEGER, DIMENSION(num_process_rows) :: diag_displ
  