Treating ties in a Gaussian Elimination process

Hey there, I have an M-File written for performing Gaussian elimination (shown below) however I can't figure out how to deal with ties in the matrix for example if there are two elements that are both the max value being called on how do I get the function to just select the first one it comes across when finding the max. Any help is appreciated.

function [Ang_update] = Gauss_elim(f, Df, k)

% [Ang_update] = Gauss_elim(f, Df, k)

% Using Gaussian Elimination to calculate the next approximation of the

% Newton Raphson iteration

% Input f is the column vector of functions containing the unknown variables

% Input Df is the Jacobian matrix of the functions in question

% Output Ang_update is the list of angles at the most recent iteration

% Version 1, 08/04/2024 by Conor Pierce, Student Number 21302251

if ~isvector(f)

error('Input f is a column vector of the system of equation to be solved.')

end

if ~isvector(k)

error('k is a column vector (8 x 1) containing the users initial estimates as entered in the Newton Raphson function.')

end

k0 = k; % Vector of initial estimates

F = f; % Initialise f(kn)

% Finding the size of the matrix

s = size(F);

n = s(1);

% Set Df as the Jacobian and formulate succesive steps

J = Df;

A = J;

B = Df*k0 - F;

% Begin set up for Gaussian elimination

A1 = A;

B1 = B;

% Set up identity matrix

Identity = eye(n);

Q_total = Identity;

% Total Pivoting

for m = 1 : n-1

active = m:n; % Identify active rows and columns

[Y, I] = max(abs(A1(active, active))); % Y is max value in column, I is position of that value

[~, I1] = max(Y); % Column to pivot is I1

row_num = I(I1);

col_num = I1;

p = 1:n;

p([m, row_num + (m-1)]) = p([row_num + (m-1), m]); % Row swap

q = 1:n;

q([m, col_num + (m-1)]) = q([col_num + (m-1), m]); % Column swap

Q_total = Q_total(:, q);

A1 = A1(p, q);

B1 = B1(p, :);

% Subtraction

I_vec = zeros(1, n);

I_vec(1, m) = 1;

L1 = Identity - ([zeros(m, 1); A1(m+1: n,m)]*I_vec)/A1(m, m);

A1 = L1*A1;

B1 = L1*B1;

end

% Set up column vector for outputs

Ang_update = zeros(n, 1);

for m = 0:n-1

% Back substitution

Ang_update(n-m, 1) = (B1(n-m, 1) - (A1(n-m, :)*Ang_update))/A1(n-m, n-m);

end

Ang_update = Q_total*Ang_update; % Reorder

end

Answer 1

[John D'Errico]

Um, max does exactly that. Where is there a problem?

x = [1 2 3 4 3 4 1];
[xmax,ind] = max(x)

xmax = 4

ind = 4

[Conor Pierce]

Yes, but does this hold for an 8 x 8 matrix?

[John D'Errico]

Ok, you are using complete pivoting, not just row pivoting as is common. Max, by default, finds the maximum in each column, and you don't want that.

A = randi(10,5,5)

A = 5x5

1 4 3 3 6 2 8 6 2 6 7 8 8 3 9 3 5 2 8 4 3 8 6 9 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Amax,ind] = max(A)

Amax = 1x5

7 8 8 9 9

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ind = 1x5

3 2 3 5 3

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

But if you do this:

[Amax,ind] = max(A,[],'all')

Amax = 9

ind = 20

It finds the largest overall element in the array. The index returned is a linear index, so the 9th element when the array is unrolled into a vector. You can recover the row and column indices of that by

[imax,jmax] = ind2sub(size(A),ind)

imax = 5

jmax = 4

Note that there were 2 instances of a 9. The first one was seen as A(5,4).