たびたび申し訳ありません

[徳竹]

昨日do文で処理すると上手く答えが出ないという質問に対してdo文の途中で変数の値が変わっているのではないか?
という指摘を受けたものです。あの時早とちりで解決したと言ってしまいましたが、まだ少し誤りがあるようです。

もしかして一回目のループで作成した配列の値が残っているのではないかと思い、（この行為にどの程度意味があるかわかりませんが）ループの最後で今まで作成した配列にすべて0を代入したところ次元が低い問題では正確に答えが出るのですが、次元が上がると答えが違ってしまいます。

そこで質問なのですがプログラムの途中で配列を初期化の様にするにはどの様にすればいいのでしょうか？またそもそもdo文においてそんなことをする必要があるのでしょうか？同じプログラムに二度も質問してしまい大変申し訳ありません、よろしくお願いします。環境はLinux,コンパイラはgccです。

[Mist]

> do文の途中で変数の値が変わっているのではないか? という指摘を受けたものです。

そんな指摘してません。
do内で使用する変数は2回目の実行前に初期状態に戻さないと同じ結果は得られないんじゃないですか、です。


> 途中で配列を初期化の様にするにはどの様にすればいいのでしょうか？

言いたいことがいまいちよくわかんないけど、プログラムの途中で宣言時のような書き方で初期化したいってことですか？

a[3] = { 1, 2, 3};

結論から言うと無理です。
初期値をセットしている配列を別に用意してmemcpyなどでコピーするのが比較的楽です。


> またそもそもdo文においてそんなことをする必要があるのでしょうか？

do文がどうとかいうことではなく、あなたがそのプログラムでどういうことをしてどういう結果を得たいかによります。
必要があるかどうかはあなたが組んだプログラムしだいです。
必要がないように考慮したコーディングであれば必要ないですし、考慮されていないコーディングであれば必要です。(つまりケースバイケース、必ず必要・不必要とは答えられません)
もう一度言いますが、同じ結果が出したいのであればdo～while内で使用する全ての変数について、1回目の実行開始前と2回目の実行開始前で同じ値になっていることを見直してください。

[Mist]

プログラムで具体的に書いてみる。

int a = 1;
int b = 0;
int loop = 2;

do {
    b = a + b;
    printf("b = %d\n", b);

    --loop;
} while (loop);

このプログラムを実行すると
b=1
b=2
となりますね。
これは一回目の実行によってbに1という値に更新されていためです。
同じ結果にしたいのであれば

do {
    b = a + b;
    printf("b = %d\n", b);

    b = 0;
    --loop;
} while (loop);

というように、一回目の実行が終わった時点でbを元に戻しておく必要があるわけです。
で、前回提示されたプログラムは一部分しかないため判りませんでしたが、上のプログラムのように変数の再初期化が不十分であるために、2回目の実行で一回目とは異なる計算・処理が行われることによって違う結果になるのではないですか？
というのが私からの回答です。

[徳竹]

すばやい回答ありがとうございます。なるほど変数の再初期化ですか...。
ただそこは前回の質問の後直して今も確かめてみたのですがおそらく大丈夫でだと思われます。やはり配列など何か別の所で少し狂いが生じていると思われます。

本来ならばプログラムを全部載せてどこがいけないのでしょうかとやりたい所ですが、200行近くありさすがにそこまでお願いするわけにはいかないので自力で頑張ります。わざわざありがとうございました。

[non]

200行ぐらいの短いプログラムなら、見た方がはやいけどな。
添付しちゃえば。

[徳竹]

やはり自力では苦しいので、そう言って下さるのなら載せてみます。

前にも書きましたが行列Aの最大固有値に対する固有ベクトルをarnoldi法を用いて求めるプログラムを書いています。arnoldi法とはまず初期値kと初期ベクトルrをとり、グラムシュミットの正規直交化を用いて行列Aを上ヘッセンベルグk+1×kの行列Hに変換する方法です。同時に求められるn×(k+1)の行列Qk+1から最初のk列を取ってきた行列Qkとの間にAQk=Qk+1×Hの関係が成立します。

その後行列Hの行番号と列番号が等しい要素だけ１を引きそれをHnewとして特異値分解します。求めた右特異ベクトルのk列目をｖとします。このｖをQkの右からかけたものをqとします。今度はqを初期値として同じ操作を繰り返しHnewの特異値が十分小さくなったときのqを固有ベクトルとして採用します。

長々と書いてしまいましたが、arnoldi法の部分や特異値分解の部分は正しい値が出ているので数学的なことではなくプログラミングの問題だとおもわれます。図々しいお願いだというのは承知しているのですが、良かったらお願いします。

