main.cpp 5.39 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
3
#include <string>
#include <iostream>

Igor Konshin's avatar
Igor Konshin committed
4
#include "../../inmost.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
5
using namespace INMOST;
Igor Konshin's avatar
Igor Konshin committed
6

Kirill Terekhov's avatar
Kirill Terekhov committed
7
8
9
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)
{
	int rank,procs;
	if( argc < 3 )
	{
18
		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] [solver_options.txt]" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
19
20
21
22
23
		return -1;
	}
	Solver::Type type;
	switch(atoi(argv[1]))
	{
24
25
26
27
28
29
30
31
32
33
34
35
		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
36
	}
Igor Konshin's avatar
Igor Konshin committed
37
	Solver::Initialize(&argc,&argv,argc > 4 ? argv[4] : NULL); // Initialize the linear solver in accordance with args
Kirill Terekhov's avatar
Kirill Terekhov committed
38
39
	{
#if defined(USE_MPI)
Igor Konshin's avatar
Igor Konshin committed
40
41
		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
42
43
44
45
46
#else
		rank = 0;
		procs = 1;
#endif
		//std::cout << rank << "/" << procs << " " << argv[0] << std::endl;
Igor Konshin's avatar
Igor Konshin committed
47
48
49
50
		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
		//std::cout << rank << " load matrix from " << std::string(argv[2]) << " ..." << std::endl;
51
		double t = Timer(), tt = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
52
53
54
55
56
		mat.Load(std::string(argv[2])); //if interval parameters not set, matrix will be divided automatically
		BARRIER
		if( !rank ) std::cout << "load matrix: " << Timer() - t << std::endl;
		//mat.Save("test.mtx");
		t = Timer();
Igor Konshin's avatar
Igor Konshin committed
57
		if( argc > 3 )
Kirill Terekhov's avatar
Kirill Terekhov committed
58
59
		{
			//std::cout << rank << " load vector from " << std::string(argv[3]) << std::endl;
Igor Konshin's avatar
Igor Konshin committed
60
			b.Load(std::string(argv[3])); // Load RHS vector
Kirill Terekhov's avatar
Kirill Terekhov committed
61
		}
Igor Konshin's avatar
Igor Konshin committed
62
		else // Set local RHS to 1 if it was not specified
Kirill Terekhov's avatar
Kirill Terekhov committed
63
64
65
66
		{
			INMOST_DATA_ENUM_TYPE mbeg,mend,k;
			mat.GetInterval(mbeg,mend);
			b.SetInterval(mbeg,mend);
Igor Konshin's avatar
Igor Konshin committed
67
			for( k = mbeg; k < mend; ++k ) b[k] = 1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
68
69
70
71
72
73
		}
		BARRIER
		if( !rank ) std::cout << "load vector: " << Timer() - t << std::endl;
		bool success = false;
		int iters;
		double resid, realresid = 0;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
74
		std::string reason;
Kirill Terekhov's avatar
Kirill Terekhov committed
75
		{
Igor Konshin's avatar
Igor Konshin committed
76
			Solver s(type); // Declare the linear solver by specified type
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
77

Kirill Terekhov's avatar
Kirill Terekhov committed
78
79
80
81
			s.SetParameterEnum("gmres_substeps",4);
			s.SetParameterReal("relative_tolerance",1.0e-9);
			s.SetParameterReal("absolute_tolerance",1.0e-16);

Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
82
83
84
85
			s.SetParameterEnum("reorder_nonzeros",0);
			s.SetParameterEnum("rescale_iterations",8);
			s.SetParameterEnum("adapt_ddpq_tolerance",0);
			
Kirill Terekhov's avatar
Kirill Terekhov committed
86
			s.SetParameterReal("drop_tolerance",1.0e-2);
Kirill Terekhov's avatar
Kirill Terekhov committed
87
88
			s.SetParameterReal("reuse_tolerance",1.0e-4);
			s.SetParameterReal("ddpq_tolerance",0.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
89
90

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

Kirill Terekhov's avatar
Kirill Terekhov committed
93
			t = Timer();
Igor Konshin's avatar
Igor Konshin committed
94
			s.SetMatrix(mat); // Compute the preconditioner for the original matrix
Kirill Terekhov's avatar
Kirill Terekhov committed
95
96
97
			BARRIER
			if( !rank ) std::cout << "preconditioner: " << Timer() - t << std::endl;
			t = Timer();
Igor Konshin's avatar
Igor Konshin committed
98
			success = s.Solve(b,x); // Solve the linear system with the previously computted preconditioner
Kirill Terekhov's avatar
Kirill Terekhov committed
99
			BARRIER
Kirill Terekhov's avatar
Kirill Terekhov committed
100
			if( !rank ) std::cout << "solver: " << Timer() - t << "\t\t\t" << std::endl;
Igor Konshin's avatar
Igor Konshin committed
101
102
			iters = s.Iterations(); // Get the number of iterations performed
			resid = s.Residual();   // Get the final residual achieved
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
103
			reason = s.GetReason();
Kirill Terekhov's avatar
Kirill Terekhov committed
104
105
106
			//x.Save("output.rhs");
		}
		tt = Timer() - tt;
Igor Konshin's avatar
Igor Konshin committed
107
108

		{ // Compute the true residual
Kirill Terekhov's avatar
Kirill Terekhov committed
109
110
111
112
113
114
115
			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
116
117
118

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

Kirill Terekhov's avatar
Kirill Terekhov committed
119
120
121
			{
				INMOST_DATA_ENUM_TYPE mbeg,mend,k;
				info.GetLocalRegion(info.GetRank(),mbeg,mend);
Igor Konshin's avatar
Igor Konshin committed
122
				for( k = mbeg; k < mend; ++k )
Kirill Terekhov's avatar
Kirill Terekhov committed
123
124
125
126
127
128
129
130
				{
					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);
Kirill Terekhov's avatar
Kirill Terekhov committed
131
			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
132
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
133
			realresid = sqrt(recv[0]/(recv[1]+1.0e-100));
Kirill Terekhov's avatar
Kirill Terekhov committed
134
			//realresid = sqrt(realresid);
Igor Konshin's avatar
Igor Konshin committed
135

Kirill Terekhov's avatar
Kirill Terekhov committed
136
137
138
			info.RestoreVector(x);
			if( !rank ) std::cout << "norms: " << Timer() - t << std::endl;
		}
Igor Konshin's avatar
Igor Konshin committed
139

Kirill Terekhov's avatar
Kirill Terekhov committed
140
		{
Igor Konshin's avatar
Igor Konshin committed
141
			if( rank == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
142
			{
Igor Konshin's avatar
Igor Konshin committed
143
				std::cout << procs << " processors";
Kirill Terekhov's avatar
Kirill Terekhov committed
144
145
				if( success )
				{
Igor Konshin's avatar
Igor Konshin committed
146
147
					std::cout << " solved in " << tt << " secs";
					std::cout << " with " << iters << " iterations to " << resid << " norm";
Kirill Terekhov's avatar
Kirill Terekhov committed
148
149
				}
				else std::cout << " failed to solve";
Igor Konshin's avatar
Igor Konshin committed
150
151
152
				std::cout  << " matrix \"" << argv[2] << "\"";
				if( argc > 3 ) std::cout  << " vector \"" << argv[3] << "\"";
				std::cout << " true residual ||Ax-b||/||b|| " << realresid;
Kirill Terekhov's avatar
Kirill Terekhov committed
153
				std::cout << std::endl;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
154
				std::cout << "reason: " << reason << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
155
156
157
			}
		}
	}
Igor Konshin's avatar
Igor Konshin committed
158
159

	Solver::Finalize(); // Finalize solver and close MPI activity
Kirill Terekhov's avatar
Kirill Terekhov committed
160
161
	return 0;
}
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
162