main.cpp 8.75 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
#include <string>
#include <iostream>
3
4
5
#include <iomanip>
#include <cmath>
#include <math.h>
Kirill Terekhov's avatar
Kirill Terekhov committed
6

Kirill Terekhov's avatar
Kirill Terekhov committed
7
#include "inmost.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
8
using namespace INMOST;
Igor Konshin's avatar
Igor Konshin committed
9

Kirill Terekhov's avatar
Kirill Terekhov committed
10
11
12
13
14
15
16
17
#if defined(USE_MPI)
#define BARRIER MPI_Barrier(MPI_COMM_WORLD);
#else
#define BARRIER
#endif

int main(int argc, char ** argv)
{
18
	if( argc < 3 || argc > 1 && ( atoi(argv[1]) < 0 || atoi(argv[1]) > 11 ) )
Kirill Terekhov's avatar
Kirill Terekhov committed
19
	{
20
21
22
		std::cout << "Usage: " << argv[0] << " method_number<0:INNER_ILU2,1:INNER_DDPQILUC,2:INNER_MPTILUC,3:INNER_MPTILU2,4:Trilinos_Aztec,5:Trilinos_Belos,6:Trilinos_ML,7:Trilinos_Ifpack,8:PETSc,9:ANI,10:FCBIILU2,11:K3BIILU2> matrix.mtx [right_hand_side.rhs] [exact_solution] [solver_options.txt]" << std::endl;
		std::cout << "Example: " << argv[0] << "  0 a.mtx b.rhs" << std::endl;
		std::cout << "or just: " << argv[0] << " 11 a.mtx - 1 database.txt" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
23
24
25
		return -1;
	}
	Solver::Type type;
26
	switch(atoi(argv[1])) // Actually: type = atoi(argv[1])
Kirill Terekhov's avatar
Kirill Terekhov committed
27
	{
28
29
30
31
32
33
34
35
36
37
38
39
		case  0: type = Solver::INNER_ILU2;      break;
		case  1: type = Solver::INNER_DDPQILUC;  break;
		case  2: type = Solver::INNER_MPTILUC;   break;
		case  3: type = Solver::INNER_MPTILU2;   break;
		case  4: type = Solver::Trilinos_Aztec;  break;
		case  5: type = Solver::Trilinos_Belos;  break;
		case  6: type = Solver::Trilinos_ML;     break;
		case  7: type = Solver::Trilinos_Ifpack; break;
		case  8: type = Solver::PETSc;           break;
		case  9: type = Solver::ANI;             break;
		case 10: type = Solver::FCBIILU2;        break;
		case 11: type = Solver::K3BIILU2;        break;
Kirill Terekhov's avatar
Kirill Terekhov committed
40
	}
41
	Solver::Initialize(&argc,&argv,argc > 5 ? argv[5] : NULL); // Initialize the linear solver in accordance with args
Kirill Terekhov's avatar
Kirill Terekhov committed
42
	{
43
		int rank, procs;
Kirill Terekhov's avatar
Kirill Terekhov committed
44
#if defined(USE_MPI)
Igor Konshin's avatar
Igor Konshin committed
45
46
		MPI_Comm_rank(MPI_COMM_WORLD,&rank);  // Get the rank of the current process
		MPI_Comm_size(MPI_COMM_WORLD,&procs); // Get the total number of processors used
Kirill Terekhov's avatar
Kirill Terekhov committed
47
48
49
50
#else
		rank = 0;
		procs = 1;
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
51
52
    if( rank == 0 )
      std::cout << "Solving with " << Solver::TypeName(type) << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
53
		//std::cout << rank << "/" << procs << " " << argv[0] << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
54
55
56
		Sparse::Matrix mat("A"); // Declare the matrix of the linear system to be solved
		Sparse::Vector b("rhs"); // Declare the right-hand side vector
		Sparse::Vector x("sol"); // Declare the solution vector
Igor Konshin's avatar
Igor Konshin committed
57
		//std::cout << rank << " load matrix from " << std::string(argv[2]) << " ..." << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
58
		double t = Timer(), tt = Timer(), tsol = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
59
60
		mat.Load(std::string(argv[2])); //if interval parameters not set, matrix will be divided automatically
		BARRIER
61
		if( !rank ) std::cout << "load matrix:    " << Timer() - t << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
62
63
		//mat.Save("test.mtx");
		t = Timer();
64
65
		if( argc > 4 && std::string(argv[3]) == std::string("-")
		             && std::string(argv[4]) == std::string("1") )
Kirill Terekhov's avatar
Kirill Terekhov committed
66
		{
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
			if( !rank ) std::cout << "[" << rank << "]: set RHS = A * [1]" << std::endl;
			Solver::Vector ex("exact");  // Declare the temporary exact solution vector
			INMOST_DATA_ENUM_TYPE mbeg,mend,k;
			mat.GetInterval(mbeg,mend);
			ex.SetInterval(mbeg,mend);
			for( k = mbeg; k < mend; ++k ) ex[k] = 1.0;
			Solver::Matrix a = mat; // Declare a temporary copy of the original matrix 
			Solver::OrderInfo info;
			info.PrepareMatrix(a,0);
			info.PrepareVector(ex);
			info.Update(ex);
			a.MatVec(1.0,ex,0.0,b); // Multiply the original matrix by a vector
			//b.Save("b.rhs");      // Save the computed RHS if required
		}
		else if( argc > 3 && std::string(argv[3]) != std::string("1") )
		{
			if( !rank ) std::cout << "[" << rank << "]: load vector from " << std::string(argv[3]) << std::endl;
Igor Konshin's avatar
Igor Konshin committed
84
			b.Load(std::string(argv[3])); // Load RHS vector
Kirill Terekhov's avatar
Kirill Terekhov committed
85
		}
Igor Konshin's avatar
Igor Konshin committed
86
		else // Set local RHS to 1 if it was not specified
Kirill Terekhov's avatar
Kirill Terekhov committed
87
		{
88
			if( !rank ) std::cout << "[" << rank << "]: set RHS=1" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
89
90
91
			INMOST_DATA_ENUM_TYPE mbeg,mend,k;
			mat.GetInterval(mbeg,mend);
			b.SetInterval(mbeg,mend);
Igor Konshin's avatar
Igor Konshin committed
92
			for( k = mbeg; k < mend; ++k ) b[k] = 1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
93
94
		}
		BARRIER
95
96
		if( !rank ) std::cout << "load vector:    " << Timer() - t << std::endl;
		tt = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
97
98
99
		bool success = false;
		int iters;
		double resid, realresid = 0;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
100
		std::string reason;
Kirill Terekhov's avatar
Kirill Terekhov committed
101
		{
Igor Konshin's avatar
Igor Konshin committed
102
			Solver s(type); // Declare the linear solver by specified type
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
103

Kirill Terekhov's avatar
Kirill Terekhov committed
104
			s.SetParameterEnum("gmres_substeps",4);
105
			s.SetParameterReal("relative_tolerance",1.0e-6);
Kirill Terekhov's avatar
Kirill Terekhov committed
106
107
			s.SetParameterReal("absolute_tolerance",1.0e-16);

Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
108
109
110
111
			s.SetParameterEnum("reorder_nonzeros",0);
			s.SetParameterEnum("rescale_iterations",8);
			s.SetParameterEnum("adapt_ddpq_tolerance",0);
			
112
113
			s.SetParameterReal("drop_tolerance",3.0e-3);
			s.SetParameterReal("reuse_tolerance",1.0e-5);
Kirill Terekhov's avatar
Kirill Terekhov committed
114
			s.SetParameterReal("ddpq_tolerance",0.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
115
116

			s.SetParameterEnum("condition_estimation",1);
Kirill Terekhov's avatar
Kirill Terekhov committed
117
      s.SetParameterEnum("schwartz_overlap",3);
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
118
119
			

Kirill Terekhov's avatar
Kirill Terekhov committed
120
			t = Timer();
Igor Konshin's avatar
Igor Konshin committed
121
			s.SetMatrix(mat); // Compute the preconditioner for the original matrix
Kirill Terekhov's avatar
Kirill Terekhov committed
122
			BARRIER
Kirill Terekhov's avatar
Kirill Terekhov committed
123
      tsol += Timer()-t;
Kirill Terekhov's avatar
Kirill Terekhov committed
124
125
			if( !rank ) std::cout << "preconditioner: " << Timer() - t << std::endl;
			t = Timer();
Igor Konshin's avatar
Igor Konshin committed
126
			success = s.Solve(b,x); // Solve the linear system with the previously computted preconditioner
Kirill Terekhov's avatar
Kirill Terekhov committed
127
			BARRIER
128
129
130
131
132
			if( !rank ) std::cout << "iterations:     " << Timer() - t << std::endl;
			iters  = s.Iterations(); // Get the number of iterations performed
			resid  = s.Residual();   // Get the final residual achieved
			reason  = s.GetReason(); // Get the convergence reason
			//x.Save("output.sol");  // Save the solution if required
Kirill Terekhov's avatar
Kirill Terekhov committed
133
134
		}
		tt = Timer() - tt;
Igor Konshin's avatar
Igor Konshin committed
135
136

		{ // Compute the true residual
Kirill Terekhov's avatar
Kirill Terekhov committed
137
			double aresid = 0, bresid = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
138
			Sparse::Vector test;
Kirill Terekhov's avatar
Kirill Terekhov committed
139
140
141
142
143
			t = Timer();
			Solver::OrderInfo info;
			info.PrepareMatrix(mat,0);
			info.PrepareVector(x);
			info.Update(x);
Igor Konshin's avatar
Igor Konshin committed
144
145
146

			mat.MatVec(1.0,x,0.0,test); // Multiply the original matrix by a vector

Kirill Terekhov's avatar
Kirill Terekhov committed
147
148
149
			{
				INMOST_DATA_ENUM_TYPE mbeg,mend,k;
				info.GetLocalRegion(info.GetRank(),mbeg,mend);
Igor Konshin's avatar
Igor Konshin committed
150
				for( k = mbeg; k < mend; ++k )
Kirill Terekhov's avatar
Kirill Terekhov committed
151
152
153
154
155
156
157
158
				{
					aresid += (test[k]-b[k])*(test[k]-b[k]);
					bresid += b[k]*b[k];
				}
			}
			double temp[2] = {aresid,bresid}, recv[2] = {aresid,bresid};
#if defined(USE_MPI)
			MPI_Reduce(temp,recv,2,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
159
			if( info.GetRank() == 0 ) std::cout << "||Ax-b||=" << sqrt(recv[0]) << " ||b||=" << sqrt(recv[1]) << " ||Ax-b||/||b||=" << sqrt(recv[0]/(recv[1]+1.0e-100)) << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
160
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
161
			realresid = sqrt(recv[0]/(recv[1]+1.0e-100));
Kirill Terekhov's avatar
Kirill Terekhov committed
162
			//realresid = sqrt(realresid);
Igor Konshin's avatar
Igor Konshin committed
163

Kirill Terekhov's avatar
Kirill Terekhov committed
164
165
166
			info.RestoreVector(x);
			if( !rank ) std::cout << "norms: " << Timer() - t << std::endl;
		}
Igor Konshin's avatar
Igor Konshin committed
167

Kirill Terekhov's avatar
Kirill Terekhov committed
168
		{
169
			if( !rank )
Kirill Terekhov's avatar
Kirill Terekhov committed
170
			{
171
				std::cout << procs << " processors for Solver::type=" << type;
Kirill Terekhov's avatar
Kirill Terekhov committed
172
173
				if( success )
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
174
					std::cout << " solved in " << tt << " secs (without reading " << tsol << " secs)";
Igor Konshin's avatar
Igor Konshin committed
175
					std::cout << " with " << iters << " iterations to " << resid << " norm";
Kirill Terekhov's avatar
Kirill Terekhov committed
176
177
				}
				else std::cout << " failed to solve";
Igor Konshin's avatar
Igor Konshin committed
178
179
				std::cout  << " matrix \"" << argv[2] << "\"";
				if( argc > 3 ) std::cout  << " vector \"" << argv[3] << "\"";
180
				std::cout << " true residual ||Ax-b||/||b||=" << realresid;
Kirill Terekhov's avatar
Kirill Terekhov committed
181
				std::cout << std::endl;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
182
				std::cout << "reason: " << reason << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
183
184
			}
		}
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238

		{ // Compare solution with exact one if available
			if( argc > 4 )
			{
				Solver::Vector ex("exact");  // Declare the exact solution vector
				Solver::Vector err("error"); // Declare the solution error vector
				INMOST_DATA_ENUM_TYPE mbeg,mend,k;
				mat.GetInterval(mbeg,mend);
				err.SetInterval(mbeg,mend);
				if( std::string(argv[4]) == std::string("1"))
				{
					ex.SetInterval(mbeg,mend);
					for( k = mbeg; k < mend; ++k ) ex[k] = 1.0;
				}
				else
				{
					//std::cout << "[" << rank << "]: load exact solution vector from " << std::string(argv[4]) << std::endl;
					ex.Load(std::string(argv[4])); // Load exact solution vector
				}
				BARRIER
				double dif1 = 0, dif2 = 0, dif8 = 0, norm = 0;
				//std::cout << "[" << rank << "]: mbeg=" << mbeg << " mend=" << mend << std::endl;
				for( k = mbeg; k < mend; ++k )
				{
					double dloc = err[k] = fabs(x[k] - ex[k]);
					dif1 += dloc;
					dif2 += dloc*dloc;
					dif8 = (dif8 > dloc) ? dif8 : dloc;
					norm += fabs(ex[k]);
				}
#if defined(USE_MPI)
				if( procs > 1 )
				{
					double temp = dif1;
					MPI_Reduce(&temp,&dif1,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
					temp = dif2;
					MPI_Reduce(&temp,&dif2,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
					temp = dif8;
					MPI_Reduce(&temp,&dif8,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
					temp = norm;
					MPI_Reduce(&temp,&norm,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
				}
#endif
				dif2 = sqrt(dif2);
				norm += 1.0e-100;
				if( !rank )
				{
					std::cout << "Difference with exact solution \"" << std::string(argv[4]) << "\": " << std::scientific << std::setprecision(6) << std::endl;
					std::cout << "dif1 = " << dif1 << "  dif2 = " << dif2 << "  dif8 = " << dif8 << "  ||ex||_1 = " << norm << std::endl;
					std::cout << "rel1 = " << dif1/norm << "  rel2 = " << dif2/norm << "  rel8 = " << dif8/norm << std::endl;
				}
				//err.Save("error.sol");
			}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
239
	}
Igor Konshin's avatar
Igor Konshin committed
240
241

	Solver::Finalize(); // Finalize solver and close MPI activity
Kirill Terekhov's avatar
Kirill Terekhov committed
242
243
	return 0;
}