C言語関係掲示板

過去ログ

No.970 複素逆行列

[戻る] [ホームページ]
No.12497

複素逆行列
投稿者---窮地(2004/02/05 18:56:22)


複素行列の取り扱いについては過去ログにもありましたが,私の疑問はまた違ったところにあるのでぜひ皆さんの力を貸してください
ソースがかなり長くなり,書き込めなかったので以下のURLにワークスペースごと圧縮したファイルをupしてますので使ってください
私は複素数をr*exp(jθ)という極座標形式で考えています.私のわからないところはmain関数内のmat->A_Amp[i][j]とmat->A_Pha[i][j]に直接数値を代入してやると逆行列が求まってもとの行列とかけると単位行列になります.
が,mat->〜にfor文で値を代入すると計算結果が違ってくることです
どこがどう間違っているのかもうさっぱりわかりません
お願いします.何とかならないでしょうか?
http://www.jupiter.sannet.ne.jp/doctor/C/index.html

No.12500

Re:複素逆行列
投稿者---たか(2004/02/05 20:11:05)


このプログラムはVC++のプロジェクト形式ですから、まずVC7.1で読み
込んで実行してみました。すると確かにおっしゃられる通りの症状が
出ました。

次に、CプログラムだけをBCC5.6.4でコンパイルして実行すると、atan2で
はDOMAIN error、powでは OVERFLOW error が発生しました。

Code Guardの出力を全部貼り付けると大変大きいので、一部を貼り付けて
みます。

Error 00001. 0x820000 (スレッド 0x0F30):
Function failure:
関数呼び出し:
atan2(-0.000000, -0.000000)=+NAN,

呼び出し履歴:
0x00403D88(=inv_mat2.exe:0x01:002D88) inv_mat2.c#364
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00402FE5(=inv_mat2.exe:0x01:001FE5) inv_mat2.c#244

------------------------------------------
Error 00002. 0x820000 (r) (スレッド 0x0F30):
Function failure:
関数呼び出し:
atan2(0.000000, 0.000000)=+NAN,

呼び出し履歴:
0x00403D88(=inv_mat2.exe:0x01:002D88) inv_mat2.c#364
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00402FE5(=inv_mat2.exe:0x01:001FE5) inv_mat2.c#244

------------------------------------------
Error 00003. 0x820000 (スレッド 0x0F30):
Function failure:
関数呼び出し:
pow(+NAN, 2.000000)=+INF,

呼び出し履歴:
0x00403D18(=inv_mat2.exe:0x01:002D18) inv_mat2.c#363
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00402FE5(=inv_mat2.exe:0x01:001FE5) inv_mat2.c#244
0x004028F4(=inv_mat2.exe:0x01:0018F4) inv_mat2.c#177

------------------------------------------
Error 00004. 0x820000 (スレッド 0x0F30):
Function failure:
関数呼び出し:
pow(+NAN, 2.000000)=+INF,

呼び出し履歴:
0x00403D3F(=inv_mat2.exe:0x01:002D3F) inv_mat2.c#363
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00403CF6(=inv_mat2.exe:0x01:002CF6) inv_mat2.c#361
0x00402FE5(=inv_mat2.exe:0x01:001FE5) inv_mat2.c#244
0x004028F4(=inv_mat2.exe:0x01:0018F4) inv_mat2.c#177

No.12501

Re:複素逆行列
投稿者---たか(2004/02/05 20:17:44)


次に、for文で代入した mat->A_Amp[][] と mat->A_Pha[i][j] の出力が
直接配列に代入したものと同一かどうか調べてみました。

--- for文版 ---
(3.000)+j(1.571) (1.000)+j(0.000) (5.385)+j(1.190) (3.606)+j(0.588) (4.123)+j(0.245) (6.403)+j(0.675) (5.000)+j(1.571)
(2.236)+j(1.107) (2.236)+j(0.464) (5.000)+j(0.927) (3.606)+j(0.983) (3.606)+j(0.588) (5.000)+j(0.644) (4.123)+j(1.326)
(2.236)+j(0.464) (3.606)+j(0.588) (5.000)+j(0.644) (4.123)+j(1.326) (3.606)+j(0.983) (3.606)+j(0.588) (3.606)+j(0.983)
(3.000)+j(0.000) (5.000)+j(0.644) (5.385)+j(0.381) (5.000)+j(1.571) (4.123)+j(1.326) (2.236)+j(0.464) (3.606)+j(0.588)
(2.236)+j(0.464) (3.606)+j(0.588) (5.000)+j(0.644) (4.123)+j(1.326) (3.606)+j(0.983) (3.606)+j(0.588) (3.606)+j(0.983)
(2.236)+j(1.107) (2.236)+j(0.464) (5.000)+j(0.927) (3.606)+j(0.983) (3.606)+j(0.588) (5.000)+j(0.644) (4.123)+j(1.326)
(3.000)+j(1.571) (1.000)+j(0.000) (5.385)+j(1.190) (3.606)+j(0.588) (4.123)+j(0.245) (6.403)+j(0.675) (5.000)+j(1.571)

--- 直接代入版 ---
(1.000)+j(0.524) (3.000)+j(0.000) (-2.000)+j(0.000) (2.000)+j(1.047) (5.000)+j(0.000) (10.000)+j(0.000) (8.000)+j(0.785)
(1.000)+j(0.000) (4.000)+j(0.000) (-4.000)+j(0.000) (2.000)+j(0.524) (4.000)+j(0.000) (-9.000)+j(0.000) (-3.000)+j(0.262)
(0.000)+j(0.000) (2.000)+j(0.349) (-3.000)+j(0.000) (2.000)+j(0.000) (1.000)+j(1.047) (0.000)+j(0.000) (0.000)+j(0.000)
(4.000)+j(0.785) (3.000)+j(0.000) (6.000)+j(0.000) (-1.000)+j(0.000) (-5.000)+j(0.000) (-3.000)+j(0.524) (13.000)+j(0.000)
(-4.000)+j(0.000) (6.000)+j(0.349) (-3.000)+j(0.000) (17.000)+j(0.000) (-7.000)+j(0.000) (-3.000)+j(0.175) (5.000)+j(0.000)
(6.000)+j(0.000) (-7.000)+j(0.698) (-2.000)+j(0.000) (0.000)+j(0.000) (3.000)+j(0.873) (9.000)+j(0.000) (6.000)+j(0.000)
(0.000)+j(0.000) (6.000)+j(0.000) (-3.000)+j(1.920) (-7.000)+j(0.000) (12.000)+j(0.000) (8.000)+j(0.349) (-7.000)+j(0.000)


No.12503

Re:複素逆行列
投稿者---たか(2004/02/05 20:24:17)


以上の結果から言える事をまとめてみたいと思います。

1.VC++またはVC.Netのpowとatan2はエラーがマスクされているか、また
 はBCCと引数、結果の扱いが異なる。OVERFLOW errorも見逃す事から、
 前者の可能性が大です。このためプログラムの間違いの発覚が困難だっ
 た。しかも実行速度が異様に遅い。コプロセッサ例外の処理にかなりの
 時間を取られている事が予想されます。

2.pow()とatan2()の引数がCodeGuardのログファイルから見て取れるよ
 うに数学的に明らかにおかしい。これはプログラムにまだバグが残って
 いる事を示しています。このバグをまず取らない限り、次の3番目の
 内容は無意味です。

3.for文で代入した配列の内容と直接代入した配列の内容が異なる。
 これでは出力結果が異なって当然です。恐らく直接代入する値が間違っ
 ていると思われます。

No.12505

Re:複素逆行列
投稿者---YuO(2004/02/05 21:43:44)


>1.VC++またはVC.Netのpowとatan2はエラーがマスクされているか、また
> はBCCと引数、結果の扱いが異なる。OVERFLOW errorも見逃す事から、
> 前者の可能性が大です。このためプログラムの間違いの発覚が困難だっ
> た。しかも実行速度が異様に遅い。コプロセッサ例外の処理にかなりの
> 時間を取られている事が予想されます。

  • 定義域エラーが起きたときは処理系定義の値を返し,errnoをEDOMに設定する。
  • 値域エラーが起きたときはHUGE_VALを返し,errnoをERANGEに設定する。

というのが標準の要求ですから,errnoを調べていない以上エラーはわからなくて当たり前です。

勝手に文字列を出力してしまうBCCのライブラリはどうかと思いますが……。


No.12510

Re:複素逆行列
投稿者---たか(2004/02/05 23:34:33)


>定義域エラーが起きたときは処理系定義の値を返し,errnoをEDOMに設定する。
>値域エラーが起きたときはHUGE_VALを返し,errnoをERANGEに設定する。

なるほど。確かにそうですね。しかし、いちいちerrnoを調べる人はそう
そういないと思います。いるとしたら数値演算ライブラリの作者ぐらい
でしょうか。BCCはそういう意味では標準の動作から外れていますね。
しかしそのお陰でたまたま今回はバグの理由がはっきりしました。

No.12511

