Осуществление алгоритмов линейной алгебры

Казахстанский филиал Московского Государственного Университета

имени М.В. Ломоносова

 

 

 

 

 

 

 

Вычислительная  практика

Студент ММ-11

Шамбасов  Мухаммед

Кафедра математики и информатики

 

 

 

 

 

 

 

Астана 2013

 

Содержание

Введение

1.Основные алгоритмы

1.1.История алгоритмов

1.2.Реализация на С

1.3.Быстродействие

Вывод.

 

Введение

  Моя отчетная работа посвящена вычислительной практике. Цель данной практики – ознакомление студентов с базовыми понятиями вычислительной математики, проверка полученных знаний студентом за 1 курс в области программирования, а также умение использовать их на практике. В течении двух недель мы изучали алгоритмы: их историю, оптимизация алгоритма и так далее. По-моему учебная практика была достаточно сложной, но в то же время полезной, так как мы наглядно увидели применение языка Си в математических задачах, а точнее в линейной алгебре.

Далее в своем отчете я  опишу некоторые алгоритмы, их историю, применение, также, по возможности, быстродействие некоторых алгоритмов, к примеру  «Умножение матриц». А также программы, реализующие эти алгоритмы.

 

 

 

1.Умножение матриц.

На первом занятии учебной  практики мы познакомились с задачей  о произведении заданных двух матриц. Перейдем к определению.

Умножение матриц — одна из основных операций над матрицами. Матрица, получаемая в результате операции умножения, называется произведением матриц.

Определение:

Пусть даны две прямоугольные  матрицы  и  размерности  и  соответственно: матрицы   и   размерности   и   соответственно:

          Тогда их произведение будет называться матрица С,

 имеющая следующий  вид :

 

 

где,

Для начала мы написали обычное  умножение двух заданных матриц и  сложность такого вычисления обычного алгоритма произведения матриц по определению составляет   Θ(n³). 

Вот код программы на СИ, реализующий перемножение матриц обычным  методом.

1. #include<stdio.h>
2. #include<stdlib.h>
3. #include<time.h>
4. #define n 3
5. int main()
6. {       int i,k,j;
7.  float a[n][n],b[n][n],c[n][n];
8.  FILE *fp;
9.  fp=fopen("textf", "r");
10.  while (!feof(fp))
11.  {
12.   for (i=0; i<n; i++)
13.    for (j=0; j<n; j++)
14.     fscanf (fp, "%f", &a[i][j]);
15.  
16.   for (i=0; i<n; i++)
17.    for (j=0; j<n; j++)
18.     fscanf (fp, "%f", &b[i][j]);
19.  }
20.  for (i=0; i<n; i++)
21.   for (j=0; j<n; j++)
22.   c[i][j]=0;
23.  for (k=0; k<n; k++)
24.   for (i=0; i<n; i++)
25.    for (j=0; j<n; j++)
26.     c[i][j]=c[i][j]+a[i][k]*b[k][j];
27.  for (i=0; i<n; i++)
28.   for (j=0; j<n; j++)
29.    if ((c[i][j]>-0.00001)&&(c[i][j]<0.00001)) c[i][j]=0;
30.  for (i=0; i<n; i++){
31.   for (j=0; j<n; j++)
32.    printf ("%f ", c[i][j]);
33.   printf ("\n");
34.  }
35.  fclose(fp);
36. }

Описание: 1-4  - вход в стандартные библиотеки.

6-7 – описание  переменных.

8-9 – открытие  файла textf, в котором записаны матрицы А и В,  для чтения. 

10-19 – чтение  матриц с файла.

20-26 – реализация  формулы, для нахождения элементов   матрицы С.

27-29 – Эпсилон  – окрестность.

30-34 – Распечатка  матрицы С.

35-36 – Закрытие  файла и конец программы.

2.Метод Штрассена

История алгоритма.

Алгоритм Штрассена предназначен для быстрого Умножения матриц. Он был разработан Фолькером Штрассеномв 1969 году и является обобщением метода умножения Карацубы на матрицы.

