solver_bcgsl.hpp 60.5 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

#include "inmost_solver.h"
11
#include <iomanip>
Kirill Terekhov's avatar
Kirill Terekhov committed
12

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

Kirill Terekhov's avatar
Kirill Terekhov committed
18

Kirill Terekhov's avatar
Kirill Terekhov committed
19

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

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

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

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

        // diagonalize the bidiagonal form
Dmitry Bagaev's avatar
Dmitry Bagaev committed
190
191
        for (k = n - 1; k >= 0; k--)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
192
            // loop over singular values
Dmitry Bagaev's avatar
Dmitry Bagaev committed
193
194
            for (its = 0; its < 30; its++)
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
195
196
                // loop over allowed iterations
                flag = 1;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
197
198
                for (l = k; l >= 0; l--)
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
199
200
                    // test for splitting
                    nm = l - 1;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
201
202
                    if (fabs(rv1[l]) + anorm == anorm)
                    {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
203
204
205
                        flag = 0;
                        break;
                    }
Kirill Terekhov's avatar
Kirill Terekhov committed
206
                    if( l == 0 ) break;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
207
208
209
                    if (fabs(w[nm]) + anorm == anorm)
                        break;
                }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
210
211
                if (flag)
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
212
213
                    c = 0.0;
                    s = 1.0;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
214
215
                    for (i = l; i <= k; i++)
                    {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
216
                        f = s * rv1[i];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
217
218
                        if (fabs(f) + anorm != anorm)
                        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
219
220
221
222
223
                            g = w[i];
                            h = PYTHAG(f, g);
                            w[i] = h;
                            h = 1.0 / h;
                            c = g * h;
224
                            s = (- f * h);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
225
226
                            for (j = nm < 0 ? 1 : 0; j < n; j++)
                            {
227
228
229
230
                                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);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
231
232
233
234
235
                            }
                        }
                    }
                }
                z = w[k];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
236
237
                if (l == k)
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
238
                    // convergence
Dmitry Bagaev's avatar
Dmitry Bagaev committed
239
240
                    if (z < 0.0)
                    {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
241
242
243
                        // make singular value nonnegative
                        w[k] = (-z);
                        for (j = 0; j < n; j++)
244
                            v[j*n+k] = (-v[j*n+k]);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
245
246
247
                    }
                    break;
                }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
248
249
                if (its >= 30)
                {
250
                    std::cerr << "No convergence after 30,000! iterations" << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
                    return 1;
                }

                // 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;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
266
267
                for (j = l; j <= nm; j++)
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
268
269
270
271
272
273
274
275
276
277
278
279
280
                    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;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
281
282
                    for (jj = 0; jj < n; jj++)
                    {
283
284
285
286
                        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);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
287
288
289
                    }
                    z = PYTHAG(f, h);
                    w[j] = z;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
290
291
                    if (z)
                    {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
292
293
294
295
296
297
                        z = 1.0 / z;
                        c = f * z;
                        s = h * z;
                    }
                    f = (c * g) + (s * y);
                    x = (c * y) - (s * g);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
298
299
                    for (jj = 0; jj < n; jj++)
                    {
300
301
302
303
                        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);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
304
305
306
307
308
309
310
311
                    }
                }
                rv1[l] = 0.0;
                rv1[k] = f;
                w[k] = x;
            }
        }

312
        for(int i = 0; i < n; i++)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
313
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
314
            INMOST_DATA_REAL_TYPE temp;
315
            for(int j = i + 1; j < n; j++)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
316
            {
317
318
319
                temp = u[i+n*j];
                u[i+n*j] = u[j+n*i];
                u[j+n*i] = temp;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
320
321
            }
        }
322
        for(i = 0; i < n; i++)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
323
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
324
            k = i;
325
326
            for(j = i+1; j < n; ++j)
                if( w[k] < w[j] ) k = j;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
327
            INMOST_DATA_REAL_TYPE temp;
328
            if( w[k] > w[i] )
Dmitry Bagaev's avatar
Dmitry Bagaev committed
329
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
330
331
332
                temp = w[k];
                w[k] = w[i];
                w[i] = temp;
333
                for(int j = 0; j < n; ++j)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
334
                {
335
336
337
338
339
340
                    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;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
341
342
343
344
345
                }
            }
        }
        return 0;
    }
Kirill Terekhov's avatar
Kirill Terekhov committed
346
#endif //USE_LAPACK_SVD
Kirill Terekhov's avatar
Kirill Terekhov committed
347
#else //PSEUDOINVERSE
Dmitry Bagaev's avatar
Dmitry Bagaev committed
348
    static int solvenxn(INMOST_DATA_REAL_TYPE * A, INMOST_DATA_REAL_TYPE * x, INMOST_DATA_REAL_TYPE * b, int n, int * order)
