import MatrixDemo;

public class Gauss
{
    public static void gaussian(MatrixDemo a, MatrixDemo b,  MatrixDemo x) 
    {
        int i, j, n;
        n= a.nrows;                      // Número de incógnitas
        MatrixDemo q= new MatrixDemo(n, n+1);
        for (i=0; i < n; i++) {
            for (j=0; j < n; j++)                   // Formar la matriz q
                q.data[i][j] = a.data[i][j]; // q(i,j)= a(i,j)
            q.data[i][n] = b.data[i][0];     // q(i,n)= b(i,0)
        }
        forward_solve(q);                  // Realizar eliminación gaussiana
        back_solve(q);                     // Realizar sustitución inversa
        
        for (i=0; i<n; i++)
            x.data[i][0] = q.data[i][n];   // x(i,0)= q(i,n)
    }
    // Función de eliminación directa
    private static void forward_solve(MatrixDemo q)
    { 
        int i, j, k, maxr, n;
        double t, pivot;
        n= q.nrows;
        
        for (i=0; i < n; i++) {		// Buscar fila con máx. elemento en esta
            maxr= i;			// columna en la diagonal o por debajo
            for (j= i+1; j < n; j++)
                if (Math.abs(q.data[j][i]) > Math.abs(q.data[maxr][i]))
                    maxr= j;
            if (maxr != i)		// Si la fila no es la fila actual, saltar
                for (k=i; k <= n; k++)
                { t= q.data[i][k];                     // t= q(i,k)
                  q.data[i][k] = q.data[maxr][k];  // q(i,k)= q(maxr, k)
                  q.data[maxr][k] =  t;                 // q(maxr, k)= t
                }
            for (j= i+1; j <n; j++)         // Calcular relación de giro
            { pivot= q.data[j][i] / q.data[i][i]; // q(j,i)/q(i,i)
              for (k= n; k >=i; k--)
                  q.data[j][k] -= q.data[i][k] * pivot;
                  // q(j,k) -= q(i,k)*pivot;	// Actualizar la fila j debajo de la diag.
            }
        }
    }
    
    private static void back_solve(MatrixDemo q)        // Función de sustitución inversa
    {                                               // a partir de la fila n-1
        int j, k, n;
        double t;                                   // t- temporal
        n= q.nrows;
        
        for (j=n-1; j >=0; j--)                     // Comenzar en la última fila
        {
            t= 0.0;
            for (k= j+1; k < n; k++)                // t += q(j,k)* q(k,n)
                t += q.data[j][k] * q.data[k][n];   
            q.data[j][n] = (q.data[j][n] - t)/q.data[j][j]; // q(j, n)= (q(j, n) -t)/q(j,j);
        }
    }
}