В первую очередь параллель это:
for (int k = 0; k < i; k++)
temp = temp + U[k, i] * U[k, j];
Просто разбей диапазон значений j на M ядер.
Проще всего так. Если у нас есть N операций, и мы хотим разбить на M ядер то это (N+M-1)/M с округлением вниз операций на ядро.
Update: Не всё так просто как мне казалось, в общем, это не ответ.
Можно другие попробовать распараллелить, но это это самая ресурсоёмкая часть.
Вообще хорошо
гуглится описание распараллеливания этого алгоритма.