Алговики - Вклад [ru]

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-12T20:43:53Z

Малков Кирилл: /* Информационный граф */

== Свойства и структура алгоритмов ==
=== Общее описание алгоритма ===
Матрица <math>B = [b_{ij}]</math> называется (верхней) двухдиагональной, или бидиагональной, если <math>b_{ij} = 0</math>, <math>i > j</math> или <math>i + 1 < j</math>. Алгоритм бидиагонализации матрицы методом отражений приводит произвольную <math>n \times n</math> матрицу A к бидиагональному виду
<math>B = PAQ</math>,
где <math>P</math> и <math>Q</math> - произведения конечного числа матриц отражения(или вращения).
=== Математическое описание алгоритма ===
Пусть дана матрица <math>A</math> размера <math>n \times n</math>. Вычисляется бидиагональная матрица <math>D</math> размера <math>n \times n</math>, такая что <math>A = PDQ</math> для некоторых матриц <math>P</math> и <math>Q</math>, являющихся произведением некоторого числа матриц отражения.Сначала мы умножаем A слева на матрицу отражения, аннулирующую все поддиагональные элементы первого столбца. Затем умножаем результат справа на матрицу отражения, аннулирующую элементы первой строки в позициях с 3-й по n-ю. Далее умножением слева аннулируем все поддиагональные элементы второго столбца, затем умножением справа аннулируем получаем нули во второй строке в позициях с 4-й по n-ю и т.д. После каждого умножения на матрицу отражения все ранее полученные нули остаются.
Формулы метода:
Пусть <math>A = [a_1, a_2, ..., a_n],</math>
<math>D_1 = A, P_1 = I, Q_1 = I.</math>
 Если <math>i</math> четное, то <math>D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*)</math>
 <math>u_i = \frac{[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]}{||[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]||_2}, Q_i = (I - 2 * u_{i} * u_{i}^*) * Q_{i-1}, P_i = P_{i-1}.</math>
 Если <math>i</math> нечетное, то <math>D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1},</math>
 <math>u_i = \frac{[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]}{||[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]||_2}, Q_i = Q_{i-1}, P_i = P_{i-1} * (I - 2 * u_{i} * u_{i}^*).</math>

=== Вычислительное ядро алгоритма ===
Вычислительное ядро алгоритма составляет нахождение на каждом этапе матрицы отражения(всего <math>2n - 2</math>).

=== Макроструктура алгоритма ===
Операции скалярного произведения, умножения матрицы на вектор
=== Схема реализации последовательного алгоритма ===
<math>D_1 = A, P_1 = I, Q_1 = I.</math>
<math>1) для четного i : D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*), Q_i = (I - 2 * u_{i} * u_{i}^*) * Q_{i-1}, P_i = P_{i-1};</math>
<math> для нечетного i : D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1}, Q_i = Q_{i-1}, P_i = P_{i-1} * (I - 2 * u_{i} * u_{i}^*).</math>
<math>2) для четного i : u_i = \frac{[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]}{||[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]||_2},</math>
<math> для нечетного i : u_i = \frac{[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]}{||[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]||_2}.</math>
=== Последовательная сложность ===
Для бидиагонализации матрицы порядка n методом отражений при последовательном варианте требуется:
<math>1) 2 * n - 2</math> вычислений квадратного корня,
<math>2) 4 * n - 4 </math>делений,
<math>3) (4 * n - 4)(1 + n + n^2)</math> умножений,
<math>4) (4 * n - 4)(1 + n + n^2)</math> сложений.
Сложность последовательного алгоритма есть <math>O(n^3)</math>.

=== Информационный граф ===
Информационный граф для количества процессов равного 4.
[[Файл:Diagram.png|Информационный граф]]
 <math>I-2uu^*</math> - вычисление матрицы отражения;
 <math>D, P, Q</math> - умножение матрицы отражения на матрицы D, P, Q;
 <math>TRANS</math> - транспонирование матрицы D;

=== Ресурс параллелизма алгоритма ===
При параллельном вычислении двухдиагональной матрицы методом отражений требуется последовательно выполнить следующие ярусы:
 -- k ярусов с вычислением <math>4 * (n/k) - 4</math> делений,
 -- k ярусов с вычислением <math>(4 * (n/k) - 4)(1 + n + n^2)</math> умножений,
 -- k ярусов с вычислением <math>(4 * (n/k) - 4)(1 + n + n^2)</math> сложений.
 Получим, что параллельная сложность алгоритма есть <math>O(n^2)</math>.

=== Входные и выходные данные алгоритма ===
На вход подается матрица <math>A</math> и ее размер <math>n</math>. Выходные данные -- матрицы <math>P, Q, D,</math> такие что <math>A = PDQ </math>

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

== Программная реализация алгоритма ==
<source lang="cpp">
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <mkl.h>
#include <time.h>
#include <omp.h>

void matrixA_gen(double *A, int num, int n, int m)
{
int r1 = 0, r2 = 0;
for(int i = 0; i < n; i++)
for(int j = 0; j < m; j++)
{
//A[i * m + j] = -1.0 * num * rand() / (double) RAND_MAX * i * j + 2.0 * rand() / (double) RAND_MAX;
A[i*m + j] = i + 1 + num;
}
}

void matrixP_gen(double *P, int num, int n, int m)
{
for(int i = 0; i < n; i++)
{
for(int j = 0; j < m; j++)
{
if((i % m) == j && i >= num * m && i < (num + 1) * m)
{
P[i * m + j] = 1;
} else {
P[i * m + j] = 0;
}
}
}
}

double sc_mult(double *a, double *b, int n)
{
double sum = 0;
for(int i = 0; i < n; i++)
sum += a[i] * b[i];
return sum;
};

double norm(double *a, int n)
{
return sqrt(sc_mult(a, a, n));
};

void mult_vec(double *A, double *u, int n, int m)
{
char t = 'n';
int int_one = 1;
double *a, one = 1.0, zero = 0.0, min_two = -2.0;
a = (double *)calloc(m, sizeof(*a));
dgemv(&t, &m, &n, &one, A, &m, u, &int_one, &zero, a, &int_one);
dger(&m, &n, &min_two, a, &int_one, u, &int_one, A, &m);
};

