sparse.cpp 56.9 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
3
4
5
6
7
#define _CRT_SECURE_NO_WARNINGS
#include "inmost_sparse.h"
#include <fstream>
#include <sstream>

namespace INMOST
{
8
9
    namespace Sparse
    {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
10

Kirill Terekhov's avatar
Kirill Terekhov committed
11
#if defined(USE_SOLVER) || defined(USE_AUTODIFF)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
12
13
14
15
        bool _hasRowEntryType = false;

        INMOST_MPI_Type RowEntryType = INMOST_MPI_DATATYPE_NULL;

16
        INMOST_MPI_Type GetRowEntryType() {return RowEntryType;}
Dmitry Bagaev's avatar
Dmitry Bagaev committed
17

18
        bool HaveRowEntryType() {return _hasRowEntryType;}
Dmitry Bagaev's avatar
Dmitry Bagaev committed
19

20
21
        void CreateRowEntryType()
        {
22
#if defined(USE_MPI)
23
24
            if( !HaveRowEntryType() )
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
25
                int ierr;
26
27
                MPI_Datatype type[3] = { INMOST_MPI_DATA_ENUM_TYPE, INMOST_MPI_DATA_REAL_TYPE, MPI_UB};
                int blocklen[3] = { 1, 1, 1 };
Dmitry Bagaev's avatar
Dmitry Bagaev committed
28
                MPI_Aint disp[3];
29
30
                disp[0] = offsetof(Sparse::Row::entry,first);
                disp[1] = offsetof(Sparse::Row::entry,second);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
31
32
                disp[2] = sizeof(Sparse::Row::entry);
                ierr = MPI_Type_create_struct(3, blocklen, disp, type, &RowEntryType);
33
34
                if( ierr != MPI_SUCCESS )
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
35
36
37
                    std::cout << __FILE__ << ":" << __LINE__ << "problem in MPI_Type_create_struct" << std::endl;
                }
                ierr = MPI_Type_commit(&RowEntryType);
38
39
                if( ierr != MPI_SUCCESS )
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
40
41
42
                    std::cout << __FILE__ << ":" << __LINE__ << "problem in MPI_Type_commit" << std::endl;
                }
            }
43
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
44
45
46
            _hasRowEntryType = true;
        }

47
48
        void DestroyRowEntryType()
        {
49
#if defined(USE_MPI)
50
51
            if( HaveRowEntryType() )
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
52
53
54
                MPI_Type_free(&RowEntryType);
                RowEntryType = INMOST_MPI_DATATYPE_NULL;
            }