Re:複素逆行列
投稿者---たか(2004/02/05 23:39:25)


ちなみに本題とは関係ありませんが、以下のようにBCBでは標準と動作
が異なる事を承知の上でわざと動作を変えている事がわかります。

>C++Builder の _matherr および _matherrl の動作は,UNIX 形式の
>_matherr および _matherrl の動作(メッセージを表示して終了する)
>とは異なります。UNIX 形式のルーチンを使いたい場合は,C++Builder
>のマスターディスクに入っている MATHERR.C および MATHERRL.C を使っ
>てください。

No.12523

Re:複素逆行列
投稿者---たか(2004/02/06 14:29:01)


仕事中ですのでほとんど時間が取れませんが、途中経過だけ書きます。
まず配列に値を代入したものの逆行列を取ろうとすると必ずエラーを出
して止まります。

次に配列に代入した物を極座標表示に直したものを普通の複素数として
扱い(これは無意味なのですがとにかくどんなデータでもいい)逆行列を
出して元の行列と乗算すると、単位行列になりません。これは私のミス
のはずです。

取り敢えず元の配列の行列式は0らしいです。後はランダムにデータを
生成してその逆行列を取り、ちゃんと単位行列になってきたら、極座標
表示での逆行列と行列式を求めるようにしようと思います。

ちなみに計算時間は瞬間的です。

No.12545

Re:複素逆行列
投稿者---たか(2004/02/07 13:51:53)


複素数を扱うので申し訳ありませんがC++で書いています。C++だと行数
がかなり減らせてしかも楽です。但し必要とあらばCに直せますのでご安心
を。

元のプログラムは余因子行列を用いて行列式及び逆行列を求めていまし
たが、ここでは簡単のためガウス・ジョルダン法を使用しています。

結論から申しますと、頂いたデータは複素数型と極座標型の行列データ
の両方の行列式が著しく小さい、もしくは0のため、計算途中でQuiet
NANやINFを発生し続行不能になっています。

このプログラムは読み込んだデータを無視して乱数の複素行列を発生し
ていますが、逆行列及び単位行列がちゃんと求まります。乱数ですから
ごくまれに行列式が0になるとやはり計算不能になりますがその確率は
低いです。

#include <iostream>
#include <fstream>
#include <vector>
#include <complex>
#include <cstdlib>
#include <iomanip>
#include <ctime>

#ifdef __BORLANDC__
#pragma option -w-8012
#endif

typedef std::complex<double> Cd;

const double EPS = 1e-14; // これより小さい行列要素は無視する(現在行列式との比較のみ使用)

class Cmatrix {
  static bool first;
  std::vector<std::vector<Cd> > cp;
public:
  explicit Cmatrix(const char* filename = 0);
  Cmatrix(const Cmatrix& p) : cp(p.cp) {}
  Cmatrix& operator=(const Cmatrix& p) {
    if (this != &p)
      cp = p.cp;
    return *this;
  }
  bool read(const char* filename);
  void write(const char* filename = 0, const char* comment = 0) const; // 二回目以降の呼び出しは追加モード
  Cd inverse(); // 返り値は行列式
  Cmatrix operator*(const Cmatrix&) const;
  void randomize(); // 乱数の複素配列を代入する
};

Cmatrix::Cmatrix(const char* filename)
{
  if (filename)
    if (!read(filename)) {
      std::cerr << "ファイル '" << filename << "' が開けません" << std::endl;
      std::exit(1);
    }
}

bool Cmatrix::first = true;

bool Cmatrix::read(const char* filename)
{
  int m, n;
  std::ifstream ifs(filename);
  Cd tmp;

  if (!ifs) return false;

  ifs >> m >> n;
  if (m != n) {
    std::cerr << "正方行列のデータではありません" << std::endl;
    std::exit(1);
  }
  cp.resize(m);
  for (int i = 0; i < cp.size(); i++) {
    cp[i].resize(n);
    for (int j = 0; j < cp[i].size(); j++) {
      ifs >> tmp;
#ifdef __BORLANDC__ // BCC5.6〜のバグを回避
#if (__BORLANDC__ >= 0x560)
      cp[i][j] = Cd(_STL::abs(tmp), std::arg(tmp));
#endif
#else
      cp[i][j] = Cd(std::abs(tmp), std::arg(tmp));
#endif
    }
  }
  return true;
}

void Cmatrix::write(const char* filename, const char* comment) const
{
  std::ostream* os;

  if (filename)
    os = new std::ofstream(filename, (first) ? std::ios::out : (std::ios::out | std::ios::app));
  else
    os = &std::cout;
  
  if (comment)
    (*os) << comment << std::endl;

  for (int i = 0; i < cp.size(); i++) {
    for (int j = 0; j < cp[i].size(); j++)
      (*os) << std::scientific << std::setprecision(6) << cp[i][j] << ' ';
    (*os) << std::endl;
  }
  if (os != &std::cout) delete os;

  first = false;
}

Cd Cmatrix::inverse()
{
  Cd det = 1.;

  for (int k = 0; k < cp.size(); k++) {
    Cd p = cp[k][k];
    det *= p;
    for (int i = 0; i < cp.size(); i++)
      cp[k][i] /= p;
    cp[k][k] = Cd(1.) / p;
    for (int j = 0; j < cp.size(); j++)
      if (j != k) {
         Cd d = cp[j][k];
         for (int i = 0; i < cp.size(); i++)
           if (i != k)
             cp[j][i] -= cp[k][i] * d;
           else 
             cp[j][i] = -d / p;
      }
  }
  return det;
}

Cmatrix Cmatrix::operator*(const Cmatrix& cp2) const
{
  if (cp.size() != cp2.cp.size()) std::exit(1);

  Cmatrix res(*this);
  Cd s;
  
  for (int i = 0; i < cp.size(); i++) {
    for (int j = 0; j < cp[i].size(); j++) {
      s = 0.;
      for (int k = 0; k < cp.size(); k++)
        s += cp[i][k] * cp2.cp[k][j];
      res.cp[i][j] = s;
    }
  }

  return res;
}

void Cmatrix::randomize()
{
  std::srand(static_cast<unsigned>(std::time(0)));

  for (int i = 0; i < cp.size(); i++)
    for (int j = 0; j < cp.size(); j++)
      cp[i][j] = Cd((static_cast<double>(std::rand()) / RAND_MAX) * 100 - 50,
                    (static_cast<double>(std::rand()) / RAND_MAX) * 100 - 50);
}

int main()
{
  Cmatrix cm("complex1.txt");

  cm.randomize(); // 読み込んだデータを無視して乱数の複素データを代入する
  
  cm.write("result.txt", "複素行列");
  
  Cmatrix icm(cm);
  Cd t = icm.inverse();
  std::cout << "det = " << t << std::endl;
  
#ifdef __BORLANDC__
#if (__BORLANDC__ >= 0x560)
  if (_STL::abs(t) > EPS) // 単に複素数の大きさで正則かどうかを判断
#endif
#else
  if (std::abs(t) > EPS)
#endif
    icm.write("result.txt", "複素逆行列");
  else {
    std::cout << "逆行列の作成に失敗しました。" << std::endl;
    std::exit(1);
  }

  Cmatrix res(cm * icm);
  res.write("result.txt", "単位行列");
}


No.12546

Re:複素逆行列
投稿者---たか(2004/02/07 13:53:15)


あ、このプログラムだと#ifディレクティブがおかしくてBCC5.5.1などを
使用している人はうまくコンパイルできません。#ifを取って_STL::で
なくstd::の方だけ残して下さい。

No.12547

それと
投稿者---たか(2004/02/07 14:05:26)


このプログラムは通常の複素数型しか扱えませんので極座標形式の複素数
は正しく計算できません。

極座標形式用のcomplexクラスを新たに作るか、データ読み込み時に
std::polarを使って通常の複素数に変換して計算する必要があります。
後者の方が楽でしょう。

No.12557

Re:それと
投稿者---窮地(2004/02/07 22:16:48)


複素逆行列は適応アルゴリズムの中で使用するので,パソコンではまわしません.SunOSでまわすのですがC++コンパイラが残念ながら入っていないのでなんとしてもCで書かないといけないんです.

No.12558

Re:それと
投稿者---たか(2004/02/07 22:59:41)


>複素逆行列は適応アルゴリズムの中で使用するので,パソコンではまわしません.SunOSでまわすのですがC++コンパイラが残念ながら入っていないのでなんとしてもCで書かないといけないんです.

わかりました。でも一番最初に頂いたデータではCでも動きません。
それだけはご了承下さい。ピボット選択すれば動くかもしれませんので
もう少し試してみます。

No.12559

もう一つファイルをアップします
投稿者---窮地(2004/02/07 23:03:20)


mat->A_Amp,mat->A_Phaに代入される値はX_ReとX_Imから計算さした値で,よってmat->に代入される値はfor文で代入する値も,直接代入する値も一緒です.

ただしmat->A_Ampに関してはfor文を使って値を代入しても乗算すれば単位行列になります.
mat->A_Phaはやっぱり直接値を代入しないといけないみたいです

