（汎用的に）行列の積を求めるプログラム

質問

L*M実行列A , M*N実行列B , M*N実行列C
を引数として、A*BをCに代入するプログラムを作りたいと思っています。
具体的には、エラーがでるのはわかっていますがイメージとして

void multiply(int L,int M,int N,double A[L][M],double B[M][N],double C[L][N])
{
　　　　　　for(int l=0;l<L;l++)　　　　for(int n=0;n<N;n++)　　　　C[l][n]=0;

　　　　　　for(int l=0;l<L;l++){
　　　　　　　　　　　for(int m=0;m<M;m++){
　　　　　　　　　　　　　　　　　　for(int n=0;n<N;n++)　　　　　　　C[l][n]+=A[l][m]*B[m][n];
　　　　　　　　　　　}
　　　　　}
}

という感じです。Cを自分なりに学習したところ、「2次元配列を渡すには行or列数が少なくとも既知でないといけない」のは知識として得ているのですが、ではこの問題を回避して所望のプログラムを書け、と言われると行き詰ってしまいました。

卒論をするにあたり、種々のサイズの行列演算が必要なので、できるだけ効率的な関数を書いておきたいのですが、どうすれば一般の行列の積の演算が可能になるのでしょうか？

回答No.11

麻野 なぎ（@AsanoNagi） · ベストアンサー

No.2 です。

とりあえず、（インクルードフォルダだけではなく）多分同じ所にある、「ライブラリのサーチフォルダ」（のような名前の所）に、該当の lib が存在するフォルダを追加すればいいのかなと思います。

maxtrix くらいなら、ヘッダファイルだけをインクルードすれば使えると思っていたのですが、なぜに、ライブラリをリンクしようとするのか？　という気はしますが。

もしも、明示的に、libboost_serialization-vc100-mt-gd-1_44.lib をプロジェクトに加えていたのなら、このファイルをプロジェクトから外してみるのもいいかなとは思います。

お礼

ありがとうございます！
プロパティマネージャのmicrosoft.cpp.win32のプロパティ内から、リンカー→全般→追加のライブラリディレクトリに該当のlibフォルダを追加したところ解決しました。

>>ヘッダファイルだけをインクルードすれば使えると思っていたのですが、なぜに、ライブラリをリンクしようとするのか？　という気はしますが。
すいません、、、gogle先生で探した導入マニュアルページに従っていく以外の知識が無いものでして・・・

回答No.12

noname#137556

ANo.10 のお礼について

multiply 関数などの定義に合わせると，全要素分をまとめて確保しないといけません。
　　double *pa = new double[L*M];

> MATLABは知人の研究室で採用しているらしいですが、
> このソフトは1コマンド1応答？形式と聞いています。

おっしゃる通り対話式に計算させることも出来ますし，
プログラム（M ファイル）を作っておいて実行させることも出来ます。

もちろん，for や if などの制御構造も使えますし，独自の関数を定義することも可能です。
http://www.mathworks.co.jp/help/ja_JP/techdoc/ref/for.html

回答No.10

noname#137556

#4 です。

行列を作る部分も関数化してこんなのをイメージしてました。

　MATRIX a, b, c;　　　　// MATRIX 構造体
　Eye(&a, L, M);　　　　　// a を単位行列で初期化
　Eye(&b, M, N);　　　　　// b を単位行列で初期化
　Product(&c, &a, &b);　// 積の計算 c = a b

計算が目的なら，MATLAB を使えば #4 の補足の計算は，
　a = eye(L, M)
　b = eye(M, N)
　c = a * b
これだけで済みます。

分野によっては当たり前のように使われるので，覚えておいて損はないと思いますよ。

よく似た Scilab は無料で入手できるので，遊んでみてください。
http://www.scilab.org/

お礼

再度の回答ありがとうございます。
MATLABは知人の研究室で採用しているらしいですが、このソフトは1コマンド1応答？形式と聞いています。
例えばfor文のように、パラメータを更新しつつ繰り返し繰り返し演算させることはできるのでしょうか？

先の1次配列として計算させるというお言葉をヒントにして再度プログラムを書いたのですが、ビルドエラーは出ませんが異常終了してしまいました。。。申し訳ないですがご指導頂けないでしょうか？
今までポインタを使わず、なるべく配列のみで処理してきたのでこのあたりの理解が劣っているのは、恥ずかしながら自覚しておりますが。

#include<iostream>
#include<cmath>
using namespace std;

void initialize(int L,int M,double* pa);
void multiply(int L,int M,int N,double* pa,double* pb,double* pc);
void print(int L,int M,double* pa,char* ch);

