C..7 AL-LIN

AL-LIN3.C

/*
 * AL-LIN3.cpp
 * g++ -I/usr/local/include/profil -o AL-LIN3 AL-LIN3.cpp -lProfil -lBias -llr
 */

/*
   Ax=b の A,近似逆行列R,右辺b から精度保証つきの解の検証 (中尾・山本 p.14)。
   教科書中の点行列、ベクトルの B,Z,x を区間で考えた (IB,IZ,xt)。
*/

#include <iostream>
#include <iomanip>           /* setprecision */
#include "Constants.h"       /* Machine::Epsilon */
#include "IntervalMatrix.h"  /* 区間行列 */
using namespace std;

//#define PRINT
#define ZPRINT
#define LOOP 300
#define DELTA Machine::Epsilon

int main(void)
{
  int i, j, N, k;

  INTERVAL delta(1.0 - DELTA, 1.0 + DELTA);
  cout << setprecision(16) ; 
  cout << "Machine ε= " << Machine::Epsilon << endl;
  cout << "        δ= " << DELTA << endl;

  cout << "何次元ですか ?" << endl;
  cin >> N;

  VECTOR x(N), b(N), Z(N);
  MATRIX A(N,N), R(N,N), I(N,N);
  INTERVAL_VECTOR xt(N), IZ(N), Y(N), X(N);
  INTERVAL_MATRIX IR(N,N), IA(N,N), IB(N,N);

  cout << "行列A" << endl;
  cin >> A;
  cout << A << endl;
  /* IA = Hull(A) */
  IA = A;

  cout << "近似逆行列R" << endl;
  cin >> R;
  cout << R << endl;
  /* IR = Hull(R); */
  IR = R;

  cout << "右辺b" << endl;
  cin >> b;
  cout << b << endl;

  /* 単位行列 */
  Clear(I);
  for (i = 1; i <= N; i++)
    I(i,i) = 1.0;

  IB = I - IR * IA;
  xt = IR * b;
  IZ = IR * (b - IA * xt);

  X = Z;

#ifdef PRINT
  cout << "区間行列 B" << endl << IB << endl;
  cout << "区間ベクトル ~x" << endl << xt << endl;
  cout << "区間行列 Z" << endl << IZ << endl;
  cout << "区間ベクトル X" << endl << X << endl;
  cout << "δ" << endl << delta << endl;
#endif

  for (k = 1; k <= LOOP; k++) {
    /* X のデルタ拡大 */
    Y = delta * X;
    X = IZ + IB * Y;

#ifdef ZPRINT
    cout << k << "回目" << endl ;
    cout << "X="  << X << endl;
    cout << "Y= " << Y << endl;
#endif

    /* (Inf(Y) < Inf(X) && Sup(X) < Sup(Y)) としていた */
    if (X < Y) {
      cout << endl << "反復終了" << endl;
      cout <<  xt  << endl << "+" << endl <<  X  << endl << "=" << endl
           << xt+X << endl
           << "の中に解がある。(反復 " << k << "回)" << endl
           << "Xの幅は" << endl << Diam(X) << endl;
      exit(0);
    }
  }
  cout << endl << "反復回数オーバー" << endl
       <<  xt << endl << "+" << endl <<  X  << endl
       << "=" << endl << xt+X << endl
       << "の中には解はありません。" << endl;

  return 0;
}