anchor.diff:      dense_a = 0
       CALL ConstructDenseFromSparse(sparse_a, dense_a)
   
  -    !! Construct the vector holding the accumulated diagonal values
  +    !! Extract the diagonal
       ALLOCATE(diag(sparse_a%columns))
  -    diag = 0
  -    fill_counter = 0
  -    DO JJ = AMat%start_column, AMat%end_column - 1
  -       IF (JJ .GE. AMat%start_row .AND. JJ .LT. AMat%end_row) THEN
  -          local_JJ = JJ - AMat%start_column + 1
  -          local_II = JJ - AMat%start_row + 1
  -          diag(local_JJ) = dense_a(local_II, local_JJ)
  -          fill_counter = fill_counter + 1
  -       END IF
  -    END DO
  -    diags_per_proc(my_row+1) = fill_counter
  -    !! Duplicate the diagonal entries along the process column (across rows)
  -    CALL MPI_Allgather(MPI_IN_PLACE, 1, MPI_INTEGER, diags_per_proc, 1, &
  -         & MPI_INTEGER, column_comm, grid_error)
  -    diag_displ(1) = 0
  -    DO II = 2, num_process_rows
  -       diag_displ(II) = diag_displ(II-1) + diags_per_proc(II-1)
  -    END DO
  -    CALL MPI_Allgatherv(MPI_IN_PLACE, diags_per_proc(my_row+1), MPINTREAL, &
  -         & diag, diags_per_proc, diag_displ, MPINTREAL, column_comm, grid_error)
  +    CALL ConstructDiag(AMat, dense_a, diag)
   
       !! 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
  +    CALL ConstructRankLookup(AMat, row_root_lookup, col_root_lookup)
   
       !! Allocate space for L
       ALLOCATE(values_per_column_l(sparse_a%columns))
  
positive.code:     DEALLOCATE(recv_index)
      DEALLOCATE(recv_values)
      DEALLOCATE(dot_values)
      CALL DestructSparseMatrix(sparse_a)
    END SUBROUTINE CholeskyDecomposition
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
negative.code:     INTEGER :: II, JJ, local_JJ, local_II, local_row
      REAL(NTREAL) :: Aval, insert_value, inverse_factor
      REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: dot_values
      INTEGER, DIMENSION(:), ALLOCATABLE :: col_root_lookup
      INTEGER :: col_root
  
      !! Optional Parameters
      IF (PRESENT(solver_parameters_in)) THEN
  