void mult(double *A, double *B, int l, int n)
{
double *G;
G = (double *)calloc(n * n, sizeof(*G));
if(l)
{
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
G[i * n + j] += A[i * n + k] * B[k * n + j];
} else {
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
G[i * n + j] += B[i * n + k] * A[k * n + j];
}
for (int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
A[i*n + j] = G[i*n +j];
free(G);
};

void reflect(double *u, double *a, int n, int j)
{
for (int i = 0; i < n; i++)
{
if (i < j)
{
u[i] = 0;
} else {
if (i == j){
u[i] = a[i] - norm(a, n);
} else {
u[i] = a[i];
}
}
}
double sc = norm(u, n);
if(sc != 0)
{
for(int i = 0; i < n; i++)
{
u[i] /= sc;
}
}
};

void trans(double *A, int n)
{
double *B;
B = (double *) calloc(n * n, sizeof(*B));
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
B[i * n + j] = A[j * n + i];
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
A[i * n + j] = B[i * n + j];
free(B);
return;
}

int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
int size, num, error, n, m, i = 0;
double d_norm = 0;
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &num);
n = 64;
srand(time(NULL));
if(size != 0)
{
m = n/size;
}
MPI_Request req;
MPI_Status status;
double *M, *A, *P, *Q, *C, *D, *Df, *Pf, *Qf, *V, *A1, *W;
double *u, *a;
A1 = (double *) calloc(m * n, sizeof(*A1));
V = (double *)calloc(n * n, sizeof(*V));
W = (double *)calloc(n * n, sizeof(*W));
Qf = (double *)calloc(n * n, sizeof(*Qf));
M = (double *)calloc(n * n, sizeof(*M));
Pf = (double *)calloc(n * n, sizeof(*Pf));
A = (double *)calloc(n * m, sizeof(*A));
P = (double *)calloc(n * m, sizeof(*P));
Q = (double *)calloc(n * m, sizeof(*Q));
C = (double *)calloc(m * m, sizeof(*C));
D = (double *)calloc(n * n, sizeof(*D));
Df = (double *)calloc(n * n, sizeof(*Df));
u = (double *)calloc(n, sizeof(*u));
a = (double *)calloc(n, sizeof(*a));
matrixA_gen(A, num, n, m);
matrixP_gen(P, num, n, m);
matrixP_gen(Q, num, n, m);
MPI_Gather(A, n * m, MPI_DOUBLE, D, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if(num == 0)
{
for(int k = 0; k < size; k++)
for(int j = 0; j < n; j++)
for(int h = 0; h < m; h++)
{
M[h + k * m + j * n] = D[k * m * n + j * m + h];
}
}

double t1 = MPI_Wtime();

for(int f = 0; f < 2 * n - 2; f++)
{
i = f / 2;
if(i < num * m || i >= (num + 1) * m )
{
MPI_Recv(u, n, MPI_DOUBLE, (i / m), 0, MPI_COMM_WORLD, &status);
mult_vec(A, u, n, m);
if(!(f%2))
mult_vec(P, u, n, m);
else
mult_vec(Q, u, n, m);
} else {
if(!(f%2))
{
for(int k = 0; k < n; k++)
{
if(k < i)
{
a[k] = 0;
} else {
a[k] = A[(i % m) + k * m];
}
}
reflect(u, a, n, i);
for(int l = 0; l < size; l++)
{
if(l != num)
{
MPI_Isend(u, n, MPI_DOUBLE, l, 0, MPI_COMM_WORLD, &req);
}
}
} else {
for(int k = 0; k < n; k++)
{
if(k < i+1)
{
a[k] = 0;
} else {
a[k] = A[(i % m) + k * m];
}
}
reflect(u, a, n, i+1);
for(int l = 0; l < size; l++)
{
if(l != num)
{
MPI_Isend(u, n, MPI_DOUBLE, l, 0, MPI_COMM_WORLD, &req);
}
}
}
mult_vec(A, u, n, m);
if(!(f%2))
mult_vec(P, u, n, m);
else
mult_vec(Q, u, n, m);
}
for(int k = 0; k < n; k++)
for(int j = 0; j < m; j++)
A1[k * m + j] = A[k * m + j];

for(int j = 0; j < size; j++)
{
if(j != num)
{
for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
C[h * m + l] = A1[j * m * m + h + l * m];
MPI_Isend(C, m * m, MPI_DOUBLE, j, 0, MPI_COMM_WORLD, &req);
} else {
for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
C[h * m + l] = A1[j * m * m + h + l * m];

for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
A1[j * m * m + h + l * m] = C[h + l * m];
}
}

for(int k = 0; k < size; k++)
{
if(k != num)
{
MPI_Recv(C, m * m, MPI_DOUBLE, k, 0, MPI_COMM_WORLD, &status);

for(int l = 0; l < m; l++)
{
for(int h = 0; h < m; h++)
{
A1[k * m * m + l + h * m] = C[h * m + l];
}
}
}
}
for(int k = 0; k < n; k++)
for(int j = 0; j < m; j++)
A[k * m + j] = A1[k * m + j];
}
double t2 = MPI_Wtime();
MPI_Gather(P, n * m, MPI_DOUBLE, D, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(Q, n * m, MPI_DOUBLE, W, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(A, n * m, MPI_DOUBLE, V, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
return 0;
}
</source>

[[Файл:Безымянный777.png|Результаты тестирования алгоритма на кластере ИВМ РАН]]
=== Особенности реализации последовательного алгоритма ===
Для оптимальности работы алгоритма следует реализовать функцию умножения матрицы отражения матрицы отражения на произвольную матрицу, это заметно сократит время время его работы.
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
Алгоритм бидиагонализации матрицы довольно последователен: на каждом шаге считается матрица отражения, и распараллелить прямым образом алгоритм нельзя, поскольку каждый раз матрица отражения строится к новой матрице. Однако можно разбить Матрицу <math> D</math> на блоки, и каждый процесс будет занят умножением матрицы отражения на соответствующий блок.

=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
Параллельная реализация алгоритма наиболее эффективна на машинах с быстрым обменом данных меду процессами
=== Существующие реализации алгоритма ===

== Литература ==
Е. Е. Тыртышников, "Методы численного анализа".

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-12T20:42:08Z

Малков Кирилл: /* Программная реализация алгоритма */

== Свойства и структура алгоритмов ==
=== Общее описание алгоритма ===
Матрица <math>B = [b_{ij}]</math> называется (верхней) двухдиагональной, или бидиагональной, если <math>b_{ij} = 0</math>, <math>i > j</math> или <math>i + 1 < j</math>. Алгоритм бидиагонализации матрицы методом отражений приводит произвольную <math>n \times n</math> матрицу A к бидиагональному виду
<math>B = PAQ</math>,
где <math>P</math> и <math>Q</math> - произведения конечного числа матриц отражения(или вращения).
=== Математическое описание алгоритма ===
Пусть дана матрица <math>A</math> размера <math>n \times n</math>. Вычисляется бидиагональная матрица <math>D</math> размера <math>n \times n</math>, такая что <math>A = PDQ</math> для некоторых матриц <math>P</math> и <math>Q</math>, являющихся произведением некоторого числа матриц отражения.Сначала мы умножаем A слева на матрицу отражения, аннулирующую все поддиагональные элементы первого столбца. Затем умножаем результат справа на матрицу отражения, аннулирующую элементы первой строки в позициях с 3-й по n-ю. Далее умножением слева аннулируем все поддиагональные элементы второго столбца, затем умножением справа аннулируем получаем нули во второй строке в позициях с 4-й по n-ю и т.д. После каждого умножения на матрицу отражения все ранее полученные нули остаются.
Формулы метода:
Пусть <math>A = [a_1, a_2, ..., a_n],</math>
<math>D_1 = A, P_1 = I, Q_1 = I.</math>
 Если <math>i</math> четное, то <math>D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*)</math>
 <math>u_i = \frac{[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]}{||[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]||_2}, Q_i = (I - 2 * u_{i} * u_{i}^*) * Q_{i-1}, P_i = P_{i-1}.</math>
 Если <math>i</math> нечетное, то <math>D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1},</math>
 <math>u_i = \frac{[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]}{||[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]||_2}, Q_i = Q_{i-1}, P_i = P_{i-1} * (I - 2 * u_{i} * u_{i}^*).</math>

=== Вычислительное ядро алгоритма ===
Вычислительное ядро алгоритма составляет нахождение на каждом этапе матрицы отражения(всего <math>2n - 2</math>).

=== Макроструктура алгоритма ===
Операции скалярного произведения, умножения матрицы на вектор
=== Схема реализации последовательного алгоритма ===
<math>D_1 = A, P_1 = I, Q_1 = I.</math>
<math>1) для четного i : D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*), Q_i = (I - 2 * u_{i} * u_{i}^*) * Q_{i-1}, P_i = P_{i-1};</math>
<math> для нечетного i : D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1}, Q_i = Q_{i-1}, P_i = P_{i-1} * (I - 2 * u_{i} * u_{i}^*).</math>
<math>2) для четного i : u_i = \frac{[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]}{||[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]||_2},</math>
<math> для нечетного i : u_i = \frac{[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]}{||[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]||_2}.</math>
=== Последовательная сложность ===
Для бидиагонализации матрицы порядка n методом отражений при последовательном варианте требуется:
<math>1) 2 * n - 2</math> вычислений квадратного корня,
<math>2) 4 * n - 4 </math>делений,
<math>3) (4 * n - 4)(1 + n + n^2)</math> умножений,
<math>4) (4 * n - 4)(1 + n + n^2)</math> сложений.
Сложность последовательного алгоритма есть <math>O(n^3)</math>.

