並列化検証
実際に扱っているソースを上げ、ここに効果のあった場所などを明記していく。
なにも並列化しないままでの実行時間は
データ量100000で
time =0.00558114
なお、OpenMPのみの効果を検証していく。
おそらく、MPIを使用した方が速くなると思われる（俺の直感）
- 時間的な都合により、効果の検証を十分に行えませんでした。
cca.cpp
#include<iostream>
#include<string>
#include<list>
#include<fstream>
#include<stdlib.h>
#include<math.h>
#include<sys/time.h>
using namespace std;
list<string> split(string,string);
double todouble(string);
int toint(string);
int cca(int, int, int, double[][1000], double[], double[][4], double, double[], int, int);
int gauss(double[][4], int, int, double);
int power(int, int, int, int, double, double[][2], double[][2], double[][4], double[],
double[][4], double[], double[], double[]);
Dataクラス
/**
* データ形式
* 1行目 変数数 属性数 データ数
* 2行目以降 第2変数の属性１ 第1変数の属性２ 第2変数の属性１ 第2変数の属性２
*/
class Data{
public:
string line; //1行目で設定を読み込む
int r,s,n; //各設定
list<string> data; //データ
string sep; //データ間の分割記号
Data(char filename[], string separater); //コンストラクタ
void init(); //初期化
int getData(double x[][1000]); //データ取得
};
Dataのコンストラクタ
Data::Data(char filename[],string separater){
ifstream fs; //ファイルストリーム
string tmp; //データ一時保存用変数
fs.open(filename); //ファイルオープン
getline(fs, line, '\n'); //lineに1行目を読み込み(設定読み込み)
//データの取り込み
while(true){
getline(fs, tmp, '\n');
//cout << tmp << endl;
data.push_back(tmp);
if(fs.eof())
break;
}
fs.close(); //ファイルオープン
sep = separater;
}
データの初期化（設定読み込み）
void Data::init(){
list<string> strs = split(line, sep); //1行目を空白で分析
list<string>::iterator it = strs.begin(); //イテレータの取得
//各設定の読み込み
r = toint(*it);
++it;
s = toint(*it);
++it;
n = toint(*it);
}
データの取得
int Data::getData(double x[][1000]){
//データ用変数の初期化
list<string>::iterator it_d = data.begin(); //1行でのデータのリスト
list<string> tmp_d; //1行のデータの分割したものを格納するリスト
list<string>::iterator it_tmp; //tmp_dを参照するためのイテレータ
//データの配列への取り込み
for(int i = 0; i < n; i++){
tmp_d = split(*it_d, sep); //1行を分割
it_tmp = tmp_d.begin(); //イテレータの取得
for(int j = 0; j < r+s; j++){
x[j][i] = todouble(*it_tmp);//配列に格納
++it_tmp;
}
/*for(int k = 0; k < s; k++){
y[k][i] =toint(*it_tmp);
++it_tmp;
}*/
++it_d;
}
return 0;
}
split関数
list<string> split(string str, string separater)
{
list<string> result; //１つ１つの文字列を持つlist
int i = 0; //切り出す文字列の先頭の番号
int j = str.find_first_of(separater); //切り出す文字列の末尾の番号
while(j != -1){
result.push_back(str.substr(i,j-i)); //文字列を切り出しlistに追加
i = j+1; //先頭の番号を更新
j = str.find_first_of(separater,i); //末尾の番号を更新
//最後尾の文字列の取り込み
if(j == -1)
result.push_back(str.substr(i));
}
return result;
}
stringをintに
int toint(string str){
return atoi(str.c_str());
}
stringをdoubleに
double todouble(string str){
return atof(str.c_str());
}
時間計測
static double elapseClock()
{
struct timeval tv;
gettimeofday(&tv, NULL);
return tv.tv_sec + (double)tv.tv_usec*1e-6;
}
メイン関数
int main() {
char filename[] = "Data.txt";
Data data(filename," ");
data.init();
double x[data.r+data.s][1000];
data.getData(x);
for(int a = 0; a < 100; a++){
for(int b = 0; b < data.r+data.s; b++)
cout << x[b][a] << " ";
cout << endl;
}
double overSt = elapseClock();
double overEd = elapseClock();
double overhead = overEd-overSt;
double ryz[data.r];
double v0[data.r];
double ab[data.r][4];
for(int i = 0; i < data.r; i++)
v0[i] = 1.0;
double before = elapseClock();
int sw = cca(data.r, data.s, data.n, x, ryz, ab, 1.0e-10, v0, 15, 2000);
double after = elapseClock();
if(sw > 0){
for(int j = 0; j < sw; j++){
cout << "相関係数 " << ryz[j] << endl;
cout << " a";
for(int k = 0; k < data.r; k++)
cout << " " << ab[j][k];
cout << endl;
cout << " b";
for(int l = 0; l < data.s; l++)
cout << " " << ab[j][data.s+l];
cout << endl;
}
}
else
cout << "***error 解を求めることができませんでした" << endl;
cout << "time =" << after-before-overhead << endl;
return 0;
}
正準相関分析
int cca(int r, int s, int n, double x[][1000], double ryz[],
double ab[][4], double eps, double v0[], int m, int ct){
int q = r + s, i1, i2, i3, sw = 0, sw1;
double vv, len;
double mean[q];
double w1[s];
double v1[r];
double v2[r];
double A[s][2];
double B[r][2];
double C[q][4];
double C11[r][r];
double C11i[r][r];
double C12[r][s];
double C21[s][r];
double C22[s][s];
double C22i[s][s];
double w[s][4];
// #pragma omp parallel for private(i2) share(i1) を使ってみたが若干速度ダウン。
// 平均値の計算（各データ列ごとに）
for (i1 = 0; i1 < q; i1++) {
mean[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
mean[i1] += x[i1][i2];
mean[i1] /= n;
}
// 分散共分散行列の計算
for (i1 = 0; i1 < q; i1++) {
for (i2 = i1; i2 < q; i2++) {
vv = 0.0;
for (i3 = 0; i3 < n; i3++)
vv += (x[i1][i3] - mean[i1]) * (x[i2][i3] - mean[i2]);
vv /= (n - 1);
C[i1][i2] = vv;
if (i1 != i2)
C[i2][i1] = vv;
}
}
/////////////若干の効果あり。要素数が多くないだけに実行したときによっては若干遅くなることも。
#pragma omp parallel
{ // C11, C12, C21, C22の設定
// C12
#pragma omp parallel for private(i2) shared(i1) nowait
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < s; i2++)
C12[i1][i2] = C[i1][i2+r];
}
// C21
#pragma omp parallel for private(i4) shared(i3) nowait
for (i3 = 0; i3 < s; i3++) {
for (i4 = 0; i4 < r; i4++)
C21[i3][i4] = C[i3+r][i4];
}
// C11とその逆行列
#pragma omp parallel for private(i6) shared(i5) nowait
for (i5 = 0; i5 < r; i5++) {
for (i6 = 0; i6 < r; i6++) {
w[i5][i6] = C[i5][i6];
C11[i5][i6] = C[i5][i6];
}
}
#pragma omp parallel for private(i8) shared(i7) nowait
for (i7 = 0; i7 < r; i7++){
for(i8 = r; i8 < 2*r; i8++)
w[i7][i8] = 0.0;
w[i7][i7+r] = 1.0;
}
}
sw1 = gauss(w, r, r, 1.0e-10);
if (sw1 == 0) {
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++)
C11i[i1][i2] = w[i1][i2+r];
}
}
else
sw = 1;
// C22とその逆行列
/////////若干だが効果あり。大枠のループの中に２つのループがあっても大丈夫のよう。
/////////ただ、変数はきっちりわけてprivateにすることが必要なよう。
#pragma omp parallel for private(i1,i2,i3)
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < s; i2++) {
w[i1][i2] = C[i1+r][i2+r];
C22[i1][i2] = C[i1+r][i2+r];
}
for (i3 = s; i3 < 2*s; i3++)
w[i1][i3] = 0.0;
w[i1][i1+s] = 1.0;
}
sw1 = gauss(w, s, s, eps);
if (sw1 == 0) {
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < s; i2++)
C22i[i1][i2] = w[i1][i2+s];
}
}
else
sw = 1;
// 固有値λ及び固有ベクトルaの計算
if (sw == 0) {
// 行列の計算
for (i1 = 0; i1 < s; i1++) {
for (i2 = 0; i2 < r; i2++) {
A[i1][i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
A[i1][i2] += C22i[i1][i3] * C21[i3][i2];
}
}
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < s; i2++) {
w[i1][i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
w[i1][i2] += C12[i1][i3] * A[i3][i2];
}
}
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++) {
A[i1][i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
A[i1][i2] += C11i[i1][i3] * w[i3][i2];
}
}
// 固有値と固有ベクトル（べき乗法）
for (i1 = 0; i1 < r; i1++) {
for (i2 = 0; i2 < r; i2++)
B[i1][i2] = 0.0;
B[i1][i1] = 1.0;
}
sw = power(0, r, m, ct, eps, A, B, C, ryz, ab, v0, v1, v2);
if (sw > 0) {
for (i1 = 0; i1 < r; i1++) {
// 相関係数
ryz[i1] = sqrt(ryz[i1]);
// 大きさの調整(a)
for (i2 = 0; i2 < r; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
w1[i2] += C11[i2][i3] * ab[i1][i3];
}
len = 0.0;
for (i2 = 0; i2 < r; i2++)
len += ab[i1][i2] * w1[i2];
len = sqrt(len);
for (i2 = 0; i2 < r; i2++)
ab[i1][i2] /= len;
// bの計算
for (i2 = 0; i2 < s; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < r; i3++)
w1[i2] += C21[i2][i3] * ab[i1][i3];
}
for (i2 = 0; i2 < s; i2++) {
ab[i1][i2+r] = 0.0;
for (i3 = 0; i3 < s; i3++)
ab[i1][i2+r] += C22i[i2][i3] * w1[i3];
}
for (i2 = 0; i2 < s; i2++)
ab[i1][i2+r] /= ryz[i1];
// 大きさの調整(b)
for (i2 = 0; i2 < s; i2++) {
w1[i2] = 0.0;
for (i3 = 0; i3 < s; i3++)
w1[i2] += C22[i2][i3] * ab[i1][i3+r];
}
len = 0.0;
for (i2 = 0; i2 < s; i2++)
len += ab[i1][i2+r] * w1[i2];
len = sqrt(len);
for (i2 = 0; i2 < s; i2++)
ab[i1][i2+r] /= len;
}
}
}
else
sw = 0;
return sw;
}
計算
int gauss(double w[][4], int n, int m, double eps){
double y1,y2;
int ind = 0, nm, m1, m2, i1, i2, i3;
nm = n + m;
for (i1 = 0; i1 < n && ind == 0; i1++) {
y1 = .0;
m1 = i1 + 1;
m2 = 0;
// ピボット要素の選択
for (i2 = i1; i2 < n; i2++) {
y2 = fabs(w[i2][i1]);
if (y1 < y2) {
y1 = y2;
m2 = i2;
}
}
// 逆行列が存在しない
if (y1 < eps)
ind = 1;
// 逆行列が存在する
else {
// 行の入れ替え
for (i2 = i1; i2 < nm; i2++) {
y1 = w[i1][i2];
w[i1][i2] = w[m2][i2];
w[m2][i2] = y1;
}
// 掃き出し操作
y1 = 1.0 / w[i1][i1];
for (i2 = m1; i2 < nm; i2++)
w[i1][i2] *= y1;
for (i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
for (i3 = m1; i3 < nm; i3++)
w[i2][i3] -= w[i2][i1] * w[i1][i3];
}
}
}
}
return(ind);
}
べき乗法による固有ベクトルの計算
int power(int i, int n, int m, int ct, double eps, double A[][2], double B[][2],
double C[][4], double r[], double v[][4], double 0[], double v1[], double v2[])
{
double l1, l2, x;
int i1, i2, i3, ind, k, k1;
// 初期設定
ind = i;
k = 0;
l1 = 0.0;
for (i1 = 0; i1 < n; i1++) {
l1 += v0[i1] * v0[i1];
v1[i1] = v0[i1];
}
l1 = sqrt(l1);
// 繰り返し計算
while (k < ct) {
// 丸め誤差の修正
if (k%m == 0) {
l2 = 0.0;
for (i1 = 0; i1 < n; i1++) {
v2[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v2[i1] += B[i1][i2] * v1[i2];
l2 += v2[i1] * v2[i1];
}
l2 = sqrt(l2);
for (i1 = 0; i1 < n; i1++)
v1[i1] = v2[i1] / l2;
}
// 次の近似
l2 = 0.0;
for (i1 = 0; i1 < n; i1++) {
v2[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v2[i1] += A[i1][i2] * v1[i2];
l2 += v2[i1] * v2[i1];
}
l2 = sqrt(l2);
for (i1 = 0; i1 < n; i1++)
v2[i1] /= l2;
// 収束判定
// 収束した場合
if (fabs((l2-l1)/l1) < eps) {
k1 = -1;
for (i1 = 0; i1 < n && k1 < 0; i1++) {
if (fabs(v2[i1]) > 0.001) {
k1 = i1;
if (v2[k1]*v1[k1] < 0.0)
l2 = -l2;
}
}
k = ct;
r[i] = l2;
for (i1 = 0; i1 < n; i1++)
v[i][i1] = v2[i1];
if (i == n-1)
ind = i + 1;
else {
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
C[i1][i2] = 0.0;
for (i3 = 0; i3 < n; i3++) {
x = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
C[i1][i2] += x * B[i3][i2];
}
}
}
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++)
B[i1][i2] = C[i1][i2];
}
for (i1 = 0; i1 < n; i1++) {
v1[i1] = 0.0;
for (i2 = 0; i2 < n; i2++)
v1[i1] += B[i1][i2] * v0[i2];
}
for (i1 = 0; i1 < n; i1++)
v0[i1] = v1[i1];
ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
}
}
// 収束しない場合
else {
for (i1 = 0; i1 < n; i1++)
v1[i1] = v2[i1];
l1 = l2;
k++;
}
}
return ind;
}
最終更新：2009年07月27日 14:15