それと行列式の値が0に近いから単位行列が出てこないというのは,直接値を代入すれば単位行列になることを考えれば,筋の通らない話ではないでしょうか・・・
http://www.jupiter.sannet.ne.jp/doctor/C/index.html

それと私はあまりC言語には詳しくないのであまり難しいことを書かれても理解できないんです(笑) 残念ながら・・・

No.12561

Re:もう一つファイルをアップします
投稿者---たか(2004/02/08 00:00:45)


余因子行列を使って逆行列を求める時は行列式が分母に入りますよね?
それで、行列式が0、または≒0の時はOVERFLOWになるのだと思います。
手書きなら計算できるかもしれませんが、コプロセッサの限界なのだと
思います。(IEEE754)

そして例えVC++で単位行列が得られても、BCB5.6.4でエラーが出て落ちる
という事は潜在的なエラーが残っている事を意味しています。このまま
ではSunOSで動かしても正常に動くという保証はありません。

しかもあの遅さは異様です。通常9×9程度の逆行列及び乗算なら今の
コンピュータではほぼ瞬間的に計算が終わるはずです。どこかで異常が
起きている可能性が大きいのです。

従ってVC++で単位行列が求まった事実はSunOSで動く事を保証するもので
は今の所ないのです。

No.12562

Re:もう一つファイルをアップします
投稿者---たか(2004/02/08 00:08:59)


>ただしmat->A_Ampに関してはfor文を使って値を代入しても乗算すれば単位行列になります.
>mat->A_Phaはやっぱり直接値を代入しないといけないみたいです

これはもしかしたら添え字がずれているか、または演算時にInf及びQNaN
を出すタイミングが違うため、結果が違う事が考えられます。その箇所を
特定する事はとても困難です。ですから行列式の値が極端に小さい時は
逆行列の作成は諦めた方が適切だと言えます。

コンピュータから見れば同じ計算をすれば必ず同じ結果が得られるはず
です。結果が違うという事は計算の仕方が違うという事です。

No.12563

実行結果
投稿者---たか(2004/02/08 00:25:41)


VC++とBCC5.6.4の実行結果が異なるので、私はBCC5.6.4の実行結果を
書く事にします。ついでにMinGW(gcc3.3.1)の結果も書きます。

新しくいただいたデータのBCC5.6.4の結果は約20秒の計算時間の後、
単位行列を出力し、その後DOMAIN error、OVERFLOW errorが40個ほど
出ます。(WinXP、Pen4 2.8G)

MinGWでは
C:/MinGW/inv_matx1.c:489: error: conflicting types for `inverse_matrix'
C:/MinGW/inv_matx1.c:7: error: previous declaration of `inverse_matrix'
C:/MinGW/inv_matx1.c:578: error: conflicting types for `determinant'
C:/MinGW/inv_matx1.c:8: error: previous declaration of `determinant'

と4つのエラーが出てコンパイルできませんでした。
ああこんな時間に。また明日。

No.12565

Re:実行結果
投稿者---窮地(2004/02/08 00:47:52)


確かに計算が異様に遅いのはおかしいとは感じていましたが,VCとBCCで処理系が異なったら結果も違うというのはおかしいですね

話はがらりと変わりますが
C言語に大変お詳しいのでそういった関係のご職業でもなさっているんでしょうか?

No.12576

Re:実行結果
投稿者---たか(2004/02/09 12:16:27)


>確かに計算が異様に遅いのはおかしいとは感じていましたが,VCとBCCで処理系が異なったら結果も違うというのはおかしいですね

一番最初に頂いたプログラムに試しに乱数を突っ込んで計算してみたら、
正しい逆行列と単位行列が結果として得られました。実行時エラーもあり
ませんでした。という事は余因子行列を求める過程でどこかすごく無駄な
事をしているはずです。但しCodeGuardをかけてみたら、どちらもメモリ
リークはなく、正しく計算されているようです。

私は以前に余因子行列でやはり逆行列を求めるプログラムを書いた事が
ありますが、実行速度はガウス・ジョルダン消去法と大差なかった事を
覚えています。

これで、最初にいただいたデータはどちらも行列式が≒0で元々逆行列
を求めるのが困難あるいは不可能である事がわかりました。一番最後に
頂いたデータはやはり単位行列を出力した後DOMAIN errorやOVERFLOW error
を吐きます。従ってプログラムの実行速度が遅い事を除けば、一番最初
のプログラムを使っていただいても差し支えない事になります。実行時
エラーが出る時は入力したデータが不適切、と言えると思います。

>話はがらりと変わりますが
>C言語に大変お詳しいのでそういった関係のご職業でもなさっているんでしょうか?

私もそうですが、ここで多量の質問に素早く答えられる方は皆その方面か
近い方面で活躍されている方々だと思います。

No.12578

Re:実行結果
投稿者---窮地(2004/02/09 15:25:49)



>一番最初に頂いたプログラムに試しに乱数を突っ込んで計算してみたら、
>正しい逆行列と単位行列が結果として得られました。実行時エラーもあり
>ませんでした。という事は余因子行列を求める過程でどこかすごく無駄な
>事をしているはずです。但しCodeGuardをかけてみたら、どちらもメモリ
>リークはなく、正しく計算されているようです。
>これで、最初にいただいたデータはどちらも行列式が≒0で元々逆行列
>を求めるのが困難あるいは不可能である事がわかりました。

一番最初のプログラム(ワークスペースごとUPしてたやつ)ですが,位相が全部0になってたと思うんですがその行列式を計算してやると(-1847107.000)+j(0.000)となりました
>一番最後に頂いたデータはやはり単位行列を出力した後DOMAINerrorOVERFLOW error
>を吐きます。
一番最後のプログラムは行列式がほぼゼロになっていました

>従ってプログラムの実行速度が遅い事を除けば、一番最初
>のプログラムを使っていただいても差し支えない事になります。実行時
>エラーが出る時は入力したデータが不適切、と言えると思います。

HPにアップしている逆行列を求めるプログラムはすべて同じで異なるのは
mat->に代入する値だけなんですが・・・

もう私には手におえなくなってきてます(悲)


No.12579

Re:実行結果
投稿者---窮地(2004/02/09 15:31:24)


逆行列を求めるには余因子行列から求めるよりガウス・ジョルダン消去法を使ったほうが無難なんでしょうか?

ガウス・ジョルダン消去法はピボットが0になったりしたら計算できなくなるのでそれを心配してるのですが・・・
もちろん行を入れ替えればすむことなんですけど一介の学生にそこまでのプログラムを作るのは難しいような気もします


No.12580

Re:実行結果
投稿者---たか(2004/02/09 15:34:55)


>ガウス・ジョルダン消去法はピボットが0になったりしたら計算できなくなるのでそれを心配してるのですが・・・
>もちろん行を入れ替えればすむことなんですけど一介の学生にそこまでのプログラムを作るのは難しいような気もします

それなら部分ピボット選択を行って分母が0にならないようにするだけで
かなりの確率で0除算は回避できます。しかも余因子行列を使うプログラ
ムならいいかというと分母に行列式が入るので問題は消えません。

No.12582

Re:実行結果
投稿者---たか(2004/02/09 15:36:13)


とにかく一度、データに乱数を入れて、エラー無しにちゃんと正しい答え
が得られる事を確認してみて下さい。その後でどうすればいいかを考える
事にしましょう。

No.12583

Re:実行結果
投稿者---窮地(2004/02/09 15:38:24)


>とにかく一度、データに乱数を入れて、エラー無しにちゃんと正しい答え
>が得られる事を確認してみて下さい。その後でどうすればいいかを考える
>事にしましょう。

for文中に乱数を発生させてmat->に代入さしてやればいいんでしょうか

No.12584

Re:実行結果
投稿者---たか(2004/02/09 15:44:46)


>for文中に乱数を発生させてmat->に代入さしてやればいいんでしょうか

X_Re、X_Imに乱数を代入した方が安全だと思います。そこから極座標形式
に変換してmatに代入しているわけですから。

No.12585

Re:実行結果
投稿者---窮地(2004/02/09 15:56:40)


>>for文中に乱数を発生させてmat->に代入さしてやればいいんでしょうか
>
>X_Re、X_Imに乱数を代入した方が安全だと思います。そこから極座標形式
>に変換してmatに代入しているわけですから。

確かに実部と虚部に乱数を発生させて見るとちゃんと単位行列になりました
ということは関数inverse_matrixとdeterminantはアルゴリズムは間違ってないということでしょうか


No.12586

Re:実行結果
投稿者---たか(2004/02/09 16:27:01)


>確かに実部と虚部に乱数を発生させて見るとちゃんと単位行列になりました
>ということは関数inverse_matrixとdeterminantはアルゴリズムは間違ってないということでしょうか

間違ってないと断定はできません(あの遅さから)が、一応正しい答えは
出せると思います。ちゃんと逆行列及び行列式は出るわけですから。

