1 // Copyright (c) 2018 ERGO-Code. See license.txt for license.
2
3 #include "normal_matrix.h"
4 #include <cassert>
5 #include "timer.h"
6 #include "utils.h"
7
8 namespace ipx {
9
10 // Method for computing matrix-vector products:
11 // 1 one pass over A
12 // 2 two passes over A, first by cols, second by rows
13 // 3 two passes over A, first by cols, second by cols
14 // Method 2 and 3 are implemented only for W != NULL.
15 // For W == NULL method 1 is used always.
16 //
17 // The three variants of matrix-vector products were implemented and profiled
18 // during the development phase of IPX. It turned out that the one-pass variant
19 // is the fastest on average (about 20% better than the best two-pass variant),
20 // and also the fastest on most LP models. Therefore, it is used for
21 // matrix-vector products of the form AA' here and in SplittedNormalMatrix.
22 #define MATVECMETHOD 1
23
NormalMatrix(const Model & model)24 NormalMatrix::NormalMatrix(const Model& model) : model_(model) {
25 #if MATVECMETHOD > 1
26 // The two-pass variants require n+m workspace to store the intermediate
27 // result W*AI'*rhs.
28 work_.resize(model.rows() + model.cols());
29 #endif
30 }
31
Prepare(const double * W)32 void NormalMatrix::Prepare(const double* W) {
33 W_ = W;
34 prepared_ = true;
35 }
36
time() const37 double NormalMatrix::time() const {
38 return time_;
39 }
40
reset_time()41 void NormalMatrix::reset_time() {
42 time_ = 0.0;
43 }
44
_Apply(const Vector & rhs,Vector & lhs,double * rhs_dot_lhs)45 void NormalMatrix::_Apply(const Vector& rhs, Vector& lhs,
46 double* rhs_dot_lhs) {
47 const Int m = model_.rows();
48 const Int n = model_.cols();
49 const Int* Ap = model_.AI().colptr();
50 const Int* Ai = model_.AI().rowidx();
51 const double* Ax = model_.AI().values();
52 #if MATVECMETHOD == 2
53 const Int* Atp = model_.AIt().colptr();
54 const Int* Ati = model_.AIt().rowidx();
55 const double* Atx = model_.AIt().values();
56 #endif
57 Timer timer;
58
59 assert(prepared_);
60 assert(lhs.size() == m);
61 assert(rhs.size() == m);
62
63 if (W_) {
64 #if MATVECMETHOD == 1
65 for (Int i = 0; i < m; i++)
66 lhs[i] = rhs[i] * W_[n+i];
67 for (Int j = 0; j < n; j++) {
68 Int begin = Ap[j], end = Ap[j+1];
69 double d = 0.0;
70 for (Int p = begin; p < end; p++)
71 d += rhs[Ai[p]] * Ax[p];
72 d *= W_[j];
73 for (Int p = begin; p < end; p++)
74 lhs[Ai[p]] += d * Ax[p];
75 }
76 #elif MATVECMETHOD == 2
77 for (Int j = 0; j < n; j++) {
78 Int begin = Ap[j], end = Ap[j+1];
79 double d = 0.0;
80 for (Int p = begin; p < end; p++)
81 d += rhs[Ai[p]] * Ax[p];
82 work_[j] = d * W_[j];
83 }
84 for (Int i = 0; i < m; i++)
85 lhs[i] = rhs[i] * W_[n+i];
86 for (Int i = 0; i < m; i++) {
87 Int begin = Atp[i], end = Atp[i+1]-1; // skip identity entry
88 double d = 0.0;
89 for (Int p = begin; p < end; p++)
90 d += work_[Ati[p]] * Atx[p];
91 lhs[i] += d;
92 }
93 #elif MATVECMETHOD == 3
94 for (Int j = 0; j < n; j++) {
95 Int begin = Ap[j], end = Ap[j+1];
96 double d = 0.0;
97 for (Int p = begin; p < end; p++)
98 d += rhs[Ai[p]] * Ax[p];
99 work_[j] = d * W_[j];
100 }
101 for (Int i = 0; i < m; i++)
102 lhs[i] = rhs[i] * W_[n+i];
103 for (Int j = 0; j < n; j++) {
104 Int begin = Ap[j], end = Ap[j+1];
105 double d = work_[j];
106 for (Int p = begin; p < end; p++)
107 lhs[Ai[p]] += d * Ax[p];
108 }
109 #else
110 #error "invalid MATHVECMETHOD"
111 #endif
112 } else {
113 lhs = 0.0;
114 for (Int j = 0; j < n; j++) {
115 Int begin = Ap[j], end = Ap[j+1];
116 double d = 0.0;
117 for (Int p = begin; p < end; p++)
118 d += rhs[Ai[p]] * Ax[p];
119 for (Int p = begin; p < end; p++)
120 lhs[Ai[p]] += d * Ax[p];
121 }
122 }
123 if (rhs_dot_lhs)
124 *rhs_dot_lhs = Dot(rhs,lhs);
125 time_ += timer.Elapsed();
126 }
127
128 } // namespace ipx
129