アルゴリズムも一応書いておきます。
arnoldi法のアルゴリズム
１：初期ベクトルqを自身の２－ノルムで割ったものをq[1]とする
2：for j=1to k
3;z=Aq[j]
4;for i=1 to j
5; h[j]=q*z
6:z=z-h[j]*q
7:end for
8:h[j+1][j]=zの２－norm
9:if h[j+1][j]=0,then stop
10: q[j+1]=z/h[j+1][j]
end for

その後は配列h[k+1][k]を行列H,ベクトルq(1)からq(k+1)を列とした行列をQk+1として上記の操作を行います。
まとめると

Repeat
[Qk+1,H]=arnoldi(A,q,k)
H-I~=Uσ(Vの転置行列)　（但しI~は単位行列の一番下に全て要素が0の行を加えたもの）
v=Vのk列目
q=Qk*v
Until σmin(H-I~)<ε

N次元の行列Aは対角成分がすべて０、それ以外はランダムで0または1としそれを列和が1となるようした
行列をA’とし、
その行列とすべての要素が1/Nの行列をDとして、A=0.85*A’+0.15*Dとして行列Aを設定します。初期ベクトルはすべて要素が２。またk=4としました。

肝心の結果に関してなのですがn=4として読み込むと二回だ特異値が０となり,正しい値に収束します。しかしn=5とした場合
０．６４５４５７、
０．１４５８２１、
０．４２１３６８、
０．５９８１０３、
０．１６３８３７という値が得られます。
しかしこのプログラムを一回のループでbreakしさらに初期ベクトルをキーボードから読み込むようにかえます。そして初期ベクトル２を読み込み得られたベクトルを再び読み込ませると値が
０．６４５４３３、
０．１８５２０２、
０．３８８１１３、
０．５５７６６６、
０．２９５７９２となり値が異なります（ちなみにこれを繰り返すと特異値が収束して正しい答えが得られます）

一体何がマズイでしょうか？というのが質問内容です。一ヶ月程前に初めてC言語を始めネットや本を参考にして何とか書いた素人なので前回指摘していただいた点で見落としがあるかもしれません。
拙い文章で申し訳ありません、図々しい質問だとは思いますがよかったらお願いいたします。ちなみにプログラムの途中でclapackを用いました。

[non]

これって、別の言語で書かれていたものを、移植している？
配列はC言語ではindexが０から始まりますが、このプログラムは１から始まるように
作られています。０を使わないでこのように作るのはよくありますが・・・

私が前のスレッドで
>２　宣言した配列の範囲を超えた場所に書き込んでいる。
>　例えば、

> for(l=1;l<=n;l++){
> w[[/url][j]-=h[j]*v[[/url];}


>などで、double w[n][n];で宣言されていたとか

って書きましたよね。

徳竹さんのプログラムはちらっとしか見ていないので、見間違いかもしれませんが
（字下げをもっとわかりやすくしましょう）

5行目　#define COL 4
74行目　k=COL;　
137行目　static double vv[CO[/url]
142行目　for(i=1;i<=k;i++)
143行目　　vv=vt[-1+i*CO[/url]

つまり、vv[0]～vv[3]までしかないのに
vv[4]に書き込んだことになります。
C言語はエラーを出さないから、プログラマは注意が必要です。
static領域に取られたこの配列の前か後にある（処理系によるので）配列のデータを壊します。
ぱっと見ただけなので、ここだけしか発見しませんでしたが（っていうより１カ所見つければ良いかと）
他にもあるかもしれません。あとは、ご自分で。

[徳竹]

見づらいプログラムで申し訳ありません。確かにVV[4]に書き込んでますね。まったく気づきませんでした。普段matlabを使っているのでその時の癖が出てしまいました。
最後にひとつだけ質問したいのですがこのような事態を避けるためにあらかじめ配列の大きさに十分な余裕を持たせておく（上の例ではvv[100]などとする）ことをしていいのでしょうか？
本当は自分で確認するべきですが、今日明日と出来る環境にいないのでそれだけ質問させてください。後は自分で頑張りますありがとうございました。

[Mist]

メモリに余裕のある環境であれば許されるんじゃないでしょうか。
ただ、そういった悪い癖をつけてしまうと後々苦労しますよ。

for (i=0;i<N;i++)

こういうふうに書く癖付けたほうがいいと思いますが。

[non]

メモリさえあればかまわないと思いますが、データについては、自分で把握することが必要ですね。

[徳竹]

分かりました、やはり書き直してみます。いろいろとありがとうございました。