Hola,
la factorización la pueden hacer ustedes. En Matlab/Octave, por ejemplo, tienen la función "lu" y "chol".
Sobre las matrices dispersas, el tiempo de ejecución va a depender de la potencia de tu máquina, la cantidad de no ceros y de cómo están dispuestos en la matriz, por lo que es posible que hayas elegido matrices que a partir de ese tamaño demoran demasiado. La idea es ver justamente cómo escala el tiempo de ejecución de los solvers directos vs. iterativos. Quizás podés intentar con matrices de un poco menor orden de nnz para que en la más grande tengas un tiempo razonable... por ejemplo nnz ~[5n,25n,50n]..?