今日のYDL on PS3

とりあえずSPEを沢山使ってみたいなあと思い、並列計算の真骨頂である積分計算をしようと。まあ、なんだ、科学技術関係の計算なんて半分は積分計算だし。残りの半分は線形代数の計算だし。線形代数なんてLAPACKに文投げれば良いし。
そういうことで積分。積分は楽しいです。久しぶりにΓ関数の本とか読んでみたくなったぞ。

区分求積法とかとか

とりあえずサブルーチンから手を付ける感じです。あとは、あれだ、積分区間を区切るだとかは、いろいろとテクニックがあるのですよね。それはまた後でな感じで。
基本、足し算なんだけど、ベクトル計算させるとはえーらしいので、ベクトル計算させるライブラリの仕様を知ってる必要があります。ということで、

vec_madd(a, b, c) ベクタa、bの各要素を乗算し、ベクタcの各要素を加算します。
http://cell.fixstars.com/ps3linux/index.php/2.3%E3%80%80%E7%B0%A1%E5%8D%98%E3%81%AASIMD%E6%BC%94%E7%AE%97

を使います。
そして、被積分関数は数表で与えられると仮定するっぽい。

#include <stdio.h>
#include <spu_intrinsics.h>
#include <spu_mfcio.h>

//  積分計算関係のパラメータの構造体を定義します
typedef struct {
  float xInf, xSup, dx, y;
} IntegralParameter;

//  JJという積分計算のパラメータの構造体を宣言します
IntegralParameter JJ[NUM_SPE] __attribute__((aligned(16)));

//  一次精度の積分計算をさせます
float IntegralO1 (float xInf, float xSup, float dx) {
  int i = 0;
  float x = 0.0, *JJ;      //  JJは返すべき値の先頭アドレスをぶち込むっぽいです

  vector float Sum = 0.0;   //  ベクトル計算をさせるための変数を定義します
  vector float i   = 0.0;   //  ベクトル計算をさせるための変数を定義します

  for (x = xInf; x < xSup; x += dx) {
    Sum += spu_madd (YY [i], dx, Sum);  //  足し算します
    i    = (int) (x / dx);              //  連続量xを離散量iにすげーいい換言に変換します

  }

  JJ = (float *) %Sum;   //  それぞれのSPEが持ってる値の入ってるアドレスをJJに渡します


  return JJ     //  結局、結果が入っているアドレスを返します
}

これであってるかって、すげー微妙だ。動くのかな。