FORTRANの行列作成プログラム

    2003/09/08 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------

1. 概要

 実密行列を係数とする連立一次方程式Ax=bをブロックガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求める行列作成プログラム。
IDAは行列の種類でIDA=0は非対角行列が全て1で対角が1+SIG(入力値)となる。
IDA=1なら全要素を0.0〜1.0の一様乱数で作成する。
IDBは右辺行列の指定であるが、現在は解xが1.0となるようにbを作成。

2. プログラム

C=================================================================C
      SUBROUTINE DRMAT(A,B,N,ND,IDA,IDB,SIG)
C=================================================================C
C  Real-Dense Matrix Generator                                    C
C-----------------------------------------------------------------C
C    A(ND,N)  R*8, Out, Output Matrix A                           C
C    B(N)     R*8, Out, Output Right Vector B                     C
C    N        I*4, In,  Matrix Order of A                         C
C    ND       I*4, In,  Array Size of A ( ND >= N )               C
C    IDA      I*4, In,  Identify of Matrix A                      C
C                       IDA=0 ; Diag=1+Sig, Nondiag=1             C
C                       IDA=1 ; Randum number                     C
C    IDB      I*4, In,  Identify of right Vector B                C
C                       IDB=0 ; Set x=1                           C
C    SIG      R*8, In,  Value of Generation                       C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro ,   2003/08/17                    C
C        ( Hitachi Ltd. and Waseda University )                   C
C=================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(ND,N), B(N)
      REAL*4 RSIG
C  Set Matrix A
      if(IDA.eq.1) then
        RSIG = abs(mod(SIG,1.0d0))
        CALL RANDOM(RSIG)
        do j=1,N
          do i=1,N
            CALL RANDOM(RSIG)
            A(i,j) = RSIG
          end do
        end do
      else  
        do j=1,N
          do i=1,N
            A(i,j) = 1.0
          end do
          A(j,j) = 1.0 + SIG
        end do
      end if
C  Set Right Vector B
      do i=1,N
        B(i) = 0
      end do
      do j=1,N
        do i=1,N
          B(i) = B(i) + A(i,j)
        end do
      end do
C
      RETURN
      END