sparse.cpp 55.6 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
        }

665
666
        void     Vector::Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg, INMOST_DATA_ENUM_TYPE mend)
        {
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
683
684
685
686
687
688
689
            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
690
                    case 0:
691
692
693
694
695
696
697
                        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
698
                        }
699
                        SetInterval(mbeg,mend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
700
701
702
                        break;
                    case 1:
                        istr >> val;
703
                        if( ind >= mbeg && ind < mend ) data[ind] = val;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
704
705
706
707
708
709
710
711
                        ind++;
                        break;
                }
            }
            input.close();
        }


712
713
        void     Vector::Save(std::string file)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
714
715
            INMOST_DATA_ENUM_TYPE vecsize = Size();

Kirill Terekhov's avatar
Kirill Terekhov committed
716
#if defined(USE_MPI)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
717
718
            int rank = 0, size = 1;
            {
719
720
                MPI_Comm_rank(GetCommunicator(),&rank);
                MPI_Comm_size(GetCommunicator(),&size);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
721
                INMOST_DATA_ENUM_TYPE temp = vecsize;
722
                MPI_Allreduce(&temp,&vecsize,1,INMOST_MPI_DATA_ENUM_TYPE,MPI_SUM,GetCommunicator());
Dmitry Bagaev's avatar
Dmitry Bagaev committed
723
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
724
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
725
726
727
            std::stringstream rhs(std::ios::in | std::ios::out);
            rhs << std::scientific;
            rhs.precision(15);
728
            for(iterator it = Begin(); it != End(); ++it) rhs << *it << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
729
#if defined(USE_MPI) && defined(USE_MPI_FILE) // Use mpi files
Dmitry Bagaev's avatar
Dmitry Bagaev committed
730
731
732
733
            {
                int ierr;
                MPI_File fh;
                MPI_Status stat;
734
735
                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
736
                ierr = MPI_File_close(&fh);
737
738
739
740
741
                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
742
743
744
745
746
                    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;
747
748
                    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
749
                }
750
751
                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
752
                ierr = MPI_File_close(&fh);
753
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
754
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
755
#elif defined(USE_MPI) //USE_MPI alternative
Dmitry Bagaev's avatar
Dmitry Bagaev committed
756
            std::string senddata = rhs.str(), recvdata;
Kirill Terekhov's avatar
Kirill Terekhov committed
757
758
759
760
761
762
			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
763

Kirill Terekhov's avatar
Kirill Terekhov committed
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
				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
780
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
781
782
783
784
785
786
			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
787
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
788
789
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
790
////////class Matrix
791
        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
792
793
794
        {
            INMOST_DATA_ENUM_TYPE mbeg, mend;
            INMOST_DATA_INTEGER_TYPE ind, imbeg, imend;
795
796
797
798
799
            if( out.Empty() )
            {
                INMOST_DATA_ENUM_TYPE vbeg,vend;
                GetInterval(vbeg,vend);
                out.SetInterval(vbeg,vend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
800
801
802
803
            }
            //CHECK SOMEHOW FOR DEBUG THAT PROVIDED VECTORS ARE OK
            //~ assert(GetFirstIndex() == out.GetFirstIndex());
            //~ assert(Size() == out.Size());
804
            GetInterval(mbeg,mend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
805
806
            imbeg = mbeg;
            imend = mend;
Kirill Terekhov's avatar
Kirill Terekhov committed
807
808
809
#if defined(USE_OMP)
#pragma omp for private(ind)
#endif
810
            for(ind = imbeg; ind < imend; ++ind) //iterate rows of matrix
Dmitry Bagaev's avatar
Dmitry Bagaev committed
811
812
813
814
815
                out[ind] = beta * out[ind] + alpha * (*this)[ind].RowVec(x);
            // outer procedure should update out vector, if needed
        }


816
        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
817
818
819
        {
            INMOST_DATA_ENUM_TYPE mbeg, mend;
            INMOST_DATA_INTEGER_TYPE ind, imbeg, imend;
820
821
822
823
824
            if( out.Empty() )
            {
                INMOST_DATA_ENUM_TYPE vbeg,vend;
                GetInterval(vbeg,vend);
                out.SetInterval(vbeg,vend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
825
826
827
828
            }
            //CHECK SOMEHOW FOR DEBUG THAT PROVIDED VECTORS ARE OK
            //~ assert(GetFirstIndex() == out.GetFirstIndex());
            //~ assert(Size() == out.Size());
829
            GetInterval(mbeg,mend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
830
831
            imbeg = mbeg;
            imend = mend;
832
            if( beta ) for(Vector::iterator it = out.Begin(); it != out.End(); ++it) (*it) *= beta;
Kirill Terekhov's avatar
Kirill Terekhov committed
833
834
835
#if defined(USE_OMP)
#pragma omp for private(ind)
#endif
836
837
838
            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
839
840
841
842
843
844
845
                    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)
846
847
                :data(start,end)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
848
849
            is_parallel = false;
            comm = _comm;
850
            SetInterval(start,end);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
851
852
853
            name = _name;
        }

854
855
        Matrix::Matrix(const Matrix & other) :data(other.data)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
856
857
858
859
860
            comm = other.comm;
            name = other.name;
        }


861
862
        Matrix & Matrix::operator =(Matrix const & other)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
863
864
865
866
867
868
869
870
            comm = other.comm;
            data = other.data;
            name = other.name;
            return *this;
        }

        Matrix::~Matrix() {}

871
872
        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
873
            INMOST_DATA_ENUM_TYPE i = to + size, j = from + size;
874
875
            if( size > 0 && to != from )
                while( j != from ) data[--j].MoveRow(data[--i]);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
876
877
        }

878
879
        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
880
            INMOST_DATA_ENUM_TYPE i = to + size, j = from + size;
881
882
            if( size > 0 && to != from )
                while( j != from ) data[--j].MoveRow(data[--i]);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
883
884
885
        }