=== Информационный граф ===
[[Файл:Diagram.png|Информационный граф]]
 <math>I-2uu^*</math> - вычисление матрицы отражения;
 <math>D, P, Q</math> - умножение матрицы отражения на матрицы D, P, Q;
 <math>TRANS</math> - транспонирование матрицы D;

=== Ресурс параллелизма алгоритма ===
При параллельном вычислении двухдиагональной матрицы методом отражений требуется последовательно выполнить следующие ярусы:
 -- k ярусов с вычислением <math>4 * (n/k) - 4</math> делений,
 -- k ярусов с вычислением <math>(4 * (n/k) - 4)(1 + n + n^2)</math> умножений,
 -- k ярусов с вычислением <math>(4 * (n/k) - 4)(1 + n + n^2)</math> сложений.
 Получим, что параллельная сложность алгоритма есть <math>O(n^2)</math>.

=== Входные и выходные данные алгоритма ===
На вход подается матрица <math>A</math> и ее размер <math>n</math>. Выходные данные -- матрицы <math>P, Q, D,</math> такие что <math>A = PDQ </math>

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

== Программная реализация алгоритма ==
<source lang="cpp">
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <mkl.h>
#include <time.h>
#include <omp.h>

void matrixA_gen(double *A, int num, int n, int m)
{
int r1 = 0, r2 = 0;
for(int i = 0; i < n; i++)
for(int j = 0; j < m; j++)
{
//A[i * m + j] = -1.0 * num * rand() / (double) RAND_MAX * i * j + 2.0 * rand() / (double) RAND_MAX;
A[i*m + j] = i + 1 + num;
}
}

void matrixP_gen(double *P, int num, int n, int m)
{
for(int i = 0; i < n; i++)
{
for(int j = 0; j < m; j++)
{
if((i % m) == j && i >= num * m && i < (num + 1) * m)
{
P[i * m + j] = 1;
} else {
P[i * m + j] = 0;
}
}
}
}

double sc_mult(double *a, double *b, int n)
{
double sum = 0;
for(int i = 0; i < n; i++)
sum += a[i] * b[i];
return sum;
};

double norm(double *a, int n)
{
return sqrt(sc_mult(a, a, n));
};

void mult_vec(double *A, double *u, int n, int m)
{
char t = 'n';
int int_one = 1;
double *a, one = 1.0, zero = 0.0, min_two = -2.0;
a = (double *)calloc(m, sizeof(*a));
dgemv(&t, &m, &n, &one, A, &m, u, &int_one, &zero, a, &int_one);
dger(&m, &n, &min_two, a, &int_one, u, &int_one, A, &m);
};

void mult(double *A, double *B, int l, int n)
{
double *G;
G = (double *)calloc(n * n, sizeof(*G));
if(l)
{
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
G[i * n + j] += A[i * n + k] * B[k * n + j];
} else {
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
G[i * n + j] += B[i * n + k] * A[k * n + j];
}
for (int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
A[i*n + j] = G[i*n +j];
free(G);
};

void reflect(double *u, double *a, int n, int j)
{
for (int i = 0; i < n; i++)
{
if (i < j)
{
u[i] = 0;
} else {
if (i == j){
u[i] = a[i] - norm(a, n);
} else {
u[i] = a[i];
}
}
}
double sc = norm(u, n);
if(sc != 0)
{
for(int i = 0; i < n; i++)
{
u[i] /= sc;
}
}
};

void trans(double *A, int n)
{
double *B;
B = (double *) calloc(n * n, sizeof(*B));
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
B[i * n + j] = A[j * n + i];
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
A[i * n + j] = B[i * n + j];
free(B);
return;
}

int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
int size, num, error, n, m, i = 0;
double d_norm = 0;
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &num);
n = 64;
srand(time(NULL));
if(size != 0)
{
m = n/size;
}
MPI_Request req;
MPI_Status status;
double *M, *A, *P, *Q, *C, *D, *Df, *Pf, *Qf, *V, *A1, *W;
double *u, *a;
A1 = (double *) calloc(m * n, sizeof(*A1));
V = (double *)calloc(n * n, sizeof(*V));
W = (double *)calloc(n * n, sizeof(*W));
Qf = (double *)calloc(n * n, sizeof(*Qf));
M = (double *)calloc(n * n, sizeof(*M));
Pf = (double *)calloc(n * n, sizeof(*Pf));
A = (double *)calloc(n * m, sizeof(*A));
P = (double *)calloc(n * m, sizeof(*P));
Q = (double *)calloc(n * m, sizeof(*Q));
C = (double *)calloc(m * m, sizeof(*C));
D = (double *)calloc(n * n, sizeof(*D));
Df = (double *)calloc(n * n, sizeof(*Df));
u = (double *)calloc(n, sizeof(*u));
a = (double *)calloc(n, sizeof(*a));
matrixA_gen(A, num, n, m);
matrixP_gen(P, num, n, m);
matrixP_gen(Q, num, n, m);
MPI_Gather(A, n * m, MPI_DOUBLE, D, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if(num == 0)
{
for(int k = 0; k < size; k++)
for(int j = 0; j < n; j++)
for(int h = 0; h < m; h++)
{
M[h + k * m + j * n] = D[k * m * n + j * m + h];
}
}

double t1 = MPI_Wtime();

for(int f = 0; f < 2 * n - 2; f++)
{
i = f / 2;
if(i < num * m || i >= (num + 1) * m )
{
MPI_Recv(u, n, MPI_DOUBLE, (i / m), 0, MPI_COMM_WORLD, &status);
mult_vec(A, u, n, m);
if(!(f%2))
mult_vec(P, u, n, m);
else
mult_vec(Q, u, n, m);
} else {
if(!(f%2))
{
for(int k = 0; k < n; k++)
{
if(k < i)
{
a[k] = 0;
} else {
a[k] = A[(i % m) + k * m];
}
}
reflect(u, a, n, i);
for(int l = 0; l < size; l++)
{
if(l != num)
{
MPI_Isend(u, n, MPI_DOUBLE, l, 0, MPI_COMM_WORLD, &req);
}
}
} else {
for(int k = 0; k < n; k++)
{
if(k < i+1)
{
a[k] = 0;
} else {
a[k] = A[(i % m) + k * m];
}
}
reflect(u, a, n, i+1);
for(int l = 0; l < size; l++)
{
if(l != num)
{
MPI_Isend(u, n, MPI_DOUBLE, l, 0, MPI_COMM_WORLD, &req);
}
}
}
mult_vec(A, u, n, m);
if(!(f%2))
mult_vec(P, u, n, m);
else
mult_vec(Q, u, n, m);
}
for(int k = 0; k < n; k++)
for(int j = 0; j < m; j++)
A1[k * m + j] = A[k * m + j];

for(int j = 0; j < size; j++)
{
if(j != num)
{
for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
C[h * m + l] = A1[j * m * m + h + l * m];
MPI_Isend(C, m * m, MPI_DOUBLE, j, 0, MPI_COMM_WORLD, &req);
} else {
for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
C[h * m + l] = A1[j * m * m + h + l * m];

for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
A1[j * m * m + h + l * m] = C[h + l * m];
}
}

for(int k = 0; k < size; k++)
{
if(k != num)
{
MPI_Recv(C, m * m, MPI_DOUBLE, k, 0, MPI_COMM_WORLD, &status);

for(int l = 0; l < m; l++)
{
for(int h = 0; h < m; h++)
{
A1[k * m * m + l + h * m] = C[h * m + l];
}
}
}
}
for(int k = 0; k < n; k++)
for(int j = 0; j < m; j++)
A[k * m + j] = A1[k * m + j];
}
double t2 = MPI_Wtime();
MPI_Gather(P, n * m, MPI_DOUBLE, D, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(Q, n * m, MPI_DOUBLE, W, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(A, n * m, MPI_DOUBLE, V, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
return 0;
}
</source>

[[Файл:Безымянный777.png|Результаты тестирования алгоритма на кластере ИВМ РАН]]
=== Особенности реализации последовательного алгоритма ===
Для оптимальности работы алгоритма следует реализовать функцию умножения матрицы отражения матрицы отражения на произвольную матрицу, это заметно сократит время время его работы.
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
Алгоритм бидиагонализации матрицы довольно последователен: на каждом шаге считается матрица отражения, и распараллелить прямым образом алгоритм нельзя, поскольку каждый раз матрица отражения строится к новой матрице. Однако можно разбить Матрицу <math> D</math> на блоки, и каждый процесс будет занят умножением матрицы отражения на соответствующий блок.

=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
Параллельная реализация алгоритма наиболее эффективна на машинах с быстрым обменом данных меду процессами
=== Существующие реализации алгоритма ===

== Литература ==
Е. Е. Тыртышников, "Методы численного анализа".

Файл:Безымянный777.png

2021-12-12T20:41:09Z

Малков Кирилл:

Результаты тестирования алгоритма на кластере ИВМ РАН

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-12T20:30:58Z

Малков Кирилл: /* Информационный граф */

== Свойства и структура алгоритмов ==
=== Общее описание алгоритма ===
Матрица <math>B = [b_{ij}]</math> называется (верхней) двухдиагональной, или бидиагональной, если <math>b_{ij} = 0</math>, <math>i > j</math> или <math>i + 1 < j</math>. Алгоритм бидиагонализации матрицы методом отражений приводит произвольную <math>n \times n</math> матрицу A к бидиагональному виду
<math>B = PAQ</math>,
где <math>P</math> и <math>Q</math> - произведения конечного числа матриц отражения(или вращения).
=== Математическое описание алгоритма ===
Пусть дана матрица <math>A</math> размера <math>n \times n</math>. Вычисляется бидиагональная матрица <math>D</math> размера <math>n \times n</math>, такая что <math>A = PDQ</math> для некоторых матриц <math>P</math> и <math>Q</math>, являющихся произведением некоторого числа матриц отражения.Сначала мы умножаем A слева на матрицу отражения, аннулирующую все поддиагональные элементы первого столбца. Затем умножаем результат справа на матрицу отражения, аннулирующую элементы первой строки в позициях с 3-й по n-ю. Далее умножением слева аннулируем все поддиагональные элементы второго столбца, затем умножением справа аннулируем получаем нули во второй строке в позициях с 4-й по n-ю и т.д. После каждого умножения на матрицу отражения все ранее полученные нули остаются.
Формулы метода:
Пусть <math>A = [a_1, a_2, ..., a_n],</math>
<math>D_1 = A, P_1 = I, Q_1 = I.</math>
 Если <math>i</math> четное, то <math>D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*)</math>
 <math>u_i = \frac{[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]}{||[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]||_2}, Q_i = (I - 2 * u_{i} * u_{i}^*) * Q_{i-1}, P_i = P_{i-1}.</math>
 Если <math>i</math> нечетное, то <math>D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1},</math>
 <math>u_i = \frac{[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]}{||[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]||_2}, Q_i = Q_{i-1}, P_i = P_{i-1} * (I - 2 * u_{i} * u_{i}^*).</math>

=== Вычислительное ядро алгоритма ===
Вычислительное ядро алгоритма составляет нахождение на каждом этапе матрицы отражения(всего <math>2n - 2</math>).

=== Макроструктура алгоритма ===
Операции скалярного произведения, умножения матрицы на вектор
=== Схема реализации последовательного алгоритма ===
<math>D_1 = A, P_1 = I, Q_1 = I.</math>
<math>1) для четного i : D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*), Q_i = (I - 2 * u_{i} * u_{i}^*) * Q_{i-1}, P_i = P_{i-1};</math>
<math> для нечетного i : D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1}, Q_i = Q_{i-1}, P_i = P_{i-1} * (I - 2 * u_{i} * u_{i}^*).</math>
<math>2) для четного i : u_i = \frac{[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]}{||[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]||_2},</math>
<math> для нечетного i : u_i = \frac{[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]}{||[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]||_2}.</math>
=== Последовательная сложность ===
Для бидиагонализации матрицы порядка n методом отражений при последовательном варианте требуется:
<math>1) 2 * n - 2</math> вычислений квадратного корня,
<math>2) 4 * n - 4 </math>делений,
<math>3) (4 * n - 4)(1 + n + n^2)</math> умножений,
<math>4) (4 * n - 4)(1 + n + n^2)</math> сложений.
Сложность последовательного алгоритма есть <math>O(n^3)</math>.