55
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
56
57
58
            _hasRowEntryType = false;
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
59
////////class RowMerger
60
61
62
63
64
65
66
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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

        RowMerger::RowMerger() : Sorted(true), Nonzeros(0), IntervalBeg(0), IntervalEnd(0) {}

        INMOST_DATA_REAL_TYPE & RowMerger::operator[] (INMOST_DATA_ENUM_TYPE pos)
        {
            pos = MapIndex(pos);
            if( LinkedList[pos+1].first != UNDEF ) return LinkedList[pos+1].second;
            else
            {
                INMOST_DATA_ENUM_TYPE index = LinkedList.get_interval_beg(), next;
                if( Sorted )
                {
                    next = index;
                    while(next < pos+1)
                    {
                        index = next;
                        next = LinkedList[index].first;
                    }
                    assert(index < pos+1);
                    assert(pos+1 < next);
                    ++Nonzeros;
                    LinkedList[index].first = pos+1;
                    LinkedList[pos+1].first = next;
                    return LinkedList[pos+1].second;
                }
                else
                {
                    INMOST_DATA_ENUM_TYPE index = LinkedList.get_interval_beg();
                    ++Nonzeros;
                    LinkedList[pos+1].first = LinkedList[index].first;
                    LinkedList[index].first = pos+1;
                    return LinkedList[pos+1].second;
                }
            }
        }

        INMOST_DATA_REAL_TYPE RowMerger::operator[] (INMOST_DATA_ENUM_TYPE pos) const
        {
            pos = MapIndex(pos);
            if( LinkedList[pos+1].first != UNDEF ) return LinkedList[pos+1].second;
            else throw -1;
        }

        RowMerger::RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted)
                : Sorted(Sorted), Nonzeros(0), IntervalBeg(interval_begin), IntervalEnd(interval_end), LinkedList(interval_begin,interval_end+1,Row::make_entry(UNDEF,0.0))
        {
            LinkedList.begin()->first = EOL;
        }

        RowMerger::RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, const std::vector<INMOST_DATA_ENUM_TYPE> & Pre, const std::vector<INMOST_DATA_ENUM_TYPE> & Post, bool Sorted)
                : Sorted(Sorted), Nonzeros(0), IntervalBeg(interval_begin), IntervalEnd(interval_end), NonlocalPre(Pre), NonlocalPost(Post), LinkedList(interval_begin-Pre.size(),interval_end+Post.size()+1,Row::make_entry(UNDEF,0.0))
        {
            //assert(std::is_sorted(Pre.begin(),Pre.end()));
            //assert(std::is_sorted(Post.begin(),Post.end()));
            LinkedList.begin()->first = EOL;
        }

        void RowMerger::Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool _Sorted)
        {
            LinkedList.set_interval_beg(interval_begin-NonlocalPre.size());
            LinkedList.set_interval_end(interval_end+1+NonlocalPost.size());
            IntervalBeg = interval_begin;
            IntervalEnd = interval_end;
            std::fill(LinkedList.begin(),LinkedList.end(),Row::make_entry(UNDEF,0.0));
            LinkedList.begin()->first = EOL;
            Nonzeros = 0;
            Sorted = _Sorted;
        }

        void RowMerger::Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, const std::vector<INMOST_DATA_ENUM_TYPE> & Pre, const std::vector<INMOST_DATA_ENUM_TYPE> & Post, bool _Sorted)
        {
            NonlocalPre = Pre;
            NonlocalPost = Post;
            Resize(interval_begin,interval_end,_Sorted);
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
135
#if defined(USE_SOLVER)
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
        void RowMerger::Resize(const Matrix & A, bool _Sorted)
        {
            INMOST_DATA_ENUM_TYPE mbeg, mend, k, l, ind;
            //retrive interval of indices
            A.GetInterval(mbeg,mend);
            //gather non-local mapping from matrix
            std::set<INMOST_DATA_ENUM_TYPE> Pre, Post;
            for(k = mbeg; k < mend; ++k)
            {
                for(l = 0; l < A[k].Size(); ++l)
                {
                    ind = A[k].GetIndex(l);
                    if( ind < mbeg ) Pre.insert(ind);
                    else if( ind >= mend ) Post.insert(ind);
                }
            }
            //sort to prepare for binary search
            NonlocalPre.clear();
            NonlocalPre.insert(NonlocalPre.end(),Pre.begin(),Pre.end());
            NonlocalPost.clear();
            NonlocalPre.insert(NonlocalPost.end(),Post.begin(),Post.end());
            Resize(mbeg,mend,_Sorted);
        }

        RowMerger::RowMerger(const Matrix & A, bool Sorted) : Sorted(Sorted), Nonzeros(0)
        {
            Resize(A,Sorted);
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
164
#endif
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
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

        RowMerger::~RowMerger() {}

        void RowMerger::SetNonlocal(const std::vector<INMOST_DATA_ENUM_TYPE> & Pre, const std::vector<INMOST_DATA_ENUM_TYPE> & Post)
        {
            //assert(std::is_sorted(Pre.begin(),Pre.end()));
            //assert(std::is_sorted(Post.begin(),Post.end()));
            NonlocalPre = Pre;
            NonlocalPost = Post;
            Resize(IntervalBeg,IntervalEnd,Sorted);
        }


        INMOST_DATA_ENUM_TYPE RowMerger::MapIndex(INMOST_DATA_ENUM_TYPE pos) const
        {
            if( pos < IntervalBeg )
            {
                assert(!NonlocalPre.empty()); //there are indices provided
                std::vector< INMOST_DATA_ENUM_TYPE >::const_iterator search = std::lower_bound(NonlocalPre.begin(),NonlocalPre.end(),pos);
                assert(*search == pos); //is there such index?
                return IntervalBeg - NonlocalPre.size() + static_cast<INMOST_DATA_ENUM_TYPE>(search-NonlocalPre.begin());
            }
            else if( pos >= IntervalEnd )
            {
                assert(!NonlocalPost.empty()); //there are indices provided
                std::vector< INMOST_DATA_ENUM_TYPE >::const_iterator search = std::lower_bound(NonlocalPost.begin(),NonlocalPost.end(),pos);
                assert(*search == pos); //is there such index?
                return IntervalEnd + static_cast<INMOST_DATA_ENUM_TYPE>(search-NonlocalPost.begin());
            }
            return pos;
        }

        INMOST_DATA_ENUM_TYPE RowMerger::UnmapIndex(INMOST_DATA_ENUM_TYPE pos) const
        {
            if( pos < IntervalBeg )
                return NonlocalPre[pos+NonlocalPre.size()-IntervalBeg];
            else if( pos >= IntervalEnd )
                return NonlocalPost[pos-IntervalEnd];
            return pos;
        }


        void RowMerger::Clear()
        {
            INMOST_DATA_ENUM_TYPE i = LinkedList.begin()->first, j;
            LinkedList.begin()->first = EOL;
            while( i != EOL )
            {
                j = LinkedList[i].first;
                LinkedList[i].first = UNDEF;
                LinkedList[i].second = 0.0;
                i = j;
            }
            Nonzeros = 0;
        }


        void RowMerger::PushRow(INMOST_DATA_REAL_TYPE coef, Row & r, bool PreSortRow)
        {
            if( Sorted && PreSortRow ) std::sort(r.Begin(),r.End());
            PushRow(coef,r);
        }

        void RowMerger::PushRow(INMOST_DATA_REAL_TYPE coef, const Row & r)
        {
            assert(Nonzeros == 0); //Linked list should be empty
            assert(LinkedList.begin()->first == EOL); //again check that list is empty
            INMOST_DATA_ENUM_TYPE index = LinkedList.get_interval_beg();
            Row::const_iterator it = r.Begin(), jt;
            while( it != r.End() )
            {
                INMOST_DATA_ENUM_TYPE pos = it->first;
                pos = MapIndex(pos);
                LinkedList[index].first = pos+1;
                LinkedList[pos+1].first = EOL;
                LinkedList[pos+1].second = it->second*coef;
                index = pos+1;
                ++Nonzeros;
                jt = it;
                ++it;
                assert(!Sorted || it == r.End() || jt->first < it->first);
            }
        }

        void RowMerger::AddRow(INMOST_DATA_REAL_TYPE coef, Row & r, bool PreSortRow)
        {
            if( Sorted && PreSortRow ) std::sort(r.Begin(),r.End());
            AddRow(coef,r);
        }

        void RowMerger::AddRow(INMOST_DATA_REAL_TYPE coef, const Row & r)
        {
            INMOST_DATA_ENUM_TYPE index = LinkedList.get_interval_beg(), next;
            Row::const_iterator it = r.Begin(), jt;
            while( it != r.End() )
            {
                INMOST_DATA_ENUM_TYPE pos = it->first;
                pos = MapIndex(pos);
                if( LinkedList[pos+1].first != UNDEF )
                    LinkedList[pos+1].second += coef*it->second;
                else if( Sorted )
                {
                    next = index;
                    while(next < pos+1)
                    {
                        index = next;
                        next = LinkedList[index].first;
                    }
                    assert(index < pos+1);
                    assert(pos+1 < next);
                    LinkedList[index].first = pos+1;
                    LinkedList[pos+1].first = next;
                    LinkedList[pos+1].second = coef*it->second;
                    ++Nonzeros;
                }
                else
                {
                    LinkedList[pos+1].first = LinkedList[index].first;
                    LinkedList[pos+1].second = coef*it->second;
                    LinkedList[index].first = pos+1;
                    ++Nonzeros;
                }
                jt = it;
                ++it;
                assert(!Sorted || it == r.End() || jt->first < it->first);
            }
        }


        void RowMerger::RetriveRow(Row & r)
        {
            r.Resize(Nonzeros);
            INMOST_DATA_ENUM_TYPE i = LinkedList.begin()->first, k = 0;
            while( i != EOL )
            {
                if( LinkedList[i].second )
                {
                    r.GetIndex(k) = UnmapIndex(i-1);
                    r.GetValue(k) = LinkedList[i].second;
                    ++k;
                }
                i = LinkedList[i].first;
            }
            r.Resize(k);
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
310
////////class HessianRow
Dmitry Bagaev's avatar
Dmitry Bagaev committed
311

312
313
        void   HessianRow::RowVec(INMOST_DATA_REAL_TYPE alpha, const Row & rU, INMOST_DATA_REAL_TYPE beta, Row & rJ) const
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
314
            INMOST_DATA_ENUM_TYPE end = rU.Size();
315
316
317
318
319
320
321
322
323
324
            for(INMOST_DATA_ENUM_TYPE i = 0; i < end; i++) rJ.GetValue(i) *= beta;
            for(INMOST_DATA_ENUM_TYPE i = 0; i < end; i++)
            {
                const index & ind = GetIndex(i);
                if( ind.first == ind.second )
                    rJ[ind.first] += alpha*rU.get_safe(ind.first)*GetValue(i);
                else
                {
                    rJ[ind.first ] += alpha*rU.get_safe(ind.second)*GetValue(i);
                    rJ[ind.second] += alpha*rU.get_safe(ind.first )*GetValue(i);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
325
326
327
328
                }
            }
        }

329
330
331
332
        bool HessianRow::isSorted() const
        {
            for(INMOST_DATA_ENUM_TYPE k = 1; k < Size(); ++k)
                if( GetIndex(k) < GetIndex(k-1) ) return false;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
333
334
335
            return true;
        }

336
337
        void HessianRow::MergeSortedRows(INMOST_DATA_REAL_TYPE alpha, const HessianRow & left, INMOST_DATA_REAL_TYPE beta, const HessianRow & right, HessianRow & output)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
338
339
            assert(left.isSorted());
            assert(right.isSorted());
340
            output.Resize(left.Size()+right.Size());
Dmitry Bagaev's avatar
Dmitry Bagaev committed
341
            INMOST_DATA_ENUM_TYPE i = 0, j = 0, q = 0;
342
343
344
345
            while( i < left.Size() && j < right.Size() )
            {
                if( left.GetIndex(i) < right.GetIndex(j) )
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
346
                    output.GetIndex(q) = left.GetIndex(i);
347
                    output.GetValue(q) = alpha*left.GetValue(i);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
348
349
                    ++q;
                    ++i;
350
351
352
                }
                else if( left.GetIndex(i) == right.GetIndex(j) )
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
353
                    output.GetIndex(q) = left.GetIndex(i);
354
                    output.GetValue(q) = alpha*left.GetValue(i) + beta*right.GetValue(j);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
355
356
357
                    ++q;
                    ++i;
                    ++j;
358
359
                }
                else //right is smaller
Dmitry Bagaev's avatar
Dmitry Bagaev committed
360
361
                {
                    output.GetIndex(q) = right.GetIndex(j);
362
                    output.GetValue(q) = beta*right.GetValue(j);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
363
364
365
366
                    ++q;
                    ++j;
                }
            }
