solver_bcgsl.hpp 46.9 KB
Newer Older
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1

Kirill Terekhov's avatar
Kirill Terekhov committed
2 3 4
#ifndef __SOLVER_BCGS__
#define __SOLVER_BCGS__

Kirill Terekhov's avatar
Kirill Terekhov committed
5 6 7 8
//\todo
// 1. comply solvers with Method prototype, after TODO in solver_prototypes.hpp is done
// 2. Implement tricks from Read/solver/bcgsl/download.pdf with convex update and true residual correction
// 3. Detect numerical accuracy breakdown - when preconditioned residual is too far from true residual (probably 2 will fix).
Kirill Terekhov's avatar
Kirill Terekhov committed
9 10 11

#include "inmost_solver.h"

Kirill Terekhov's avatar
Kirill Terekhov committed
12
#define PERTURBATE_RTILDE
Kirill Terekhov's avatar
Kirill Terekhov committed
13 14
//#define CONVEX_COMBINATION
#define PSEUDOINVERSE  // same trick as in petsc with pseudoinverse
Kirill Terekhov's avatar
Kirill Terekhov committed
15 16
//#define USE_LAPACK_SVD // use lapack's dgesvd routine instead of built-in svdnxn

Kirill Terekhov's avatar
Kirill Terekhov committed
17
//#if !defined(NDEBUG)
Kirill Terekhov's avatar
Kirill Terekhov committed
18
#define REPORT_RESIDUAL
Kirill Terekhov's avatar
Kirill Terekhov committed
19
//#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
20 21
//#define USE_OMP

Kirill Terekhov's avatar
Kirill Terekhov committed
22 23 24

