main.cpp 8.56 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

Igor Konshin's avatar
Igor Konshin 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
51
52

		//std::cout << "[" << rank << "/" << procs << "]: " << argv[0] << std::endl;
Igor Konshin's avatar
Igor Konshin committed
53
54
55
		Solver::Matrix mat("A"); // Declare the matrix of the linear system to be solved
		Solver::Vector b("rhs"); // Declare the right-hand side vector
		Solver::Vector x("sol"); // Declare the solution vector
56
57
		//std::cout << "[" << rank << "]: load matrix from " << std::string(argv[2]) << " ..." << std::endl;
		double t = Timer(), tt;
Kirill Terekhov's avatar
Kirill Terekhov committed
58
59
		mat.Load(std::string(argv[2])); //if interval parameters not set, matrix will be divided automatically
		BARRIER
60
		if( !rank ) std::cout << "load matrix:    " << Timer() - t << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
61
62
		//mat.Save("test.mtx");
		t = Timer();
63
64
		if( argc > 4 && std::string(argv[3]) == std::string("-")
		             && std::string(argv[4]) == std::string("1") )
Kirill Terekhov's avatar
Kirill Terekhov committed
65
		{
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
			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
83
			b.Load(std::string(argv[3])); // Load RHS vector
Kirill Terekhov's avatar
Kirill Terekhov committed
84
		}
Igor Konshin's avatar
Igor Konshin committed
85
		else // Set local RHS to 1 if it was not specified
Kirill Terekhov's avatar
Kirill Terekhov committed
86
		{
87
			if( !rank ) std::cout << "[" << rank << "]: set RHS=1" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
88
89
90
			INMOST_DATA_ENUM_TYPE mbeg,mend,k;
			mat.GetInterval(mbeg,mend);
			b.SetInterval(mbeg,mend);
Igor Konshin's avatar
Igor Konshin committed
91
			for( k = mbeg; k < mend; ++k ) b[k] = 1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
92
93
		}
		BARRIER
94
95
		if( !rank ) std::cout << "load vector:    " << Timer() - t << std::endl;
		tt = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
96
97
98
		bool success = false;
		int iters;
		double resid, realresid = 0;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
99
		std::string reason;
Kirill Terekhov's avatar
Kirill Terekhov committed
100
		{
Igor Konshin's avatar
Igor Konshin committed
101
			Solver s(type); // Declare the linear solver by specified type
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
102

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

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

			s.SetParameterEnum("condition_estimation",1);
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
116
117
			

Kirill Terekhov's avatar
Kirill Terekhov committed
118
			t = Timer();
Igor Konshin's avatar
Igor Konshin committed
119
			s.SetMatrix(mat); // Compute the preconditioner for the original matrix
Kirill Terekhov's avatar
Kirill Terekhov committed
120
121
122
			BARRIER
			if( !rank ) std::cout << "preconditioner: " << Timer() - t << std::endl;
			t = Timer();
Igor Konshin's avatar
Igor Konshin committed
123
			success = s.Solve(b,x); // Solve the linear system with the previously computted preconditioner
Kirill Terekhov's avatar
Kirill Terekhov committed
124
			BARRIER
125
126
127
128
129
			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
130
131
		}
		tt = Timer() - tt;
Igor Konshin's avatar
Igor Konshin committed
132
133

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

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

Kirill Terekhov's avatar
Kirill Terekhov committed
144
145
146
			{
				INMOST_DATA_ENUM_TYPE mbeg,mend,k;
				info.GetLocalRegion(info.GetRank(),mbeg,mend);
Igor Konshin's avatar
Igor Konshin committed
147
				for( k = mbeg; k < mend; ++k )
Kirill Terekhov's avatar
Kirill Terekhov committed
148
149
150
151
152
153
154
155
				{
					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);
156
			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
157
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
158
			realresid = sqrt(recv[0]/(recv[1]+1.0e-100));
Kirill Terekhov's avatar
Kirill Terekhov committed
159
			//realresid = sqrt(realresid);
Igor Konshin's avatar
Igor Konshin committed
160

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

Kirill Terekhov's avatar
Kirill Terekhov committed
165
		{
166
			if( !rank )
Kirill Terekhov's avatar
Kirill Terekhov committed
167
			{
168
				std::cout << procs << " processors for Solver::type=" << type;
Kirill Terekhov's avatar
Kirill Terekhov committed
169
170
				if( success )
				{
Igor Konshin's avatar
Igor Konshin committed
171
172
					std::cout << " solved in " << tt << " secs";
					std::cout << " with " << iters << " iterations to " << resid << " norm";
Kirill Terekhov's avatar
Kirill Terekhov committed
173
174
				}
				else std::cout << " failed to solve";
Igor Konshin's avatar
Igor Konshin committed
175
176
				std::cout  << " matrix \"" << argv[2] << "\"";
				if( argc > 3 ) std::cout  << " vector \"" << argv[3] << "\"";
177
				std::cout << " true residual ||Ax-b||/||b||=" << realresid;
Kirill Terekhov's avatar
Kirill Terekhov committed
178
				std::cout << std::endl;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
179
				std::cout << "reason: " << reason << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
180
181
			}
		}
182
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

		{ // 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
236
	}
Igor Konshin's avatar
Igor Konshin committed
237
238

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