■この門くぐるべからず(その8)
QR法は行列の固有値を計算するための反復計算法です.QR分解はグラム・シュミット法やハウスホルダー法で行列の直交分解できますが,以下のサブルーチンではグラム・シュミット法を用いて,配列CZに格納された行列の非規格直交ベクトルを規格直交ベクトルに変換し,配列BZに格納します.
DSC法のプログラムでは,正規直交ベクトルを求めるのにこのサブルーチンを利用しています.小生は数値計算の専門家ではないので詳細はさておきますが,QR分解の具体的な手続きを知りたい場合や計算誤差に興味のある場合は,数値解析に関する成書が多数出版されていますので,それらを参考にしてください.
===================================
【1】QR分解のサブルーチン
11490 '
11500 ' ** ORTHOGONALIZATION ( CZ==>BZ ) **
11510 '
11520 *ORTHOGONAL:' [GRAM-SCHMIDT]
11530 FOR I=1 TO NZ
11540 FOR J=1 TO NZ
11550 H1(I,J)=0:H2(I,J)=0:H3(I,J)=0
11560 NEXT J
11570 NEXT I
11580 '
11590 S=0
11600 FOR I=1 TO NZ:S=S+CZ(I,1)*CZ(I,1):NEXT I
11610 H1(1,1)=SQR(S)
11620 FOR I=1 TO NZ
11630 H2(I,1)=CZ(I,1)/H1(1,1)
11640 NEXT I
11650 '
11660 FOR K=2 TO NZ
11670 FOR I=1 TO K-1
11680 S=0
11690 FOR J=1 TO NZ:S=S+H2(J,I)*CZ(J,K):NEXT J
11700 H1(I,K)=S
11710 NEXT I
11720 FOR I=1 TO NZ
11730 S=0
11740 FOR J=1 TO K-1:S=S+H1(J,K)*H2(I,J):NEXT J
11750 H3(I,K)=CZ(I,K)-S
11760 NEXT I
11770 S=0
11780 FOR I=1 TO NZ:S=S+H3(I,K)*H3(I,K):NEXT I
11790 H1(K,K)=SQR(S)
11800 FOR I=1 TO NZ
11810 H2(I,K)=H3(I,K)/H1(K,K)
11820 NEXT I
11830 NEXT K
11840 '
11850 FOR I=1 TO NZ
11860 FOR J=1 TO NZ:BZ(I,J)=H2(I,J):NEXT J
11870 NEXT I
11880 RETURN
===================================