← Back to branch summary

~chaffra/+junk/trilinos

« back to all changes in this revision

Viewing changes to packages/belos/epetra/test/TFQMR/test_tfqmr_hb.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme, Johannes Ring
  • Date: 2009-12-13 12:53:22 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20091213125322-in0nrdjc55deqsw9
Tags: 10.0.3.dfsg-1
[Christophe Prud'homme]
* New upstream release

[Johannes Ring]
* debian/patches/libname.patch: Add prefix 'libtrilinos_' to all
  libraries. 
* debian/patches/soname.patch: Add soversion to libraries.
* debian/watch: Update download URL.
* debian/control:
  - Remove python-numeric from Build-Depends (virtual package).
  - Remove automake and autotools from Build-Depends and add cmake to
    reflect switch to CMake.
  - Add python-support to Build-Depends.
* debian/rules: 
  - Cleanup and updates for switch to CMake.

Show diffs side-by-side

added added

removed removed

Lines of Context:
packages/belos/epetra/test/TFQMR/test_tfqmr_hb.cpp
 
1
// @HEADER
 
2
// ***********************************************************************
 
3
//
 
4
//                 Belos: Block Linear Solvers Package
 
5
//                 Copyright (2004) Sandia Corporation
 
6
//
 
7
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
 
8
// license for use of this work by or on behalf of the U.S. Government.
 
9
//
 
10
// This library is free software; you can redistribute it and/or modify
 
11
// it under the terms of the GNU Lesser General Public License as
 
12
// published by the Free Software Foundation; either version 2.1 of the
 
13
// License, or (at your option) any later version.
 
14
//
 
15
// This library is distributed in the hope that it will be useful, but
 
16
// WITHOUT ANY WARRANTY; without even the implied warranty of
 
17
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
18
// Lesser General Public License for more details.
 
19
//
 
20
// You should have received a copy of the GNU Lesser General Public
 
21
// License along with this library; if not, write to the Free Software
 
22
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 
23
// USA
 
24
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
 
25
//
 
26
// ***********************************************************************
 
27
// @HEADER
 
28
//
 
29
// This driver reads a problem from a Harwell-Boeing (HB) file.
 
30
// The right-hand-side from the HB file is used instead of random vectors.
 
31
// The initial guesses are all set to zero. 
 
32
//
 
33
// NOTE: No preconditioner is used in this case. 
 
34
//
 
35
#include "BelosConfigDefs.hpp"
 
36
#include "BelosLinearProblem.hpp"
 
37
#include "BelosEpetraAdapter.hpp"
 
38
#include "BelosTFQMRSolMgr.hpp"
 
39
#include "createEpetraProblem.hpp"
 
40
#include "Trilinos_Util.h"
 
41
#include "Epetra_CrsMatrix.h"
 
42
#include "Epetra_Map.h"
 
43
#include "Teuchos_CommandLineProcessor.hpp"
 
44
#include "Teuchos_ParameterList.hpp"
 
45
 
 
46
int main(int argc, char *argv[]) {
 
47
  //
 
48
#ifdef EPETRA_MPI       
 
49
  // Initialize MPI     
 
50
  MPI_Init(&argc,&argv);        
 
51
  Belos::MPIFinalize mpiFinalize; // Will call finalize with *any* return
 
52
  (void)mpiFinalize;
 
53
#endif
 
54
  //
 
55
  using Teuchos::RCP;
 
56
  using Teuchos::rcp;
 
57
  //
 
58
  // Get test parameters from command-line processor
 
59
  //  
 
60
  bool verbose = false, proc_verbose = false;
 
61
  int frequency = -1;  // how often residuals are printed by solver
 
62
  int numrhs = 1;  // total number of right-hand sides to solve for
 
63
  int blocksize = 1;  // blocksize used by solver
 
64
  int maxiters = -1;  // maximum number of iterations for solver to use
 
65
  std::string filename("orsirr1.hb");
 
66
  double tol = 1.0e-5;  // relative residual tolerance
 
67
 
 
68
  Teuchos::CommandLineProcessor cmdp(false,true);
 
69
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
 
70
  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
 
71
  cmdp.setOption("tol",&tol,"Relative residual tolerance used by TFQMR solver.");
 
72
  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
 
73
  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
 
74
  cmdp.setOption("block-size",&blocksize,"Block size to be used by TFQMR solver.");
 
75
  cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 := adapted to problem/block size).");
 
76
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
 
77
    return -1;
 
78
  }
 
79
  if (!verbose)
 
80
    frequency = -1;  // reset frequency if test is not verbose
 
81
  //
 
82
  // Get the problem
 
83
  //
 
84
  int MyPID;
 
85
  RCP<Epetra_CrsMatrix> A;
 
86
  RCP<Epetra_MultiVector> B, X;
 
87
  int return_val =Belos::createEpetraProblem(filename,NULL,&A,&B,&X,&MyPID);
 
88
  if(return_val != 0) return return_val;
 
89
  proc_verbose = ( verbose && (MyPID==0) );
 
90
  //
 
91
  // Solve using Belos
 
92
  //
 
93
  typedef double                          ST;
 