886
        void     Matrix::Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg, INMOST_DATA_ENUM_TYPE mend, std::string file_ord)
887
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
888
889
            char str[16384];
            std::ifstream input(file.c_str());
890
            if( input.fail() ) throw -1;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
891
892
893
894
            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
895
#if defined(USE_MPI)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
896
897
            int flag = 0;
            MPI_Initialized(&flag);
898
899
900
901
            if( flag && mend == ENUMUNDEF && mbeg == ENUMUNDEF )
            {
                MPI_Comm_rank(GetCommunicator(),&rank);
                MPI_Comm_size(GetCommunicator(),&size);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
902
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
903
#endif
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
            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
928
            int line = 0;
929
930
            while( !input.getline(str,16384).eof() )
            {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
931
                line++;
932
933
934
935
936
                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
937
                    case 0:
938
939
940
941
942
943
944
                        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
945
                        }
946
                        SetInterval(mbeg,mend);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
947
948
949
950
                        //~ std::cout << rank << " my interval " << mbeg << ":" << mend << std::endl;
                        break;
                    case 1:
                        istr >> row >> col >> val;
951
                        row--; col--;
952
                        if( ord ) { row = ord[row]; col = ord[col]; }
953
                        if( row >= mbeg && row < mend ) data[row][col] = val;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
954
955
956
957
                        break;
                }
            }
            int nonzero = 0;
958
            for(iterator it = Begin(); it != End(); ++it) nonzero += it->Size();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
959
960
961
962
            //~ std::cout << rank << " total nonzero " << max_lines << " my nonzero " << nonzero << std::endl;
            input.close();
        }

963
964
965
        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
966
            bool have_annotation = false;
967
968
969
            if( text && !text->Empty() )
            {
                if( text->GetFirstIndex() == GetFirstIndex() &&
Dmitry Bagaev's avatar
Dmitry Bagaev committed
970
971
                    text->GetLastIndex() == GetLastIndex())
                    have_annotation = true;
972
973
974
975
                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
976
977
                }
            }
