anchor.diff:      DEALLOCATE(recv_index)
       DEALLOCATE(recv_values)
       DEALLOCATE(dot_values)
  -    DEALLOCATE(diag_correction_index)
  -    DEALLOCATE(diag_correction_values)
       DEALLOCATE(a_buf)
  +    DEALLOCATE(row_root_lookup)
  +    DEALLOCATE(col_root_lookup)
       CALL DestructSparseMatrix(sparse_a)
     END SUBROUTINE PivotedCholeskyDecomposition
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
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:           recv_index(:recv_num_values) = index_l(:recv_num_values,local_pi_j)
            recv_values(:recv_num_values) = values_l(:recv_num_values, local_pi_j)
         ELSE
            inverse_factor = 0
         END IF
         col_root = col_root_lookup(pi_j)
  
         !! Compute A Buffer AMat[:,pi_j]
         local_pi_j = pi_j - AMat%start_row + 1
         IF (pi_j .GE. AMat%start_row .AND. pi_j .LT. AMat%end_row) THEN
            DO II = 1, num_local_pivots
               local_pi_i = local_pivots(II)
               a_buf(II) = dense_a(local_pi_j,local_pi_i)
            END DO
         END IF
         row_root = row_root_lookup(pi_j)
  
         !! Broadcast column JJ, and Inverse Factor
         CALL BroadcastVector(recv_num_values, recv_index, recv_values, &
              & col_root, row_comm)
         CALL MPI_Bcast(inverse_factor, 1, MPINTREAL, col_root, row_comm, &
              & grid_error)
  
         !! Broadcast A values
         CALL BroadcastARow(num_local_pivots,a_buf,row_root,column_comm)
  
         !! Compute Dot Products
         CALL DotAllPivoted(recv_num_values, recv_index, recv_values, &
              & values_per_column_l, index_l, values_l, local_pivots, &
              & num_local_pivots, dot_values, column_comm)
  
         !! Loop over other columns
         DO II = 1, num_local_pivots
            !! Insert Into L
            Aval = a_buf(II)
            insert_value = inverse_factor * (Aval - dot_values(II))
            local_pi_i = local_pivots(II)
            IF (JJ .GE. AMat%start_row .AND. JJ .LT. AMat%end_row) THEN
               CALL AppendToVector(values_per_column_l(local_pi_i), &
                    & index_l(:,local_pi_i), values_l(:, local_pi_i), &
  