strange bicg vs. bicgstab behaviour for same matrix

//==========================================================================
//       CVS Information:
//
//	     $RCSfile: parser.cpp,v $  $Revision: 1.3 $  $State: Exp $ 	
//	     $Author: llee $  $Date: 2001/10/26 14:28:24 $ 	
//	     $Locker:  $ 	
//---------------------------------------------------------------------------
//	                                                                     
//	 DESCRIPTION                                                         
//	                                                                     
//---------------------------------------------------------------------------
//
//       LICENSE AGREEMENT
//	     $COPYRIGHT$	
//---------------------------------------------------------------------------
//	                                                                   
//	 REVISION HISTORY:                                                 
//	                                                                   
//	 $Log: parser.cpp,v $
//	 Revision 1.3  2001/10/26 14:28:24  llee
//	 *** empty log message ***
//	
//	 Revision 1.2  2001/07/05 22:28:58  llee1
//	 gcc 3.0 fix
//	
//	 Revision 1.1  2000/07/26 21:49:58  llee1
//	 change file extension from .cc to .cpp
//	
//	 Revision 1.3  2000/07/18 14:30:44  llee1
//	 *** empty log message ***
//	
//	 Revision 1.2  2000/07/17 15:44:06  llee1
//	 *** empty log message ***
//		
//	                                                                    
//===========================================================================
#include <iostream>
#include "parser.h"

#if defined ITL_NO_SSTREAM
#include <stdio.h>
#else
#include <sstream>
#endif

#include <assert.h>

namespace std {

  parser::parser(int argc, char* argv[], const string& delim) {
    string raw_data;
    //copy argv to raw_data;
    for (int i=1; i<argc; ++i) {
      raw_data += argv[i];
      raw_data += " ";
    }

    remove_special_character(raw_data, '=');
    split(delim, raw_data);
  }
  
  void parser::register_flag (const string& flag, int i, const string& mesg) {
    registers[flag] = i;
    messages[flag] = mesg;
  }

  bool parser::flags_coexist (const string& flag1, const string& flag2, 
			      bool coexist) {
    bool found1 = seg.find(flag1) != seg.end();
    bool found2 = seg.find(flag2) != seg.end();
    
    if ( coexist ) {
      if ( found1 != found2 ) {
	cerr << flag1 << " and " << flag2 
	     << " must be both in the comand line. Exiting ..." << endl;
	assert(false);
      }
    } else if ( found1 && found2 ) {
      cerr << flag1 << " and " << flag2 
	   << " cannot be both in the comand line. Exiting ..." << endl;
      //exit
      assert(false);
    }

    return found1 && found2 != coexist;
  }
  
#if defined ITL_NO_SSTREAM
  void parser::get(const string& flag, double& value, double d) { 
    if ( registers.find(flag) == registers.end() ) {
      cerr << flag << " is not registered! Exiting..." << endl;
      assert(false);
    }

    map<string,string>::iterator i = seg.find(flag);

    if ( i == seg.end() ) 
      value = d;
    else {
      if ( registers[flag] == 0 )
	value = ! d; //not the default
      else if ((*i).second.size()==0) {
	cerr << flag << " requires a value." << endl 
	     << " Please provid it right after the flag with whitespace " 
	     << endl << " or equal sign between flag and value." << endl;
	assert(false);
      } else {
	sscanf((*i).second.data(), "%lf", &value);
      }	  
    }
  }

  void parser::get(const string& flag, int& value, int d) { 
    if ( registers.find(flag) == registers.end() ) {
      cerr << flag << " is not registered! Exiting..." << endl;
      assert(false);
    }

    map<string,string>::iterator i = seg.find(flag);

    if ( i == seg.end() ) 
      value = d;
    else {
      if ( registers[flag] == 0 )
	value = ! d; //not the default
      else if ((*i).second.size()==0) {
	cerr << flag << " requires a value." << endl 
	     << " Please provid it right after the flag with whitespace " 
	     << endl << " or equal sign between flag and value." << endl;
	assert(false);
      } else {
	sscanf((*i).second.data(), "%d", &value);
      }	  
    }
  }