367
368
369
370
371
372
            while( i < left.Size() )
            {
                if( q > 0 && (output.GetIndex(q-1) == left.GetIndex(i)) )
                    output.GetValue(q-1) += alpha*left.GetValue(i);
                else
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
373
                    output.GetIndex(q) = left.GetIndex(i);
374
                    output.GetValue(q) = alpha*left.GetValue(i);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
375
376
377
378
                    ++q;
                }
                ++i;
            }
379
380
381
382
383
384
            while( j < right.Size() )
            {
                if( q > 0 && (output.GetIndex(q-1) == right.GetIndex(j)) )
                    output.GetValue(q-1) += beta*right.GetValue(j);
                else
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
385
                    output.GetIndex(q) = right.GetIndex(j);
386
                    output.GetValue(q) = beta*right.GetValue(j);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
387
388
389
390
391
392
393
                    ++q;
                }
                ++j;
            }
            output.Resize(q);
        }

394
395
        void HessianRow::MergeJacobianHessian(INMOST_DATA_REAL_TYPE a, const Row & JL, const Row & JR, INMOST_DATA_REAL_TYPE b, const HessianRow & HL, INMOST_DATA_REAL_TYPE c, const HessianRow & HR, HessianRow & output)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
396
397
            // merge three sorted arrays at once
            // one of the array is synthesized from JL and JR on the go
398
            static const entry stub_entry = make_entry(make_index(ENUMUNDEF,ENUMUNDEF),0.0);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
399
400
401
402
            assert(JL.isSorted());
            assert(JR.isSorted());
            assert(HL.isSorted());
            assert(HR.isSorted());
403
404
405
406
407
408
409
410
411
412
413
414
415
            output.Resize(HL.Size()+HR.Size()+JL.Size()*JR.Size());
            INMOST_DATA_ENUM_TYPE i = 0, j = 0, k1 = 0, l1 = 0, k2 = 0, l2 = 0, q = 0, kk1 = 0, ll2 = 0, r;
            entry candidate[4] = {stub_entry,stub_entry,stub_entry,stub_entry};
            if( i < HL.Size() )
                candidate[0] = make_entry(HL.GetIndex(i),b*HL.GetValue(i));
            if( j < HR.Size() )
                candidate[1] = make_entry(HR.GetIndex(j),c*HR.GetValue(j));
            if( k1 < JL.Size() && l1 < JR.Size() )
                candidate[2] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),0.5*a*JL.GetValue(k1)*JR.GetValue(l1));
            if( k2 < JL.Size() && l2 < JR.Size() )
                candidate[3] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),0.5*a*JL.GetValue(k2)*JR.GetValue(l2));
            do
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
416
417
                //pick smallest
                r = 0;
418
419
420
                if( candidate[1].first < candidate[r].first ) r = 1;
                if( candidate[2].first < candidate[r].first ) r = 2;
                if( candidate[3].first < candidate[r].first ) r = 3;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
421
                //all candidates are stub - exit
422
                if( candidate[r].first == stub_entry.first ) break;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
423
                //record selected entry
424
425
426
427
                if( q > 0 && (output.GetIndex(q-1) == candidate[r].first) )
                    output.GetValue(q-1) += candidate[r].second;
                else
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
428
429
430
431
                    output.GetIndex(q) = candidate[r].first;
                    output.GetValue(q) = candidate[r].second;
                    ++q;
                }
432
                if( r == 0 ) //update left hessian index
Dmitry Bagaev's avatar
Dmitry Bagaev committed
433
                {
434
435
                    if( ++i < HL.Size() )
                        candidate[0] = make_entry(HL.GetIndex(i),b*HL.GetValue(i));
Dmitry Bagaev's avatar
Dmitry Bagaev committed
436
                    else candidate[0] = stub_entry;
437
438
                }
                else if( r == 1 ) //update right hessian index
Dmitry Bagaev's avatar
Dmitry Bagaev committed
439
                {
440
441
                    if( ++j < HR.Size() )
                        candidate[1] = make_entry(HR.GetIndex(j),c*HR.GetValue(j));
Dmitry Bagaev's avatar
Dmitry Bagaev committed
442
                    else candidate[1] = stub_entry;
443
444
                }
                else if( r == 2 ) //update jacobians indexes
