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

Свойства и структура алгоритмов

Общее описание алгоритма

Матрица [math]\displaystyle{ B = [b_{ij}] }[/math] называется (верхней) двухдиагональной, или бидиагональной, если [math]\displaystyle{ b_{ij} = 0 }[/math], [math]\displaystyle{ i \gt j }[/math] или [math]\displaystyle{ i + 1 \lt j }[/math]. Алгоритм бидиагонализации матрицы методом отражений приводит произвольную [math]\displaystyle{ n \times n }[/math] матрицу A к бидиагональному виду [math]\displaystyle{ B = PAQ }[/math], где [math]\displaystyle{ P }[/math] и [math]\displaystyle{ Q }[/math] - произведения конечного числа матриц отражения(или вращения).

Математическое описание алгоритма

Пусть дана матрица [math]\displaystyle{ A }[/math] размера [math]\displaystyle{ n \times n }[/math]. Вычисляется бидиагональная матрица [math]\displaystyle{ D }[/math] размера [math]\displaystyle{ n \times n }[/math], такая что [math]\displaystyle{ A = PDQ }[/math] для некоторых матриц [math]\displaystyle{ P }[/math] и [math]\displaystyle{ Q }[/math], являющихся произведением некоторого числа матриц отражения.Сначала мы умножаем A слева на матрицу отражения, аннулирующую все поддиагональные элементы первого столбца. Затем умножаем результат справа на матрицу отражения, аннулирующую элементы первой строки в позициях с 3-й по n-ю. Далее умножением слева аннулируем все поддиагональные элементы второго столбца, затем умножением справа аннулируем получаем нули во второй строке в позициях с 4-й по n-ю и т.д. После каждого умножения на матрицу отражения все ранее полученные нули остаются. Формулы метода: Пусть [math]\displaystyle{ A = [a_1, a_2, ..., a_n], }[/math] [math]\displaystyle{ D_1 = A, P_1 = I, Q_1 = I. }[/math]
Если [math]\displaystyle{ i }[/math] четное, то [math]\displaystyle{ D_i = D_{i-1} * (I - 2 * u_{i-1} * u_{i-1}^*) }[/math]
[math]\displaystyle{ 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]\displaystyle{ i }[/math] нечетное, то [math]\displaystyle{ D_i = (I - 2 * u_{i-1} * u_{i-1}^*) * D_{i-1}, }[/math]
[math]\displaystyle{ 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]\displaystyle{ 2n - 2 }[/math]).

Макроструктура алгоритма

Операции скалярного произведения, умножения матрицы на вектор

Схема реализации последовательного алгоритма

[math]\displaystyle{ D_1 = A, P_1 = I, Q_1 = I. }[/math] [math]\displaystyle{ 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]\displaystyle{ для нечетного 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]\displaystyle{ 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]\displaystyle{ для нечетного 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]\displaystyle{ 1) 2 * n - 2 }[/math] вычислений квадратного корня, [math]\displaystyle{ 2) 4 * n - 4 }[/math]делений, [math]\displaystyle{ 3) (4 * n - 4)(1 + n + n^2) }[/math] умножений, [math]\displaystyle{ 4) (4 * n - 4)(1 + n + n^2) }[/math] сложений. Сложность последовательного алгоритма есть [math]\displaystyle{ O(n^3) }[/math].

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

Информационный граф для количества процессов равного 4.
[math]\displaystyle{ I-2uu^* }[/math] - вычисление матрицы отражения;
[math]\displaystyle{ D, P, Q }[/math] - умножение матрицы отражения на матрицы D, P, Q;
[math]\displaystyle{ TRANS }[/math] - транспонирование матрицы D;

Ресурс параллелизма алгоритма

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

Входные и выходные данные алгоритма

На вход подается матрица [math]\displaystyle{ A }[/math] и ее размер [math]\displaystyle{ n }[/math]. Выходные данные -- матрицы [math]\displaystyle{ P, Q, D, }[/math] такие что [math]\displaystyle{ A = PDQ }[/math]

Свойства алгоритма

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

Программная реализация алгоритма

#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;
}

Особенности реализации последовательного алгоритма

Для оптимальности работы алгоритма следует реализовать функцию умножения матрицы отражения матрицы отражения на произвольную матрицу, это заметно сократит время время его работы.

Возможные способы и особенности параллельной реализации алгоритма

Алгоритм бидиагонализации матрицы довольно последователен: на каждом шаге считается матрица отражения, и распараллелить прямым образом алгоритм нельзя, поскольку каждый раз матрица отражения строится к новой матрице. Однако можно разбить Матрицу [math]\displaystyle{ D }[/math] на блоки, и каждый процесс будет занят умножением матрицы отражения на соответствующий блок.

Выводы для классов архитектур

Параллельная реализация алгоритма наиболее эффективна на машинах с быстрым обменом данных меду процессами

Литература

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