int main()
{
const int L=2;
const int M=3;
const int N=4;

double** A = new double* [L];
for(int l=0;l<L;l++) A[l] = new double [M];//A[L][M]
double* pa=&A[0][0];

double** B = new double* [M];
for(int m=0;m<M;m++) B[m] = new double [N];//B[M][N]
double* pb=&B[0][0];

double** C = new double* [L];
for(int l=0;l<L;l++) C[l] = new double [N];//C[L][N]
double* pc=&C[0][0];

initialize(L,M,pa);
initialize(M,N,pb);
initialize(L,N,pc);

multiply(L,M,N,pa,pb,pc);
print(L,M,pa,"A");
print(M,N,pb,"B");
print(L,N,pc,"C");

for( int l=0; l<L; l++) delete [] A[l];
delete [] A;

for( int m=0; m<M; m++) delete [] B[m];
delete [] B;

for( int l=0; l<L; l++) delete [] C[l];
delete [] C;

return 0;
}

void initialize(int L,int M,double* pa)
{
for(int l=0;l<L;l++){
for(int m=0;m<M;m++){
//if(l!=m) pa[l*M+m]=0;
//else pa[l*M+m]=1;
pa[l*M+m]=1;
}
}
}

void multiply(int L,int M,int N,double* pa,double* pb,double* pc)
{
for(int l=0;l<L;l++) for(int m=0;m<M;m++) for(int n=0;n<N;n++) pc[l*N+n]=pa[l*M+m]*pb[m*N+n];
}

void print(int L,int M,double* pa,char *ch)
{
for(int l=0;l<L;l++){
for(int m=0;m<M;m++){
cout<<ch<<"["<<l<<"]["<<m<<"]is "<<"="<<pa[l*M+m]<<endl;
}
}
}

回答No.9

Interest

私も、ANo.7 の boostライブラリに1票。

目的は卒研用のプログラムを動かすことであって、行列演算のライブラリを作ることではありませんよね。
世の中で広く使われているライブラリがあるのですから、そちらを使ったほうが時間も労力も大幅に(数週間分以上）節約できると思います。なにしろ、自分で作るとバグも大量に作りこむことになるので。:-p)

C++やboostに抵抗があるなら、matlibという手もあります。
http://robot.cs.kobe-u.ac.jp/contents/staff/hanahara/programs/matlib.html

いずれにしても、先人が苦労して作ってくれたものがあるのに、似たようなものをまた作るというのは無駄な努力です。

回答No.8

正親町（@Ohgimachi）

>L*M実行列A , M*N実行列B , M*N実行列C
C=A*BとするとCはL*N行列になるんじゃないの？

for(i=0; i<L; i++){
for(n=0;n<N;n++)
for(m=0;m<M;m++){
C[i][n]+=A[i][m]*B[m][n];
}
}

Cにはポインターという便利なものがあるので、
配列のアドレス計算をポインターで実現すればよいのです。
したがって一次元配列のアドレスとして関数に渡すのが正解です。

Cの要素を計算するこの部分は
for(m=0;m<M;m++){
C[i][n]+=A[i][m]*B[m][n];

ポインターで書くとこうなります。
duble *a,*b,*c;

a=&A[i*M],b=&B[n];
for(m=0;m<M;m++){
*c+=(*a)*(*b);
a++,b+=N;
}
c++;

まとめると、こんな感じ

duble *a,*b,*c;
int i,n;

c=C;
for(i=0; i<L; i++){
□for(n=0;n<N;n++){
□□double *a, *b, *amax;
□□a=&A[i*M]; //aの先頭アドレスを計算する
□□amax=a+M; //aの終端アドレスを計算する
□□b=&B[n]; //bの先頭アドレスを計算する
□□for(;a<amax;a++){
□□□*c+=(*a)*(*b);
□□□b+=N;
□□}
□□c++;
□}
}

回答No.7

麻野 なぎ（@AsanoNagi）

Ｃ＋＋の範疇で解決すれば良くて、行列の計算が「道具として」必要ということなら、boost::matrix

http://www.boost.org/doc/libs/1_39_0/libs/numeric/ublas/doc/matrix.htm

がお勧めですが。
ソースも見ることができるので、ソースを見て勉強するのもありだと思います。

お礼

回答ありがとうございます

早速boostのインストールを
http://www.kmonos.net/alang/boost/install.html
を参考にしつつ行ってみましたが、いざhttp://www.page.sannet.ne.jp/d_takahashi/boost/ublas/index.html
の見本プログラムをビルトしてみると
>LINK : fatal error LNK1104: ファイル 'libboost_serialization-vc100-mt-gd-1_44.lib' を開くことができません。
と出てしまい実行ができません。
本質問の趣旨からは外れてしまうのですが、googleで検索しても全くこのエラーの情報が無いので、よろしければアドバイス頂けないでしょうか？

私の行った具体的なインストール手順は
http://www.boostpro.com/download/
からBoostPro 1.44.0 Installer を実行し、チェックボックスで選択すべき個所では
visual studio 2010を、また、default variantsは下2つ（single thread)を除き全てを選択しました。
更に、visual studio2010において
http://msdn.microsoft.com/query/dev10.query?appId=Dev10IDEF1&l=JA-JP&k=k(VS.TOOLSOPTIONSPAGES.PROJECTS.VCDIRECTORIES)&rd=true
を参考にし、VC++ ディレクトリのインクルードディレクトリにC:\Program Files(x86)\boost
を追加指定した次第です。