Dmitry Bagaev's avatar
Dmitry Bagaev committed
445
                {
446
                    if( ++k1 == JL.Size() )
Dmitry Bagaev's avatar
Dmitry Bagaev committed
447
                    {
448
449
450
451
452
                        ++l1;
                        if( l1 < JR.Size() )
                        {
                            while(kk1 < JL.Size() && JL.GetIndex(kk1) < JR.GetIndex(l1) ) kk1++;
                            k1 = kk1;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
453
454
                        }
                    }
455
456
                    if( k1 < JL.Size() && l1 < JR.Size() )
                        candidate[2] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),(k1==l1?0.5:1)*a*JL.GetValue(k1)*JR.GetValue(l1));
Dmitry Bagaev's avatar
Dmitry Bagaev committed
457
458
459
                    else
                        candidate[2] = stub_entry;
                }
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
                else //update jacobians indexes
                {
                    if( ++l2 == JR.Size() )
                    {
                        ++k2;
                        if( k2 < JL.Size() )
                        {
                            while(ll2 < JR.Size() && JL.GetIndex(k2) > JR.GetIndex(ll2) ) ll2++;
                            l2 = ll2;
                        }
                    }
                    if( k2 < JL.Size() && l2 < JR.Size() )
                        candidate[3] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),(k2==l2?0.5:1)*a*JL.GetValue(k2)*JR.GetValue(l2));
                    else
                        candidate[3] = stub_entry;
                }
            }
            while(true);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
478
479
480
            output.Resize(q);
        }

481
482
        void HessianRow::MergeJacobianHessian(INMOST_DATA_REAL_TYPE a, const Row & JL, const Row & JR, INMOST_DATA_REAL_TYPE b, const HessianRow & H, HessianRow & output)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
483
484
            // merge three sorted arrays at once
            // one of the array is synthesized from JL and JR on the go
485
            static const entry stub_entry = make_entry(make_index(ENUMUNDEF,ENUMUNDEF),0.0);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
486
487
488
            assert(JL.isSorted());
            assert(JR.isSorted());
            assert(H.isSorted());
489
490
491
492
493
494
495
496
497
498
499
            output.Resize(H.Size()+JL.Size()*JR.Size());
            INMOST_DATA_ENUM_TYPE i = 0, k1 = 0, l1 = 0, q = 0, k2 = 0, l2 = 0, r, ll2 = 0, kk1 = 0;
            entry candidate[3] = {stub_entry,stub_entry,stub_entry};
            if( i < H.Size() )
                candidate[0] = make_entry(H.GetIndex(i),b*H.GetValue(i));
            if( k1 < JL.Size() && l1 < JR.Size() )
                candidate[1] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),0.5*a*JL.GetValue(k1)*JR.GetValue(l1));
            if( k2 < JL.Size() && l2 < JR.Size() )
                candidate[2] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),0.5*a*JL.GetValue(k2)*JR.GetValue(l2));
            do
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
500
501
                //pick smallest
                r = 0;
502
503
                if( candidate[1].first < candidate[r].first ) r = 1;
                if( candidate[2].first < candidate[r].first ) r = 2;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
504
                //all candidates are stub - exit
505
                if( candidate[r].first == stub_entry.first ) break;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
506
                //record selected entry
507
508
509
510
511
512
                if( q > 0 && (output.GetIndex(q-1) == candidate[r].first) )
                {
                    output.GetValue(q-1) += candidate[r].second;
                }
                else
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
513
514
515
516
                    output.GetIndex(q) = candidate[r].first;
                    output.GetValue(q) = candidate[r].second;
                    ++q;
                }
517
                if( r == 0 ) //update left hessian index
Dmitry Bagaev's avatar
Dmitry Bagaev committed
518
                {
519
520
                    if( ++i < H.Size() )
                        candidate[0] = make_entry(H.GetIndex(i),b*H.GetValue(i));
Dmitry Bagaev's avatar
Dmitry Bagaev committed
521
                    else candidate[0] = stub_entry;
522
523
                }
                else if( r == 1 )
Dmitry Bagaev's avatar
Dmitry Bagaev committed
524
                {
525
                    if( ++k1 == JL.Size() )
Dmitry Bagaev's avatar
Dmitry Bagaev committed
526
                    {
527
528
529
530
531
                        ++l1;
                        if( l1 < JR.Size() )
                        {
                            while(kk1 < JL.Size() && JL.GetIndex(kk1) < JR.GetIndex(l1) ) kk1++;
                            k1 = kk1;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
532
533
                        }
                    }
534
535
                    if( k1 < JL.Size() && l1 < JR.Size() )
                        candidate[1] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),(k1==l1?0.5:1)*a*JL.GetValue(k1)*JR.GetValue(l1));
Dmitry Bagaev's avatar
Dmitry Bagaev committed
536
537
538
                    else
                        candidate[1] = stub_entry;
                }
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
                else //update jacobians indexes
                {
                    if( ++l2 == JR.Size() )
                    {
                        ++k2;
                        if( k2 < JL.Size() )
                        {
                            while(ll2 < JR.Size() && JL.GetIndex(k2) > JR.GetIndex(ll2) ) ll2++;
                            l2 = ll2;
                        }
                    }
                    if( k2 < JL.Size() && l2 < JR.Size() )
                        candidate[2] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),(k2==l2?0.5:1)*a*JL.GetValue(k2)*JR.GetValue(l2));
                    else
                        candidate[2] = stub_entry;
                }
            }
            while(true);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
557
558
559
            output.Resize(q);
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
560
561
////////class Row
#if defined(USE_SOLVER)
562
563
        INMOST_DATA_REAL_TYPE   Row::RowVec(Vector & x) const
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
564
565
            INMOST_DATA_REAL_TYPE ret = 0;
            INMOST_DATA_ENUM_TYPE end = Size();
566
            for(INMOST_DATA_ENUM_TYPE i = 0; i < end; i++) ret = ret + x[GetIndex(i)]*GetValue(i);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
567
568
            return ret;
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
569
#endif //USE_SOLVER
570
571
572
573
        bool Row::isSorted() const
        {
            for(INMOST_DATA_ENUM_TYPE k = 1; k < Size(); ++k)
                if( GetIndex(k) < GetIndex(k-1) ) return false;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
574
575
            return true;
        }
576
577
        void Row::MergeSortedRows(INMOST_DATA_REAL_TYPE alpha, const Row & left, INMOST_DATA_REAL_TYPE beta, const Row & right, Row & output)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
578
579
            assert(left.isSorted());
            assert(right.isSorted());
580
            output.Resize(left.Size()+right.Size());
Dmitry Bagaev's avatar
Dmitry Bagaev committed
581
            INMOST_DATA_ENUM_TYPE i = 0, j = 0, q = 0;
