@@ -44,6 +44,7 @@ static void staticReg(double& pivot, Int sign, const FHoptions& options,
4444 // apply static regularisation
4545
4646 double old_pivot = pivot;
47+ if (sign == 0 ) sign = pivot > 0 ? 1 : -1 ;
4748 if (sign > 0 )
4849 pivot += options.reg_p ;
4950 else
@@ -85,29 +86,24 @@ static bool blockBunchKaufman(const Int j, Int n, double* A, Int lda,
8586 auto res = maxInCol (j, n, j, A, lda);
8687 double gamma_j = res.second ;
8788 const Int r = res.first ;
88- double Ajj = sign[j] > 0 ? A[j + lda * j] + options.reg_d
89- : A[j + lda * j] - options.reg_p ;
89+ const Int sign_j = sign[j] != 0 ? sign[j] : (A[j + lda * j] > 0 ? 1 : -1 );
90+ double Ajj = sign_j > 0 ? A[j + lda * j] + options.reg_d
91+ : A[j + lda * j] - options.reg_p ;
9092
91- if (std::max (std::abs (Ajj), gamma_j) <= thresh || sign[j] * Ajj < 0 ||
93+ if (std::max (std::abs (Ajj), gamma_j) <= thresh || sign_j * Ajj < 0 ||
9294 j == n - 1 ) {
9395 // Must accept current pivot
9496 double old_pivot = A[j + lda * j];
95- staticReg (A[j + lda * j], sign[j] , options, totalreg[j]);
97+ staticReg (A[j + lda * j], sign_j , options, totalreg[j]);
9698
97- if (sign[j] * A[j + lda * j] < 0 ) {
99+ if (sign_j * A[j + lda * j] < 0 ) {
98100 data.setWrongSign (A[j + lda * j]);
99- // A[j + lda * j] = sign[j] * thresh;
101+ // A[j + lda * j] = sign_j * thresh;
100102 }
101103
102104 if (std::max (std::abs (Ajj), gamma_j) < thresh) {
103105 // perturbe pivot
104- double temp_sign = sign[j];
105- if (temp_sign == 0 && A[j + lda * j] > 0 )
106- temp_sign = 1 ;
107- else if (temp_sign == 0 )
108- temp_sign = -1 ;
109-
110- A[j + lda * j] = temp_sign * thresh;
106+ A[j + lda * j] = sign_j * thresh;
111107 data.countRegPiv ();
112108 }
113109 totalreg[j] = A[j + lda * j] - old_pivot;
@@ -118,35 +114,38 @@ static bool blockBunchKaufman(const Int j, Int n, double* A, Int lda,
118114 assert (r >= 0 );
119115 res = maxInCol (j, n, r, A, lda);
120116 double gamma_r = res.second ;
121- double Arr = sign[r] > 0 ? A[r + lda * r] + options.reg_p
122- : A[r + lda * r] - options.reg_d ;
117+ const Int sign_r = sign[r] != 0 ? sign[r] : (A[r + lda * r] > 0 ? 1 : -1 );
118+ double Arr = sign_r > 0 ? A[r + lda * r] + options.reg_p
119+ : A[r + lda * r] - options.reg_d ;
123120
124121 if ((std::abs (Ajj) >= kAlphaBK * gamma_j ||
125122 std::abs (Ajj) * gamma_r >= kAlphaBK * gamma_j * gamma_j)) {
126123 // Accept current pivot
127- staticReg (A[j + lda * j], sign[j] , options, totalreg[j]);
124+ staticReg (A[j + lda * j], sign_j , options, totalreg[j]);
128125
129- if (sign[j] * A[j + lda * j] < 0 ) data.setWrongSign (A[j + lda * j]);
126+ if (sign_j * A[j + lda * j] < 0 ) data.setWrongSign (A[j + lda * j]);
130127
131128 } else if (std::abs (Arr) >= kAlphaBK * gamma_r) {
132129 // Use pivot r
133130
134131 swapCols (' U' , n, A, lda, j, r, swaps, sign, data);
135- staticReg (A[j + lda * j], sign[j] , options, totalreg[j]);
132+ staticReg (A[j + lda * j], sign_j , options, totalreg[j]);
136133
137- if (sign[j] * A[j + lda * j] < 0 ) data.setWrongSign (A[j + lda * j]);
134+ if (sign_j * A[j + lda * j] < 0 ) data.setWrongSign (A[j + lda * j]);
138135
139136 } else {
140137 // Use 2x2 pivot (j,r)
141138
139+ const Int sign_jp1 = sign[j + 1 ] != 0
140+ ? sign[j + 1 ]
141+ : (A[j + 1 + lda * (j + 1 )] > 0 ? 1 : -1 );
142142 swapCols (' U' , n, A, lda, j + 1 , r, swaps, sign, data);
143143 flag_2x2 = true ;
144- staticReg (A[j + lda * j], sign[j], options, totalreg[j]);
145- staticReg (A[j + 1 + lda * (j + 1 )], sign[j + 1 ], options,
146- totalreg[j + 1 ]);
144+ staticReg (A[j + lda * j], sign_j, options, totalreg[j]);
145+ staticReg (A[j + 1 + lda * (j + 1 )], sign_jp1, options, totalreg[j + 1 ]);
147146
148- if (sign[j] * A[j + lda * j] < 0 ) data.setWrongSign (A[j + lda * j]);
149- if (sign[j + 1 ] * A[j + 1 + lda * (j + 1 )] < 0 )
147+ if (sign_j * A[j + lda * j] < 0 ) data.setWrongSign (A[j + lda * j]);
148+ if (sign_jp1 * A[j + 1 + lda * (j + 1 )] < 0 )
150149 data.setWrongSign (A[j + 1 + lda * (j + 1 )]);
151150 }
152151 }
0 commit comments