Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_Superlumt_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// Amesos2: Templated Direct Sparse Solver Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//
42// @HEADER
43
52#ifndef AMESOS2_SUPERLUMT_DEF_HPP
53#define AMESOS2_SUPERLUMT_DEF_HPP
54
55#include <Teuchos_Tuple.hpp>
56#include <Teuchos_ParameterList.hpp>
57#include <Teuchos_StandardParameterEntryValidators.hpp>
58
60#include "Amesos2_Util.hpp"
61
62
63namespace SLUMT {
64 /*
65 * We have to declare this extern function because of a bug in the
66 * SuperLU_MT header files. In each header is declared a function
67 * "extern int xsp_colorder(...)" where x is in {'s','d','c','z'},
68 * but this function is never defined anywhere, so if you use the
69 * function as defined, it will compile fine, but will choke during
70 * linking. No other code in SuperLU_MT actually calls these
71 * functions. Instead, all other SuperLU_MT function just call
72 * "sp_colorder", which is defined within the SuperLU_MT library,
73 * but not declared.
74 */
75 extern "C" {
76 int sp_colorder(SuperMatrix*,int*,superlumt_options_t*,SuperMatrix*);
77 }
78} // end namespace SLUMT
79
80
81namespace Amesos2 {
82
83 /*
84 * Note: Although many of the type definitions for SuperLU_MT are
85 * identical to those of SuperLU, we do not mix the definitions so
86 * that we do not introduce messy coupling between the two
87 * interfaces. Also, there exist enough differences between the two
88 * packages to merit dedicated utility functions.
89 *
90 * We have also taken a different approach to interfacing with
91 * SuperLU_MT than with SuperLU which I believe leads to a better
92 * seperation of the 4 solution steps. We may in the future adopt a
93 * similar strategy for SuperLU.
94 */
95
96 template <class Matrix, class Vector>
97 Superlumt<Matrix,Vector>::Superlumt(Teuchos::RCP<const Matrix> A,
98 Teuchos::RCP<Vector> X,
99 Teuchos::RCP<const Vector> B )
100 : SolverCore<Amesos2::Superlumt,Matrix,Vector>(A, X, B)
101 , nzvals_() // initialize to empty arrays
102 , rowind_()
103 , colptr_()
104 , is_contiguous_(true)
105 {
106 Teuchos::RCP<Teuchos::ParameterList> default_params
107 = Teuchos::parameterList( *(this->getValidParameters()) );
108 this->setParameters(default_params);
109
110 data_.options.lwork = 0; // Use system memory for factorization
111
112 data_.perm_r.resize(this->globalNumRows_);
113 data_.perm_c.resize(this->globalNumCols_);
114 data_.options.perm_r = data_.perm_r.getRawPtr();
115 data_.options.perm_c = data_.perm_c.getRawPtr();
116
117 data_.R.resize(this->globalNumRows_);
118 data_.C.resize(this->globalNumCols_);
119
120 data_.options.refact = SLUMT::NO; // initially we are not refactoring
121 data_.equed = SLUMT::NOEQUIL; // No equilibration has yet been performed
122 data_.rowequ = false;
123 data_.colequ = false;
124
125 data_.A.Store = NULL;
126 data_.AC.Store = NULL;
127 data_.BX.Store = NULL;
128 data_.L.Store = NULL;
129 data_.U.Store = NULL;
130
131 data_.stat.ops = NULL;
132 }
133
134
135 template <class Matrix, class Vector>
137 {
138 /* Free SuperLU_MT data_types
139 * - Matrices
140 * - Vectors
141 * - Stat object
142 */
143
144 // Storage is initialized in numericFactorization_impl()
145 if ( data_.A.Store != NULL ){
146 // Our Teuchos::Array's will destroy rowind, colptr, and nzval for us
147 SLUMT::D::Destroy_SuperMatrix_Store( &(data_.A) );
148 }
149
150 // Cleanup memory allocated during last call to sp_colorder if needed
151 if( data_.AC.Store != NULL ){
152 SLUMT::D::Destroy_CompCol_Permuted( &(data_.AC) ); // free's colbeg, colend, and Store
153 }
154
155 if ( data_.L.Store != NULL ){ // will only ever be true for this->root_
156 SLUMT::D::Destroy_SuperNode_SCP( &(data_.L) );
157 SLUMT::D::Destroy_CompCol_NCP( &(data_.U) );
158
159 // memory alloc'd in sp_colorder
160 free( data_.options.etree );
161 free( data_.options.colcnt_h );
162 free( data_.options.part_super_h );
163 }
164
165
166 // Storage is initialized in solve_impl()
167 if ( data_.BX.Store != NULL ){
168 /* Cannot use SLU::Destroy_Dense_Matrix routine here, since it attempts to
169 * free the array of non-zero values, but that array has already been
170 * deallocated by the MultiVector object. So we release just the Store
171 * instead.
172 */
173 SLUMT::D::Destroy_SuperMatrix_Store( &(data_.BX) );
174 }
175
176 if ( data_.stat.ops != NULL )
177 SLUMT::D::StatFree( &(data_.stat) );
178 }
179
180 template<class Matrix, class Vector>
181 int
183 {
184 // Use either the column-ordering found in the users perm_c or the requested computed ordering
185 int perm_spec = data_.options.ColPerm;
186 if( perm_spec != SLUMT::MY_PERMC && this->root_ ){
187#ifdef HAVE_AMESOS2_TIMERS
188 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
189#endif
190
191 SLUMT::S::get_perm_c(perm_spec, &(data_.A), data_.perm_c.getRawPtr());
192 }
193 // Ordering will be applied directly before numeric factorization
194
195 return(0);
196 }
197
198
199
200 template <class Matrix, class Vector>
201 int
203 {
204 // We assume that a call to symbolicFactorization is indicating that
205 // the structure of the matrix has changed. SuperLU doesn't allow
206 // us to do a symbolic factorization directly, but we can leave a
207 // flag that tells it it needs to symbolically factor later.
208 data_.options.refact = SLUMT::NO;
209
210 if( this->status_.numericFactorizationDone() && this->root_ ){
211 // If we've done a numeric factorization already, then we need to
212 // cleanup the old L and U. Stores and other data will be
213 // allocated during numeric factorization. Only rank 0 has valid
214 // pointers
215 SLUMT::D::Destroy_SuperNode_Matrix( &(data_.L) );
216 SLUMT::D::Destroy_CompCol_NCP( &(data_.U) );
217 data_.L.Store = NULL;
218 data_.U.Store = NULL;
219 }
220
221 return(0);
222 }
223
224
225 template <class Matrix, class Vector>
226 int
228 {
229 using Teuchos::as;
230
231#ifdef HAVE_AMESOS2_DEBUG
232 const int nprocs = data_.options.nprocs;
233 TEUCHOS_TEST_FOR_EXCEPTION( nprocs <= 0,
234 std::invalid_argument,
235 "The number of threads to spawn should be greater than 0." );
236#endif
237
238 int info = 0;
239
240 if ( this->root_ ) {
241
242 if( data_.options.fact == SLUMT::EQUILIBRATE ){
243 magnitude_type rowcnd, colcnd, amax;
244 int info;
245
246 function_map::gsequ(&(data_.A), data_.R.getRawPtr(),
247 data_.C.getRawPtr(), &rowcnd, &colcnd,
248 &amax, &info);
249 TEUCHOS_TEST_FOR_EXCEPTION( info != 0,
250 std::runtime_error,
251 "SuperLU_MT gsequ returned with status " << info );
252
253 function_map::laqgs(&(data_.A), data_.R.getRawPtr(),
254 data_.C.getRawPtr(), rowcnd, colcnd,
255 amax, &(data_.equed));
256
257 data_.rowequ = (data_.equed == SLUMT::ROW) || (data_.equed == SLUMT::BOTH);
258 data_.colequ = (data_.equed == SLUMT::COL) || (data_.equed == SLUMT::BOTH);
259
260 data_.options.fact = SLUMT::DOFACT;
261 }
262
263 // Cleanup memory allocated during last call to sp_colorder if needed
264 if( data_.AC.Store != NULL ){
265 SLUMT::D::Destroy_CompCol_Permuted( &(data_.AC) ); // free's colbeg, colend, and Store
266 if( data_.options.refact == SLUMT::NO ){ // then we recompute etree; free the old one
267 free( data_.options.etree );
268 free( data_.options.colcnt_h );
269 free( data_.options.part_super_h );
270 }
271 data_.AC.Store = NULL;
272 }
273
274 // Apply the column ordering, so that AC is the column-permuted A, and compute etree
275 SLUMT::sp_colorder(&(data_.A), data_.perm_c.getRawPtr(),
276 &(data_.options), &(data_.AC));
277
278
279 // Allocate and initialize status variable
280 const int n = as<int>(this->globalNumCols_); // n is the number of columns in A
281 if( data_.stat.ops != NULL ){ SLUMT::D::StatFree( &(data_.stat) ); data_.stat.ops = NULL; }
282 SLUMT::D::StatAlloc(n, data_.options.nprocs, data_.options.panel_size,
283 data_.options.relax, &(data_.stat));
284 SLUMT::D::StatInit(n, data_.options.nprocs, &(data_.stat));
285
286
287 { // Do factorization
288#ifdef HAVE_AMESOS2_TIMERS
289 Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
290#endif
291
292#ifdef HAVE_AMESOS2_VERBOSE_DEBUG
293 std::cout << "SuperLU_MT:: Before numeric factorization" << std::endl;
294 std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
295 std::cout << "rowind_ : " << rowind_.toString() << std::endl;
296 std::cout << "colptr_ : " << colptr_.toString() << std::endl;
297#endif
298
299 function_map::gstrf(&(data_.options), &(data_.AC),
300 data_.perm_r.getRawPtr(), &(data_.L), &(data_.U),
301 &(data_.stat), &info);
302 }
303
304 // Set the number of non-zero values in the L and U factors
305 this->setNnzLU( as<size_t>(((SLUMT::SCformat*)data_.L.Store)->nnz +
306 ((SLUMT::NCformat*)data_.U.Store)->nnz) );
307 }
308
309 // Check output
310 const global_size_type info_st = as<global_size_type>(info);
311 TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st <= this->globalNumCols_),
312 std::runtime_error,
313 "Factorization complete, but matrix is singular. Division by zero eminent");
314 TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st > this->globalNumCols_),
315 std::runtime_error,
316 "Memory allocation failure in SuperLU_MT factorization");
317 // The other option, that info_st < 0 denotes invalid parameters to
318 // the function, but we'll assume for now that that won't happen.
319
320 data_.options.fact = SLUMT::FACTORED;
321 data_.options.refact = SLUMT::YES;
322
323 /* All processes should return the same error code */
324 Teuchos::broadcast(*(this->getComm()),0,&info);
325 return(info);
326 }
327
328
329 template <class Matrix, class Vector>
330 int
331 Superlumt<Matrix,Vector>::solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
332 const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
333 {
334 using Teuchos::as;
335
336 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
337 const size_t nrhs = X->getGlobalNumVectors();
338
339 Teuchos::Array<slu_type> bxvals_(ld_rhs * nrhs);
340 size_t ldbx_ = as<size_t>(ld_rhs);
341
342 { // Get values from B
343#ifdef HAVE_AMESOS2_TIMERS
344 Teuchos::TimeMonitor convTimer( this->timers_.vecConvTime_ );
345 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
346#endif
347
348 if ( is_contiguous_ == true ) {
350 slu_type>::do_get(B, bxvals_(),
351 ldbx_,
352 ROOTED, this->rowIndexBase_);
353 }
354 else {
355 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
356 "NON_CONTIGUOUS_AND_ROOTED not supported.");
357 }
358 }
359
360 int info = 0; // returned error code (0 = success)
361 if ( this->root_ ) {
362 // Clean up old B stores if it has already been created
363 if( data_.BX.Store != NULL ){
364 SLUMT::D::Destroy_SuperMatrix_Store( &(data_.BX) );
365 data_.BX.Store = NULL;
366 }
367
368 {
369#ifdef HAVE_AMESOS2_TIMERS
370 Teuchos::TimeMonitor convTimer( this->timers_.vecConvTime_ );
371#endif
372 SLUMT::Dtype_t dtype = type_map::dtype;
373 function_map::create_Dense_Matrix(&(data_.BX), as<int>(ld_rhs), as<int>(nrhs),
374 bxvals_.getRawPtr(), as<int>(ldbx_),
375 SLUMT::SLU_DN, dtype, SLUMT::SLU_GE);
376 }
377
378 if( data_.options.trans == SLUMT::NOTRANS ){
379 if( data_.rowequ ){ // row equilibration has been done on AC
380 // scale bxvals_ by diag(R)
381 Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.R(),
382 SLUMT::slu_mt_mult<slu_type,magnitude_type>());
383 }
384 } else if( data_.colequ ){ // column equilibration has been done on AC
385 // scale bxvals_ by diag(C)
386 Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.C(),
387 SLUMT::slu_mt_mult<slu_type,magnitude_type>());
388 }
389
390
391 {
392#ifdef HAVE_AMESOS2_TIMERS
393 Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
394#endif
395
396 function_map::gstrs(data_.options.trans, &(data_.L),
397 &(data_.U), data_.perm_r.getRawPtr(),
398 data_.perm_c.getRawPtr(), &(data_.BX),
399 &(data_.stat), &info);
400 }
401 } // end block for solve time
402
403 /* All processes should have the same error code */
404 Teuchos::broadcast(*(this->getComm()),0,&info);
405
406 TEUCHOS_TEST_FOR_EXCEPTION( info < 0,
407 std::runtime_error,
408 "Argument " << -info << " to gstrs had an illegal value" );
409
410 // "Un-scale" the solution so that it is a solution of the original system
411 if( data_.options.trans == SLUMT::NOTRANS ){
412 if( data_.colequ ){ // column equilibration has been done on AC
413 // scale bxvals_ by diag(C)
414 Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.C(),
415 SLUMT::slu_mt_mult<slu_type,magnitude_type>());
416 }
417 } else if( data_.rowequ ){ // row equilibration has been done on AC
418 // scale bxvals_ by diag(R)
419 Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.R(),
420 SLUMT::slu_mt_mult<slu_type,magnitude_type>());
421 }
422
423 /* Update X's global values */
424 {
425#ifdef HAVE_AMESOS2_TIMERS
426 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
427#endif
428
429 if ( is_contiguous_ == true ) {
431 MultiVecAdapter<Vector>, slu_type>::do_put(X, bxvals_(), ldbx_, ROOTED);
432 }
433 else {
434 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
435 "NON_CONTIGUOUS_AND_ROOTED not supported.");
436 }
437 }
438
439 return(info);
440 }
441
442
443 template <class Matrix, class Vector>
444 bool
446 {
447 // The SuperLU_MT factorization routines can handle square as well as
448 // rectangular matrices, but SuperLU_MT can only apply the solve routines to
449 // square matrices, so we check the matrix for squareness.
450 return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
451 }
452
453
454 template <class Matrix, class Vector>
455 void
456 Superlumt<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
457 {
458 using Teuchos::as;
459 using Teuchos::RCP;
460 using Teuchos::getIntegralValue;
461 using Teuchos::ParameterEntryValidator;
462
463 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
464
465 int default_nprocs = Kokkos::DefaultHostExecutionSpace().concurrency();
466 data_.options.nprocs = parameterList->get<int>("nprocs", default_nprocs);
467
468 data_.options.trans = this->control_.useTranspose_ ? SLUMT::TRANS : SLUMT::NOTRANS;
469 // SuperLU_MT "trans" option can override the Amesos2 option
470 if( parameterList->isParameter("trans") ){
471 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("trans").validator();
472 parameterList->getEntry("trans").setValidator(trans_validator);
473
474 data_.options.trans = getIntegralValue<SLUMT::trans_t>(*parameterList, "trans");
475 }
476
477 data_.options.panel_size = parameterList->get<int>("panel_size", SLUMT::D::sp_ienv(1));
478
479 data_.options.relax = parameterList->get<int>("relax", SLUMT::D::sp_ienv(2));
480
481 const bool equil = parameterList->get<bool>("Equil", true);
482 data_.options.fact = equil ? SLUMT::EQUILIBRATE : SLUMT::DOFACT;
483
484 const bool symmetric_mode = parameterList->get<bool>("SymmetricMode", false);
485 data_.options.SymmetricMode = symmetric_mode ? SLUMT::YES : SLUMT::NO;
486
487 const bool print_stat = parameterList->get<bool>("PrintStat", false);
488 data_.options.PrintStat = print_stat ? SLUMT::YES : SLUMT::NO;
489
490 data_.options.diag_pivot_thresh = parameterList->get<double>("diag_pivot_thresh", 1.0);
491
492 if( parameterList->isParameter("ColPerm") ){
493 RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry("ColPerm").validator();
494 parameterList->getEntry("ColPerm").setValidator(colperm_validator);
495
496 data_.options.ColPerm = getIntegralValue<SLUMT::colperm_t>(*parameterList, "ColPerm");
497 }
498
499 // TODO: until we support retrieving col/row permutations from the user
500 data_.options.usepr = SLUMT::NO;
501
502 if( parameterList->isParameter("IsContiguous") ){
503 is_contiguous_ = parameterList->get<bool>("IsContiguous");
504 }
505 }
506
507
508 template <class Matrix, class Vector>
509 Teuchos::RCP<const Teuchos::ParameterList>
511 {
512 using std::string;
513 using Teuchos::tuple;
514 using Teuchos::ParameterList;
515 using Teuchos::EnhancedNumberValidator;
516 using Teuchos::setStringToIntegralParameter;
517 using Teuchos::stringToIntegralParameterEntryValidator;
518
519 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
520
521 if( is_null(valid_params) ){
522 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
523
524 int default_nprocs = Kokkos::DefaultHostExecutionSpace().concurrency();
525 Teuchos::RCP<EnhancedNumberValidator<int> > nprocs_validator
526 = Teuchos::rcp( new EnhancedNumberValidator<int>() );
527 nprocs_validator->setMin(default_nprocs);
528 pl->set("nprocs", default_nprocs, "The number of processors to factorize with", nprocs_validator);
529
530 setStringToIntegralParameter<SLUMT::trans_t>("trans", "NOTRANS",
531 "Solve for the transpose system or not",
532 tuple<string>("TRANS","NOTRANS","CONJ"),
533 tuple<string>("Solve with transpose",
534 "Do not solve with transpose",
535 "Solve with the conjugate transpose"),
536 tuple<SLUMT::trans_t>(SLUMT::TRANS,
537 SLUMT::NOTRANS,
538 SLUMT::CONJ),
539 pl.getRawPtr());
540
541 Teuchos::RCP<EnhancedNumberValidator<int> > panel_size_validator
542 = Teuchos::rcp( new EnhancedNumberValidator<int>() );
543 panel_size_validator->setMin(1);
544 pl->set("panel_size", SLUMT::D::sp_ienv(1),
545 "Specifies the number of consecutive columns to be treated as a unit of task.",
546 panel_size_validator);
547
548 Teuchos::RCP<EnhancedNumberValidator<int> > relax_validator
549 = Teuchos::rcp( new EnhancedNumberValidator<int>() );
550 relax_validator->setMin(1);
551 pl->set("relax", SLUMT::D::sp_ienv(2),
552 "Specifies the number of columns to be grouped as a relaxed supernode",
553 relax_validator);
554
555 pl->set("Equil", true, "Whether to equilibrate the system before solve");
556
557 Teuchos::RCP<EnhancedNumberValidator<double> > diag_pivot_thresh_validator
558 = Teuchos::rcp( new EnhancedNumberValidator<double>(0.0, 1.0) );
559 pl->set("diag_pivot_thresh", 1.0,
560 "Specifies the threshold used for a diagonal entry to be an acceptable pivot",
561 diag_pivot_thresh_validator); // partial pivoting
562
563 // Note: MY_PERMC not yet supported
564 setStringToIntegralParameter<SLUMT::colperm_t>("ColPerm", "COLAMD",
565 "Specifies how to permute the columns of the "
566 "matrix for sparsity preservation",
567 tuple<string>("NATURAL",
568 "MMD_AT_PLUS_A",
569 "MMD_ATA",
570 "COLAMD"),
571 tuple<string>("Natural ordering",
572 "Minimum degree ordering on A^T + A",
573 "Minimum degree ordering on A^T A",
574 "Approximate minimum degree column ordering"),
575 tuple<SLUMT::colperm_t>(SLUMT::NATURAL,
576 SLUMT::MMD_AT_PLUS_A,
577 SLUMT::MMD_ATA,
578 SLUMT::COLAMD),
579 pl.getRawPtr());
580
581 pl->set("SymmetricMode", false, "Specifies whether to use the symmetric mode");
582
583 // TODO: until we support getting row/col permutations from user
584 // pl->set("usepr", false);
585
586 pl->set("PrintStat", false, "Specifies whether to print the solver's statistics");
587
588 pl->set("IsContiguous", true, "Whether GIDs contiguous");
589
590 valid_params = pl;
591 }
592
593 return valid_params;
594 }
595
596
597 template <class Matrix, class Vector>
598 bool
600 {
601 using Teuchos::as;
602
603#ifdef HAVE_AMESOS2_TIMERS
604 Teuchos::TimeMonitor convTimer( this->timers_.mtxConvTime_ );
605#endif
606
607 if( current_phase == SYMBFACT ) return false;
608
609 // Store is allocated on create_CompCol_Matrix
610 if( data_.A.Store != NULL ){
611 SLUMT::D::Destroy_SuperMatrix_Store( &(data_.A) );
612 data_.A.Store = NULL;
613 }
614
615 if( this->root_ ){
616 nzvals_.resize(this->globalNumNonZeros_);
617 rowind_.resize(this->globalNumNonZeros_);
618 colptr_.resize(this->globalNumCols_ + 1);
619 }
620
621 int nnz_ret = 0;
622 {
623#ifdef HAVE_AMESOS2_TIMERS
624 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
625#endif
626
627 // Use compressed-column store for SuperLU_MT
628 if ( is_contiguous_ == true ) {
629 Util::get_ccs_helper<
630 MatrixAdapter<Matrix>,slu_type,int,int>::do_get(this->matrixA_.ptr(),
631 nzvals_, rowind_, colptr_,
632 nnz_ret, ROOTED, ARBITRARY, this->rowIndexBase_);
633 }
634 else {
635 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
636 "NON_CONTIGUOUS_AND_ROOTED not supported.");
637 }
638 }
639
640 if( this->root_ ){
641 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<int>(this->globalNumNonZeros_),
642 std::runtime_error,
643 "rank=0 failed to get all nonzero values in getCcs()");
644
645 SLUMT::Dtype_t dtype = type_map::dtype;
646 function_map::create_CompCol_Matrix(&(data_.A),
647 as<int>(this->globalNumRows_),
648 as<int>(this->globalNumCols_),
649 nnz_ret,
650 nzvals_.getRawPtr(),
651 rowind_.getRawPtr(),
652 colptr_.getRawPtr(),
653 SLUMT::SLU_NC,
654 dtype, SLUMT::SLU_GE);
655
656 TEUCHOS_TEST_FOR_EXCEPTION( data_.A.Store == NULL,
657 std::runtime_error,
658 "Failed to allocate matrix store" );
659 }
660
661 return true;
662 }
663
664
665 template<class Matrix, class Vector>
666 const char* Superlumt<Matrix,Vector>::name = "SuperLU_MT";
667
668
669} // end namespace Amesos2
670
671#endif // AMESOS2_SUPERLUMT_DEF_HPP
@ ROOTED
Definition Amesos2_TypeDecl.hpp:127
@ ARBITRARY
Definition Amesos2_TypeDecl.hpp:143
Utility functions for Amesos2.
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers.
Definition Amesos2_SolverCore_decl.hpp:106
super_type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Set/update internal variables and solver options.
Definition Amesos2_SolverCore_def.hpp:518
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a const parameter list of all of the valid parameters that this->setParameterList(....
Definition Amesos2_SolverCore_def.hpp:550
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition Amesos2_SolverCore_decl.hpp:479
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition Amesos2_SolverCore_decl.hpp:476
Amesos2 interface to the Multi-threaded version of SuperLU.
Definition Amesos2_Superlumt_decl.hpp:78
Superlumt(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition Amesos2_Superlumt_def.hpp:97
~Superlumt()
Destructor.
Definition Amesos2_Superlumt_def.hpp:136
EPhase
Used to indicate a phase in the direct solution.
Definition Amesos2_TypeDecl.hpp:65
Helper class for getting 1-D copies of multivectors.
Definition Amesos2_MultiVecAdapter_decl.hpp:267
Helper class for putting 1-D data arrays into multivectors.
Definition Amesos2_MultiVecAdapter_decl.hpp:373