582
583
584
585
            while( i < left.Size() && j < right.Size() )
            {
                if( left.GetIndex(i) < right.GetIndex(j) )
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
586
                    output.GetIndex(q) = left.GetIndex(i);
587
                    output.GetValue(q) = alpha*left.GetValue(i);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
588
589
                    ++q;
                    ++i;
590
591
592
                }
                else if( left.GetIndex(i) == right.GetIndex(j) )
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
593
                    output.GetIndex(q) = left.GetIndex(i);
594
                    output.GetValue(q) = alpha*left.GetValue(i) + beta*right.GetValue(j);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
595
596
597
                    ++q;
                    ++i;
                    ++j;
598
599
                }
                else //right is smaller
Dmitry Bagaev's avatar
Dmitry Bagaev committed
600
601
                {
                    output.GetIndex(q) = right.GetIndex(j);
602
                    output.GetValue(q) = beta*right.GetValue(j);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
603
604
605
606
                    ++q;
                    ++j;
                }
            }
607
608
609
610
611
612
            while( i < left.Size() )
            {
                if( q > 0 && (output.GetIndex(q-1) == left.GetIndex(i)) )
                    output.GetValue(q-1) += alpha*left.GetValue(i);
                else
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
613
                    output.GetIndex(q) = left.GetIndex(i);
614
                    output.GetValue(q) = alpha*left.GetValue(i);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
615
616
617
618
                    ++q;
                }
                ++i;
            }
619
620
621
622
623
624
            while( j < right.Size() )
            {
                if( q > 0 && (output.GetIndex(q-1) == right.GetIndex(j)) )
                    output.GetValue(q-1) += beta*right.GetValue(j);
                else
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
625
                    output.GetIndex(q) = right.GetIndex(j);
626
                    output.GetValue(q) = beta*right.GetValue(j);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
627
628
629
630
631
632
                    ++q;
                }
                ++j;
            }
            output.Resize(q);
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
633
#endif //defined(USE_SOLVER) || defined(USE_AUTODIFF)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
634
635


Kirill Terekhov's avatar
Kirill Terekhov committed
636
637
#if defined(USE_SOLVER)
////////class Vector
638
639
        Vector::Vector(std::string _name, INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end, INMOST_MPI_Comm _comm) :data(start,end)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
640
641
642
643
644
            comm = _comm;
            name = _name;
            is_parallel = false;
        }

645
646
        Vector::Vector(const Vector & other) : data(other.data)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
647
648
649
650
651
            comm = other.comm;
            name = other.name;
            is_parallel = other.is_parallel;
        }

652
653
        Vector & Vector::operator =(Vector const & other)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
654
655
656
657
658
659
660
            comm = other.comm;
            data = other.data;
            name = other.name;
            is_parallel = other.is_parallel;
            return *this;
        }

661
662
        Vector::~Vector()
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
663
664
        }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
665
        void     Vector::Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg, INMOST_DATA_ENUM_TYPE mend, std::string file_ord)
666
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
667
668
            char str[16384];
            std::ifstream input(file.c_str());
669
            if( input.fail() ) throw -1;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
670
671
672
673
            int state = 0, k;
            INMOST_DATA_ENUM_TYPE vec_size, vec_block, ind = 0;
            INMOST_DATA_REAL_TYPE val;
            int size = 1, rank = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
674
#if defined(USE_MPI)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
675
676
            int flag = 0;
            MPI_Initialized(&flag);
