FEI Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
fei_Aztec_LinSysCore.cpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ************************************************************************
4// FEI: Finite Element Interface to Linear Solvers
5// Copyright (2005) Sandia Corporation.
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8// U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Alan Williams (william@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41*/
42
43
45#include <fei_sstream.hpp>
46
47#ifdef HAVE_FEI_AZTECOO
48
49#include <stdlib.h>
50#include <string.h>
51#include <stdio.h>
52#include <stdexcept>
53
54#include "fei_defs.h"
55#include "fei_Data.hpp"
56#include "fei_Lookup.hpp"
58#include "snl_fei_Utils.hpp"
59#include "fei_ArrayUtils.hpp"
60
61#undef fei_file
62#define fei_file "fei_Aztec_LinSysCore.cpp"
63#include "fei_ErrMacros.hpp"
64
65#include "fei_mpi.h"
66
67#ifndef FEI_SER
68
69#define AZTEC_MPI AZTEC_MPI
70#define AZ_MPI AZ_MPI
71#ifndef MPI
72#define MPI MPI
73#endif
74
75#endif
76
77#include <az_aztec.h>
78
79#include "fei_Aztec_Map.hpp"
84
86
87namespace fei_trilinos {
88
89static int azlsc_solveCounter_ = 0;
90static int azlsc_debugFileCounter_ = 0;
91
92static std::map<std::string, unsigned> fei_aztec_named_solve_counter;
93
94//=========CONSTRUCTOR==========================================================
96 : comm_(comm),
97 lookup_(NULL),
98 haveLookup_(false),
99 update_(NULL),
100 map_(),
101 A_(NULL),
102 A_ptr_(NULL),
103 x_(NULL),
104 b_(NULL),
105 bc_(NULL),
106 essBCindices_(NULL),
107 numEssBCs_(0),
108 bcsLoaded_(false),
109 explicitDirichletBCs_(false),
110 BCenforcement_no_column_mod_(false),
111 b_ptr_(NULL),
112 matrixAllocated_(false),
113 vectorsAllocated_(false),
114 blkMatrixAllocated_(false),
115 matrixLoaded_(false),
116 rhsLoaded_(false),
117 needNewPreconditioner_(false),
118 tooLateToChooseBlock_(false),
119 blockMatrix_(false),
120 blkMap_(),
121 blkA_(NULL),
122 blkA_ptr_(NULL),
123 blkUpdate_(NULL),
124 azA_(NULL),
125 azP_(NULL),
126 precondCreated_(false),
127 azS_(NULL),
128 scalingCreated_(false),
129 aztec_options_(NULL),
130 aztec_params_(NULL),
131 aztec_status_(NULL),
132 tmp_x_(NULL),
133 tmp_x_touched_(false),
134 tmp_b_(NULL),
135 tmp_bc_(NULL),
136 tmp_b_allocated_(false),
137 ML_Vanek_(false),
138 numLevels_(9), //arbitrary default...
139 rhsIDs_(NULL),
140 numRHSs_(0),
141 currentRHS_(-1),
142 numGlobalEqns_(0),
143 localOffset_(0),
144 numLocalEqns_(0),
145 numGlobalEqnBlks_(0),
146 localBlkOffset_(0),
147 numLocalEqnBlks_(0),
148 localBlockSizes_(NULL),
149 outputLevel_(0),
150 numParams_(0),
151 paramStrings_(NULL),
152 name_("dbg"),
153 debugOutput_(0),
154 debugFileCounter_(0),
155 debugPath_(NULL),
156 debugFileName_(NULL),
157 debugFile_(NULL),
158 named_solve_counter_(fei_aztec_named_solve_counter)
159{
160 masterProc_ = 0;
161 numProcs_ = 1;
162 thisProc_ = 0;
163#ifndef FEI_SER
164 MPI_Comm_size(comm_, &numProcs_);
165 MPI_Comm_rank(comm_, &thisProc_);
166#endif
167
168 aztec_options_ = new int[AZ_OPTIONS_SIZE];
169 aztec_params_ = new double[AZ_PARAMS_SIZE];
170
171 AZ_defaults(aztec_options_, aztec_params_);
172
173 aztec_status_ = new double[AZ_STATUS_SIZE];
174 for(int i=0; i<AZ_STATUS_SIZE; i++) aztec_status_[i] = 0.0;
175
176 numRHSs_ = 1;
177 rhsIDs_ = new int[numRHSs_];
178 rhsIDs_[0] = 0;
179
180 std::map<std::string,unsigned>::iterator
181 iter = named_solve_counter_.find(name_);
182 if (iter == named_solve_counter_.end()) {
183 named_solve_counter_.insert(std::make_pair(name_, 0));
184 }
185}
186
187struct free_any_remaining_az_memory {
188 free_any_remaining_az_memory(){}
189 ~free_any_remaining_az_memory()
190 {
191 AZ_manage_memory(0, AZ_CLEAR_ALL, 0, NULL, NULL);
192 }
193};
194
195//========DESTRUCTOR============================================================
196Aztec_LinSysCore::~Aztec_LinSysCore() {
197 static free_any_remaining_az_memory when_program_exits;
198
199 if (blkMatrixAllocated_) {
200 delete blkA_;
201 delete [] blkUpdate_;
202 delete [] localBlockSizes_;
203 blkMatrixAllocated_ = false;
204 }
205 if (matrixAllocated_) {
206 delete A_;
207 matrixAllocated_ = false;
208 }
209
210 if (precondCreated_) {
211 AZ_precond_destroy(&azP_);
212 precondCreated_ = false;
213 }
214
215 if (scalingCreated_) {
216 AZ_scaling_destroy(&azS_);
217 scalingCreated_ = false;
218 }
219
220 if (vectorsAllocated_) delete x_;
221 delete [] tmp_x_;
222 delete [] tmp_bc_;
223
224 for(int i=0; i<numRHSs_; i++) {
225 if (vectorsAllocated_) delete b_[i];
226 if (tmp_b_allocated_) delete [] tmp_b_[i];
227 }
228 if (vectorsAllocated_) delete [] b_;
229 if (tmp_b_allocated_) delete [] tmp_b_;
230 tmp_b_allocated_ = false;
231 delete [] update_;
232
233 delete [] aztec_options_;
234 delete [] aztec_params_;
235 delete [] aztec_status_;
236
237 delete [] rhsIDs_;
238 numRHSs_ = 0;
239
240 for(int j=0; j<numParams_; j++) {
241 delete [] paramStrings_[j];
242 }
243 delete [] paramStrings_;
244 numParams_ = 0;
245
246 if (debugOutput_) {
247 debugOutput_ = 0;
248 fclose(debugFile_);
249 delete [] debugPath_;
250 delete [] debugFileName_;
251 }
252
253 delete [] essBCindices_;
254 essBCindices_ = NULL;
255 numEssBCs_ = 0;
256 delete bc_;
257}
258
259//==============================================================================
260LinearSystemCore* Aztec_LinSysCore::clone() {
261 return(new Aztec_LinSysCore(comm_));
262}
263
264//==============================================================================
265int Aztec_LinSysCore::parameters(int numParams, const char*const * params) {
266//
267// this function takes parameters for setting internal things like solver
268// and preconditioner choice, etc.
269//
270 debugOutput("parameters");
271
272 if (numParams == 0 || params == NULL) {
273 debugOutput("--- no parameters");
274 return(0);
275 }
276
277 const char* param = NULL;
278
279 snl_fei::mergeStringLists(paramStrings_, numParams_,
280 params, numParams);
281
282 param = snl_fei::getParamValue("outputLevel",numParams,params);
283 if (param != NULL){
284 sscanf(param,"%d", &outputLevel_);
285 }
286
287 param = snl_fei::getParamValue("AZ_matrix_type", numParams, params);
288 if (param != NULL){
289 setMatrixType(param);
290 }
291
292 param = snl_fei::getParam("EXPLICIT_BC_ENFORCEMENT", numParams, params);
293 if (param != NULL){
294 explicitDirichletBCs_ = true;
295 }
296 else {
297 explicitDirichletBCs_ = false;
298 }
299
300 param = snl_fei::getParam("BC_ENFORCEMENT_NO_COLUMN_MOD",numParams,params);
301 if (param != NULL){
302 BCenforcement_no_column_mod_ = true;
303 }
304 else {
305 BCenforcement_no_column_mod_ = false;
306 }
307
308 param = snl_fei::getParamValue("numLevels", numParams, params);
309 if (param != NULL){
310 sscanf(param,"%d", &numLevels_);
311 }
312
313 bool dbgOutputParam = false;
314 char* dbgFileName = NULL;
315 char* dbgPath = NULL;
316
317 bool output_level_on = false;
318 param = snl_fei::getParamValue("FEI_OUTPUT_LEVEL",numParams,params);
319 if (param != NULL) {
320 std::string str(param);
321 if (str == "ALL" || str == "MATRIX_FILES" || str == "FULL_LOGS") {
322 output_level_on = true;
323 }
324 }
325
326 param = snl_fei::getParamValue("debugOutput",numParams,params);
327 if (param != NULL || output_level_on){
328 dbgOutputParam = true;
329 dbgFileName = new char[128];
330 dbgPath = param==NULL ? new char[3] : new char[strlen(param)+1];
331 if (param == NULL) sprintf(dbgPath, "./");
332 else strcpy(dbgPath, param);
333
334 sprintf(dbgFileName, "AZLSC.%d.%d.%d",
335 numProcs_, thisProc_, azlsc_debugFileCounter_);
336 }
337
338 param = snl_fei::getParamValue("name",numParams, params);
339 if(param != NULL){
340 name_ = param;
341
342 std::map<std::string,unsigned>::iterator
343 iter = named_solve_counter_.find(name_);
344 if (iter == named_solve_counter_.end()) {
345 named_solve_counter_.insert(std::make_pair(name_, 0));
346 }
347
348 if (dbgOutputParam) {
349 if (dbgFileName != NULL) delete [] dbgFileName;
350 dbgFileName = new char[256];
351 sprintf(dbgFileName, "AZLSC_%s.%d.%d.%d", name_.c_str(),
352 numProcs_, thisProc_, azlsc_debugFileCounter_);
353 }
354 }
355
356 if (dbgOutputParam) {
357 if (azlsc_debugFileCounter_ == azlsc_solveCounter_) {
358 setDebugOutput(dbgPath,dbgFileName);
359 ++azlsc_debugFileCounter_;
360 }
361 delete [] dbgFileName;
362 delete [] dbgPath;
363 }
364
365 if (debugOutput_) {
366 fprintf(debugFile_,"--- numParams %d\n",numParams);
367 for(int i=0; i<numParams; i++){
368 fprintf(debugFile_,"------ paramStrings[%d]: %s\n",i,
369 paramStrings_[i]);
370 }
371 }
372
373 debugOutput("leaving parameters function");
374 return(0);
375}
376
377//==============================================================================
378int Aztec_LinSysCore::setLookup(Lookup& lookup)
379{
380 lookup_ = &lookup;
381 haveLookup_ = true;
382 return(0);
383}
384
385//==============================================================================
386int Aztec_LinSysCore::setGlobalOffsets(int len, int* nodeOffsets,
387 int* eqnOffsets, int* blkEqnOffsets)
388{
389 localOffset_ = eqnOffsets[thisProc_];
390 numLocalEqns_ = eqnOffsets[thisProc_+1]-localOffset_;
391 numGlobalEqns_ = eqnOffsets[numProcs_];
392
393 localBlkOffset_ = blkEqnOffsets[thisProc_];
394 numLocalEqnBlks_ = blkEqnOffsets[thisProc_+1]-localBlkOffset_;
395 numGlobalEqnBlks_ = blkEqnOffsets[numProcs_];
396
397 int err = createMiscStuff();
398 if (err) return(err);
399
400 if (debugOutput_) {
401 fprintf(debugFile_, "setGlobalNodeOffsets, len: %d\n", len);
402 for(int i=0; i<len; i++) {
403 fprintf(debugFile_, " nodeOffsets[%d]: %d, eqnOffsets[%d]: %d, blkEqnOffsets[%d], %d\n", i, nodeOffsets[i], i, eqnOffsets[i], i, blkEqnOffsets[i]);
404 }
405 fflush(debugFile_);
406 }
407 return(0);
408}
409
410//==============================================================================
411int Aztec_LinSysCore::setConnectivities(GlobalID elemBlock,
412 int numElements,
413 int numNodesPerElem,
414 const GlobalID* elemIDs,
415 const int* const* connNodes)
416{
417 (void) elemBlock;
418 (void) numElements;
419 (void) numNodesPerElem;
420 (void) elemIDs;
421 (void) connNodes;
422 return(0);
423}
424
425//==============================================================================
426int Aztec_LinSysCore::setMatrixType(const char* name) {
427
428 if (!strcmp(name, "AZ_VBR_MATRIX")) {
429 if (!tooLateToChooseBlock_) {
430 blockMatrix_ = true;
431 debugOutput("setMatrixType: chose block matrix");
432 }
433 else {
434 fei::console_out() << "Aztec_LinSysCore::setMatrixType: WARNING: Too late to choose"
435 << " the DVBR matrix; make this choice before calling "
436 << "setMatrixStructure. DMSR will be used." << FEI_ENDL;
437 }
438 }
439 else if (!strcmp(name, "AZ_MSR_MATRIX")) {
440 //do nothing, this is the default case.
441 }
442 else {
443 if (thisProc_ == 0) {
444 FEI_COUT << "Aztec_LinSysCore: Warning: requested matrix-type <"<<name <<"> not recognized." << FEI_ENDL;
445 FEI_COUT << "Aztec_LinSysCore: valid matrix-types are: 'AZ_MSR_MATRIX', 'AZ_VBR_MATRIX'" << FEI_ENDL;
446 }
447 }
448 return(0);
449}
450
451//==============================================================================
452int Aztec_LinSysCore::setMatrixStructure(int** ptColIndices,
453 int* ptRowLengths,
454 int** blkColIndices,
455 int* blkRowLengths,
456 int* ptRowsPerBlkRow)
457{
458 if (debugOutput_) {
459 fprintf(debugFile_, "setMatrixStructure\n");
460 for(int i=0; i<numLocalEqnBlks_; i++) {
461 int blkEqn = localBlkOffset_+i;
462 fprintf(debugFile_, " blkRow %d, ptRowsPerBlkRow %d\n",
463 blkEqn, ptRowsPerBlkRow[i]);
464 }
465 fflush(debugFile_);
466 }
467
468 int err = allocateMatrix(ptColIndices, ptRowLengths, blkColIndices,
469 blkRowLengths, ptRowsPerBlkRow);
470 return(err);
471}
472
473//==============================================================================
474int Aztec_LinSysCore::createMiscStuff()
475{
476//
477//This function is where we establish the structures/objects associated
478//with the linear algebra library. i.e., do initial allocations, etc.
479//
480 if (debugOutput_) {
481 fprintf(debugFile_,
482 "createMiscStuff: numRHSs_: %d\n", numRHSs_);
483 fflush(debugFile_);
484 }
485
486 if (numLocalEqns_ > numGlobalEqns_)
487 messageAbort("createMiscStuff: numLocalEqns > numGlobalEqns");
488
489 if (0 > localOffset_)
490 messageAbort("createMiscStuff: localOffset_ out of range");
491
492 update_ = numLocalEqns_ > 0 ? new int[numLocalEqns_] : NULL;
493 if (numLocalEqns_ > 0 && update_ == NULL) {
494 return(-1);
495 }
496
497 int i, j;
498 for(j=0; j<numLocalEqns_; j++) update_[j] = localOffset_+j;
499
500 azS_ = AZ_scaling_create();
501 if (azS_ != NULL) scalingCreated_ = true;
502 else {
503 fei::console_out() << "Aztec_LinSysCore::createMiscStuff ERROR: failed to create scaling"
504 << FEI_ENDL;
505 return(-1);
506 }
507
508 if (numRHSs_ <= 0)
509 messageAbort("numRHSs_==0. Out of scope or destructor already called?");
510
511 tmp_x_ = numLocalEqns_ > 0 ? new double[numLocalEqns_] : NULL;
512 if (numLocalEqns_ > 0 && tmp_x_ == NULL) return(-1);
513
514 tmp_bc_ = numLocalEqns_ > 0 ? new double[numLocalEqns_] : NULL;
515 if (numLocalEqns_ > 0 && tmp_bc_ == NULL) return(-1);
516
517 for(i=0; i<numLocalEqns_; i++){
518 tmp_x_[i] = 0.0;
519 tmp_bc_[i] = 0.0;
520 }
521
522 if (!tmp_b_allocated_) {
523 tmp_b_ = new double*[numRHSs_];
524
525 for(j=0; j<numRHSs_; j++) {
526 tmp_b_[j] = new double[numLocalEqns_];
527 for(i=0; i<numLocalEqns_; i++) {
528 tmp_b_[j][i] = 0.0;
529 }
530 }
531
532 tmp_b_allocated_ = true;
533 }
534
535 if (currentRHS_ < 0) currentRHS_ = 0;
536 return(0);
537}
538
539//==============================================================================
540int Aztec_LinSysCore::allocateMatrix(int** ptColIndices,
541 int* ptRowLengths,
542 int** blkColIndices,
543 int* blkRowLengths,
544 int* ptRowsPerBlkRow)
545{
546 int err;
547 if (blockMatrix_) {
548 err = createBlockMatrix(blkColIndices, blkRowLengths, ptRowsPerBlkRow);
549 return(err);
550 }
551
552 tooLateToChooseBlock_ = true;
553
554 map_.reset( new Aztec_Map(numGlobalEqns_, numLocalEqns_, update_, localOffset_, comm_));
555
556 A_ = new AztecDMSR_Matrix(map_);
557
558 if (A_ptr_ == NULL) A_ptr_ = A_;
559
560 //
561 //we're going to have to look through the colIndices and see if any
562 //of the rows do NOT have an entry on the diagonal. For each row that
563 //does have an entry on the diagonal, we subtract 1 from the row-length
564 //since the AztecDMSR_Matrix::allocate function wants the length of each
565 //row NOT including any entry on the diagonal.
566 //
567
568 int* row_lengths = numLocalEqns_ > 0 ? new int[numLocalEqns_] : NULL;
569
570 for (int i = 0; i < numLocalEqns_; i++) {
571 if (fei::searchList(localOffset_+i,
572 ptColIndices[i], ptRowLengths[i]) >= 0) {
573 row_lengths[i] = ptRowLengths[i] - 1;
574 }
575 else {
576 row_lengths[i] = ptRowLengths[i];
577 }
578 }
579
580 // so now we know all the row lengths, and can configure our matrix.
581
582 A_ptr_->allocate(row_lengths, ptColIndices);
583
584 delete [] row_lengths;
585
586 matrixAllocated_ = true;
587 return(0);
588}
589
590//==============================================================================
591int Aztec_LinSysCore::createBlockMatrix(int** blkColIndices,
592 int* blkRowLengths,
593 int* ptRowsPerBlkRow)
594{
595 int i;
596
597 blkUpdate_ = new int[numLocalEqnBlks_];
598
599 for(i=0; i<numLocalEqnBlks_; i++) {
600 blkUpdate_[i] = localBlkOffset_ + i;
601 }
602
603 blkMap_.reset(new Aztec_BlockMap(numGlobalEqns_, numLocalEqns_,
604 update_, localOffset_, comm_,
605 numGlobalEqnBlks_, numLocalEqnBlks_,
606 blkUpdate_,
607 localBlkOffset_, ptRowsPerBlkRow));
608
609 blkA_ = new AztecDVBR_Matrix(blkMap_);
610
611 if (blkA_ptr_ == NULL) blkA_ptr_ = blkA_;
612
613 localBlockSizes_ = new int[numLocalEqnBlks_];
614
615 //now we need to count how many non-zero blocks there are.
616 numNonzeroBlocks_ = 0;
617 for(i=0; i<numLocalEqnBlks_; i++) {
618
619 numNonzeroBlocks_ += blkRowLengths[i];
620
621 localBlockSizes_[i] = ptRowsPerBlkRow[i];
622 }
623
624 //and now we need to flatten the list of block-column-indices into a 1-D
625 //list to give to the AztecDVBR_Matrix::allocate function.
626 int* blk_col_inds = new int[numNonzeroBlocks_];
627
628 int offset = 0;
629 for(i=0; i<numLocalEqnBlks_; i++) {
630 for(int j=0; j<blkRowLengths[i]; j++) {
631 blk_col_inds[offset++] = blkColIndices[i][j];
632 }
633 }
634
635 //finally we're ready to allocate our matrix.
636 blkA_->allocate(blkRowLengths, blk_col_inds);
637
638 delete [] blk_col_inds;
639
640 blkMatrixAllocated_ = true;
641 return(0);
642}
643
644//==============================================================================
645int Aztec_LinSysCore::resetMatrixAndVector(double s)
646{
647 if (blkMatrixAllocated_) blkA_ptr_->put(s);
648 if (matrixAllocated_) A_ptr_->put(s);
649
650 if (rhsLoaded_) {
651 for(int i=0; i<numRHSs_; i++){
652 b_[i]->put(s);
653 }
654 }
655 if (b_ptr_ != NULL) b_ptr_->put(s);
656
657 if (bc_ != NULL) bc_->put(s);
658
659 if (!tmp_b_allocated_) {
660 for(int j=0; j<numRHSs_; j++) {
661 tmp_b_[j] = new double[numLocalEqns_];
662 }
663 }
664 matrixLoaded_ = false;
665 rhsLoaded_ = false;
666 bcsLoaded_ = false;
667
668 delete [] tmp_bc_;
669 tmp_bc_ = new double[numLocalEqns_];
670
671 for(int ii=0; ii<numLocalEqns_; ii++) tmp_bc_[ii] = s;
672
673 delete [] essBCindices_;
674 essBCindices_ = NULL;
675 numEssBCs_ = 0;
676
677 for(int j=0; j<numRHSs_; j++) {
678 for(int i=0; i<numLocalEqns_; i++) {
679 tmp_b_[j][i] = s;
680 }
681 }
682 return(0);
683}
684
685//==============================================================================
686int Aztec_LinSysCore::resetMatrix(double s) {
687 if (blkMatrixAllocated_) blkA_ptr_->put(s);
688 if (matrixAllocated_) A_ptr_->put(s);
689
690 matrixLoaded_ = false;
691 bcsLoaded_ = false;
692
693 if (bc_ != NULL) bc_->put(s);
694
695 delete [] tmp_bc_;
696 tmp_bc_ = new double[numLocalEqns_];
697
698 for(int ii=0; ii<numLocalEqns_; ii++) tmp_bc_[ii] = s;
699
700 delete [] essBCindices_;
701 essBCindices_ = NULL;
702 numEssBCs_ = 0;
703
704 return(0);
705}
706
707//==============================================================================
708int Aztec_LinSysCore::resetRHSVector(double s) {
709
710 if (rhsLoaded_) {
711 for(int i=0; i<numRHSs_; i++){
712 b_[i]->put(s);
713 }
714 }
715
716 if (b_ptr_ != NULL) b_ptr_->put(s);
717
718 if (!tmp_b_allocated_) {
719 for(int j=0; j<numRHSs_; j++) {
720 tmp_b_[j] = new double[numLocalEqns_];
721 }
722 }
723 rhsLoaded_ = false;
724 bcsLoaded_ = false;
725
726 for(int j=0; j<numRHSs_; j++) {
727 double* cur_b = tmp_b_[j];
728 for(int i=0; i<numLocalEqns_; i++) {
729 cur_b[i] = s;
730 }
731 }
732 return(0);
733}
734
735//==============================================================================
736int Aztec_LinSysCore::sumIntoSystemMatrix(int numPtRows, const int* ptRows,
737 int numPtCols, const int* ptCols,
738 const double* const* values)
739{
740 matrixLoaded_ = false;
741
742 return(sumIntoPointRow(numPtRows, ptRows, numPtCols, ptCols, values, false));
743}
744
745//==============================================================================
746int Aztec_LinSysCore::sumIntoSystemMatrix(int numPtRows, const int* ptRows,
747 int numPtCols, const int* ptCols,
748 int numBlkRows, const int* blkRows,
749 int numBlkCols, const int* blkCols,
750 const double* const* values)
751{
752 matrixLoaded_ = false;
753
754 if ((A_ptr_ == NULL) && (blkA_ptr_ == NULL))
755 messageAbort("sumIntoSystemMatrix: matrix is NULL.");
756
757 if (blockMatrix_) {
758 return(sumIntoBlockRow(numBlkRows, blkRows, numBlkCols, blkCols,
759 values, numPtCols, false));
760 }
761 else {
762 return( sumIntoPointRow(numPtRows, ptRows, numPtCols, ptCols, values,
763 false));
764 }
765}
766
767//==============================================================================
768int Aztec_LinSysCore::putIntoSystemMatrix(int numPtRows, const int* ptRows,
769 int numPtCols, const int* ptCols,
770 const double* const* values)
771{
772 matrixLoaded_ = false;
773
774 return(sumIntoPointRow(numPtRows, ptRows, numPtCols, ptCols, values, true));
775}
776
777//==============================================================================
778int Aztec_LinSysCore::getMatrixRowLength(int row, int& length)
779{
780 if (!blockMatrix_) {
781 if (row >= localOffset_ && row < (localOffset_+numLocalEqns_)) {
782 length = A_ptr_->rowLength(row);
783 return(0);
784 }
785 }
786 else {
787 if (!haveLookup_) return(-1);
788 int blkRow = lookup_->ptEqnToBlkEqn(row);
789 if (blkRow < 0) return(-1);
790 return( blkA_ptr_->getNumNonzerosPerRow(blkRow, length) );
791 }
792
793 return(-1);
794}
795
796//==============================================================================
797int Aztec_LinSysCore::getMatrixRow(int row, double* coefs,
798 int* indices,
799 int len, int& rowLength)
800{
801 int err = getMatrixRowLength(row, rowLength);
802 if (err != 0) return(-1);
803
804 int* tmpIndices = indices;
805 double* tmpCoefs = coefs;
806
807 if (len < rowLength) {
808 tmpIndices = new int[rowLength];
809 tmpCoefs = new double[rowLength];
810 }
811
812 if (!blockMatrix_) {
813 A_ptr_->getRow(row, rowLength, tmpCoefs, tmpIndices);
814 }
815 else {
816 fei::console_out() << "Aztec_LinSysCore::getMatrixRow: not implemented." << FEI_ENDL;
817 return(-1);
818 }
819
820 if (len < rowLength) {
821 for(int i=0; i<len; i++) {
822 indices[i] = tmpIndices[i];
823 coefs[i] = tmpCoefs[i];
824 }
825 delete [] tmpIndices;
826 delete [] tmpCoefs;
827 }
828
829 return(0);
830}
831
832//==============================================================================
833int Aztec_LinSysCore::sumIntoBlockRow(int numBlkRows, const int* blkRows,
834 int numBlkCols, const int* blkCols,
835 const double* const* values,
836 int numPtCols,
837 bool overwriteInsteadOfAccumulate)
838{
839 int i;
840 //first, let's figure out which of the incoming blkRows is the biggest --
841 //i.e., which contains the most point-rows.
842 int maxBlkSize = 0;
843
844 for(i=0; i<numBlkRows; i++) {
845 int thisSize = lookup_->getBlkEqnSize(blkRows[i]);
846
847 if (maxBlkSize < thisSize) maxBlkSize = thisSize;
848 }
849
850 //now we can allocate an array to hold the values for passing each block
851 //row into the block-matrix.
852 double* coefs = new double[numPtCols*maxBlkSize];
853
854 int rowOffset = 0;
855 for(i=0; i<numBlkRows; i++) {
856 //now, copy the values for this block row into the coefs array.
857 copyBlockRow(i, blkRows, numBlkCols, blkCols, &(values[rowOffset]), coefs);
858
859 //now shove it into the matrix.
860 if (overwriteInsteadOfAccumulate) {
861 int err = blkA_ptr_->putBlockRow(blkRows[i], coefs,
862 (int*)blkCols, numBlkCols);
863 if (err != 0) {
864 fei::console_out() << thisProc_ << " DVBR::putBlockRow failed." << FEI_ENDL;
865 return(err);
866 }
867 }
868 else {
869 int err = blkA_ptr_->sumIntoBlockRow(blkRows[i], coefs,
870 (int*)blkCols, numBlkCols);
871 if (err != 0) {
872 fei::console_out() << thisProc_ << " DVBR::sumIntoBlockRow failed." << FEI_ENDL;
873 return(err);
874 }
875 }
876
877 rowOffset += lookup_->getBlkEqnSize(blkRows[i]);
878 }
879
880 delete [] coefs;
881 return(0);
882}
883
884//==============================================================================
885int Aztec_LinSysCore::copyBlockRow(int i, const int* blkRows,
886 int numBlkCols, const int* blkCols,
887 const double* const* values,
888 double* coefs){
889
890 int rowSize = 0, colSize = 0;
891 //coefs needs to contain the entries for each block together, and those
892 //entries need to be in column-major order
893 int colOffset = 0;
894 int coefOffset = 0;
895 for(int b=0; b<numBlkCols; b++) {
896
897 int err = blkA_ptr_->getBlockSize(blkRows[i], blkCols[b],
898 rowSize, colSize);
899 if (err) {
900 fei::console_out() << "Aztec_LSC:copyBlockRow: ERROR in getBlockSize" << FEI_ENDL;
901 return(-1);
902 }
903
904 for(int j=colOffset; j<colOffset+colSize; j++) {
905 for(int r=0; r<rowSize; r++) {
906 coefs[coefOffset++] = values[r][j];
907 }
908 }
909 colOffset += colSize;
910 }
911 return(0);
912}
913
914//==============================================================================
915int Aztec_LinSysCore::getBlockSize(int blkInd) {
916
917 int localBlkRow = blkInd - localBlkOffset_;
918 if (localBlkRow >= 0 && localBlkRow < numLocalEqnBlks_) {
919 //if this is a local blkInd, we have its size in an array.
920 return(localBlockSizes_[localBlkRow]);
921 }
922 else {
923 //else it ain't local, so we'll get its size from the matrix because
924 //we know that the matrix obtained it when it allocated itself.
925
926 int numRemoteBlocks = blkA_ptr_->getNumRemoteBlocks();
927 int* remoteInds = blkA_ptr_->getRemoteBlockIndices();
928 int* remoteSizes = blkA_ptr_->getRemoteBlockSizes();
929
930 int ins;
931 int index = fei::binarySearch(blkInd, remoteInds, numRemoteBlocks, ins);
932 if (index < 0)
933 messageAbort("getBlockSize: can't find blkInd.");
934
935 return(remoteSizes[index]);
936 }
937}
938
939//==============================================================================
940int Aztec_LinSysCore::sumIntoPointRow(int numPtRows, const int* ptRows,
941 int numPtCols, const int* ptColIndices,
942 const double* const* values,
943 bool overwriteInsteadOfAccumulate)
944{
945 int i, j;
946 if (debugOutput_) {
947 fprintf(debugFile_,"sumIntoPointRow, %d rows\n", numPtRows);
948 for(i=0; i<numPtRows; ++i) {
949 for(j=0; j<numPtCols; ++j) {
950 fprintf(debugFile_," sipr row %d, col %d, value: %e\n", ptRows[i],
951 ptColIndices[j], values[i][j]);
952 }
953 }
954 fflush(debugFile_);
955 }
956
957 if (!blockMatrix_) {
958 if (overwriteInsteadOfAccumulate) {
959 for(i=0; i<numPtRows; ++i) {
960 CHK_ERR( A_ptr_->putRow(ptRows[i], numPtCols, values[i], ptColIndices) );
961 }
962 }
963 else {
964 int err = A_ptr_->sumIntoRow(numPtRows, ptRows, numPtCols, ptColIndices,
965 values);
966 if (err != 0) {
967 FEI_OSTRINGSTREAM osstr;
968 osstr << "Aztec_LinSysCore::sumIntoPointRow ERROR calling A_ptr->sumIntoRow";
969 throw std::runtime_error(osstr.str());
970 }
971 }
972
973 return(0);
974 }
975
976 if (!haveLookup_) {
977 messageAbort("sumIntoPointRow: need lookup object, don't have it.");
978 }
979
980 int* blkIntData = new int[numPtRows*2 + numPtCols*2];
981 int* blkRows = blkIntData;
982 int* blkRowOffsets = blkIntData+numPtRows;
983 int* blkCols = blkIntData+2*numPtRows;
984 int* blkColOffsets = blkIntData+2*numPtRows+numPtCols;;
985
986 //now fill the blkRows and blkCols lists.
987
988 for(i=0; i<numPtRows; i++) {
989 blkRows[i] = lookup_->ptEqnToBlkEqn(ptRows[i]);
990 blkRowOffsets[i] = lookup_->getOffsetIntoBlkEqn(blkRows[i], ptRows[i]);
991 }
992
993 for(i=0; i<numPtCols; i++) {
994 blkCols[i] = lookup_->ptEqnToBlkEqn(ptColIndices[i]);
995 if (blkCols[i] < 0) {
996 fei::console_out() << thisProc_ << " lookup ptEqnToBlkEqn("<<ptColIndices[i]<<"): "
997 << blkCols[i] << FEI_ENDL;
998 messageAbort("negative blk-col");
999 }
1000
1001 blkColOffsets[i] = lookup_->getOffsetIntoBlkEqn(blkCols[i], ptColIndices[i]);
1002 }
1003
1004 for(i=0; i<numPtRows; i++) {
1005 for(j=0; j<numPtCols; j++) {
1006 sumPointIntoBlockRow(blkRows[i], blkRowOffsets[i],
1007 blkCols[j], blkColOffsets[j], values[i][j]);
1008 }
1009 }
1010
1011 delete [] blkIntData;
1012 return(0);
1013}
1014
1015//==============================================================================
1016int Aztec_LinSysCore::sumPointIntoBlockRow(int blkRow, int rowOffset,
1017 int blkCol, int colOffset,
1018 double value)
1019{
1020 int rowSize = lookup_->getBlkEqnSize(blkRow);
1021 int colSize = lookup_->getBlkEqnSize(blkCol);
1022
1023 int len = rowSize*colSize;
1024 if (len <= 0) {
1025 fei::console_out() << thisProc_ << ", ALSC::spibr: blkRow: " << blkRow << ", blkCol: " << blkCol << ", rowSize: " << rowSize << ", colSize: " << colSize << FEI_ENDL;
1026 messageAbort("sumPointIntoBlockRow: len <= 0");
1027 }
1028
1029 double* val = new double[rowSize*colSize];
1030
1031 for(int i=0; i<len; i++) val[i] = 0.0;
1032
1033 val[colOffset*rowSize + rowOffset] = value;
1034
1035 int err = blkA_ptr_->sumIntoBlockRow(blkRow, val, &blkCol, 1);
1036 if (err != 0) {
1037 fei::console_out() << thisProc_ << " DVBR::sumIntoBlockRow failed" << FEI_ENDL;
1038 return(err);
1039 }
1040
1041 delete [] val;
1042 return(0);
1043}
1044
1045//==============================================================================
1046int Aztec_LinSysCore::sumIntoRHSVector(int num,
1047 const double* values,
1048 const int* indices)
1049{
1050 //
1051 //This function scatters (accumulates) values into the linear-system's
1052 //currently selected RHS vector.
1053 //
1054 // num is how many values are being passed,
1055 // indices holds the global 'row-numbers' into which the values go,
1056 // and values holds the actual coefficients to be scattered.
1057 //
1058 rhsLoaded_ = false;
1059 double* cur_b = tmp_b_[currentRHS_];
1060
1061 if (debugOutput_) {
1062 for(int i=0; i<num; ++i) {
1063 fprintf(debugFile_,"sumIntoRHS %d, %e\n", indices[i], values[i]);
1064 fflush(debugFile_);
1065 }
1066 }
1067
1068 for(int i=0; i<num; i++){
1069 int localRow = indices[i] - localOffset_;
1070 if (localRow < 0 || localRow > numLocalEqns_) continue;
1071
1072 cur_b[localRow] += values[i];
1073 }
1074 return(0);
1075}
1076
1077//==============================================================================
1078int Aztec_LinSysCore::putIntoRHSVector(int num, const double* values,
1079 const int* indices)
1080{
1081//
1082//This function scatters (puts) values into the linear-system's
1083//currently selected RHS vector.
1084//
1085// num is how many values are being passed,
1086// indices holds the global 'row-numbers' into which the values go,
1087// and values holds the actual coefficients to be scattered.
1088//
1089 rhsLoaded_ = false;
1090
1091 for(int i=0; i<num; i++){
1092 int localRow = indices[i] - localOffset_;
1093 if (localRow < 0 || localRow > numLocalEqns_) continue;
1094
1095 if (debugOutput_) {
1096 fprintf(debugFile_,"putIntoRHS %d, %e\n", indices[i], values[i]);
1097 fflush(debugFile_);
1098 }
1099
1100 tmp_b_[currentRHS_][localRow] = values[i];
1101 }
1102 return(0);
1103}
1104
1105//==============================================================================
1106int Aztec_LinSysCore::getFromRHSVector(int num, double* values,
1107 const int* indices)
1108{
1109 //
1110 //This function retrieves values from the linear-system's
1111 //currently selected RHS vector.
1112 //
1113 // num is how many values are being requested,
1114 // indices holds the global 'row-numbers' for which values are requested,
1115 // and values holds the actual coefficients to be returned.
1116 //
1117 //if non-local indices are supplied, the corresponding positions in the values
1118 //array are not referenced.
1119 //
1120
1121 for(int i=0; i<num; i++){
1122 int localRow = indices[i] - localOffset_;
1123 if (localRow < 0 || localRow > numLocalEqns_) continue;
1124
1125 values[i] = tmp_b_[currentRHS_][localRow];
1126 }
1127 return(0);
1128}
1129
1130//==============================================================================
1131int Aztec_LinSysCore::matrixLoadComplete() {
1132
1133 if (debugOutput_) {
1134 fprintf(debugFile_,"matrixLoadComplete\n");
1135 fflush(debugFile_);
1136 }
1137
1138 if (matrixLoaded_ && rhsLoaded_) return(0);
1139
1141 int* data_org = NULL;
1142
1143 if (blockMatrix_) {
1144 tmpMap = blkMap_;
1145 }
1146 else {
1147 tmpMap = map_;
1148 }
1149
1150 if (!matrixLoaded_) {
1151 if (blockMatrix_) {
1152 if (!blkA_ptr_->isLoaded()) blkA_ptr_->loadComplete();
1153 }
1154 else {
1155 if (!A_ptr_->isFilled()) {
1156 A_ptr_->fillComplete();
1157 }
1158 }
1159
1160 matrixLoaded_ = true;
1161
1162 needNewPreconditioner_ = true;
1163 }
1164
1165 if (blockMatrix_) data_org = blkA_ptr_->getData_org();
1166 else data_org = A_ptr_->getAZ_MATRIX_PTR()->data_org;
1167
1168 Aztec_LSVector* tmp = NULL;
1169
1170 //if x_ is not null, then it probably has a previous solution in it that
1171 //we might as well keep for the next initial guess unless the user has
1172 //specifically set an initial guess.
1173 if (x_ != NULL) {
1174 tmp = new Aztec_LSVector(*x_);
1175 *tmp = *x_;
1176 }
1177
1178 //if x_ hasn't been allocated yet, we better do that now.
1179 if (x_ == NULL) x_ = new Aztec_LSVector(tmpMap, data_org);
1180 if (bc_ == NULL) bc_ = new Aztec_LSVector(tmpMap, data_org);
1181
1182 //if we did save off a copy of x_ above, let's put it back now.
1183 if (tmp != NULL) {
1184 *x_ = *tmp;
1185 delete tmp;
1186 }
1187 //if we've put boundary-condition data in tmp_bc_ then we better move it
1188 //into bc_ now.
1189 if (tmp_bc_ != NULL) {
1190 for(int j=0; j<numEssBCs_; j++) {
1191 int index = essBCindices_[j];
1192 (*bc_)[index] = tmp_bc_[index-localOffset_];
1193 }
1194 }
1195
1196 if (tmp_x_touched_) {
1197 //and now we can fill x_ with the stuff we've been holding in
1198 //tmp_x_, if tmp_x_ has been touched... i.e., if the user loaded an
1199 //initial guess into it.
1200 for(int i=0; i<numLocalEqns_; i++){
1201 (*x_)[i+localOffset_] = tmp_x_[i];
1202 }
1203 }
1204
1205 //now we're going to get the AZ_MATRIX ptr out of the AztecDMSR matrix
1206 //wrapper class.
1207
1208 if (blockMatrix_) azA_ = blkA_ptr_->getAZ_MATRIX_Ptr();
1209 else azA_ = A_ptr_->getAZ_MATRIX_PTR();
1210
1211 if (rhsLoaded_) return(0);
1212
1213 if (b_ != NULL) {
1214 for(int i=0; i<numRHSs_; i++) delete b_[i];
1215 delete [] b_;
1216 }
1217
1218 b_ = new Aztec_LSVector*[numRHSs_];
1219 for(int j=0; j<numRHSs_; j++) {
1220 b_[j] = new Aztec_LSVector(tmpMap, data_org);
1221 if (b_[j] == NULL) return(-1);
1222
1223 //now fill b_[j] with the stuff we've been holding in tmp_b_[j].
1224 for(int i=0; i<numLocalEqns_; i++){
1225 (*(b_[j]))[i+localOffset_] = tmp_b_[j][i];
1226 }
1227 }
1228
1229 b_ptr_ = b_[currentRHS_];
1230 vectorsAllocated_ = true;
1231
1232 if (!bcsLoaded_) modifyRHSforBCs();
1233 bcsLoaded_ = false;
1234
1235 rhsLoaded_ = true;
1236
1237 if (debugOutput_) {
1238 fprintf(debugFile_,"leaving matrixLoadComplete\n");
1239 fflush(debugFile_);
1240 }
1241 return(0);
1242}
1243
1244//==============================================================================
1245int Aztec_LinSysCore::getBlkEqnsAndOffsets(int* ptEqns, int* blkEqns,
1246 int* blkOffsets, int numEqns)
1247{
1248 if (!haveLookup_) messageAbort("getBlkEqnsAndOffsets: don't have Lookup");
1249
1250 for(int i=0; i<numEqns; i++) {
1251 blkEqns[i] = lookup_->ptEqnToBlkEqn(ptEqns[i]);
1252 if (blkEqns[i] < 0) {
1253 fei::console_out() << thisProc_ << "ptEqn: " << ptEqns[i] << ", blkEqn: " << blkEqns[i]
1254 << FEI_ENDL;
1255 messageAbort("getBlkEqnsAndOffsets: blk-eqn lookup failed.");
1256 }
1257
1258 blkOffsets[i] = lookup_->getOffsetIntoBlkEqn(blkEqns[i], ptEqns[i]);
1259 if (blkOffsets[i] < 0) {
1260 messageAbort("getBlkEqnsAndOffsets: blk-offset lookup failed.");
1261 }
1262 }
1263 return(0);
1264}
1265
1266//==============================================================================
1267int Aztec_LinSysCore::enforceEssentialBC(int* globalEqn,
1268 double* alpha,
1269 double* gamma, int len)
1270{
1271 //
1272 //This function must enforce an essential boundary condition on each local
1273 //equation in 'globalEqn'. This means that the following modification
1274 //should be made to A and b, for each globalEqn:
1275 //
1276 //for(each locally-owned equation i in the globalEqn array){
1277 // zero row i and put 1.0 on the diagonal
1278 // for(each row j in column i) {
1279 // if (i==j) b[i] = gamma/alpha;
1280 // else b[j] -= (gamma/alpha)*A[j,i];
1281 // A[j,i] = 0.0;
1282 // }
1283 //}
1284 //
1285
1286 if (len == 0) return(0);
1287
1288 std::vector<int> bcEqns; bcEqns.reserve(len);
1289 std::vector<int> indirect; indirect.reserve(len);
1290 for(int i=0; i<len; ++i) {
1291 bcEqns.push_back(globalEqn[i]);
1292 indirect.push_back(i);
1293 }
1294
1295 fei::insertion_sort_with_companions<int>(len, &bcEqns[0], &indirect[0]);
1296
1297 if (debugOutput_) {
1298 fprintf(debugFile_,"numEssBCs: %d\n", len);
1299 for(int ii=0; ii<len; ++ii) {
1300 fprintf(debugFile_, " EssBC eqn %d, alpha %e gamma %e\n",
1301 bcEqns[ii], alpha[indirect[ii]], gamma[indirect[ii]]);
1302 }
1303 fflush(debugFile_);
1304 }
1305
1306 bcsLoaded_ = true;
1307
1308 int localEnd = localOffset_ + numLocalEqns_ - 1;
1309
1310 if (debugOutput_) {
1311 fprintf(debugFile_,"localOffset_: %d, localEnd: %d\n", localOffset_, localEnd);
1312 fflush(debugFile_);
1313 }
1314
1315 {
1316 int* newBCindices = new int[numEssBCs_+len];
1317 int ii;
1318 for(ii=0; ii<numEssBCs_; ii++) newBCindices[ii] = essBCindices_[ii];
1319 int offset = 0;
1320 for(ii=0; ii<len; ii++) {
1321 if ((localOffset_ <= globalEqn[ii]) && (globalEqn[ii] <= localEnd)){
1322 newBCindices[numEssBCs_+offset++] = globalEqn[ii];
1323 }
1324 }
1325
1326 if (offset > 0) {
1327 delete [] essBCindices_;
1328 essBCindices_ = newBCindices;
1329 numEssBCs_ += offset;
1330 }
1331 else {
1332 delete [] newBCindices;
1333 }
1334 }
1335
1336 if (blockMatrix_) {
1337 int* blkIntData = new int[len*2];
1338 int* blkEqns = blkIntData;
1339 int* blkOffsets = blkIntData+len;
1340
1341 getBlkEqnsAndOffsets(globalEqn, blkEqns, blkOffsets, len);
1342
1343 enforceBlkEssentialBC(blkEqns, blkOffsets, alpha, gamma, len);
1344
1345 delete [] blkIntData;
1346
1347 return(0);
1348 }
1349
1350 for(int i=0; i<len; i++) {
1351
1352 int globalEqn_i = bcEqns[i];
1353
1354 //if globalEqn[i] isn't local, then the processor that owns it
1355 //should be running this code too. Otherwise there's trouble...
1356
1357 if ((localOffset_ > globalEqn_i) || (globalEqn_i > localEnd)) continue;
1358
1359 //zero this row, except for the diagonal coefficient.
1360 A_ptr_->setDiagEntry(globalEqn_i, 1.0);
1361 int* offDiagIndices = NULL;
1362 double* offDiagCoefs = NULL;
1363 int offDiagLength = 0;
1364 A_ptr_->getOffDiagRowPointers(globalEqn_i, offDiagIndices,
1365 offDiagCoefs, offDiagLength);
1366
1367 for(int jjj=0; jjj<offDiagLength; jjj++) offDiagCoefs[jjj] = 0.0;
1368
1369 //also, make the rhs modification here.
1370 double bc_term = gamma[indirect[i]]/alpha[indirect[i]];
1371 double rhs_term = bc_term;
1372 if (explicitDirichletBCs_) rhs_term = 0.0;
1373
1374 if (rhsLoaded_) {
1375 (*b_ptr_)[globalEqn_i] = rhs_term;
1376 (*bc_)[globalEqn_i] = bc_term;
1377 }
1378 else {
1379 tmp_b_[currentRHS_][globalEqn_i -localOffset_] = rhs_term;
1380 tmp_bc_[globalEqn_i -localOffset_] = bc_term;
1381 }
1382 }
1383
1384 if (BCenforcement_no_column_mod_ == true) {
1385 return(0);
1386 }
1387
1388 for(int row=localOffset_; row<=localEnd; row++) {
1389
1390 int insertPoint = -1;
1391 int index = fei::binarySearch(row, &bcEqns[0], len, insertPoint);
1392 if (index >= 0) continue;
1393
1394 int* offDiagIndices2 = NULL;
1395 double* offDiagCoefs2 = NULL;
1396 int offDiagLength2 = 0;
1397 A_ptr_->getOffDiagRowPointers(row, offDiagIndices2,
1398 offDiagCoefs2, offDiagLength2);
1399
1400 //look through this row to find the non-zeros in positions that
1401 //correspond to eqns in bcEqns and make the appropriate modification.
1402
1403 for(int j=0; j<offDiagLength2; j++) {
1404
1405 int col_index = A_ptr_->getAztec_Map()->getTransformedEqn(offDiagIndices2[j]);
1406
1407 int idx = fei::binarySearch(col_index, &bcEqns[0], len, insertPoint);
1408 if (idx < 0) continue;
1409
1410 double bc_term = gamma[indirect[idx]]/alpha[indirect[idx]];
1411
1412 double value = offDiagCoefs2[j]*bc_term;
1413
1414 if (rhsLoaded_) {
1415 (*b_ptr_)[row] -= value;
1416 }
1417 else {
1418 tmp_b_[currentRHS_][row-localOffset_] -= value;
1419 }
1420
1421 if (debugOutput_) {
1422 fprintf(debugFile_,"BC mod, rhs %d -= %e\n", row, value);
1423 fprintf(debugFile_,"BC, set A(%d,%d)==%e, to 0.0\n",
1424 row, bcEqns[idx], offDiagCoefs2[j]);
1425 }
1426
1427 offDiagCoefs2[j] = 0.0;
1428
1429 }//for offDiagLength2
1430
1431 }//for localEnd
1432
1433 return(0);
1434}
1435
1436//==============================================================================
1437int Aztec_LinSysCore::enforceBlkEssentialBC(int* blkEqn,
1438 int* blkOffset,
1439 double* alpha,
1440 double* gamma,
1441 int len)
1442{
1443//
1444//This function must enforce an essential boundary condition on each local
1445//equation specified by the pair 'blkEqn' and 'blkOffset'. A point-equation
1446//resides within a block-equation. The blkOffset gives the distance from the
1447//beginning of the block-equation to the point-equation.
1448
1449//The following modification should be made to A and b, for each incoming
1450//point-equation:
1451//
1452//for(each local equation i){
1453// for(each column j in row i) {
1454// if (i==j) b[i] = gamma/alpha;
1455// else b[j] -= (gamma/alpha)*A[j,i];
1456// }
1457//}
1458//
1459//all of row and column 'i' in A should be zeroed,
1460//except for 1.0 on the diagonal.
1461//
1462 int val_length = 0;
1463 double* val = NULL;
1464 int colInd_length = 0;
1465 int* blkCols = NULL;
1466 int rowNNZs = 0, numCols = 0, err = 0;
1467
1468 double* val2 = NULL;
1469 int val2Len = val_length;
1470 int* cols2 = NULL;
1471 int cols2Len = colInd_length;
1472
1473 int localEnd = localBlkOffset_ + numLocalEqnBlks_ - 1;
1474
1475 for(int i=0; i<len; i++) {
1476
1477 //if blkEqn[i] isn't local, then the processor that owns it
1478 //should be running this code too. Otherwise there's trouble...
1479
1480 if ((localBlkOffset_ > blkEqn[i]) || (blkEqn[i] > localEnd)){
1481 continue;
1482 }
1483
1484 err = getBlockRow(blkEqn[i], val, val_length, blkCols, colInd_length,
1485 numCols, rowNNZs);
1486 if (err) {
1487 fei::console_out() << "Aztec_LSC: ERROR in getBlockRow" << FEI_ENDL;
1488 return(-1);
1489 }
1490
1491 //now let's do the BC modification to this block-row.
1492
1493 err = blkRowEssBCMod(blkEqn[i], blkOffset[i], val, blkCols, numCols,
1494 rowNNZs, alpha[i], gamma[i]);
1495
1496 blkA_ptr_->putBlockRow(blkEqn[i], val, blkCols, numCols);
1497
1498 //now let's take advantage of the symmetry of element-contributions to
1499 //do the column-wise part of the BC modification. Since the structure of
1500 //the matrix arising from element contributions is symmetric, we know that
1501 //the column indices in row 'blkEqn' correspond to rows which have entries
1502 //in column 'blkEqn'. So we can modify the column in those rows rather
1503 //than searching all rows in the matrix looking for that column.
1504
1505 //so let's loop over the block-columns and do the column-wise essBC mod
1506 //to those rows.
1507
1508 for(int j=0; j<numCols; j++) {
1509
1510 int col_row = blkCols[j];
1511
1512 //if this column-index doesn't correspond to a local row, skip it.
1513 if ((localOffset_ > col_row) || (col_row > localEnd)) continue;
1514
1515 //if this column is the diagonal column, skip it.
1516 if (col_row == blkEqn[i]) continue;
1517
1518 int thisNNZ = 0;
1519 int thisNumBlks = 0;
1520 err = getBlockRow(col_row, val2, val2Len, cols2, cols2Len,
1521 thisNumBlks, thisNNZ);
1522
1523 err = blkColEssBCMod(col_row, blkEqn[i], blkOffset[i], val2, cols2,
1524 thisNumBlks, thisNNZ, alpha[i], gamma[i]);
1525
1526 blkA_ptr_->putBlockRow(col_row, val2, cols2, thisNumBlks);
1527
1528 }// end for(j<rowLength) loop
1529 }
1530
1531 delete [] val;
1532 delete [] blkCols;
1533 delete [] val2;
1534 delete [] cols2;
1535 return(0);
1536}
1537
1538//==============================================================================
1539int Aztec_LinSysCore::blkRowEssBCMod(int blkEqn, int blkOffset, double* val,
1540 int* blkCols, int numCols, int numPtNNZ,
1541 double alpha, double gamma)
1542{
1543//
1544//This function performs an essential boundary-condition modification for a
1545//single block-row.
1546//
1547 //we need to know which point-row this block-row corresponds to, so
1548 //we can do stuff to the rhs vector.
1549
1550 int pointRow = blockRowToPointRow(blkEqn);
1551
1552 int offset = 0;
1553
1554 for(int j=0; j<numCols; j++) {
1555
1556 int err, ptRows = 0, ptCols = 0;
1557 err = blkA_ptr_->getBlockSize(blkEqn, blkCols[j], ptRows, ptCols);
1558
1559 if (err) {
1560 fei::console_out() << "Aztec_LSC::blkRowEssBCMod: error in getBlockSize" << FEI_ENDL;
1561 return(1);
1562 }
1563
1564 if (blkCols[j] == blkEqn) {
1565 //this is the diagonal block, so we need to diagonalize the
1566 //'blkOffset'th point-row and point-column, leaving a 1.0 on the
1567 //diagonal.
1568 double bc_term = gamma/alpha;
1569 double rhs_term = bc_term;
1570 if (explicitDirichletBCs_) rhs_term = 0.0;
1571
1572 int thisOffset = offset;
1573
1574 for(int jj=0; jj<ptCols; jj++) {
1575 if (jj==blkOffset) {
1576 //if this is the point-column of interest, we move the
1577 //contents over to the rhs in the appropriate way, except
1578 //for the entry on the row-of-interest, which we just
1579 //set to 1.0.
1580
1581 for(int row=0; row<ptRows; row++) {
1582 if (row==blkOffset) {
1583 if (rhsLoaded_) {
1584 (*b_ptr_)[pointRow+row] = rhs_term;
1585 (*bc_)[pointRow+row] = bc_term;
1586 }
1587 else {
1588 tmp_b_[currentRHS_][pointRow+row-localOffset_] = rhs_term;
1589 tmp_bc_[pointRow+row-localOffset_] = bc_term;
1590 }
1591 val[thisOffset+row] = 1.0;
1592 }
1593 else {
1594 if (rhsLoaded_) {
1595 (*b_ptr_)[pointRow+row] -= val[thisOffset+row]
1596 * bc_term;
1597 }
1598 else {
1599 tmp_b_[currentRHS_][pointRow+row-localOffset_] -= val[thisOffset+row]*bc_term;
1600 }
1601 val[thisOffset+row] = 0.0;
1602 }
1603 }
1604 }
1605 else {
1606 val[thisOffset+blkOffset] = 0.0;
1607 }
1608
1609 thisOffset += ptRows;
1610 }
1611 }
1612 else {
1613 //this isn't the diagonal block, so we just want to zero the
1614 //whole 'blkOffset'th point-row in this block.
1615 int thisOffset = offset + blkOffset;
1616 for(int ii=0; ii<ptCols; ii++) {
1617 val[thisOffset] = 0.0;
1618 thisOffset += ptRows;
1619 }
1620 }
1621
1622 offset += ptRows*ptCols;
1623 }
1624
1625 return(0);
1626}
1627
1628//==============================================================================
1629int Aztec_LinSysCore::blkColEssBCMod(int blkRow, int blkEqn, int blkOffset,
1630 double* val, int* blkCols, int numCols,
1631 int numPtNNZ, double alpha, double gamma)
1632{
1633//
1634//This function does the column-wise modification for an Essential BC, to
1635//block column 'blkEqn', to the point-equation 'blkOffset' equations into the
1636//block, for the block-row contained in val,blkCols.
1637//
1638//NOTE: blkEqn is a 0-based equation number.
1639
1640 int thisPtRow = blockRowToPointRow(blkRow);
1641 int thisRowSize = 0, thisColSize = 0;
1642
1643 //look through the row to find the non-zero blk in position
1644 //blkEqn and make the appropriate modification.
1645 int err, offset = 0;
1646 for(int j=0; j<numCols; j++) {
1647 err = blkA_ptr_->getBlockSize(blkRow, blkCols[j],
1648 thisRowSize, thisColSize);
1649
1650 if (err) {
1651 fei::console_out() << "Aztec_LinSysCore::blkColEssBCMod: ERROR in getBlockSize" << FEI_ENDL;
1652 return(1);
1653 }
1654
1655 if (blkCols[j] == blkEqn) {
1656 double bc_term = gamma/alpha;
1657
1658 int thisOffset = offset + blkOffset*thisRowSize;
1659
1660 for(int row=0; row<thisRowSize; row++) {
1661 if (rhsLoaded_) {
1662 (*b_ptr_)[thisPtRow+row] -= val[thisOffset+row] * bc_term;
1663 }
1664 else {
1665 tmp_b_[currentRHS_][thisPtRow+row-localOffset_] -=
1666 val[thisOffset+row]*bc_term;
1667 }
1668 val[thisOffset+row] = 0.0;
1669 }
1670
1671 break;
1672 }
1673
1674 offset += thisRowSize*thisColSize;
1675 }// end for(j<numCols) loop
1676
1677 return(0);
1678}
1679
1680//==============================================================================
1681int Aztec_LinSysCore::blockRowToPointRow(int blkRow) {
1682//
1683//This function returns a (global 0-based) point-equation which corresponds to
1684//the first point-equation in block-row 'blkRow'.
1685//
1686 int localBlkRow = blkRow - localBlkOffset_;
1687
1688 if (localBlkRow < 0 || localBlkRow > numLocalEqnBlks_) return(-1);
1689
1690 int localPtRow = 0;
1691 for(int i=0; i<localBlkRow; i++) {
1692 localPtRow += localBlockSizes_[i];
1693 }
1694
1695 int pointRow = localPtRow + localOffset_;
1696 return(pointRow);
1697}
1698
1699//==============================================================================
1700int Aztec_LinSysCore::getBlockRow(int blkRow, double*& val, int& valLen,
1701 int*& blkColInds, int& blkColIndLen,
1702 int& numNzBlks, int& numNNZ) {
1703//
1704//This function gets a block row from the VBR matrix. The val array and the
1705//blkColInds array are of lengths valLen and blkColIndLen, respectively.
1706//On output, val and blkColInds are lengthened if necessary, with the new
1707//lengths updated in valLen and blkColIndLen. The actual number of nonzero
1708//blocks and nonzero point-entries are returned in numNzBlks and numNNZ.
1709//
1710 numNNZ = 0;
1711 int err = blkA_ptr_->getNumNonzerosPerRow(blkRow, numNNZ);
1712 if (err) {
1713 fei::console_out() << "Aztec_LSC::getBlockRow: ERROR in getNumNonzerosPerRow" << FEI_ENDL;
1714 return(1);
1715 }
1716
1717 numNzBlks = 0;
1718 err = blkA_ptr_->getNumBlocksPerRow(blkRow, numNzBlks);
1719 if (err) {
1720 fei::console_out() << "Aztec_LSC::getBlockRow: ERROR in getNumBlocksPerRow" << FEI_ENDL;
1721 return(1);
1722 }
1723
1724 if (numNNZ > valLen) {
1725 double* newVals = new double[numNNZ];
1726 delete [] val;
1727 val = newVals;
1728 valLen = numNNZ;
1729 }
1730
1731 if (numNzBlks > blkColIndLen) {
1732 int* newCols = new int[numNzBlks];
1733 delete [] blkColInds;
1734 blkColInds = newCols;
1735 blkColIndLen = numNzBlks;
1736 }
1737
1738 err = blkA_ptr_->getBlockRow(blkRow, val, blkColInds, numNzBlks);
1739
1740 return(err);
1741}
1742
1743//==============================================================================
1744int Aztec_LinSysCore::enforceRemoteEssBCs(int numEqns, int* globalEqns,
1745 int** colIndices, int* colIndLen,
1746 double** coefs)
1747{
1748 bcsLoaded_ = true;
1749
1750 //writeA("preRemBC");
1751
1752 if (debugOutput_) {
1753 for(int i=0; i<numEqns; ++i) {
1754 fprintf(debugFile_,"remBC row %d, (cols,coefs): ", globalEqns[i]);
1755 for(int j=0; j<colIndLen[i]; ++j) {
1756 fprintf(debugFile_, "(%d,%e) ",colIndices[i][j], coefs[i][j]);
1757 }
1758 fprintf(debugFile_,"\n");
1759 }
1760 fflush(debugFile_);
1761 }
1762
1763 if (blockMatrix_) {
1764 int* blkEqns = new int[numEqns];
1765 int* blkOffsets = new int[numEqns];
1766
1767 getBlkEqnsAndOffsets(globalEqns, blkEqns, blkOffsets, numEqns);
1768
1769 int** blkCols = new int*[numEqns];
1770 int** blkColOffsets = new int*[numEqns];
1771
1772 for(int i=0; i<numEqns; i++) {
1773 blkCols[i] = new int[colIndLen[i]];
1774 blkColOffsets[i] = new int[colIndLen[i]];
1775
1776 getBlkEqnsAndOffsets(colIndices[i], blkCols[i], blkColOffsets[i],
1777 colIndLen[i]);
1778 }
1779
1780 enforceBlkRemoteEssBCs(numEqns, blkEqns, blkCols, blkColOffsets,
1781 colIndLen, coefs);
1782
1783 delete [] blkEqns;
1784 delete [] blkOffsets;
1785
1786 for(int j=0; j<numEqns; j++) {
1787 delete [] blkCols[j];
1788 delete [] blkColOffsets[j];
1789 }
1790 delete [] blkCols;
1791 delete [] blkColOffsets;
1792
1793 return(0);
1794 }
1795
1796 int localEnd = localOffset_ + numLocalEqns_ - 1;
1797
1798 for(int i=0; i<numEqns; i++) {
1799 int globalEqn_i = globalEqns[i];
1800
1801 if ((localOffset_ > globalEqn_i) || (globalEqn_i > localEnd)){
1802 continue;
1803 }
1804
1805 int rowLen = 0;
1806 int* AcolInds = NULL;
1807 double* Acoefs = NULL;
1808
1809 A_ptr_->getOffDiagRowPointers(globalEqn_i, AcolInds, Acoefs, rowLen);
1810
1811 for(int j=0; j<colIndLen[i]; j++) {
1812 for(int k=0; k<rowLen; k++) {
1813 if (A_ptr_->getAztec_Map()->getTransformedEqn(AcolInds[k]) == colIndices[i][j]) {
1814 double value = Acoefs[k]*coefs[i][j];
1815
1816 double old_rhs_val = 0.0;
1817 if (rhsLoaded_) {
1818 old_rhs_val = (*b_ptr_)[globalEqn_i];
1819 (*b_ptr_)[globalEqn_i] -= value;
1820 }
1821 else {
1822 old_rhs_val = tmp_b_[currentRHS_][globalEqn_i -localOffset_];
1823 tmp_b_[currentRHS_][globalEqn_i -localOffset_] -= value;
1824 }
1825
1826 if (debugOutput_) {
1827 fprintf(debugFile_,"remBC mod, rhs %d (%e) -= %e\n",
1828 globalEqn_i, old_rhs_val, value);
1829 fprintf(debugFile_,"remBC, set A(%d,%d)==%e, to 0.0\n",
1830 globalEqn_i, AcolInds[k], Acoefs[k]);
1831 }
1832
1833 Acoefs[k] = 0.0;
1834 }
1835 }
1836 }
1837
1838 }
1839
1840 return(0);
1841}
1842
1843//==============================================================================
1844int Aztec_LinSysCore::enforceBlkRemoteEssBCs(int numEqns, int* blkEqns,
1845 int** blkColInds, int** blkColOffsets,
1846 int* blkColLens, double** remEssBCCoefs) {
1847 int val_length = 0;
1848 double* val = NULL;
1849 int colInd_length = 0;
1850 int* blkCols = NULL;
1851 int rowNNZs = 0, numCols = 0, err = 0;
1852
1853 int localEnd = localBlkOffset_ + numLocalEqnBlks_ - 1;
1854
1855 for(int i=0; i<numEqns; i++) {
1856 if ((localBlkOffset_ > blkEqns[i]) || (blkEqns[i] > localEnd)){
1857 continue;
1858 }
1859
1860 err = getBlockRow(blkEqns[i], val, val_length, blkCols, colInd_length,
1861 numCols, rowNNZs);
1862 if (err) {
1863 fei::console_out() << "Aztec_LSC:enforceBlkRemoteEssBC ERROR in getBlockRow" << FEI_ENDL;
1864 return(-1);
1865 }
1866
1867 //we need to know which point-row this block-row corresponds to, so
1868 //we can do stuff to the rhs vector.
1869
1870 int pointRow = blockRowToPointRow(blkEqns[i]);
1871
1872 int offset = 0;
1873
1874 for(int j=0; j<numCols; j++) {
1875 int ptRows = 0, ptCols = 0;
1876 err = blkA_ptr_->getBlockSize(blkEqns[i], blkCols[j],
1877 ptRows, ptCols);
1878 if (err) {
1879 fei::console_out() << "Aztec_LSC::enforceBlkRemoteEssBCs: error in getBlockSize"
1880 << FEI_ENDL;
1881 return(-1);
1882 }
1883
1884 for(int k=0; k<blkColLens[i]; k++) {
1885 if (blkColInds[i][k] == blkCols[j]) {
1886 int thisOffset = offset + blkColOffsets[i][k] * ptRows;
1887 double rhsTerm = remEssBCCoefs[i][k];
1888
1889 double* bvec = &(tmp_b_[currentRHS_][pointRow-localOffset_]);
1890 double* bcvec = &(tmp_bc_[pointRow-localOffset_]);
1891 for(int row=0; row<ptRows; row++) {
1892 double& coef = val[thisOffset+row];
1893 bvec[row] -= coef*rhsTerm;
1894 bcvec[row] -= coef*rhsTerm;
1895 coef = 0.0;
1896 }
1897
1898 blkA_ptr_->putBlockRow(blkEqns[i], &val[offset],
1899 &blkCols[j], 1);
1900 }
1901 }
1902
1903 offset += ptRows*ptCols;
1904 }
1905 }
1906
1907 delete [] val;
1908 delete [] blkCols;
1909 return(0);
1910}
1911
1912//==============================================================================
1913int Aztec_LinSysCore::getMatrixPtr(Data& data)
1914{
1915 if (!matrixLoaded_) matrixLoadComplete();
1916
1917 if (blockMatrix_) {
1918 data.setTypeName("AztecDVBR_Matrix");
1919 data.setDataPtr((void*)blkA_ptr_);
1920 }
1921 else {
1922 data.setTypeName("AztecDMSR_Matrix");
1923 data.setDataPtr((void*)A_ptr_);
1924 }
1925 return(0);
1926}
1927
1928//==============================================================================
1929int Aztec_LinSysCore::copyInMatrix(double scalar, const Data& data) {
1930//
1931//Overwrites the current internal matrix with a scaled copy of the
1932//input argument.
1933//
1934 if (blockMatrix_) {
1935 if (strcmp("AztecDVBR_Matrix", data.getTypeName()))
1936 messageAbort("copyInMatrix: data type string not 'AztecDVBR_Matrix'.");
1937
1938 AztecDVBR_Matrix* source = (AztecDVBR_Matrix*)data.getDataPtr();
1939
1940 if (blkA_ != NULL) delete blkA_;
1941 blkA_ = new AztecDVBR_Matrix(*source);
1942 blkA_ptr_ = blkA_;
1943
1944 VBRmatPlusScaledMat(blkA_ptr_, scalar, source);
1945
1946 azA_ = blkA_ptr_->getAZ_MATRIX_Ptr();
1947 }
1948 else {
1949 if (strcmp("AztecDMSR_Matrix", data.getTypeName()))
1950 messageAbort("copyInMatrix: data type string not 'AztecDMSR_Matrix'.");
1951
1952 AztecDMSR_Matrix* source = (AztecDMSR_Matrix*)data.getDataPtr();
1953 A_ptr_->copyStructure(*source);
1954
1955 MSRmatPlusScaledMat(A_ptr_, scalar, source);
1956
1957 azA_ = A_ptr_->getAZ_MATRIX_PTR();
1958 }
1959
1960 needNewPreconditioner_ = true;
1961 return(0);
1962}
1963
1964//==============================================================================
1965int Aztec_LinSysCore::copyOutMatrix(double scalar, Data& data) {
1966//
1967//Passes out a scaled copy of the current internal matrix.
1968//
1969
1970 if (!matrixLoaded_) matrixLoadComplete();
1971
1972 if (blockMatrix_) {
1973 AztecDVBR_Matrix* outmat = new AztecDVBR_Matrix(*blkA_ptr_);
1974
1975 //the copy-constructor above took all structural info from blkA_ptr_,
1976 //but not the data coefficients.
1977
1978 VBRmatPlusScaledMat(outmat, scalar, blkA_ptr_);
1979
1980 data.setTypeName("AztecDVBR_Matrix");
1981 data.setDataPtr((void*)outmat);
1982 }
1983 else {
1984 AztecDMSR_Matrix* outmat = new AztecDMSR_Matrix(*A_ptr_);
1985
1986 outmat->scale(scalar);
1987
1988 data.setTypeName("AztecDMSR_Matrix");
1989 data.setDataPtr((void*)outmat);
1990 }
1991 return(0);
1992}
1993
1994//==============================================================================
1995int Aztec_LinSysCore::sumInMatrix(double scalar, const Data& data) {
1996
1997 if (!matrixLoaded_) matrixLoadComplete();
1998
1999 if (blockMatrix_) {
2000 if (strcmp("AztecDVBR_Matrix", data.getTypeName())) {
2001 fei::console_out() << "Aztec_LinSysCore::sumInMatrix ERROR, incoming type-string: "
2002 << data.getTypeName() << ", expected AztecDVBR_Matrix." << FEI_ENDL;
2003 messageAbort("Aborting.");
2004 }
2005 AztecDVBR_Matrix* source = (AztecDVBR_Matrix*)data.getDataPtr();
2006
2007 VBRmatPlusScaledMat(blkA_ptr_, scalar, source);
2008 }
2009 else {
2010 if (strcmp("AztecDMSR_Matrix", data.getTypeName()))
2011 messageAbort("sumInMatrix: data type string not 'AztecDMSR_Matrix'.");
2012
2013 AztecDMSR_Matrix* source = (AztecDMSR_Matrix*)data.getDataPtr();
2014
2015 MSRmatPlusScaledMat(A_ptr_, scalar, source);
2016 }
2017
2018 needNewPreconditioner_ = true;
2019 return(0);
2020}
2021
2022//==============================================================================
2023int Aztec_LinSysCore::getRHSVectorPtr(Data& data) {
2024
2025 if (!matrixLoaded_) matrixLoadComplete();
2026
2027 data.setTypeName("Aztec_LSVector");
2028 data.setDataPtr((void*)b_ptr_);
2029 return(0);
2030}
2031
2032//==============================================================================
2033int Aztec_LinSysCore::copyInRHSVector(double scalar, const Data& data) {
2034
2035 if (!rhsLoaded_) matrixLoadComplete();
2036
2037 if (strcmp("Aztec_LSVector", data.getTypeName()))
2038 messageAbort("copyInRHSVector: data's type string not 'Aztec_LSVector'.");
2039
2040 Aztec_LSVector* sourcevec = (Aztec_LSVector*)data.getDataPtr();
2041
2042 *b_ptr_ = *sourcevec;
2043
2044 if (scalar != 1.0) b_ptr_->scale(scalar);
2045 return(0);
2046}
2047
2048//==============================================================================
2049int Aztec_LinSysCore::copyOutRHSVector(double scalar, Data& data) {
2050
2051 if (!rhsLoaded_) matrixLoadComplete();
2052
2053 Aztec_LSVector* outvec = new Aztec_LSVector(*b_ptr_);
2054
2055 outvec->put(0.0);
2056
2057 outvec->addVec(scalar, *b_ptr_);
2058
2059 data.setTypeName("Aztec_LSVector");
2060 data.setDataPtr((void*)outvec);
2061 return(0);
2062}
2063
2064//==============================================================================
2065int Aztec_LinSysCore::sumInRHSVector(double scalar, const Data& data) {
2066
2067 if (!rhsLoaded_) matrixLoadComplete();
2068
2069 if (strcmp("Aztec_LSVector", data.getTypeName()))
2070 messageAbort("sumInRHSVector: data's type string not 'Aztec_LSVector'.");
2071
2072 Aztec_LSVector* source = (Aztec_LSVector*)data.getDataPtr();
2073
2074 b_ptr_->addVec(scalar, *source);
2075 return(0);
2076}
2077
2078//==============================================================================
2079int Aztec_LinSysCore::destroyMatrixData(Data& data) {
2080
2081 if (blockMatrix_) {
2082 if (strcmp("AztecDVBR_Matrix", data.getTypeName()))
2083 messageAbort("destroyMatrixData: data doesn't contain a AztecDVBR_Matrix.");
2084
2085 AztecDVBR_Matrix* mat = (AztecDVBR_Matrix*)data.getDataPtr();
2086 delete mat;
2087 }
2088 else {
2089 if (strcmp("AztecDMSR_Matrix", data.getTypeName()))
2090 messageAbort("destroyMatrixData: data doesn't contain a AztecDMSR_Matrix.");
2091
2092 AztecDMSR_Matrix* mat = (AztecDMSR_Matrix*)data.getDataPtr();
2093 delete mat;
2094 }
2095 return(0);
2096}
2097
2098//==============================================================================
2099int Aztec_LinSysCore::destroyVectorData(Data& data) {
2100
2101 if (strcmp("Aztec_LSVector", data.getTypeName()))
2102 messageAbort("destroyVectorData: data doesn't contain a Aztec_LSVector.");
2103
2104 Aztec_LSVector* vec = (Aztec_LSVector*)data.getDataPtr();
2105 delete vec;
2106 return(0);
2107}
2108
2109//==============================================================================
2110int Aztec_LinSysCore::setNumRHSVectors(int numRHSs, const int* rhsIDs) {
2111
2112 if (numRHSs < 0)
2113 messageAbort("setNumRHSVectors: numRHSs < 0.");
2114
2115 if (numRHSs == 0) return(0);
2116
2117 delete [] rhsIDs_;
2118 numRHSs_ = numRHSs;
2119 rhsIDs_ = new int[numRHSs_];
2120
2121 for(int i=0; i<numRHSs_; i++) rhsIDs_[i] = rhsIDs[i];
2122 return(0);
2123}
2124
2125//==============================================================================
2126int Aztec_LinSysCore::setRHSID(int rhsID) {
2127
2128 for(int i=0; i<numRHSs_; i++){
2129 if (rhsIDs_[i] == rhsID){
2130 currentRHS_ = i;
2131 if (rhsLoaded_) b_ptr_ = b_[currentRHS_];
2132 return(0);
2133 }
2134 }
2135
2136 messageAbort("setRHSID: rhsID not found.");
2137 return(-1);
2138}
2139
2140//==============================================================================
2141int Aztec_LinSysCore::putInitialGuess(const int* eqnNumbers,
2142 const double* values,
2143 int len) {
2144//
2145//This function scatters (puts) values into the linear-system's soln vector.
2146//
2147// num is how many values are being passed,
2148// indices holds the global 'row-numbers' into which the values go,
2149// and values holds the actual coefficients to be scattered.
2150//
2151
2152 int localEnd = localOffset_ + numLocalEqns_ -1;
2153 if (matrixLoaded_) {
2154 for(int i=0; i<len; i++){
2155 if ((localOffset_ > eqnNumbers[i]) || (eqnNumbers[i] > localEnd))
2156 continue;
2157
2158 (*x_)[eqnNumbers[i]] = values[i];
2159 }
2160 }
2161 else {
2162 tmp_x_touched_ = true;
2163 for(int i=0; i<len; i++){
2164 if ((localOffset_ > eqnNumbers[i]) || (eqnNumbers[i] > localEnd))
2165 continue;
2166 tmp_x_[eqnNumbers[i]-localOffset_] = values[i];
2167 }
2168 }
2169 return(0);
2170}
2171
2172//==============================================================================
2173int Aztec_LinSysCore::getSolution(double* answers, int len) {
2174//
2175//The caller must allocate the memory for 'answers',
2176//and len must be set to the right value -- i.e., len should equal
2177//numLocalEqns_.
2178//
2179 if (len != numLocalEqns_)
2180 messageAbort("getSolution: len != numLocalEqns_.");
2181
2182 for(int i=0; i<numLocalEqns_; i++) {
2183 answers[i] = (*x_)[localOffset_ + i];
2184 }
2185 return(0);
2186}
2187
2188//==============================================================================
2189int Aztec_LinSysCore::getSolnEntry(int eqnNumber, double& answer) {
2190//
2191//This function returns a single solution entry, the coefficient for
2192//equation number eqnNumber.
2193//
2194 int localEnd = localOffset_ + numLocalEqns_ -1;
2195 if ((localOffset_ > eqnNumber) || (eqnNumber > localEnd))
2196 messageAbort("getSolnEntry: eqnNumber out of range.");
2197
2198 answer = (*x_)[eqnNumber];
2199 return(0);
2200}
2201
2202//==============================================================================
2203int Aztec_LinSysCore::formResidual(double* values, int len)
2204{
2205 if (len != numLocalEqns_) {
2206 messageAbort("formResidual: len != numLocalEqns_.");
2207 }
2208
2209 if (!matrixLoaded_ || !rhsLoaded_) {
2210 matrixLoadComplete();
2211 }
2212
2213 //after the solve x_ and b_ptr_ were 'inv_order'd back into user-ordering.
2214 //Before we can do this residual calculation, we have to put them into Aztec
2215 //ordering again.
2216
2217 int* update_index = NULL;
2218 if (blockMatrix_) {
2219 update_index = blkA_ptr_->getUpdate_index();
2220 }
2221 else {
2222 update_index = A_ptr_->getAztec_Map()->update_index;
2223 }
2224
2225 AZ_reorder_vec((double*)(x_->startPointer()), azA_->data_org, update_index,
2226 azA_->rpntr);
2227
2228 AZ_reorder_vec((double*)(b_ptr_->startPointer()), azA_->data_org,
2229 update_index, azA_->rpntr);
2230
2231 Aztec_LSVector* r = new Aztec_LSVector(*x_);
2232
2233 if (blockMatrix_) blkA_ptr_->matvec(*x_, *r); //form r = A*x
2234 else A_ptr_->matvec(*x_, *r);
2235
2236 r->addVec(-1.0, *b_ptr_); //form r = A*x - b
2237
2238 r->scale(-1.0); //form r = b - A*x (not really necessary, but...)
2239
2240 //now let's get the residual r into user ordering...
2241
2242 Aztec_LSVector* rtmp = new Aztec_LSVector(*x_);
2243
2244 AZ_invorder_vec((double*)(r->startPointer()), azA_->data_org, update_index,
2245 azA_->rpntr, (double*)rtmp->startPointer());
2246
2247 *r = *rtmp;
2248
2249 const double* rptr = r->startPointer();
2250
2251 for(int i=0; i<numLocalEqns_; i++) {
2252 values[i] = rptr[i];
2253 }
2254
2255 //finally, let's put x_ and b_ptr_ back into user ordering before we exit...
2256 AZ_invorder_vec((double*)(x_->startPointer()), azA_->data_org, update_index,
2257 azA_->rpntr, (double*)rtmp->startPointer());
2258
2259 *x_ = *rtmp;
2260
2261 AZ_invorder_vec((double*)(b_ptr_->startPointer()), azA_->data_org,
2262 update_index, azA_->rpntr, (double*)rtmp->startPointer());
2263
2264 *b_ptr_ = *rtmp;
2265
2266 delete rtmp;
2267 delete r;
2268 return(0);
2269}
2270
2271//==============================================================================
2272int Aztec_LinSysCore::selectSolver(const char* name) {
2273//
2274// This function takes a string naming the solver and sets the solver choice
2275// in the options array accordingly.
2276//
2277 char msg[64];
2278
2279 if (!strcmp(name, "AZ_gmres")) {
2280 aztec_options_[AZ_solver] = AZ_gmres;
2281 sprintf(msg, "AZ_gmres solver.");
2282 }
2283 else if (!strcmp(name, "AZ_cg")) {
2284 aztec_options_[AZ_solver] = AZ_cg;
2285 sprintf(msg, "AZ_cg solver.");
2286 }
2287 else if (!strcmp(name, "AZ_bicgstab")) {
2288 aztec_options_[AZ_solver] = AZ_bicgstab;
2289 sprintf(msg, "AZ_bicgstab solver.");
2290 }
2291 else if (!strcmp(name, "AZ_cgs")) {
2292 aztec_options_[AZ_solver] = AZ_cgs;
2293 sprintf(msg, "AZ_cgs solver.");
2294 }
2295 else if (!strcmp(name, "AZ_tfqmr")) {
2296 aztec_options_[AZ_solver] = AZ_tfqmr;
2297 sprintf(msg, "AZ_tfqmr solver.");
2298 }
2299 else if (!strcmp(name, "AZ_lu")) {
2300 aztec_options_[AZ_solver] = AZ_lu;
2301 sprintf(msg, "AZ_lu solver.");
2302 }
2303 else {
2304 aztec_options_[AZ_solver] = AZ_gmres;
2305 sprintf(msg, "AZ_gmres default solver.");
2306 if (thisProc_ == 0) {
2307 FEI_COUT << "Aztec_LinSysCore: Warning: requested solver <" << name << "> not recognized, defaulting to AZ_gmres." << FEI_ENDL;
2308 }
2309 }
2310
2311 debugOutput(msg);
2312
2313 return(0);
2314}
2315
2316//==============================================================================
2317int Aztec_LinSysCore::selectPreconditioner(const char* name) {
2318
2319 char msg[128];
2320 sprintf(msg, "selectPreconditioner(%s)", name);
2321 debugOutput(msg);
2322
2323 if (!strcmp(name, "AZ_none")) {
2324 aztec_options_[AZ_precond] = AZ_none;
2325 sprintf(msg, " -- selected: AZ_none.");
2326 }
2327 else if (!strcmp(name, "AZ_Jacobi")) {
2328 aztec_options_[AZ_precond] = AZ_Jacobi;
2329 sprintf(msg, " -- selected: AZ_Jacobi.");
2330 }
2331 else if (!strcmp(name, "AZ_Neumann")) {
2332 aztec_options_[AZ_precond] = AZ_Neumann;
2333 sprintf(msg, " -- selected: AZ_Neumann.");
2334 }
2335 else if (!strcmp(name, "AZ_ls")) {
2336 aztec_options_[AZ_precond] = AZ_ls;
2337 sprintf(msg, " -- selected: AZ_ls.");
2338 }
2339 else if (!strcmp(name, "AZ_sym_GS")) {
2340 aztec_options_[AZ_precond] = AZ_sym_GS;
2341 sprintf(msg, " -- selected: AZ_sym_GS.");
2342 }
2343 else if (!strcmp(name, "AZ_dom_decomp")) {
2344 aztec_options_[AZ_precond] = AZ_dom_decomp;
2345 sprintf(msg, " -- selected: AZ_dom_decomp.");
2346 }
2347//#ifndef NOT_USING_ML
2348#ifdef USING_ML
2349 else if (!strcmp(name, "ML_Vanek")) {
2350 ML_Vanek_ = true;
2351 aztec_options_[AZ_precond] = AZ_user_precond;
2352 sprintf(msg, " -- selected: AZ_user_precond.");
2353 }
2354#endif
2355 else {
2356 aztec_options_[AZ_precond] = AZ_none;
2357 sprintf(msg," -- selected: Default, AZ_none.\n");
2358 if (thisProc_ == 0) {
2359 FEI_COUT << "Aztec_LinSysCore: Warning: requested preconditioner <" << name << "> not recognized, defaulting to AZ_none." << FEI_ENDL;
2360 }
2361 }
2362
2363 debugOutput(msg);
2364
2365 return(0);
2366}
2367
2368//==============================================================================
2369void Aztec_LinSysCore::setSubdomainSolve(const char* name) {
2370
2371 char msg[128];
2372 sprintf(msg, "setSubdomainSolve(%s)", name);
2373 debugOutput(msg);
2374
2375 if (!strcmp(name, "AZ_lu")) {
2376 aztec_options_[AZ_subdomain_solve] = AZ_lu;
2377 sprintf(msg, " -- selected AZ_lu");
2378 }
2379 else if (!strcmp(name, "AZ_ilu")) {
2380 aztec_options_[AZ_subdomain_solve] = AZ_ilu;
2381 sprintf(msg, " -- selected AZ_ilu");
2382 }
2383 else if (!strcmp(name, "AZ_ilut")) {
2384 aztec_options_[AZ_subdomain_solve] = AZ_ilut;
2385 sprintf(msg, " -- selected AZ_ilut");
2386 }
2387 else if (!strcmp(name, "AZ_rilu")) {
2388 aztec_options_[AZ_subdomain_solve] = AZ_rilu;
2389 sprintf(msg, " -- selected AZ_rilu");
2390 }
2391 else if (!strcmp(name, "AZ_bilu")) {
2392 aztec_options_[AZ_subdomain_solve] = AZ_bilu;
2393 sprintf(msg, " -- selected AZ_bilu");
2394 }
2395 else if (!strcmp(name, "AZ_icc")) {
2396 aztec_options_[AZ_subdomain_solve] = AZ_icc;
2397 sprintf(msg, " -- selected AZ_icc");
2398 }
2399 else {
2400 if (thisProc_ == 0) {
2401 FEI_COUT << "Aztec_LinSysCore: Warning: requested subdomain-solve <" << name << "> not recognized." << FEI_ENDL;
2402 }
2403 }
2404
2405 debugOutput(msg);
2406}
2407
2408//==============================================================================
2409void Aztec_LinSysCore::setTypeOverlap(const char* name) {
2410
2411 char msg[128];
2412 sprintf(msg, "setTypeOverlap(%s)", name);
2413 debugOutput(msg);
2414
2415 if (!strcmp(name, "AZ_standard")) {
2416 aztec_options_[AZ_type_overlap] = AZ_standard;
2417 sprintf(msg, " -- selected AZ_standard");
2418 }
2419 else if (!strcmp(name, "AZ_ilu")) {
2420 aztec_options_[AZ_type_overlap] = AZ_symmetric;
2421 sprintf(msg, " -- selected AZ_symmetric");
2422 }
2423 else {
2424 if (thisProc_ == 0) {
2425 FEI_COUT << "Aztec_LinSysCore: Warning: requested type-overlap <" << name << "> not recognized." << FEI_ENDL;
2426 }
2427 }
2428
2429 debugOutput(msg);
2430}
2431
2432//==============================================================================
2433int Aztec_LinSysCore::writeSystem(const char* name)
2434{
2435 writeA(name);
2436 return(0);
2437}
2438
2439//==============================================================================
2440void Aztec_LinSysCore::recordUserParams()
2441{
2442 checkForParam("AZ_tol", numParams_, paramStrings_,
2443 aztec_params_[AZ_tol]);
2444
2445 checkForParam("AZ_drop", numParams_, paramStrings_,
2446 aztec_params_[AZ_drop]);
2447
2448 checkForParam("AZ_ilut_fill", numParams_, paramStrings_,
2449 aztec_params_[AZ_ilut_fill]);
2450
2451 checkForParam("AZ_omega", numParams_, paramStrings_,
2452 aztec_params_[AZ_omega]);
2453
2454 checkForParam("AZ_weights", numParams_, paramStrings_,
2455 aztec_params_[AZ_weights]);
2456}
2457
2458//==============================================================================
2459void Aztec_LinSysCore::recordUserOptions()
2460{
2461//
2462//Private function, to be called from launchSolver.
2463//
2464 const char* param = NULL;
2465
2466 param = snl_fei::getParamValue("AZ_solver", numParams_, paramStrings_);
2467 if (param != NULL){
2468 selectSolver(param);
2469 }
2470
2471 param = snl_fei::getParamValue("AZ_precond", numParams_, paramStrings_);
2472 if (param != NULL){
2473 selectPreconditioner(param);
2474 }
2475
2476 param = snl_fei::getParamValue("AZ_subdomain_solve",
2477 numParams_, paramStrings_);
2478 if (param != NULL){
2479 setSubdomainSolve(param);
2480 }
2481
2482 param = snl_fei::getParamValue("AZ_scaling", numParams_, paramStrings_);
2483 if (param != NULL){
2484 setScalingOption(param);
2485 }
2486
2487 param = snl_fei::getParamValue("AZ_conv", numParams_, paramStrings_);
2488 if (param != NULL){
2489 setConvTest(param);
2490 }
2491
2492 param = snl_fei::getParamValue("AZ_pre_calc",
2493 numParams_, paramStrings_);
2494 if (param != NULL){
2495 setPreCalc(param);
2496 }
2497
2498 param = snl_fei::getParamValue("AZ_overlap", numParams_, paramStrings_);
2499 if (param != NULL){
2500 setOverlap(param);
2501 }
2502
2503 param = snl_fei::getParamValue("AZ_type_overlap",
2504 numParams_, paramStrings_);
2505 if (param != NULL){
2506 setTypeOverlap(param);
2507 }
2508
2509 param = snl_fei::getParamValue("AZ_orthog", numParams_, paramStrings_);
2510 if (param != NULL){
2511 setOrthog(param);
2512 }
2513
2514 param = snl_fei::getParamValue("AZ_aux_vec", numParams_, paramStrings_);
2515 if (param != NULL){
2516 setAuxVec(param);
2517 }
2518
2519 param = snl_fei::getParamValue("AZ_output", numParams_, paramStrings_);
2520 if (param != NULL){
2521 setAZ_output(param);
2522 }
2523 else aztec_options_[AZ_output] = outputLevel_;
2524
2525 checkForOption("AZ_poly_ord", numParams_, paramStrings_,
2526 aztec_options_[AZ_poly_ord]);
2527
2528 checkForOption("AZ_kspace", numParams_, paramStrings_,
2529 aztec_options_[AZ_kspace]);
2530
2531 checkForOption("AZ_max_iter", numParams_, paramStrings_,
2532 aztec_options_[AZ_max_iter]);
2533
2534 checkForOption("AZ_reorder", numParams_, paramStrings_,
2535 aztec_options_[AZ_reorder]);
2536
2537 checkForOption("AZ_graph_fill", numParams_, paramStrings_,
2538 aztec_options_[AZ_graph_fill]);
2539
2540 checkForOption("AZ_keep_info", numParams_, paramStrings_,
2541 aztec_options_[AZ_keep_info]);
2542}
2543
2544//==============================================================================
2545int Aztec_LinSysCore::writeA(const char* name)
2546{
2547 if (name == NULL) {
2548 return(-1);
2549 }
2550
2551 FEI_OSTRINGSTREAM osstr;
2552
2553 if (debugPath_ != NULL) {
2554 osstr << debugPath_;
2555 }
2556 else osstr << ".";
2557
2558 osstr << "/A_" << name << ".mtx";
2559
2560 if (blockMatrix_) osstr << ".vbr";
2561
2562 std::string str = osstr.str();
2563 const char* matname = str.c_str();
2564
2565 if (blockMatrix_) {
2566 blkA_ptr_->writeToFile(matname);
2567 }
2568 else {
2569 A_ptr_->writeToFile(matname);
2570 }
2571
2572 return(0);
2573}
2574
2575//==============================================================================
2576int Aztec_LinSysCore::writeVec(Aztec_LSVector* v, const char* name)
2577{
2578 if (name == NULL || v == NULL) {
2579 return(-1);
2580 }
2581
2582 FEI_OSTRINGSTREAM osstr;
2583
2584 if (debugPath_ != NULL) {
2585 osstr << debugPath_;
2586 }
2587 else osstr << ".";
2588
2589 osstr << "/" << name << ".vec";
2590
2591 std::string str = osstr.str();
2592
2593 const char* vecname = str.c_str();
2594
2595 v->writeToFile(vecname);
2596
2597 return(0);
2598}
2599
2600//==============================================================================
2601int Aztec_LinSysCore::modifyRHSforBCs()
2602{
2603 if (explicitDirichletBCs_) {
2604 for(int j=0; j<numEssBCs_; j++) {
2605 int index = essBCindices_[j];
2606 (*b_ptr_)[index] = 0.0;
2607 }
2608 }
2609 else {
2610 for(int j=0; j<numEssBCs_; j++) {
2611 int index = essBCindices_[j];
2612 (*b_ptr_)[index] = tmp_bc_[index-localOffset_];
2613 }
2614 }
2615
2616 return(0);
2617}
2618
2619//==============================================================================
2620int Aztec_LinSysCore::explicitlySetDirichletBCs()
2621{
2622 for(int j=0; j<numEssBCs_; j++) {
2623 int index = essBCindices_[j];
2624 if (rhsLoaded_) {
2625 (*b_ptr_)[index] = (*bc_)[index];
2626 (*x_)[index] = (*bc_)[index];
2627 }
2628 else {
2629 (*b_ptr_)[index] = tmp_bc_[index-localOffset_];
2630 (*x_)[index] = tmp_bc_[index-localOffset_];
2631 }
2632 }
2633 return(0);
2634}
2635
2636//==============================================================================
2637int Aztec_LinSysCore::launchSolver(int& solveStatus, int& iterations) {
2638//
2639//This function does any last-second setup required for the
2640//linear solver, then goes ahead and launches the solver to get
2641//the solution vector.
2642//Also, if possible, the number of iterations that were performed
2643//is stored in the iterations_ variable.
2644//
2645
2646 unsigned counter = 0;
2647 std::map<std::string,unsigned>::iterator
2648 iter = named_solve_counter_.find(name_);
2649 if (iter == named_solve_counter_.end()) {
2650 fei::console_out() << "fei: Aztec_LinSysCore::LaunchSolver: internal error."
2651 << FEI_ENDL;
2652 }
2653 else {
2654 counter = iter->second++;
2655 }
2656
2657 if (debugOutput_ && outputLevel_ > 1) {
2658 FEI_OSTRINGSTREAM osstr;
2659 osstr << name_ << "_Aztec.np"<<numProcs_<<".slv"<<counter;
2660 std::string str = osstr.str();
2661
2662 writeA(str.c_str());
2663
2664 FEI_OSTRINGSTREAM x0_osstr;
2665 x0_osstr << "x0_" << str;
2666 std::string x0_str = x0_osstr.str();
2667
2668 writeVec(x_, x0_str.c_str());
2669
2670 FEI_OSTRINGSTREAM b_osstr;
2671 b_osstr << "b_" << str;
2672 std::string b_str = b_osstr.str();
2673
2674 writeVec(b_ptr_, b_str.c_str());
2675 }
2676
2677 if (needNewPreconditioner_) {
2678 if (precondCreated_) {
2679 AZ_precond_destroy(&azP_);
2680 }
2681
2682//#ifndef NOT_USING_ML
2683#ifdef USING_ML
2684 ML* ml = NULL;
2685
2686 if (ML_Vanek_) {
2687 int numFineSweeps = 2;
2688 int numCoarseSweeps = 2;
2689 double omega = 0.67;
2690
2691 initialize_ML(azA_, &azP_, numLevels_,
2692 numFineSweeps, numCoarseSweeps, omega,
2693 map_->getProcConfig(), &ml);
2694 }
2695 else {
2696 //set the preconditioning matrix Pmat to point to azA_ (Amat).
2697 azA_->data_org[AZ_name] = 0;
2698 azP_ = AZ_precond_create(azA_, AZ_precondition, NULL);
2699 }
2700#else
2701 azA_->data_org[AZ_name] = 0;
2702 azP_ = AZ_precond_create(azA_, AZ_precondition, NULL);
2703#endif
2704 needNewPreconditioner_ = false;
2705 aztec_options_[AZ_pre_calc] = AZ_calc;
2706 aztec_options_[AZ_conv] = AZ_rhs;
2707 aztec_options_[AZ_orthog] = AZ_modified;
2708
2709 recordUserParams();
2710 recordUserOptions();
2711
2712 aztec_options_[AZ_keep_info] = 1;
2713 }
2714 else {
2715 aztec_options_[AZ_pre_calc] = AZ_reuse;
2716 }
2717
2718 precondCreated_ = true;
2719
2720 int* proc_config = NULL;
2721 int* update_index = NULL;
2722 if (blockMatrix_) {
2723 proc_config = blkMap_->getProcConfig();
2724 update_index = blkA_ptr_->getUpdate_index();
2725 }
2726 else {
2727 proc_config = map_->getProcConfig();
2728 update_index = A_ptr_->getAztec_Map()->update_index;
2729 }
2730
2731 AZ_reorder_vec((double*)(x_->startPointer()), azA_->data_org, update_index,
2732 azA_->rpntr);
2733
2734 AZ_reorder_vec((double*)(b_ptr_->startPointer()), azA_->data_org,
2735 update_index, azA_->rpntr);
2736
2737 AZ_iterate((double*)(x_->startPointer()),
2738 (double*)(b_ptr_->startPointer()),
2739 aztec_options_, aztec_params_, aztec_status_,
2740 proc_config, azA_, azP_, azS_);
2741
2742 iterations = (int)aztec_status_[AZ_its];
2743
2744 solveStatus = (int)aztec_status_[AZ_why];
2745
2746 azlsc_solveCounter_++;
2747
2748 Aztec_LSVector* xtmp = new Aztec_LSVector(*x_);
2749
2750 //now we need to put x_ back into user-ordering for when we're asked to
2751 //hand out solution entries.
2752
2753 AZ_invorder_vec((double*)(x_->startPointer()), azA_->data_org, update_index,
2754 azA_->rpntr, (double*)xtmp->startPointer());
2755
2756 *x_ = *xtmp;
2757
2758 //let's put b back into user-ordering too...
2759 AZ_invorder_vec((double*)(b_ptr_->startPointer()), azA_->data_org,
2760 update_index, azA_->rpntr, (double*)xtmp->startPointer());
2761
2762 *b_ptr_ = *xtmp;
2763
2764 delete xtmp;
2765
2766 if (explicitDirichletBCs_) explicitlySetDirichletBCs();
2767 else {
2768 }
2769
2770 if (debugOutput_) {
2771 FEI_OSTRINGSTREAM osstr;
2772 osstr << name_ << "_Aztec.np"<<numProcs_<<".slv"<<counter;
2773 std::string str = osstr.str();
2774
2775 FEI_OSTRINGSTREAM x_osstr;
2776 x_osstr << "x_" << str;
2777 std::string x_str = x_osstr.str();
2778
2779 writeVec(x_, x_str.c_str());
2780 }
2781 for(int j=0; j<numEssBCs_; j++) {
2782 int index = essBCindices_[j];
2783 if (rhsLoaded_) {
2784 (*bc_)[index] = 0.0;
2785 }
2786 else {
2787 tmp_bc_[index-localOffset_] = 0.0;
2788 }
2789 }
2790 delete [] essBCindices_;
2791 essBCindices_ = NULL;
2792 numEssBCs_ = 0;
2793
2794 return(0);
2795}
2796
2797//==============================================================================
2798void Aztec_LinSysCore::setScalingOption(const char* param) {
2799
2800 char* msg = new char[128];
2801 for(int i=0; i<128; i++) msg[i] = '\0';
2802
2803 if (!strcmp(param, "AZ_none")) {
2804 aztec_options_[AZ_scaling] = AZ_none;
2805 sprintf(msg, "No scaling");
2806 }
2807 else if (!strcmp(param, "AZ_Jacobi")) {
2808 aztec_options_[AZ_scaling] = AZ_Jacobi;
2809 sprintf(msg, "AZ_Jacobi scaling");
2810 }
2811 else if (!strcmp(param, "AZ_BJacobi")) {
2812 aztec_options_[AZ_scaling] = AZ_BJacobi;
2813 sprintf(msg, "AZ_BJacobi scaling");
2814 }
2815 else if (!strcmp(param, "AZ_row_sum")) {
2816 aztec_options_[AZ_scaling] = AZ_row_sum;
2817 sprintf(msg, "AZ_row_sum scaling");
2818 }
2819 else if (!strcmp(param, "AZ_sym_diag")) {
2820 aztec_options_[AZ_scaling] = AZ_sym_diag;
2821 sprintf(msg, "AZ_sym_diag scaling");
2822 }
2823 else if (!strcmp(param, "AZ_sym_row_sum")) {
2824 aztec_options_[AZ_scaling] = AZ_sym_row_sum;
2825 sprintf(msg, "AZ_sym_row_sum scaling");
2826 }
2827 else {
2828 if (thisProc_ == 0) {
2829 FEI_COUT << "Aztec_LinSysCore: Warning: requested scaling <" << param << "> not recognized." << FEI_ENDL;
2830 }
2831 }
2832
2833 debugOutput(msg);
2834
2835 delete [] msg;
2836
2837 return;
2838}
2839
2840//==============================================================================
2841void Aztec_LinSysCore::setConvTest(const char* param) {
2842
2843 char msg[64];
2844
2845 if (!strcmp(param, "AZ_r0")) {
2846 aztec_options_[AZ_conv] = AZ_r0;
2847 sprintf(msg, "AZ_conv AZ_r0");
2848 }
2849 else if (!strcmp(param, "AZ_rhs")) {
2850 aztec_options_[AZ_conv] = AZ_rhs;
2851 sprintf(msg, "AZ_conv AZ_rhs");
2852 }
2853 else if (!strcmp(param, "AZ_Anorm")) {
2854 aztec_options_[AZ_conv] = AZ_Anorm;
2855 sprintf(msg, "AZ_conv AZ_Anorm");
2856 }
2857 else if (!strcmp(param, "AZ_sol")) {
2858 aztec_options_[AZ_conv] = AZ_sol;
2859 sprintf(msg, "AZ_conv AZ_sol");
2860 }
2861 else if (!strcmp(param, "AZ_weighted")) {
2862 aztec_options_[AZ_conv] = AZ_weighted;
2863 sprintf(msg, "AZ_conv AZ_weighted");
2864 }
2865 else if (!strcmp(param, "AZ_noscaled")) {
2866 aztec_options_[AZ_conv] = AZ_noscaled;
2867 sprintf(msg, "AZ_conv AZ_noscaled");
2868 }
2869 else {
2870 if (thisProc_ == 0) {
2871 FEI_COUT << "Aztec_LinSysCore: Warning: requested convergence test <" << param << "> not recognized." << FEI_ENDL;
2872 }
2873 }
2874
2875 debugOutput(msg);
2876
2877 return;
2878}
2879
2880//==============================================================================
2881void Aztec_LinSysCore::setPreCalc(const char* param)
2882{
2883 char msg[64];
2884
2885 if (!strcmp(param, "AZ_calc")) {
2886 aztec_options_[AZ_pre_calc] = AZ_calc;
2887 sprintf(msg, "AZ_pre_calc AZ_calc");
2888 }
2889 else if (!strcmp(param, "AZ_recalc")) {
2890 aztec_options_[AZ_pre_calc] = AZ_recalc;
2891 sprintf(msg, "AZ_pre_calc AZ_recalc");
2892 }
2893 else if (!strcmp(param, "AZ_reuse")) {
2894 aztec_options_[AZ_pre_calc] = AZ_reuse;
2895 sprintf(msg, "AZ_pre_calc AZ_reuse");
2896 }
2897 else {
2898 if (thisProc_ == 0) {
2899 FEI_COUT << "Aztec_LinSysCore: Warning: requested pre_calc <" << param << "> not recognized." << FEI_ENDL;
2900 }
2901 }
2902
2903 debugOutput(msg);
2904
2905 return;
2906}
2907
2908//==============================================================================
2909void Aztec_LinSysCore::setOverlap(const char* param)
2910{
2911 char msg[64];
2912
2913 if (!strcmp(param, "AZ_none")) {
2914 aztec_options_[AZ_overlap] = AZ_none;
2915 sprintf(msg, "AZ_overlap AZ_none");
2916 }
2917 else if (!strcmp(param, "AZ_diag")) {
2918 aztec_options_[AZ_overlap] = AZ_diag;
2919 sprintf(msg, "AZ_overlap AZ_diag");
2920 }
2921 else if (!strcmp(param, "AZ_full")) {
2922 aztec_options_[AZ_overlap] = AZ_full;
2923 sprintf(msg, "AZ_overlap AZ_full");
2924 }
2925 else {
2926 checkForOption("AZ_overlap", numParams_, paramStrings_,
2927 aztec_options_[AZ_overlap]);
2928 }
2929
2930 debugOutput(msg);
2931
2932 return;
2933}
2934
2935//==============================================================================
2936void Aztec_LinSysCore::setOrthog(const char* param)
2937{
2938 char msg[64];
2939
2940 if (!strcmp(param, "AZ_classic")) {
2941 aztec_options_[AZ_orthog] = AZ_classic;
2942 sprintf(msg, "AZ_orthog AZ_classic");
2943 }
2944 else if (!strcmp(param, "AZ_modified")) {
2945 aztec_options_[AZ_orthog] = AZ_modified;
2946 sprintf(msg, "AZ_orthog AZ_modified");
2947 }
2948 else {
2949 if (thisProc_ == 0) {
2950 FEI_COUT << "Aztec_LinSysCore: Warning: requested orthog. <" << param << "> not recognized." << FEI_ENDL;
2951 }
2952 }
2953
2954 debugOutput(msg);
2955
2956 return;
2957}
2958
2959//==============================================================================
2960void Aztec_LinSysCore::setAuxVec(const char* param)
2961{
2962 char msg[64];
2963
2964 if (!strcmp(param, "AZ_resid")) {
2965 aztec_options_[AZ_aux_vec] = AZ_resid;
2966 sprintf(msg, "AZ_aux_vec AZ_resid");
2967 }
2968 else if (!strcmp(param, "AZ_rand")) {
2969 aztec_options_[AZ_aux_vec] = AZ_rand;
2970 sprintf(msg, "AZ_aux_vec AZ_rand");
2971 }
2972 else {
2973 if (thisProc_ == 0) {
2974 FEI_COUT << "Aztec_LinSysCore: Warning: requested aux_vec <" << param << "> not recognized." << FEI_ENDL;
2975 }
2976 }
2977
2978 debugOutput(msg);
2979
2980 return;
2981}
2982
2983//==============================================================================
2984void Aztec_LinSysCore::setAZ_output(const char* param)
2985{
2986 char msg[64];
2987 int out = -1;
2988 int num = sscanf(param, "%d", &out);
2989 if (num == 1 && out > -1) {
2990 sprintf(msg, "AZ_output %d", out);
2991 aztec_options_[AZ_output] = out;
2992 }
2993 else if (!strcmp(param, "AZ_all")) {
2994 aztec_options_[AZ_output] = AZ_all;
2995 sprintf(msg, "AZ_output AZ_all");
2996 }
2997 else if (!strcmp(param, "AZ_none")) {
2998 aztec_options_[AZ_output] = AZ_none;
2999 sprintf(msg, "AZ_output AZ_none");
3000 }
3001 else if (!strcmp(param, "AZ_warnings")) {
3002 aztec_options_[AZ_output] = AZ_warnings;
3003 sprintf(msg, "AZ_output AZ_warnings");
3004 }
3005 else if (!strcmp(param, "AZ_last")) {
3006 aztec_options_[AZ_output] = AZ_last;
3007 sprintf(msg, "AZ_output AZ_last");
3008 }
3009 else {
3010 if (thisProc_ == 0) {
3011 FEI_COUT << "Aztec_LinSysCore: Warning: requested AZ_output <" << param << "> not recognized." << FEI_ENDL;
3012 }
3013 }
3014
3015 debugOutput(msg);
3016
3017 return;
3018}
3019
3020//==============================================================================
3021void Aztec_LinSysCore::checkForParam(const char* paramName,
3022 int numParams, char** paramStrings,
3023 double& param) {
3024 const char* parameter =
3025 snl_fei::getParamValue(paramName, numParams, paramStrings);
3026 if (parameter != NULL) {
3027 sscanf(parameter, "%le", &param);
3028 }
3029}
3030
3031//==============================================================================
3032void Aztec_LinSysCore::checkForOption(const char* paramName,
3033 int numParams, char** paramStrings,
3034 int& param) {
3035 const char* parameter =
3036 snl_fei::getParamValue(paramName, numParams, paramStrings);
3037 if (parameter != NULL) {
3038 sscanf(parameter, "%d", &param);
3039 }
3040}
3041
3042//==============================================================================
3043void Aztec_LinSysCore::setDebugOutput(const char* path, const char* name){
3044//
3045//This function turns on debug output, and opens a file to put it in.
3046//
3047 if (debugOutput_) {
3048 fprintf(debugFile_,"setDebugOutput closing this file.");
3049 fflush(debugFile_);
3050 fclose(debugFile_);
3051 debugFile_ = NULL;
3052 }
3053
3054 int pathLength = strlen(path);
3055 if (path != debugPath_) {
3056 delete [] debugPath_;
3057 debugPath_ = new char[pathLength + 1];
3058 sprintf(debugPath_, "%s", path);
3059 }
3060
3061 int nameLength = strlen(name);
3062 if (name != debugFileName_) {
3063 delete [] debugFileName_;
3064 debugFileName_ = new char[nameLength + 1];
3065 sprintf(debugFileName_,"%s",name);
3066 }
3067
3068 char* dbFileName = new char[pathLength + nameLength + 3];
3069
3070 sprintf(dbFileName, "%s/%s", path, name);
3071
3072 debugOutput_ = 1;
3073 debugFile_ = fopen(dbFileName, "w");
3074
3075 if (!debugFile_){
3076 fei::console_out() << "couldn't open debug output file: " << dbFileName << FEI_ENDL;
3077 debugOutput_ = 0;
3078 delete [] debugPath_;
3079 debugPath_ = NULL;
3080 delete [] debugFileName_;
3081 debugFileName_ = NULL;
3082 }
3083
3084 delete [] dbFileName;
3085}
3086
3087//==============================================================================
3088int Aztec_LinSysCore::VBRmatPlusScaledMat(AztecDVBR_Matrix* A,
3089 double scalar,
3090 AztecDVBR_Matrix* source)
3091{
3092 int* nnz = new int[numLocalEqnBlks_];
3093 int* nblk = new int[numLocalEqnBlks_];
3094 int* src_nnz = new int[numLocalEqnBlks_];
3095 int* src_nblk = new int[numLocalEqnBlks_];
3096
3097 if (nnz == NULL || nblk == NULL || src_nnz==NULL || src_nblk==NULL) {
3098 messageAbort("VBRMatPlusScaledMat: allocation failed");
3099 }
3100
3101 A->getNumNonzerosPerRow(nnz);
3102 A->getNumBlocksPerRow(nblk);
3103 source->getNumNonzerosPerRow(src_nnz);
3104 source->getNumBlocksPerRow(src_nblk);
3105
3106 int i, max_nnz = 0, max_nblk = 0;
3107 for(i=0; i<numLocalEqnBlks_; i++) {
3108 if (nnz[i] != src_nnz[i] || nblk[i] != src_nblk[i]) {
3109 messageAbort("VBRmatPlusScaledMat: matrix sizes don't match.");
3110 }
3111 if (max_nnz < nnz[i]) max_nnz = nnz[i];
3112 if (max_nblk < nblk[i]) max_nblk = nblk[i];
3113 }
3114
3115 delete [] nnz;
3116 delete [] nblk;
3117 delete [] src_nnz;
3118 delete [] src_nblk;
3119
3120 double* val = new double[max_nnz];
3121 int* colInds = new int[max_nblk];
3122 if (val==NULL || colInds==NULL) {
3123 messageAbort("VBRmatPlusScaledMat: allocation failed");
3124 }
3125 int len, nnzBlks;
3126
3127 for(i=0; i<numLocalEqnBlks_; i++) {
3128 int row = localBlkOffset_+i;
3129 int err = source->getNumBlocksPerRow(row, nnzBlks);
3130 err += source->getNumNonzerosPerRow(row, len);
3131 err += source->getBlockRow(row, val, colInds, nnzBlks);
3132
3133 if (err) messageAbort("VBRmatPlusScaledMat: error getting src row");
3134
3135 for(int j=0; j<len; j++) val[j] *= scalar;
3136
3137 err = A->sumIntoBlockRow(row, val, colInds, nnzBlks);
3138 if (err) messageAbort("VBRmatPlusScaledMat: error summing in row");
3139 }
3140
3141 delete [] val;
3142 delete [] colInds;
3143 return(0);
3144}
3145
3146//==============================================================================
3147int Aztec_LinSysCore::MSRmatPlusScaledMat(AztecDMSR_Matrix* A,
3148 double scalar,
3149 AztecDMSR_Matrix* source)
3150{
3151 return(A->addScaledMatrix(scalar, *source));
3152}
3153
3154//==============================================================================
3155void Aztec_LinSysCore::debugOutput(const char* msg) const {
3156 if (debugOutput_) {
3157 fprintf(debugFile_, "%s\n", msg);
3158 fflush(debugFile_);
3159 }
3160}
3161
3162//==============================================================================
3163int Aztec_LinSysCore::messageAbort(const char* msg) const {
3164 fei::console_out() << "Aztec_LinSysCore: " << msg << " Aborting." << FEI_ENDL;
3165#ifndef FEI_SER
3166 MPI_Abort(comm_, -1);
3167#else
3168 abort();
3169#endif
3170 return(-1);
3171}
3172
3173}//namespace fei_trilinos
3174
3175#endif
3176//HAVE_FEI_AZTECOO
const char * getTypeName() const
Definition: fei_Data.hpp:32
void setDataPtr(void *ptr)
Definition: fei_Data.hpp:35
void * getDataPtr() const
Definition: fei_Data.hpp:38
void setTypeName(const char *name)
Definition: fei_Data.hpp:28
#define CHK_ERR(a)
int GlobalID
Definition: fei_defs.h:60
#define FEI_ENDL
#define FEI_COUT
#define MPI_Comm
Definition: fei_mpi.h:56
#define MPI_Abort(a, b)
Definition: fei_mpi.h:59
#define FEI_OSTRINGSTREAM
Definition: fei_sstream.hpp:32
int searchList(const T &item, const T *list, int len)
int binarySearch(const T &item, const T *list, int len)
std::ostream & console_out()
int mergeStringLists(char **&strings, int &numStrings, const char *const *stringsToMerge, int numStringsToMerge)
const char * getParam(const char *key, int numParams, const char *const *paramStrings)
const char * getParamValue(const char *key, int numParams, const char *const *paramStrings, char separator=' ')