=== Информационный граф ===
[[Файл:Diagram.png|Информационный граф]]
 <math>I-2uu^*</math> - вычисление матрицы отражения;
 <math>D, P, Q</math> - умножение матрицы отражения на матрицы D, P, Q;
 <math>TRANS</math> - транспонирование матрицы D;

=== Ресурс параллелизма алгоритма ===
При параллельном вычислении двухдиагональной матрицы методом отражений требуется последовательно выполнить следующие ярусы:
 -- k ярусов с вычислением <math>4 * (n/k) - 4</math> делений,
 -- k ярусов с вычислением <math>(4 * (n/k) - 4)(1 + n + n^2)</math> умножений,
 -- k ярусов с вычислением <math>(4 * (n/k) - 4)(1 + n + n^2)</math> сложений.
 Получим, что параллельная сложность алгоритма есть <math>O(n^2)</math>.

=== Входные и выходные данные алгоритма ===
На вход подается матрица <math>A</math> и ее размер <math>n</math>. Выходные данные -- матрицы <math>P, Q, D,</math> такие что <math>A = PDQ </math>

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

== Программная реализация алгоритма ==
<source lang="cpp">
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <mkl.h>
#include <time.h>
#include <omp.h>

void matrixA_gen(double *A, int num, int n, int m)
{
int r1 = 0, r2 = 0;
for(int i = 0; i < n; i++)
for(int j = 0; j < m; j++)
{
//A[i * m + j] = -1.0 * num * rand() / (double) RAND_MAX * i * j + 2.0 * rand() / (double) RAND_MAX;
A[i*m + j] = i + 1 + num;
}
}

void matrixP_gen(double *P, int num, int n, int m)
{
for(int i = 0; i < n; i++)
{
for(int j = 0; j < m; j++)
{
if((i % m) == j && i >= num * m && i < (num + 1) * m)
{
P[i * m + j] = 1;
} else {
P[i * m + j] = 0;
}
}
}
}

double sc_mult(double *a, double *b, int n)
{
double sum = 0;
for(int i = 0; i < n; i++)
sum += a[i] * b[i];
return sum;
};

double norm(double *a, int n)
{
return sqrt(sc_mult(a, a, n));
};

void mult_vec(double *A, double *u, int n, int m)
{
char t = 'n';
int int_one = 1;
double *a, one = 1.0, zero = 0.0, min_two = -2.0;
a = (double *)calloc(m, sizeof(*a));
dgemv(&t, &m, &n, &one, A, &m, u, &int_one, &zero, a, &int_one);
dger(&m, &n, &min_two, a, &int_one, u, &int_one, A, &m);
};