まず行列式だけ求めておき、その値が(複素数なら大きさでよいでしょ
う)0に近かったら計算を回避するのが良いと思います。

No.12602

Re:実行結果
投稿者---たか(2004/02/09 21:55:54)


余因子行列を使ったプログラムを作ってみました。これだとガウス・ジョ
ルダン法ではエラーが出る場合でも行列式及び逆行列を求められる事が
あります。最初に頂いたデータでもぎりぎりで結果が求まりました。

計算速度はガウス・ジョルダン法より少し遅いですが1秒はかかりませ
ん。

もしこのプログラムを使ってみたいとおっしゃられるならC表現に直し
ます。(これが大変なんですが・・・・)

//
// 余因子行列を使って行列式及び逆行列を求める
//

#include <iostream>
#include <fstream>
#include <vector>
#include <complex>
#include <cstdlib>
#include <iomanip>
#include <ctime>

#ifdef __BORLANDC__
#pragma option -w-8012
#endif

typedef std::complex<double> Cd;

const double EPS = 1e-14; // これより小さい行列要素は無視する(現在行列式との比較のみ使用)

class Cmatrix {
  static bool first;
  std::vector<std::vector<Cd> > cp;
public:
  explicit Cmatrix(const char* filename = 0);
  explicit Cmatrix(int m, int n) {
    cp.resize(m);
    for (int i = 0; i < cp.size(); i++)
      cp[i].resize(n);
  }
  Cmatrix(const Cmatrix& p) : cp(p.cp) {}
  Cmatrix& operator=(const Cmatrix& p) {
    if (this != &p)
      cp = p.cp;
    return *this;
  }
  bool read(const char* filename);
  void write(const char* filename = 0, const char* comment = 0) const; // 二回目以降の呼び出しは追加モード
  Cd det() const;
  Cd inverse();
  Cd Cmatrix::cofactor(int i, int j) const;
  Cd Cmatrix::minorDet(int i, int j) const;
  Cmatrix operator*(const Cmatrix&) const;
  void randomize(); // 乱数の複素配列を代入する
};

Cmatrix::Cmatrix(const char* filename)
{
  if (filename)
    if (!read(filename)) {
      std::cerr << "ファイル '" << filename << "' が開けません" << std::endl;
      std::exit(1);
    }
}

bool Cmatrix::first = true;

bool Cmatrix::read(const char* filename)
{
  int m, n;
  std::ifstream ifs(filename);
  Cd tmp;

  if (!ifs) return false;

  ifs >> m >> n;
  if (m != n) {
    std::cerr << "正方行列のデータではありません" << std::endl;
    std::exit(1);
  }
  cp.resize(m);
  for (int i = 0; i < cp.size(); i++) {
    cp[i].resize(n);
    for (int j = 0; j < cp[i].size(); j++) {
      ifs >> tmp;
#ifdef __BORLANDC__ // BCC5.6〜のバグを回避
#if (__BORLANDC__ >= 0x560)
      cp[i][j] = Cd(_STL::abs(tmp), std::arg(tmp));
#endif
#else
      cp[i][j] = Cd(std::abs(tmp), std::arg(tmp));
#endif
    }
  }
  return true;
}

void Cmatrix::write(const char* filename, const char* comment) const
{
  std::ostream* os;

  if (filename)
    os = new std::ofstream(filename, (first) ? std::ios::out : (std::ios::out | std::ios::app));
  else
    os = &std::cout;
  
  if (comment)
    (*os) << comment << std::endl;

  for (int i = 0; i < cp.size(); i++) {
    for (int j = 0; j < cp[i].size(); j++)
      (*os) << std::scientific << std::setprecision(6) << cp[i][j] << ' ';
    (*os) << std::endl;
  }
  if (os != &std::cout) delete os;

  first = false;
}

Cd Cmatrix::det() const
{
  Cd value;
  int j;

  // 余因子展開による行列式の計算
  value = 0;

  if (cp[0].size() == 1)
    value = cp[0][0];
  else 
    for (j = 0; j < cp[0].size(); j++)
      value += cp[0][j] * cofactor(0, j);

  return value;
}

Cd Cmatrix::inverse()
{
  Cmatrix invp(cp.size(), cp[0].size());    // 解となる逆行列
  Cd determinant;  // 行列式の値

  determinant = det();

  // 余因子展開による逆行列の計算
  for (int i = 0; i < invp.cp[0].size(); i++)
    for (int j = 0; j < invp.cp.size(); j++)
      invp.cp[j][i] = cofactor(i, j) / determinant;

  cp = invp.cp;

  return determinant;
}

Cd Cmatrix::cofactor(int i, int j) const
{
  return ((i + j) % 2 == 0) ? minorDet(i, j) : -minorDet(i, j);
}

Cd Cmatrix::minorDet(int i, int j) const
{
  Cmatrix md(cp.size() - 1, cp[0].size() - 1);
  int k, l, kk, ll;

  for (k = 0; k < cp.size() - 1; k++) {
    kk = (k < i) ? k : k + 1;
    for (l = 0; l < cp[0].size() - 1; l++) {
      ll = (l < j) ? l : l + 1;
      md.cp[k][l] = cp[kk][ll];
    }
  }
  return md.det();
}

Cmatrix Cmatrix::operator*(const Cmatrix& cp2) const
{
  Cmatrix res(*this);
  Cd s;
  
  for (int i = 0; i < cp.size(); i++) {
    for (int j = 0; j < cp[i].size(); j++) {
      s = 0.;
      for (int k = 0; k < cp.size(); k++)
        s += cp[i][k] * cp2.cp[k][j];
      res.cp[i][j] = s;
    }
  }

  return res;
}

void Cmatrix::randomize()
{
  std::srand(static_cast<unsigned>(std::time(0)));

  for (int i = 0; i < cp.size(); i++)
    for (int j = 0; j < cp.size(); j++)
      cp[i][j] = Cd((static_cast<double>(std::rand()) / RAND_MAX) * 100 - 50,
                    (static_cast<double>(std::rand()) / RAND_MAX) * 100 - 50);
}

int main()
{
  Cmatrix cm("complex1.txt");

//  cm.randomize(); // 読み込んだデータを無視して乱数の複素データを代入する
  
  cm.write("result.txt", "複素行列");
  
  Cmatrix icm(cm);
  Cd t = icm.inverse();
  std::cout << "det = " << t << std::endl;

#ifdef __BORLANDC__
#if (__BORLANDC__ >= 0x560)
  if (_STL::abs(t) > EPS) // 単に複素数の大きさで正則かどうかを判断
#endif
#else
  if (std::abs(t) > EPS)
#endif
    icm.write("result.txt", "複素逆行列");
  else {
    std::cout << "逆行列の作成に失敗しました。" << std::endl;
    std::exit(1);
  }

  Cmatrix res(cm * icm);
  res.write("result.txt", "単位行列");
}


No.12663

Re:実行結果
投稿者---窮地(2004/02/11 22:33:33)


>余因子行列を使ったプログラムを作ってみました。これだとガウス・ジョ
>ルダン法ではエラーが出る場合でも行列式及び逆行列を求められる事が
>あります。最初に頂いたデータでもぎりぎりで結果が求まりました。
>
>計算速度はガウス・ジョルダン法より少し遅いですが1秒はかかりませ
>ん。
>
>もしこのプログラムを使ってみたいとおっしゃられるならC表現に直し
>ます。(これが大変なんですが・・・・)
>
少し忙しくてHPをのぞけなかったんですが,できればC表現でお願いします.
話は戻るんですがRLS.cに乱数を発生さしてみたんですけど単位行列が求まりませんでした.
http://www.jupiter.sannet.ne.jp/doctor/C/index.html
一度実行して試してみていただけませんか?

No.12666

Re:実行結果
投稿者---たか(2004/02/11 23:51:03)


>少し忙しくてHPをのぞけなかったんですが,できればC表現でお願いします.
>話は戻るんですがRLS.cに乱数を発生さしてみたんですけど単位行列が求まりませんでした.
>http://www.jupiter.sannet.ne.jp/doctor/C/index.html
>一度実行して試してみていただけませんか?

C表現に直すと関数がどこからも見えてしまい危険になりますが、使い方
さえ間違わなければOKです。ぼちぼち修正してみます。

それから頂いたURLのプログラムですが、こちらの環境(VC7.1もあります
が敢えて使いません。BCC5.6.4で)DOMAIN error、OVERFLOW errorが出て
しまいます。配列への代入の仕方で結果が異なるという事は、余因子を
求める過程のどこかで添え字ミスか何かあるのかもしれません。
CodeGuardでは確認できませんでしたが。

No.12673

Re:実行結果
投稿者---窮地(2004/02/12 10:43:10)


>配列への代入の仕方で結果が異なるという事は、
乱数を発生させる場所によって結果が違ってきています
X_Re,X_Imに乱数を発生させてAmpとPhaに直した後にさらに行列計算をさせるからでしょうか


No.12692

Re:実行結果
投稿者---たか(2004/02/12 17:33:38)


