Home > K25104 > LinearSystems > Gaussian_total.m

Gaussian_total

PURPOSE ^

Performs Gaussian elimination with total pivoting on A to

SYNOPSIS ^

function [ U, c , P ] = Gaussian_total( A,b )

DESCRIPTION ^

 Performs Gaussian elimination with total pivoting on A to 
   transform the system Ax=b to U(Px)=c.
   Note the permutation matrix P stores the column swaps.
 Input arguments:
   A, square matrix
   b, column vector
 Output arguments:
   U, square upper triangular matrix
   c, column vector
   P, square permutaion matrix

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ U, c , P ] = Gaussian_total( A,b )
0002 % Performs Gaussian elimination with total pivoting on A to
0003 %   transform the system Ax=b to U(Px)=c.
0004 %   Note the permutation matrix P stores the column swaps.
0005 % Input arguments:
0006 %   A, square matrix
0007 %   b, column vector
0008 % Output arguments:
0009 %   U, square upper triangular matrix
0010 %   c, column vector
0011 %   P, square permutaion matrix
0012 
0013 U=A; c=b;       % initialize U,c
0014 [m,n]=size(A); % finding the size of A
0015 P=eye(n);   % initialize permutation matrix to store column permutations
0016 % check user inputs
0017 if n~= m
0018     error('input is not a square matrix');
0019 elseif size(b,1) ~= n
0020     error('input dimensions do not match');
0021 end
0022 
0023 for i=1:n;  % considering the rows in turn
0024     % choose the element with the greatest absolute value over the bottom
0025     % right (n-i)x(n-i) submatrix
0026     [c_value,c_ind]=max(abs(U(i:n,i:n))); 
0027     [~,r_ind]=max(abs(c_value));
0028     c_ind=c_ind(r_ind);
0029     r_ind=r_ind+i-1;
0030     c_ind=c_ind+i-1;  % index of largest element is (r_ind,c_ind)
0031     
0032     if abs(U(r_ind,c_ind))<1e-12;
0033         error('matrix is singular, there is no unique solution');
0034     else
0035         U([r_ind,i],:)=U([i,r_ind],:);  % interchange rows
0036         c([r_ind,i],1)=c([i,r_ind],1);  % perform same row swap on the
0037                                         % right hand side
0038     
0039         U(:,[c_ind,i])=U(:,[i,c_ind]);  % interchange columns
0040         P(:,[c_ind,i])=P(:,[i,c_ind]);  % perform same column swap on P
0041     end
0042     
0043     for j=i+1:n;
0044         temp=(U(j,i)/U(i,i));       % eliminate elements below the diagonal
0045         U(j,:)=U(j,:)-U(i,:)*temp;  % by subtracting a linear combination
0046                                     % of rows
0047         
0048         c(j)=c(j)-c(i)*temp;        % perform the same row operation on the
0049                                     % right hand side
0050     end
0051 end
0052 
0053 end
0054

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005