void mult(double *A, double *B, int l, int n)
{
double *G;
G = (double *)calloc(n * n, sizeof(*G));
if(l)
{
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
G[i * n + j] += A[i * n + k] * B[k * n + j];
} else {
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
G[i * n + j] += B[i * n + k] * A[k * n + j];
}
for (int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
A[i*n + j] = G[i*n +j];
free(G);
};

void reflect(double *u, double *a, int n, int j)
{
for (int i = 0; i < n; i++)
{
if (i < j)
{
u[i] = 0;
} else {
if (i == j){
u[i] = a[i] - norm(a, n);
} else {
u[i] = a[i];
}
}
}
double sc = norm(u, n);
if(sc != 0)
{
for(int i = 0; i < n; i++)
{
u[i] /= sc;
}
}
};

void trans(double *A, int n)
{
double *B;
B = (double *) calloc(n * n, sizeof(*B));
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
B[i * n + j] = A[j * n + i];
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
A[i * n + j] = B[i * n + j];
free(B);
return;
}

int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
int size, num, error, n, m, i = 0;
double d_norm = 0;
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &num);
n = 64;
srand(time(NULL));
if(size != 0)
{
m = n/size;
}
MPI_Request req;
MPI_Status status;
double *M, *A, *P, *Q, *C, *D, *Df, *Pf, *Qf, *V, *A1, *W;
double *u, *a;
A1 = (double *) calloc(m * n, sizeof(*A1));
V = (double *)calloc(n * n, sizeof(*V));
W = (double *)calloc(n * n, sizeof(*W));
Qf = (double *)calloc(n * n, sizeof(*Qf));
M = (double *)calloc(n * n, sizeof(*M));
Pf = (double *)calloc(n * n, sizeof(*Pf));
A = (double *)calloc(n * m, sizeof(*A));
P = (double *)calloc(n * m, sizeof(*P));
Q = (double *)calloc(n * m, sizeof(*Q));
C = (double *)calloc(m * m, sizeof(*C));
D = (double *)calloc(n * n, sizeof(*D));
Df = (double *)calloc(n * n, sizeof(*Df));
u = (double *)calloc(n, sizeof(*u));
a = (double *)calloc(n, sizeof(*a));
matrixA_gen(A, num, n, m);
matrixP_gen(P, num, n, m);
matrixP_gen(Q, num, n, m);
MPI_Gather(A, n * m, MPI_DOUBLE, D, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if(num == 0)
{
for(int k = 0; k < size; k++)
for(int j = 0; j < n; j++)
for(int h = 0; h < m; h++)
{
M[h + k * m + j * n] = D[k * m * n + j * m + h];
}
}

double t1 = MPI_Wtime();

for(int f = 0; f < 2 * n - 2; f++)
{
i = f / 2;
if(i < num * m || i >= (num + 1) * m )
{
MPI_Recv(u, n, MPI_DOUBLE, (i / m), 0, MPI_COMM_WORLD, &status);
mult_vec(A, u, n, m);
if(!(f%2))
mult_vec(P, u, n, m);
else
mult_vec(Q, u, n, m);
} else {
if(!(f%2))
{
for(int k = 0; k < n; k++)
{
if(k < i)
{
a[k] = 0;
} else {
a[k] = A[(i % m) + k * m];
}
}
reflect(u, a, n, i);
for(int l = 0; l < size; l++)
{
if(l != num)
{
MPI_Isend(u, n, MPI_DOUBLE, l, 0, MPI_COMM_WORLD, &req);
}
}
} else {
for(int k = 0; k < n; k++)
{
if(k < i+1)
{
a[k] = 0;
} else {
a[k] = A[(i % m) + k * m];
}
}
reflect(u, a, n, i+1);
for(int l = 0; l < size; l++)
{
if(l != num)
{
MPI_Isend(u, n, MPI_DOUBLE, l, 0, MPI_COMM_WORLD, &req);
}
}
}
mult_vec(A, u, n, m);
if(!(f%2))
mult_vec(P, u, n, m);
else
mult_vec(Q, u, n, m);
}
for(int k = 0; k < n; k++)
for(int j = 0; j < m; j++)
A1[k * m + j] = A[k * m + j];

for(int j = 0; j < size; j++)
{
if(j != num)
{
for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
C[h * m + l] = A1[j * m * m + h + l * m];
MPI_Isend(C, m * m, MPI_DOUBLE, j, 0, MPI_COMM_WORLD, &req);
} else {
for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
C[h * m + l] = A1[j * m * m + h + l * m];

for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
A1[j * m * m + h + l * m] = C[h + l * m];
}
}

for(int k = 0; k < size; k++)
{
if(k != num)
{
MPI_Recv(C, m * m, MPI_DOUBLE, k, 0, MPI_COMM_WORLD, &status);

for(int l = 0; l < m; l++)
{
for(int h = 0; h < m; h++)
{
A1[k * m * m + l + h * m] = C[h * m + l];
}
}
}
}
for(int k = 0; k < n; k++)
for(int j = 0; j < m; j++)
A[k * m + j] = A1[k * m + j];
}
double t2 = MPI_Wtime();
MPI_Gather(P, n * m, MPI_DOUBLE, D, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(Q, n * m, MPI_DOUBLE, W, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(A, n * m, MPI_DOUBLE, V, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
return 0;
}
</source>
=== Особенности реализации последовательного алгоритма ===
Для оптимальности работы алгоритма следует реализовать функцию умножения матрицы отражения матрицы отражения на произвольную матрицу, это заметно сократит время время его работы.
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
Алгоритм бидиагонализации матрицы довольно последователен: на каждом шаге считается матрица отражения, и распараллелить прямым образом алгоритм нельзя, поскольку каждый раз матрица отражения строится к новой матрице. Однако можно разбить Матрицу <math> D</math> на блоки, и каждый процесс будет занят умножением матрицы отражения на соответствующий блок.

=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
Параллельная реализация алгоритма наиболее эффективна на машинах с быстрым обменом данных меду процессами
=== Существующие реализации алгоритма ===

== Литература ==
Е. Е. Тыртышников, "Методы численного анализа".

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-12T20:30:34Z

Малков Кирилл: /* Информационный граф */

== Свойства и структура алгоритмов ==
=== Общее описание алгоритма ===
Матрица <math>B = [b_{ij}]</math> называется (верхней) двухдиагональной, или бидиагональной, если <math>b_{ij} = 0</math>, <math>i > j</math> или <math>i + 1 < j</math>. Алгоритм бидиагонализации матрицы методом отражений приводит произвольную <math>n \times n</math> матрицу A к бидиагональному виду
<math>B = PAQ</math>,
где <math>P</math> и <math>Q</math> - произведения конечного числа матриц отражения(или вращения).
=== Математическое описание алгоритма ===
Пусть дана матрица <math>A</math> размера <math>n \times n</math>. Вычисляется бидиагональная матрица <math>D</math> размера <math>n \times n</math>, такая что <math>A = PDQ</math> для некоторых матриц <math>P</math> и <math>Q</math>, являющихся произведением некоторого числа матриц отражения.Сначала мы умножаем A слева на матрицу отражения, аннулирующую все поддиагональные элементы первого столбца. Затем умножаем результат справа на матрицу отражения, аннулирующую элементы первой строки в позициях с 3-й по n-ю. Далее умножением слева аннулируем все поддиагональные элементы второго столбца, затем умножением справа аннулируем получаем нули во второй строке в позициях с 4-й по n-ю и т.д. После каждого умножения на матрицу отражения все ранее полученные нули остаются.
Формулы метода:
Пусть <math>A = [a_1, a_2, ..., a_n],</math>
<math>D_1 = A, P_1 = I, Q_1 = I.</math>
 Если <math>i</math> четное, то <math>D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*)</math>
 <math>u_i = \frac{[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]}{||[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]||_2}, Q_i = (I - 2 * u_{i} * u_{i}^*) * Q_{i-1}, P_i = P_{i-1}.</math>
 Если <math>i</math> нечетное, то <math>D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1},</math>
 <math>u_i = \frac{[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]}{||[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]||_2}, Q_i = Q_{i-1}, P_i = P_{i-1} * (I - 2 * u_{i} * u_{i}^*).</math>

=== Вычислительное ядро алгоритма ===
Вычислительное ядро алгоритма составляет нахождение на каждом этапе матрицы отражения(всего <math>2n - 2</math>).

=== Макроструктура алгоритма ===
Операции скалярного произведения, умножения матрицы на вектор
=== Схема реализации последовательного алгоритма ===
<math>D_1 = A, P_1 = I, Q_1 = I.</math>
<math>1) для четного i : D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*), Q_i = (I - 2 * u_{i} * u_{i}^*) * Q_{i-1}, P_i = P_{i-1};</math>
<math> для нечетного i : D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1}, Q_i = Q_{i-1}, P_i = P_{i-1} * (I - 2 * u_{i} * u_{i}^*).</math>
<math>2) для четного i : u_i = \frac{[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]}{||[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]||_2},</math>
<math> для нечетного i : u_i = \frac{[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]}{||[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]||_2}.</math>
=== Последовательная сложность ===
Для бидиагонализации матрицы порядка n методом отражений при последовательном варианте требуется:
<math>1) 2 * n - 2</math> вычислений квадратного корня,
<math>2) 4 * n - 4 </math>делений,
<math>3) (4 * n - 4)(1 + n + n^2)</math> умножений,
<math>4) (4 * n - 4)(1 + n + n^2)</math> сложений.
Сложность последовательного алгоритма есть <math>O(n^3)</math>.

