О Терафлопсах

во-первых, конечные элементы и конечные разности дают разряженные матрицы

зависит от задачи, но да, в большинстве случаев МКЭ это разряженные матрицы.

а во-вторых, решать уравнения точно (до double точности) почти никогда не нужно. а значит нам нужны эффективные итеративные метод
Эффективные итерационные методы типа CG, для получения приемлемой точности требуют double вычислений. CG отличается от простых методов тем, что у него есть "контекст", а именно есть большая разница как мы попали в точку x1 (например это первое приближение, или допустим на некотором большом шаге n/1000), и точность вычислений при CG определяет скорость сходимости.
Я пробовал решать решать систему МКЭ на GPU(не суть), с помощью CG (+ предобуславливатель Якоби) - так вот, на некоторых задачах single-precision вообще не сходилось к заданной точности, а на некоторых сходилось через количество итераций большее чем для double-precision.
Далее, нужно ещё учитывать то, что когда мы уменьшаем точность в матрице коэффициентов (допустим используем float вместо double) и используем float для вычислений, то мы не просто уменьшаем точность результата каждой итерации, а ещё и рассматриваем немного другую задачу.
упираются в память
да, верно - упирается в память. Причём для основных форматов разряженных матриц примерно половина bandwidth это коэффициенты матрицы, другая половина - индекс колонки для каждого коэфф. матрицы (то есть обычно uint32, и от точности вычислений не зависят) (ещё есть несколько других векторов, но они значительно меньше). То есть, при переходе float->double, необходимый bandwidth не увеличится в два раза, поэтому большого смысла в использовании float нет.

а значит нам нужны эффективные итеративные методы. для них главные операции - умножение матрицы на вектор и применение прекондишинера имеют сложность n^2 и упираются в память, гигафлопсов особо не достичь.

Уравнения получаемые в результате применения МКЭ, решаются не только с помощью итерационных методов, но и с помощью прямых.
И у каждых своя область применения. Прямые обычно применяют тогда, когда возможно (хватает памяти) - потому что получается намного быстрее чем итерационными. итерационные в остальных случаях (так как требуют гораздо меньше памяти. при итерационных методах, можно даже не хранить матрицу коэффициентов явно).
Так вот, прямые методы упираются обычно не в память, а во флопсы.
И у меня есть предположение, что коды используемые в hpl, отчасти могут быть использованы при решении разряженных слау прямыми методами.

для них главные операции - умножение матрицы на вектор и применение прекондишинера имеют сложность n^2 и упираются в память, гигафлопсов особо не достичь.
Ну вот вы сами себе и ответили - HPL решает Ax=b с плотными матрицами, чтобы как раз померятся flops'ами. Причём я бы сказал даже не flops'ами, а dgemm/s'ами. Делать тупо dgemm не интересно, мерить тупо flops'ы - тоже, поэтому и выбрали что-то более менее практическое.