1 /* 2 HMat-OSS (HMatrix library, open source software) 3 4 Copyright (C) 2014-2015 Airbus Group SAS 5 6 This program is free software; you can redistribute it and/or 7 modify it under the terms of the GNU General Public License 8 as published by the Free Software Foundation; either version 2 9 of the License, or (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with this program; if not, write to the Free Software 18 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 19 20 http://github.com/jeromerobert/hmat-oss 21 */ 22 23 #include "config.h" 24 25 /*! \file 26 \ingroup HMatrix 27 \brief Templated Matrix class for recursive algorithms used by HMatrix. 28 */ 29 30 #include "recursion.hpp" 31 #include "h_matrix.hpp" 32 #include "common/context.hpp" 33 34 namespace hmat { 35 36 template<typename T, typename Mat> recursiveLdltDecomposition(hmat_progress_t * progress)37 void RecursionMatrix<T, Mat>::recursiveLdltDecomposition(hmat_progress_t * progress) { 38 39 // Recursive LDLT factorization: 40 // 41 // [ h11 |th21 ] [ L11 | 0 ] [ D1 | 0 ] [ tL11 | tL21 ] 42 // [ ----+---- ] = [ ----+---- ] * [----+----] * [------+------] 43 // [ h21 | h22 ] [ L21 | L22 ] [ 0 | D2 ] [ 0 | tL22 ] 44 // 45 // h11 = L11 * D1 * tL11 (gives L11 and D1 by LDLT decomposition) 46 // h21 = L21 * D1 * tL11 (gives L21 when L11 and D1 are known) 47 // h22 = L21 * D1 * tL21 + L22 * D2 * tL22 (gives L22 and D2 by LDLT decomposition of h22-L21*D1*tL21) 48 49 // for all i, j<=i : hij = sum_k Lik Dk t{Ljk} k<=i,j 50 // The algorithm loops over 3 steps: for all k=1, 2, ..., n 51 // - We factorize the element (k,k) 52 // - We use "solve" to compute the rest of the column 'k' 53 // - We update the rest of the matrix [k+1, .., n]x[k+1, .., n] (below diag) 54 55 HMAT_ASSERT_MSG(me()->nrChildRow()==me()->nrChildCol(), 56 "RecursionMatrix<T, Mat>::recursiveLdltDecomposition: case not allowed " 57 "Nr Child A[%d, %d] Dimensions A=%s ", 58 me()->nrChildRow(), me()->nrChildCol(), me()->description().c_str()); 59 60 for (int k=0 ; k<me()->nrChildRow() ; k++) { 61 // Hkk <- Lkk * Dk * tLkk 62 me()->get(k,k)->ldltDecomposition(progress); 63 // Solve the rest of column k: solve Lik Dk tLkk = Hik and get Lik 64 for (int i=k+1 ; i<me()->nrChildRow() ; i++) { 65 if (!me()->get(i,k)) 66 continue; 67 me()->get(k,k)->solveUpperTriangularRight(me()->get(i,k), Factorization::LDLT, Diag::NONUNIT, Uplo::LOWER); 68 me()->get(i,k)->multiplyWithDiag(me()->get(k,k), Side::RIGHT, true); 69 } 70 // update the rest of the matrix [k+1, .., n]x[k+1, .., n] (below diag) 71 for (int i=k+1 ; i<me()->nrChildRow() ; i++) { 72 if (!me()->get(i,k)) 73 continue; 74 // Hij <- Hij - Lik Dk tLjk 75 // if i=j, we can use mdmtProduct, otherwise we use mdntProduct 76 for (int j=k+1 ; j<i ; j++) 77 if (me()->get(i,j) && me()->get(j,k)) 78 me()->get(i,j)->mdntProduct(me()->get(i,k), me()->get(k,k), me()->get(j,k)); // hij -= Lik.Dk.tLjk 79 me()->get(i,i)->mdmtProduct(me()->get(i,k), me()->get(k,k)); // hii -= Lik.Dk.tLik 80 } 81 } 82 83 } 84 85 template<typename T, typename Mat> recursiveSolveUpperTriangularRight(Mat * b,Factorization algo,Diag diag,Uplo uplo) const86 void RecursionMatrix<T, Mat>::recursiveSolveUpperTriangularRight(Mat* b, Factorization algo, Diag diag, Uplo uplo) const { 87 88 // Forward substitution: 89 // [ X11 | X12 ] [ U11 | U12 ] [ b11 | b12 ] 90 // [ ----+---- ] * [-----+-----] = [ ----+---- ] 91 // [ X21 | X22 ] [ 0 | U22 ] [ b21 | b22 ] 92 // 93 // X11 * U11 = b11 (by recursive forward substitution) 94 // X21 * U11 = b21 (by recursive forward substitution) 95 // X11 * U12 + X12 * U22 = b12 (forward substitution of X12*U22=b12-X11*U12) 96 // X21 * U12 + X22 * U22 = b22 (forward substitution of X22*U22=b22-X21*U12) 97 // X and b are not necessarily square 98 99 // First we handle the general case, where dimensions (in terms of number of children) are compatible 100 if (me()->nrChildRow() == b->nrChildCol()) { 101 102 for (int k=0 ; k<b->nrChildRow() ; k++) // loop on the lines of b 103 for (int i=0 ; i<me()->nrChildRow() ; i++) { 104 // Update b[k,i] with the contribution of the solutions already computed b[k,j] j<i 105 if (!b->get(k, i)) continue; 106 for (int j=0 ; j<i ; j++) { 107 const Mat* u_ji = (uplo == Uplo::LOWER ? me()->get(i, j) : me()->get(j, i)); 108 if (b->get(k, j) && u_ji) 109 b->get(k, i)->gemm('N', uplo == Uplo::LOWER ? 'T' : 'N', Constants<T>::mone, b->get(k, j), u_ji, Constants<T>::pone); 110 } 111 // Solve the i-th diagonal system 112 me()->get(i, i)->solveUpperTriangularRight(b->get(k,i), algo, diag, uplo); 113 } 114 115 } else if (me()->nrChildRow()>1 && b->nrChildCol()==1 && b->nrChildRow()>1) { 116 // Then we handle the specific case where me() is divided in rows, and b is not divided in cols: we recurse on b's children only 117 // [ X11 ] [ U11 | U12 ] [ b11 ] 118 // [ --------- ] * [ ----+---- ] = [ --------- ] 119 // [ X21 ] [ 0 | U22 ] [ b21 ] 120 for (int k=0 ; k<b->nrChildRow() ; k++) // loop on the rows of b 121 me()->recursiveSolveUpperTriangularRight(b->get(k,0), algo, diag, uplo); 122 123 } else { 124 HMAT_ASSERT_MSG(false, "RecursionMatrix<T, Mat>::recursiveSolveUpperTriangularRight: case not yet handled " 125 "Nr Child A[%d, %d] b[%d, %d] " 126 "Dimensions A=%s b=%s", 127 me()->nrChildRow(), me()->nrChildCol(), b->nrChildRow(), b->nrChildCol(), 128 me()->description().c_str(), b->description().c_str()); 129 130 } 131 132 } 133 134 template<typename T, typename Mat> recursiveMdmtProduct(const Mat * m,const Mat * d)135 void RecursionMatrix<T, Mat>::recursiveMdmtProduct(const Mat* m, const Mat* d) { 136 137 // 138 // [ h11 |th21 ] [ M11 | M12 ] [ D1 | 0 ] [ tM11 | tM21 ] 139 // [ ----+---- ] -= [ ----+---- ] * [----+----] * [------+------] 140 // [ h21 | h22 ] [ M21 | M22 ] [ 0 | D2 ] [ tM12 | tM22 ] 141 // 142 // h11 -= M11.D1.tM11 + M12.D2.tM12 143 // h21 -= M21.D1.tM11 + M22.D2.tM12 144 // h22 -= M21.D1.tM21 + M22.D2.tM22 145 // 146 // hij -= sum_k Mik.Dk.tMjk 147 148 // First we handle the general case, where dimensions (in terms of number of children) are compatible. 149 // recursiveMdmtProduct is called when neither this and m are leaves, but d may be a leaf. 150 const int block_row_d = d->isLeaf() ? 1 : d->nrChildRow(); 151 const int block_col_d = d->isLeaf() ? 1 : d->nrChildCol(); 152 if (me()->nrChildRow()==me()->nrChildCol() && block_row_d==block_col_d && me()->nrChildRow()==m->nrChildRow() && m->nrChildCol()==block_row_d) { 153 if (d->isLeaf()) { 154 // hij -= Mik.Dk.tMjk : if i=j, we use mdmtProduct. Otherwise, mdntProduct 155 for(int i=0 ; i<me()->nrChildRow() ; i++) { 156 if (!m->get(i,0)) 157 continue; 158 for (int j=0 ; j<i ; j++) 159 if (me()->get(i,j) && m->get(j,0)) 160 me()->get(i,j)->mdntProduct(m->get(i,0), d, m->get(j,0)); // hij -= Mi0.D.tMj0 161 me()->get(i,i)->mdmtProduct(m->get(i,0), d); // hii -= Mi0.D.tMi0 162 } 163 } else { 164 // hij -= Mik.Dk.tMjk : if i=j, we use mdmtProduct. Otherwise, mdntProduct 165 for(int i=0 ; i<me()->nrChildRow() ; i++) 166 for(int k=0 ; k<m->nrChildCol() ; k++) { 167 const Mat* m_ik = m->get(i,k); 168 if (!m_ik) 169 continue; 170 const Mat* d_k = d->get(k,k); 171 for (int j=0 ; j<i ; j++) 172 if (me()->get(i,j) && m->get(j,k)) 173 me()->get(i,j)->mdntProduct(m_ik, d_k, m->get(j,k)); // hij -= Mik.Dk.tMjk 174 me()->get(i,i)->mdmtProduct(m_ik, d_k); // hii -= Mik.Dk.tMik 175 } 176 } 177 } else { 178 HMAT_ASSERT_MSG(false, "RecursionMatrix<T, Mat>::recursiveMdmtProduct: case not yet handled " 179 "Nr Child this[%d, %d] m[%d, %d] d[%d, %d]" 180 "Dimensions this=%s m=%s d=%s", 181 me()->nrChildRow(), me()->nrChildCol(), m->nrChildRow(), m->nrChildCol(), 182 d->nrChildRow(), d->nrChildCol(), 183 me()->description().c_str(), m->description().c_str(), d->description().c_str()); 184 185 } 186 } 187 188 template<typename T, typename Mat> recursiveSolveLowerTriangularLeft(Mat * b,Factorization algo,Diag diag,Uplo uplo,MainOp mainOp) const189 void RecursionMatrix<T, Mat>::recursiveSolveLowerTriangularLeft(Mat* b, Factorization algo, Diag diag, Uplo uplo, MainOp mainOp) const { 190 191 // Forward substitution: 192 // [ L11 | 0 ] [ X11 | X12 ] [ b11 | b12 ] 193 // [ ----+---- ] * [-----+-----] = [ ----+---- ] 194 // [ L21 | L22 ] [ X21 | X22 ] [ b21 | b22 ] 195 // 196 // L11 * X11 = b11 (by recursive forward substitution) 197 // L11 * X12 = b12 (by recursive forward substitution) 198 // L21 * X11 + L22 * X21 = b21 (forward substitution of L22*X21=b21-L21*X11) 199 // L21 * X12 + L22 * X22 = b22 (forward substitution of L22*X22=b22-L21*X12) 200 // X and b are not necessarily square 201 202 // First we handle the general case, where dimensions (in terms of number of children) are compatible 203 if (me()->nrChildCol() == b->nrChildRow()) { 204 205 for (int k=0 ; k<b->nrChildCol() ; k++) // loop on the column of b 206 for (int i=0 ; i<me()->nrChildRow() ; i++) { 207 // Update b[i,k] with the contribution of the solutions already computed b[j,k] j<i 208 if (!b->get(i, k)) continue; 209 for (int j=0 ; j<i ; j++) 210 if (me()->get(i,j) && b->get(j,k)) 211 b->get(i, k)->gemm('N', 'N', Constants<T>::mone, me()->get(i, j), b->get(j,k), 212 Constants<T>::pone, mainOp); 213 // Solve the i-th diagonal system 214 me()->get(i, i)->solveLowerTriangularLeft(b->get(i,k), algo, diag, uplo, mainOp); 215 } 216 217 } else if (me()->nrChildCol()>1 && b->nrChildRow()==1 && b->nrChildCol()>1) { 218 // Then we handle the specific case where me() is divided in columns, and b is not divided in rows: we recurse on b's children only 219 // [ L11 | 0 ] [ | ] [ | ] 220 // [ ----+---- ] * [ X11 | X12 ] = [ b11 | b12 ] 221 // [ L21 | U22 ] [ | ] [ | ] 222 for (int k=0 ; k<b->nrChildCol() ; k++) // loop on the column of b 223 me()->recursiveSolveLowerTriangularLeft(b->get(0,k), algo, diag, uplo, mainOp); 224 } else { 225 HMAT_ASSERT_MSG(false, "RecursionMatrix<T, Mat>::recursiveSolveLowerTriangularLeft: case not yet handled " 226 "Nr Child A[%d, %d] b[%d, %d] " 227 "Dimensions A=%s b=%s", 228 me()->nrChildRow(), me()->nrChildCol(), b->nrChildRow(), b->nrChildCol(), 229 me()->description().c_str(), b->description().c_str()); 230 231 } 232 } 233 234 template<typename T, typename Mat> recursiveLuDecomposition(hmat_progress_t * progress)235 void RecursionMatrix<T, Mat>::recursiveLuDecomposition(hmat_progress_t * progress) { 236 237 // | | | | | | | | | 238 // | h11 | h12 | | L11 | | | U11 | U12 | 239 // |-----|-----| = |-----|-----| * |-----|-----| 240 // | h21 | h22 | | L21 | L22 | | | U22 | 241 // | | | | | | | | | 242 // 243 // h11 = L11 * U11 => (L11,U11) = h11.luDecomposition 244 // h12 = L11 * U12 => trsm L 245 // h21 = L21 * U11 => trsm R 246 // h22 = L21 * U12 + L22 * U22 => (L22,U22) = (h22 - L21*U12).luDecomposition() 247 // 248 // hij = sum_k Lik Ukj k<=i,j 249 // The algorithm loops over 3 steps: for all k=1, 2, ..., n 250 // - We factorize the element (k,k) 251 // - We use "solve" to compute the rest of the columns and rows 'k' 252 // - We update the rest of the matrix [k+1, .., n]x[k+1, .., n] 253 254 HMAT_ASSERT_MSG(me()->nrChildRow()==me()->nrChildCol(), 255 "RecursionMatrix<T, Mat>::recursiveLuDecomposition: case not allowed " 256 "Nr Child A[%d, %d] Dimensions A=%s ", 257 me()->nrChildRow(), me()->nrChildCol(), me()->description().c_str()); 258 259 for (int k=0 ; k<me()->nrChildRow() ; k++) { 260 if(me()->get(k,k) == nullptr) 261 // inert diagonal block. The associated row & column are considered as also inert. 262 continue; 263 // Hkk <- Lkk * Ukk 264 me()->get(k,k)->luDecomposition(progress); 265 // Solve the rest of line k: solve Lkk Uki = Hki and get Uki 266 for (int i=k+1 ; i<me()->nrChildRow() ; i++) 267 if (me()->get(k,i)) 268 me()->get(k,k)->solveLowerTriangularLeft(me()->get(k,i), Factorization::LU, Diag::UNIT, Uplo::LOWER); 269 // Solve the rest of column k: solve Lik Ukk = Hik and get Lik 270 for (int i=k+1 ; i<me()->nrChildRow() ; i++) 271 if (me()->get(i,k)) 272 me()->get(k,k)->solveUpperTriangularRight(me()->get(i,k), Factorization::LU, Diag::NONUNIT, Uplo::UPPER); 273 // update the rest of the matrix starting at (k+1, k+1) 274 for (int i=k+1 ; i<me()->nrChildRow() ; i++) { 275 if (!me()->get(i,k)) 276 continue; 277 for (int j=k+1 ; j<me()->nrChildRow() ; j++) 278 // Hij <- Hij - Lik Ukj 279 if (me()->get(i,j) && me()->get(k,j)) 280 me()->get(i,j)->gemm('N', 'N', Constants<T>::mone, me()->get(i,k), me()->get(k,j), Constants<T>::pone); 281 } 282 } 283 284 } 285 286 template<typename T, typename Mat> recursiveInverseNosym()287 void RecursionMatrix<T, Mat>::recursiveInverseNosym() { 288 289 // Matrix inversion: 290 // The idea to inverse M is to consider the extended matrix obtained by putting Identity next to M : 291 // 292 // [ M11 | M12 | I | 0 ] 293 // [ ----+-----+-----+---- ] 294 // [ M21 | M22 | 0 | I ] 295 // 296 // We then apply operations on the line of this matrix (matrix multiplication of an entire line, 297 // linear combination of lines) 298 // to transform the 'M' part into identity. Doing so, the identity part will at the end contain M-1. 299 // We loop on the column of M. 300 // At the end of loop 'k', the 'k' first columns of 'M' are now identity, 301 // and the 'k' first columns of Identity have changed (it's no longer identity, it's not yet M-1). 302 // The matrix 'this' stores the first 'k' block of the identity part of the extended matrix, and the last n-k blocks of the M part 303 // At the end, 'this' contains M-1 304 305 HMAT_ASSERT_MSG(me()->nrChildRow()==me()->nrChildCol(), 306 "RecursionMatrix<T, Mat>::recursiveInverseNosym: case not allowed " 307 "Nr Child A[%d, %d] Dimensions A=%s ", 308 me()->nrChildRow(), me()->nrChildCol(), me()->description().c_str()); 309 310 for (int k=0 ; k<me()->nrChildRow() ; k++){ 311 // Inverse M_kk 312 me()->get(k,k)->inverse(); 313 // Update line 'k' = left-multiplied by M_kk-1 314 for (int j=0 ; j<me()->nrChildCol() ; j++) 315 if (j!=k) { 316 // Mkj <- Mkk^-1 Mkj we use a temp matrix X because this type of product is not allowed with gemm (beta=0 erases Mkj before using it !) 317 Mat* x = me()->get(k,j)->copy(); 318 me()->get(k,j)->gemm('N', 'N', Constants<T>::pone, me()->get(k,k), x, Constants<T>::zero); 319 x->destroy(); 320 } 321 // Update the rest of matrix M 322 for (int i=0 ; i<me()->nrChildRow() ; i++) 323 // line 'i' -= Mik x line 'k' (which has just been multiplied by Mkk-1) 324 for (int j=0 ; j<me()->nrChildCol() ; j++) 325 if (i!=k && j!=k) 326 // Mij <- Mij - Mik Mkk^-1 Mkj (with Mkk-1.Mkj allready stored in Mkj) 327 me()->get(i,j)->gemm('N', 'N', Constants<T>::mone, me()->get(i,k), me()->get(k,j), Constants<T>::pone); 328 // Update column 'k' = right-multiplied by -M_kk-1 329 for (int i=0 ; i<me()->nrChildRow() ; i++) 330 if (i!=k) { 331 // Mik <- - Mik Mkk^-1 332 Mat* x = me()->get(i,k)->copy(); 333 me()->get(i,k)->gemm('N', 'N', Constants<T>::mone, x, me()->get(k,k), Constants<T>::zero); 334 x->destroy(); 335 } 336 } 337 338 } 339 340 template<typename T, typename Mat> recursiveLltDecomposition(hmat_progress_t * progress)341 void RecursionMatrix<T, Mat>::recursiveLltDecomposition(hmat_progress_t * progress) { 342 343 // | | | | | | | | | 344 // | h11 | h21 | | L1 | | | L1t | Lt | 345 // |-----|-----| = |-----|-----| * |-----|-----| 346 // | h21 | h22 | | L | L2 | | | L2t | 347 // | | | | | | | | | 348 // 349 // h11 = L1 * L1t => L1 = h11.lltDecomposition 350 // h21 = L*L1t => L = L1t.solve(h21) => trsm R with L1 lower stored 351 // h22 = L*Lt + L2 * L2t => L2 = (h22 - L*Lt).lltDecomposition() 352 // 353 // 354 // for all i, j<=i : hij = sum_k Lik t{Ljk} k<=i,j 355 // The algorithm loops over 3 steps: for all k=1, 2, ..., n 356 // - We factorize the element (k,k) 357 // - We use "solve" to compute the rest of the column 'k' 358 // - We update the rest of the matrix [k+1, .., n]x[k+1, .., n] (below diag) 359 360 HMAT_ASSERT_MSG(me()->nrChildRow()==me()->nrChildCol(), 361 "RecursionMatrix<T, Mat>::recursiveLltDecomposition: case not allowed " 362 "Nr Child A[%d, %d] Dimensions A=%s ", 363 me()->nrChildRow(), me()->nrChildCol(), me()->description().c_str()); 364 365 for (int k=0 ; k<me()->nrChildRow() ; k++) { 366 // Hkk <- Lkk * tLkk 367 me()->get(k,k)->lltDecomposition(progress); 368 // Solve the rest of column k: solve Lik tLkk = Hik and get Lik 369 for (int i=k+1 ; i<me()->nrChildRow() ; i++) 370 if (me()->get(i,k)) 371 me()->get(k,k)->solveUpperTriangularRight(me()->get(i,k), Factorization::LLT, Diag::NONUNIT, Uplo::LOWER); 372 // update the rest of the matrix [k+1, .., n]x[k+1, .., n] (below diag) 373 for (int i=k+1 ; i<me()->nrChildRow() ; i++) { 374 if (!me()->get(i,k)) 375 continue; 376 for (int j=k+1 ; j<=i ; j++) 377 // Hij <- Hij - Lik tLjk 378 if (me()->get(i,j) && me()->get(j,k)) 379 me()->get(i,j)->gemm('N', 'T', Constants<T>::mone, me()->get(i,k), me()->get(j,k), Constants<T>::pone); 380 } 381 } 382 } 383 384 template<typename T, typename Mat> recursiveSolveUpperTriangularLeft(Mat * b,Factorization algo,Diag diag,Uplo uplo,MainOp mainOp) const385 void RecursionMatrix<T, Mat>::recursiveSolveUpperTriangularLeft(Mat* b, 386 Factorization algo, Diag diag, Uplo uplo, MainOp mainOp) const { 387 388 // Backward substitution: 389 // [ U11 | U12 ] [ X11 | X12 ] [ b11 | b12 ] 390 // [ ----+---- ] * [-----+-----] = [ ----+---- ] 391 // [ 0 | U22 ] [ X21 | X22 ] [ b21 | b22 ] 392 // 393 // U22 * X21 = b21 (by recursive backward substitution) 394 // U22 * X22 = b22 (by recursive backward substitution) 395 // U11 * X12 + U12 * X22 = b12 (backward substitution of U11*X12=b12-U12*X22) 396 // U11 * X11 + U12 * X21 = b11 (backward substitution of U11*X11=b11-U12*X21) 397 // X and b are not necessarily square 398 399 // First we handle the general case, where dimensions (in terms of number of children) are compatible 400 if (me()->nrChildCol() == b->nrChildRow()) { 401 402 for (int k=0 ; k<b->nrChildCol() ; k++) { // Loop on the column of the RHS 403 for (int i=me()->nrChildRow()-1 ; i>=0 ; i--) { 404 if (!b->get(i,k)) 405 continue; 406 // Solve the i-th diagonal system 407 me()->get(i, i)->solveUpperTriangularLeft(b->get(i,k), algo, diag, uplo, mainOp); 408 // Update b[j,k] j<i with the contribution of the solutions just computed b[i,k] 409 for (int j=0 ; j<i ; j++) { 410 const Mat* u_ji = (uplo == Uplo::LOWER ? me()->get(i, j) : me()->get(j, i)); 411 if(u_ji && b->get(j,k)) 412 b->get(j,k)->gemm(uplo == Uplo::LOWER ? 'T' : 'N', 'N', Constants<T>::mone, u_ji, 413 b->get(i,k), Constants<T>::pone, mainOp); 414 } 415 } 416 } 417 418 } else if (me()->nrChildCol()>1 && b->nrChildRow()==1 && b->nrChildCol()>1) { 419 // Then we handle the specific case where me() is divided in columns, and b is not divided in rows: we recurse on b's children only 420 // [ U11 | U12 ] [ | ] [ | ] 421 // [ ----+---- ] * [ X11 | X12 ] = [ b11 | b12 ] 422 // [ 0 | U22 ] [ | ] [ | ] 423 for (int k=0 ; k<b->nrChildCol() ; k++) // loop on the column of b 424 me()->recursiveSolveUpperTriangularLeft(b->get(0,k), algo, diag, uplo); 425 426 } else { 427 HMAT_ASSERT_MSG(false, "RecursionMatrix<T, Mat>::recursiveSolveUpperTriangularLeft: case not yet handled " 428 "Nr Child A[%d, %d] b[%d, %d] " 429 "Dimensions A=%s b=%s", 430 me()->nrChildRow(), me()->nrChildCol(), b->nrChildRow(), b->nrChildCol(), 431 me()->description().c_str(), b->description().c_str()); 432 433 } 434 435 } 436 transposeMeta(bool temporaryOnly)437 template <typename T, typename Mat> void RecursionMatrix<T, Mat>::transposeMeta(bool temporaryOnly) { 438 if (!me()->isLeaf()) { 439 // We cannot not, in general, transpose in-place, so we need a backup of 'children' 440 std::vector<Mat*> children_bak(me()->nrChild()); 441 for(int i = 0; i < me()->nrChild(); i++) 442 children_bak[i] = me()->getChild(i); 443 // and finally we fill 'children' 444 int k = 0; 445 for (int j = 0; j < me()->nrChildRow(); j++) 446 for (int i = 0; i < me()->nrChildCol(); i++) 447 me()->getChild(j + i * me()->nrChildRow()) = children_bak[k++]; 448 for (int i = 0; i < me()->nrChild(); i++) 449 if (me()->getChild(i)) 450 me()->getChild(i)->transposeMeta(temporaryOnly); 451 } 452 } 453 } // end namespace hmat 454