=== Информационный граф ===
[[Файл:Diagram.png|Информационный граф]]
 <math>I-2uu*</math> - вычисление матрицы отражения;
 <math>D, P, Q</math> - умножение матрицы отражения на матрицы D, P, Q;
 <math>TRANS</math> - транспонирование матрицы D;

=== Ресурс параллелизма алгоритма ===
При параллельном вычислении двухдиагональной матрицы методом отражений требуется последовательно выполнить следующие ярусы:
 -- k ярусов с вычислением <math>4 * (n/k) - 4</math> делений,
 -- k ярусов с вычислением <math>(4 * (n/k) - 4)(1 + n + n^2)</math> умножений,
 -- k ярусов с вычислением <math>(4 * (n/k) - 4)(1 + n + n^2)</math> сложений.
 Получим, что параллельная сложность алгоритма есть <math>O(n^2)</math>.

=== Входные и выходные данные алгоритма ===
На вход подается матрица <math>A</math> и ее размер <math>n</math>. Выходные данные -- матрицы <math>P, Q, D,</math> такие что <math>A = PDQ </math>

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

== Программная реализация алгоритма ==
<source lang="cpp">
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <mkl.h>
#include <time.h>
#include <omp.h>

void matrixA_gen(double *A, int num, int n, int m)
{
int r1 = 0, r2 = 0;
for(int i = 0; i < n; i++)
for(int j = 0; j < m; j++)
{
//A[i * m + j] = -1.0 * num * rand() / (double) RAND_MAX * i * j + 2.0 * rand() / (double) RAND_MAX;
A[i*m + j] = i + 1 + num;
}
}

void matrixP_gen(double *P, int num, int n, int m)
{
for(int i = 0; i < n; i++)
{
for(int j = 0; j < m; j++)
{
if((i % m) == j && i >= num * m && i < (num + 1) * m)
{
P[i * m + j] = 1;
} else {
P[i * m + j] = 0;
}
}
}
}

double sc_mult(double *a, double *b, int n)
{
double sum = 0;
for(int i = 0; i < n; i++)
sum += a[i] * b[i];
return sum;
};

double norm(double *a, int n)
{
return sqrt(sc_mult(a, a, n));
};

void mult_vec(double *A, double *u, int n, int m)
{
char t = 'n';
int int_one = 1;
double *a, one = 1.0, zero = 0.0, min_two = -2.0;
a = (double *)calloc(m, sizeof(*a));
dgemv(&t, &m, &n, &one, A, &m, u, &int_one, &zero, a, &int_one);
dger(&m, &n, &min_two, a, &int_one, u, &int_one, A, &m);
};

