1+ #include " HybridPackedFormatHandler.h"
2+
3+ #include " CallAndTimeBlas.h"
4+ #include " DataCollector.h"
5+ #include " DenseFact.h"
6+ #include " ipm/hipo/auxiliary/Auxiliary.h"
7+
8+ namespace hipo {
9+
10+ HybridPackedFormatHandler::HybridPackedFormatHandler (
11+ const Symbolic& S, int sn, const Regul& regul, DataCollector& data,
12+ std::vector<double >& frontal)
13+ : FormatHandler(S, sn, regul, frontal), data_{data} {
14+ // initialize frontal and clique
15+ initFrontal ();
16+ initClique ();
17+ }
18+
19+ void HybridPackedFormatHandler::initFrontal () {
20+ const int n_blocks = (sn_size_ - 1 ) / nb_ + 1 ;
21+ diag_start_.resize (n_blocks);
22+ int frontal_size = getDiagStart (ldf_, sn_size_, nb_, n_blocks, diag_start_);
23+ frontal_.assign (frontal_size + extra_space, 0.0 );
24+ // NB: the plus 10 is not needed, but it avoids weird problems later on.
25+ }
26+
27+ void HybridPackedFormatHandler::initClique () {
28+ clique_.resize (S_->cliqueSize (sn_));
29+ }
30+
31+ void HybridPackedFormatHandler::assembleFrontal (int i, int j, double val) {
32+ int block = j / nb_;
33+ int ldb = ldf_ - block * nb_;
34+ int ii = i - block * nb_;
35+ int jj = j - block * nb_;
36+ frontal_[diag_start_[block] + ii + ldb * jj] = val;
37+ }
38+
39+ void HybridPackedFormatHandler::assembleFrontalMultiple (
40+ int num, const std::vector<double >& child, int nc, int child_sn, int row,
41+ int col, int i, int j) {
42+ const int jblock = col / nb_;
43+ row -= jblock * nb_;
44+ col -= jblock * nb_;
45+ const int start_block = S_->cliqueBlockStart (child_sn, jblock);
46+ const int ld = nc - nb_ * jblock;
47+
48+ int block = j / nb_;
49+ int ldb = ldf_ - block * nb_;
50+ int ii = i - block * nb_;
51+ int jj = j - block * nb_;
52+
53+ callAndTime_daxpy (num, 1.0 , &child[start_block + row + ld * col], 1 ,
54+ &frontal_[diag_start_[block] + ii + ldb * jj], 1 , data_);
55+ }
56+
57+ int HybridPackedFormatHandler::denseFactorise (double reg_thresh) {
58+ int status;
59+
60+ status = denseFactFP2FH (frontal_.data (), ldf_, sn_size_, nb_, data_);
61+ if (status) return status;
62+
63+ // find the position within pivot_sign corresponding to this supernode
64+ int sn_start = S_->snStart (sn_);
65+ const int * pivot_sign = &S_->pivotSign ().data ()[sn_start];
66+
67+ status =
68+ denseFactFH (' P' data (), clique_.data (),
69+ pivot_sign, reg_thresh, regul_, local_reg_.data (),
70+ swaps_.data (), pivot_2x2_.data (), S_->parNode (), data_);
71+
72+ return status;
73+ }
74+
75+ void HybridPackedFormatHandler::assembleClique (const std::vector<double >& child,
76+ int nc, int child_sn) {
77+ // go through the columns of the contribution of the child
78+ for (int col = 0 ; col < nc; ++col) {
79+ // relative index of column in the frontal matrix
80+ int j = S_->relindClique (child_sn, col);
81+
82+ if (j >= sn_size_) {
83+ // assemble into clique
84+
85+ // adjust relative index to access clique
86+ j -= sn_size_;
87+
88+ // go through the rows of the contribution of the child
89+ int row = col;
90+ while (row < nc) {
91+ // relative index of the entry in the matrix clique
92+ const int i = S_->relindClique (child_sn, row) - sn_size_;
93+
94+ // how many entries to sum
95+ const int consecutive = S_->consecutiveSums (child_sn, row);
96+
97+ // use daxpy_ for summing consecutive entries
98+
99+ const int jblock_c = col / nb_;
100+ const int jb_c = std::min (nb_, nc - nb_ * jblock_c);
101+ const int row_c = row - jblock_c * nb_;
102+ const int col_c = col - jblock_c * nb_;
103+ const int start_block_c = S_->cliqueBlockStart (child_sn, jblock_c);
104+ const int ld_c = nc - nb_ * jblock_c;
105+
106+ const int jblock = j / nb_;
107+ const int jb = std::min (nb_, ldc_ - nb_ * jblock);
108+ const int ii = i - jblock * nb_;
109+ const int jj = j - jblock * nb_;
110+ const int start_block = S_->cliqueBlockStart (sn_, jblock);
111+ const int ld = ldc_ - nb_ * jblock;
112+
113+ callAndTime_daxpy (consecutive, 1.0 ,
114+ &child[start_block_c + row_c + ld_c * col_c], 1 ,
115+ &clique_[start_block + ii + ld * jj], 1 , data_);
116+
117+ row += consecutive;
118+ }
119+ }
120+ }
121+ }
122+
123+ void HybridPackedFormatHandler::extremeEntries () {
124+ #ifdef HIPO_COLLECT_EXPENSIVE_DATA
125+ double minD = std::numeric_limits<double >::max ();
126+ double maxD = 0.0 ;
127+ double minoffD = std::numeric_limits<double >::max ();
128+ double maxoffD = 0.0 ;
129+
130+ // number of blocks of columns
131+ const int n_blocks = (sn_size_ - 1 ) / nb_ + 1 ;
132+
133+ // index to access frontal
134+ int index{};
135+
136+ // go through blocks of columns for this supernode
137+ for (int j = 0 ; j < n_blocks; ++j) {
138+ // number of columns in the block
139+ const int jb = std::min (nb_, sn_size_ - nb_ * j);
140+
141+ for (int k = 0 ; k < jb; ++k) {
142+ // off diagonal entries
143+ for (int i = 0 ; i < k; ++i) {
144+ if (frontal_[index] != 0.0 ) {
145+ minoffD = std::min (minoffD, std::abs (frontal_[index]));
146+ maxoffD = std::max (maxoffD, std::abs (frontal_[index]));
147+ }
148+ index++;
149+ }
150+
151+ // diagonal entry
152+ minD = std::min (minD, std::abs (1.0 / frontal_[index]));
153+ maxD = std::max (maxD, std::abs (1.0 / frontal_[index]));
154+
155+ index += jb - k;
156+ }
157+
158+ const int entries_left = (ldf_ - nb_ * j - jb) * jb;
159+
160+ for (int i = 0 ; i < entries_left; ++i) {
161+ if (frontal_[index] != 0.0 ) {
162+ minoffD = std::min (minoffD, std::abs (frontal_[index]));
163+ maxoffD = std::max (maxoffD, std::abs (frontal_[index]));
164+ }
165+ index++;
166+ }
167+ }
168+
169+ data_.setExtremeEntries (minD, maxD, minoffD, maxoffD);
170+ #endif
171+ }
172+
173+ } // namespace hipo
0 commit comments