В отличие от традиционного  алгоритма умножения матриц (по формуле c_ik = Σa_ijb_jk), работающего за время Θ(n³) = Θ(n^log₂ ⁸), алгоритм Штрассена умножает матрицы за время Θ(n^log₂ ⁷) = Θ(n^2.81), что даёт выигрыш на больших плотных матрицах начиная, примерно, от 64×64.

Несмотря на то, что алгоритм Штрассена является не самым быстрым  из существующих алгоритмов быстрого умножения матриц, он проще программируется, поэтому именно он чаще используется на практике.

Штрассен был первым, кто  показал возможность умножения  матриц более эффективным способом, чем стандартный. После публикации его работы в 1969 начались активные поиски более быстрого алгоритма. Самым  асимптотически быстрым алгоритмом на сегодняшний день является алгоритм Копперсмита - Винограда, позволяющий  перемножать матрицы за Θ(n^2.376) операций, предложенный в 1987 году и усовершенствованный в 2011 году до уровня Θ(n^2.3727). Этот алгоритм не представляет практического интереса в силу астрономически большой константы в оценке арифметической сложности. Вопрос о предельной в асимптотике скорости перемножения матриц не решен. Существует гипотеза Штрассена о том, что для достаточно больших n существует алгоритм перемножения двух матриц размера n х n за Θ(n^2+ε) операций, где ε сколь угодно малое наперед заданное положительное число. Эта гипотеза имеет сугубо теоретический интерес, так как размер матриц для которых ε действительно мало по всей видимости очень велик.

Вопрос о построении наиболее быстрого и устойчивого практического алгоритма умножения больших матриц также остается нерешенным.

 

Алгоритм

Пусть A, B — две квадратные матрицы над кольцом R. Матрица C получается по формуле:

Если размер умножаемых матриц n не является натуральной степенью двойки, мы дополняем исходные матрицы дополнительными нулевыми строками и столбцами. При этом мы получаем удобные для рекурсивного умножения размеры, но теряем в эффективности за счёт дополнительных ненужных умножений.

Разделим матрицы A, B и C на равные по размеру блочные матрицы

 

где

тогда

Однако с помощью этой процедуры нам не удалось уменьшить  количество умножений. Как и в  обычном методе, нам требуется 8 умножений.

Теперь определим новые  элементы

которые затем используются для выражения Ci, j. Таким образом, нам нужно всего 7 умножений на каждом этапе рекурсии. Элементы матрицы C выражаются из Pk по формулам 

Итерационный процесс  продолжается n раз, до тех пор пока матрицы Ci, j не выродятся в числа (элементы кольца R). На практике итерации останавливают при размере матриц от 32 до 128 и далее используют обычный метод умножения матриц. Это делают из-за того, что алгоритм Штрассена теряет эффективность по сравнению с обычным на малых матрицах в силу большего числа сложений.

 

 

Код на Си

Весь смысл метода Штрассена  содержится в функции strassen.

int strassen(int** A,int** B,int** C,int N)