void mult(double *A, double *B, int l, int n)
{
double *G;
G = (double *)calloc(n * n, sizeof(*G));
if(l)
{
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
G[i * n + j] += A[i * n + k] * B[k * n + j];
} else {
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
G[i * n + j] += B[i * n + k] * A[k * n + j];
}
for (int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
A[i*n + j] = G[i*n +j];
free(G);
};

void reflect(double *u, double *a, int n, int j)
{
for (int i = 0; i < n; i++)
{
if (i < j)
{
u[i] = 0;
} else {
if (i == j){
u[i] = a[i] - norm(a, n);
} else {
u[i] = a[i];
}
}
}
double sc = norm(u, n);
if(sc != 0)
{
for(int i = 0; i < n; i++)
{
u[i] /= sc;
}
}
};

void trans(double *A, int n)
{
double *B;
B = (double *) calloc(n * n, sizeof(*B));
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
B[i * n + j] = A[j * n + i];
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
A[i * n + j] = B[i * n + j];
free(B);
return;
}

int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
int size, num, error, n, m, i = 0;
double d_norm = 0;
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &num);
n = 64;
srand(time(NULL));
if(size != 0)
{
m = n/size;
}
MPI_Request req;
MPI_Status status;
double *M, *A, *P, *Q, *C, *D, *Df, *Pf, *Qf, *V, *A1, *W;
double *u, *a;
A1 = (double *) calloc(m * n, sizeof(*A1));
V = (double *)calloc(n * n, sizeof(*V));
W = (double *)calloc(n * n, sizeof(*W));
Qf = (double *)calloc(n * n, sizeof(*Qf));
M = (double *)calloc(n * n, sizeof(*M));
Pf = (double *)calloc(n * n, sizeof(*Pf));
A = (double *)calloc(n * m, sizeof(*A));
P = (double *)calloc(n * m, sizeof(*P));
Q = (double *)calloc(n * m, sizeof(*Q));
C = (double *)calloc(m * m, sizeof(*C));
D = (double *)calloc(n * n, sizeof(*D));
Df = (double *)calloc(n * n, sizeof(*Df));
u = (double *)calloc(n, sizeof(*u));
a = (double *)calloc(n, sizeof(*a));
matrixA_gen(A, num, n, m);
matrixP_gen(P, num, n, m);
matrixP_gen(Q, num, n, m);
MPI_Gather(A, n * m, MPI_DOUBLE, D, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if(num == 0)
{
for(int k = 0; k < size; k++)
for(int j = 0; j < n; j++)
for(int h = 0; h < m; h++)
{
M[h + k * m + j * n] = D[k * m * n + j * m + h];
}
}

double t1 = MPI_Wtime();

for(int f = 0; f < 2 * n - 2; f++)
{
i = f / 2;
if(i < num * m || i >= (num + 1) * m )
{
MPI_Recv(u, n, MPI_DOUBLE, (i / m), 0, MPI_COMM_WORLD, &status);
mult_vec(A, u, n, m);
if(!(f%2))
mult_vec(P, u, n, m);
else
mult_vec(Q, u, n, m);
} else {
if(!(f%2))
{
for(int k = 0; k < n; k++)
{
if(k < i)
{
a[k] = 0;
} else {
a[k] = A[(i % m) + k * m];
}
}
reflect(u, a, n, i);
for(int l = 0; l < size; l++)
{
if(l != num)
{
MPI_Isend(u, n, MPI_DOUBLE, l, 0, MPI_COMM_WORLD, &req);
}
}
} else {
for(int k = 0; k < n; k++)
{
if(k < i+1)
{
a[k] = 0;
} else {
a[k] = A[(i % m) + k * m];
}
}
reflect(u, a, n, i+1);
for(int l = 0; l < size; l++)
{
if(l != num)
{
MPI_Isend(u, n, MPI_DOUBLE, l, 0, MPI_COMM_WORLD, &req);
}
}
}
mult_vec(A, u, n, m);
if(!(f%2))
mult_vec(P, u, n, m);
else
mult_vec(Q, u, n, m);
}
for(int k = 0; k < n; k++)
for(int j = 0; j < m; j++)
A1[k * m + j] = A[k * m + j];

for(int j = 0; j < size; j++)
{
if(j != num)
{
for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
C[h * m + l] = A1[j * m * m + h + l * m];
MPI_Isend(C, m * m, MPI_DOUBLE, j, 0, MPI_COMM_WORLD, &req);
} else {
for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
C[h * m + l] = A1[j * m * m + h + l * m];

for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
A1[j * m * m + h + l * m] = C[h + l * m];
}
}

for(int k = 0; k < size; k++)
{
if(k != num)
{
MPI_Recv(C, m * m, MPI_DOUBLE, k, 0, MPI_COMM_WORLD, &status);

for(int l = 0; l < m; l++)
{
for(int h = 0; h < m; h++)
{
A1[k * m * m + l + h * m] = C[h * m + l];
}
}
}
}
for(int k = 0; k < n; k++)
for(int j = 0; j < m; j++)
A[k * m + j] = A1[k * m + j];
}
double t2 = MPI_Wtime();
MPI_Gather(P, n * m, MPI_DOUBLE, D, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(Q, n * m, MPI_DOUBLE, W, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(A, n * m, MPI_DOUBLE, V, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
return 0;
}
</source>
=== Особенности реализации последовательного алгоритма ===
Для оптимальности работы алгоритма следует реализовать функцию умножения матрицы отражения матрицы отражения на произвольную матрицу, это заметно сократит время время его работы.
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
Алгоритм бидиагонализации матрицы довольно последователен: на каждом шаге считается матрица отражения, и распараллелить прямым образом алгоритм нельзя, поскольку каждый раз матрица отражения строится к новой матрице. Однако можно разбить Матрицу <math> D</math> на блоки, и каждый процесс будет занят умножением матрицы отражения на соответствующий блок.

=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
Параллельная реализация алгоритма наиболее эффективна на машинах с быстрым обменом данных меду процессами
=== Существующие реализации алгоритма ===

== Литература ==
Е. Е. Тыртышников, "Методы численного анализа".

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-12T20:28:53Z

Малков Кирилл:

== Свойства и структура алгоритмов ==
=== Общее описание алгоритма ===
Матрица <math>B = [b_{ij}]</math> называется (верхней) двухдиагональной, или бидиагональной, если <math>b_{ij} = 0</math>, <math>i > j</math> или <math>i + 1 < j</math>. Алгоритм бидиагонализации матрицы методом отражений приводит произвольную <math>n \times n</math> матрицу A к бидиагональному виду
<math>B = PAQ</math>,
где <math>P</math> и <math>Q</math> - произведения конечного числа матриц отражения(или вращения).
=== Математическое описание алгоритма ===
Пусть дана матрица <math>A</math> размера <math>n \times n</math>. Вычисляется бидиагональная матрица <math>D</math> размера <math>n \times n</math>, такая что <math>A = PDQ</math> для некоторых матриц <math>P</math> и <math>Q</math>, являющихся произведением некоторого числа матриц отражения.Сначала мы умножаем A слева на матрицу отражения, аннулирующую все поддиагональные элементы первого столбца. Затем умножаем результат справа на матрицу отражения, аннулирующую элементы первой строки в позициях с 3-й по n-ю. Далее умножением слева аннулируем все поддиагональные элементы второго столбца, затем умножением справа аннулируем получаем нули во второй строке в позициях с 4-й по n-ю и т.д. После каждого умножения на матрицу отражения все ранее полученные нули остаются.
Формулы метода:
Пусть <math>A = [a_1, a_2, ..., a_n],</math>
<math>D_1 = A, P_1 = I, Q_1 = I.</math>
 Если <math>i</math> четное, то <math>D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*)</math>
 <math>u_i = \frac{[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]}{||[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]||_2}, Q_i = (I - 2 * u_{i} * u_{i}^*) * Q_{i-1}, P_i = P_{i-1}.</math>
 Если <math>i</math> нечетное, то <math>D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1},</math>
 <math>u_i = \frac{[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]}{||[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]||_2}, Q_i = Q_{i-1}, P_i = P_{i-1} * (I - 2 * u_{i} * u_{i}^*).</math>

=== Вычислительное ядро алгоритма ===
Вычислительное ядро алгоритма составляет нахождение на каждом этапе матрицы отражения(всего <math>2n - 2</math>).

=== Макроструктура алгоритма ===
Операции скалярного произведения, умножения матрицы на вектор
=== Схема реализации последовательного алгоритма ===
<math>D_1 = A, P_1 = I, Q_1 = I.</math>
<math>1) для четного i : D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*), Q_i = (I - 2 * u_{i} * u_{i}^*) * Q_{i-1}, P_i = P_{i-1};</math>
<math> для нечетного i : D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1}, Q_i = Q_{i-1}, P_i = P_{i-1} * (I - 2 * u_{i} * u_{i}^*).</math>
<math>2) для четного i : u_i = \frac{[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]}{||[0, 0, ..., 0, d_{i+1} - ||d_{i+1}||_2, d_{i+2}, ..., d_n]||_2},</math>
<math> для нечетного i : u_i = \frac{[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]}{||[0, 0, ..., 0, d_i - ||d_i||_2, d_{i+1}, ..., d_n]||_2}.</math>
=== Последовательная сложность ===
Для бидиагонализации матрицы порядка n методом отражений при последовательном варианте требуется:
<math>1) 2 * n - 2</math> вычислений квадратного корня,
<math>2) 4 * n - 4 </math>делений,
<math>3) (4 * n - 4)(1 + n + n^2)</math> умножений,
<math>4) (4 * n - 4)(1 + n + n^2)</math> сложений.
Сложность последовательного алгоритма есть <math>O(n^3)</math>.

=== Информационный граф ===
[[Файл:Diagram.png|Информационный граф]]
<math>I-2uu*</math> - вычисление матрицы отражения;
<math>D, P, Q</math> - умножение матрицы отражения на матрицы D, P, Q;
<math>TRANS</math> - транспонирование матрицы D;
=== Ресурс параллелизма алгоритма ===
При параллельном вычислении двухдиагональной матрицы методом отражений требуется последовательно выполнить следующие ярусы:
 -- k ярусов с вычислением <math>4 * (n/k) - 4</math> делений,
 -- k ярусов с вычислением <math>(4 * (n/k) - 4)(1 + n + n^2)</math> умножений,
 -- k ярусов с вычислением <math>(4 * (n/k) - 4)(1 + n + n^2)</math> сложений.
 Получим, что параллельная сложность алгоритма есть <math>O(n^2)</math>.

=== Входные и выходные данные алгоритма ===
На вход подается матрица <math>A</math> и ее размер <math>n</math>. Выходные данные -- матрицы <math>P, Q, D,</math> такие что <math>A = PDQ </math>

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

== Программная реализация алгоритма ==
<source lang="cpp">
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <mkl.h>
#include <time.h>
#include <omp.h>

void matrixA_gen(double *A, int num, int n, int m)
{
int r1 = 0, r2 = 0;
for(int i = 0; i < n; i++)
for(int j = 0; j < m; j++)
{
//A[i * m + j] = -1.0 * num * rand() / (double) RAND_MAX * i * j + 2.0 * rand() / (double) RAND_MAX;
A[i*m + j] = i + 1 + num;
}
}

void matrixP_gen(double *P, int num, int n, int m)
{
for(int i = 0; i < n; i++)
{
for(int j = 0; j < m; j++)
{
if((i % m) == j && i >= num * m && i < (num + 1) * m)
{
P[i * m + j] = 1;
} else {
P[i * m + j] = 0;
}
}
}
}

double sc_mult(double *a, double *b, int n)
{
double sum = 0;
for(int i = 0; i < n; i++)
sum += a[i] * b[i];
return sum;
};

double norm(double *a, int n)
{
return sqrt(sc_mult(a, a, n));
};

void mult_vec(double *A, double *u, int n, int m)
{
char t = 'n';
int int_one = 1;
double *a, one = 1.0, zero = 0.0, min_two = -2.0;
a = (double *)calloc(m, sizeof(*a));
dgemv(&t, &m, &n, &one, A, &m, u, &int_one, &zero, a, &int_one);
dger(&m, &n, &min_two, a, &int_one, u, &int_one, A, &m);
};