978
            for(iterator it = Begin(); it != End(); ++it) nonzero += it->Size();
Kirill Terekhov's avatar
Kirill Terekhov committed
979
#if defined(USE_MPI)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
980
981
            int rank = 0, size = 1;
            {
982
983
984
985
                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
986
987
988
                matsize = two[0];
                nonzero = two[1];
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
989
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
990
991
992
            std::stringstream mtx(std::ios::in | std::ios::out);
            mtx << std::scientific;
            mtx.precision(15);
993
994
995
996
997
998
            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
999
                }
1000
1001
                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
1002
1003
                ++row;
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
1004
#if defined(USE_MPI) && defined(USE_MPI_FILE) // USE_MPI2?
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1005
1006
1007
1008
1009
            {
                char estring[MPI_MAX_ERROR_STRING];
                int ierr, len;
                MPI_File fh;
                MPI_Status stat;
1010
1011
                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
1012
                ierr = MPI_File_close(&fh);
1013
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1014
                ierr = MPI_Barrier(GetCommunicator());
1015
1016
1017
1018
1019
                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
1020
                    std::cout << estring << std::endl;
1021
                    MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1022
                }
1023
1024
                if( rank == 0 )
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1025
1026
1027
1028
1029
1030
1031
                    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());
1032
1033
                    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
1034
                }
1035
1036
                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
1037
                ierr = MPI_File_close(&fh);
1038
                if( ierr != MPI_SUCCESS ) MPI_Abort(GetCommunicator(),__LINE__);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1039
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
1040
#elif defined(USE_MPI)//USE_MPI alternative
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1041
            std::string senddata = mtx.str(), recvdata;
Kirill Terekhov's avatar
Kirill Terekhov committed
1042
1043
1044
1045
1046
1047
			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
1048

Kirill Terekhov's avatar
Kirill Terekhov committed
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
				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
1069
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
1070
1071
1072
1073
1074
1075
1076
			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
1077
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1078
1079
        }

1080
1081
        void Matrix::Swap(Matrix & other)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1082
1083
            data.swap(other.data);
            name.swap(other.name);
1084
1085
            std::swap(comm,other.comm);
            std::swap(is_parallel,other.is_parallel);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1086
1087
        }

Kirill Terekhov's avatar
Kirill Terekhov committed
1088
#if defined(USE_OMP)
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1089
        bool LockService::HaveLocks() const
Kirill Terekhov's avatar
Kirill Terekhov committed
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
		{
			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
1137
        bool LockService::HaveLocks() const {return false;}
1138
1139
1140
1141
        bool LockService::Lock(INMOST_DATA_ENUM_TYPE row) {(void)row; return true;}
        bool LockService::TestLock(INMOST_DATA_ENUM_TYPE row) {(void)row; return true;}
        bool LockService::UnLock(INMOST_DATA_ENUM_TYPE row) {(void)row; return true;}
        void LockService::SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end) {(void)beg; (void)end;}
1142
1143
1144
1145
1146
        void LockService::DestroyLocks() {}
        INMOST_DATA_ENUM_TYPE  LockService::GetFirstIndex() const {return 0;}
        INMOST_DATA_ENUM_TYPE  LockService::GetLastIndex() const {return 0;}
        bool LockService::Empty() const {return true;}
        void LockService::GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const {start = 0; end = 0;}
Kirill Terekhov's avatar
Kirill Terekhov committed
1147
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1148
1149


Kirill Terekhov's avatar
Kirill Terekhov committed
1150
////////class HessianMatrix
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1151

1152
1153
        void     HessianMatrix::Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg, INMOST_DATA_ENUM_TYPE mend)
        {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1154
1155
            char str[16384];
            std::ifstream input(file.c_str());
1156
            if( input.fail() ) throw -1;