{

int half=N/2, i, j, full=N;

if (N<=512)

{

multi(A,B,C,N);

}

else

{

int full = half*2,N=full;

 

int** A11; int** A12; int** A21; int** A22;

int** B11; int** B12; int** B21; int** B22;

int** C11; int** C12; int** C21; int** C22;

int** P1; int** P2; int** P3; int** P4; int** P5; int** P6; int** P7;

int** Result1; int** Result2;

int i,j;

 

A11=(int**)malloc(half*sizeof(int*));A12=(int**)malloc(half*sizeof(int*));

A21=(int**)malloc(half*sizeof(int*));A22=(int**)malloc(half*sizeof(int*));

B11=(int**)malloc(half*sizeof(int*));B12=(int**)malloc(half*sizeof(int*));

B21=(int**)malloc(half*sizeof(int*));B22=(int**)malloc(half*sizeof(int*));

C11=(int**)malloc(half*sizeof(int*));C12=(int**)malloc(half*sizeof(int*));

C21=(int**)malloc(half*sizeof(int*));C22=(int**)malloc(half*sizeof(int*));

P1=(int**)malloc(half*sizeof(int*));P2=(int**)malloc(half*sizeof(int*));

P3=(int**)malloc(half*sizeof(int*));P4=(int**)malloc(half*sizeof(int*));

P5=(int**)malloc(half*sizeof(int*));P6=(int**)malloc(half*sizeof(int*));

P7=(int**)malloc(half*sizeof(int*));Result1=(int**)malloc(half*sizeof(int*)); Result2=(int**)malloc(half*sizeof(int*));

 

for (i=0; i<half; i++)

{

A11[i]=(int*)malloc(half*sizeof(int));A12[i]=(int*)malloc(half*sizeof(int));

A21[i]=(int*)malloc(half*sizeof(int));A22[i]=(int*)malloc(half*sizeof(int));

B11[i]=(int*)malloc(half*sizeof(int)); B12[i]=(int*)malloc(half*sizeof(int));

B21[i]=(int*)malloc(half*sizeof(int));B22[i]=(int*)malloc(half*sizeof(int));

C11[i]=(int*)malloc(half*sizeof(int)); C12[i]=(int*)malloc(half*sizeof(int));

C21[i]=(int*)malloc(half*sizeof(int));C22[i]=(int*)malloc(half*sizeof(int));

P1[i]=(int*)malloc(half*sizeof(int));P2[i]=(int*)malloc(half*sizeof(int));

P3[i]=(int*)malloc(half*sizeof(int));P4[i]=(int*)malloc(half*sizeof(int));

P5[i]=(int*)malloc(half*sizeof(int));P6[i]=(int*)malloc(half*sizeof(int));

P7[i]=(int*)malloc(half*sizeof(int)); Result1[i]=(int*)malloc(half*sizeof(int));

Result2[i]=(int*)malloc(half*sizeof(int));

}

for (i=0;i<half;i++)

{

for(j=0;j<half;j++)

{

A11[i][j] = A[i][j]; A12[i][j] = A[i][j + half];

A21[i][j] = A[i + half][j];A22[i][j] = A[i + half][j + half];

            

B11[i][j] = B[i][j]; B12[i][j] = B[i][j + half];

B21[i][j] = B[i + half][j]; B22[i][j] = B[i + half][j + half];

}}

add(A11,A22,Result1,half);

add(B11,B22,Result2,half);

strassen(Result1,Result2,P1,half);

 

add(A21,A22,Result1,half);

strassen(Result1,B11,P2,half);

 

sub(B12,B22,Result1,half);

strassen(A11,Result1,P3,half);

 

sub(B21,B11,Result1,half);

strassen(A22,Result1,P4,half);

 

add(A11,A12,Result1,half);

strassen(Result1,B22,P5,half);

 

sub(A21,A11,Result1,half);

add(B11,B12,Result2,half);

strassen(Result1,Result2,P6,half);

 

sub(A12,A22,Result1,half);

add(B21,B22,Result2,half);

strassen(Result1,Result2,P7,half);

 

add(P3,P5,C12,half);

 

add(P2,P4,C21,half);

 

add(P1,P7,Result1,half);

sub(P4,P5,Result2,half);

add(Result1,Result2,C11,half);

 

sub(P1,P2,Result1,half);

add(P3,P6,Result2,half);

add(Result1,Result2,C22,half);

 

for (i=0;i<half;i++)

{

for(j=0;j<half;j++)

{

C[i][j] = C11[i][j]; C[i][j+half] = C12[i][j];

C[i+half][j] = C21[i][j]; C[i+half][j+half] = C22[i][j];

}}

 

for ( i = 0; i < half; i++)

{

free(A11[i]);free(A12[i]);free(A21[i]);free(A22[i]);

free(B11[i]);free(B12[i]);free(B21[i]);free(B22[i]);

free(C11[i]);free(C12[i]);free(C21[i]);free(C22[i]);

free(P1[i]);free(P2[i]);free(P3[i]);free(P4[i]);free(P5[i]); free(P6[i]);free(P7[i]);

free(Result1[i]);free(Result2[i]);

}

free(A11);free(A12);free(A21);free(A22);

free(B11);free(B12);free(B21);free(B22);

free(C11);free(C12);free(C21);free(C22);

free(P1);free(P2);free(P3);free(P4);free(P5);free(P6);free(P7);

free(Result1);free(Result2);

 

}

return 0;

}