>>配列への代入の仕方で結果が異なるという事は、
>乱数を発生させる場所によって結果が違ってきています
>X_Re,X_Imに乱数を発生させてAmpとPhaに直した後にさらに行列計算をさせるからでしょうか

余因子行列を求める場所に何らかのバグが残っているかも知れません。
しかしネストしている上、malloc()がちりばめてあり、とても読めません。

現在の行列式と逆行列のプログラムはきっぱりと捨てられる方が良いかと
思います。バグ付きのソースをいじくり回して余計な時間を消費するより
も、動くとわかっているプログラムを利用した方が結局お得です。

No.12693

Re:実行結果
投稿者---たか(2004/02/12 17:38:57)


それで、C++のソースをCにコンバートしたものを貼り付けておきますので
よろしかったらご利用下さい。正しく動く事が確認してあります。

注意していただきたいのは、このプログラムの行列式と逆行列を求める
プログラムは極座標形式は扱えませんので、呼び出す前に普通の複素数
表現に変換してから呼び出して下さい。

読み込むデータ形式は

2 2
(0,1)(2,3)
(4,5)(6,7)

のような形をしています。これはC++の複素数クラスの入出力に合わせた
ものです。

No.12694

C言語のソースです
投稿者---たか(2004/02/12 17:41:55)


//
// 余因子行列を使って行列式及び逆行列を求める(C言語版)
//

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define EPS  1e-14 // これより小さい行列要素は無視する(現在行列式との比較のみ使用)

#define ABS(x) (sqrt((x.re) * (x.re) + (x.im) * (x.im)))
#define ARG(x) (atan2((x.im), (x.re)))

int first = 1;

struct Complex {
  double re, im;
};

struct Cmatrix {
  int m, n;
  struct Complex **cp;
};

void make_matrix(int m, int n, struct Cmatrix **cmp);
void erase_matrix(struct Cmatrix *cmp);
int read_matrix(const char *filename, struct Cmatrix **cmp);
void write(const char* filename, const char* comment, struct Cmatrix *cmp);
struct Complex det(const struct Cmatrix *cmp);
struct Complex inverse(struct Cmatrix **cmp);
struct Complex cofactor(int i, int j, const struct Cmatrix *cmp);
struct Complex minorDet(int i, int j, const struct Cmatrix *cmp);
struct Complex cmpadd(const struct Complex a, const struct Complex b);
struct Complex cmpsub(const struct Complex a, const struct Complex b);
struct Complex cmpmul(const struct Complex a, const struct Complex b);
struct Complex cmpdiv(const struct Complex a, const struct Complex b);
void randomize_matrix(struct Cmatrix *cmp);
void cmatmul(const struct Cmatrix *a, const struct Cmatrix *b, struct Cmatrix *c);
void cmpcopy(struct Cmatrix *a, const struct Cmatrix *b);

void make_matrix(int m, int n, struct Cmatrix **cmp)
{
  int i;
  
  *cmp = (struct Cmatrix *)malloc(sizeof(struct Cmatrix));
  
  (*cmp)->m = m;
  (*cmp)->n = n;

  (*cmp)->cp = (struct Complex **)malloc(m * sizeof(struct Complex *));
  for (i = 0; i < m; i++)
    (*cmp)->cp[i] = (struct Complex *)malloc(n * sizeof(struct Complex));
}

void erase_matrix(struct Cmatrix *cmp)
{
  int i;
  
  for (i = 0; i < cmp->m; i++)
    free(cmp->cp[i]);
  free(cmp->cp);
  free(cmp);
}

int read_matrix(const char *filename, struct Cmatrix **cmp)
{
  int i, j, m, n;
  FILE *fi;
  struct Complex t;

  if ((fi = fopen(filename, "r")) == NULL) return 0;

  fscanf(fi, "%d %d ", &m, &n);
  
  if (m != n) {
    fprintf(stderr, "正方行列のデータではありません\n");
    exit(1);
  }
  
  make_matrix(m, n, cmp);
  
  for (i = 0; i < m; i++)
    for (j = 0; j < n; j++) {
      fscanf(fi, "%*c%lf%*c%lf%*c ", &t.re, &t.im);
      (*cmp)->cp[i][j].re = t.re; //ABS(t);
      (*cmp)->cp[i][j].im = t.im; //ARG(t);
    }
  
  fclose(fi);
  return 1;
}

void write(const char* filename, const char* comment, struct Cmatrix *cmp)
{
  int i, j;
  FILE *fo;

  if (filename)
    fo = fopen(filename, (first) ? "w" : "a");
  else
    fo = stdout;
  
  if (comment)
    fprintf(fo, "%s\n", comment);

  for (i = 0; i < cmp->m; i++) {
    for (j = 0; j < cmp->n; j++)
      fprintf(fo, "(%e,%e)", cmp->cp[i][j].re, cmp->cp[i][j].im);
    fprintf(fo, "\n");
  }
  
  fclose(fo);
  first = 0;
}

struct Complex det(const struct Cmatrix *cmp)
{
  struct Complex value = {0, 0};
  int j;

  // 余因子展開による行列式の計算
  if (cmp->n == 1)
    value = cmp->cp[0][0];
  else if (cmp->n == 2)
    value = cmpsub(cmpmul(cmp->cp[0][0], cmp->cp[1][1]), cmpmul(cmp->cp[0][1], cmp->cp[1][0]));
  else
    for (j = 0; j < cmp->n; j++)
      value = cmpadd(value, cmpmul(cmp->cp[0][j], cofactor(0, j, cmp)));

  return value;
}

struct Complex inverse(struct Cmatrix **cmp)
{
  int i, j;
  struct Cmatrix *invp;
  struct Complex determinant;         // 行列式の値
  
  make_matrix((*cmp)->m, (*cmp)->n, &invp);  // 解となる逆行列

  determinant = det(*cmp);
  
  // 余因子展開による逆行列の計算
  for (i = 0; i < invp->n; i++)
    for (j = 0; j < invp->n; j++)
      invp->cp[j][i] = cmpdiv(cofactor(i, j, *cmp), determinant);

  erase_matrix(*cmp);
  *cmp = invp;
  
  return determinant;
}

struct Complex cofactor(int i, int j, const struct Cmatrix *cmp)
{
  struct Complex t = {0, 0};
  
  return ((i + j) % 2 == 0) ? minorDet(i, j, cmp) : cmpsub(t, minorDet(i, j, cmp));
}

struct Complex minorDet(int i, int j, const struct Cmatrix *cmp)
{
  struct Cmatrix *md;
  struct Complex res;
  int k, l, kk, ll;
  
  make_matrix(cmp->m - 1, cmp->n - 1, &md);

  for (k = 0; k < md->m; k++) {
    kk = (k < i) ? k : k + 1;
    for (l = 0; l < md->n; l++) {
      ll = (l < j) ? l : l + 1;
      md->cp[k][l] = cmp->cp[kk][ll];
    }
  }
  
  res = det(md);
  erase_matrix(md);

  return res;
}

struct Complex cmpadd(const struct Complex a, const struct Complex b)
{
  struct Complex t;
  
  t.re = a.re + b.re;
  t.im = a.im + b.im;
  
  return t;
}

struct Complex cmpsub(const struct Complex a, const struct Complex b)
{
  struct Complex t;
  
  t.re = a.re - b.re;
  t.im = a.im - b.im;
  
  return t;
}

struct Complex cmpmul(const struct Complex a, const struct Complex b)
{
  struct Complex t;
  
  t.re = a.re * b.re - a.im * b.im;
  t.im = a.re * b.im + a.im * b.re;
  
  return t;
}

struct Complex cmpdiv(const struct Complex a, const struct Complex b)
{
  struct Complex t;
  
  t.re = (a.re * b.re + a.im * b.im) / (b.re * b.re + b.im * b.im);
  t.im = -(a.re * b.im - b.re * a.im) / (b.re * b.re + b.im * b.im);
  
  return t;
}


No.12695

C言語のソースです(続き)
投稿者---たか(2004/02/12 17:42:47)


void cmatmul(const struct Cmatrix *a, const struct Cmatrix *b, struct Cmatrix *c)
{
  int i, j, k;
  struct Complex s;
  
  for (i = 0; i < a->m; i++)
    for (j = 0; j < a->n; j++) {
      s.re = s.im = 0;
      for (k = 0; k < a->m; k++)
        s = cmpadd(s, cmpmul(a->cp[i][k], b->cp[k][j]));
      c->cp[i][j] = s;
    }
}

void randomize_matrix(struct Cmatrix *cmp)
{
  int i, j;
  
  srand((unsigned)time(NULL));
  for (i = 0; i < cmp->m; i++)
    for (j = 0; j < cmp->n; j++) {
      cmp->cp[i][j].re = (double)rand() / RAND_MAX * 100 - 50;
      cmp->cp[i][j].im = (double)rand() / RAND_MAX * 100 - 50;
    }
}

