1 /* 2 * 3 * Copyright (c) Kresimir Fresl 2003 4 * 5 * Distributed under the Boost Software License, Version 1.0. 6 * (See accompanying file LICENSE_1_0.txt or copy at 7 * http://www.boost.org/LICENSE_1_0.txt) 8 * 9 * Author acknowledges the support of the Faculty of Civil Engineering, 10 * University of Zagreb, Croatia. 11 * 12 */ 13 14 /* for UMFPACK Copyright, License and Availability see umfpack_inc.hpp */ 15 16 17 #ifndef BOOST_NUMERIC_BINDINGS_UMFPACK_HPP 18 #define BOOST_NUMERIC_BINDINGS_UMFPACK_HPP 19 20 21 #include <boost/noncopyable.hpp> 22 #include <boost/numeric/bindings/umfpack/umfpack_overloads.hpp> 23 #include <boost/numeric/bindings/value_type.hpp> 24 #include <boost/numeric/bindings/begin.hpp> 25 #include <boost/numeric/bindings/end.hpp> 26 #include <boost/numeric/bindings/size.hpp> 27 #include <boost/numeric/bindings/data_order.hpp> 28 #include <boost/numeric/bindings/index_base.hpp> 29 30 31 namespace boost { namespace numeric { namespace bindings { namespace umfpack { 32 33 template <typename MatrA> check_umfpack_structure()34 void check_umfpack_structure() 35 { 36 BOOST_STATIC_ASSERT((boost::is_same< 37 typename bindings::detail::property_at< MatrA, tag::matrix_type >::type, 38 tag::general 39 >::value)); 40 BOOST_STATIC_ASSERT((boost::is_same< 41 typename bindings::result_of::data_order<MatrA>::type, 42 tag::column_major 43 >::value)); 44 typedef typename bindings::result_of::index_base<MatrA>::type index_b; 45 BOOST_STATIC_ASSERT(index_b::value == 0); 46 typedef typename bindings::detail::property_at< 47 MatrA, tag::data_structure >::type storage_f; 48 BOOST_STATIC_ASSERT( 49 (boost::is_same<storage_f, tag::compressed_sparse>::value || 50 boost::is_same<storage_f, tag::coordinate_sparse>::value )); 51 } 52 53 template <typename T = double> 54 struct symbolic_type : private noncopyable { 55 void *ptr; symbolic_typeboost::numeric::bindings::umfpack::symbolic_type56 symbolic_type():ptr(0){} ~symbolic_typeboost::numeric::bindings::umfpack::symbolic_type57 ~symbolic_type() { 58 if (ptr) 59 detail::free_symbolic (T(), 0, &ptr); 60 } freeboost::numeric::bindings::umfpack::symbolic_type61 void free() { 62 if (ptr) 63 detail::free_symbolic (T(), 0, &ptr); 64 ptr = 0; 65 } 66 }; 67 68 template <typename T> free(symbolic_type<T> & s)69 void free (symbolic_type<T>& s) { s.free(); } 70 71 template <typename T = double> 72 struct numeric_type : private noncopyable { 73 void *ptr; numeric_typeboost::numeric::bindings::umfpack::numeric_type74 numeric_type():ptr(0){} ~numeric_typeboost::numeric::bindings::umfpack::numeric_type75 ~numeric_type() { 76 if (ptr) 77 detail::free_numeric (T(), 0, &ptr); 78 } freeboost::numeric::bindings::umfpack::numeric_type79 void free() { 80 if (ptr) 81 detail::free_numeric (T(), 0, &ptr); 82 ptr = 0; 83 } 84 }; 85 86 template <typename T> free(numeric_type<T> & n)87 void free (numeric_type<T>& n) { n.free(); } 88 89 90 template <typename T = double> 91 struct control_type : private noncopyable { 92 double ptr[UMFPACK_CONTROL]; control_typeboost::numeric::bindings::umfpack::control_type93 control_type() { detail::defaults (T(), 0, ptr); } operator []boost::numeric::bindings::umfpack::control_type94 double operator[] (int i) const { return ptr[i]; } operator []boost::numeric::bindings::umfpack::control_type95 double& operator[] (int i) { return ptr[i]; } defaultsboost::numeric::bindings::umfpack::control_type96 void defaults() { detail::defaults (T(), 0, ptr); } 97 }; 98 99 template <typename T> defaults(control_type<T> & c)100 void defaults (control_type<T>& c) { c.defaults(); } 101 102 template <typename T = double> 103 struct info_type : private noncopyable { 104 double ptr[UMFPACK_INFO]; operator []boost::numeric::bindings::umfpack::info_type105 double operator[] (int i) const { return ptr[i]; } operator []boost::numeric::bindings::umfpack::info_type106 double& operator[] (int i) { return ptr[i]; } 107 }; 108 109 110 ///////////////////////////////////// 111 // solving system of linear equations 112 ///////////////////////////////////// 113 114 115 // symbolic 116 /* 117 * Given nonzero pattern of a sparse matrix A in column-oriented form, 118 * umfpack_*_symbolic performs a column pre-ordering to reduce fill-in 119 * (using COLAMD or AMD) and a symbolic factorisation. This is required 120 * before the matrix can be numerically factorised with umfpack_*_numeric. 121 */ 122 namespace detail { 123 124 template <typename MatrA> 125 inline symbolic(tag::compressed_sparse,MatrA const & A,void ** Symbolic,double const * Control=0,double * Info=0)126 int symbolic (tag::compressed_sparse, 127 MatrA const& A, void **Symbolic, 128 double const* Control = 0, double* Info = 0) 129 { 130 return detail::symbolic (bindings::size_row (A), 131 bindings::size_column (A), 132 bindings::begin_compressed_index_major (A), 133 bindings::begin_index_minor (A), 134 bindings::begin_value (A), 135 Symbolic, Control, Info); 136 } 137 138 template <typename MatrA, typename QVec> 139 inline symbolic(tag::compressed_sparse,MatrA const & A,QVec const & Qinit,void ** Symbolic,double const * Control=0,double * Info=0)140 int symbolic (tag::compressed_sparse, 141 MatrA const& A, QVec const& Qinit, void **Symbolic, 142 double const* Control = 0, double* Info = 0) 143 { 144 #ifdef CHECK_TEST_COVERAGE 145 typedef typename MatrA::not_yet_tested i_m_still_here; 146 #endif 147 return detail::qsymbolic (bindings::size_row (A), 148 bindings::size_column (A), 149 bindings::begin_compressed_index_major (A), 150 bindings::begin_index_minor (A), 151 bindings::begin_value (A), 152 bindings::begin_value (Qinit), 153 Symbolic, Control, Info); 154 } 155 156 template <typename MatrA> 157 inline symbolic(tag::coordinate_sparse,MatrA const & A,void ** Symbolic,double const * Control=0,double * Info=0)158 int symbolic (tag::coordinate_sparse, 159 MatrA const& A, void **Symbolic, 160 double const* Control = 0, double* Info = 0) 161 { 162 int n_row = bindings::size_row (A); 163 int n_col = bindings::size_column (A); 164 int nnz = bindings::end_value (A) - bindings::begin_value (A); 165 166 typedef typename bindings::value_type<MatrA>::type val_t; 167 168 int const* Ti = bindings::begin_index_minor (A); 169 int const* Tj = bindings::begin_index_major (A); 170 bindings::detail::array<int> Ap (n_col+1); 171 if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory; 172 bindings::detail::array<int> Ai (nnz); 173 if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory; 174 175 int status = detail::triplet_to_col (n_row, n_col, nnz, 176 Ti, Tj, static_cast<val_t*> (0), 177 Ap.storage(), Ai.storage(), 178 static_cast<val_t*> (0), 0); 179 if (status != UMFPACK_OK) return status; 180 181 return detail::symbolic (n_row, n_col, 182 Ap.storage(), Ai.storage(), 183 bindings::begin_value (A), 184 Symbolic, Control, Info); 185 } 186 187 template <typename MatrA, typename QVec> 188 inline symbolic(tag::coordinate_sparse,MatrA const & A,QVec const & Qinit,void ** Symbolic,double const * Control=0,double * Info=0)189 int symbolic (tag::coordinate_sparse, 190 MatrA const& A, QVec const& Qinit, void **Symbolic, 191 double const* Control = 0, double* Info = 0) 192 { 193 #ifdef CHECK_TEST_COVERAGE 194 typedef typename MatrA::not_yet_tested i_m_still_here; 195 #endif 196 int n_row = bindings::size_row (A); 197 int n_col = bindings::size_column (A); 198 int nnz = bindings::end_value (A) - bindings::begin_value (A); 199 200 typedef typename bindings::value_type<MatrA>::type val_t; 201 202 int const* Ti = bindings::begin_index_minor (A); 203 int const* Tj = bindings::begin_index_major (A); 204 bindings::detail::array<int> Ap (n_col+1); 205 if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory; 206 bindings::detail::array<int> Ai (nnz); 207 if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory; 208 209 int status = detail::triplet_to_col (n_row, n_col, nnz, 210 Ti, Tj, static_cast<val_t*> (0), 211 Ap.storage(), Ai.storage(), 212 static_cast<val_t*> (0), 0); 213 if (status != UMFPACK_OK) return status; 214 215 return detail::qsymbolic (n_row, n_col, 216 Ap.storage(), Ai.storage(), 217 bindings::begin_value (A), 218 bindings::begin_value (Qinit), 219 Symbolic, Control, Info); 220 } 221 222 } // detail 223 224 template <typename MatrA> 225 inline symbolic(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,double const * Control=0,double * Info=0)226 int symbolic (MatrA const& A, 227 symbolic_type< 228 typename bindings::value_type<MatrA>::type 229 >& Symbolic, 230 double const* Control = 0, double* Info = 0) 231 { 232 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 233 check_umfpack_structure<MatrA>(); 234 #endif 235 typedef typename bindings::detail::property_at< 236 MatrA, tag::data_structure >::type storage_f; 237 238 return detail::symbolic (storage_f(), A, &Symbolic.ptr, Control, Info); 239 } 240 241 template <typename MatrA> 242 inline symbolic(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)243 int symbolic (MatrA const& A, 244 symbolic_type< 245 typename bindings::value_type<MatrA>::type 246 >& Symbolic, 247 control_type< 248 typename bindings::value_type<MatrA>::type 249 > const& Control, 250 info_type< 251 typename bindings::value_type<MatrA>::type 252 >& Info) 253 { 254 return symbolic (A, Symbolic, Control.ptr, Info.ptr); 255 } 256 257 template <typename MatrA> 258 inline symbolic(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,control_type<typename bindings::value_type<MatrA>::type> const & Control)259 int symbolic (MatrA const& A, 260 symbolic_type< 261 typename bindings::value_type<MatrA>::type 262 >& Symbolic, 263 control_type< 264 typename bindings::value_type<MatrA>::type 265 > const& Control) 266 { 267 return symbolic (A, Symbolic, Control.ptr); 268 } 269 270 template <typename MatrA, typename QVec> 271 inline symbolic(MatrA const & A,QVec const & Qinit,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,double const * Control=0,double * Info=0)272 int symbolic (MatrA const& A, QVec const& Qinit, 273 symbolic_type< 274 typename bindings::value_type<MatrA>::type 275 >& Symbolic, 276 double const* Control = 0, double* Info = 0) 277 { 278 #ifdef CHECK_TEST_COVERAGE 279 typedef typename MatrA::not_yet_tested i_m_still_here; 280 #endif 281 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 282 check_umfpack_structure<MatrA>(); 283 #endif 284 typedef typename bindings::detail::property_at< 285 MatrA, tag::data_structure >::type storage_f; 286 287 assert (bindings::size_column (A) == bindings::size (Qinit)); 288 289 return detail::symbolic (storage_f(), A, Qinit, 290 &Symbolic.ptr, Control, Info); 291 } 292 293 template <typename MatrA, typename QVec> 294 inline symbolic(MatrA const & A,QVec const & Qinit,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)295 int symbolic (MatrA const& A, QVec const& Qinit, 296 symbolic_type< 297 typename bindings::value_type<MatrA>::type 298 >& Symbolic, 299 control_type< 300 typename bindings::value_type<MatrA>::type 301 > const& Control, 302 info_type< 303 typename bindings::value_type<MatrA>::type 304 >& Info) 305 { 306 return symbolic (A, Qinit, Symbolic, Control.ptr, Info.ptr); 307 } 308 309 template <typename MatrA, typename QVec> 310 inline symbolic(MatrA const & A,QVec const & Qinit,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,control_type<typename bindings::value_type<MatrA>::type> const & Control)311 int symbolic (MatrA const& A, QVec const& Qinit, 312 symbolic_type< 313 typename bindings::value_type<MatrA>::type 314 >& Symbolic, 315 control_type< 316 typename bindings::value_type<MatrA>::type 317 > const& Control) 318 { 319 return symbolic (A, Qinit, Symbolic, Control.ptr); 320 } 321 322 323 // numeric 324 /* 325 * Given a sparse matrix A in column-oriented form, and a symbolic analysis 326 * computed by umfpack_*_*symbolic, the umfpack_*_numeric routine performs 327 * the numerical factorisation, PAQ=LU, PRAQ=LU, or P(R\A)Q=LU, where P 328 * and Q are permutation matrices (represented as permutation vectors), 329 * R is the row scaling, L is unit-lower triangular, and U is upper 330 * triangular. This is required before the system Ax=b (or other related 331 * linear systems) can be solved. 332 */ 333 namespace detail { 334 335 template <typename MatrA> 336 inline numeric(tag::compressed_sparse,MatrA const & A,void * Symbolic,void ** Numeric,double const * Control=0,double * Info=0)337 int numeric (tag::compressed_sparse, MatrA const& A, 338 void *Symbolic, void** Numeric, 339 double const* Control = 0, double* Info = 0) 340 { 341 return detail::numeric (bindings::size_row (A), 342 bindings::size_column (A), 343 bindings::begin_compressed_index_major (A), 344 bindings::begin_index_minor (A), 345 bindings::begin_value (A), 346 Symbolic, Numeric, Control, Info); 347 } 348 349 template <typename MatrA> 350 inline numeric(tag::coordinate_sparse,MatrA const & A,void * Symbolic,void ** Numeric,double const * Control=0,double * Info=0)351 int numeric (tag::coordinate_sparse, MatrA const& A, 352 void *Symbolic, void** Numeric, 353 double const* Control = 0, double* Info = 0) 354 { 355 int n_row = bindings::size_row (A); 356 int n_col = bindings::size_column (A); 357 int nnz = bindings::end_value (A) - bindings::begin_value (A); 358 359 typedef typename bindings::value_type<MatrA>::type val_t; 360 361 int const* Ti = bindings::begin_index_minor (A); 362 int const* Tj = bindings::begin_index_major (A); 363 bindings::detail::array<int> Ap (n_col+1); 364 if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory; 365 bindings::detail::array<int> Ai (nnz); 366 if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory; 367 368 int status = detail::triplet_to_col (n_row, n_col, nnz, 369 Ti, Tj, static_cast<val_t*> (0), 370 Ap.storage(), Ai.storage(), 371 static_cast<val_t*> (0), 0); 372 if (status != UMFPACK_OK) return status; 373 374 return detail::numeric (n_row, n_col, 375 Ap.storage(), Ai.storage(), 376 bindings::begin_value (A), 377 Symbolic, Numeric, Control, Info); 378 } 379 380 } // detail 381 382 template <typename MatrA> 383 inline numeric(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> const & Symbolic,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,double const * Control=0,double * Info=0)384 int numeric (MatrA const& A, 385 symbolic_type< 386 typename bindings::value_type<MatrA>::type 387 > const& Symbolic, 388 numeric_type< 389 typename bindings::value_type<MatrA>::type 390 >& Numeric, 391 double const* Control = 0, double* Info = 0) 392 { 393 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 394 check_umfpack_structure<MatrA>(); 395 #endif 396 typedef typename bindings::detail::property_at< 397 MatrA, tag::data_structure >::type storage_f; 398 399 return detail::numeric (storage_f(), A, 400 Symbolic.ptr, &Numeric.ptr, Control, Info); 401 } 402 403 template <typename MatrA> 404 inline numeric(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> const & Symbolic,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)405 int numeric (MatrA const& A, 406 symbolic_type< 407 typename bindings::value_type<MatrA>::type 408 > const& Symbolic, 409 numeric_type< 410 typename bindings::value_type<MatrA>::type 411 >& Numeric, 412 control_type< 413 typename bindings::value_type<MatrA>::type 414 > const& Control, 415 info_type< 416 typename bindings::value_type<MatrA>::type 417 >& Info) 418 419 { 420 // g++ (3.2) is unable to distinguish 421 // function numeric() and namespace boost::numeric ;o) 422 return umfpack::numeric (A, Symbolic, Numeric, Control.ptr, Info.ptr); 423 } 424 425 template <typename MatrA> 426 inline numeric(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> const & Symbolic,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control)427 int numeric (MatrA const& A, 428 symbolic_type< 429 typename bindings::value_type<MatrA>::type 430 > const& Symbolic, 431 numeric_type< 432 typename bindings::value_type<MatrA>::type 433 >& Numeric, 434 control_type< 435 typename bindings::value_type<MatrA>::type 436 > const& Control) 437 { 438 return umfpack::numeric (A, Symbolic, Numeric, Control.ptr); 439 } 440 441 442 // factor 443 /* 444 * symbolic and numeric 445 */ 446 namespace detail { 447 448 template <typename MatrA> 449 inline factor(tag::compressed_sparse,MatrA const & A,void ** Numeric,double const * Control=0,double * Info=0)450 int factor (tag::compressed_sparse, MatrA const& A, 451 void** Numeric, double const* Control = 0, double* Info = 0) 452 { 453 #ifdef CHECK_TEST_COVERAGE 454 typedef typename MatrA::not_yet_tested i_m_still_here; 455 #endif 456 symbolic_type<typename bindings::value_type<MatrA>::type> 457 Symbolic; 458 459 int status; 460 status = detail::symbolic (bindings::size_row (A), 461 bindings::size_column (A), 462 bindings::begin_compressed_index_major (A), 463 bindings::begin_index_minor (A), 464 bindings::begin_value (A), 465 &Symbolic.ptr, Control, Info); 466 if (status != UMFPACK_OK) return status; 467 468 return detail::numeric (bindings::size_row (A), 469 bindings::size_column (A), 470 bindings::begin_compressed_index_major (A), 471 bindings::begin_index_minor (A), 472 bindings::begin_value (A), 473 Symbolic.ptr, Numeric, Control, Info); 474 } 475 476 template <typename MatrA> 477 inline factor(tag::coordinate_sparse,MatrA const & A,void ** Numeric,double const * Control=0,double * Info=0)478 int factor (tag::coordinate_sparse, MatrA const& A, 479 void** Numeric, double const* Control = 0, double* Info = 0) 480 { 481 #ifdef CHECK_TEST_COVERAGE 482 typedef typename MatrA::not_yet_tested i_m_still_here; 483 #endif 484 int n_row = bindings::size_row (A); 485 int n_col = bindings::size_column (A); 486 int nnz = bindings::end_value (A) - bindings::begin_value (A); 487 488 typedef typename bindings::value_type<MatrA>::type val_t; 489 490 int const* Ti = bindings::begin_index_minor (A); 491 int const* Tj = bindings::begin_index_major (A); 492 bindings::detail::array<int> Ap (n_col+1); 493 if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory; 494 bindings::detail::array<int> Ai (nnz); 495 if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory; 496 497 int status = detail::triplet_to_col (n_row, n_col, nnz, 498 Ti, Tj, static_cast<val_t*> (0), 499 Ap.storage(), Ai.storage(), 500 static_cast<val_t*> (0), 0); 501 if (status != UMFPACK_OK) return status; 502 503 symbolic_type<typename bindings::value_type<MatrA>::type> 504 Symbolic; 505 506 status = detail::symbolic (n_row, n_col, 507 Ap.storage(), Ai.storage(), 508 bindings::begin_value (A), 509 &Symbolic.ptr, Control, Info); 510 if (status != UMFPACK_OK) return status; 511 512 return detail::numeric (n_row, n_col, 513 Ap.storage(), Ai.storage(), 514 bindings::begin_value (A), 515 Symbolic.ptr, Numeric, Control, Info); 516 } 517 518 } // detail 519 520 template <typename MatrA> 521 inline factor(MatrA const & A,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,double const * Control=0,double * Info=0)522 int factor (MatrA const& A, 523 numeric_type< 524 typename bindings::value_type<MatrA>::type 525 >& Numeric, 526 double const* Control = 0, double* Info = 0) 527 { 528 #ifdef CHECK_TEST_COVERAGE 529 typedef typename MatrA::not_yet_tested i_m_still_here; 530 #endif 531 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 532 check_umfpack_structure<MatrA>(); 533 #endif 534 typedef typename bindings::detail::property_at< 535 MatrA, tag::data_structure >::type storage_f; 536 537 return detail::factor (storage_f(), A, &Numeric.ptr, Control, Info); 538 } 539 540 template <typename MatrA> 541 inline factor(MatrA const & A,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)542 int factor (MatrA const& A, 543 numeric_type< 544 typename bindings::value_type<MatrA>::type 545 >& Numeric, 546 control_type< 547 typename bindings::value_type<MatrA>::type 548 > const& Control, 549 info_type< 550 typename bindings::value_type<MatrA>::type 551 >& Info) 552 { 553 return factor (A, Numeric, Control.ptr, Info.ptr); 554 } 555 556 template <typename MatrA> 557 inline factor(MatrA const & A,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control)558 int factor (MatrA const& A, 559 numeric_type< 560 typename bindings::value_type<MatrA>::type 561 >& Numeric, 562 control_type< 563 typename bindings::value_type<MatrA>::type 564 > const& Control) 565 { 566 return factor (A, Numeric, Control.ptr); 567 } 568 569 570 // solve 571 /* 572 * Given LU factors computed by umfpack_*_numeric and the right-hand-side, 573 * B, solve a linear system for the solution X. Iterative refinement is 574 * optionally performed. Only square systems are handled. 575 */ 576 namespace detail { 577 578 template <typename MatrA, typename VecX, typename VecB> 579 inline solve(tag::compressed_sparse,int sys,MatrA const & A,VecX & X,VecB const & B,void * Numeric,double const * Control=0,double * Info=0)580 int solve (tag::compressed_sparse, int sys, 581 MatrA const& A, VecX& X, VecB const& B, 582 void *Numeric, double const* Control = 0, double* Info = 0) 583 { 584 return detail::solve (sys, bindings::size_row (A), 585 bindings::begin_compressed_index_major (A), 586 bindings::begin_index_minor (A), 587 bindings::begin_value (A), 588 bindings::begin_value (X), 589 bindings::begin_value (B), 590 Numeric, Control, Info); 591 } 592 593 template <typename MatrA, typename VecX, typename VecB> 594 inline solve(tag::coordinate_sparse,int sys,MatrA const & A,VecX & X,VecB const & B,void * Numeric,double const * Control=0,double * Info=0)595 int solve (tag::coordinate_sparse, int sys, 596 MatrA const& A, VecX& X, VecB const& B, 597 void *Numeric, double const* Control = 0, double* Info = 0) 598 { 599 600 int n = bindings::size_row (A); 601 int nnz = bindings::end_value (A) - bindings::begin_value (A); 602 603 typedef typename bindings::value_type<MatrA>::type val_t; 604 605 int const* Ti = bindings::begin_index_minor (A); 606 int const* Tj = bindings::begin_index_major (A); 607 bindings::detail::array<int> Ap (n+1); 608 if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory; 609 bindings::detail::array<int> Ai (nnz); 610 if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory; 611 612 int status = detail::triplet_to_col (n, n, nnz, 613 Ti, Tj, static_cast<val_t*> (0), 614 Ap.storage(), Ai.storage(), 615 static_cast<val_t*> (0), 0); 616 if (status != UMFPACK_OK) return status; 617 618 return detail::solve (sys, n, Ap.storage(), Ai.storage(), 619 bindings::begin_value (A), 620 bindings::begin_value (X), 621 bindings::begin_value (B), 622 Numeric, Control, Info); 623 } 624 625 } // detail 626 627 template <typename MatrA, typename VecX, typename VecB> 628 inline solve(int sys,MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,double const * Control=0,double * Info=0)629 int solve (int sys, MatrA const& A, VecX& X, VecB const& B, 630 numeric_type< 631 typename bindings::value_type<MatrA>::type 632 > const& Numeric, 633 double const* Control = 0, double* Info = 0) 634 { 635 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 636 check_umfpack_structure<MatrA>(); 637 #endif 638 typedef typename bindings::detail::property_at< 639 MatrA, tag::data_structure >::type storage_f; 640 641 assert (bindings::size_row (A) == bindings::size_row (A)); 642 assert (bindings::size_column (A) == bindings::size (X)); 643 assert (bindings::size_column (A) == bindings::size (B)); 644 645 return detail::solve (storage_f(), sys, A, X, B, 646 Numeric.ptr, Control, Info); 647 } 648 649 template <typename MatrA, typename VecX, typename VecB> 650 inline solve(int sys,MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)651 int solve (int sys, MatrA const& A, VecX& X, VecB const& B, 652 numeric_type< 653 typename bindings::value_type<MatrA>::type 654 > const& Numeric, 655 control_type< 656 typename bindings::value_type<MatrA>::type 657 > const& Control, 658 info_type< 659 typename bindings::value_type<MatrA>::type 660 >& Info) 661 { 662 return solve (sys, A, X, B, Numeric, Control.ptr, Info.ptr); 663 } 664 665 template <typename MatrA, typename VecX, typename VecB> 666 inline solve(int sys,MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control)667 int solve (int sys, MatrA const& A, VecX& X, VecB const& B, 668 numeric_type< 669 typename bindings::value_type<MatrA>::type 670 > const& Numeric, 671 control_type< 672 typename bindings::value_type<MatrA>::type 673 > const& Control) 674 { 675 return solve (sys, A, X, B, Numeric, Control.ptr); 676 } 677 678 template <typename MatrA, typename VecX, typename VecB> 679 inline solve(MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,double const * Control=0,double * Info=0)680 int solve (MatrA const& A, VecX& X, VecB const& B, 681 numeric_type< 682 typename bindings::value_type<MatrA>::type 683 > const& Numeric, 684 double const* Control = 0, double* Info = 0) 685 { 686 return solve (UMFPACK_A, A, X, B, Numeric, Control, Info); 687 } 688 689 template <typename MatrA, typename VecX, typename VecB> 690 inline solve(MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)691 int solve (MatrA const& A, VecX& X, VecB const& B, 692 numeric_type< 693 typename bindings::value_type<MatrA>::type 694 > const& Numeric, 695 control_type< 696 typename bindings::value_type<MatrA>::type 697 > const& Control, 698 info_type< 699 typename bindings::value_type<MatrA>::type 700 >& Info) 701 { 702 return solve (UMFPACK_A, A, X, B, Numeric, 703 Control.ptr, Info.ptr); 704 } 705 706 template <typename MatrA, typename VecX, typename VecB> 707 inline solve(MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control)708 int solve (MatrA const& A, VecX& X, VecB const& B, 709 numeric_type< 710 typename bindings::value_type<MatrA>::type 711 > const& Numeric, 712 control_type< 713 typename bindings::value_type<MatrA>::type 714 > const& Control) 715 { 716 return solve (UMFPACK_A, A, X, B, Numeric, Control.ptr); 717 } 718 719 720 // umf_solve 721 /* 722 * symbolic, numeric and solve 723 */ 724 namespace detail { 725 726 template <typename MatrA, typename VecX, typename VecB> 727 inline umf_solve(tag::compressed_sparse,MatrA const & A,VecX & X,VecB const & B,double const * Control=0,double * Info=0)728 int umf_solve (tag::compressed_sparse, 729 MatrA const& A, VecX& X, VecB const& B, 730 double const* Control = 0, double* Info = 0) 731 { 732 #ifdef CHECK_TEST_COVERAGE 733 typedef typename MatrA::not_yet_tested i_m_still_here; 734 #endif 735 symbolic_type<typename bindings::value_type<MatrA>::type> 736 Symbolic; 737 numeric_type<typename bindings::value_type<MatrA>::type> 738 Numeric; 739 740 int status; 741 status = detail::symbolic (bindings::size_row (A), 742 bindings::size_column (A), 743 bindings::begin_compressed_index_major (A), 744 bindings::begin_index_minor (A), 745 bindings::begin_value (A), 746 &Symbolic.ptr, Control, Info); 747 if (status != UMFPACK_OK) return status; 748 749 status = detail::numeric (bindings::size_row (A), 750 bindings::size_column (A), 751 bindings::begin_compressed_index_major (A), 752 bindings::begin_index_minor (A), 753 bindings::begin_value (A), 754 Symbolic.ptr, &Numeric.ptr, Control, Info); 755 if (status != UMFPACK_OK) return status; 756 757 return detail::solve (UMFPACK_A, bindings::size_row (A), 758 bindings::begin_compressed_index_major (A), 759 bindings::begin_index_minor (A), 760 bindings::begin_value (A), 761 bindings::begin_value (X), 762 bindings::begin_value (B), 763 Numeric.ptr, Control, Info); 764 } 765 766 template <typename MatrA, typename VecX, typename VecB> 767 inline umf_solve(tag::coordinate_sparse,MatrA const & A,VecX & X,VecB const & B,double const * Control=0,double * Info=0)768 int umf_solve (tag::coordinate_sparse, 769 MatrA const& A, VecX& X, VecB const& B, 770 double const* Control = 0, double* Info = 0) 771 { 772 #ifdef CHECK_TEST_COVERAGE 773 typedef typename MatrAA::not_yet_tested i_m_still_here; 774 #endif 775 int n_row = bindings::size_row (A); 776 int n_col = bindings::size_column (A); 777 int nnz = bindings::end_value (A) - bindings::begin_value (A); 778 779 typedef typename bindings::value_type<MatrA>::type val_t; 780 781 int const* Ti = bindings::begin_index_minor (A); 782 int const* Tj = bindings::begin_index_major (A); 783 bindings::detail::array<int> Ap (n_col+1); 784 if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory; 785 bindings::detail::array<int> Ai (nnz); 786 if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory; 787 788 int status = detail::triplet_to_col (n_row, n_col, nnz, 789 Ti, Tj, static_cast<val_t*> (0), 790 Ap.storage(), Ai.storage(), 791 static_cast<val_t*> (0), 0); 792 if (status != UMFPACK_OK) return status; 793 794 symbolic_type<typename bindings::value_type<MatrA>::type> 795 Symbolic; 796 numeric_type<typename bindings::value_type<MatrA>::type> 797 Numeric; 798 799 status = detail::symbolic (n_row, n_col, 800 Ap.storage(), Ai.storage(), 801 bindings::begin_value (A), 802 &Symbolic.ptr, Control, Info); 803 if (status != UMFPACK_OK) return status; 804 805 status = detail::numeric (n_row, n_col, 806 Ap.storage(), Ai.storage(), 807 bindings::begin_value (A), 808 Symbolic.ptr, &Numeric.ptr, Control, Info); 809 if (status != UMFPACK_OK) return status; 810 811 return detail::solve (UMFPACK_A, n_row, Ap.storage(), Ai.storage(), 812 bindings::begin_value (A), 813 bindings::begin_value (X), 814 bindings::begin_value (B), 815 Numeric.ptr, Control, Info); 816 } 817 818 } // detail 819 820 template <typename MatrA, typename VecX, typename VecB> 821 inline umf_solve(MatrA const & A,VecX & X,VecB const & B,double const * Control=0,double * Info=0)822 int umf_solve (MatrA const& A, VecX& X, VecB const& B, 823 double const* Control = 0, double* Info = 0) 824 { 825 #ifdef CHECK_TEST_COVERAGE 826 typedef typename MatrA::not_yet_tested i_m_still_here; 827 #endif 828 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 829 check_umfpack_structure<MatrA>(); 830 #endif 831 typedef typename bindings::detail::property_at< 832 MatrA, tag::data_structure >::type storage_f; 833 834 assert (bindings::size_row (A) == bindings::size_row (A)); 835 assert (bindings::size_column (A) == bindings::size (X)); 836 assert (bindings::size_column (A) == bindings::size (B)); 837 838 return detail::umf_solve (storage_f(), A, X, B, Control, Info); 839 } 840 841 template <typename MatrA, typename VecX, typename VecB> 842 inline umf_solve(MatrA const & A,VecX & X,VecB const & B,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)843 int umf_solve (MatrA const& A, VecX& X, VecB const& B, 844 control_type< 845 typename bindings::value_type<MatrA>::type 846 > const& Control, 847 info_type< 848 typename bindings::value_type<MatrA>::type 849 >& Info) 850 { 851 return umf_solve (A, X, B, Control.ptr, Info.ptr); 852 } 853 854 template <typename MatrA, typename VecX, typename VecB> 855 inline umf_solve(MatrA const & A,VecX & X,VecB const & B,control_type<typename bindings::value_type<MatrA>::type> const & Control)856 int umf_solve (MatrA const& A, VecX& X, VecB const& B, 857 control_type< 858 typename bindings::value_type<MatrA>::type 859 > const& Control) 860 { 861 return umf_solve (A, X, B, Control.ptr); 862 } 863 864 865 /////////////////////// 866 // matrix manipulations 867 /////////////////////// 868 869 870 // scale 871 872 template <typename VecX, typename VecB> 873 inline scale(VecX & X,VecB const & B,numeric_type<typename bindings::value_type<VecB>::type> const & Numeric)874 int scale (VecX& X, VecB const& B, 875 numeric_type< 876 typename bindings::value_type<VecB>::type 877 > const& Numeric) 878 { 879 return detail::scale (bindings::size (B), 880 bindings::begin_value (X), 881 bindings::begin_value (B), 882 Numeric.ptr); 883 } 884 885 886 //////////// 887 // reporting 888 //////////// 889 890 891 // report status 892 893 template <typename T> 894 inline report_status(control_type<T> const & Control,int status)895 void report_status (control_type<T> const& Control, int status) { 896 detail::report_status (T(), 0, Control.ptr, status); 897 } 898 899 #if 0 900 template <typename T> 901 inline 902 void report_status (int printing_level, int status) { 903 control_type<T> Control; 904 Control[UMFPACK_PRL] = printing_level; 905 detail::report_status (T(), 0, Control.ptr, status); 906 } 907 template <typename T> 908 inline 909 void report_status (int status) { 910 control_type<T> Control; 911 detail::report_status (T(), 0, Control.ptr, status); 912 } 913 #endif 914 915 916 // report control 917 918 template <typename T> 919 inline report_control(control_type<T> const & Control)920 void report_control (control_type<T> const& Control) { 921 detail::report_control (T(), 0, Control.ptr); 922 } 923 924 925 // report info 926 927 template <typename T> 928 inline report_info(control_type<T> const & Control,info_type<T> const & Info)929 void report_info (control_type<T> const& Control, info_type<T> const& Info) { 930 detail::report_info (T(), 0, Control.ptr, Info.ptr); 931 } 932 933 #if 0 934 template <typename T> 935 inline 936 void report_info (int printing_level, info_type<T> const& Info) { 937 control_type<T> Control; 938 Control[UMFPACK_PRL] = printing_level; 939 detail::report_info (T(), 0, Control.ptr, Info.ptr); 940 } 941 template <typename T> 942 inline 943 void report_info (info_type<T> const& Info) { 944 control_type<T> Control; 945 detail::report_info (T(), 0, Control.ptr, Info.ptr); 946 } 947 #endif 948 949 950 // report matrix (compressed column and coordinate) 951 952 namespace detail { 953 954 template <typename MatrA> 955 inline report_matrix(tag::compressed_sparse,MatrA const & A,double const * Control)956 int report_matrix (tag::compressed_sparse, MatrA const& A, 957 double const* Control) 958 { 959 return detail::report_matrix (bindings::size_row (A), 960 bindings::size_column (A), 961 bindings::begin_compressed_index_major (A), 962 bindings::begin_index_minor (A), 963 bindings::begin_value (A), 964 1, Control); 965 } 966 967 template <typename MatrA> 968 inline report_matrix(tag::coordinate_sparse,MatrA const & A,double const * Control)969 int report_matrix (tag::coordinate_sparse, MatrA const& A, 970 double const* Control) 971 { 972 return detail::report_triplet (bindings::size_row (A), 973 bindings::size_column (A), 974 bindings::end_value (A) - bindings::begin_value (A), 975 bindings::begin_index_major (A), 976 bindings::begin_index_minor (A), 977 bindings::begin_value (A), 978 Control); 979 } 980 981 } // detail 982 983 template <typename MatrA> 984 inline report_matrix(MatrA const & A,control_type<typename bindings::value_type<MatrA>::type> const & Control)985 int report_matrix (MatrA const& A, 986 control_type< 987 typename bindings::value_type<MatrA>::type 988 > const& Control) 989 { 990 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 991 check_umfpack_structure<MatrA>(); 992 #endif 993 typedef typename bindings::detail::property_at< 994 MatrA, tag::data_structure >::type storage_f; 995 996 return detail::report_matrix (storage_f(), A, Control.ptr); 997 } 998 999 1000 // report vector 1001 1002 template <typename VecX> 1003 inline report_vector(VecX const & X,control_type<typename bindings::value_type<VecX>::type> const & Control)1004 int report_vector (VecX const& X, 1005 control_type< 1006 typename bindings::value_type<VecX>::type 1007 > const& Control) 1008 { 1009 return detail::report_vector (bindings::size (X), 1010 bindings::begin_value (X), 1011 Control.ptr); 1012 } 1013 1014 1015 // report numeric 1016 1017 template <typename T> 1018 inline report_numeric(numeric_type<T> const & Numeric,control_type<T> const & Control)1019 int report_numeric (numeric_type<T> const& Numeric, 1020 control_type<T> const& Control) 1021 { 1022 return detail::report_numeric (T(), 0, Numeric.ptr, Control.ptr); 1023 } 1024 1025 1026 // report symbolic 1027 1028 template <typename T> 1029 inline report_symbolic(symbolic_type<T> const & Symbolic,control_type<T> const & Control)1030 int report_symbolic (symbolic_type<T> const& Symbolic, 1031 control_type<T> const& Control) 1032 { 1033 return detail::report_symbolic (T(), 0, Symbolic.ptr, Control.ptr); 1034 } 1035 1036 1037 // report permutation vector 1038 1039 template <typename VecP, typename T> 1040 inline report_permutation(VecP const & Perm,control_type<T> const & Control)1041 int report_permutation (VecP const& Perm, control_type<T> const& Control) { 1042 #ifdef CHECK_TEST_COVERAGE 1043 typedef typename T::not_yet_tested i_m_still_here; 1044 #endif 1045 return detail::report_perm (T(), 0, 1046 bindings::begin_value (Perm), 1047 Control.ptr); 1048 } 1049 1050 1051 }}}} 1052 1053 #endif // BOOST_NUMERIC_BINDINGS_UMFPACK_HPP 1054