2次元FDM用SOR法(C)プログラム

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

1. 概要

 2次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してSOR法で
反復解xを求める。
加速係数(ω)を与える必要がある。収束速度はωにより大きく左右される。

2. プログラム

#include <stdio.h>
#include <math.h>
// Global Define
#define NDX  301
#define NDY  301
#define A(i,j,l)  A[l-1][j][i]
#define B(i,j)    B[j][i]
#define X(i,j)    X[j][i]
extern double A[5][NDY][NDX], B[NDY][NDX], X[NDY][NDX] ;
extern FILE *FT1 ;
//=================================================================C
int SOR2D(int NX,int NY,double Omega,double EPS,int *ITER,double *ERR)
//=================================================================C
//  Solve Ax=b by SOR with 2 dimensional FDM                       C
//    Given Omega ( Acceleration factor )                          C
//-----------------------------------------------------------------C
//    A(0:NX,0:NY,5)  R*8, In, A Coefficient Matrix                C
//    B(0:NX,0:NY) R*8, In,  A Right-hand Vector(b)                C
//    NX           I*4, In,  Grid Numbers on X-axis                C
//    NY           I*4, In,  Grid Numbers on Y-axis                C
//    X(0:NX,0:NY) R*8, I/O, Initial and Solution vector           C
//    Omega        R*8, In,  SOR Acceleration factor               C
//    EPS          R*8, In,  if ||r||/||b|| <= EPS --> return      C
//    ITER         I*4, I/O, Number of Iteration                   C
//    ERR          R*8, Out, ERR=||r||/||b||                       C
//    IERR         I*4, Out, IERR=0,  Normal Return                C
//                               =1,  No Convergent                C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2003/08/21                     C
//        ( Hitachi Ltd. and Waseda University )                   C
//=================================================================C
{ int     i, j, k, IERR ;
  double  R, BN, RN, RA ; 
//  Initial
  IERR = 0 ;
  BN   = 0.0 ;
  for (j=1; j<=NY-1; j++) {
    for (i=1; i<=NX-1; i++) {
      BN = BN + B(i,j)*B(i,j) ; }
    }
//  Main Loop
  for (k=1; k<=*ITER; k++) { 
    RN = 0.0 ;
    for (j=1; j<=NY-1; j++) {
      for (i=1; i<=NX-1; i++) {
        R      = (B(i,j) - A(i,j,1)*X(i,j-1) - A(i,j,2)*X(i-1,j)
                         - A(i,j,4)*X(i+1,j) - A(i,j,5)*X(i,j+1) )
                         / A(i,j,3) - X(i,j) ;
        X(i,j) = X(i,j) + Omega*R ;
        RA     = R*A(i,j,3) ;
        RN     = RN + RA*RA ;  }
     }
//   if(ERR <= EPS) Convergent
    *ERR = sqrt(RN/BN) ;
    if(*ERR <= EPS) goto M100 ;
   }
  IERR = 1 ;
 M100: *ITER = k ;
 return (IERR) ; }