fortran 固有値

プログラミング初心者です。
http://www.netlib.org/toms/496
↑このサブルーチン（複素数からなる行列の固有値と固有ベクトルを求めるプログラム）を使いたいのですが、勉強不足のため、それすらできません。多分、配列の宣言の仕方がまずいと思うのですが…
mainプログラムはどのように書けばいいのでしょうか？
よろしくお願いします。
行列は3*3で試しています。
↓これが今のプログラムです。

complex A,B,X,EIGA,EIGB
integer N,NA,NB,NX,ITER
dimension A(N,N),B(N,N),X(N,N),EIGA(N),EIGB(N),ITER(N)
logical WANTX
N=3
NA=3
NB=3
NX=3

A(1,1)=
：
A(3,3)=

call LZHES(N,A,NA,B,NB,X,NX,WANTX)
call LZIT(N,A NA,B,NB,X,NX,ANTX,ITER,EIGA,EIGB)

end

No.4 ベストアンサー 回答者： ultraCS 回答日時： 2006/12/15 04:37

サブルーチンのソース拝見しました。

早起きしてしまった、というより昼間寝過ぎで寝られない(実は療養中)ので、飲みながらということで勘違いがあったら失礼

ここ

IF (ANORM.EQ.0.) ANORM = 1.0
IF (BNORM.EQ.0.) BNORM = 1.0
EPSB = BNORM
EPSA = ANORM
40 EPSA = EPSA/2.0
EPSB = EPSB/2.0
C = ANORM + EPSA
IF (C.GT.ANORM) GO TO 40

ここで、ANORM=EPSA=0.はあり得ないので、EPSAが浮動小数点で桁落ちするまでループが回ります。FLOATINGだとかなり時間がかかると思う(とんだと思うくらいにはかかるはず)。
このループはEPSAが0.にならない限り(桁落ちしない限り)、抜けません。BNORMがANORMより小さければ、そのときはEPSBも桁落ちするので全く意味のないループです。大体、EPSAと0.を比較すりゃすむのに、何故、Cを計算する。
どうも、ここは、間違っているような・・・

あと、お節介だけど

・WANTXが初期化されていません。多分0クリアされているので、実際は.FALSE.と思いますが、定義しておいたほうがいいでしょう。

・これはサブルーチン作者に言いたいが、D0という変数は使っちゃ駄目です。デバッグ時にD0はDOと見間違うのでね。ステートメントやよく使う関数と紛らわしい変数定義はしないように

この回答へのお礼

回答ありがとうございます。
ご指摘いただいたとおり、このループから抜け出せなくなっているようです。
しかし、わたしには固有値問題に関する知識があまりなく、そもそもこのループが何のためになるのかわからないため、どこをどう直せばいいのかわかりません。（一応、いろいろ試してみましたがうまくいきませんでした。）
また、サブルーチンLZHESの出力も明らかにおかしいので、こちらもどこかおかしいのかもしれません。
そこで、時間はかかりますが固有値問題に関する理論、fortranについて勉強しようとおもいます。

朝早くから、わざわざありがとうございました。大変参考になりました。

お礼日時：2006/12/15 12:04

No.3 回答者： ultraCS 回答日時： 2006/12/13 18:30

どこで停まるのか、少なくとも、CALLの前後にデバッグ用の出力は入れないとわからないですよ。

リンク先のライブラリも見えないし

で、もし、このソースをそのまま流しているのだとすると

call LZHES(N,A,NA,B,NB,X,NX,WANTX)
call LZIT(N,A NA,B,NB,X,NX,ANTX,ITER,EIGA,EIGB)
　　　　　　　　↑ここにカンマがないんだけど

どこのFORTRANかわからないので、こういうときにどういう動作をするかは断言しないけど(リンカで教えてくれる親切な処理系もあるけど)。LZITの中では嵐が荒れ狂っていることでしょうね。

実はカンマがあったなら、与えた引数では解けないで途方に暮れているのかも(サブルーチンが参照できないので自信なし)

あとは、つまらんアドバイスだが、可読性を上げるために

complex A,B,X,EIGA,EIGB
integer N,NA,NB,NX,ITER
dimension A(N,N),B(N,N),X(N,N),EIGA(N),EIGB(N),ITER(N)
logical WANTX
　　↓
parameter (N=3, NA=3, NB=3, NX=3)
integer ITER(N)
complex A(N,N),B(N,N),X(N,N),EIGA(N),EIGB(N)
logical WANTX