94
  typedef Epetra_Operator                 OP;
 
95
  typedef Epetra_MultiVector              MV;
 
96
  typedef Belos::OperatorTraits<ST,MV,OP> OPT;
 
97
  typedef Belos::MultiVecTraits<ST,MV>    MVT;
 
98
  //
 
99
  // *****Construct initial guess and random right-hand-sides *****
 
100
  //
 
101
  if (numrhs == 1) {
 
102
    MVT::MvInit( *X, 1.0 );    
 
103
  } else {
 
104
    X = rcp( new Epetra_MultiVector( A->Map(), numrhs ) );
 
105
    MVT::MvRandom( *X );
 
106
    B = rcp( new Epetra_MultiVector( A->Map(), numrhs ) );
 
107
  }
 
108
  OPT::Apply( *A, *X, *B );
 
109
  MVT::MvInit( *X, 0.0 );
 
110
  //
 
111
  // ********Other information used by block solver***********
 
112
  // *****************(can be user specified)******************
 
113
  //
 
114
  const int NumGlobalElements = B->GlobalLength();
 
115
  if (maxiters == -1)
 
116
    maxiters = NumGlobalElements/blocksize - 1; // maximum number of iterations to run
 
117
  //
 
118
  ParameterList belosList;
 
119
  belosList.set( "Block Size", blocksize );              // Blocksize to be used by iterative solver
 
120
  belosList.set( "Maximum Iterations", maxiters );       // Maximum number of iterations allowed
 
121
  belosList.set( "Convergence Tolerance", tol );         // Relative convergence tolerance requested
 
122
  if (verbose) {
 
123
    belosList.set( "Verbosity", Belos::Errors + Belos::Warnings + 
 
124
                   Belos::TimingDetails + Belos::FinalSummary + Belos::StatusTestDetails );
 
125
    if (frequency > 0)
 
126
      belosList.set( "Output Frequency", frequency );
 
127
  }
 
128
  else
 
129
    belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
 
130
  //
 
131
  // Construct an unpreconditioned linear problem instance.
 
132
  //
 
133
  Belos::LinearProblem<double,MV,OP> problem( A, X, B );
 
134
  bool set = problem.setProblem();
 
135
  if (set == false) {
 
136
    if (proc_verbose)
 
137
      std::cout << std::endl << "ERROR:  Belos::LinearProblem failed to set up correctly!" << std::endl;
 
138
    return -1;
 
139
  }
 
140
  //
 
141
  // Create an iterative solver manager.
 
142
  //
 
143
  RCP< Belos::SolverManager<double,MV,OP> > solver
 
144
    = rcp( new Belos::TFQMRSolMgr<double,MV,OP>(rcp(&problem,false), rcp(&belosList,false)) );
 
145
  //
 
146
  // **********Print out information about problem*******************
 
147
  //
 
148
  if (proc_verbose) {
 
149
    std::cout << std::endl << std::endl;
 
150
    std::cout << "Dimension of matrix: " << NumGlobalElements << std::endl;
 
151
    std::cout << "Number of right-hand sides: " << numrhs << std::endl;
 
152
    std::cout << "Block size used by solver: " << blocksize << std::endl;
 
153
    std::cout << "Max number of TFQMR iterations: " << maxiters << std::endl; 
 
154
    std::cout << "Relative residual tolerance: " << tol << std::endl;
 
155
    std::cout << std::endl;
 
156
  }
 
157
  //
 
158
  // Perform solve
 
159
  //
 
160
  Belos::ReturnType ret = solver->solve();
 
161
  //
 
162
  // Compute actual residuals.
 
163
  //
 
164
  bool badRes = false;
 
165
  std::vector<double> actual_resids( numrhs );
 
166
  std::vector<double> rhs_norm( numrhs );
 
167
  Epetra_MultiVector resid(A->Map(), numrhs);
 
168
  OPT::Apply( *A, *X, resid );
 
169
  MVT::MvAddMv( -1.0, resid, 1.0, *B, resid ); 
 
170
  MVT::MvNorm( resid, actual_resids );
 
171
  MVT::MvNorm( *B, rhs_norm );
 
172
  if (proc_verbose) {
 
173
    std::cout<< "---------- Actual Residuals (normalized) ----------"<<std::endl<<std::endl;
 
174
    for ( int i=0; i<numrhs; i++) {
 
175
      double actRes = actual_resids[i]/rhs_norm[i];
 
176
      std::cout<<"Problem "<<i<<" : \t"<< actRes <<std::endl;
 
177
      if (actRes > tol) badRes = true;
 
178
    }
 
179
  }
 
180
 
 
181
  if (ret!=Belos::Converged || badRes) {
 
182
    if (proc_verbose)
 
183
      std::cout << std::endl << "End Result: TEST FAILED" << std::endl; 
 
184
    return -1;
 
185
  }
 
186
  //
 
187
  // Default return value
 
188
  //
 
189
  if (proc_verbose)
 
190
    std::cout << std::endl << "End Result: TEST PASSED" << std::endl;
 
191
  return 0;
 
192
  //
 
193
} // end test_bl_cg_hb.cpp
 
194
 
 
195
 
 
196
 
 
197