more bicgstab information

              Matrix #  	 Matrix Name  	 nrows  	 ncols  	 nnz 
 is_symmetric 	 iter 	 True Residual 	 Converged 	 error-code
0 	bp__1600.mtx 	822 	822 	4841 	N 	0 	1.#QNAN0 	N  	0
2 	gemat12.mtx 	4929 	4929 	33111 	N 	0 	1.#QNAN0 	N  	0
5 	tub1000.mtx 	1000 	1000 	3996 	N 	9 	0.000005 	N  	0
6 	utm300.mtx 	300 	300 	3155 	N 	15 	0.000012 	N  	0
8 	bcsstm20.mtx 	485 	485 	970 	Y 	8 	0.000021 	N  	0
12 	hor__131.mtx 	434 	434 	4710 	N 	21 	0.00002 	N  	0
13 	psmigr_1.mtx 	3140 	3140 	543162 	N 	0 	1.#QNAN0 	N  	0
18 	bp___800.mtx 	822 	822 	4534 	N 	0 	1.#QNAN0 	N  	0
1 	cavity01.mtx 	317 	317 	7327 	N 	60 	131.827524 	Y  	1
3 	s1rmq4m1.mtx 	5489 	5489 	286600 	Y 	60 	4117.281278 	Y  	1
4 	shl__400.mtx 	663 	663 	1712 	N 	60 	7876.350876 	Y  	1
7 	west2021.mtx 	2021 	2021 	7353 	N 	60 	3344.341296 	Y  	1
9 	e20r5000.mtx 	4241 	4241 	131556 	N 	60 	5618.114418 	Y  	1
10 	fidap037.mtx 	3565 	3565 	67591 	N 	60 	16.887688 	Y  	1
11 	fidapm37.mtx 	9152 	9152 	765944 	N 	60 	95844.2395 	Y  	1
14 	s1rmq4m1.mtx 	5489 	5489 	286600 	Y 	60 	4117.281278 	Y  	1
15 	shl____0.mtx 	663 	663 	1687 	N 	60 	1948.317709 	Y  	1
16 	west2021.mtx 	2021 	2021 	7353 	N 	60 	3344.341296 	Y  	1
17 	af23560.mtx 	23560 	23560 	484256 	N 	60 	52209.26516 	Y  	1
19 	dw8192.mtx 	8192 	8192 	41746 	N 	60 	1563.168486 	Y  	1
20 	memplus.mtx 	17758 	17758 	126150 	N 	60 	11.250001 	Y  	1
21 	orsreg_1.mtx 	2205 	2205 	14133 	N 	60 	0.068578 	Y  	1


Here is a more complete list of matrices from the Matrix-Market 
collection for which I obtain strange results using bicgstab.  As you 
can see bicgstab reports with an error-code of 0 (non-convergence 
according to the ITL documentation).

As you can also see 4 of the candidate non-converging solutions produce 
QNaNs.  The convergence behaviour does not appear to be influenced by 
how the b vector is initialised (lower or higher than 1.0 for each 
element in b).

Regards,

- David


[email protected] wrote:

>Send mtl-devel mailing list submissions to
>	[email protected]
>
>To subscribe or unsubscribe via the World Wide Web, visit
>	http://www.osl.iu.edu/mailman/listinfo.cgi/mtl-devel
>or, via email, send a message with subject or body 'help' to
>	[email protected]
>
>You can reach the person managing the list at
>	[email protected]
>
>When replying, please edit your Subject line so it is more specific
>than "Re: Contents of mtl-devel digest..."
>  
>
>------------------------------------------------------------------------
>
>Today's Topics:
>
>   1. strange bicg vs. bicgstab behaviour for same matrix
>      (david moloney)
>  
>
>
> ------------------------------------------------------------------------
>
> Subject:
> MTL: strange bicg vs. bicgstab behaviour for same matrix
> From:
> david moloney <[email protected]>
> Date:
> Tue, 09 Nov 2004 15:02:49 +0000
> To:
> [email protected], David Moloney <[email protected]>
>
> To:
> [email protected], David Moloney <[email protected]>
>
>
> Hi all,
>
> I am using a number of Matrix-Market matrices as test cases for some 
> work I am doing using mtl and itl.  I have found that a number of 
> these matrices produce QNaNs from the two_norm() function in mtl when 
> bicgstab is used and behave correctly when the bicg solver from itl is 
> used.
> The only way that a QNaN can be produced from mtl::two_norm() 
> according to the IEEE standard, is if the square-root of a negative 
> number is taken.  This however is impossible as two_norm() is the sum 
> of the squares of the elements of a vector b and all of the elements 
> of b are set to 1 initially.  On the first iteration the result 
> returned by two_norm should be sqrt(822), where 822 is the numbers of 
> columns and rows in the A matrix.  Incidentally sqrt(822) = 28.67054 
> as can be seen in the first iteration of both versions of the program 
> below.
>
> The summation on the other hand can produce a QNaN if one or other 
> element of b when squared is infinity, but the fact that bicgstab 
> succeeds suggests this is not the case.
>
> The matrix used is bp__1600.mtx from the matrix-market test-suite 
> (http://math.nist.gov/MatrixMarket/data/Harwell-Boeing/smtape/bp__1600.html) 
> and the code (modified from itl example code) is attached.
> Any opinions welcome,
>
> - David
>
>
> ************************************************************************************** 
>
> bicg results
> ************************************************************************************** 
>
>
> iteration 0: resid 28.6705
> finished! error code = 2
> 0 iterations
> 28.6705 is actual final residual.
> 1 is actual relative tolerance achieved.
> Relative tol: 1e-006  Absolute tol: 0
>
> Solving bp__1600.mtx using itl/bicg
>
> Summary of results for bp__1600.mtx
> -----------------------------------------
> converged                   Y
> # iterations                0
> normb                       28.6705
> final residual              28.6705
> True Residual               28.6705
> Target relative tolerance   1e-006
> Target absolute tolerance   0
> Achieved relative tolerance 1
> -----------------------------------------
>
> ************************************************************************************** 
>
> bicgstab results
> ************************************************************************************** 
>
>
> normr_ = 28.6705
> resid_[1] = r[28.6705]/normb[28.6705]
> iteration 0: resid 28.6705
> normr_ = 1.#QNAN
> resid_[1.#QNAN] = r[1.#QNAN]/normb[28.6705]
> iteration 0: resid 1.#QNAN
> finished! error code = 0
> 0 iterations
> 1.#QNAN is actual final residual.
> 1.#QNAN is actual relative tolerance achieved.
> Relative tol: 1e-006  Absolute tol: 0
>
> Solving bp__1600.mtx using itl/bicgstab
>
> Summary of results for bp__1600.mtx
> ----------------------------------------
> converged                   N
> # iterations                0
> normb                       28.6705
> final residual              1.#QNAN
> True Residual               1.#QNAN
> Target relative tolerance   1e-006
> Target absolute tolerance   0
> Achieved relative tolerance 1.#QNAN
> -----------------------------------------
>
>
>
>------------------------------------------------------------------------
>
>//==========================================================================
>//       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;
>}
>
>
>  
>
>------------------------------------------------------------------------
>
>_______________________________________________
>This list is archived at http://www.lsc.nd.edu/MailArchives/mtl-devel/
>

_______________________________________________
This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/