Нет, одного LU-разложения (хоть мыслите то Вы в правильном направлении) недостаточно. После того как оно (LU) выполнено, надо попросить соответствующий пакет решить систему линейных уражнений с матрицей правых векторов, представляющих из себя единичную матрицу. Ответ и будет обратной матрицей.

Должен признаться что сам с этими пакетами не работал, но принимал участие в написании подобных. Так в нашем пакете была возможность вызова процедуры решения системы линейных алгебраических уравнений методом полного (бывает ещё и incomplete) LU разложения с возможностью передачи в процедуру не одного вектора правой части, а целой матрицы правых частей. Физически это соответствует тому, что Вы хотите решить не одну, а сразу несколько систем линейных уравнений с одной и той же матрицей, но разными правыми частями. Вообще то я бы ожидал что такая возможность существует и в BLAS и в LAPACK. Посмотрите внимательнее. Если это так, то подавайте на вход в соответствующую процедуру в качестве матрицы правых частей единичную матрицу - тогда на выходе в качестве решения получите матрицу обратную Вашей матрице. Если такой возможности в этих пакетах нет (что было бы странно), то обязательно должен быть режим в котором процедура выполняет только LU-разложение. Вызывайте процедуру в этом режиме один раз - она разложит матрицу на произведение нижнедиагональной и верхнедиагональной. Затем вызываете эту же процедуру в цикле n раз, но уже в режиме только решения СЛАУ и в качестве вектора правой части на i-том цикле итерации передаёте в процедуру вектор все элементы которого за исключением i-того нули, а i-тый равен единице. Это другая альтернатива. Посмотрите повмнимательней на интерфеис процедуры и почитайте описание входных/выходных параметров - Вам должно стать всё понятно.

Да