void cmpcopy(struct Cmatrix *a, const struct Cmatrix *b)
{
  int i, j;
  
  for (i = 0; i < b->m; i++)
    for (j = 0; j < b->n; j++)
      a->cp[i][j] = b->cp[i][j];
}

int main(void)
{
  struct Cmatrix *cm, *icm, *res;
  struct Complex t;

  read_matrix("complex.txt", &cm);
//  make_matrix(10, 10, &cm);
//  randomize_matrix(cm); // 読み込んだデータを無視して乱数の複素データを代入する

  write("result.txt", "複素行列", cm);

  make_matrix(cm->m, cm->n, &icm);
  cmpcopy(icm, cm);
  t = inverse(&icm);
  printf("det = (%e,%e)\n", t.re, t.im);

  if (ABS(t) > EPS) // 単に複素数の大きさで正則かどうかを判断
    write("result.txt", "複素逆行列", icm);
  else {
    fprintf(stderr, "逆行列の作成に失敗しました\n");
    exit(1);
  }

  make_matrix(cm->m, cm->n, &res);
  cmatmul(cm, icm, res);
  write("result.txt", "単位行列", res);
  
  erase_matrix(res);
  erase_matrix(icm);
  erase_matrix(cm);
  return 0;
}


No.12772

Re:C言語のソースです(続き)
投稿者---窮地(2004/02/15 11:16:19)


ソースいただきました
私から見れば大変難しそうなんですが現在読解中です

No.12776

Re:C言語のソースです(続き)
投稿者---たか(2004/02/15 18:06:54)


>ソースいただきました
>私から見れば大変難しそうなんですが現在読解中です

元に頂いたプログラムから見て、大きく三つの点が違います。

1.配列は普通の二次元配列ではなく動的二次元配列を採用している。
  これはC++のソースをCに直した際の副作用でもあるが、配列の
  大きさが変化しても適応できる利点がある。欠点は使用後は必ず
  解放しなければならない事、使用前に初期化make_matrix()を呼び
  出さなければならない事。(C++ではこの欠点が隠される)

2.極座標形式でなくて普通の複素数形式で扱っている。これは配列
  への代入時に変換する方法、またはcmp_add()など4つの四則演算
  を極座標形式に対応させる方法のいずれかで回避できる。

3.余因子行列の生成に再帰を用いた。これは一見するとプログラム
  が見づらく感じるが、元のプログラムの行列式と逆行列を求める
  プログラムはほとんど同じ計算の部分を2度書いていたのに対し、
  重複部分を一本化したので実は見通しが良い。

という感じです。実際に適応されてみたいデータがありましたらテキス
ト形式でUPされてみてください。使い方をmain()で示そうと思います。

No.12793

使い方を教えてください
投稿者---窮地(2004/02/17 00:20:16)


RLS.cの中にRxx_Re[][]とRxx_Im[][]がありますが,どうやって取り込めばいいんでしょうか?
また,求まった逆行列を例えばinvRxx_Re[][],invRxx_Im[][]という配列を作った場合,どうやってを代入すればいいんでしょうか?

いただいたソースは私にとってはかなり難しいです・・・
とりあえず使い方だけでも教えてください

No.12811

Re:使い方を教えてください
投稿者---たか(2004/02/17 14:49:12)


>RLS.cの中にRxx_Re[][]とRxx_Im[][]がありますが,どうやって取り込めばいいんでしょうか?
>また,求まった逆行列を例えばinvRxx_Re[][],invRxx_Im[][]という配列を作った場合,どうやってを代入すればいいんでしょうか?
>
>いただいたソースは私にとってはかなり難しいです・・・
>とりあえず使い方だけでも教えてください

例えば次のようなコードを挟めばいいでしょう。

  struct Cmatrix *Rxx, *iRxx;

  make_matrix(7, 7, &Rxx);

  for (i = 0; i < Rxx->m; i++)
    for (j = 0; j < Rxx->n; j++) {
      Rxx.cp[i][j].re = Rxx_Re[i][j] * cos(Rxx_Im[i][j]);
      Rxx.cp[i][j].im = Rxx_Re[i][j] * sin(Rxx_Im[i][j]);
    }

  make_matrix(Rxx->m, Rxx->n, &iRxx);
  cmpcopy(iRxx, Rxx);
  inverse(&iRxx); // 逆行列を求める

  for (i = 0; i < iRxx->m; i++)
    for (j = 0; j < iRxx->n; j++) {
      invRxx_Re[i][j] = ABS(iRxx.cp[i][j]);
      invRxx_Im[i][j] = ARG(iRxx.cp[i][j]);
    }
  
  erase_matrix(Rxx);
  erase_matrix(iRxx);


No.12502

Re:複素逆行列
投稿者---YuO(2004/02/05 20:20:59)


>ソースがかなり長くなり,書き込めなかったので以下のURLにワークスペースごと圧縮したファイルをupしてますので使ってください

ソースファイルだけがあれば十分です。
それよりも,ソースファイルにコメントをつけて下さい。
#ファイル名は英数字と-_.のみからなるようにしておくのが無難。

ちなみに,異常系の処理を先に書くようにすると,
インデントが少なくなることが多いです。


>私は複素数をr*exp(jθ)という極座標形式で考えています.

jってことは電気屋さんかな?
世間一般では虚数単位にiを使うので,iと書いた方が伝わりやすいですよ。


>私のわからないところはmain関数内のmat->A_Amp[i][j]とmat->A_Pha[i][j]に直接数値を代入してやると逆行列が求まってもとの行列とかけると単位行列になります.
>が,mat->〜にfor文で値を代入すると計算結果が違ってくることです

そもそも,直接代入している数値が間違っていませんか?
j3 = 3 * ejπ/2ですから,
<blcokquote>
mat->A_Amp[0][0]=1.0;
mat->A_Pha[0][0]=30.0*PI/180;
</blockquote>

<pre>mat-&gt;A_Amp[0][0] = 3.0;
mat-&gt;A_Pha[0][0] = 90.0 * PI / 180;</pre>
ではないでしょうか?


>どこがどう間違っているのかもうさっぱりわかりません

たぶん,そこら中が……。
大域変数まであるので全部読む気になりません。

せっかくVC++という優秀なデバッガを持つシステムを利用しているのだから,
ちゃんとデバッガ使って値を追いかけてみたらどうでしょう。

No.12509

Re:複素逆行列
投稿者---窮地(2004/02/05 23:26:42)


説明不足でした
mat->に直接代入している値はA_ReとA_Imから計算さして代入している値ではありません.私が勝手においた値です

for文を使ってmat->に代入している数値はA_ReとA_Imから計算さして出した値です

ですので逆行列自体が違うのはもちろんなんですがmat*imatが単位行列にならないのが不思議というか困っています
誤解がないようにもう一度2ファイルUPします
私自身もそんなにC言語に詳しい人間ではありませんので,皆さんのご指摘を理解できるかどうかわかりませんがよろしくお願いします
http://www.jupiter.sannet.ne.jp/doctor/C/index.html

No.12513

Re:複素逆行列
投稿者---たか(2004/02/05 23:48:37)


>ですので逆行列自体が違うのはもちろんなんですがmat*imatが単位行列にならないのが不思議というか困っています

いや、ですから、BCBでDOMAIN error、OVERFLOW errorが出ますから、
inverse_matrix()かdeterminant()のどちらかが間違っているか、元の
行列に逆行列が存在しないかのどちらかです。プログラムが一見正常に
終了するように見えますが、処理系によっては演算不能でアボートして
しまいます。

No.12516

Re:複素逆行列
投稿者---窮地(2004/02/05 23:55:12)


>>ですので逆行列自体が違うのはもちろんなんですがmat*imatが単位行列にならないのが不思議というか困っています
>
>いや、ですから、BCBでDOMAIN error、OVERFLOW errorが出ますから、
>inverse_matrix()かdeterminant()のどちらかが間違っているか、元の
>行列に逆行列が存在しないかのどちらかです。プログラムが一見正常に
>終了するように見えますが、処理系によっては演算不能でアボートして
>しまいます。

VC++ではエラーでないし,直接値を代入したときにはちゃんと単位行列になるのでinverse_matrix()とdeterminant()には問題ないのかなぁと思っているのですが・・・
本当に”窮地”です


No.12518

Re:複素逆行列
投稿者---たか(2004/02/05 23:58:28)


>VC++ではエラーでないし,直接値を代入したときにはちゃんと単位行列になるのでinverse_matrix()とdeterminant()には問題ないのかなぁと思っているのですが・・・
>本当に”窮地”です

しかし、処理系を変えた途端に動かなくなってしまうという事は、VC++で
は本当に”たまたま”単位行列が結果として得られているとしか思えません。

No.12512

Re:複素逆行列
投稿者---窮地(2004/02/05 23:45:17)


>>12503
pow()とatan2()の引数がCodeGuardのログファイルから見て取れるよ
 うに数学的に明らかにおかしい。これはプログラムにまだバグが残って
 いる事を示しています。このバグをまず取らない限り、次の3番目の
 内容は無意味です。

