function Demidopivot(a,b)
%% THE METHOF of PRINCIPAL ELEMENTS
%
% See, Computational Mathematics
% B.P. Demidivich, I.A. Maron
% Third Printing, 1981, MIR.
% page 278.
%
% Implemented by Prof. Diomar Cesar Lobo PhD
% UFF - Volta Redonda, RJ, Brazil
% April 16th, 2013
%OBS: Matrix A must be Augmented matriz [a b]
%%
A = [a b]; atemp=A;
disp(A)
[n,m]=size(A);
%
%Initialize vector to keep track of column swaps.
for i = 1:n
col(i) = i;
end
%
%Do Demidovich Method of Principal Elements, pg. 278.
for k = 1:n-1
[n,m]=size(A);
max = 0;
isav = k;
jsav = k;
for i = k:n % Find pivot a(p,q)
for j=k:n
mtest = A(i,j);
if (abs(mtest)) > max
max = abs(A(i,j));
isav = i;
jsav = j;
end
end
end
out=sprintf('valor amax=%f, imax=%d, jmax=%d\n',max,isav,jsav);disp(out);
%
% calculate the multipliers mi=-a(i,q)/a(p,q)
for i =k:n
if (i ~= isav)
factorm(i)=-A(i,jsav)/A(isav,jsav);
else
if (i == isav)
factorm(i)=-A(isav,jsav)/A(isav,jsav);
end
end
end
% calculate the M matrix
for i = k:n
if(i ~= isav)
for j = k:m
atemp(i,j)=A(i,j)+A(isav,j)*factorm(i);
end
end
end
%
atemp
A = atemp;
% rebuild matrix in upper triangular form
if (isav ~= k)
%switch rows.
for j = 1:m
temp = A(k,j);
temp1 = A(isav,j);
A(k,j) = temp1;
A(isav,j) = temp;
end
disp('swap rows');
disp(A);
end
if (jsav ~= k)
%switch columns.
tempcol = col(jsav);
col(jsav) = col(k);
col(k) = tempcol;
for i=1:n
temp = A(i,k);
temp1 = A(i,jsav);
A(i,jsav) = temp;
A(i,k) = temp1;
end
disp('swap columns');
disp(A);
disp('variables now in order: ');
disp(col');
end
atemp = A;
%
end %End of matrix M sequence
%
if (A(n,n) == 0) %Check for singular matrix.
disp('matrix of coefficients is singular')
disp(A);
return;
end
%
%Do back substitution.
for k = 1:m-n
x(k,n) = A(n,n+k)/A(n,n);
end
for k = 1:m-n
for i = n-1:-1:1
sum(k) = 0;
for j = i+1:n
sum(k) = sum(k) + A(i,j)*x(k,j);
end
x(k,i) = A(i,n+k) - sum(k);
if (A(i,i) == 0) %Check for singular matrix.
disp('matrix of coefficients is singular');
disp(A);
break;
end
x(k,i) = x(k,i)/A(i,i);
end
end
%
%Display upper triangular form.
disp('upper triangular form:');
disp(A);
%
%Swap solution vector if columns have been exchanged.
%
for i = 1:n
if(col(i) ~= i)
disp('current order of variables: ');
disp(col');
pause;
disp('switching back to original order ');
disp(' ');
break;
end
end
%
%Display solution.
disp('solution:');
%
for j = 1:m-n
for i = 1:n
xo(j,col(i)) = x(j,i);
end
end
disp(xo');
return;
%
