Cの行列作成プログラム

    2007/05/30 東京工芸大学 後 保範 ( 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. プログラム

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//================================================================C
void DRMAT(double *A, double *B, int N, int ND, int IDA, 
           int IDB, double SIG)
//================================================================C
//  Real-Dense Matrix Generator                                   C
//----------------------------------------------------------------C
//    A(ND*N)  R*8, Out, Output Matrix A                          C
//    B(N)     R*8, Out, Output Right Vector B                    C
//    N        I*4, In,  Matrix Order of A                        C
//    ND       I*4, In,  Array Size of A ( ND >= N )              C
//    IDA      I*4, In,  Identify of Matrix A                     C
//                       IDA=0 ; Diag=1+Sig, Nondiag=1            C
//                       IDA=1 ; Randum number                    C
//    IDB      I*4, In,  Identify of right Vector B               C
//                       IDB=0 ; Set x=1                          C
//    SIG      R*8, In,  Value of Generation                      C
//----------------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2007/05/30                    C
//        ( Tokyo Polytechnic University )                        C
//      Address is equal to Fortran Version                       C
//================================================================C
{ int    i, j, IRAND;
  double V; 
//  Set Matrix A
  V     = 32767.0 ;
  if(IDA == 1) {
    IRAND = fabs(SIG)*V ;
    srand(IRAND) ;
    for (j=0; j<N; j++) {
      for (i=0; i<N; i++)
       { A[i+j*ND] = rand()/V ;  }
     } }
  else {  
    for (j=0; j<N; j++) {
      for (i=0; i<N; i++)
       { A[i+j*ND] = 1.0 ; }
      A[j+j*ND] = 1.0 + SIG ;
   } }
//  Set Right Vector B
  for (i=0; i<N; i++)
   { B[i] = 0 ; }
  for (j=0; j<N; j++) {
    for (i=0; i<N; i++)
     { B[i] = B[i] + A[i+j*ND] ; }
  }
 }