349
350
351
352
353
354
355
356
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
422
423
424
425
426
427
428
429
430
431
	{
		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++)
		{
			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];
		}

		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];

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

Dmitry Bagaev's avatar
Dmitry Bagaev committed
631
632

            bool halt = false;
633
            if( last_resid < atol || last_resid < rtol*resid0 )
Dmitry Bagaev's avatar
Dmitry Bagaev committed
634
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
635
636
637
                reason = "initial solution satisfy tolerances";
                halt = true;
            }
638

Kirill Terekhov's avatar
Kirill Terekhov committed
639
640
641
#if defined(USE_OMP)
#pragma omp parallel
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
642
643
            {
                INMOST_DATA_ENUM_TYPE i = 0;
644
                while( !halt )
Dmitry Bagaev's avatar
Dmitry Bagaev committed
645
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
646

Kirill Terekhov's avatar
Kirill Terekhov committed
647
648
649
#if defined(USE_OMP)
#pragma omp single
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
650
                    {
651
                        rho0 = -omega*rho0;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
652
653
654
                    }


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

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

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

Kirill Terekhov's avatar
Kirill Terekhov committed
769

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

783
                        if( fabs(eta) < 1.0e-50 )
Dmitry Bagaev's avatar
Dmitry Bagaev committed
784
                        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
785
786
787
                            //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

Kirill Terekhov's avatar
Kirill Terekhov committed
789
790
791
#if defined(USE_OMP)
#pragma omp single
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
792
                        alpha = rho0 / eta;
Kirill Terekhov's avatar
Kirill Terekhov committed
793

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

#if defined(USE_OMP)
#pragma omp for
#endif
820
821
                        for(INMOST_DATA_INTEGER_TYPE k = ivbeg; k < ivend; ++k)
                            SOL[k] += alpha*u[0][k];
Kirill Terekhov's avatar
Kirill Terekhov committed
822

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

Kirill Terekhov's avatar
Kirill Terekhov committed
832

Kirill Terekhov's avatar
Kirill Terekhov committed
833
834
835
#if defined(USE_OMP)
#pragma omp single
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
836
                        resid = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
837
838
839
#if defined(USE_OMP)
#pragma omp for reduction(+:resid)
#endif
840
841
842
                        for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k)
                            resid += r[0][k]*r[0][k];
                        info->Integrate(&resid,1);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
843
                        //info->ScalarProd(r[0],r[0],vlocbeg,vlocend,resid); // resid = dot(r[j],r[j])
Kirill Terekhov's avatar
Kirill Terekhov committed
844
845
846
#if defined(USE_OMP)
#pragma omp single
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
847
                        {
848
                            resid = sqrt(resid/rhs_norm); // resid = sqrt(dot(r[j],r[j]))
Dmitry Bagaev's avatar
Dmitry Bagaev committed
849
850
851
                            last_it++;
                        }

852
                        if( resid < atol || resid < rtol*resid0 )
Dmitry Bagaev's avatar
Dmitry Bagaev committed
853
                        {
Kirill Terekhov's avatar
Kirill Terekhov committed
854
855
856
#if defined(USE_OMP)
#pragma omp single
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
857
858
859
860
861
                            reason = "early exit in bi-cg block";
                            last_resid = resid;
                            halt = true;
                            break;
                        }
862
                        ApplyOperator(r[j],r[j+1]); // r[j+1] = A*R*r[j]
Dmitry Bagaev's avatar
Dmitry Bagaev committed
863
                    }
Kirill Terekhov's avatar
Kirill Terekhov committed
864

865
                    if( halt ) break;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
866
                    INMOST_DATA_ENUM_TYPE size = l;
Kirill Terekhov's avatar
Kirill Terekhov committed
867
#if defined(CONVEX_COMBINATION)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
868
                    size = l+1;
Kirill Terekhov's avatar
Kirill Terekhov committed
869
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
870
                    // Penalization for convex combination for update below
871
                    for(INMOST_DATA_ENUM_TYPE j = 1; j < l+1; j++)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
872
                    {
873
                        for(INMOST_DATA_ENUM_TYPE m = 1; m < j+1; m++)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
874
                        {
Kirill Terekhov's avatar
Kirill Terekhov committed
875
876
877
#if defined(USE_OMP)
#pragma omp single
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
878
                            tau_sum = 0.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
879
880
881
#if defined(USE_OMP)
#pragma omp for reduction(+:tau_sum)
#endif
882
883
                            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
884

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

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

950
                    if( halt ) break;
Kirill Terekhov's avatar
Kirill Terekhov committed
951

Kirill Terekhov's avatar
Kirill Terekhov committed
952
#if defined(CONVEX_COMBINATION)
Kirill Terekhov's avatar
Kirill Terekhov committed
953
954
955
#if defined(USE_OMP)
#pragma omp single
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
956
957
958
959
960
961
962
963
964
965
966
967
                    {
                        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