  //
  // Added by David Moloney to allow string arguments to be extracted
  // from argv[] ... for example filenames
  //
  void parser::get(const string& flag, char*& value, char* d) { 
    if ( registers.find(flag) == registers.end() ) {
      cerr << flag << " is not registered! Exiting..." << endl;
      assert(false);
    }
    map<string,string>::iterator i = seg.find(flag);
    if (i==seg.end()) value = d;
    else {
      if (registers[flag]==0)
	    value = d; //not the default
      else if ((*i).second.size()==0) {
	    cerr << flag << " requires a value." << endl 
	         << " Please provid it right after the flag with whitespace " 
	         << endl << " or equal sign between flag and value." << endl;
	    assert(false);
      } 
      else {
	    sscanf((*i).second.data(), "%s", value);
      }	  
    }
  }

#else
  template <typename T>
  void parser::get(const string& flag, T& value, const T& d) { 
    if ( registers.find(flag) == registers.end() ) {
      cerr << flag << " is not registered! Exiting..." << endl;
      assert(false);
    }

    map<string,string>::iterator i = seg.find(flag);

    if ( i == seg.end() ) 
      value = d;
    else {
      if ( registers[flag] == 0 )
	value = ! d; //not the default
      else if ((*i).second.size()==0) {
	cerr << flag << " requires a value." << endl 
	     << " Please provid it right after the flag with whitespace " << endl
	     << " or equal sign between flag and value." << endl;
	assert(false);
      } else {
	istringstream iss((*i).second);
	iss >> value;
	
	if (iss.fail()) {
	  cerr << flag << " has an invalid value. Exit..." << endl;
	  //exit
	  assert(false);
	}
      }	  
    }
  }
#endif

  void parser::help(char *msg) {
    string h("--help");
    bool found = seg.find(h) != seg.end();
    if ( found ) {
      cout << endl << endl;
      cout << msg << " Usage of flags: " << endl << endl;
      for ( map<string,int>::iterator i = registers.begin();
	  i != registers.end(); ++i) {
	    cout << "  " << (*i).first;
	    if ( (*i).second ) 
	    cout << " <value> ";
	    cout << ": " << messages[(*i).first] << endl;
      }
      cout << endl << endl;
    }
  }
  
  void parser::remove_special_character(string& raw_data, char c) {
    char * temp = const_cast<char*>(raw_data.c_str());
    for (string::size_type i=0; i<raw_data.size(); ++i)
      if ( temp[i] == c )
	temp[i] = ' ';
  }
  
  void parser::split(const string& delim, const string& raw_data) {

    string::size_type pos = 0;
    bool cond = true;
    
    while ( cond ) {
      string::size_type prev = pos+2;
      pos = raw_data.find(delim, prev);
      
      if ( pos == string::npos ) {
	pos = raw_data.size();
	cond = false;
      }

      string::size_type s = pos - prev + 2;
      
      halve( raw_data.substr(prev-2, s) );
    }
  }
  
  //str is one unit include a flag and its (possible) value(s)
  //it is not allowed that str is started " "
  void parser::halve(const string& str) {
    string::size_type pos = str.find(" ");
    if ( pos == string::npos )
      seg[str] = "";
    else
      seg[str.substr(0,pos)] = str.substr(pos, str.size());
  }

#if !defined ITL_NO_SSTREAM
  template 
  void parser::get(const string& flag, int& value, const int& d);
  template 
  void parser::get(const string& flag, double& value, const double& d);
  template
  void parser::get(const string& flag, bool& value, const bool& d);
#endif  
}

#include "mtl/matrix.h"
#include "mtl/mtl.h"
#include "mtl/utils.h"
#include <itl/interface/mtl.h>

#define __DEBUG__ // itl debug on

#define __BICGSTAB__ // Bi-Conjugate Gradient Stablilised Method
#ifdef __BICGSTAB__
  #include "itl/preconditioner/ilu.h"
  #include "itl/krylov/bicgstab.h"
#endif

//#define __BICG__ // Bi-Conjugate Gradient Method
#ifdef __BICG__
  #include "itl/preconditioner/ilu.h"
  #include "itl/krylov/bicg.h"
#endif

//#define __CG__ // Conjugate Gradient Method
#ifdef __CG__
  #include "itl/preconditioner/cholesky.h"
  #include "itl/krylov/cg.h"
#endif

#include "laplacian.h"
#include "parser.h"

