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)
{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
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
			if( !rank ) std::cout << "[" << rank << "]: set RHS = A * [1]" << std::endl;
68
			Sparse::Vector ex("exact");  // Declare the temporary exact solution vector
69 70 71 72
			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;
73
			Sparse::Matrix a = mat; // Declare a temporary copy of the original matrix 
74 75 76 77 78 79 80 81 82 83
			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

		{ // Compare solution with exact one if available
			if( argc > 4 )
			{
189 190
				Sparse::Vector ex("exact");  // Declare the exact solution vector
				Sparse::Vector err("error"); // Declare the solution error vector
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
				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;
}