はどういう意味でしょうか?
pow()とatan2()はどこが間違ってるんでしょうか?

No.12515

Re:複素逆行列
投稿者---たか(2004/02/05 23:54:18)


atan2(-0.000000, -0.000000)=+NAN,
pow(+NAN, 2.000000)=+INF,

atan2()はQNaN、pow()は+Infを生成しています。それぞれ不定数、無限大
数を表します。

No.12519

Re:複素逆行列
投稿者---窮地(2004/02/06 00:01:05)


>atan2(-0.000000, -0.000000)=+NAN,
>pow(+NAN, 2.000000)=+INF,
>
>atan2()はQNaN、pow()は+Infを生成しています。それぞれ不定数、無限大
>数を表します。

atan2(-0.000000, -0.000000)は0になりました

No.12520

Re:複素逆行列
投稿者---たか(2004/02/06 00:02:20)


>atan2(-0.000000, -0.000000)は0になりました

errnoを確かめてみてください。
また明日書きます。


No.12514

Re:複素逆行列
投稿者---窮地(2004/02/05 23:49:43)


>>12502
電気屋さんです

j3 = 3 * ejπ/2ですから,
<blcokquote>
mat->A_Amp[0][0]=1.0;
mat->A_Pha[0][0]=30.0*PI/180;
</blockquote>

<pre>mat-&gt;A_Amp[0][0] = 3.0;
mat-&gt;A_Pha[0][0] = 90.0 * PI / 180;</pre>
ではないでしょうか?

mat-&gt;A_Ampとはどういうことなんでしょうか

No.12517

Re:複素逆行列
投稿者---たか(2004/02/05 23:55:33)


>mat-&gt;A_Ampとはどういうことなんでしょうか

mat->A_AmpをHTMLツールで変換して貼り付けたけれども変換する場所が
まずくてそのまま表示されてしまったのでしょう。

No.12868

ツリー構造が長くなったので改行します
投稿者---窮地(2004/02/20 21:47:33)


まずは,忙しい中時間を割いてプログラムを作っていただいたことに大変感謝いたします.
せっかくいただいたプログラムソースですが,私には読解不可能でした(私の力不足)
そこで,もう一度自分で作ったプログラムを改良して逆行列を求められるようにしました.
以前のように異様に時間がかかるということはありません.
ただVC++ではエラーは出ませんでしたが違う処理系(SunOS)でまわすとなにやらエラーらしきものが出るのですが,一体どこが悪いか教えていただけませんか?
できれば私の作ったプログラムをベースにバグフィックスしていただきたいのですが・・・
http://www.jupiter.sannet.ne.jp/doctor/C/index.html


No.12873

Re:ツリー構造が長くなったので改行します
投稿者---たか(2004/02/21 17:08:50)


>ただVC++ではエラーは出ませんでしたが違う処理系(SunOS)でまわすとなにやらエラーらしきものが出るのですが,一体どこが悪いか教えていただけませんか?
>できれば私の作ったプログラムをベースにバグフィックスしていただきたいのですが・・・

まず、SunOSで回した時の「なにやらエラーらしきもの」のメッセージを
教えていただけないでしょうか。デバッグにはそういう些細な情報が
大変重要な役目を果たします。

それから、頂いたプログラムをBCCでCodeGuardを掛けると、メモリリーク
がいくつか見受けられます。この手のプログラムは修正が大変困難です。
malloc()した領域がfree()されてないのです。私だったらそのバグが出た
関数をデバッグする事はせず、ばっさり切り捨てて新しく作り直します。
その方が早いからです。

一応CodeGuardが吐いたエラーログを貼り付けておきます。

No.12874

メモリリークが起きています
投稿者---たか(2004/02/21 17:10:56)


Error 00001. 0x300010 (スレッド 0x0FC8):
リソースリーク: メモリブロック (0xF0DD70) が解放されていません。

メモリブロック(0x00F0DD70) [長さ: 16 バイト] は malloc
によって確保されました。
呼び出し履歴:
0x004026F1(=inv_mat1.exe:0x01:0016F1) inv_mat1.c#266
0x00401F0C(=inv_mat1.exe:0x01:000F0C) inv_mat1.c#168
0x00401562(=inv_mat1.exe:0x01:000562) inv_mat1.c#67
0x0040B902(=inv_mat1.exe:0x01:00A902)

------------------------------------------
Error 00002. 0x300010 (スレッド 0x0FC8):
リソースリーク: メモリブロック (0xF03D58) が解放されていません。

メモリブロック(0x00F03D58) [長さ: 16 バイト] は malloc
によって確保されました。
呼び出し履歴:
0x004026F1(=inv_mat1.exe:0x01:0016F1) inv_mat1.c#266
0x00401F0C(=inv_mat1.exe:0x01:000F0C) inv_mat1.c#168
0x00401562(=inv_mat1.exe:0x01:000562) inv_mat1.c#67
0x0040B902(=inv_mat1.exe:0x01:00A902)

------------------------------------------
Error 00003. 0x300010 (スレッド 0x0FC8):
リソースリーク: メモリブロック (0xF0DCF8) が解放されていません。

メモリブロック(0x00F0DCF8) [長さ: 16 バイト] は malloc
によって確保されました。
呼び出し履歴:
0x004026F1(=inv_mat1.exe:0x01:0016F1) inv_mat1.c#266
0x00401F0C(=inv_mat1.exe:0x01:000F0C) inv_mat1.c#168
0x00401562(=inv_mat1.exe:0x01:000562) inv_mat1.c#67
0x0040B902(=inv_mat1.exe:0x01:00A902)

------------------------------------------
Error 00004. 0x300010 (スレッド 0x0FC8):
リソースリーク: メモリブロック (0xF057A0) が解放されていません。

メモリブロック(0x00F057A0) [長さ: 16 バイト] は malloc
によって確保されました。
呼び出し履歴:
0x004026F1(=inv_mat1.exe:0x01:0016F1) inv_mat1.c#266
0x00401F0C(=inv_mat1.exe:0x01:000F0C) inv_mat1.c#168
0x00401562(=inv_mat1.exe:0x01:000562) inv_mat1.c#67
0x0040B902(=inv_mat1.exe:0x01:00A902)

------------------------------------------
Error 00005. 0x300010 (スレッド 0x0FC8):
リソースリーク: メモリブロック (0xF07600) が解放されていません。

メモリブロック(0x00F07600) [長さ: 16 バイト] は malloc
によって確保されました。
呼び出し履歴:
0x004026F1(=inv_mat1.exe:0x01:0016F1) inv_mat1.c#266
0x00401F0C(=inv_mat1.exe:0x01:000F0C) inv_mat1.c#168
0x00401562(=inv_mat1.exe:0x01:000562) inv_mat1.c#67
0x0040B902(=inv_mat1.exe:0x01:00A902)

------------------------------------------
Error 00006. 0x300010 (スレッド 0x0FC8):
リソースリーク: メモリブロック (0xF0B538) が解放されていません。

メモリブロック(0x00F0B538) [長さ: 16 バイト] は malloc
によって確保されました。
呼び出し履歴:
0x004026F1(=inv_mat1.exe:0x01:0016F1) inv_mat1.c#266
0x00401F0C(=inv_mat1.exe:0x01:000F0C) inv_mat1.c#168
0x00401562(=inv_mat1.exe:0x01:000562) inv_mat1.c#67
0x0040B902(=inv_mat1.exe:0x01:00A902)

------------------------------------------
Error 00007. 0x300010 (スレッド 0x0FC8):
リソースリーク: メモリブロック (0xF09570) が解放されていません。

メモリブロック(0x00F09570) [長さ: 16 バイト] は malloc
によって確保されました。
呼び出し履歴:
0x004026F1(=inv_mat1.exe:0x01:0016F1) inv_mat1.c#266
0x00401F0C(=inv_mat1.exe:0x01:000F0C) inv_mat1.c#168
0x00401562(=inv_mat1.exe:0x01:000562) inv_mat1.c#67
0x0040B902(=inv_mat1.exe:0x01:00A902)

------------------------------------------
Error 00008. 0x300010 (スレッド 0x0FC8):
リソースリーク: メモリブロック (0xF0B574) が解放されていません。

メモリブロック(0x00F0B574) [長さ: 16 バイト] は malloc
によって確保されました。
呼び出し履歴:
0x004026F1(=inv_mat1.exe:0x01:0016F1) inv_mat1.c#266
0x00401F0C(=inv_mat1.exe:0x01:000F0C) inv_mat1.c#168
0x00401562(=inv_mat1.exe:0x01:000562) inv_mat1.c#67
0x0040B902(=inv_mat1.exe:0x01:00A902)

------------------------------------------
Error 00009. 0x300010 (スレッド 0x0FC8):
リソースリーク: メモリブロック (0xF035B4) が解放されていません。

メモリブロック(0x00F035B4) [長さ: 16 バイト] は malloc
によって確保されました。
呼び出し履歴:
0x004026F1(=inv_mat1.exe:0x01:0016F1) inv_mat1.c#266
0x00401F0C(=inv_mat1.exe:0x01:000F0C) inv_mat1.c#168
0x00401562(=inv_mat1.exe:0x01:000562) inv_mat1.c#67
0x0040B902(=inv_mat1.exe:0x01:00A902)