677
678
679
680
            if( flag && mend == ENUMUNDEF && mbeg == ENUMUNDEF )
            {
                MPI_Comm_rank(GetCommunicator(),&rank);
                MPI_Comm_size(GetCommunicator(),&size);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
681
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
682
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
            int * ord = NULL;
            if (file_ord != "")
            {
                std::ifstream input_ord;
                input_ord.open(file_ord.c_str(), std::ifstream::in);
                if( input_ord.fail() ) throw -2;
                int n;
                input_ord >> n;
                ord = (int *) malloc(sizeof(int) * n);
                for (int i=0; i<n; i++) input_ord >> ord[i];
                int nbl;
                input_ord >> nbl;
                if( nbl != size ) throw -3;
                int * ibl;
                ibl = (int *) malloc(sizeof(int) * (nbl+1));
                for (int i=0; i<nbl+1; i++) input_ord >> ibl[i];
                if( mbeg == ENUMUNDEF ) mbeg = ibl[rank];
                if( mend == ENUMUNDEF ) mend = ibl[rank+1];
                for (int i=0; i<n; i++) ord[i] -= 1;
                free(ibl);
                input_ord.close();
                //std::cout<<"Vector::Load(): n="<<n<<" nbl="<<nbl<<" np="<<size<<" id="<<rank<<" mbeg="<<mbeg<<" mend="<<mend<<std::endl;//db
            }

707
708
709
710
711
712
713
            while( !input.getline(str,16384).eof() )
            {
                k = 0; while( isspace(str[k]) ) k++;
                if( str[k] == '%' || str[k] == '\0' ) continue;
                std::istringstream istr(str+k);
                switch(state)
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
714
                    case 0:
715
716
717
718
719
720
721
                        istr >> vec_size; state = 1;
                        vec_block = vec_size/size;
                        if( mbeg == ENUMUNDEF ) mbeg = rank*vec_block;
                        if( mend == ENUMUNDEF )
                        {
                            if( rank == size-1 ) mend = vec_size;
                            else mend = mbeg+vec_block;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
722
                        }
723
                        SetInterval(mbeg,mend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
724
725
726
                        break;
                    case 1:
                        istr >> val;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
727
728
729
730
731
                        if( ord ) {
                            if( ord[ind] >= mbeg && ord[ind] < mend ) data[ord[ind]] = val;
                        } else {
                            if( ind >= mbeg && ind < mend ) data[ind] = val;
                        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
732
733
734
735
736
                        ind++;
                        break;
                }
            }
            input.close();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
737
            if (file_ord != "") free(ord);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
738
739
740
        }


741
742
        void     Vector::Save(std::string file)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
743
744
            INMOST_DATA_ENUM_TYPE vecsize = Size();

Kirill Terekhov's avatar
Kirill Terekhov committed
745
#if defined(USE_MPI)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
746
747
            int rank = 0, size = 1;
            {
748
749
                MPI_Comm_rank(GetCommunicator(),&rank);
                MPI_Comm_size(GetCommunicator(),&size);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
750
                INMOST_DATA_ENUM_TYPE temp = vecsize;
751
                MPI_Allreduce(&temp,&vecsize,1,INMOST_MPI_DATA_ENUM_TYPE,MPI_SUM,GetCommunicator());
Dmitry Bagaev's avatar
Dmitry Bagaev committed
752
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
753
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
754
755
756
            std::stringstream rhs(std::ios::in | std::ios::out);
            rhs << std::scientific;
            rhs.precision(15);
757
            for(iterator it = Begin(); it != End(); ++it) rhs << *it << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
758
#if defined(USE_MPI) && defined(USE_MPI_FILE) // Use mpi files
Dmitry Bagaev's avatar
Dmitry Bagaev committed
759
760
761
762
            {
                int ierr;
                MPI_File fh;
                MPI_Status stat;
763
764
                ierr = MPI_File_open(GetCommunicator(),const_cast<char *>(file.c_str()), MPI_MODE_CREATE | MPI_MODE_DELETE_ON_CLOSE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh);
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
765
                ierr = MPI_File_close(&fh);
766
767
768
769
770
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
                ierr = MPI_File_open(GetCommunicator(),const_cast<char *>(file.c_str()),MPI_MODE_WRONLY | MPI_MODE_CREATE,MPI_INFO_NULL,&fh);
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
                if( rank == 0 )
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
771
772
773
774
775
                    std::stringstream header;
                    //header << "% vector " << name << std::endl;
                    //header << "% is written by INMOST" << std::endl;
                    //header << "% by MPI_File_* api" << std::endl;
                    header << vecsize << std::endl;
776
777
                    ierr = MPI_File_write_shared(fh,const_cast<char *>(header.str().c_str()),static_cast<int>(header.str().size()),MPI_CHAR,&stat);
                    if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
778
                }
779
780
                ierr = MPI_File_write_ordered(fh,const_cast<char *>(rhs.str().c_str()),static_cast<int>(rhs.str().size()),MPI_CHAR,&stat);
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
781
                ierr = MPI_File_close(&fh);
782
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
783
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
784
#elif defined(USE_MPI) //USE_MPI alternative
Dmitry Bagaev's avatar
Dmitry Bagaev committed
785
            std::string senddata = rhs.str(), recvdata;
Kirill Terekhov's avatar
Kirill Terekhov committed
786
787
788
789
790
791
			int sendsize = static_cast<int>(senddata.size());
			std::vector<int> recvsize(size), displ(size);
			MPI_Gather(&sendsize,1,MPI_INT,&recvsize[0],1,MPI_INT,0,GetCommunicator());
			if( rank == 0 )
			{
				int totsize = recvsize[0];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
792

Kirill Terekhov's avatar
Kirill Terekhov committed
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
				displ[0] = 0;
				for(int i = 1; i < size; i++)
				{
					totsize += recvsize[i];
					displ[i] = displ[i-1]+recvsize[i-1];
				}
				recvdata.resize(totsize);
			}
			else recvdata.resize(1); //protect from dereferencing null
			MPI_Gatherv(&senddata[0],sendsize,MPI_CHAR,&recvdata[0],&recvsize[0],&displ[0],MPI_CHAR,0,GetCommunicator());
			if( rank == 0 )
			{
				std::fstream output(file.c_str(),std::ios::out);
				output << vecsize << std::endl;
				output << recvdata;
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
809
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
810
811
812
813
814
815
			std::fstream output(file.c_str(),std::ios::out);
			//output << "% vector " << name << std::endl;
			//output << "% is written by INMOST" << std::endl;
			//output << "% by sequential write" << std::endl;
			output << vecsize << std::endl;
			output << rhs.rdbuf();
Kirill Terekhov's avatar
Kirill Terekhov committed
816
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
817
818
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
819
////////class Matrix
820
        void Matrix::MatVec(INMOST_DATA_REAL_TYPE alpha, Vector & x, INMOST_DATA_REAL_TYPE beta, Vector & out) const //y = alpha*A*x + beta * y
Dmitry Bagaev's avatar
Dmitry Bagaev committed
821
822
823
        {
            INMOST_DATA_ENUM_TYPE mbeg, mend;
            INMOST_DATA_INTEGER_TYPE ind, imbeg, imend;
824
825
826
827
828
            if( out.Empty() )
            {
                INMOST_DATA_ENUM_TYPE vbeg,vend;
                GetInterval(vbeg,vend);
                out.SetInterval(vbeg,vend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
829
830
831
832
            }
            //CHECK SOMEHOW FOR DEBUG THAT PROVIDED VECTORS ARE OK
            //~ assert(GetFirstIndex() == out.GetFirstIndex());
            //~ assert(Size() == out.Size());
833
            GetInterval(mbeg,mend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
834
835
            imbeg = mbeg;
            imend = mend;
Kirill Terekhov's avatar
Kirill Terekhov committed
836
837
838
#if defined(USE_OMP)
#pragma omp for private(ind)
#endif
839
            for(ind = imbeg; ind < imend; ++ind) //iterate rows of matrix
Dmitry Bagaev's avatar
Dmitry Bagaev committed
840
841
842
843
844
                out[ind] = beta * out[ind] + alpha * (*this)[ind].RowVec(x);
            // outer procedure should update out vector, if needed
        }


845
        void Matrix::MatVecTranspose(INMOST_DATA_REAL_TYPE alpha, Vector & x, INMOST_DATA_REAL_TYPE beta, Vector & out) const //y = alpha*A*x + beta * y
Dmitry Bagaev's avatar
Dmitry Bagaev committed
846
847
848
        {
            INMOST_DATA_ENUM_TYPE mbeg, mend;
            INMOST_DATA_INTEGER_TYPE ind, imbeg, imend;
849
850
851
852
853
            if( out.Empty() )
            {
                INMOST_DATA_ENUM_TYPE vbeg,vend;
                GetInterval(vbeg,vend);
                out.SetInterval(vbeg,vend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
854
855
856
857
            }
            //CHECK SOMEHOW FOR DEBUG THAT PROVIDED VECTORS ARE OK
            //~ assert(GetFirstIndex() == out.GetFirstIndex());
            //~ assert(Size() == out.Size());
858
            GetInterval(mbeg,mend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
859
860
            imbeg = mbeg;
            imend = mend;
861
            if( beta ) for(Vector::iterator it = out.Begin(); it != out.End(); ++it) (*it) *= beta;
Kirill Terekhov's avatar
Kirill Terekhov committed
862
863
864
#if defined(USE_OMP)
#pragma omp for private(ind)
#endif
865
866
867
            for(ind = imbeg; ind < imend; ++ind)
            {
                for(Row::const_iterator it = (*this)[ind].Begin(); it != (*this)[ind].End(); ++it)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
868
869
870
871
872
873
874
                    out[it->first] += alpha * x[ind] * it->second;
            }
            // outer procedure should update out vector, if needed
        }


        Matrix::Matrix(std::string _name, INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end, INMOST_MPI_Comm _comm)
875
876
                :data(start,end)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
877
878
            is_parallel = false;
            comm = _comm;
879
            SetInterval(start,end);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
880
881
882
            name = _name;
        }

883
884
        Matrix::Matrix(const Matrix & other) :data(other.data)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
885
886
887
888
889
            comm = other.comm;
            name = other.name;
        }


890
891
        Matrix & Matrix::operator =(Matrix const & other)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
892
893
894
895
896
897
898
899
            comm = other.comm;
            data = other.data;
            name = other.name;
            return *this;
        }

        Matrix::~Matrix() {}

900
901
        void      Matrix::MoveRows(INMOST_DATA_ENUM_TYPE from, INMOST_DATA_ENUM_TYPE to, INMOST_DATA_ENUM_TYPE size)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
902
            INMOST_DATA_ENUM_TYPE i = to + size, j = from + size;
903
904
            if( size > 0 && to != from )
                while( j != from ) data[--j].MoveRow(data[--i]);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
905
906
        }

907
908
        void      HessianMatrix::MoveRows(INMOST_DATA_ENUM_TYPE from, INMOST_DATA_ENUM_TYPE to, INMOST_DATA_ENUM_TYPE size)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
909
            INMOST_DATA_ENUM_TYPE i = to + size, j = from + size;
910
911
            if( size > 0 && to != from )
                while( j != from ) data[--j].MoveRow(data[--i]);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
912
913
914
        }


915
        void     Matrix::Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg, INMOST_DATA_ENUM_TYPE mend, std::string file_ord)
916
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
917
918
            char str[16384];
            std::ifstream input(file.c_str());
919
            if( input.fail() ) throw -1;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
920
921
922
923
            int state = 0, k;
            INMOST_DATA_ENUM_TYPE mat_size, max_lines, row, col, mat_block;
            INMOST_DATA_REAL_TYPE val;
            int size = 1, rank = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
924
#if defined(USE_MPI)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
925
926
            int flag = 0;
            MPI_Initialized(&flag);
927
928
929
930
            if( flag && mend == ENUMUNDEF && mbeg == ENUMUNDEF )
            {
                MPI_Comm_rank(GetCommunicator(),&rank);
                MPI_Comm_size(GetCommunicator(),&size);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
931
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
932
#endif
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
            int * ord = NULL;
            if (file_ord != "")
            {
                std::ifstream input_ord;
                input_ord.open(file_ord.c_str(), std::ifstream::in);
                if( input_ord.fail() ) throw -2;
                int n;
                input_ord >> n; // check if( n == line )
                ord = (int *) malloc(sizeof(int) * n);
                for (int i=0; i<n; i++) input_ord >> ord[i];
                int nbl;
                input_ord >> nbl;
                if( nbl != size ) throw -3;
                int * ibl;
                ibl = (int *) malloc(sizeof(int) * (nbl+1));
                for (int i=0; i<nbl+1; i++) input_ord >> ibl[i];
                if( mbeg == ENUMUNDEF ) mbeg = (unsigned int) ibl[rank];
                if( mend == ENUMUNDEF ) mend = (unsigned int) ibl[rank + 1];
                for (int i=0; i<n; i++) ord[i] -= 1;
                free(ibl);
                input_ord.close();
                //std::cout<<"Matrix::Load(): n="<<n<<" nbl="<<nbl<<" np="<<size<<" id="<<rank<<" mbeg="<<mbeg<<" mend="<<mend<<std::endl;//db
            }

Dmitry Bagaev's avatar
Dmitry Bagaev committed
957
            int line = 0;
958
959
            while( !input.getline(str,16384).eof() )
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
960
                line++;
961
962
963
964
965
                k = 0; while( isspace(str[k]) ) k++;
                if( str[k] == '%' || str[k] == '\0' ) continue;
                std::istringstream istr(str+k);
                switch(state)
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
966
                    case 0:
967
968
969
970
971
972
973
                        istr >> mat_size >> mat_size >> max_lines; state = 1;
                        mat_block = mat_size/size;
                        if( mbeg == ENUMUNDEF ) mbeg = rank*mat_block;
                        if( mend == ENUMUNDEF )
                        {
                            if( rank == size-1 ) mend = mat_size;
                            else mend = mbeg+mat_block;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
974
                        }
975
                        SetInterval(mbeg,mend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
976
977
978
979
                        //~ std::cout << rank << " my interval " << mbeg << ":" << mend << std::endl;
                        break;
                    case 1:
                        istr >> row >> col >> val;
980
                        row--; col--;
981
                        if( ord ) { row = ord[row]; col = ord[col]; }
982
                        if( row >= mbeg && row < mend ) data[row][col] = val;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
983
984
985
986
                        break;
                }
            }
            int nonzero = 0;
987
            for(iterator it = Begin(); it != End(); ++it) nonzero += it->Size();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
988
989
            //~ std::cout << rank << " total nonzero " << max_lines << " my nonzero " << nonzero << std::endl;
            input.close();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
990
            if (file_ord != "") free(ord);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
991
992
        }