namespace INMOST
{
Kirill Terekhov's avatar
Kirill Terekhov committed
25 26 27 28 29 30 31
	//lapack svd
#if defined(PSEUDOINVERSE)
#if defined(USE_LAPACK_SVD)
	extern "C"
	{
		void dgesvd_(char*,char*,int*,int*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*);
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
32 33 34 35 36
#else
	//SVD adopted from http://www.public.iastate.edu/~dicook/JSS/paper/code/svd.c
	static INMOST_DATA_REAL_TYPE SIGNFUNC(INMOST_DATA_REAL_TYPE a, INMOST_DATA_REAL_TYPE b) { return b >= 0.0 ? fabs(a) : -fabs(a); }
	static INMOST_DATA_REAL_TYPE MAXFUNC(INMOST_DATA_REAL_TYPE x,INMOST_DATA_REAL_TYPE y) { return x > y? x : y; }
	static INMOST_DATA_REAL_TYPE PYTHAG(INMOST_DATA_REAL_TYPE a, INMOST_DATA_REAL_TYPE b)
Kirill Terekhov's avatar
Kirill Terekhov committed
37
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
38 39 40 41 42
		INMOST_DATA_REAL_TYPE at = fabs(a), bt = fabs(b), ct, result;
		if (at > bt)       { ct = bt / at; result = sqrt(at) * sqrt(at + ct * bt); }
		else if (bt > 0.0) { ct = at / bt; result = sqrt(bt) * sqrt(bt + ct * at); }
		else result = 0.0;
		return(result);
Kirill Terekhov's avatar
Kirill Terekhov committed
43
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
44
	static int svdnxn(INMOST_DATA_REAL_TYPE * pa,INMOST_DATA_REAL_TYPE * pu, INMOST_DATA_REAL_TYPE *pw, INMOST_DATA_REAL_TYPE * pv, const int n)
Kirill Terekhov's avatar
Kirill Terekhov committed
45
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
46 47 48 49
		shell<INMOST_DATA_REAL_TYPE> a(pa,n*n);
		shell<INMOST_DATA_REAL_TYPE> u(pu,n*n);
		shell<INMOST_DATA_REAL_TYPE> v(pv,n*n);
		shell<INMOST_DATA_REAL_TYPE> w(pw,n);
Kirill Terekhov's avatar
Kirill Terekhov committed
50
		//std::copy(a.begin(),a.end(),u.begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
51
		//memcpy(u,a,sizeof(INMOST_DATA_REAL_TYPE)*n*n);
Kirill Terekhov's avatar
Kirill Terekhov committed
52 53 54
		int flag, i, its, j, jj, k, l, nm;
		INMOST_DATA_REAL_TYPE c, f, h, s, x, y, z;
		INMOST_DATA_REAL_TYPE anorm = 0.0, g = 0.0, scale = 0.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
55 56
		dynarray<INMOST_DATA_REAL_TYPE,64> rv1;
    rv1.resize(n);
Kirill Terekhov's avatar
Kirill Terekhov committed
57
		// Householder reduction to bidiagonal form
Kirill Terekhov's avatar
Kirill Terekhov committed
58
		for (i = 0; i < n*n; ++i) u[i] = a[i];
Kirill Terekhov's avatar
Kirill Terekhov committed
59
		for (i = 0; i < n; i++) 
Kirill Terekhov's avatar
Kirill Terekhov committed
60
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
61 62 63 64 65
			// left-hand reduction
			l = i + 1;
			rv1[i] = scale * g;
			g = s = scale = 0.0;
			if (i < n) 
Kirill Terekhov's avatar
Kirill Terekhov committed
66
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
				for (k = i; k < n; k++) 
					scale += fabs(u[k*n+i]);
				if (scale) 
				{
					for (k = i; k < n; k++) 
					{
						u[k*n+i] = u[k*n+i]/scale;
						s += (u[k*n+i] * u[k*n+i]);
					}
					f = u[i*n+i];
					g = -SIGNFUNC(sqrt(s), f);
					h = f * g - s;
					u[i*n+i] = (f - g);
					if (i != n - 1) 
					{
						for (j = l; j < n; j++) 
						{
							for (s = 0.0, k = i; k < n; k++) 
								s += (u[k*n+i] * u[k*n+j]);
							f = s / h;
							for (k = i; k < n; k++) 
								u[k*n+j] += (f * u[k*n+i]);
						}
					}
					for (k = i; k < n; k++) 
						u[k*n+i] = (u[k*n+i]*scale);
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
94
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
95
			w[i] = (scale * g);
Kirill Terekhov's avatar
Kirill Terekhov committed
96

Kirill Terekhov's avatar
Kirill Terekhov committed
97 98 99
			// right-hand reduction 
			g = s = scale = 0.0;
			if (i < n && i != n - 1) 
Kirill Terekhov's avatar
Kirill Terekhov committed
100
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
101 102 103
				for (k = l; k < n; k++) 
					scale += fabs(u[i*n+k]);
				if (scale) 
Kirill Terekhov's avatar
Kirill Terekhov committed
104
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
105
					for (k = l; k < n; k++) 
Kirill Terekhov's avatar
Kirill Terekhov committed
106
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
107 108
						u[i*n+k] = (u[i*n+k]/scale);
						s += (u[i*n+k] * u[i*n+k]);
Kirill Terekhov's avatar
Kirill Terekhov committed
109
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
110 111 112 113 114 115 116
					f = u[i*n+l];
					g = -SIGNFUNC(sqrt(s), f);
					h = f * g - s;
					u[i*n+l] = (f - g);
					for (k = l; k < n; k++) 
						rv1[k] = u[i*n+k] / h;
					if (i != n - 1) 
Kirill Terekhov's avatar
Kirill Terekhov committed
117
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
118 119 120 121 122 123 124
						for (j = l; j < n; j++) 
						{
							for (s = 0.0, k = l; k < n; k++) 
								s += (u[j*n+k] * u[i*n+k]);
							for (k = l; k < n; k++) 
								u[j*n+k] += (s * rv1[k]);
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
125
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
126 127
					for (k = l; k < n; k++) 
						u[i*n+k] = (u[i*n+k]*scale);
Kirill Terekhov's avatar
Kirill Terekhov committed
128
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
129 130 131
			}
			anorm = MAXFUNC(anorm, (fabs(w[i]) + fabs(rv1[i])));
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
132

Kirill Terekhov's avatar
Kirill Terekhov committed
133 134 135 136 137 138
		// accumulate the right-hand transformation 
		for (i = n - 1; i >= 0; i--) 
		{
			if (i < n - 1) 
			{
				if (g) 
Kirill Terekhov's avatar
Kirill Terekhov committed
139
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
140 141 142 143
					for (j = l; j < n; j++)
						v[j*n+i] = ((u[i*n+j] / u[i*n+l]) / g);
					// double division to avoid underflow 
					for (j = l; j < n; j++) 
Kirill Terekhov's avatar
Kirill Terekhov committed
144
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
145 146 147 148
						for (s = 0.0, k = l; k < n; k++) 
							s += (u[i*n+k] * v[k*n+j]);
						for (k = l; k < n; k++) 
							v[k*n+j] += (s * v[k*n+i]);
Kirill Terekhov's avatar
Kirill Terekhov committed
149
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
150
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
151 152
				for (j = l; j < n; j++) 
					v[i*n+j] = v[j*n+i] = 0.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
153
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
154 155 156
			v[i*n+i] = 1.0;
			g = rv1[i];
			l = i;
Kirill Terekhov's avatar
Kirill Terekhov committed
157
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
158

Kirill Terekhov's avatar
Kirill Terekhov committed
159 160
		// accumulate the left-hand transformation
		for (i = n - 1; i >= 0; i--) 
Kirill Terekhov's avatar
Kirill Terekhov committed
161
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
162 163 164 165 166 167 168 169 170
			l = i + 1;
			g = w[i];
			if (i < n - 1) 
				for (j = l; j < n; j++) 
					u[i*n+j] = 0.0;
			if (g) 
			{
				g = 1.0 / g;
				if (i != n - 1) 
Kirill Terekhov's avatar
Kirill Terekhov committed
171
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
172 173 174 175 176 177 178 179
					for (j = l; j < n; j++) 
					{
						for (s = 0.0, k = l; k < n; k++) 
							s += (u[k*n+i] * u[k*n+j]);
						f = (s / u[i*n+i]) * g;
						for (k = i; k < n; k++) 
							u[k*n+j] += (f * u[k*n+i]);
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
180
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
181 182
				for (j = i; j < n; j++) 
					u[j*n+i] = (u[j*n+i]*g);
Kirill Terekhov's avatar
Kirill Terekhov committed
183
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
184
			else 
Kirill Terekhov's avatar
Kirill Terekhov committed
185
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
186 187 188 189 190
				for (j = i; j < n; j++) 
					u[j*n+i] = 0.0;
			}
			++u[i*n+i];
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
191

Kirill Terekhov's avatar
Kirill Terekhov committed
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
		// diagonalize the bidiagonal form
		for (k = n - 1; k >= 0; k--) 
		{                             
			// loop over singular values 
			for (its = 0; its < 30; its++) 
			{                         
				// loop over allowed iterations
				flag = 1;
				for (l = k; l >= 0; l--) 
				{                     
					// test for splitting 
					nm = l - 1;
					if (fabs(rv1[l]) + anorm == anorm) 
					{
						flag = 0;
						break;
					}
					if (fabs(w[nm]) + anorm == anorm) 
						break;
				}
				if (flag) 
				{
					c = 0.0;
					s = 1.0;
					for (i = l; i <= k; i++) 
					{
						f = s * rv1[i];
						if (fabs(f) + anorm != anorm) 
						{
							g = w[i];
							h = PYTHAG(f, g);
							w[i] = h; 
							h = 1.0 / h;
							c = g * h;
							s = (- f * h);
Kirill Terekhov's avatar
Kirill Terekhov committed
227
							for (j = nm < 0 ? 1 : 0; j < n; j++) 
Kirill Terekhov's avatar
Kirill Terekhov committed
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
							{
								y = u[j*n+nm];
								z = u[j*n+i];
								u[j*n+nm] = (y * c + z * s);
								u[j*n+i] = (z * c - y * s);
							}
						}
					}
				}
				z = w[k];
				if (l == k) 
				{                  
					// convergence
					if (z < 0.0) 
					{              
						// make singular value nonnegative
						w[k] = (-z);
						for (j = 0; j < n; j++) 
							v[j*n+k] = (-v[j*n+k]);
					}
					break;
				}
				if (its >= 30) 
				{
					fprintf(stderr, "No convergence after 30,000! iterations \n");
					return 1;
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
255

Kirill Terekhov's avatar
Kirill Terekhov committed
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
				// shift from bottom 2 x 2 minor
				x = w[l];
				nm = k - 1;
				y = w[nm];
				g = rv1[nm];
				h = rv1[k];
				f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
				g = PYTHAG(f, 1.0);
				f = ((x - z) * (x + z) + h * ((y / (f + SIGNFUNC(g, f))) - h)) / x;

				// next QR transformation
				c = s = 1.0;
				for (j = l; j <= nm; j++) 
				{
					i = j + 1;
					g = rv1[i];
					y = w[i];
					h = s * g;
					g = c * g;
					z = PYTHAG(f, h);
					rv1[j] = z;
					c = f / z;
					s = h / z;
					f = x * c + g * s;
					g = g * c - x * s;
					h = y * s;
					y = y * c;
					for (jj = 0; jj < n; jj++) 
					{
						x = v[jj*n+j];
						z = v[jj*n+i];
						v[jj*n+j] = (x * c + z * s);
						v[jj*n+i] = (z * c - x * s);
					}
					z = PYTHAG(f, h);
					w[j] = z;
					if (z) 
					{
						z = 1.0 / z;
						c = f * z;
						s = h * z;
					}
					f = (c * g) + (s * y);
					x = (c * y) - (s * g);
					for (jj = 0; jj < n; jj++) 
					{
						y = u[jj*n+j];
						z = u[jj*n+i];
						u[jj*n+j] = (y * c + z * s);
						u[jj*n+i] = (z * c - y * s);
					}
				}
				rv1[l] = 0.0;
				rv1[k] = f;
				w[k] = x;
Kirill Terekhov's avatar
Kirill Terekhov committed
311 312
			}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
313 314

		for(int i = 0; i < n; i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
315
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
316
			INMOST_DATA_REAL_TYPE temp;
Kirill Terekhov's avatar
Kirill Terekhov committed
317
			for(int j = i + 1; j < n; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
318
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
319 320 321
				temp = u[i+n*j];
				u[i+n*j] = u[j+n*i];
				u[j+n*i] = temp;
Kirill Terekhov's avatar
Kirill Terekhov committed
322 323
			}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
324
		for(i = 0; i < n; i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
325
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
326 327 328 329 330
			k = i;
			for(j = i+1; j < n; ++j)
				if( w[k] < w[j] ) k = j;
			INMOST_DATA_REAL_TYPE temp;
			if( w[k] > w[i] )
Kirill Terekhov's avatar
Kirill Terekhov committed
331
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
332 333 334 335 336 337 338 339 340 341 342 343
				temp = w[k];
				w[k] = w[i];
				w[i] = temp;
				for(int j = 0; j < n; ++j)
				{
					temp = u[k*n+j];
					u[k*n+j] = u[i*n+j];
					u[i*n+j] = temp;
					temp = v[j*n+k];
					v[j*n+k] = v[j*n+i];
					v[j*n+i] = temp;
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
344
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
345
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
346
		return 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
347
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
348
#endif //USE_LAPACK_SVD
Kirill Terekhov's avatar
Kirill Terekhov committed
349
#else //PSEUDOINVERSE
Kirill Terekhov's avatar
Kirill Terekhov committed
350
	static int solvenxn(INMOST_DATA_REAL_TYPE * A, INMOST_DATA_REAL_TYPE * x, INMOST_DATA_REAL_TYPE * b, int n, int * order)
Kirill Terekhov's avatar
Kirill Terekhov committed
351 352 353 354 355
	{
		INMOST_DATA_REAL_TYPE temp, max;
		int temp2;
		for(int i = 0; i < n; i++) order[i] = i;
		for(int i = 0; i < n; i++)
Kirill Terekhov's avatar
Kirill Terekhov committed
356
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421
			int maxk = i, maxq = i;
			max = fabs(A[maxk*n+maxq]);
			//Find best pivot
			for(int q = i; q < n; q++) // over columns
			{
				for(int k = i; k < n; k++) // over rows
				{
					if( fabs(A[k*n+q]) > max )
					{
						max = fabs(A[k*n+q]);
						maxk = k;
						maxq = q;
					}
				}
			}
			//Exchange rows
			if( maxk != i ) 
			{
				for(int q = 0; q < n; q++)
				{
					temp = A[maxk*n+q];
					A[maxk*n+q] = A[i*n+q];
					A[i*n+q] = temp;
				}
				//exchange rhs
				{
					temp = b[maxk];
					b[maxk] = b[i];
					b[i] = temp;
				}
			}
			//Exchange columns
			if( maxq != i ) 
			{
				for(int k = 0; k < n; k++)
				{
					temp = A[k*n+maxq];
					A[k*n+maxq] = A[k*n+i];
					A[k*n+i] = temp;
				}
				//remember order in sol
				{
					temp2 = order[maxq];
					order[maxq] = order[i];
					order[i] = temp2;
				}
			}
			if( fabs(b[i]/A[i*n+i]) > 1.0e+100 )
				return i+1;
		
			for(int k = i+1; k < n; k++)
			{
				A[i*n+k] /= A[i*n+i];
				A[k*n+i] /= A[i*n+i];
			}
			for(int k = i+1; k < n; k++)
			for(int q = i+1; q < n; q++)
			{
				A[k*n+q] -= A[k*n+i] * A[i*n+i] * A[i*n+q];
			}
			for(int j = i+1; j < n; j++) //iterate over columns of L
			{
				b[j] -= b[i] * A[j*n+i];
			}
			b[i] /= A[i*n+i];
Kirill Terekhov's avatar
Kirill Terekhov committed
422
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
423 424 425 426 427 428 429 430

		for(int i = n-1; i >= 0; i--) //iterate over rows of U
			for(int j = i+1; j < n; j++) 
			{
				b[i] -= b[j] * A[i*n+j];
			}
		for(int i = 0; i < n; i++)
			x[order[i]] = b[i];
Kirill Terekhov's avatar
Kirill Terekhov committed
431
	
Kirill Terekhov's avatar
Kirill Terekhov committed
432 433
		return 0;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
434
#endif //PSEUDOINVERSE
Kirill Terekhov's avatar
Kirill Terekhov committed
435 436
	class BCGSL_solver : public IterativeMethod
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
437
    INMOST_DATA_REAL_TYPE length;
Kirill Terekhov's avatar
Kirill Terekhov committed
438 439 440 441
		INMOST_DATA_REAL_TYPE rtol, atol, divtol, last_resid;
		INMOST_DATA_ENUM_TYPE iters, maxits, l, last_it;
		INMOST_DATA_REAL_TYPE resid;
		INMOST_DATA_REAL_TYPE * tau, * sigma, * gamma, *theta1, * theta2, * theta3;
Kirill Terekhov's avatar
Kirill Terekhov committed
442 443
		Sparse::Vector r_tilde, x0, t, * u, * r;
		Sparse::Matrix * Alink;
Kirill Terekhov's avatar
Kirill Terekhov committed
444
		Method * prec;
Kirill Terekhov's avatar
Kirill Terekhov committed
445
		std::string reason;
Kirill Terekhov's avatar
Kirill Terekhov committed
446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469
		Solver::OrderInfo * info;
		bool init;
	public:
		INMOST_DATA_ENUM_TYPE GetIterations() {return last_it;}
		INMOST_DATA_REAL_TYPE GetResidual() {return last_resid;}
		INMOST_DATA_REAL_TYPE & RealParameter(std::string name)
		{
			if (name[0] == ':')
			{
				if (prec != NULL) return prec->RealParameter(name.substr(1, name.size() - 1));
			}
			if (name == "rtol") return rtol;
			else if (name == "atol") return atol;
			else if (name == "divtol") return divtol;
			else if (prec != NULL) return prec->RealParameter(name);
			throw - 1;
		}
		INMOST_DATA_ENUM_TYPE & EnumParameter(std::string name)
		{
			if (name[0] == ':')
			{
				if (prec != NULL) return prec->EnumParameter(name.substr(1, name.size() - 1));
			}
			if (name == "maxits") return maxits;
470 471 472 473 474
			else if (name == "levels") 
			{
				if( init ) throw - 1; //solver was already initialized, value should not be changed
				return l;
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
475 476 477 478 479 480 481 482 483 484 485 486 487
			else if (prec != NULL) return prec->EnumParameter(name);
			throw - 1;
		}
		BCGSL_solver(Method * prec, Solver::OrderInfo & info)
			:rtol(1e-8), atol(1e-9), divtol(1e+40), maxits(1500),l(2),prec(prec),info(&info)
		{
			Alink = NULL;
			init = false;
		}
		bool Initialize()
		{
			if (isInitialized()) Finalize();
			if (prec != NULL && !prec->isInitialized()) prec->Initialize();
Kirill Terekhov's avatar
Kirill Terekhov committed
488 489
			info->PrepareVector(r_tilde);
			info->PrepareVector(x0);
Kirill Terekhov's avatar
Kirill Terekhov committed
490
			info->PrepareVector(t);
Kirill Terekhov's avatar
Kirill Terekhov committed
491 492 493 494
			tau = new INMOST_DATA_REAL_TYPE[l * 3 + (l+1)*(l+1) + (l+1)*2];
			sigma = tau + (l+1)*(l+1);
			gamma = sigma + l+1;
			theta1 = gamma + l+1;
Kirill Terekhov's avatar
Kirill Terekhov committed
495 496
			theta2 = theta1 + l;
			theta3 = theta2 + l;
Kirill Terekhov's avatar
Kirill Terekhov committed
497
			u = new Sparse::Vector[l * 2 + 2];
Kirill Terekhov's avatar
Kirill Terekhov committed
498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542
			r = u + l + 1;
			for (INMOST_DATA_ENUM_TYPE k = 0; k < l + 1; k++)
			{
				info->PrepareVector(r[k]);
				info->PrepareVector(u[k]);
			}
			init = true;
			return true;
		}
		bool isInitialized() { return init && (prec == NULL || prec->isInitialized()); }
		bool Finalize()
		{
			if (isInitialized())
			{
				if (!prec->isFinalized()) prec->Finalize();
				delete[] u;
				delete[] tau;
				init = false;
			}
			return true;
		}
		bool isFinalized() { return !init && (prec == NULL || prec->isFinalized()); }
		void Copy(const Method * other)
		{
			const BCGSL_solver * b = dynamic_cast<const BCGSL_solver *>(other);
			assert(b != NULL);
			rtol = b->rtol;
			atol = b->atol;
			divtol = b->divtol;
			last_resid = b->last_resid;
			iters = b->iters;
			maxits = b->maxits;
			l = b->l;
			last_it = b->last_it;
			resid = b->resid;
			Alink = b->Alink;
			info = b->info;
			if (init) Finalize();
			if (b->prec != NULL)
			{
				if (prec == NULL) prec = b->prec->Duplicate();
				else prec->Copy(b->prec);
			}
			if (b->init) Initialize();
		}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
543
		BCGSL_solver(const BCGSL_solver & other) :IterativeMethod(other)
Kirill Terekhov's avatar
Kirill Terekhov committed
544 545 546 547 548 549 550 551 552 553 554 555 556
		{
			Copy(&other);
		}
		BCGSL_solver & operator =(BCGSL_solver const & other)
		{
			Copy(&other);
			return *this;
		}
		~BCGSL_solver()
		{
			if (!isFinalized()) Finalize();
			if (prec != NULL) delete prec;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
557
		void ApplyOperator(Sparse::Vector & Input, Sparse::Vector & Output)
Kirill Terekhov's avatar
Kirill Terekhov committed
558 559 560 561 562 563 564 565 566 567
		{
			if (prec != NULL) //right preconditioning here! for left preconditioner have to reverse order
			{
				prec->Solve(Input, t); 
				info->Update(t);
				Alink->MatVec(1.0,t,0,Output);
				info->Update(Output);
			}
			else
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
568
				Alink->MatVec(1.0,Input,0,Output);
Kirill Terekhov's avatar
Kirill Terekhov committed
569 570 571
				info->Update(Output);
			}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
572
		bool Solve(Sparse::Vector & RHS, Sparse::Vector & SOL)
Kirill Terekhov's avatar
Kirill Terekhov committed
573 574 575
		{
			assert(isInitialized());
			INMOST_DATA_ENUM_TYPE vbeg,vend, vlocbeg, vlocend;
Kirill Terekhov's avatar
Kirill Terekhov committed
576
			INMOST_DATA_INTEGER_TYPE ivbeg, ivend, ivlocbeg, ivlocend;
Kirill Terekhov's avatar
Kirill Terekhov committed
577
			INMOST_DATA_REAL_TYPE rho0 = 1, rho1 = 1, alpha = 0, beta = 0, omega = 1, eta = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
578
			INMOST_DATA_REAL_TYPE resid0, resid, rhs_norm, tau_sum = 0, sigma_sum = 0,r_tilde_norm = 0;//, temp[2];
Kirill Terekhov's avatar
Kirill Terekhov committed
579 580 581
			iters = 0;
			info->PrepareVector(SOL);
			info->PrepareVector(RHS);
Kirill Terekhov's avatar
Kirill Terekhov committed
582 583
			info->Update(SOL);
			info->Update(RHS);
Kirill Terekhov's avatar
Kirill Terekhov committed
584 585 586 587
			if( prec != NULL ) prec->ReplaceSOL(SOL);
			if( prec != NULL ) prec->ReplaceRHS(RHS);
			info->GetLocalRegion(info->GetRank(),vlocbeg,vlocend);
			info->GetVectorRegion(vbeg,vend);
Kirill Terekhov's avatar
Kirill Terekhov committed
588 589 590 591 592
			ivbeg = vbeg;
			ivend = vend;
			ivlocbeg = vlocbeg;
			ivlocend = vlocend;
			//info->ScalarProd(RHS,RHS,vlocbeg,vlocend,rhs_norm);
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
593
			rhs_norm = 1;
Kirill Terekhov's avatar
Kirill Terekhov committed
594
			//r[0] = b
Kirill Terekhov's avatar
Kirill Terekhov committed
595
			std::copy(RHS.Begin(),RHS.End(),r[0].Begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
596
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
597
				// r[0] = r[0] - A x
Kirill Terekhov's avatar
Kirill Terekhov committed
598
				Alink->MatVec(-1,SOL,1,r[0]); //global multiplication, r probably needs an update
Kirill Terekhov's avatar
Kirill Terekhov committed
599
				info->Update(r[0]); // r is good
600
				std::copy(SOL.Begin(),SOL.End(),x0.Begin()); //x0 = x
Kirill Terekhov's avatar
Kirill Terekhov committed
601
				std::fill(SOL.Begin(),SOL.End(),0.0); //x = 0
Kirill Terekhov's avatar
Kirill Terekhov committed
602
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
603 604
			std::copy(r[0].Begin(),r[0].End(),r_tilde.Begin()); // r_tilde = r[0]
			std::fill(u[0].Begin(),u[0].End(),0); // u[0] = 0
Kirill Terekhov's avatar
Kirill Terekhov committed
605 606 607 608 609 610 611 612
			resid = 0;
#if defined(USE_OMP)
#pragma omp parallel for reduction(+:resid)
#endif
			for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k)
				resid += r[0][k]*r[0][k];
			info->Integrate(&resid,1);
			//info->ScalarProd(r[0],r[0],vlocbeg,vlocend,resid); //resid = dot(r[0],r[0])
Kirill Terekhov's avatar
Kirill Terekhov committed
613
      last_resid = resid = resid0 = sqrt(resid/rhs_norm); //resid = sqrt(dot(r[0],r[0])
Kirill Terekhov's avatar
Kirill Terekhov committed
614 615 616 617
#if defined(USE_OMP)
#pragma omp parallel for
#endif
			for(INMOST_DATA_INTEGER_TYPE k = ivbeg; k < ivend; ++k) // r_tilde = r[0] / dot(r[0],r[0])
Kirill Terekhov's avatar
Kirill Terekhov committed
618
				r_tilde[k] /= resid;
Kirill Terekhov's avatar
Kirill Terekhov committed
619 620 621 622 623 624 625
			last_it = 0;
#if defined(REPORT_RESIDUAL)
			if( info->GetRank() == 0 ) 
			{
				//std::cout << "iter " << last_it << " residual " << resid << std::endl;
				//std::cout << "iter " << last_it << " resid " << resid << "\r";
				//printf("iter %3d resid %12g | %12g relative %12g | %12g\r", last_it, resid, atol, resid / resid0, rtol);
Kirill Terekhov's avatar
Kirill Terekhov committed
626
				printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
Kirill Terekhov's avatar
Kirill Terekhov committed
627 628 629
				fflush(stdout);
			}
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
630
			
Kirill Terekhov's avatar
Kirill Terekhov committed
631
			bool halt = false;
Kirill Terekhov's avatar
Kirill Terekhov committed
632 633 634
			if( last_resid < atol || last_resid < rtol*resid0 ) 
			{
				reason = "initial solution satisfy tolerances";
Kirill Terekhov's avatar
Kirill Terekhov committed
635
				halt = true;
Kirill Terekhov's avatar
Kirill Terekhov committed
636
			}
637

Kirill Terekhov's avatar
Kirill Terekhov committed
638 639 640
#if defined(USE_OMP)
#pragma omp parallel
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
641
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
642
				INMOST_DATA_ENUM_TYPE i = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
643
				while( !halt )
Kirill Terekhov's avatar
Kirill Terekhov committed
644
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
645
					
Kirill Terekhov's avatar
Kirill Terekhov committed
646 647 648
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
649
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
650
						rho0 = -omega*rho0;
Kirill Terekhov's avatar
Kirill Terekhov committed
651
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
652
				
Kirill Terekhov's avatar
Kirill Terekhov committed
653
					
Kirill Terekhov's avatar
Kirill Terekhov committed
654
					for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
655
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
656 657 658
            if( false )
            {
perturbate:
Kirill Terekhov's avatar
Kirill Terekhov committed
659
#if defined(PERTURBATE_RTILDE)
Kirill Terekhov's avatar
Kirill Terekhov committed
660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684
              //std::cout << "Rescue solver by perturbing orthogonal direction!" << std::endl;
              //Compute length of the vector
#if defined(USE_OMP)
#pragma omp single
#endif
              {
                length = static_cast<INMOST_DATA_REAL_TYPE>(ivlocend-ivlocbeg);
                r_tilde_norm = 0.0;
              }
              info->Integrate(&length,1);
              //perform perturbation (note that rand() with openmp may give identical sequences of random values
#if defined(USE_OMP)
#pragma omp for
#endif
              for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
              {
                INMOST_DATA_REAL_TYPE unit = 2*static_cast<INMOST_DATA_REAL_TYPE>(rand())/static_cast<INMOST_DATA_REAL_TYPE>(RAND_MAX)-1.0;
                r_tilde[k] += unit/length;
              }
              //compute norm for orthogonal vector
#if defined(USE_OMP)
#pragma omp for reduction(+:r_tilde_norm)
#endif
              for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
                r_tilde_norm += r_tilde[k]*r_tilde[k];
Kirill Terekhov's avatar
Kirill Terekhov committed
685 686 687 688 689 690
#if defined(USE_OMP)
#pragma omp single
#endif
              {
                r_tilde_norm = sqrt(r_tilde_norm);
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
691 692 693 694 695 696 697 698
              info->Integrate(&r_tilde_norm,1);
              //normalize orthogonal vector to unity
#if defined(USE_OMP)
#pragma omp for
#endif
              for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
                r_tilde[k] /= r_tilde_norm;
              //Recompute rho1
Kirill Terekhov's avatar
Kirill Terekhov committed
699 700 701 702 703
#else
              reason = "r_tilde perturbation algorithm is disabled";
              halt = true;
              break;
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
704
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
705 706 707 708 709 710 711 712 713 714
#if defined(USE_OMP)
#pragma omp single
#endif
						rho1 = 0;
#if defined(USE_OMP)
#pragma omp for reduction(+:rho1)
#endif
						for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
							rho1+= r[j][k]*r_tilde[k];
						info->Integrate(&rho1,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
715 716 717 718 719
            if( fabs(rho1) < 1.0e-50 )
            {
              //std::cout << "Asking to perturbate r_tilde since rho1 is too small " << rho1 << std::endl;
              goto perturbate;
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
720 721 722 723 724 725 726
						//info->ScalarProd(r[j],r_tilde,vlocbeg,vlocend,rho1); // rho1 = dot(r[j],r_tilde)
#if defined(USE_OMP)
#pragma omp single
#endif
						{
							beta = alpha * (rho1/rho0);
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
727

Kirill Terekhov's avatar
Kirill Terekhov committed
728 729 730 731 732 733
						if( fabs(beta) > 1.0e+100 ) 
						{
							//std::cout << "alpha " << alpha << " rho1 " << rho1 << " rho0 " << rho0 << " beta " << beta << std::endl;
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
734 735 736 737
              {
							  reason = "multiplier(1) is too large";
							  halt = true;
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
738 739
							break;
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
740

Kirill Terekhov's avatar
Kirill Terekhov committed
741 742
						if( beta != beta )
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
743
              //std::cout << "alpha " << alpha << " rho1 " << rho1 << " rho0 " << rho0 << " beta " << beta << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
744 745 746
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
747 748 749 750
              {
							  reason = "multiplier(1) is NaN";
							  halt = true;
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
751 752 753 754 755 756 757 758
							break;
						}
#if defined(USE_OMP)
#pragma omp single
#endif
						{
							rho0 = rho1;
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
759
						for(INMOST_DATA_ENUM_TYPE m = 0; m < j+1; m++)
Kirill Terekhov's avatar
Kirill Terekhov committed
760 761 762 763 764
						{
#if defined(USE_OMP)
#pragma omp for
#endif
							for(INMOST_DATA_INTEGER_TYPE k = ivbeg; k < ivend; ++k)
Kirill Terekhov's avatar
Kirill Terekhov committed
765
								u[m][k] = r[m][k] - beta*u[m][k];
Kirill Terekhov's avatar
Kirill Terekhov committed
766
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
767

Kirill Terekhov's avatar
Kirill Terekhov committed
768

Kirill Terekhov's avatar
Kirill Terekhov committed
769 770 771 772 773 774 775 776 777 778 779 780
						ApplyOperator(u[j],u[j+1]); // u[j+1] = A*R*u[j]
#if defined(USE_OMP)
#pragma omp single
#endif
						eta = 0;
#if defined(USE_OMP)
#pragma omp for reduction(+:eta)
#endif
						for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
							eta += u[j+1][k]*r_tilde[k];
						info->Integrate(&eta,1);
						//info->ScalarProd(u[j+1],r_tilde,vlocbeg,vlocend,eta); //eta = dot(u[j+1],r_tilde)
Kirill Terekhov's avatar
Kirill Terekhov committed
781

Kirill Terekhov's avatar
Kirill Terekhov committed
782 783 784 785 786 787
            if( fabs(eta) < 1.0e-50 )
            {
              //std::cout << "Asking to perturbate r_tilde since eta is too small " << eta << std::endl;
              goto perturbate;
            }

Kirill Terekhov's avatar
Kirill Terekhov committed
788 789 790 791 792 793 794 795 796 797
#if defined(USE_OMP)
#pragma omp single
#endif
						alpha = rho0 / eta;

						if( fabs(alpha) > 1.0e+100 ) 
						{
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
798 799 800 801
              {
							  reason = "multiplier(2) is too large";
							  halt = true;
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
802 803 804 805 806 807 808
							break;
						}
						if( alpha != alpha )
						{
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
809 810 811 812
              {
							  reason = "multiplier(2) is NaN";
							  halt = true;
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
813 814 815 816 817 818 819 820 821
							break;
						}

#if defined(USE_OMP)
#pragma omp for
#endif
						for(INMOST_DATA_INTEGER_TYPE k = ivbeg; k < ivend; ++k)
							SOL[k] += alpha*u[0][k];

Kirill Terekhov's avatar
Kirill Terekhov committed
822
						for(INMOST_DATA_ENUM_TYPE m = 0; m < j+1; m++)
Kirill Terekhov's avatar
Kirill Terekhov committed
823 824 825 826 827
						{
#if defined(USE_OMP)
#pragma omp for
#endif
							for(INMOST_DATA_INTEGER_TYPE k = ivbeg; k < ivend; ++k) //r[i] = r[i] - alpha * u[i+1]
Kirill Terekhov's avatar
Kirill Terekhov committed
828
								r[m][k] -= alpha*u[m+1][k];
Kirill Terekhov's avatar
Kirill Terekhov committed
829
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
830 831

					
Kirill Terekhov's avatar
Kirill Terekhov committed
832 833 834 835 836 837 838 839 840 841 842 843 844 845
#if defined(USE_OMP)
#pragma omp single
#endif
						resid = 0;
#if defined(USE_OMP)
#pragma omp for reduction(+:resid)
#endif
						for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
							resid += r[0][k]*r[0][k];
						info->Integrate(&resid,1);
						//info->ScalarProd(r[0],r[0],vlocbeg,vlocend,resid); // resid = dot(r[j],r[j])
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
846 847 848 849
            {
						  resid = sqrt(resid/rhs_norm); // resid = sqrt(dot(r[j],r[j]))
              last_it++;
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
850
					
Kirill Terekhov's avatar
Kirill Terekhov committed
851 852 853 854 855 856 857 858 859 860
						if( resid < atol || resid < rtol*resid0 ) 
						{
#if defined(USE_OMP)
#pragma omp single
#endif
							reason = "early exit in bi-cg block";
							last_resid = resid;
							halt = true;
							break;
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
861
						ApplyOperator(r[j],r[j+1]); // r[j+1] = A*R*r[j]
Kirill Terekhov's avatar
Kirill Terekhov committed
862
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
863

Kirill Terekhov's avatar
Kirill Terekhov committed
864 865 866 867 868 869 870
					if( halt ) break;
					INMOST_DATA_ENUM_TYPE size = l;
#if defined(CONVEX_COMBINATION)
					size = l+1;
#endif
					// Penalization for convex combination for update below
					for(INMOST_DATA_ENUM_TYPE j = 1; j < l+1; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
871
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
872 873 874 875 876 877 878 879 880 881 882
						for(INMOST_DATA_ENUM_TYPE m = 1; m < j+1; m++)
						{
#if defined(USE_OMP)
#pragma omp single
#endif
							tau_sum = 0.0;
#if defined(USE_OMP)
#pragma omp for reduction(+:tau_sum)
#endif
							for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k)
								tau_sum += r[j][k]*r[m][k];
Kirill Terekhov's avatar
Kirill Terekhov committed
883 884 885 886 887 888

              if( fabs(tau_sum) > 1.0e+100 )
              {
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
889 890 891 892
                {
							    reason = "multiplier(tau) is too large";
                  halt = true;
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
893 894 895
              }
              if( tau_sum != tau_sum )
              {
Kirill Terekhov's avatar
Kirill Terekhov committed
896 897 898
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
899 900 901 902
                {
							    reason = "multiplier(tau) is NaN";
                  halt = true;
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
903 904 905 906 907 908 909
              }
#if defined(USE_OMP)
#pragma omp single
#endif
              {
							  tau[(j-1) + (m-1)*size] = tau[(m-1) + (j-1)*size] = tau_sum;
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
910 911 912 913 914 915 916 917 918 919
						}
#if defined(USE_OMP)
#pragma omp single
#endif
						sigma_sum = 0;
#if defined(USE_OMP)
#pragma omp for reduction(+:sigma_sum)
#endif
						for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k)
							sigma_sum += r[0][k]*r[j][k];
Kirill Terekhov's avatar
Kirill Terekhov committed
920 921 922

            if( fabs(sigma_sum) > 1.0e+100 )
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
923 924 925
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
926 927 928 929
              {
							  reason = "multiplier(sigma) is too large";
                halt = true;
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
930 931 932 933 934 935
            }
            if( sigma_sum != sigma_sum )
            {
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
936 937 938 939
              {
							  reason = "multiplier(sigma) is NaN";
                halt = true;
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
940 941 942 943 944 945 946
            }
#if defined(USE_OMP)
#pragma omp single
#endif
            {
						  sigma[j-1] = sigma_sum;
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
947
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
948 949 950

          if( halt ) break;

Kirill Terekhov's avatar
Kirill Terekhov committed
951
#if defined(CONVEX_COMBINATION)
Kirill Terekhov's avatar
Kirill Terekhov committed
952 953 954
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
955
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
956 957 958 959 960 961 962 963 964
						INMOST_DATA_REAL_TYPE lagrangian = 0.0;
						for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++) lagrangian += tau[j+size*j];
						sigma[l] = lagrangian;
						tau[(l+1)*(l+1)-1] = 0.0;
						for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++)
						{
							tau[l + j*(l+1)] = -lagrangian;
							tau[l*(l+1)+j] = lagrangian;
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
965 966
					}
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
967 968 969
					info->Integrate(tau,(l+2)*(l+1)); //sigma is updated with tau
          //info->Integrate(tau,size*size);
          //info->Integrate(sigma,size);
Kirill Terekhov's avatar
Kirill Terekhov committed
970

Kirill Terekhov's avatar
Kirill Terekhov committed
971 972 973
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
974
#if defined(PSEUDOINVERSE)
Kirill Terekhov's avatar
Kirill Terekhov committed
975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033
					{
						int dgesvd_info = 0;
#if defined(USE_LAPACK_SVD)
						char c = 'A';
						INMOST_DATA_REAL_TYPE U[128*128], V[128*128], w[128];
						INMOST_DATA_REAL_TYPE work[5*128];
						int lwork = 5*128;
						int n = static_cast<int>(size);
						dgesvd_(&c,&c,&n,&n,tau,&n,w,U,&n,V,&n,work,&lwork,&dgesvd_info);
#else
						/*
						char c = 'A';
						INMOST_DATA_REAL_TYPE U2[128*128], V2[128*128], w2[128], tau2[128*128];
						INMOST_DATA_REAL_TYPE work[5*128];
						int lwork = 5*128;
						int n = l;
						memcpy(tau2,tau,sizeof(INMOST_DATA_REAL_TYPE)*l*l);
						dgesvd_(&c,&c,&n,&n,tau2,&n,w2,U2,&n,V2,&n,work,&lwork,&dgesvd_info);
						printf("dgesvd\n");
						printf("w\n");
						for(int q = 0; q < l; ++q) printf("%g ",w2[q]);
						printf("\nU\n");
						for(int q = 0; q < l*l; ++q) 
						{
							printf("%g ",U2[q]);
							if( (q+1)%l == 0 ) printf("\n");
						}
						printf("V\n");
						for(int q = 0; q < l*l; ++q) 
						{
							printf("%g ",V2[q]);
							if( (q+1)%l == 0 ) printf("\n");
						}
						*/
						INMOST_DATA_REAL_TYPE U[128*128], V[128*128], w[128];
						dgesvd_info = svdnxn(tau,U,w,V,size);
						//for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++) w[j] = S[j*l+j];
						/*
						printf("svdnxn\n");
						printf("w\n");
						for(int q = 0; q < l; ++q) printf("%g ",w[q]);
						printf("\nU\n");
						for(int q = 0; q < l*l; ++q) 
						{
							printf("%g ",U[q]);
							if( (q+1)%l == 0 ) printf("\n");
						}
						printf("V\n");
						for(int q = 0; q < l*l; ++q) 
						{
							printf("%g ",V[q]);
							if( (q+1)%l == 0 ) printf("\n");
						}
						*/
#endif		
						/*
						printf("w ");
						for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++) printf("%20g ",w[j]);
						printf("\n");
Kirill Terekhov's avatar
Kirill Terekhov committed
1034

Kirill Terekhov's avatar
Kirill Terekhov committed
1035 1036 1037 1038 1039 1040 1041
						printf("U\n");
						for(INMOST_DATA_ENUM_TYPE j = 0; j < l*l; j++) 
						{
							printf("%20g ",U[j]);
							if( (j+1) % l == 0 ) printf("\n");
						}
						printf("\n");
Kirill Terekhov's avatar
Kirill Terekhov committed
1042

Kirill Terekhov's avatar
Kirill Terekhov committed
1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077
						printf("VT\n");
						for(INMOST_DATA_ENUM_TYPE j = 0; j < l*l; j++) 
						{
							printf("%20g ",V[j]);
							if( (j+1) % l == 0 ) printf("\n");
						}
						printf("\n");
						*/
						if( dgesvd_info != 0 )
						{
							printf("(%s:%d) dgesvd %d\n",__FILE__,__LINE__,dgesvd_info);
							exit(-1);
						}
					
						INMOST_DATA_REAL_TYPE maxw = w[0], tol;
						for(INMOST_DATA_ENUM_TYPE j = 1; j < size; j++) if(w[j]>maxw) maxw = w[j];
						tol = size*maxw*1.0e-14;
						memset(gamma,0,sizeof(INMOST_DATA_REAL_TYPE)*size);
						for(INMOST_DATA_ENUM_TYPE j = 0; j < size; j++)
						{
							if( w[j] > tol )
							{
								INMOST_DATA_REAL_TYPE sum = 0;
								for(INMOST_DATA_ENUM_TYPE k = 0; k < size; ++k)
									sum += sigma[k]*U[j*size+k];
								for(INMOST_DATA_ENUM_TYPE k = 0; k < size; ++k)
									gamma[k] += sum/w[j]*V[k*size+j];
							}
						}
					}

					//svdnxn(tau,U,S,V,l);
					//INMOST_DATA_REAL_TYPE inv_tau[64];
					//pseudoinverse(tau,inv_tau,l);
					//matmul(inv_tau,sigma,gamma,l,l,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
1078
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
1079
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098
						int order[128];
						int row = solvenxn(tau,gamma,sigma,size,order);
						/*
						double sum = 0.0;
						for(int j = 0; j < l; ++j) 
						{
							sum += gamma[j];
							std::cout << gamma[j] << " ";
						}
						std::cout << "sum: " << sum;
						//std::cout << " lagrangian: " << gamma[l];
						std::cout << std::endl;
						*/
						if( row != 0 )
						{
							std::cout << "breakdown on row " << row << std::endl;
							reason = "breakdown in matrix inversion in polynomial part";
							break;
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
1099
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1100 1101 1102 1103 1104 1105
#endif
#if defined(USE_OMP)
#pragma omp single
#endif
					omega = gamma[l-1];
					if( fabs(omega) > 1.0e+100 )
Kirill Terekhov's avatar
Kirill Terekhov committed
1106
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1107 1108 1109
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
1110 1111 1112 1113
            {
						  reason = "multiplier(3) is too large";
						  halt = true;
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
1114
						break;
Kirill Terekhov's avatar
Kirill Terekhov committed
1115
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1116
					if( omega != omega )
Kirill Terekhov's avatar
Kirill Terekhov committed
1117
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1118 1119 1120
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
1121 1122 1123 1124
            {
						  reason = "multiplier(3) is NaN";
						  halt = true;
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
1125
						break;
Kirill Terekhov's avatar
Kirill Terekhov committed
1126
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1127
					for(INMOST_DATA_ENUM_TYPE j = 1; j < l+1; ++j)
Kirill Terekhov's avatar
Kirill Terekhov committed
1128
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1129 1130 1131 1132 1133 1134 1135 1136 1137
#if defined(USE_OMP)
#pragma omp for
#endif
						for(INMOST_DATA_INTEGER_TYPE k = ivbeg; k < ivend; ++k)
						{
							u[0][k] -= gamma[j-1]*u[j][k];
							SOL[k]  += gamma[j-1]*r[j-1][k];
							r[0][k] -= gamma[j-1]*r[j][k];
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
1138
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1139 1140
				
				
Kirill Terekhov's avatar
Kirill Terekhov committed
1141
					/*
Kirill Terekhov's avatar
Kirill Terekhov committed
1142
					for(INMOST_DATA_ENUM_TYPE j = 1; j < l+1; j++)
Kirill Terekhov's avatar
Kirill Terekhov committed
1143
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162
						for(INMOST_DATA_ENUM_TYPE i = 1; i < j; i++)
						{
							tau[i-1 + (j-1)*l] = 0;
							for(INMOST_DATA_ENUM_TYPE k = vlocbeg; k < vlocend; ++k)
								tau[i-1 + (j-1)*l] += r[j][k]*r[i][k];
							info->Integrate(&tau[i-1 + (j-1)*l],1);
							tau[i-1 + (j-1)*l] /= sigma[i-1];
							for(INMOST_DATA_ENUM_TYPE k = vbeg; k < vend; ++k)
								r[j][k] -= tau[i-1 + (j-1)*l]*r[i][k];
						}
						INMOST_DATA_REAL_TYPE temp[2] = {0,0};
						for(INMOST_DATA_ENUM_TYPE k = vlocbeg; k < vlocend; ++k)
						{
							temp[0] += r[j][k]*r[j][k];
							temp[1] += r[0][k]*r[j][k];
						}
						info->Integrate(temp,2);
						sigma[j-1] = temp[0];//+1.0e-35; //REVIEW
						theta2[j-1] = temp[1]/sigma[j-1];
Kirill Terekhov's avatar
Kirill Terekhov committed
1163
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179
					omega = theta1[l-1] = theta2[l-1];
					for(INMOST_DATA_ENUM_TYPE j = l-1; j > 0; j--)
					{
						eta = 0;
						for(INMOST_DATA_ENUM_TYPE i = j+1; i < l+1; i++)
							eta += tau[j-1 + (i-1)*l] * theta1[i-1];
						theta1[j-1] = theta2[j-1] - eta;
					}
					for(INMOST_DATA_ENUM_TYPE j = 1; j < l; j++)
					{
						eta = 0;
						for(INMOST_DATA_ENUM_TYPE i = j+1; i < l; i++)
							eta += tau[j-1 + (i-1)*l] * theta1[i];
						theta3[j-1] = theta1[j] + eta;
					}
					for(INMOST_DATA_ENUM_TYPE k = vbeg; k < vend; ++k)
Kirill Terekhov's avatar
Kirill Terekhov committed
1180
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192
						SOL[k] += theta1[0]*r[0][k];
						r[0][k] -= theta2[l-1]*r[l][k];
						u[0][k] -= theta1[l-1]*u[l][k];
					}
					for(INMOST_DATA_ENUM_TYPE j = 1; j < l; j++)
					{
						for(INMOST_DATA_ENUM_TYPE k = vbeg; k < vend; ++k)
						{
							u[0][k] -= theta1[j-1]*u[j][k];
							SOL[k] += theta3[j-1]*r[j][k];
							r[0][k] -= theta2[j-1]*r[j][k];
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
1193 1194
					}
					*/
Kirill Terekhov's avatar
Kirill Terekhov committed
1195
					//last_it = l+1;
Kirill Terekhov's avatar
Kirill Terekhov committed
1196
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211
#if defined(USE_OMP)
#pragma omp single
#endif
						resid = 0;
#if defined(USE_OMP)
#pragma omp for reduction(+:resid)
#endif
						for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k)
							resid += r[0][k]*r[0][k];
						info->Integrate(&resid,1);
						//info->ScalarProd(r[0],r[0],vlocbeg,vlocend,resid);
#if defined(USE_OMP)
#pragma omp single
#endif
						resid = sqrt(resid/rhs_norm);
Kirill Terekhov's avatar
Kirill Terekhov committed
1212
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1213
					
Kirill Terekhov's avatar
Kirill Terekhov committed
1214 1215
#if defined(REPORT_RESIDUAL)
					if( info->GetRank() == 0 ) 
Kirill Terekhov's avatar
Kirill Terekhov committed
1216
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1217 1218 1219 1220 1221 1222
						//std::cout << "iter " << last_it << " residual " << resid << " time " << tt << " matvec " << ts*0.5/l << " precond " << tp*0.5/l << std::endl;
						//std::cout << "iter " << last_it << " resid " << resid << "\r";
						//printf("iter %3d resid %12g | %12g relative %12g | %12g\r", last_it, resid, atol, resid / resid0, rtol);
#if defined(USE_OMP)
#pragma omp single
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
1223
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
1224 1225
							printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
							fflush(stdout);
Kirill Terekhov's avatar
Kirill Terekhov committed
1226 1227
						}
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1228 1229 1230
#endif
#if defined(USE_OMP)
#pragma omp single
Kirill Terekhov's avatar
Kirill Terekhov committed
1231
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
1232 1233
					last_resid = resid;
					if( resid != resid )
Kirill Terekhov's avatar
Kirill Terekhov committed
1234
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1235 1236 1237 1238 1239
#if defined(USE_OMP)
#pragma omp single
#endif
						reason = "residual is NAN";
						break;
Kirill Terekhov's avatar
Kirill Terekhov committed
1240
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1241
					if( resid < atol )
Kirill Terekhov's avatar
Kirill Terekhov committed
1242
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1243 1244 1245 1246 1247
#if defined(USE_OMP)
#pragma omp single
#endif
						reason = "converged due to absolute tolerance";
						break;
Kirill Terekhov's avatar
Kirill Terekhov committed
1248
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1249
					if( resid < rtol*resid0 )
Kirill Terekhov's avatar
Kirill Terekhov committed
1250
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1251 1252 1253 1254 1255
#if defined(USE_OMP)
#pragma omp single
#endif
						reason = "converged due to relative tolerance";
						break;
Kirill Terekhov's avatar
Kirill Terekhov committed
1256
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1257
					if( resid > divtol )
Kirill Terekhov's avatar
Kirill Terekhov committed
1258
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
1259 1260 1261 1262 1263
#if defined(USE_OMP)
#pragma omp single
#endif
						reason = "diverged due to divergence tolerance";
						break;
Kirill Terekhov's avatar
Kirill Terekhov committed
1264
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
1265 1266 1267 1268
					if( i == maxits )
					{
#if defined(USE_OMP)
#pragma omp single
Kirill Terekhov's avatar
Kirill Terekhov committed
1269
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
1270 1271 1272 1273
						reason = "reached maximum iteration number";
						break;
					}
					i++;
Kirill Terekhov's avatar
Kirill Terekhov committed
1274
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
1275 1276 1277
			}
			if (prec != NULL)
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1278 1279
				prec->Solve(SOL, r_tilde); //undo right preconditioner
				std::copy(r_tilde.Begin(), r_tilde.End(), SOL.Begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1280
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
1281 1282 1283 1284
#if defined(USE_OMP)
#pragma omp parallel for
#endif
			for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) //undo shift
Kirill Terekhov's avatar
Kirill Terekhov committed
1285
				SOL[k] += x0[k];
Kirill Terekhov's avatar
Kirill Terekhov committed
1286 1287 1288 1289 1290 1291
			//info->RestoreMatrix(A);
			info->RestoreVector(SOL);
			info->RestoreVector(RHS);
			if( last_resid < atol || last_resid < rtol*resid0 ) return true;
			return false;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
1292 1293 1294
		bool ReplaceMAT(Sparse::Matrix & A) { if (isInitialized()) Finalize(); if (prec != NULL) prec->ReplaceMAT(A);  Alink = &A; return true; }
		bool ReplaceRHS(Sparse::Vector & RHS) { (void) RHS; return true; }
		bool ReplaceSOL(Sparse::Vector & SOL) { (void) SOL; return true; }
Kirill Terekhov's avatar
Kirill Terekhov committed
1295
		Method * Duplicate() { return new BCGSL_solver(*this);}
Kirill Terekhov's avatar
Kirill Terekhov committed
1296
		std::string GetReason() {return reason;}
Kirill Terekhov's avatar
Kirill Terekhov committed
1297 1298 1299 1300 1301 1302 1303 1304
	};


	class BCGS_solver : public IterativeMethod
	{
		INMOST_DATA_REAL_TYPE rtol, atol, divtol, last_resid;
		INMOST_DATA_ENUM_TYPE iters, maxits, last_it;
		INMOST_DATA_REAL_TYPE resid;
Kirill Terekhov's avatar
Kirill Terekhov committed
1305 1306
		Sparse::Vector r0, p, y, s, t, z, r, v;
		Sparse::Matrix * Alink;
Kirill Terekhov's avatar
Kirill Terekhov committed
1307 1308 1309
		Method * prec;
		Solver::OrderInfo * info;
		bool init;
Kirill Terekhov's avatar
Kirill Terekhov committed
1310
		std::string reason;
Kirill Terekhov's avatar
Kirill Terekhov committed
1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385
	public:
		INMOST_DATA_ENUM_TYPE GetIterations() {return last_it;}
		INMOST_DATA_REAL_TYPE GetResidual() {return last_resid;}
		INMOST_DATA_REAL_TYPE & RealParameter(std::string name)
		{
			if (name[0] == ':')
			{
				if (prec != NULL) return prec->RealParameter(name.substr(1, name.size() - 1));
			}
			if (name == "rtol") return rtol;
			else if (name == "atol") return atol;
			else if (name == "divtol") return divtol;
			else if( prec != NULL ) return prec->RealParameter(name);
			throw - 1;
		}
		INMOST_DATA_ENUM_TYPE & EnumParameter(std::string name)
		{
			if (name[0] == ':')
			{
				if (prec != NULL) return prec->EnumParameter(name.substr(1, name.size() - 1));
			}
			if (name == "maxits") return maxits;
			else if (prec != NULL) return prec->EnumParameter(name);
			throw - 1;
		}
		BCGS_solver(Method * prec, Solver::OrderInfo & info)
			:rtol(1e-8), atol(1e-11), divtol(1e+40), iters(0), maxits(1500),prec(prec),info(&info)
		{
			init = false;
		}
		bool Initialize()
		{
			assert(Alink != NULL);
			if (isInitialized()) Finalize();
			if (prec != NULL && !prec->isInitialized()) prec->Initialize();
			info->PrepareVector(r);
			info->PrepareVector(v);
			info->PrepareVector(p);
			info->PrepareVector(y);
			info->PrepareVector(s);
			info->PrepareVector(t);
			info->PrepareVector(z);
			info->PrepareVector(r0);
			init = true;
			return true;
		}
		bool isInitialized() { return init && (prec == NULL || prec->isInitialized()); }
		bool Finalize()
		{
			if (prec != NULL && !prec->isFinalized()) prec->Finalize();
			init = false;
			return true;
		}
		bool isFinalized() { return !init && (prec == NULL || prec->isFinalized()); }
		void Copy(const Method * other)
		{
			const BCGS_solver * b = dynamic_cast<const BCGS_solver *>(other);
			assert(b != NULL);
			info = b->info;
			rtol = b->rtol;
			atol = b->atol;
			divtol = b->divtol;
			maxits = b->maxits;
			last_resid = b->last_resid;
			iters = b->iters;
			last_it = b->last_it;
			resid = b->resid;
			Alink = b->Alink;
			if (b->prec != NULL)
			{
				if (prec == NULL) prec = b->prec->Duplicate();
				else prec->Copy(b->prec);
			}
			if (b->init) Initialize();
		}
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1386
		BCGS_solver(const BCGS_solver & other) : IterativeMethod(other)
Kirill Terekhov's avatar
Kirill Terekhov committed
1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399
		{
			Copy(&other);
		}
		BCGS_solver & operator =(BCGS_solver const & other)
		{
			Copy(&other);
			return *this;
		}
		~BCGS_solver()
		{
			if (!isFinalized()) Finalize();
			if (prec != NULL) delete prec;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
1400
		bool Solve(Sparse::Vector & RHS, Sparse::Vector & SOL)
Kirill Terekhov's avatar
Kirill Terekhov committed
1401 1402
		{
			assert(isInitialized());
Kirill Terekhov's avatar
Kirill Terekhov committed
1403
			INMOST_DATA_REAL_TYPE tempa = 0.0, tempb=0.0, r0_norm = 0, length;
Kirill Terekhov's avatar
Kirill Terekhov committed
1404
			INMOST_DATA_ENUM_TYPE vbeg,vend, vlocbeg, vlocend;
Kirill Terekhov's avatar
Kirill Terekhov committed
1405
			INMOST_DATA_INTEGER_TYPE ivbeg,ivend, ivlocbeg, ivlocend;
Kirill Terekhov's avatar
Kirill Terekhov committed
1406 1407 1408
			INMOST_DATA_REAL_TYPE rho = 1, rho1 = 0, alpha = 1, beta = 0, omega = 1;
      INMOST_DATA_REAL_TYPE resid0, resid, temp[2] = {0,0};
      iters = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1409 1410
			info->PrepareVector(SOL);
			info->PrepareVector(RHS);
Kirill Terekhov's avatar
Kirill Terekhov committed
1411 1412
			info->Update(SOL);
			info->Update(RHS);
Kirill Terekhov's avatar
Kirill Terekhov committed
1413 1414 1415 1416
			if (prec != NULL)prec->ReplaceSOL(SOL);
			if (prec != NULL)prec->ReplaceRHS(RHS);
			info->GetLocalRegion(info->GetRank(),vlocbeg,vlocend);
			info->GetVectorRegion(vbeg,vend);
Kirill Terekhov's avatar
Kirill Terekhov committed
1417 1418 1419 1420
			ivbeg = vbeg;
			ivend = vend;
			ivlocbeg = vlocbeg;
			ivlocend = vlocend;
Kirill Terekhov's avatar
Kirill Terekhov committed
1421
      
Kirill Terekhov's avatar
Kirill Terekhov committed
1422
			std::copy(RHS.Begin(),RHS.End(),r.Begin());
Kirill Terekhov's avatar
Kirill Terekhov committed
1423 1424 1425 1426
			{
				Alink->MatVec(-1,SOL,1,r); //global multiplication, r probably needs an update
				info->Update(r); // r is good
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
1427
			std::copy(r.Begin(),r.End(),r0.Begin());
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1428 1429
			std::fill(v.Begin(),v.End(),0.0);
			std::fill(p.Begin(),p.End(),0.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
1430 1431
			{
				resid = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1432
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
1433
#pragma omp parallel for reduction(+:resid)
Kirill Terekhov's avatar
Kirill Terekhov committed
1434 1435
#endif
				for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; k++) 
Kirill Terekhov's avatar
Kirill Terekhov committed
1436
					resid += r0[k]*r0[k];
Kirill Terekhov's avatar
Kirill Terekhov committed
1437
				info->Integrate(&resid,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
1438 1439
			}
			last_resid = resid = resid0 = sqrt(resid);
Kirill Terekhov's avatar
Kirill Terekhov committed
1440 1441 1442 1443 1444
#if defined(USE_OMP)
#pragma omp parallel for
#endif
			for(INMOST_DATA_INTEGER_TYPE k = ivbeg; k < ivend; ++k) // r_tilde = r[0] / dot(r[0],r[0])
				r0[k] /= resid;
Kirill Terekhov's avatar
Kirill Terekhov committed
1445 1446 1447 1448 1449 1450 1451
			last_it = 0;
#if defined(REPORT_RESIDUAL)
			if( info->GetRank() == 0 ) 
			{
				//std::cout << "iter " << last_it << " residual " << resid << std::endl;
				//std::cout << "iter " << last_it << " resid " << resid << "\r";
				//printf("iter %3d resid %12g | %12g relative %12g | %12g\r",last_it,resid,atol,resid/resid0,rtol);
Kirill Terekhov's avatar
Kirill Terekhov committed
1452
				printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
Kirill Terekhov's avatar
Kirill Terekhov committed
1453 1454 1455
				fflush(stdout);
			}
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
1456 1457 1458 1459 1460 1461 1462 1463

      bool halt = false;
			if( last_resid < atol || last_resid < rtol*resid0 ) 
			{
				reason = "initial solution satisfy tolerances";
				halt = true;
			}

Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1464
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
1465
#pragma omp parallel
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1466
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
1467
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
1468
				INMOST_DATA_ENUM_TYPE i = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
1469
				while( !halt )
Kirill Terekhov's avatar
Kirill Terekhov committed
1470
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
1471
					
1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516
            if( false )
            {
perturbate:
#if defined(PERTURBATE_RTILDE)
              //std::cout << "Rescue solver by perturbing orthogonal direction!" << std::endl;
              //Compute length of the vector
#if defined(USE_OMP)
#pragma omp single
#endif
              {
                length = static_cast<INMOST_DATA_REAL_TYPE>(ivlocend-ivlocbeg);
                r0_norm = 0.0;
              }
              info->Integrate(&length,1);
              //perform perturbation (note that rand() with openmp may give identical sequences of random values
#if defined(USE_OMP)
#pragma omp for
#endif
              for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
              {
                INMOST_DATA_REAL_TYPE unit = 2*static_cast<INMOST_DATA_REAL_TYPE>(rand())/static_cast<INMOST_DATA_REAL_TYPE>(RAND_MAX)-1.0;
                r0[k] += unit/length;
              }
              //compute norm for orthogonal vector