インストールの結果、C:Program Files(x86)にboostフォルダが作成され、上記のfatal error LNK1104:
で指摘されているlibファイルはC:\Program Files(x86)\boost\boost_1_44\libフォルダ内に存在していたのですが、上記fatal errorが現れました。

回答No.6

D-Matsu

引数の渡し方がわからないなら構造体にパックしてしまうのが楽でしょう。
構造体に行列サイズ情報を含めておけば引数ひとつで完全な行列の情報を渡すことができます。
あと、CなのかC++なのかははっきりしておいてください。
Cならnewは使えません。

回答No.5

Tacosan

行列を
double** p=new double* [L];
for(int l=0;l<L;l++)p[l]=new double [M];
のような形で確保するなら, 引数は double ** じゃないとダメだよね.

回答No.4

noyuo

関数としてもいいですけど、１次元の配列を２次元の配列と扱うための
アドレス計算のマクロを定義します。
/* 長さR の Y行　X列のアドレス求める */
#define REC_POS(R,X,Y) (Y*R+X)

として、C[REC_POS(L,n,l)] などとして参照します。
C,A,Bともに１つのマクロで使えるかなぁ？

回答No.3

noname#137556

基本は，全要素分のメモリを確保(malloc)して，要素を一列に並べておけば良いです。
どういう順番で列べるか（行優先か/列優先か）さえ決めておけば，欲しい要素がどこにあるか計算で求められますよね？

ただ，計算するのに，行数や列数といった情報が必要なので，まとめて構造体にして，
その構造体を引数で受け取って計算をする関数をどんどん定義していけば良いです。

C++ なら，クラスにして，演算子を再定義して，ベクトルや部分行列を継承でうまく表現して・・・
と妄想が膨らみますが，

行列の演算を実装すること自体が目的でなければ，実績のある既存のライブラリを使ったり，
そもそも，MATLAB や Scilab，R 言語といった行列の取り扱いに長けた処理系を使ったほうが良いかもしれません。

お礼

回答ありがとうございます
>>要素を一列に並べておけば
参考になりました。一応こうすれば所望のプログラムは、動作するという意味では形にはなりましたが・・・以下のように前処理をしないと動かないものしか作れず、不格好な仕上がりとなってしまいました。
どうにかして引数に2次元配列を組み込み、(1)直接演算するか、(2)前処理部分を1つの関数に纏めるかしたいのですが、やはり多次元配列の関数への渡し方が分かりません。

一応ですが、クラス、演算子オーバーロード、vectorの基礎は習いましたので、これらを用いて楽？ができるのでしたらそちらについてもアドバイスを頂ければありがたいです。

double A[L][M],B[M][N],C[L][N];

for(int l=0;l<L;l++) for(int m=0;m<M;m++) {if(l!=m) A[l][m]=0; else A[l][m]=1;}
for(int m=0;m<M;m++) for(int n=0;n<N;n++) {if(m!=n) B[m][n]=0; else B[m][n]=1;}

//行列の積を無理やり１次元配列を使って解く。まず前処理
double a[L*M],b[M*N],c[L*N];
for(int l=0;l<L;l++) for(int m=0;m<M;m++) a[M*l+m]=A[l][m];
for(int m=0;m<M;m++) for(int n=0;n<N;n++) b[N*m+n]=B[m][n];

multiply(L,M,N,a,b,c);//積計算
for(int l=0;l<L;l++) for(int n=0;n<N;n++) C[l][n]=c[N*l+n];//結果代入

なお、multiplyの定義は

void multiply(int L,int M,int N,double* A,double* B,double* C)//A[L][M] B[M][N] C[L][N]
{
for(int ln=0;ln<L*N;ln++) C[ln]=0;

for(int l=0;l<L;l++){
for(int m=0;m<M;m++){
for(int n=0;n<N;n++) C[N*l+n]+=A[M*l+m]*B[N*m+n];
}
}
}

回答No.2

k_kota

普通のCってそういう感じで文字を定義していいんでしたっけ？関数の最初でやっておくのが普通だと思うのですけど。

とりあえず、配列の参照渡しとポインタの概念くらいは分からないとどうにもならないと思います。
関数の定義はダメだと思います。

関数にはもろもろのサイズと2次配列をポインタ渡ししてやればいいです。
そういう意味では考え方はOKです。
あとはサイズ既知なのですから、計算結果のすべての要素について
それぞれ計算してやればいいです。
行列の積の計算方法は分かりますかね。

ループは代入先の行と列のループと計算時のループで３重になります。

まあ、探せば行列の基礎なのでそれくらいあると思いますけどね。
自分で探して調べたり、理解できないならプログラムするほどのものでは無いような気がします。
エクセルでいいんじゃないの？
もっと高度なことがしたいなら、３日程度できっちりプログラミングの基礎をやらないと後々苦労すると思います。