using namespace mtl;
using namespace itl;

#include "mtl/matrix_market_stream.h"
#include <iostream>
#include <fstream>
#include <stdio.h>

using std::cin;
using std::ios;
using std::ofstream;
using std::ifstream;
using std::cerr;
using std::cout;
using std::endl;
using std::parser;

typedef double Type;
typedef matrix< Type, rectangle<>, array< compressed<> >, row_major >::type Matrix;

int main (int argc, char* argv[]) {

  int    converged;
  int    mx, my, max_it;
  double ksp_atol, ksp_rtol;
  int    N;
  int    jj;
  char   *mtxfilename;

  if ( argc == 1 ) {
    cout << argv[0] << " --help will get your the usage of flags." << endl;
  }

  parser myparser(argc, argv);

  myparser.register_flag("-file",     1, "name of Matrix Market format input file");
  myparser.register_flag("-mx",       1, "number of points in x axis in the mesh");
  myparser.register_flag("-my",       1, "number of points in y axis in the mesh");
  myparser.register_flag("-max_it",   1, "maximal iteration allowed");  
  myparser.register_flag("-ksp_atol", 1, "Absolute tolerance in KSP iteration");  
  myparser.register_flag("-ksp_rtol", 1, "Relative tolerance in KSP iteration");  

  myparser.help(argv[0]);

  myparser.get("-file", mtxfilename, "bp__1600.mtx");
  myparser.get("-mx",       mx,       16           );
  myparser.get("-my",       my,       mx           );
  myparser.get("-max_it",   max_it,   60           );
  myparser.get("-ksp_atol", ksp_atol               );
  myparser.get("-ksp_rtol", ksp_rtol, 1.0e-6       );

  N = mx * my;

    matrix_market_stream<Type> mms(mtxfilename);
    Matrix A(mms);

    generate_laplacian_2D(A, mx, my);
    dense1D<Type> x(A.nrows(), Type(0));
    dense1D<Type> b(A.ncols());
    for (dense1D<Type>::iterator i=b.begin(); i!=b.end(); i++) 
      *i = 1; // initialize all elements of b vector to 1

    #ifdef __BICG__
      ILU<Matrix> precond(A);
      noisy_iteration<Type> iter(b, max_it, ksp_rtol, ksp_atol);
      converged = bicg(A, x, b, precond(), iter); //gmres algorithm
      cout << endl << " Solving " << mtxfilename << " using itl/bicg" << endl;
    #endif // __BICG__

    #ifdef __CG__
      cholesky<Matrix> precond(A);
      cholesky<Matrix>::Precond p = precond(); // identity_preconditioner p;
      noisy_iteration<double> iter(b, max_it, ksp_rtol, ksp_atol);
      converged = cg(A, x, b, p, iter); //gmres algorithm
      cout << endl << " Solving " << mtxfilename << " using itl/cg" << endl;
    #endif // __CG__

    #ifdef __BICGSTAB__
      ILU<Matrix> precond(A);
      noisy_iteration<double> iter(b, max_it, ksp_rtol, ksp_atol);
      converged = bicgstab(A, x, b, precond(), iter); //gmres algorithm
      cout << endl << " Solving " << mtxfilename << " using itl/bicgstab" << endl;
    #endif // __BICGSTAB__

    //verify the result
    //
    dense1D<Type> b1(A.ncols());
    itl::mult(A, x, b1);
    itl::add(b1, itl::scaled(b, -1.), b1);

    // report results on VDU
    //
    cout << endl;
    cout << " Summary of results for " << mtxfilename << endl;
    cout << "-------------------------------------------------------- "  << endl;
    cout << " converged                   " << (converged ? "Y" : "N")   << endl;
    cout << " # iterations                " << iter.iterations()         << endl;
    cout << " normb                       " << iter.normb()              << endl;
    cout << " final residual              " << iter.resid()              << endl;
    cout << " True Residual               " << itl::two_norm(b1)         << endl;
    cout << " Target relative tolerance   " << ksp_rtol                  << endl;
    cout << " Target absolute tolerance   " << ksp_atol                  << endl;
    cout << " Achieved relative tolerance " << iter.resid()/iter.normb() << endl;
    cout << "-------------------------------------------------------- "  << endl;
    cout << endl;

  return 0;
}