ガウス・ザイデル法

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57

static bool GaussSeidel(vector< vector<double> > &A, vector<double> &b, vector<double> &x, const int n, int max_iter = 1000)
{
    vector<double> x_old;
    x_old.resize(n);
 
    double denom, dev, max_dev;
    int i, j, iter = 0;
      
        for(i = 0; i < n; i++){
        denom = A[i][i];
        if(denom < RX_FEQ_EPS) return false;
 
        b[i] /= denom;
        for(j = 0; j < n; j++) A[i][j] /= denom;
    }
 
        while(true){
        for(i = 0; i < n; i++){
            x_old[i] = x[i];
            x[i] = 0;
            for(j = 0; j < n; j++){
                if(j != i){
                    x[i] -= A[i][j]*x[j];
                }
            }
            x[i] += b[i];
        }
 
                max_dev = fabs(x_old[0]-x[0])/x[0];
        for(i = 1; i < n; i++){
            dev = fabs(x_old[i]-x[i])/x[i];
            max_dev = (dev > max_dev) ? dev : max_dev;
        }
 
        if(max_dev <= RX_FEQ_EPS){
            return true;
        }
        else{
            iter++;
            if(iter > max_iter){
                return false;
            }
        }
    }
}