void mult(double *A, double *B, int l, int n)
{
double *G;
G = (double *)calloc(n * n, sizeof(*G));
if(l)
{
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
G[i * n + j] += A[i * n + k] * B[k * n + j];
} else {
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
G[i * n + j] += B[i * n + k] * A[k * n + j];
}
for (int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
A[i*n + j] = G[i*n +j];
free(G);
};

void reflect(double *u, double *a, int n, int j)
{
for (int i = 0; i < n; i++)
{
if (i < j)
{
u[i] = 0;
} else {
if (i == j){
u[i] = a[i] - norm(a, n);
} else {
u[i] = a[i];
}
}
}
double sc = norm(u, n);
if(sc != 0)
{
for(int i = 0; i < n; i++)
{
u[i] /= sc;
}
}
};

void trans(double *A, int n)
{
double *B;
B = (double *) calloc(n * n, sizeof(*B));
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
B[i * n + j] = A[j * n + i];
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
A[i * n + j] = B[i * n + j];
free(B);
return;
}

int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
int size, num, error, n, m, i = 0;
double d_norm = 0;
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &num);
n = 64;
srand(time(NULL));
if(size != 0)
{
m = n/size;
}
MPI_Request req;
MPI_Status status;
double *M, *A, *P, *Q, *C, *D, *Df, *Pf, *Qf, *V, *A1, *W;
double *u, *a;
A1 = (double *) calloc(m * n, sizeof(*A1));
V = (double *)calloc(n * n, sizeof(*V));
W = (double *)calloc(n * n, sizeof(*W));
Qf = (double *)calloc(n * n, sizeof(*Qf));
M = (double *)calloc(n * n, sizeof(*M));
Pf = (double *)calloc(n * n, sizeof(*Pf));
A = (double *)calloc(n * m, sizeof(*A));
P = (double *)calloc(n * m, sizeof(*P));
Q = (double *)calloc(n * m, sizeof(*Q));
C = (double *)calloc(m * m, sizeof(*C));
D = (double *)calloc(n * n, sizeof(*D));
Df = (double *)calloc(n * n, sizeof(*Df));
u = (double *)calloc(n, sizeof(*u));
a = (double *)calloc(n, sizeof(*a));
matrixA_gen(A, num, n, m);
matrixP_gen(P, num, n, m);
matrixP_gen(Q, num, n, m);
MPI_Gather(A, n * m, MPI_DOUBLE, D, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if(num == 0)
{
for(int k = 0; k < size; k++)
for(int j = 0; j < n; j++)
for(int h = 0; h < m; h++)
{
M[h + k * m + j * n] = D[k * m * n + j * m + h];
}
}

double t1 = MPI_Wtime();

for(int f = 0; f < 2 * n - 2; f++)
{
i = f / 2;
if(i < num * m || i >= (num + 1) * m )
{
MPI_Recv(u, n, MPI_DOUBLE, (i / m), 0, MPI_COMM_WORLD, &status);
mult_vec(A, u, n, m);
if(!(f%2))
mult_vec(P, u, n, m);
else
mult_vec(Q, u, n, m);
} else {
if(!(f%2))
{
for(int k = 0; k < n; k++)
{
if(k < i)
{
a[k] = 0;
} else {
a[k] = A[(i % m) + k * m];
}
}
reflect(u, a, n, i);
for(int l = 0; l < size; l++)
{
if(l != num)
{
MPI_Isend(u, n, MPI_DOUBLE, l, 0, MPI_COMM_WORLD, &req);
}
}
} else {
for(int k = 0; k < n; k++)
{
if(k < i+1)
{
a[k] = 0;
} else {
a[k] = A[(i % m) + k * m];
}
}
reflect(u, a, n, i+1);
for(int l = 0; l < size; l++)
{
if(l != num)
{
MPI_Isend(u, n, MPI_DOUBLE, l, 0, MPI_COMM_WORLD, &req);
}
}
}
mult_vec(A, u, n, m);
if(!(f%2))
mult_vec(P, u, n, m);
else
mult_vec(Q, u, n, m);
}
for(int k = 0; k < n; k++)
for(int j = 0; j < m; j++)
A1[k * m + j] = A[k * m + j];

for(int j = 0; j < size; j++)
{
if(j != num)
{
for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
C[h * m + l] = A1[j * m * m + h + l * m];
MPI_Isend(C, m * m, MPI_DOUBLE, j, 0, MPI_COMM_WORLD, &req);
} else {
for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
C[h * m + l] = A1[j * m * m + h + l * m];

for(int l = 0; l < m; l++)
for(int h = 0; h < m; h++)
A1[j * m * m + h + l * m] = C[h + l * m];
}
}

for(int k = 0; k < size; k++)
{
if(k != num)
{
MPI_Recv(C, m * m, MPI_DOUBLE, k, 0, MPI_COMM_WORLD, &status);

for(int l = 0; l < m; l++)
{
for(int h = 0; h < m; h++)
{
A1[k * m * m + l + h * m] = C[h * m + l];
}
}
}
}
for(int k = 0; k < n; k++)
for(int j = 0; j < m; j++)
A[k * m + j] = A1[k * m + j];
}
double t2 = MPI_Wtime();
MPI_Gather(P, n * m, MPI_DOUBLE, D, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(Q, n * m, MPI_DOUBLE, W, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(A, n * m, MPI_DOUBLE, V, n * m, MPI_DOUBLE, 0, MPI_COMM_WORLD);
return 0;
}
</source>
=== Особенности реализации последовательного алгоритма ===
Для оптимальности работы алгоритма следует реализовать функцию умножения матрицы отражения матрицы отражения на произвольную матрицу, это заметно сократит время время его работы.
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
Алгоритм бидиагонализации матрицы довольно последователен: на каждом шаге считается матрица отражения, и распараллелить прямым образом алгоритм нельзя, поскольку каждый раз матрица отражения строится к новой матрице. Однако можно разбить Матрицу <math> D</math> на блоки, и каждый процесс будет занят умножением матрицы отражения на соответствующий блок.

=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
Параллельная реализация алгоритма наиболее эффективна на машинах с быстрым обменом данных меду процессами
=== Существующие реализации алгоритма ===

== Литература ==
Е. Е. Тыртышников, "Методы численного анализа".

Файл:Diagram.png

2021-12-12T20:10:35Z

Малков Кирилл:

Информационный граф

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T22:11:24Z

Малков Кирилл: /* Возможные способы и особенности параллельной реализации алгоритма */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T22:10:43Z

Малков Кирилл: /* Входные и выходные данные алгоритма */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T22:09:03Z

Малков Кирилл: /* Вычислительное ядро алгоритма */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T22:07:57Z

Малков Кирилл: /* Математическое описание алгоритма */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T22:07:30Z

Малков Кирилл: /* Математическое описание алгоритма */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T22:03:10Z

Малков Кирилл: /* Свойства и структура алгоритмов */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T22:02:11Z

Малков Кирилл: /* Информационный граф */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T22:01:20Z

Малков Кирилл: /* Информационный граф */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T22:00:10Z

Малков Кирилл: /* Информационный граф */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T21:59:24Z

Малков Кирилл: /* Информационный граф */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T21:58:36Z

Малков Кирилл: /* Последовательная сложность */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T21:57:22Z

Малков Кирилл:

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T21:55:38Z

Малков Кирилл: /* Последовательная сложность */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T21:54:40Z

Малков Кирилл: /* Информационный граф */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T21:53:24Z

Малков Кирилл: /* Информационный граф */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T03:24:56Z

Малков Кирилл: /* Математическое описание алгоритма */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T03:24:01Z

Малков Кирилл:

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-08T03:14:08Z

Малков Кирилл:

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-06T23:57:47Z

Малков Кирилл:

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-06T23:00:33Z

Малков Кирилл: /* Свойства и структура алгоритмов */

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-12-06T22:42:31Z

Малков Кирилл:

Участник:Малков Кирилл/Алгоритм бидиагонализации матрицы методом отражений

2021-10-26T19:07:42Z

Малков Кирилл: Новая страница: «== Свойства и структура алгоритмов == === Общее описание алгоритма === Матрица <math>B = [b_{ij}]</math>...»

Алгоритм бидиагонализации матрицы методом отражений

2021-10-26T19:05:00Z