Hola,
yo tampoco tengo muy claro a que se tiene que llegar o cual es el óptimo, pero el procedimiento que hice fue:
Tomar las primeras n-filas de la matriz (que tiene dimensiones 2n y 2n). Teniendo una matriz de n filas por 2n columnas por el vector de dimension 2n de la forma (x_1,..x_n,y_1,...,y_n).
Al tener una matriz L1 inferior no singular me queda un sistema compatible determinado para las y_i que resuelvo mediante sustitución hacia adelante. En este punto ya pude encontrar el valor para todas las y_i.
Luego, tome las segundas n-filas y use el mismo procedimiento pero ahora cuando haga el producto me van a quedar los elementos de la matriz L2 por las variables x's y los elementos de B por los y's (que ya los tengo calculados del paso anterior). Paso ese producto de la matriz B por los y's para el lado del "termino independiente" (con los c's). Y vuelvo a hacer una sustitucion para adelante con los elementos de L2 usando los x_j como incognita y en el vector solución el c_j inicial restado por los B_ij * y_j.
Este procedimiento requiere 2 sustituciones hacia adelante, el cálculo del producto de B por el sub-vector "y" y las operaciones requeridas para "recalcular" el vector de termino independiente que se usa al resolver L2.
No se si es lo mejor, pero que la letra de matrices L no-singulares me dio la idea de arrancar por ese lado.
Saludos.
Agustin.
