Внимание! Об обновлениях смотри внизу этой страницы. Последнее обновление – от 18.04.2018.
Введение
Класс DMatrix разработан на языке C++ (в среде Borl a nd Builder 6) и предназначен для встраивания в исходный код с целью упрощения программирования операций с матрицами.
Класс позволяет использовать при программировании переопределенные операции: присвоение, сложение матриц, умножение матриц, умножение матрицы на число (справа). Например, код C++, использующий объекты данного класса, может выглядеть так:
A = B;
A = B + C;
A = B * C;
A = B * c;
где A , B и C – объекты класса, с – переменная типа float , double или long double .
Кроме того, класс содержит функции обращения матрицы, вычисления определителя и транспонирования:
A = B. Inverse (); -обращение матрицы B ;
d = B . det ( ); -вычисление определителя матрицы B ;
A = B. T ( ); -транспонирование матрицы B .
Значения ячеек матрицы имеют тип long double , что позволяет производить операции с большой точностью. В некоторых достаточно сложных прикладных задачах, построенных на итеративных алгоритмах, увеличение точности вычислений приводит к ускорению сходимости алгоритма, то есть даже ускоряет работу программы.
Важное преимущество нашего класса заключается в реализации подхода «динамическая матрица».
В прикладных задачах, связанных с обработкой потоков информации, встает задача оперативного формирования матриц на основе данных (сигналов), поступающих по нескольким каналам. При этом, во-первых, данные не должны пропадать из-за нехватки времени, они должны без временных задержек наполнять матрицы, позволяя математическому алгоритму быстро приступать к их обработке. Во-вторых, память компьютера должна использоваться рационально, чтобы избежать переполнения памяти при больших потоках информации.
Объект класса DMatrix можно изобразить в виде стакана:
Матрица, хранящаяся в данном объекте, имеет размерность m*n . Как правило, в случае работы с потоком данных число n – это количество переменных, содержащихся в потоке, и/или каких-то рассчитанных величин, например, производных по времени от переменных потока.
Рациональное использование памяти при работе с динамической матрицей подразумевает динамическое выделение памяти для строк матрицы (низ стакана) и динамическое освобождение памяти от использованных строк (верх стакана).
При операциях с матрицами, таких, как сложение, умножение, вычисление определителя и др., обработке подвергаются только актуальные данные, для указания на которые служит переменная k.
При обработке сигналов от устройств часто возникают задача формирования и решения в реальном времени систем уравнений. Для этой цели сделаны классы LSM и PROJECTION, работающие с нашими динамическими матрицами. В них реализованы алгоритмы решения систем линейных алгебраических уравнений с гибкой системой настройки параметров. Алгоритмы базируются на методе наименьших квадратов с «коэффициентом экспоненциального забывания» и проекционном решающем алгоритме. Настройки этих алгоритмов позволяют адаптировать их как для чистки входящих сигналов, так и для быстрого реагирования решения на скачки сигналов на входе. Таким образом, программист (или эксплуатант), подбирая эти настройки, должен найти оптимальный баланс между точностью решения и быстротой реагирования.
Переменные и функции класса DMatrix
Класс содержит переменные:
int m ; - количество строк матрицы
int n ; - количество столбцов матрицы
int k; - указатель на последнюю строку матрицы. При операциях в качестве актуальных строк матрицы используются строки массива с индексами(k-m), … , (k -1).
int M; - количество указателей на строки матрицы
long double **data; - указатель на двумерный массив значений ячеек матрицы
booloblom; - признак аварийного результата выполнения операции
Класс содержит функции:
long double __ fastcall det ( void ); - расчет определителя матрицы
DMatrix __ fastcall T ( void ); - транспонирование матрицы
DMatrix __ fastcall Inverse ( void ); - расчет обратной матрицы
void __ fastcall Ini ( m 0, n 0, M 0, k 0); - инициализация переменных объекта-матрицы и выделение памяти
void __ fastcall de _ allocate ( void ); - чистка памяти: удаление m строк матрицы из диапазона (k-m), … , (k -1) и указателей наэти строки
void __ fastcall Allocate ( int k 0); - выделение памяти для указанной строки
void __ fastcall Delete ( int k 0); - чистка памяти: удаление указанной строки
Работа с матрицами
Пример 1.
Рассмотрим простейший пример – сложение двух матриц:
Сначала нужно добавить в проект наш класс. Если используется среда Borland Builder , то необходимо добавить в проект файл dmatrix . cpp и подключить его в unit ’е, в котором мы будем работать с матрицами; для этого напишем в header -файле директиву#include " d matrix.h " .
Добавим на форму кнопку Button 1 и надпись Label1. В функцию Button1Click поместим следующий код.
Объявим 3 матрицы:
DMatrixA, B, C; (слагаемые и результат, C = A + B ),
опишем свойства этих матриц и выделим память для строк и столбцов матриц:
A.m = 2;
A.n = 2;
A.k = 2;
A.M = 2;
A.oblom = false;
A.data = new long double*;
for (inti = 0; i < A.k; i++)A.data[i] = new long double;
То же самое можно написать более компактным образом, с помощью функции Ini:
A . Ini (2, 2, 2, 2);
В результате выполнения функции Ini (m 0, n 0, M 0, k 0) задаются значения переменных A.m, A.n, A.M, A.k, выделяется память для A.M указателей на строки и выделяется память для A.k строк.
Примечание по работе Ini :
Иногда нет смысла сразу выделять память под строки. В этом случае можно написать:A.Ini(2, 2, 2).После выполнения функции Ini с тремя аргументами будет присвоено A.k = 0, о чем не стоит забывать, т.к. A.k = 0 указывает на матрицу из 0 строк; если мы хотим указать на матрицу из A.m строк, то A.k должен быть не меньше A . m .
Можно также сделать двухместный вызов этой функции:A.Ini(2, 2), в этом случае не будет выделена память для указателей на строки и будет присвоено A. M = 0.
Так же поступим с матрицами B и C :
B . Ini (2, 2, 2, 2);
C . Ini (2, 2, 2, 2);
Зададим значения ячеек матриц A и B:
A.data = 1;
A.data = 2;
A.data = 3;
A.data = 4;
B.data = 5;
B.data = 6;
B . data = 7;
B . data = 8;
и, наконец, сложим матрицы:
C = A + B ;
В результате сложения в массиве C.data появится требуемый результат, который можно вывести на экран. Например:
Label 1-> Caption = FloatToStr ( C . data ); - на форме отобразится значение ячейки (1-я строка, 1-й столбец), равное 6.
После получения результата и его выгрузки освободим память:
if(A.data)A.de_allocate();
if(B.data)B.de_allocate();
if(C.data)C.de_allocate();
Важное замечание по поводу работы функции de_allocate():эта функция удаляет только строки в диапазоне (k-m), … , (k -1)(имеются в виду индексы массива data в C++ коде, то есть нумерация идет с 0). То есть «рабочую» часть нашего стакана.Перед запуском de_allocate() надо быть уверенным, что все прочие строки, лежащие вне данного диапазона, были удалены ранее. В нашем примере k = m, поэтому функцию выполнять можно; вопрос постепенного удаления «отслуживших» строк рассмотрим в следующем примере.
Исходный код рассмотренного примера – .
Пример 2.
Предположим, что требуется оценивать некоторую характеристику потока измерений, поступающих в компьютер с внешних устройств. Например, это данные с N датчиков.
В каждый момент времени имеем N значений, из которых можно составить строку матрицы, то есть можно сказать, что на вход приходят строки. Пусть характеристика, которую надо вычислять – детерминант. Для того, чтобы его вычислить, нам нужно иметь N строк, т.к. матрица должна быть квадратной.
Опишем матрицу и выделим память (только для указателей на строки):
X . Ini ( N , N , 10000);// Не указали значение k - не выделилось место под строки
X . k = 1;
while ( X . k <= X . M )
{
X . Allocate ( X . k - 1);// Выделяем память для строки
if ( X . k > X . m ) X . Delete ( X . k - X . m - 1);// Удаляем последнюю использованную, но пока не удаленную строку
// Заполняем новую строку (симуляция сигналов)
for(i = 0; i < N; i++)
X.data[i] = (((rand() % RAND_MAX)+1.0)/RAND_MAX);// Равномерно распределенная на случайная величина
if ( X . k >= X . m ) det = X . det ();// Появилось m строк - можно считать определитель
X . k ++;
}
Полный текст пр ограммы - . В исходном коде добавлено умножение матрицы на себя, дискриминант берется от матрицы-произведения. Отображение результатов на форме сделано с помощью таймера.
Исходя из значения X.M можно оценить, сколько раз в секунду Ваш компьютер находит значение определителя матрицы размера N * N (всего определитель вычисляется (X.M- N +1) раз). Увеличивая размер матрицы, можно получить представление о возможностях алгоритма.
Важное замечание: в этом примере мы освобождали и выделяли память для строк на каждом шаге. Для повышения быстродействия иногда имеет смысл перераспределять память с интервалом в несколько шагов, большими кусками. То есть через несколько шагов добавлять / удалять сразу по несколько строк матрицы.
Так как функция de_allocate() освобождает память от строк с индексами из диапазона(k-m), … , (k -1), важно перед запуском функции указать правильное значение k. Для этого в исходном коде после главного цикла указано:X.k--; , чтобы сбросить последнее увеличение X.k на единицу.
Динамическая идентификация состояния системы
Для демонстрации возможностей работы с классом предлагается пример (с исходниками) .
При запуске программы (matrix.exe) открывается форма, на которой можно задавать матрицы и производить простейшие операции с ними. В исходном тексте можно увидеть, каким образом используется наш класс при программировании на C++.
При нажатии на кнопку «Динамические операции» открывается форма для решения следующей задачи.
Пусть у нас есть 4 сигнала, поступающих на вход:X 1 , X 2 , X 3 и Y . Мы предполагаем, что эти 4 переменные связаны между собой линейным уравнением:
C 1 * X 1 + C 2 * X 2 + C 3 * X 3 = Y
В качестве одного или нескольких входных параметров может использоваться как физический сигнал, так и расчетные величины, например, производные или отфильтрованные данные.
Требуется оценить в каждый момент времени коэффициенты состояния системы C 1 , C 2 , C 3 , то есть решить уравнение относительно C i .
В идеальных условиях, при отсутствии шума и при постоянных X 1 , X 2 , X 3 и Y в каждый момент времени мы получаем одно и то же уравнение, имеющее бесконечное множество решений. На практике мы имеем шум и изменяющиеся значения сигналов, в результате чего получается система из большого количества несовместных уравнений.
Для приближенного решения несовместной системы существует несколько вычислительных подходов; в нашей программе реализована комбинация 2-х методов:МНК (Метод наименьших квадратов) с «коэффициентом экспоненциального забывания» и обобщенный проекционный алгоритм.
Использование подхода МНК обусловлено необходимостью получения на каждом шаге невырожденной матрицы для последующего решения системы уравнений проекционным алгоритмом. Для увеличения гибкости настройки вводится коэффициент экспоненциального забывания W 2 , смысл которого состоит в усреднении вновь полученной методом МНК системы с аналогичной системой, полученной на предыдущем шаге. Чем больше значение W 2 , тем больший вес имеют предыдущие значения. Коэффициент принадлежит интервалу ;
Если процесс – единственный, этого делать не надо, по умолчанию Thr = 0.
То же касается и инициализации объектов LSM и PROJECTION. В этих классах в функциях инициализации добавлена еще одна переменная (по умолчанию также равная 0):
LSMLsm;
Lsm .Ini (…, …, [ номер процесса ]);
PROJECTIONProjection ;
Projection . Ini 0(…, [номер процесса]);
Все исходники, доступные к скачиванию с данной страницы, обновлены с учетом этих изменений.
Безусловно, по входным данным можно восстановить матрицы, тривиальным образом их перемножить и вывести результат в заданном виде. Как пример плохого стиля программирования: наивный алгоритм приводит к почти двукратному перерасходу потребляемой памяти и существенно уступает в производительности. Существует другое решение, прийти к которому можно путем последовательного приспособления стандартных матричных операций к формату входных данных.
Наибольшую сложность, как вычислительную, так и идейную, представляет умножение матриц, в данном случае — возведение в квадрат. Нижеприведенная последовательность шагов позволит перемножать симметричные матрицы, используя их линейное представление в формате исходных данных.
Замечание 1
Подпространство симметричных матриц незамкнуто относительно умножения: произведением двух симметричных матриц может быть несимметричная матрица.
Пример
\left(\begin{array}{ccc} 1 & 1 & 3 \\ 1 & 4 & 5 \\ 3 & 5 & 0\end{array} \right) \cdot \left(\begin{array}{ccc} 4 & 0 & 0 \\ 0 & 4 & 3 \\ 0 & 3 & 2 \end{array} \right) = \left(\begin{array}{ccc} 4 & 13 & 9 \\ 4 & 31 & 22 \\ 12 & 20 & 15 \end{array} \right)
Замечание 2
Степень симметрической матрицы также является симметрической матрицей. Доказательство основано на представлении матрицы как представления линейного оператора и на свойствах эрмитовых операторов.
Во всех дальнейших выкладках подразумевается, что матрица представлена линейным массивом из \frac{n(n+1)}{2} элементов.
Для начала, заметим, что элемент c_{i,j} матрицы C=A^{2}, равен скалярному произведению (как векторов в стандартном базисе) i-ой строки матрицы A на j-ую её строку (в силу того, что в симметричной матрице j-ая строка совпадает с j-м столбцом). Следовательно, для возведения в степень симметричной матрицы необходимо и достаточно реализовать операцию скалярного перемножения двух её строк.
Тогда следует понять, как по данному представлению матрицы получить её i-ую строку. Для удобства, выпишем имеющиеся элементы в виде полной матрицы. Заметим, что первым элементом i-ой строки будет i-ый элемент первой строки, и обобщим это наблюдение. Обозначим позицию текущего интересующего нас элемента i-ой строки как j. Если j < i, то следует выбрать i-ый элемент j-ой строки, иначе следует выбрать все элементы j-ой строки, лежащие правее данного. Графически можно интерпретировать алгоритм таким образом: начиная с i-го элемента первой строки, спускаться вертикально вниз по матрице до тех пор, пока не будет достигнута i-ая строка, далее — двигаться вправо до конца строки. Наглядность алгоритма позволяет легко воплотить его программно.
Следует только пронаблюдать, как изменяются расстояния между элементами, лежащими один под другим, при движении вниз по строкам матрицы, что позволит легко перенести алгоритм на линейное представление верхней треугольной матрицы. Предлагаем читателю самому проделать это несложное, но наглядное упражнение, для лучшего понимания алгоритма.
Теперь можно перейти непосредственно к реализации.
Объявление двумерного массива в СИ имеет следующий синтаксис:тип имя[размер №1][размер №2];
Размеры двумерного массива в СИ указываются в отдельных парных квадратных скобках после имени и могут быть любыми положительными целочисленными значениями. На практике принято значение первой размерности называть строками, а второй – столбцами. Как и в случае одномерного массива, в стандарте С89 регламентируется, что размеры двумерного массива должны быть целочисленными константами.
Стандарт С99 допускает объявление динамических двумерных массивов путём использования выражений при указании размеров матрицы, если в это выражение входят значения определенных ранее переменных (выражение должно иметь положительный целочисленный результат). Например:
int n,m;
printf("Введите размеры матрицы: ”);
scanf("%d %d”,&n,&m);
double a[n][m]; При объявлении двумерного массива в СИ допускается производить инициализацию значений элементов матрицы:
тип имя[размер №1][размер №2] = {
{значение № 11, ... значение № 1N},
...
{значение № M1, ... значение № MN}
}; Примеры объявлений с инициализацией:
int a = { //Объявлен двумерный массив
{1,2,3,4}, // 1 2 3 4
{5,6}}; // 5 6 0 0
Double b = { //Объявлен двумерный массив
{1.0, 2.0, 3.0, 4.0, 5.0}, // 1 2 3 4 5
{6.0, 7.0} // 6 7 0 0 0
}; // 0 0 0 0 0
int a = {{1,2,3,4,5},{6,7,8,9,0}}; Допускается не указывать количество строк в двумерном массиве (указываются пустые квадратные скобки). В таком случае размер массива будет определен по числу инициализирующих значений строк. Количество столбцов матрицы всегда необходимо указывать. Например:
double b = {{1,2,3,4},{5,6,7,8}}; Объявление константных матриц (значения их элементов изменить нельзя) начинается с ключевого слова const, за которым следует объявление матрицы с инициализацией. Пример:
const int matrix = {
{1,2,3,4,5},
{6,7,8,9}
}; Обращение к элементу матрицы осуществляется путем указания имени матрицы, а после имени в отдельных парных квадратных скобках индексы элемента (строка и столбец):
имя[строка][столбец] Индексация в языке СИ начинается с нуля, поэтому для матрицы размером, например, пять строк и десять столбцов правильными будут индексы строк от нуля до четырех, а столбцов – от нуля до девяти включительно.
Каждый отдельный элемент матрицы может рассматриваться как простая переменная и, соответственно, выступать в выражениях в качестве RValue или LValue значений.
Ввод и вывод матриц в языке СИ осуществляется поэлементно. Так как матрица имеет двойную размерность, то ввод и вывод осуществляется во вложенных циклах. Например:
double a;
for(int i=0;i<5;i++)
for(int j=0;j<10;j++)
scanf("%lf”,&a[i][j]);
...
for(int i=0;i<5;i++){
for(int j=0;j<10;j++)
printf("%8.2lf\t”,a[i][j]);
printf("\n”);
} Присвоение матрицы матрице также осуществляется поэлементно. Например, необходимо присвоить целочисленную матрицу x целочисленной матрице y. Фрагмент программы:
int x, y;
...
for(int i=0;i<5;i++)
for(int j=0;j<10;j++)
y[i][j] = x[i][j];
... В языке СИ допускается создание массивов размерностью три и более(т.е трехмерных, четырехмерных и т.д.). Например, объявление трёхмерного целочисленного массива с инициализацией будет иметь вид:
int a={ //это трехмерный массив
{{1,2},{3,4}},
{{5,6},{7,8}}
}; Ввод, вывод и прочая обработка такого массива осуществляется аналогично обработке двумерного массива, только уже в трех вложенных
До этого момента, мы рассматривали только одномерные массивы, то есть, к элементу массива мы обращались через один индекс. Однако, массивы могут быть и двумерными и трехмерными и, даже, n-мерными. Многомерные массивы — это массивы, у которых есть более одного индекса. Вместо одной строки элементов, многомерные массивы можно рассматривать как совокупность элементов, которые распределены по двум или более измерениям. Вот так, например, можно визуализировать двумерный массив:
В этом примере изображен двумерный массив размером 3*5, 3 — строки и 5 столбцов. Объявление двумерного массива почти ничем не отличается от объявления одномерного, за исключением того, что при объявлении двумерного массива, нужно указывать размер каждого измерения в квадратных скобочках. Например, давайте объявим двумерный массив размером 8*8, это размер поля для стандартных шашек — 8 строк и 8 столбцов:
Int checkers; // двумерный массив
То есть, двумерный массив хорошо подходит для хранения информации на шашечном поле. Также двумерный массив можно легко использовать для хранения информации о любой другой игре — шахматы, крестики нолики, сапер и т. д. Чтобы получить доступ к любому элементу такого массива, нужно воспользоваться двумя значениями — индексами, первый индекс — это номер строки, а второй — номер столбца. Все выше сказанное относится и к n-мерным массивам. Хотя, уже 4-х мерные массивы сложновато визуализировать. Присваивать значения элементам массива очень просто, вот пример:
// присваиваем первому элементу массива значение - 5 myArray = 5;
Пример присваивания значения элемента двумерного массива:
// присваиваем первому элементу массива значение - 5 myArray = 10;
В этом примере мы присвоили значение 10 элементу двумерного массива myArray , который находится во второй строке и в 4-м столбце. Визуально это выглядит так:
[__][__][__][__][__] [__][__][__][__] [__][__][__][__][__]Как видите, все просто, главное помните, что нумерация строк и столбцов всегда начинается с 0. То есть, я еще раз хочу вам напомнить, что вы никогда не должны пытаться записать данные после последнего элемента массива, например, когда у вас есть массив размером — 10 элементов и вы пытаетесь присвоить значение элементу с индексом . Память для массива была выделена только для десяти элементов, (индексы от 0 до 9), поэтому элемента с индексом 10 просто не существует. В таком случае, запись в оперативной памяти может привести к непредсказуемым последствиям — например, вы можете в конечном итоге испортить работу параллельно запущенной программы. Однако, как правило, операционная система не позволит такого рода безрассудное поведение и приведет к краху программы, если та попытается получить доступ к нераспределенной памяти.
Давайте рассмотрим практический пример использования массивов в программах:
#include
Сразу смотрим результат работы программы:
Вот такой массив у нас получился: =0 =0 =0 =0 =0 =0 =0 =0 =0 =1 =2 =3 =4 =5 =6 =7 =0 =2 =4 =6 =8 =10 =12 =14 =0 =3 =6 =9 =12 =15 =18 =21 =0 =4 =8 =12 =16 =20 =24 =28 =0 =5 =10 =15 =20 =25 =30 =35 =0 =6 =12 =18 =24 =30 =36 =42 =0 =7 =14 =21 =28 =35 =42 =49
В этом примере, мы сначала заполняем двумерный массив произведением его индексов, строки 8 — 11 . А потом выводим на экран его содержимое, строки 13 — 20 .
Если вы хотите объявить указатель на массив, то вы не должны использовать операцию взятия адреса — & , вот пример:
Char *ptrArray; char myString; ptrArray = myString; // указателю присваиваем адрес первого элемента массива myString без использования &
В то время как с обычными переменными, этого сделать нельзя, пример:
Int *ptrNumber; int number; ptrNumber = &number; // обязательно используем оператор - &
На первый взгляд, тот факт, что массивы работают как указатели, может запутать вас, но это специфика языка программирования, это нужно просто запомнить.
P.S.: После прочтения статьи, рекомендую вам отвлечься от программирования. Хороший способ это сделать — интересная игра. Достаньте смартфон и запустите свою любимую игру. А вот и полезная ссылка для фанатов андроида — скачать игры для андроид . Игр очень много и все интересные, а главное — бесплатные.
Аннотация: Приводятся правильные и неправильные способы реализации матриц и многомерных массивов на языке Си. Работа с матрицами иллюстрируется на примере приведения матрицы к ступенчатому виду методом Гаусса. Рассматриваются методы работы с файлами, использующие функции ввода-вывода из стандартной библиотеки ANSI. Приводятся способы работы с символами и текстовыми строками с помощью функций стандартной библиотеки. Материал иллюстрируется примерами, включающими программу "wc" подсчета символов, слов и строк в файле и программу "Записная книжка", которая позволяет находить телефон человека по его имени, а также сохранять и модифицировать содержимое книжки.
Представление матриц и многомерных массивов
Специального типа данных матрица или многомерный массив в Си нет, однако, можно использовать массив элементов типа массив . Например, переменная a представляет матрицу размера 3x3 с вещественными элементами:
Элементы матрицы располагаются в памяти последовательно по строкам: сначала идут элементы строки с индексом 0, затем строки с индексом 1, в конце строки с индексом 2 (в программировании отсчет индексов всегда начинается с нуля, а не с единицы!). При этом выражение
где i -- целая переменная , представляет собой указатель на начальный элемент i -й строки и имеет тип double* .
Для обращения к элементу матрицы надо записать его индексы в квадратных скобках, например, выражение
представляет собой элемент матрицы a в строке с индексом i и столбце с индексом j . Элемент матрицы можно использовать в любом выражении как обычную переменную (например, можно читать его значение или присваивать новое).
Такая реализация матрицы удобна и максимально эффективна с точки зрения времени доступа к элементам. У нее только один существенный недостаток: так можно реализовать только матрицу, размер которой известен заранее. Язык Си не позволяет описывать массивы переменного размера, размер массива должен быть известен до начала работы программы еще на стадии компиляции.
Пусть нужна матрица , размер которой определяется во время работы программы. Тогда пространство под нее надо захватывать в динамической памяти с помощью функции malloc языка Си или оператора new языка C++. При этом в динамической памяти захватывается линейный массив и возвращается указатель на него. Рассмотрим вещественную матрицу размером m строк на n столбцов. Захват памяти выполняется с помощью функции malloc языка Си
double *a; . . . a = (double *) malloc(m * n * sizeof(double));
или с помощью оператора new языка C++:
double *a; int m, n; . . . a = new double;
При этом считается, что элементы матрицы будут располагаться в массиве следующим образом: сначала идут элементы строки с индексом 0 , затем элементы строки с индексом 1 и т.д., последними идут элементы строки с индексом m - 1 . Каждая строка состоит из n элементов, следовательно, индекс элемента строки i и столбца j в линейном массиве равен
(действительно, поскольку индексы начинаются с нуля, то i равно количеству строк, которые нужно пропустить, i * n - суммарное количество элементов в пропускаемых строках; число j равно смещению внутри последней строки). Таким образом, элементу матрицы в строке i и столбце j соответствует выражение
Этот способ представления матрицы удобен и эффективен. Его основное преимущество состоит в том, что элементы матрицы хранятся в непрерывном отрезке памяти. Во-первых, это позволяет оптимизирующему компилятору преобразовывать текст программы, добиваясь максимального быстродействия; во-вторых, при выполнении программы максимально используется механизм кеш -памяти, сводящий к минимуму обращения к памяти и значительно ускоряющий работу программы.
В некоторых книгах по Си рекомендуется реализовывать матрицу как массив указателей на ее строки, при этом память под каждую строку захватывается отдельно в динамической памяти:
double **a; // Адрес массива указателей int m, n; // Размеры матрицы: m строк, n столбцов int i; . . . // Захватывается память под массив указателей a = (double **) malloc(m * sizeof(double *)); for (i = 0; i < m; ++i) { // Захватывается память под строку с индексом i a[i] = (double *) malloc(n * sizeof(double)); }
После этого к элементу a ij можно обращаться с помощью выражения
Несмотря на всю сложность этого решения, никакого выигрыша нет, наоборот, программа проигрывает в скорости! Причина состоит в том, что матрица не хранится в непрерывном участке памяти, это мешает как оптимизации программы, так и эффективному использованию кеш -памяти. Так что лучше не применять такой метод представления матрицы.
Многомерные массивы реализуются аналогично матрицам. Например, вещественный трехмерный массив размера 4 x 4 x 2 описывается как
обращение к его элементу с индексами x , y , z осуществляется с помощью выражения
Многомерные массивы переменного размера с числом индексов большим двух встречаются в программах довольно редко, но никаких проблем с их реализацией нет: они реализуются аналогично матрицам. Например, пусть надо реализовать трехмерный вещественный массив размера m x n x k . Захватывается линейный массив вещественных чисел размером m * n * k :
double *a; . . . a = (double *) malloc(m * n * k * sizeof(double));
Доступ к элементу с индексами x , y , z осуществляется с помощью выражения
a[(x * n + y) * k + z]
Пример: приведение матрицы к ступенчатому виду методом Гаусса
В качестве примера работы с матрицами рассмотрим алгоритм Гаусса приведения матрицы к ступенчатому виду. Метод Гаусса - один из основных результатов линейной алгебры и аналитической геометрии, к нему сводятся множество других теорем и методов линейной алгебры (теория и вычисление определителей, решение систем линейных уравнений, вычисление ранга матрицы и обратной матрицы, теория базисов конечномерных векторных пространств и т.д.).
Напомним, что матрица A с элементами a ij называется ступенчатой, если она обладает следующими двумя свойствами:
- если в матрице есть нулевая строка, то все строки ниже нее также нулевые;
- пусть a ij не равное 0 -- первый ненулевой элемент в строке с индексом i , т.е. элементы a il = 0 при l < j . Тогда все элементы в j -м столбце ниже элемента a ij равны нулю, и все элементы левее и ниже a ij также равны нулю: a kl = 0 при k > i и l =< j .
Ступенчатая матрица выглядит примерно так:
здесь тёмными квадратиками отмечены первые ненулевые элементы строк матрицы. Белым цветом изображаются нулевые элементы, серым цветом - произвольные элементы.
Алгоритм Гаусса использует элементарные преобразования матрицы двух типов.
- Преобразование первого рода :две строки матрицы меняются местами, и при этом знаки всех элементов одной из строк изменяются на противоположные.
- Преобразование второго рода :к одной строке матрицы прибавляется другая строка, умноженная на произвольное число.
Элементарные преобразования сохраняют определитель и ранг матрицы, а также множество решений линейной системы. Алгоритм Гаусса приводит произвольную матрицу элементарными преобразованиями к ступенчатому виду. Для ступенчатой квадратной матрицы определитель равен произведению диагональных элементов, а ранг - числу ненулевых строк (рангом по определению называется размерность линейной оболочки строк матрицы).
Метод Гаусса в математическом варианте состоит в следующем:
- ищем сначала ненулевой элемент в первом столбце. Если все элементы первого столбца нулевые, то переходим ко второму столбцу, и так далее. Если нашли ненулевой элемент в k -й строке, то при помощи элементарного преобразования первого рода меняем местами первую и k -ю строки, добиваясь того, чтобы первый элемент первой строки был отличен от нуля;
- используя элементарные преобразования второго рода, обнуляем все элементы первого столбца, начиная со второго элемента. Для этого от строки с номером k вычитаем первую строку, умноженную на коэффициент a k1 /a 11 .
- переходим ко второму столбцу (или j -му, если все элементы первого столбца были нулевыми), и в дальнейшем рассматриваем только часть матрицы, начиная со второй строки и ниже. Снова повторяем пункты 1) и 2) до тех пор, пока не приведем матрицу к ступенчатому виду.
Программистский вариант метода Гаусса имеет три отличия от математического:
r = -a kj /a ij . a k = a k + r * a i
Такая схема работает нормально только тогда, когда коэффициент r по абсолютной величине не превосходит единицы. В противном случае, ошибки округления умножаются на большой коэффициент и, таким образом, экспоненциально растут. Математики называют это явление неустойчивостью вычислительной схемы. Если вычислительная схема неустойчива, то полученные с ее помощью результаты не имеют никакого отношения к исходной задаче. В нашем случае схема устойчива, когда коэффициент r = -a kj /a ij не превосходит по модулю единицы. Для этого должно выполняться неравенство
|a ij | >= |a kj | при k > i
Отсюда следует, что при поиске разрешающего элемента в j -м столбце необходимо найти не первый попавшийся ненулевой элемент, а максимальный по абсолютной величине . Если он по модулю не превосходит , то считаем, что все элементы столбца нулевые; иначе меняем местами строки, ставя его на вершину столбца, и затем обнуляем столбец элементарными преобразованиями второго рода.
Ниже дан полный текст программы на Си , приводящей вещественную матрицу к ступенчатому виду. Функция , реализующая метод Гаусса , одновременно подсчитывает и ранг матрицы. Программа вводит размеры матрицы и ее элементы с клавиатуры и вызывает функцию приведения к ступенчатому виду. Затем программа печатает ступенчатый вид матрицы и ее ранг . В случае квадратной матрицы также вычисляется и печатается определитель матрицы , равный произведению диагональных элементов ступенчатой матрицы.
При реализации метода Гаусса используется схема построения цикла с помощью инварианта, см. раздел 1.5.2. В цикле меняются две переменные -- индекс строки i , 0 =< i < m - 1 , и индекс столбца j , 0 =< j < n - 1 . Инвариантом цикла является утверждение о том, что часть матрицы (математики говорят минор ) в столбцах 0,1,...j - 1 приведена к ступенчатому виду и что первый ненулевой элемент в строке i - 1 стоит в столбце с индексом меньшим j . В теле цикла рассматривается только минор матрицы в строках i,...,m - 1 и столбцах j,...,n - 1 . Сначала ищется максимальный по модулю элемент в j -м столбце. Если он по абсолютной величине не превосходит то j увеличивается на единицу (считается, что столбец нулевой). Иначе перестановкой строк разрешающий элемент ставится на вершину j -го столбца минора, и затем столбец обнуляется элементарными преобразованиями второго рода. После этого оба индекса i и j увеличиваются на единицу. Алгоритм завершается, когда либо i = m , либо j = n . По окончании алгоритма значение переменной i равно числу ненулевых строк ступенчатой матрицы, т.е. рангу исходной матрицы.
Для вычисления абсолютной величины вещественного числа x типа double мы пользуемся стандарной математической функцией fabs (x) , описанной в стандартном заголовочном файле " math .h.
#include
Загрузка...