993
994
995
        void     Matrix::Save(std::string file, const AnnotationService * text)
        {
            INMOST_DATA_ENUM_TYPE matsize = Size(), nonzero = 0, row = GetFirstIndex()+1;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
996
            bool have_annotation = false;
997
998
999
            if( text && !text->Empty() )
            {
                if( text->GetFirstIndex() == GetFirstIndex() &&
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1000
1001
                    text->GetLastIndex() == GetLastIndex())
                    have_annotation = true;
1002
1003
1004
1005
                else
                {
                    std::cout << "Size of provided annotation (" << text->GetFirstIndex() << "," << text->GetLastIndex() << ")" << std::endl;
                    std::cout << "differs from the size of the matrix (" << GetFirstIndex() << "," << GetLastIndex() << ")" << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1006
1007
                }
            }
1008
            for(iterator it = Begin(); it != End(); ++it) nonzero += it->Size();
Kirill Terekhov's avatar
Kirill Terekhov committed
1009
#if defined(USE_MPI)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1010
1011
            int rank = 0, size = 1;
            {
1012
1013
1014
1015
                MPI_Comm_rank(GetCommunicator(),&rank);
                MPI_Comm_size(GetCommunicator(),&size);
                INMOST_DATA_ENUM_TYPE temp_two[2] = {matsize,nonzero}, two[2];
                MPI_Allreduce(temp_two,two,2,INMOST_MPI_DATA_ENUM_TYPE,MPI_SUM,GetCommunicator());
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1016
1017
1018
                matsize = two[0];
                nonzero = two[1];
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
1019
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1020
1021
1022
            std::stringstream mtx(std::ios::in | std::ios::out);
            mtx << std::scientific;
            mtx.precision(15);
1023
1024
1025
1026
1027
1028
            for(INMOST_DATA_ENUM_TYPE k = GetFirstIndex(); k < GetLastIndex(); ++k)
            {
                if( have_annotation )
                {
                    const std::string & str = text->GetAnnotation(k);
                    if( !str.empty() ) mtx << "% " << str << "\n";
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1029
                }
1030
1031
                for(Row::iterator jt = (*this)[k].Begin(); jt != (*this)[k].End(); ++jt)
                    mtx << row << " " << jt->first+1 << " " << jt->second << "\n";
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1032
1033
                ++row;
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
1034
#if defined(USE_MPI) && defined(USE_MPI_FILE) // USE_MPI2?
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1035
1036
1037
1038
1039
            {
                char estring[MPI_MAX_ERROR_STRING];
                int ierr, len;
                MPI_File fh;
                MPI_Status stat;
1040
1041
                ierr = MPI_File_open(GetCommunicator(),const_cast<char *>(file.c_str()), MPI_MODE_CREATE | MPI_MODE_DELETE_ON_CLOSE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh);
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1042
                ierr = MPI_File_close(&fh);
1043
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1044
                ierr = MPI_Barrier(GetCommunicator());
1045
1046
1047
1048
1049
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
                ierr = MPI_File_open(GetCommunicator(),const_cast<char *>(file.c_str()),MPI_MODE_WRONLY | MPI_MODE_CREATE,MPI_INFO_NULL,&fh);
                if( ierr != MPI_SUCCESS )
                {
                    MPI_Error_string(ierr,estring,&len);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1050
                    std::cout << estring << std::endl;
1051
                    MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1052
                }
1053
1054
                if( rank == 0 )
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1055
1056
1057
1058
1059
1060
1061
                    std::stringstream header;
                    header << "%%MatrixMarket matrix coordinate real general" << std::endl;
                    header << "% matrix " << name << std::endl;
                    header << "% is written by INMOST" << std::endl;
                    header << "% by MPI_File_* api" << std::endl;
                    header << matsize << " " << matsize << " " << nonzero << std::endl;
                    //std::string header_data(header.str());
1062
1063
                    ierr = MPI_File_write_shared(fh,const_cast<char *>(header.str().c_str()),static_cast<int>(header.str().size()),MPI_CHAR,&stat);
                    if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1064
                }
1065
1066
                ierr = MPI_File_write_ordered(fh,const_cast<char *>(mtx.str().c_str()),static_cast<int>(mtx.str().size()),MPI_CHAR,&stat);
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1067
                ierr = MPI_File_close(&fh);
1068
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1069
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
1070
#elif defined(USE_MPI)//USE_MPI alternative
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1071
            std::string senddata = mtx.str(), recvdata;
Kirill Terekhov's avatar
Kirill Terekhov committed
1072
1073
1074
1075
1076
1077
			int sendsize = static_cast<int>(senddata.size());
			std::vector<int> recvsize(size), displ(size);
			MPI_Gather(&sendsize,1,MPI_INT,&recvsize[0],1,MPI_INT,0,GetCommunicator());
			if( rank == 0 )
			{
				int totsize = recvsize[0];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1078

Kirill Terekhov's avatar
Kirill Terekhov committed
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
				displ[0] = 0;
				for(int i = 1; i < size; i++)
				{
					totsize += recvsize[i];
					displ[i] = displ[i-1]+recvsize[i-1];
				}
				recvdata.resize(totsize);
			}
			else recvdata.resize(1); //protect from dereferencing null
			MPI_Gatherv(&senddata[0],sendsize,MPI_CHAR,&recvdata[0],&recvsize[0],&displ[0],MPI_CHAR,0,GetCommunicator());
			if( rank == 0 )
			{
				std::fstream output(file.c_str(),std::ios::out);
				output << "%%MatrixMarket matrix coordinate real general" << std::endl;
				output << "% matrix " << name << std::endl;
				output << "% is written by INMOST" << std::endl;
				output << "% by MPI_Gather* api and sequential write" << std::endl;
				output << matsize << " " << matsize << " " << nonzero << std::endl;
				output << recvdata;
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
1099
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
1100
1101
1102
1103
1104
1105
1106
			std::fstream output(file.c_str(),std::ios::out);
			output << "%%MatrixMarket matrix coordinate real general" << std::endl;
			output << "% matrix " << name << std::endl;
			output << "% is written by INMOST" << std::endl;
			output << "% by sequential write " << std::endl;
			output << matsize << " " << matsize << " " << nonzero << std::endl;
			output << mtx.rdbuf();
Kirill Terekhov's avatar
Kirill Terekhov committed
1107
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1108
1109
        }

1110
1111
        void Matrix::Swap(Matrix & other)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1112
1113
            data.swap(other.data);
            name.swap(other.name);
1114
1115
            std::swap(comm,other.comm);
            std::swap(is_parallel,other.is_parallel);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1116
1117
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
1118
#if defined(USE_OMP)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1119
        bool LockService::HaveLocks() const
Kirill Terekhov's avatar
Kirill Terekhov committed
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
		{
			return !locks.empty();
		}
		bool LockService::Lock(INMOST_DATA_ENUM_TYPE row)
		{
			assert( !locks.empty() );
			omp_set_lock(&locks[row]);
			return true;
		}
		bool LockService::TestLock(INMOST_DATA_ENUM_TYPE row)
		{
			assert( !locks.empty() );
			if( omp_test_lock(&locks[row]) )
				return true;
			else
				return false;
		}
		bool LockService::UnLock(INMOST_DATA_ENUM_TYPE row)
		{
			assert( !locks.empty() );
			omp_unset_lock(&locks[row]);
			return true;
		}
		void LockService::SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end)
		{
			if( !locks.empty() ) DestroyLocks();
			if( end != beg )
			{
				locks.set_interval_beg(beg);
				locks.set_interval_end(end);
				for(INMOST_DATA_ENUM_TYPE k = beg; k < end; ++k)
					omp_init_lock(&locks[k]);
			}
		}
		void LockService::DestroyLocks()
		{
			INMOST_DATA_ENUM_TYPE kbeg,kend;
			kbeg = locks.get_interval_beg();
			kend = locks.get_interval_end();
			for(INMOST_DATA_ENUM_TYPE k = kbeg; k < kend; ++k)
				omp_destroy_lock(&locks[k]);
		}
		INMOST_DATA_ENUM_TYPE  LockService::GetFirstIndex() const {return locks.get_interval_beg();}
		INMOST_DATA_ENUM_TYPE  LockService::GetLastIndex() const {return locks.get_interval_end();}
		bool LockService::Empty() const {return locks.empty();}
		void LockService::GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const {start = locks.get_interval_beg(); end = locks.get_interval_end();}
#else
1167
        bool LockService::HaveLocks() const {return false;}