・型宣言で配列まで宣言する。
・parameter文では右辺で型が決まってしまうので、別に型宣言はしない。

など、重複しての宣言は間違いの元なので避ける(一度で宣言してしまう)
また、暗黙の型宣言(I-NがINTEGERってやつ)には、よほどのことがなければ従う

この回答への補足

回答ありがとうございます。
コンパイラは英Salford Software社の「Salford FTN77 Personal Edition」です。
リンク先ですが、見れない原因がちょっとわかりません、すみません。
The Netlib( http://www.netlib.org/ )というサイトで「Eigenvalue」で検索して一番目にhitしたものなのですが…
↓検索結果はこんな感じです。

toms/496
keywords:eigenvalue, generalized eigenvalue problem
gams: D4b4
title: LZHES/LZIT
for:generalized eigenvalue problem for complexmatrices
alg: LZ algorithm
by: L.C. Kaufman
ref: ACM TOMS 1 (1975) 271-281


いただいたアドバイスどおり型宣言を変えてみましたが、うまくいきませんでした。

また、デバッグ用の出力をつけていろいろ試したところ、サブルーチンLZIT内のif文の後で止まっているようです。
↓これです。

　IF (ABS(REAL(W))+ABS(AIMAG(W)).LE.ABS(REAL(Z))
* +ABS(AIMAG(Z))) GO TO 210

if文の前までは出力されますが、if文の後では出力されませんでした。210の後でもだめでした。

ABS(REAL(W))+ABS(AIMAG(W))と、ABS(REAL(Z))+ABS(AIMAG(Z))
の値に問題があるということでしょうか？
だとしたら、何が問題なのでしょう？

ちなみに行列Aは
A(1,1)=4
A(1,2)=2
A(1,3)=1
A(2,1)=1
A(2,2)=5
A(2,3)=2
A(3,1)=2
A(3,2)=3
A(3,3)=3
としています。行列Aの要素が複素数なので実数だけでもいいと思ったのですが…

補足日時：2006/12/13 19:23

No.2 回答者： Tacosan 回答日時： 2006/12/13 17:12

「途中でとまる」ってのは, どういう状況なんでしょうか? 何かメッセージは出ていませんか?

あと, REAL などを使うと正確に答えがでるとは限りません. というより, REAL や COMPLEX を使うとほぼ間違いなく (ほんのわずかですが) 計算に誤差が生じると思ってください. これは実数を計算機で使う以上は避けて通れない問題です. なので .EQ. を使ったところで「正確に計算できないかもしれない (例えば, 本当なら 0 になるはずが 0 にならないかもしれない) ので .EQ. を使うのは危険だよ～」ってコンパイラが言ってます. 数値計算するときには必須の知識.

この回答への補足

WARNINGの理由がよくわかりました。ありがとうございます。

止まるということについてですが、mainプログラムで下のようにしてサブルーチンLZHESの処理が終わると end lzhes というメッセージが出力されるようにしています。しかし、出力されません。
warningは出ていますが、それは全て関係式によるものです。
ERRORSも出ていません
osはwindows、コンパイラはsalford ftn77 personal　edition、エディタはCPad for Salford FTN77を使っています。


call LZHES(N,A,NA,B,NB,X,NX,WANTX)
write(6,*) 'end lzhes'

何度も申し訳ありませんが、よろしくお願いします。

補足日時：2006/12/13 17:27

No.1 回答者： Tacosan 回答日時： 2006/12/13 12:57

これがメインプログラムだったら, N は定数にしないとダメ. これだと N の値がコンパイルするときに決まらないからエラーになってると思うんだけど....

complex ... の前にでも
parameter(N=3)
と書いてみて.

この回答への補足

回答ありがとうございます。最初の部分を次のように直してみました。しかし、途中で止まってしまいました。
integer ITER,N,NA,NB,NX
parameter (N=3,NA=3,NB=3,NX=3)
complex A,B,X,EIGA,EIGB

また、関係式があるところでは必ず以下のようなWARNINGが出ます。（サブルーチン内）
なぜでしょうか？よろしくお願いします。

0062) IF (D.EQ.0.0) GO TO 80
WARNING - The use of .EQ. or .NE. with non-integer operands can produce

補足日時：2006/12/13 13:26