FEI  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
fei_Aztec_LinSysCore.cpp
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 
44 #include <fei_trilinos_macros.hpp>
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"
57 #include "fei_LinearSystemCore.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"
80 #include "fei_Aztec_BlockMap.hpp"
81 #include "fei_Aztec_LSVector.hpp"
82 #include "fei_AztecDMSR_Matrix.hpp"
83 #include "fei_AztecDVBR_Matrix.hpp"
84 
85 #include "fei_Aztec_LinSysCore.hpp"
86 
87 namespace fei_trilinos {
88 
89 static int azlsc_solveCounter_ = 0;
90 static int azlsc_debugFileCounter_ = 0;
91 
92 static std::map<std::string, unsigned> fei_aztec_named_solve_counter;
93 
94 //=========CONSTRUCTOR==========================================================
95 Aztec_LinSysCore::Aztec_LinSysCore(MPI_Comm comm)
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 
187 struct 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============================================================
196 Aztec_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 //==============================================================================
260 LinearSystemCore* Aztec_LinSysCore::clone() {
261  return(new Aztec_LinSysCore(comm_));
262 }
263 
264 //==============================================================================
265 int 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 //==============================================================================
378 int Aztec_LinSysCore::setLookup(Lookup& lookup)
379 {
380  lookup_ = &lookup;
381  haveLookup_ = true;
382  return(0);
383 }
384 
385 //==============================================================================
386 int 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 //==============================================================================
411 int 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 //==============================================================================
426 int 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 //==============================================================================
452 int 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 //==============================================================================
474 int 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 //==============================================================================
540 int 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 //==============================================================================
591 int 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 //==============================================================================
645 int 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 //==============================================================================
686 int 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 //==============================================================================
708 int 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 //==============================================================================
736 int 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 //==============================================================================
746 int 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 //==============================================================================
768 int 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 //==============================================================================
778 int 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 //==============================================================================
797 int 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 //==============================================================================
833 int 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 //==============================================================================
885 int 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 //==============================================================================
915 int 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 //==============================================================================
940 int 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 //==============================================================================
1016 int 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 //==============================================================================
1046 int 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 //==============================================================================
1078 int 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 //==============================================================================
1106 int 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 //==============================================================================
1131 int 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 //==============================================================================
1245 int 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 //==============================================================================
1267 int 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 //==============================================================================
1437 int 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 //==============================================================================
1539 int 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 //==============================================================================
1629 int 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 //==============================================================================
1681 int 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 //==============================================================================
1700 int 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 //==============================================================================
1744 int 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 //==============================================================================
1844 int 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 //==============================================================================
1913 int 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 //==============================================================================
1929 int 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 //==============================================================================
1965 int 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 //==============================================================================
1995 int 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 //==============================================================================
2023 int 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 //==============================================================================
2033 int 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 //==============================================================================
2049 int 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 //==============================================================================
2065 int 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 //==============================================================================
2079 int 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 //==============================================================================
2099 int 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 //==============================================================================
2110 int 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 //==============================================================================
2126 int 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 //==============================================================================
2141 int 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 //==============================================================================
2173 int 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 //==============================================================================
2189 int 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 //==============================================================================
2203 int 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 //==============================================================================
2272 int 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 //==============================================================================
2317 int 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 //==============================================================================
2369 void 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 //==============================================================================
2409 void 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 //==============================================================================
2433 int Aztec_LinSysCore::writeSystem(const char* name)
2434 {
2435  writeA(name);
2436  return(0);
2437 }
2438 
2439 //==============================================================================
2440 void 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 //==============================================================================
2459 void 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 //==============================================================================
2545 int 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 //==============================================================================
2576 int 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 //==============================================================================
2601 int 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 //==============================================================================
2620 int 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 //==============================================================================
2637 int 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 //==============================================================================
2798 void 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 //==============================================================================
2841 void 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 //==============================================================================
2881 void 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 //==============================================================================
2909 void 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 //==============================================================================
2936 void 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 //==============================================================================
2960 void 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 //==============================================================================
2984 void 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 //==============================================================================
3021 void 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 //==============================================================================
3032 void 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 //==============================================================================
3043 void 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 //==============================================================================
3088 int 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 //==============================================================================
3147 int Aztec_LinSysCore::MSRmatPlusScaledMat(AztecDMSR_Matrix* A,
3148  double scalar,
3149  AztecDMSR_Matrix* source)
3150 {
3151  return(A->addScaledMatrix(scalar, *source));
3152 }
3153 
3154 //==============================================================================
3155 void Aztec_LinSysCore::debugOutput(const char* msg) const {
3156  if (debugOutput_) {
3157  fprintf(debugFile_, "%s\n", msg);
3158  fflush(debugFile_);
3159  }
3160 }
3161 
3162 //==============================================================================
3163 int 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
void setTypeName(const char *name)
Definition: fei_Data.hpp:28
int binarySearch(const T &item, const T *list, int len)
void * getDataPtr() const
Definition: fei_Data.hpp:38
std::ostream & console_out()
void setDataPtr(void *ptr)
Definition: fei_Data.hpp:35
const char * getTypeName() const
Definition: fei_Data.hpp:32
int searchList(const T &item, const T *list, int len)