Примечание: хотя в алгоритме Штрассена  обычное перемножение матриц начинают с 32, мы взяли с 512, что несколько  ускорило время исполнения алгоритма.

Ряд 1 – наивное перемножение; Ряд 2 – метод Штрассена

 

3. Метод Гаусса.

Метод Гаусса — классический метод решения системы линейных алгебраических уравнений (СЛАУ). Это метод последовательного исключения переменных, когда с помощью элементарных преобразований система уравнений приводится к равносильной системе треугольного вида, из которой последовательно, начиная с последних (по номеру) переменных, находятся все остальные переменные.

История. Хотя в настоящее время данный метод повсеместно называется методом Гаусса, он был известен и до К. Ф. Гаусса. Первое известное описание данного метода — в китайском трактате «Математика в девяти книгах», составленном между I в. до н.э. и II в. н. э.

Описание метода. Пусть исходная система выглядит следующим образом

Матрица А называется основной матрицей системы, b — столбцом свободных членов. Тогда, согласно свойству элементарных преобразований над строками, основную матрицу этой системы можно привести к ступенчатому виду (эти же преобразования нужно применять к столбцу свободных членов):

При этом будем считать, что базисный минор (ненулевой минор максимального  порядка) основной матрицы находится  в верхнем левом углу, то есть в него входят только коэффициенты при переменных x_j1, … x_jr. Тогда переменные x_j1, … x_jr называются главными переменными. Все остальные называются свободными. Если хотя бы одно число , где i>r , то рассматриваемая система несовместна, т.е. у неё нет ни одного решения. Пусть для любых i>r. Перенесём свободные переменные за знаки равенств и поделим каждое из уравнений системы на свой коэффициент при самом левом x ( , где i— номер строки):

Где

Если свободным  переменным системы (2) придавать все  возможные значения и решать новую  систему относительно главных неизвестных  снизу вверх (то есть от нижнего уравнения  к верхнему), то мы получим все

решения этой СЛАУ. Так как эта система получена путём элементарных преобразований над исходной системой (1), то по теореме  об эквивалентности при элементарных преобразованиях системы (1) и (2)

эквивалентны, то есть множества их решений совпадают.

Теорема Кронекера — Капелли.

Система совместна тогда и только тогда, когда ранг её основной матрицы  равен рангу её расширенной матрицы.

Вот псевдокод программы на СИ, реализующий перемножение матриц обычным методом.

1. for(s=0; s<n; s++){
2.   max=a[s][s];
3.   k=s;
4.   for (i=s; i<n; i++){
5.    if (abs(a[i][s])>max){
6.     max=abs(a[i][s]);
7.     k=i;
8.    }
9.   }
10.   if (a[k][s]<0) max=-max;
11.   if(max!=0){
12.    if (k!=s) {
13.     for (j=s; j<n; j++){
14.      c=a[s][j];
15.      a[s][j]=a[k][j]/max;
16.      a[k][j]=c; 
17.     }
18.     m=b[s];
19.     b[s]=b[k]/max;   
20.     b[k]=m;
21.     
22.    }else{
23.     for(j=s; j<n; j++){
24.      a[k][j]=a[k][j]/max; 
25.     }
26.     b[k]=b[k]/max;
27.    }
28.   } else {
29.    printf ("mn-vo reshenii\n");
30.    goto m10;
31.   }
32.   for (i=s+1; i<n; i++){
33.    c=a[i][s];
34.    for(j=s; j<n; j++)
35.     a[i][j]=a[i][j]-a[s][j]*c;
36.     b[i]=b[i]-b[s]*c;
37.   }   
38.  }
39.  for (i=0; i<n; i++){
40.   for(j=0;j<n; j++)
41.    printf("%f ", a[i][j]);
42.    printf("%f ", b[i]);
43.    printf("\n");
44.  }
45.  for (i=n-1; i>-1; i--){
46.    z=0;
47.    for(j=n-1; j>i-1; j--)
48.     z=z+a[i][j+1]*x[j+1];
49.    x[i]=b[i]-z;
50.   }
51.  for (i=0; i<n; i++)
52.   printf("%f\n",x[i]);