Aprendiendo SCILAB - Eliminación de Gauss con pivotación parcial

El siguiente código implementa el algoritmo de eliminación de Gauss con pivotación parcial para resolver el sistema de ecuaciones

Ax = b

// *****************************************************************************
// * FUNCION : gaussPivPar
// * PARAMS : 
// *            A - [n x n] Double
// *            b - [n x 1] Double
// * RETURN :
// *            < struct >
// *                U - [n x n] Double
// *                x - [n x 1] Double
// *
// * Calcula la solucion del sistema de ecuaciones dado por Ax = b mediante el 
// * algoritmo de Gauss con pivotacion parcial. 
// * Esta funcion no solo retorna la solucion sino que tambien retorna la matriz
// * triangular superior que se forma con este algoritmo.
// *****************************************************************************
function rpta = gaussPivPar(A, b)

    [m,n] = size(A)     // Dimesiones de la matriz de entrada
    
    if (m <> n) then    // Comprueba si la matriz es cuadrada
        error("La matriz no es cuadrada.")
    end
    
    // Inicia el algoritmo
    for i = 1:n-1
        // Pivotacion parcial
        [q,p] = max(A(i:n,i))
        p = p+i-1
        A = intercambiarFilas(A, i, p)
        b = intercambiarFilas(b, i, p)
        
        // Triangularizacion
        for j = i+1:n
            r = A(j,i)/A(i,i)
            A(j,:) = A(j,:) - r*A(i,:)
            b(j,:) = b(j,:) - r*b(i,:)
        end
    end
    
    // Sustitucion inversa
    for j = n:-1:1
        suma = 0
        for k = j+1:n
            suma = suma + A(j,k)*x(k)
        end
        x(j) = ( b(j) - suma ) / A(j,j)
    end
    
    rpta = struct("U", A, "x", x)
endfunction


// *****************************************************************************
// * FUNCION : intercambiarFila
// * PARAMS : 
// *            A - [m x n] Double
// *            i - [1 x 1] Integer
// *            j - [1 x 1] Integer
// *
// * Funcion de ayuda que intercambia las filas i y j de la matriz A.
// *****************************************************************************
function rpta = intercambiarFilas(A, i, j)
    temp = A(i,:)
    A(i,:) = A(j,:)
    A(j,:) = temp
    rpta = A
endfunction


// *****************************************************************************
// * EJEMPLO DE EJECUCION
// *
// * -->A = [0 2 1; 1 -2 -3; -1 1 2]
// *  A  =
// *  
// *     0.    2.    1.  
// *     1.  - 2.  - 3.  
// *   - 1.    1.    2.  
// *  
// * -->b = [-8; 0; 3]
// *  b  =
// *  
// *   - 8.  
// *     0.  
// *     3.  
// *  
// * -->sol = gaussPivPar(A, b)
// *  sol  =
// *  
// *    U: [3x3 constant]
// *    x: [3x1 constant]
// *  
// * -->sol.U
// *  ans  =
// *  
// *     1.  - 2.  - 3.   
// *     0.    2.    1.   
// *     0.    0.  - 0.5  
// *  
// * -->sol.x
// *  ans  =
// *  
// *   - 4.  
// *   - 5.  
// *     2. 
// *****************************************************************************