main.cpp 5.14 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 )
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
18
		std::cout << "Usage: " << argv[0] << " method_number<0:INNER_ILU2,1:INNER_MLILUC,2:PETSc,3:Trilinos_Aztec,4:Trilinos_Belos,5:Trilinos_Ifpack,6:Trilinos_ML,7:ANI> matrix.mtx [right_hand_side.rhs] [solver_options.txt]" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
19
20
21
22
23
24
		return -1;
	}
	Solver::Type type;
	switch(atoi(argv[1]))
	{
		case 0: type = Solver::INNER_ILU2; break;
Kirill Terekhov's avatar
Kirill Terekhov committed
25
		case 1: type = Solver::INNER_DDPQILUC; break;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
26
27
28
29
30
31
		case 2: type = Solver::PETSc; break;
		case 3: type = Solver::Trilinos_Aztec; break;
		case 4: type = Solver::Trilinos_Belos; break;
		case 5: type = Solver::Trilinos_Ifpack; break;
		case 6: type = Solver::Trilinos_ML; break;
    case 7: type = Solver::ANI; break;
Kirill Terekhov's avatar
Kirill Terekhov committed
32
		case 8: type = Solver::INNER_MPTILUC; break;
Kirill Terekhov's avatar
Kirill Terekhov committed
33
	}
Igor Konshin's avatar
Igor Konshin committed
34
	Solver::Initialize(&argc,&argv,argc > 4 ? argv[4] : NULL); // Initialize the linear solver in accordance with args
Kirill Terekhov's avatar
Kirill Terekhov committed
35
36
	{
#if defined(USE_MPI)
Igor Konshin's avatar
Igor Konshin committed
37
38
		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
39
40
41
42
43
#else
		rank = 0;
		procs = 1;
#endif
		//std::cout << rank << "/" << procs << " " << argv[0] << std::endl;
Igor Konshin's avatar
Igor Konshin committed
44
45
46
47
		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;
48
		double t = Timer(), tt = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
49
50
51
52
53
		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
54
		if( argc > 3 )
Kirill Terekhov's avatar
Kirill Terekhov committed
55
56
		{
			//std::cout << rank << " load vector from " << std::string(argv[3]) << std::endl;
Igor Konshin's avatar
Igor Konshin committed
57
			b.Load(std::string(argv[3])); // Load RHS vector
Kirill Terekhov's avatar
Kirill Terekhov committed
58
		}
Igor Konshin's avatar
Igor Konshin committed
59
		else // Set local RHS to 1 if it was not specified
Kirill Terekhov's avatar
Kirill Terekhov committed
60
61
62
63
		{
			INMOST_DATA_ENUM_TYPE mbeg,mend,k;
			mat.GetInterval(mbeg,mend);
			b.SetInterval(mbeg,mend);
Igor Konshin's avatar
Igor Konshin committed
64
			for( k = mbeg; k < mend; ++k ) b[k] = 1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
65
66
67
68
69
70
		}
		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
71
		std::string reason;
Kirill Terekhov's avatar
Kirill Terekhov committed
72
		{
Igor Konshin's avatar
Igor Konshin committed
73
			Solver s(type); // Declare the linear solver by specified type
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
74

Kirill Terekhov's avatar
Kirill Terekhov committed
75
76
77
78
			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
79
80
81
82
			s.SetParameterEnum("reorder_nonzeros",0);
			s.SetParameterEnum("rescale_iterations",8);
			s.SetParameterEnum("adapt_ddpq_tolerance",0);
			
Kirill Terekhov's avatar
Kirill Terekhov committed
83
84
85
86
87
			s.SetParameterReal("drop_tolerance",1.0e-4);
			s.SetParameterReal("reuse_tolerance",1.0e-8);
			s.SetParameterReal("ddpq_tolerance",0.8);

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

Kirill Terekhov's avatar
Kirill Terekhov committed
90
			t = Timer();
Igor Konshin's avatar
Igor Konshin committed
91
			s.SetMatrix(mat); // Compute the preconditioner for the original matrix
Kirill Terekhov's avatar
Kirill Terekhov committed
92
93
94
			BARRIER
			if( !rank ) std::cout << "preconditioner: " << Timer() - t << std::endl;
			t = Timer();
Igor Konshin's avatar
Igor Konshin committed
95
			success = s.Solve(b,x); // Solve the linear system with the previously computted preconditioner
Kirill Terekhov's avatar
Kirill Terekhov committed
96
			BARRIER
Kirill Terekhov's avatar
Kirill Terekhov committed
97
			if( !rank ) std::cout << "solver: " << Timer() - t << "\t\t\t" << std::endl;
Igor Konshin's avatar
Igor Konshin committed
98
99
			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
100
			reason = s.GetReason();
Kirill Terekhov's avatar
Kirill Terekhov committed
101
102
103
			//x.Save("output.rhs");
		}
		tt = Timer() - tt;
Igor Konshin's avatar
Igor Konshin committed
104
105

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

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

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

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

Kirill Terekhov's avatar
Kirill Terekhov committed
137
		{
Igor Konshin's avatar
Igor Konshin committed
138
			if( rank == 0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
139
			{
Igor Konshin's avatar
Igor Konshin committed
140
				std::cout << procs << " processors";
Kirill Terekhov's avatar
Kirill Terekhov committed
141
142
				if( success )
				{
Igor Konshin's avatar
Igor Konshin committed
143
144
					std::cout << " solved in " << tt << " secs";
					std::cout << " with " << iters << " iterations to " << resid << " norm";
Kirill Terekhov's avatar
Kirill Terekhov committed
145
146
				}
				else std::cout << " failed to solve";
Igor Konshin's avatar
Igor Konshin committed
147
148
149
				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
150
				std::cout << std::endl;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
151
				std::cout << "reason: " << reason << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
152
153
154
			}
		}
	}
Igor Konshin's avatar
Igor Konshin committed
155
156

	Solver::Finalize(); // Finalize solver and close MPI activity
Kirill Terekhov's avatar
Kirill Terekhov committed
157
158
	return 0;
}
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
159