把矩陣化成簡化行階梯型矩陣的函數(shù)rref 在北太天元中已經(jīng)實現(xiàn)為內(nèi)置函數(shù)了。下面是一個腳本實現(xiàn)的版本，可以供寫腳本的同志們參考。

為了和內(nèi)置的rref區(qū)分，這里的函數(shù)名改成了 rrefEx, 這個腳本參考了octave的rref.m的實現(xiàn)。

北太天元作簡化行階梯型矩陣的rrefEx.m 的內(nèi)容如下

% 把矩陣 A 化成 簡化行階梯形矩陣

% [A, k] = rrefEx (A, tol)

% 右邊的A 是輸入的矩陣，可以是single或者double，但是必須是二維矩陣

% tol 是忍量，矩陣值比tol小的數(shù)會被置為0

% 左邊的A 是輸出的矩陣，是簡化行階梯形矩陣

% k 是主元所在的列號

% 例:

% a = [1];

% [r k] = rrefEx (a);

% a = [1 3; 4 5];

% [r k] = rrefEx (a);

% a = [1 3; 4 5; 7 9];

% [r k] = rrefEx (a);

% a = [1 2 3; 2 4 6; 7 2 0];

% [r k] = rrefEx (a);

% a = [1 2 1; 2 4 2.01; 2 4 2.1];

% tol = 0.02;

% [r k] = rrefEx (a, tol);

% tol = 0.2;

% [r k] = rrefEx (a, tol);

function [A, k] = rrefEx (A, tol)

???%如果輸入的參數(shù)個數(shù)小于1，那么輸出幫助信息

?if (nargin < 1)

??????help rrefEx;

?end

?if (ndims (A) > 2)

???error ("rrefEx: A 必須是一個二維矩陣");

?end

?[rows, cols] = size (A);

?if (nargin < 2)

???if (isa (A, 'single'))

?????tol = eps ('single') * max (rows, cols) * norm (A, inf ('single'));

???else

?????tol = eps * max (rows, cols) * norm (A, inf);

???end

?end

?used = zeros (1, cols); %用來記錄用過的列

?r = 1;?

?for c = 1:cols

???% 第c列的 從第r行到最后一行 的主元 （絕對值最大的元素)

???[m, pivot] = max (abs (A(r:rows,c))); % m 是主元的大小

???pivot = r + pivot - 1; % pivot 是主元所在的行號

???if (m <= tol)

?????% 如果上面得到的主元m很小，那么跳過第c列，

???????% 同時把這些小元素都強制修改成真正的0

?????A(r:rows, c) = zeros (rows-r+1, 1);

???else

?????% 記錄主元所在的列號

?????used(1, c) = 1;

?????% 交換當前行和主元所在的行

?????A([pivot, r], c:cols) = A([r, pivot], c:cols);

?????% 現(xiàn)在主元所在的行已經(jīng)放在了當前行r了，再用一個常數(shù)1/A(r,c)乘以這一行

?????????% 從而把主元 標準化為 1.

?????A(r, c:cols) = A(r, c:cols) / A(r, c);

?????% 把當前列的主元上下方的元素都消成0

?????ridx = [1:r-1, r+1:rows];% 1:r-1 是主元上方的行號，r+1:rows是主元下方的行號

?????A(ridx, c:cols) = A(ridx, c:cols) - A(ridx, c) * A(r, c:cols);

?????%檢查是否完成了所有行

?????if (r++ == rows)

???????break;

?????end

???end

?end

?k = find (used);

end