------------------------------------------

以下同様、大変長く続きます。

No.12875

コンパイラ警告
投稿者---たか(2004/02/21 17:11:54)


警告 W8004 inv_mat1.c 188: 'cmat' に代入した値は使われていない(関数 inverse_matrix )
警告 W8004 inv_mat1.c 257: 'next' に代入した値は使われていない(関数 determinant )

No.12877

言い忘れてました
投稿者---たか(2004/02/21 17:23:05)


私は手元にSunOSを動かせる環境がありませんので、「VC++で動いて
SunOSで動かない」場合のデバッグはできません。SunOSを持っておられ
る方の登場を待ちましょう。

それと、BCC5.6.4ではメモリリークは検出されたものの、実行結果は
正しく、メモリリークを除けば「正常に」動いている事になります。

No.12878

Re:言い忘れてました
投稿者---窮地(2004/02/21 18:17:31)


>私は手元にSunOSを動かせる環境がありませんので、「VC++で動いて
>SunOSで動かない」場合のデバッグはできません。SunOSを持っておられ
>る方の登場を待ちましょう。
>
>それと、BCC5.6.4ではメモリリークは検出されたものの、実行結果は
>正しく、メモリリークを除けば「正常に」動いている事になります。

もちろん私もSunOSのような高価なものは家にはないので学校に行かないとだめなんですが,エラーについては詳細に見てきます
あと,conplex_matrix.cは一応単位行列を出してますが,RLS.cではやはり正常には動いていません.
RLS.cをアップします.

No.12880

Re:言い忘れてました
投稿者---たか(2004/02/21 20:59:57)


>もちろん私もSunOSのような高価なものは家にはないので学校に行かないとだめなんですが,エラーについては詳細に見てきます
>あと,conplex_matrix.cは一応単位行列を出してますが,RLS.cではやはり正常には動いていません.
>RLS.cをアップします.

malloc()を使われる時は、きっちり後でfree()するようにプログラムされ
ると直ると思いますよ。CにはガベージコレクションがないのでJavaで
作られるのも一手です。

それと全体的に関数が長いです。特にmain()が長すぎます。恐らくこれの
せいで他の方がデバッグに取り組んでみようという気が失せられているよ
うな気がします。もちろん長くても全く構わないのですが、デバッグの際
に甚だ苦痛を伴います。サブルーチン化できる所はどんどんされて、関数
を短くしておくとデバッグが格段に楽になります。

RLS.cはまだDOMAIN error や OVERFLOW error が出ます。pow()やatan2()
に不適切なパラメータが渡っているんですね。

No.12884

Re:言い忘れてました
投稿者---窮地(2004/02/21 23:52:05)


>malloc()を使われる時は、きっちり後でfree()するようにプログラムされ
>ると直ると思いますよ。CにはガベージコレクションがないのでJavaで
>作られるのも一手です。
>
RLS.cで逆行列を求める際にcomplex_matrix.cの内容を使ってるんですが・・・
ですからcomplex_matrix.cのバグも直らないとだめなんでしょうね.


>RLS.cはまだDOMAIN error や OVERFLOW error が出ます。pow()やatan2()
>に不適切なパラメータが渡っているんですね。
どういうことなんでしょうか?

No.12893

Re:言い忘れてました
投稿者---たか(2004/02/22 21:06:54)


>RLS.cで逆行列を求める際にcomplex_matrix.cの内容を使ってるんですが・・・
>ですからcomplex_matrix.cのバグも直らないとだめなんでしょうね.

バグを取れば動くような気がします。

>>RLS.cはまだDOMAIN error や OVERFLOW error が出ます。pow()やatan2()
>>に不適切なパラメータが渡っているんですね。
>どういうことなんでしょうか?

与える行列の行列式が潰れているか、計算ミスがあると思われます。

No.12894

Re:言い忘れてました
投稿者---窮地(2004/02/22 22:58:18)


>>RLS.cはまだDOMAIN error や OVERFLOW error が出ます。pow()やatan2()
>>>に不適切なパラメータが渡っているんですね。
>>どういうことなんでしょうか?
>
>与える行列の行列式が潰れているか、計算ミスがあると思われます。

>>pow()やatan2()に不適切なパラメーターといわれましても文法的にはぜんぜん問題ないように思うんですが・・・?

No.12895

Re:言い忘れてました・追加
投稿者---窮地(2004/02/22 23:01:11)


pow()やatan2()に渡されるパラメーターには行列式は関係ないように思います
逆にpow()やatan2()をつかって行列を作ってるくらいですから

No.12911

Re:言い忘れてました・追加
投稿者---たか(2004/02/23 19:12:10)


>pow()やatan2()に渡されるパラメーターには行列式は関係ないように思います
>逆にpow()やatan2()をつかって行列を作ってるくらいですから

調べてみたら確かにそのようです、早とちりして申し訳ありません。

エラーが出ている行は、
317、333、334、343、364、365です。これらはpow()、atan2()でエラー
を出しています(DOMAIN、OVERFLOW)。

このエラーについては扱うパラメータが私には不明ですので、申し訳
ありませんがそちらで修正願います。

また530行でたくさんのメモリリークを起こしています。

No.12920

Re:言い忘れてました・追加
投稿者---窮地(2004/02/24 08:57:22)


>>pow()やatan2()に渡されるパラメーターには行列式は関係ないように思います
>>逆にpow()やatan2()をつかって行列を作ってるくらいですから
>
>調べてみたら確かにそのようです、早とちりして申し訳ありません。
>
>エラーが出ている行は、
>317、333、334、343、364、365です。これらはpow()、atan2()でエラー
>を出しています(DOMAIN、OVERFLOW)。
>
>このエラーについては扱うパラメータが私には不明ですので、申し訳
>ありませんがそちらで修正願います。
>
>また530行でたくさんのメモリリークを起こしています。

エラーの出ている行を教えていただきありがとうございました.
317、333、334、343、364、365行は逆行列がちゃんと計算できているという仮定の下,作っていった部分ですからおそらく逆行列が出ていない状況ですのでエラーとして出たのだと思います.
ただ530行目のメモリリークは完全なバグとしかいえないですね.
関数determinantとinverse_matrixないでのエラーは530行目だけでしょうか?

No.12921

Re:言い忘れてました・追加
投稿者---窮地(2004/02/24 09:09:27)


聞き忘れてました
530行目というと
if(det!=NULL){
の部分でしょうか?

No.12943

Re:言い忘れてました・追加
投稿者---たか(2004/02/24 21:33:18)


>聞き忘れてました
>530行目というと
>if(det!=NULL){
>の部分でしょうか?

そうですね。ここでmalloc()したdetを解放し忘れているわけです。
言い方を変えればdetを解放せずに新しくmalloc()したアドレスを上書き
しているという事です。

バグはここだけか?というご質問ですが、それはわかりません。論理的
エラーはコンパイルが通ってしまうからです。今回はたまたまCodeGuard
でメモリリークを発見しましたが、これも本来コンパイル時には分からな
いものです。メモリリークは実行時エラーで、しかも通常の実行時エラー
とは違いじわじわとメモリを浸食する質の悪いエラーです。

長時間走らせる場合はメモリリークは絶対に取り除いておきましょう。

No.12879

Re:言い忘れてました
投稿者---窮地(2004/02/21 18:36:20)


せっかくいただいたソースプログラムだったのですがプロの方が作られたものはやはり素人には難しかったです
使いこなせなくてすいません・・・

No.12881

Re:言い忘れてました
投稿者---たか(2004/02/21 21:07:25)


>せっかくいただいたソースプログラムだったのですがプロの方が作られたものはやはり素人には難しかったです
>使いこなせなくてすいません・・・

どなたでもわかりやすいように書いたつもりだったんですが仕様書無し
で書いたので難しかったのではないかと思います。どの関数がどんな働
きをするのか、把握するのはパズルを解くようなものですから。

それとあちこち手を抜いてあります。正方行列しか扱えない(行列式や
逆行列はこれでいいけど乗算ぐらいはできた方がいいですね)、行と列
のチェックを省いてある、既にある行列をmake_matrix()に入れると
メモリリークが起きる等です。それらは簡単に手を加えられますが、
プログラム強度が今回はそんなに重要ではないと思い、使う人がこの
ライブラリの使い方を知っているという前提で書きました。

Cの複素数及び行列のライブラリならいろんな所に転がっています。その
ようなものを利用されれば、説明書がきちんとしており、安全性もより
高くて私の作ったような物より使いやすい物がきっとあるでしょう。

No.12882

例えば
投稿者---たか(2004/02/21 21:17:09)


LAPACKはいかがでしょうか。科学技術計算用ライブラリとして多分今
最も多く使われているものの一つです。他にもいろいろありますが、
大抵有料だったり移植されてなかったりします。