Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц: различия между версиями
V (обсуждение | вклад) |
V (обсуждение | вклад) |
||
(не показано 7 промежуточных версий этого же участника) | |||
Строка 42: | Строка 42: | ||
== Последовательная сложность алгоритма == | == Последовательная сложность алгоритма == | ||
− | + | Сложность последовательного алгоритма порядка <math>O(n^{3})</math> | |
== Информационный граф == | == Информационный граф == | ||
Строка 78: | Строка 78: | ||
Для исследования масштабируемости были проведены серии тестов с различным числом узлов (в степенях 2). Чем больше узлов, тем меньше загруженность каждого, если хотим получить одинаковый совокупный объем данных, а, следовательно, программа выполняется быстрее. | Для исследования масштабируемости были проведены серии тестов с различным числом узлов (в степенях 2). Чем больше узлов, тем меньше загруженность каждого, если хотим получить одинаковый совокупный объем данных, а, следовательно, программа выполняется быстрее. | ||
График для фиксированных n = 100, p = 10. | График для фиксированных n = 100, p = 10. | ||
− | [[Файл:Graphhh.jpg|center | + | [[Файл:Graphhh.jpg|center]] |
Все тесты проводились на [https://ru.wikipedia.org/wiki/Ломоносов_(суперкомпьютер) суперкомпьютере Ломоносов]. | Все тесты проводились на [https://ru.wikipedia.org/wiki/Ломоносов_(суперкомпьютер) суперкомпьютере Ломоносов]. | ||
Строка 89: | Строка 89: | ||
== Существующие реализации алгоритма == | == Существующие реализации алгоритма == | ||
− | + | <source lang="c++"> | |
+ | #include "mkl.h" | ||
+ | #include <ctime> | ||
+ | #include <cstdlib> | ||
+ | #include <iostream> | ||
+ | #include <string.h> | ||
+ | #include <sstream> | ||
+ | #include "mpi.h" | ||
+ | #include <fstream> | ||
+ | #include <algorithm> | ||
+ | #include <iterator> | ||
+ | #include <string> | ||
+ | #include <vector> | ||
+ | |||
+ | using namespace std; | ||
+ | |||
+ | void diagMatrix(int p, int k, float *m); | ||
+ | void ortBasis(int p, float *m); | ||
+ | void out(int p, float *m); | ||
+ | void choleskyFactorization(int p, float *m); | ||
+ | void copy(int p, float *from, float *to); | ||
+ | void samplesGenerate(int n, int p, float *mean, float *chol, float *res); | ||
+ | void wGenerate(int n, float* w); | ||
+ | void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result); | ||
+ | void sampleCovariance(int n, int p, float* Samples, float* Result); | ||
+ | void eigenDecomposition(int p, float* m, float *values); | ||
+ | void rProjector(int p, int r, float *vectors, float* projector); | ||
+ | float S(int n, int p, float *m1, float *m2); | ||
+ | |||
+ | void diagMatrix(int p, int k, float *m) | ||
+ | { | ||
+ | float numbers[5] = {36, 30, 25, 19, 15}; | ||
+ | int ind = 0; | ||
+ | for (int i = 0; i < p; ++i) | ||
+ | { | ||
+ | for (int j = 0; j < p; ++j) | ||
+ | m[i*p + j] = 0.0; | ||
+ | if (k > 0) | ||
+ | { | ||
+ | m[i*p + i] = numbers[ind]; | ||
+ | ind++; | ||
+ | k--; | ||
+ | } | ||
+ | else | ||
+ | m[i*p + i] = 1.0 + 5.0*rand() / float(RAND_MAX); | ||
+ | } | ||
+ | } | ||
+ | |||
+ | void ortBasis(int p, float *m) | ||
+ | { | ||
+ | for (int i = 0; i < p; i++) | ||
+ | for(int j = 0; j < p; j++) | ||
+ | m[i*p + j] = rand(); | ||
+ | float* tau = (float*) calloc(p, sizeof(float)); | ||
+ | LAPACKE_sgeqrf(LAPACK_ROW_MAJOR, p, p, m, p, tau); | ||
+ | LAPACKE_sorgqr(LAPACK_ROW_MAJOR, p, p, 1, m, p, tau); | ||
+ | free(tau); | ||
+ | } | ||
+ | |||
+ | void out(int r, int c, float *m) | ||
+ | { | ||
+ | for (int i = 0; i < r; i++) | ||
+ | { | ||
+ | cout << "{ "; | ||
+ | for (int j = 0; j < c; j++) | ||
+ | cout << m[i*c + j] << " "; | ||
+ | cout << "}" << endl; | ||
+ | } | ||
+ | } | ||
+ | |||
+ | void choleskyFactorization(int p, float *m) | ||
+ | { | ||
+ | lapack_int q = LAPACKE_spotrf(LAPACK_ROW_MAJOR, 'L', p, m, p); | ||
+ | } | ||
+ | |||
+ | void copy(int p, float *from, float *to) | ||
+ | { | ||
+ | mkl_somatcopy('R', 'N', p, p, 1.0, from, p, to, p); | ||
+ | } | ||
+ | |||
+ | void samplesGenerate(int n, int p, float *mean, float *chol, float *res) | ||
+ | { | ||
+ | VSLStreamStatePtr stream; | ||
+ | int i = (int) rand(); | ||
+ | vslNewStream(&stream,VSL_BRNG_MCG31, i); | ||
+ | vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, n, res, p, VSL_MATRIX_STORAGE_FULL, mean, chol); | ||
+ | } | ||
+ | |||
+ | void wGenerate(int n, float* w) | ||
+ | { | ||
+ | float *mean = (float *) calloc(n, sizeof(float)); | ||
+ | float *cov = (float*) calloc(n * n, sizeof(float)); | ||
+ | for (int i = 1; i < n; i++) | ||
+ | { | ||
+ | mean[i] = 1; | ||
+ | cov[i*n + i] = 1; | ||
+ | } | ||
+ | int i = (int) rand(); | ||
+ | VSLStreamStatePtr stream; | ||
+ | vslNewStream(&stream,VSL_BRNG_MCG31, i); | ||
+ | vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, 1, w, n, VSL_MATRIX_STORAGE_FULL, mean, cov); | ||
+ | free(mean); | ||
+ | free(cov); | ||
+ | } | ||
+ | |||
+ | void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result) | ||
+ | { | ||
+ | //Result p*p | ||
+ | float *W = (float *) calloc(n, sizeof(float)); | ||
+ | float *X = (float *) calloc(p, sizeof(float)); | ||
+ | wGenerate(n, W); | ||
+ | for(int k = 0; k < n; k++) | ||
+ | { | ||
+ | for (int i = 0; i < p; i++) | ||
+ | { | ||
+ | X[i] = Samples[k*p + i]; | ||
+ | } | ||
+ | cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, W[k], X, p, X, p, 1.0, Result, p); | ||
+ | } | ||
+ | free(W); | ||
+ | for(int k = 0; k < p; k++) | ||
+ | { | ||
+ | for (int i = 0; i < p; i++) | ||
+ | { | ||
+ | Result[k*p + i] /= n; | ||
+ | } | ||
+ | } | ||
+ | |||
+ | free(X); | ||
+ | } | ||
+ | |||
+ | void sampleCovariance(int n, int p, float* Samples, float* Result) | ||
+ | { | ||
+ | //Result p*p | ||
+ | float *X = (float *) calloc(p, sizeof(float)); | ||
+ | for(int k = 0; k < n; k++) | ||
+ | { | ||
+ | for (int i = 0; i < p; i++) | ||
+ | { | ||
+ | X[i] = Samples[k*p + i]; | ||
+ | } | ||
+ | cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, Result, p); | ||
+ | } | ||
+ | for(int k = 0; k < p; k++) | ||
+ | { | ||
+ | for (int i = 0; i < p; i++) | ||
+ | { | ||
+ | Result[k*p + i] /= n; | ||
+ | } | ||
+ | } | ||
+ | free(X); | ||
+ | } | ||
+ | |||
+ | void eigenDecomposition(int p, float* m, float *values) //decreasing orderб столбыц! | ||
+ | { | ||
+ | LAPACKE_ssyevd(LAPACK_ROW_MAJOR, 'V', 'U', p, m, p, values); | ||
+ | } | ||
+ | |||
+ | void rProjector(int p, int r, float *vectors, float* projector) | ||
+ | { | ||
+ | int num = p - r; | ||
+ | float *X = (float *) calloc(p, sizeof(float)); | ||
+ | for (int i = 0; i < p; i++) | ||
+ | { | ||
+ | X[i] = vectors[p*i + num]; | ||
+ | } | ||
+ | cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, projector, p); | ||
+ | free(X); | ||
+ | } | ||
+ | |||
+ | float S(int n, int p, float *m1, float *m2) | ||
+ | { | ||
+ | float *m = (float *) calloc(p * p, sizeof(float)); | ||
+ | for (int i = 0; i < p; i++) | ||
+ | for (int j = 0; j < p; j++) | ||
+ | m[i*p + j] = m1[i*p + j] - m2[i*p + j]; | ||
+ | float norm = LAPACKE_slange(LAPACK_ROW_MAJOR, 'F', p, p, m, p); | ||
+ | free(m); | ||
+ | return n*norm*norm; | ||
+ | } | ||
+ | |||
+ | |||
+ | |||
+ | int main(int argc, char **argv) | ||
+ | { | ||
+ | //задаем размерности | ||
+ | int paral = 16; | ||
+ | int p = 20; | ||
+ | int k = 3; | ||
+ | int M = 3000; | ||
+ | int r = 1; //нумерация с 1 | ||
+ | int N = 3; //45 | ||
+ | int n = 100; | ||
+ | int rank, size, prev, next; | ||
+ | |||
+ | //инициализация MPI | ||
+ | MPI_Request reqs[4]; | ||
+ | MPI_Status stats[4]; | ||
+ | MPI_Init(&argc, &argv); | ||
+ | MPI_Comm_size(MPI_COMM_WORLD, &size); | ||
+ | MPI_Comm_rank(MPI_COMM_WORLD, &rank); | ||
+ | // | ||
+ | double MPI_Wtime(void); | ||
+ | double t1, t2, dt; | ||
+ | // | ||
+ | t1 = MPI_Wtime(); | ||
+ | prev = rank - 1; | ||
+ | next = rank + 1; | ||
+ | if(rank==0) prev = size - 1; | ||
+ | if(rank==size - 1) next = 0; | ||
+ | |||
+ | //чтобы у всех сгенерировалась одинаковая матрица ковариации srand сработает одновременно | ||
+ | srand(time(NULL)); | ||
+ | //файлы, в который будем писать результаты, у каждого параллельного процесса свой (номер процесса + .txt) | ||
+ | ofstream myfile; | ||
+ | char intStr[2]; | ||
+ | if (rank < 10) | ||
+ | { | ||
+ | intStr[0] = '0'; | ||
+ | intStr[1] = '0' + rank; | ||
+ | } | ||
+ | else | ||
+ | { | ||
+ | intStr[0] = '1'; | ||
+ | intStr[1] = '0' + rank - 10; | ||
+ | } | ||
+ | char *txt = ".txt"; | ||
+ | char *str = strcat(intStr,txt); | ||
+ | |||
+ | myfile.open(str); | ||
+ | |||
+ | //файлы, в которые будем отправлять размерности и тд, за это будет отвечать нулевой процесс | ||
+ | ofstream dim; | ||
+ | ofstream val; | ||
+ | |||
+ | if (rank == 0){ | ||
+ | dim.open("dimensions.txt"); | ||
+ | val.open("values.txt"); | ||
+ | } | ||
+ | |||
+ | if (rank == 0) | ||
+ | { | ||
+ | dim << "p = " << p << endl; | ||
+ | dim << "k = " << k << endl; | ||
+ | dim << "M = " << M << endl; | ||
+ | dim << "r = " << r << endl; | ||
+ | dim << "N = " << N*(paral-1) << endl << endl; | ||
+ | dim << "n = " << n; | ||
+ | } | ||
+ | |||
+ | //генерация матрицы ковариации | ||
+ | float *COV = (float *) calloc(p * p, sizeof(float)); | ||
+ | float *MEAN = (float *) calloc(p, sizeof(float)); | ||
+ | float *COV_Chol = (float *) calloc(p * p, sizeof(float)); | ||
+ | float *Basis = (float *) calloc(p * p, sizeof(float)); | ||
+ | float *COV_L = (float *) calloc(p * p, sizeof(float)); | ||
+ | |||
+ | ortBasis(p, Basis); | ||
+ | |||
+ | diagMatrix(p, k, COV_L); | ||
+ | cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, Basis, p, COV_L, p, 0.0, COV, p); | ||
+ | cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, COV, p, Basis, p, 0.0, COV_L, p); | ||
+ | copy(p, COV_L, COV); | ||
+ | |||
+ | copy(p, COV, COV_Chol); | ||
+ | choleskyFactorization(p, COV_Chol); | ||
+ | free(Basis); | ||
+ | free(COV_L); | ||
+ | |||
+ | float *values_pr = (float *) calloc(p, sizeof(float)); | ||
+ | float *vectors_pr = (float *) calloc(p * p, sizeof(float)); | ||
+ | |||
+ | copy(p, COV, vectors_pr); | ||
+ | eigenDecomposition(p, vectors_pr, values_pr); | ||
+ | |||
+ | //для исследования нам понадобятся собственные значения, поэтому их тоже выводит нулевой процесс | ||
+ | if(rank ==0){ | ||
+ | for (int i = 0; i < p; i++) | ||
+ | { | ||
+ | val << values_pr[i] << " "; | ||
+ | } | ||
+ | val << endl << endl; | ||
+ | } | ||
+ | |||
+ | free(values_pr); | ||
+ | free(vectors_pr); | ||
+ | |||
+ | //вычислим проектор для реальной матрицы ковариации | ||
+ | float *values = (float *) calloc(p, sizeof(float)); | ||
+ | float *vectors = (float *) calloc(p * p, sizeof(float)); | ||
+ | float *ProjectorR = (float *) calloc(p * p, sizeof(float)); | ||
+ | |||
+ | copy(p, COV, vectors); | ||
+ | eigenDecomposition(p, vectors, values); | ||
+ | rProjector(p, r, vectors, ProjectorR); | ||
+ | //ProjectorR - нужный проектор | ||
+ | free(values); | ||
+ | free(vectors); | ||
+ | |||
+ | |||
+ | //для каждого процесса запустим свой рандом (зависящий от ранка), иначе все значения у каждого процесса будут одинаковые | ||
+ | srand(time(NULL) + rank); | ||
+ | |||
+ | float* NormsBootstrapWorld = (float *) calloc(N * M, sizeof(float)); | ||
+ | for (int i = 0; i < N; i++) | ||
+ | { | ||
+ | float *Samples = (float *) calloc(n * p, sizeof(float)); //Каждый Xj - строка в матрице | ||
+ | for (int j = 0; j < M; j++) | ||
+ | { | ||
+ | float *COV_SampleBootstrap = (float *) calloc(p * p, sizeof(float)); | ||
+ | float *values = (float *) calloc(p, sizeof(float)); | ||
+ | float *vectors = (float *) calloc(p * p, sizeof(float)); | ||
+ | float *ProjectorSB = (float *) calloc(p * p, sizeof(float)); | ||
+ | |||
+ | sampleCovarianceBootstrap(n, p, Samples, COV_SampleBootstrap); | ||
+ | copy(p, COV_SampleBootstrap, vectors); | ||
+ | eigenDecomposition(p, vectors, values); //вектора в столбцах | ||
+ | rProjector(p, r, vectors, ProjectorSB); | ||
+ | //ProjectorSB - нужный проектор | ||
+ | |||
+ | float s = S(n, p, ProjectorR, ProjectorSB); | ||
+ | NormsBootstrapWorld[i*M + j] = s; | ||
+ | |||
+ | free(COV_SampleBootstrap); | ||
+ | free(values); | ||
+ | free(vectors); | ||
+ | free(ProjectorSB); | ||
+ | } | ||
+ | free(Samples); | ||
+ | } | ||
+ | |||
+ | free(ProjectorR); | ||
+ | |||
+ | for (int i = 0; i < M; i++) | ||
+ | { | ||
+ | for (int j = 0; j < N; j++) | ||
+ | { | ||
+ | myfile << NormsBootstrapWorld[j*M + i] << " "; | ||
+ | } | ||
+ | myfile << endl; | ||
+ | } | ||
+ | //в myfile N столбцов | ||
+ | |||
+ | myfile << endl; | ||
+ | free(NormsBootstrapWorld); | ||
+ | |||
+ | free(COV); | ||
+ | free(MEAN); | ||
+ | free(COV_Chol); | ||
+ | t2 = MPI_Wtime(); | ||
+ | dt = t2 - t1; | ||
+ | if (rank == 0) | ||
+ | cout << dt << endl; | ||
+ | MPI_Finalize(); //конец параллеливания | ||
+ | } | ||
+ | </source> | ||
= Литература = | = Литература = | ||
− | + | * [https://en.wikipedia.org/wiki/Bootstrapping_(statistics) Boostrapping (wiki)] | |
+ | * [https://software.intel.com/sites/default/files/managed/ff/c8/mkl-2017-developer-reference-c_0.pdf Developer Reference for Intel® Math Kernel Library 2017 - C(pdf)] |
Текущая версия на 15:32, 30 ноября 2016
Основные авторы описания: В.С.Шумовская
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Пусть [math]X_{1},\dots, X_{n}[/math] -- независимые и нормально распределенные случайные вектора в [math]\R^p[/math] с нулевым средним и матрицей ковариацией [math]\Sigma[/math], она лежит в [math]\R^{pxp}[/math] и такая, что ее собственные значения быстро убывают, т.е. 3-5 больших, а остальные, к примеру, в диапазоне от 1 до 3.
К этой выборке применим бутстрэп и вычислим в мире бутстрэпа [math]M[/math] матриц ковариаций сигма [math]\Sigma^{o}_{j}, j = 1,\dots,M.[/math]
Далее фиксируем [math]r[/math], обозначим за [math]P_{r}, P^{o}_{j}, j = 1,\dots,M[/math] - проекторы на r-ое подпространство и вычислим следующие статистики:
[math]S^{o}_{j} = ||P_{r} - P^{o}_{r}||^{2}_{2}[/math]
Задача -- вычислить большое количество этих статистик для визуализации их распределения.
1.2 Математическое описание алгоритма
- Генерация матрицы ковариации.
- Возьмем нужный нам набор собственных значений и поместим их на диагонали матрицы [math]\Lambda.[/math]
- Далее генерируем базис [math]p[/math]-мерного пространства и поместим его в матрицу [math]U.[/math]
- Тогда [math]\Sigma = U\Lambda U^{T}.[/math]
- Проделываем следующий цикл N раз:
- Генерация векторов [math]X[/math] ~ [math]N(\theta, \Sigma).[/math]
- Проделываем следующий цикл M раз:
- Генерация весов [math]w[/math] ~ [math]N(1, 1).[/math]
- Вычисление матриц ковариаций в бутстрэпе [math]\Sigma^{o} = \frac{1}{n}\sum_{k = 1}^{n}w_{k}X_{k}X^{T}_{k}. [/math]
- Вычисление собственных векторов [math]u_{1},\dots,u_{p}[/math] (в порядке убывания).
- Вычисление проекторов (считаем, что одномерные) по формуле [math]u = u_{r}u^{T}_{r}.[/math]
- Вычисление статистик [math]S^{o} = ||P_{r} - P^{o}_{r}||^{2}_{2}.[/math]
1.3 Вычислительное ядро алгоритма
Основное время работы алгоритма приходится на работу с матрицами (присутствует очень много умножений).
1.4 Макроструктура алгоритма
1.5 Схема реализации последовательного алгоритма
1.6 Последовательная сложность алгоритма
Сложность последовательного алгоритма порядка [math]O(n^{3})[/math]
1.7 Информационный граф
1.8 Ресурс параллелизма алгоритма
В принципе, параллельный алгоритм не сильно отличается от схемы последовательного алгоритма, но внешний цикл будет выполняться уже не условные N раз, а N деленное на число вычислительных узлов, поэтому программа будет работать быстрее. Так же на каждом узле параллельно вычисляются внутри узла многие операции (умножение матриц и т.д.)
1.9 Входные и выходные данные алгоритма
1.9.1 Входные данные
- Размерности матриц
- Количество повторений циклов
- Собственные значения
Эти параметры меняются в коде
1.9.2 Выходные данные
От каждого вычислительного узла будет получена матрица значений статистик
1.10 Свойства алгоритма
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Для исследования масштабируемости были проведены серии тестов с различным числом узлов (в степенях 2). Чем больше узлов, тем меньше загруженность каждого, если хотим получить одинаковый совокупный объем данных, а, следовательно, программа выполняется быстрее. График для фиксированных n = 100, p = 10.
Все тесты проводились на суперкомпьютере Ломоносов. Характеристики вычислительного кластера: link.
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
#include "mkl.h"
#include <ctime>
#include <cstdlib>
#include <iostream>
#include <string.h>
#include <sstream>
#include "mpi.h"
#include <fstream>
#include <algorithm>
#include <iterator>
#include <string>
#include <vector>
using namespace std;
void diagMatrix(int p, int k, float *m);
void ortBasis(int p, float *m);
void out(int p, float *m);
void choleskyFactorization(int p, float *m);
void copy(int p, float *from, float *to);
void samplesGenerate(int n, int p, float *mean, float *chol, float *res);
void wGenerate(int n, float* w);
void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result);
void sampleCovariance(int n, int p, float* Samples, float* Result);
void eigenDecomposition(int p, float* m, float *values);
void rProjector(int p, int r, float *vectors, float* projector);
float S(int n, int p, float *m1, float *m2);
void diagMatrix(int p, int k, float *m)
{
float numbers[5] = {36, 30, 25, 19, 15};
int ind = 0;
for (int i = 0; i < p; ++i)
{
for (int j = 0; j < p; ++j)
m[i*p + j] = 0.0;
if (k > 0)
{
m[i*p + i] = numbers[ind];
ind++;
k--;
}
else
m[i*p + i] = 1.0 + 5.0*rand() / float(RAND_MAX);
}
}
void ortBasis(int p, float *m)
{
for (int i = 0; i < p; i++)
for(int j = 0; j < p; j++)
m[i*p + j] = rand();
float* tau = (float*) calloc(p, sizeof(float));
LAPACKE_sgeqrf(LAPACK_ROW_MAJOR, p, p, m, p, tau);
LAPACKE_sorgqr(LAPACK_ROW_MAJOR, p, p, 1, m, p, tau);
free(tau);
}
void out(int r, int c, float *m)
{
for (int i = 0; i < r; i++)
{
cout << "{ ";
for (int j = 0; j < c; j++)
cout << m[i*c + j] << " ";
cout << "}" << endl;
}
}
void choleskyFactorization(int p, float *m)
{
lapack_int q = LAPACKE_spotrf(LAPACK_ROW_MAJOR, 'L', p, m, p);
}
void copy(int p, float *from, float *to)
{
mkl_somatcopy('R', 'N', p, p, 1.0, from, p, to, p);
}
void samplesGenerate(int n, int p, float *mean, float *chol, float *res)
{
VSLStreamStatePtr stream;
int i = (int) rand();
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, n, res, p, VSL_MATRIX_STORAGE_FULL, mean, chol);
}
void wGenerate(int n, float* w)
{
float *mean = (float *) calloc(n, sizeof(float));
float *cov = (float*) calloc(n * n, sizeof(float));
for (int i = 1; i < n; i++)
{
mean[i] = 1;
cov[i*n + i] = 1;
}
int i = (int) rand();
VSLStreamStatePtr stream;
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, 1, w, n, VSL_MATRIX_STORAGE_FULL, mean, cov);
free(mean);
free(cov);
}
void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *W = (float *) calloc(n, sizeof(float));
float *X = (float *) calloc(p, sizeof(float));
wGenerate(n, W);
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, W[k], X, p, X, p, 1.0, Result, p);
}
free(W);
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}
free(X);
}
void sampleCovariance(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *X = (float *) calloc(p, sizeof(float));
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, Result, p);
}
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}
free(X);
}
void eigenDecomposition(int p, float* m, float *values) //decreasing orderб столбыц!
{
LAPACKE_ssyevd(LAPACK_ROW_MAJOR, 'V', 'U', p, m, p, values);
}
void rProjector(int p, int r, float *vectors, float* projector)
{
int num = p - r;
float *X = (float *) calloc(p, sizeof(float));
for (int i = 0; i < p; i++)
{
X[i] = vectors[p*i + num];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, projector, p);
free(X);
}
float S(int n, int p, float *m1, float *m2)
{
float *m = (float *) calloc(p * p, sizeof(float));
for (int i = 0; i < p; i++)
for (int j = 0; j < p; j++)
m[i*p + j] = m1[i*p + j] - m2[i*p + j];
float norm = LAPACKE_slange(LAPACK_ROW_MAJOR, 'F', p, p, m, p);
free(m);
return n*norm*norm;
}
int main(int argc, char **argv)
{
//задаем размерности
int paral = 16;
int p = 20;
int k = 3;
int M = 3000;
int r = 1; //нумерация с 1
int N = 3; //45
int n = 100;
int rank, size, prev, next;
//инициализация MPI
MPI_Request reqs[4];
MPI_Status stats[4];
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//
double MPI_Wtime(void);
double t1, t2, dt;
//
t1 = MPI_Wtime();
prev = rank - 1;
next = rank + 1;
if(rank==0) prev = size - 1;
if(rank==size - 1) next = 0;
//чтобы у всех сгенерировалась одинаковая матрица ковариации srand сработает одновременно
srand(time(NULL));
//файлы, в который будем писать результаты, у каждого параллельного процесса свой (номер процесса + .txt)
ofstream myfile;
char intStr[2];
if (rank < 10)
{
intStr[0] = '0';
intStr[1] = '0' + rank;
}
else
{
intStr[0] = '1';
intStr[1] = '0' + rank - 10;
}
char *txt = ".txt";
char *str = strcat(intStr,txt);
myfile.open(str);
//файлы, в которые будем отправлять размерности и тд, за это будет отвечать нулевой процесс
ofstream dim;
ofstream val;
if (rank == 0){
dim.open("dimensions.txt");
val.open("values.txt");
}
if (rank == 0)
{
dim << "p = " << p << endl;
dim << "k = " << k << endl;
dim << "M = " << M << endl;
dim << "r = " << r << endl;
dim << "N = " << N*(paral-1) << endl << endl;
dim << "n = " << n;
}
//генерация матрицы ковариации
float *COV = (float *) calloc(p * p, sizeof(float));
float *MEAN = (float *) calloc(p, sizeof(float));
float *COV_Chol = (float *) calloc(p * p, sizeof(float));
float *Basis = (float *) calloc(p * p, sizeof(float));
float *COV_L = (float *) calloc(p * p, sizeof(float));
ortBasis(p, Basis);
diagMatrix(p, k, COV_L);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, Basis, p, COV_L, p, 0.0, COV, p);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, COV, p, Basis, p, 0.0, COV_L, p);
copy(p, COV_L, COV);
copy(p, COV, COV_Chol);
choleskyFactorization(p, COV_Chol);
free(Basis);
free(COV_L);
float *values_pr = (float *) calloc(p, sizeof(float));
float *vectors_pr = (float *) calloc(p * p, sizeof(float));
copy(p, COV, vectors_pr);
eigenDecomposition(p, vectors_pr, values_pr);
//для исследования нам понадобятся собственные значения, поэтому их тоже выводит нулевой процесс
if(rank ==0){
for (int i = 0; i < p; i++)
{
val << values_pr[i] << " ";
}
val << endl << endl;
}
free(values_pr);
free(vectors_pr);
//вычислим проектор для реальной матрицы ковариации
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorR = (float *) calloc(p * p, sizeof(float));
copy(p, COV, vectors);
eigenDecomposition(p, vectors, values);
rProjector(p, r, vectors, ProjectorR);
//ProjectorR - нужный проектор
free(values);
free(vectors);
//для каждого процесса запустим свой рандом (зависящий от ранка), иначе все значения у каждого процесса будут одинаковые
srand(time(NULL) + rank);
float* NormsBootstrapWorld = (float *) calloc(N * M, sizeof(float));
for (int i = 0; i < N; i++)
{
float *Samples = (float *) calloc(n * p, sizeof(float)); //Каждый Xj - строка в матрице
for (int j = 0; j < M; j++)
{
float *COV_SampleBootstrap = (float *) calloc(p * p, sizeof(float));
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorSB = (float *) calloc(p * p, sizeof(float));
sampleCovarianceBootstrap(n, p, Samples, COV_SampleBootstrap);
copy(p, COV_SampleBootstrap, vectors);
eigenDecomposition(p, vectors, values); //вектора в столбцах
rProjector(p, r, vectors, ProjectorSB);
//ProjectorSB - нужный проектор
float s = S(n, p, ProjectorR, ProjectorSB);
NormsBootstrapWorld[i*M + j] = s;
free(COV_SampleBootstrap);
free(values);
free(vectors);
free(ProjectorSB);
}
free(Samples);
}
free(ProjectorR);
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
myfile << NormsBootstrapWorld[j*M + i] << " ";
}
myfile << endl;
}
//в myfile N столбцов
myfile << endl;
free(NormsBootstrapWorld);
free(COV);
free(MEAN);
free(COV_Chol);
t2 = MPI_Wtime();
dt = t2 - t1;
if (rank == 0)
cout << dt << endl;
MPI_Finalize(); //конец параллеливания
}