1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 //***************************************************************************
9 // system includes
10 //---------------------------------------------------------------------------
11
12 #include "utilities/_hypre_utilities.h"
13 #include <stdlib.h>
14 #include <string.h>
15 #include <stdio.h>
16 #include <math.h>
17
18 #if 0 /* RDF: Not sure this is really needed */
19 #ifdef WIN32
20 #define strcmp _stricmp
21 #endif
22 #endif
23
24 //***************************************************************************
25 // HYPRE includes
26 //---------------------------------------------------------------------------
27
28 #include "HYPRE.h"
29 #include "IJ_mv/HYPRE_IJ_mv.h"
30 #include "parcsr_mv/HYPRE_parcsr_mv.h"
31 #include "parcsr_ls/HYPRE_parcsr_ls.h"
32 #include "HYPRE_parcsr_bicgstabl.h"
33 #include "HYPRE_parcsr_TFQmr.h"
34 #include "HYPRE_parcsr_bicgs.h"
35 #include "HYPRE_parcsr_symqmr.h"
36 #include "HYPRE_parcsr_fgmres.h"
37 #include "HYPRE_parcsr_lsicg.h"
38 #include "HYPRE_LinSysCore.h"
39 #include "parcsr_mv/_hypre_parcsr_mv.h"
40 #include "HYPRE_LSI_schwarz.h"
41 #include "HYPRE_LSI_ddilut.h"
42 #include "HYPRE_LSI_ddict.h"
43 #include "HYPRE_LSI_poly.h"
44 #include "HYPRE_LSI_block.h"
45 #include "HYPRE_LSI_Uzawa_c.h"
46 #include "HYPRE_LSI_Dsuperlu.h"
47 #include "HYPRE_MLMaxwell.h"
48 #include "HYPRE_SlideReduction.h"
49
50 //***************************************************************************
51 // timers
52 //---------------------------------------------------------------------------
53
54 #ifdef HYPRE_SEQUENTIAL
55 #include <time.h>
56 extern "C"
57 {
LSC_Wtime()58 double LSC_Wtime()
59 {
60 clock_t ticks;
61 double seconds;
62 ticks = clock() ;
63 seconds = (double) ticks / (double) CLOCKS_PER_SEC;
64 return seconds;
65 }
66 }
67 #else
LSC_Wtime()68 double LSC_Wtime()
69 {
70 return (MPI_Wtime());
71 }
72 #endif
73 extern "C" {
74 int HYPRE_LSI_qsort1a( int *, int *, int, int );
75 int HYPRE_LSI_PartitionMatrix(int, int, int*, int**, double**, int*, int**);
76 int HYPRE_LSI_GetMatrixDiagonal(int, int, int *, int **, double **, int *,
77 int *, double *);
78 }
79
80 //***************************************************************************
81 // These are external functions needed internally here
82 //---------------------------------------------------------------------------
83
84 #include "HYPRE_LSI_mli.h"
85
86 extern "C" {
87
88 #ifdef HAVE_ML
89 int HYPRE_LSI_MLCreate( MPI_Comm, HYPRE_Solver *);
90 int HYPRE_LSI_MLDestroy( HYPRE_Solver );
91 #endif
92
93 #ifdef HAVE_MLMAXWELL
94 int HYPRE_LSI_MLMaxwellCreate(MPI_Comm, HYPRE_Solver *);
95 int HYPRE_LSI_MLMaxwellDestroy(HYPRE_Solver );
96 #endif
97
98 #ifdef HAVE_AMGE
99 int HYPRE_LSI_AMGeCreate();
100 int HYPRE_LSI_AMGeDestroy();
101 int HYPRE_LSI_AMGeSetNNodes(int);
102 int HYPRE_LSI_AMGeSetNElements(int);
103 int HYPRE_LSI_AMGeSetSystemSize(int);
104 int HYPRE_LSI_AMGePutRow(int,int,double*,int*);
105 int HYPRE_LSI_AMGeSolve( double *rhs, double *sol );
106 int HYPRE_LSI_AMGeSetBoundary( int leng, int *colInd );
107 int HYPRE_LSI_AMGeWriteToFile();
108 #endif
109
110 void hypre_qsort0(int *, int, int);
111 void hypre_qsort1(int *, double *, int, int);
112 }
113
114 #define habs(x) ( ( (x) > 0 ) ? x : -(x))
115
116 //***************************************************************************
117 // constructor
118 //---------------------------------------------------------------------------
119
HYPRE_LinSysCore(MPI_Comm comm)120 HYPRE_LinSysCore::HYPRE_LinSysCore(MPI_Comm comm) :
121 comm_(comm),
122 HYOutputLevel_(0),
123 memOptimizerFlag_(0),
124 mapFromSolnFlag_(0),
125 mapFromSolnLeng_(0),
126 mapFromSolnLengMax_(0),
127 mapFromSolnList_(NULL),
128 mapFromSolnList2_(NULL),
129 HYA_(NULL),
130 HYnormalA_(NULL),
131 HYb_(NULL),
132 HYnormalB_(NULL),
133 HYbs_(NULL),
134 HYx_(NULL),
135 HYr_(NULL),
136 HYpxs_(NULL),
137 HYpbs_(NULL),
138 numGlobalRows_(0),
139 localStartRow_(0),
140 localEndRow_(-1),
141 localStartCol_(-1),
142 localEndCol_(-1),
143 rowLengths_(NULL),
144 colIndices_(NULL),
145 colValues_(NULL),
146 reducedA_(NULL),
147 reducedB_(NULL),
148 reducedX_(NULL),
149 reducedR_(NULL),
150 HYA21_(NULL),
151 HYA12_(NULL),
152 A21NRows_(0),
153 A21NCols_(0),
154 HYinvA22_(NULL),
155 currA_(NULL),
156 currB_(NULL),
157 currX_(NULL),
158 currR_(NULL),
159 currentRHS_(0),
160 numRHSs_(1),
161 nStored_(0),
162 storedIndices_(NULL),
163 auxStoredIndices_(NULL),
164 mRHSFlag_(0),
165 mRHSNumGEqns_(0),
166 mRHSGEqnIDs_(NULL),
167 mRHSNEntries_(NULL),
168 mRHSBCType_(NULL),
169 mRHSRowInds_(NULL),
170 mRHSRowVals_(NULL),
171 matrixVectorsCreated_(0),
172 systemAssembled_(0),
173 slideReduction_(0),
174 slideReductionMinNorm_(-1.0),
175 slideReductionScaleMatrix_(0),
176 schurReduction_(0),
177 schurReductionCreated_(0),
178 projectionScheme_(0),
179 projectSize_(0),
180 projectCurrSize_(0),
181 projectionMatrix_(NULL),
182 normalEqnFlag_(0),
183 slideObj_(NULL),
184 selectedList_(NULL),
185 selectedListAux_(NULL),
186 nConstraints_(0),
187 constrList_(NULL),
188 matrixPartition_(0),
189 HYSolver_(NULL),
190 maxIterations_(1000),
191 tolerance_(1.0e-6),
192 normAbsRel_(0),
193 pcgRecomputeRes_(0),
194 HYPrecon_(NULL),
195 HYPreconReuse_(0),
196 HYPreconSetup_(0),
197 lookup_(NULL),
198 haveLookup_(0)
199 {
200 //-------------------------------------------------------------------
201 // find my processor ID
202 //-------------------------------------------------------------------
203
204 MPI_Comm_rank(comm, &mypid_);
205 MPI_Comm_size(comm, &numProcs_);
206
207 //-------------------------------------------------------------------
208 // default method = gmres
209 //-------------------------------------------------------------------
210
211 HYSolverName_ = new char[64];
212 strcpy(HYSolverName_,"gmres");
213 HYSolverID_ = HYGMRES;
214
215 //-------------------------------------------------------------------
216 // default preconditioner = identity
217 //-------------------------------------------------------------------
218
219 HYPreconName_ = new char[64];
220 strcpy(HYPreconName_,"diagonal");
221 HYPreconID_ = HYDIAGONAL;
222
223 //-------------------------------------------------------------------
224 // parameters for controlling amg, pilut, SuperLU, etc.
225 //-------------------------------------------------------------------
226
227 gmresDim_ = 100; // restart size in GMRES
228 fgmresUpdateTol_ = 0;
229
230 amgMaxLevels_ = 30; // default max number of levels
231 amgCoarsenType_ = 0; // default coarsening
232 amgMeasureType_ = 0; // local measure
233 amgSystemSize_ = 1; // system size
234 amgMaxIter_ = 1; // number of iterations
235 amgNumSweeps_[0] = 1; // no. of sweeps for fine grid
236 amgNumSweeps_[1] = 1; // no. of presmoothing sweeps
237 amgNumSweeps_[2] = 1; // no. of postsmoothing sweeps
238 amgNumSweeps_[3] = 1; // no. of sweeps for coarsest grid
239 amgRelaxType_[0] = 3; // hybrid for the fine grid
240 amgRelaxType_[1] = 3; // hybrid for presmoothing
241 amgRelaxType_[2] = 3; // hybrid for postsmoothing
242 amgRelaxType_[3] = 9; // direct for the coarsest level
243 amgGridRlxType_ = 0; // smoothes all points
244 amgStrongThreshold_ = 0.25;
245 amgSmoothType_ = 0; // default non point smoother, none
246 amgSmoothNumLevels_ = 0; // no. of levels for non point smoothers
247 amgSmoothNumSweeps_ = 1; // no. of sweeps for non point smoothers
248 amgCGSmoothNumSweeps_ = 0; // no. of sweeps for preconditioned CG smoother
249 amgSchwarzRelaxWt_ = 1.0; // relaxation weight for Schwarz smoother
250 amgSchwarzVariant_ = 0; // hybrid multiplicative Schwarz with
251 // no overlap across processor boundaries
252 amgSchwarzOverlap_ = 1; // minimal overlap
253 amgSchwarzDomainType_ = 2; // domain through agglomeration
254 amgUseGSMG_ = 0;
255 amgGSMGNSamples_ = 0;
256 amgAggLevels_ = 0;
257 amgInterpType_ = 0;
258 amgPmax_ = 0;
259
260 for (int i = 0; i < 25; i++) amgRelaxWeight_[i] = 1.0;
261 for (int j = 0; j < 25; j++) amgRelaxOmega_[j] = 1.0;
262
263 pilutFillin_ = 0; // how many nonzeros to keep in L and U
264 pilutDropTol_ = 0.0;
265 pilutMaxNnzPerRow_ = 0; // register the max NNZ/row in matrix A
266
267 ddilutFillin_ = 1.0; // additional fillin other than A
268 ddilutDropTol_ = 1.0e-8;
269 ddilutOverlap_ = 0;
270 ddilutReorder_ = 0;
271
272 ddictFillin_ = 1.0; // additional fillin other than A
273 ddictDropTol_ = 1.0e-8;
274
275 schwarzFillin_ = 1.0; // additional fillin other than A
276 schwarzNblocks_ = 1;
277 schwarzBlksize_ = 0;
278
279 polyOrder_ = 8; // order of polynomial preconditioner
280
281 parasailsSym_ = 0; // default is nonsymmetric
282 parasailsThreshold_ = 0.1;
283 parasailsNlevels_ = 1;
284 parasailsFilter_ = 0.05;
285 parasailsLoadbal_ = 0.0;
286 parasailsReuse_ = 0; // reuse pattern if nonzero
287
288 euclidargc_ = 2; // parameters information for Euclid
289 euclidargv_ = new char*[euclidargc_*2];
290 for (int k = 0; k < euclidargc_*2; k++) euclidargv_[k] = new char[50];
291 strcpy(euclidargv_[0], "-level");
292 strcpy(euclidargv_[1], "0");
293 strcpy(euclidargv_[2], "-sparseA");
294 strcpy(euclidargv_[3], "0.0");
295
296 superluOrdering_ = 0; // natural ordering in SuperLU
297 superluScale_[0] = 'N'; // no scaling in SuperLUX
298
299 mlMethod_ = 1;
300 mlNumPreSweeps_ = 1;
301 mlNumPostSweeps_ = 1;
302 mlPresmootherType_ = 1; // default Gauss-Seidel
303 mlPostsmootherType_ = 1; // default Gauss-Seidel
304 mlRelaxWeight_ = 0.5;
305 mlStrongThreshold_ = 0.08; // one suggested by Vanek/Brezina/Mandel
306 mlCoarseSolver_ = 0; // default coarse solver = SuperLU
307 mlCoarsenScheme_ = 1; // default coarsening scheme = uncoupled
308 mlNumPDEs_ = 3; // default block size
309
310 truncThresh_ = 0.0;
311 rnorm_ = 0.0;
312 rhsIDs_ = new int[1];
313 rhsIDs_[0] = 0;
314 feData_ = NULL;
315 haveFEData_ = 0;
316 feData_ = NULL;
317 MLI_NumNodes_ = 0;
318 MLI_FieldSize_ = 0;
319 MLI_NodalCoord_ = NULL;
320 MLI_EqnNumbers_ = NULL;
321 MLI_Hybrid_GSA_ = 0;
322 MLI_Hybrid_NSIncr_ = 2;
323 MLI_Hybrid_MaxIter_ = 100;
324 MLI_Hybrid_ConvRate_ = 0.95;
325 MLI_Hybrid_NTrials_ = 5;
326 AMSData_.numNodes_ = 0;
327 AMSData_.numLocalNodes_ = 0;
328 AMSData_.EdgeNodeList_ = NULL;
329 AMSData_.NodeNumbers_ = NULL;
330 AMSData_.NodalCoord_ = NULL;
331 amsX_ = NULL;
332 amsY_ = NULL;
333 amsZ_ = NULL;
334 amsG_ = NULL;
335 amsD0_ = NULL;
336 amsD1_ = NULL;
337 amsNumPDEs_ = 3;
338 amsMaxIter_ = 1;
339 amsTol_ = 0.0;
340 amsCycleType_ = 1;
341 amsRelaxType_ = 2;
342 amsRelaxTimes_ = 1;
343 amsRelaxWt_ = 1.0;
344 amsRelaxOmega_ = 1.0;
345 amsBetaPoisson_ = NULL;
346 amsPrintLevel_ = 0;
347 amsAlphaCoarsenType_ = 10;
348 amsAlphaAggLevels_ = 1;
349 amsAlphaRelaxType_ = 6;
350 amsAlphaStrengthThresh_ = 0.25;
351 amsBetaCoarsenType_ = 10;
352 amsBetaAggLevels_ = 1;
353 amsBetaRelaxType_ = 6;
354 amsBetaStrengthThresh_ = 0.25;
355 FEI_mixedDiagFlag_ = 0;
356 FEI_mixedDiag_ = NULL;
357 sysPDEMethod_ = -1;
358 sysPDEFormat_ = -1;
359 sysPDETol_ = 0.0;
360 sysPDEMaxIter_ = -1;
361 sysPDENumPre_ = -1;
362 sysPDENumPost_ = -1;
363 sysPDENVars_ = 3;
364
365 //-------------------------------------------------------------------
366 // parameters ML Maxwell solver
367 //-------------------------------------------------------------------
368
369 maxwellANN_ = NULL;
370 maxwellGEN_ = NULL;
371 }
372
373 //***************************************************************************
374 // destructor
375 //---------------------------------------------------------------------------
376
~HYPRE_LinSysCore()377 HYPRE_LinSysCore::~HYPRE_LinSysCore()
378 {
379 int i;
380 //-------------------------------------------------------------------
381 // diagnostic message
382 //-------------------------------------------------------------------
383
384 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
385 printf("%4d : HYPRE_LSC::entering destructor.\n",mypid_);
386
387 //-------------------------------------------------------------------
388 // clean up the allocated matrix and vectors
389 //-------------------------------------------------------------------
390
391 if ( HYA_ != NULL ) {HYPRE_IJMatrixDestroy(HYA_); HYA_ = NULL;}
392 if ( HYx_ != NULL ) {HYPRE_IJVectorDestroy(HYx_); HYx_ = NULL;}
393 if ( HYr_ != NULL ) {HYPRE_IJVectorDestroy(HYr_); HYr_ = NULL;}
394 if ( HYbs_ != NULL )
395 {
396 for ( i = 0; i < numRHSs_; i++ )
397 if ( HYbs_[i] != NULL ) HYPRE_IJVectorDestroy(HYbs_[i]);
398 delete [] HYbs_;
399 HYbs_ = NULL;
400 }
401 if ( HYpbs_ != NULL )
402 {
403 for ( i = 0; i <= projectSize_; i++ )
404 if ( HYpbs_[i] != NULL ) HYPRE_IJVectorDestroy(HYpbs_[i]);
405 delete [] HYpbs_;
406 HYpbs_ = NULL;
407 }
408 if ( HYpxs_ != NULL )
409 {
410 for ( i = 0; i <= projectSize_; i++ )
411 if ( HYpxs_[i] != NULL ) HYPRE_IJVectorDestroy(HYpxs_[i]);
412 delete [] HYpxs_;
413 HYpxs_ = NULL;
414 }
415 if (HYnormalA_!= NULL) {HYPRE_IJMatrixDestroy(HYnormalA_);HYnormalA_ = NULL;}
416 if (HYnormalB_!= NULL) {HYPRE_IJVectorDestroy(HYnormalB_);HYnormalB_ = NULL;}
417 if (reducedA_ != NULL) {HYPRE_IJMatrixDestroy(reducedA_); reducedA_ = NULL;}
418 if (reducedB_ != NULL) {HYPRE_IJVectorDestroy(reducedB_); reducedB_ = NULL;}
419 if (reducedX_ != NULL) {HYPRE_IJVectorDestroy(reducedX_); reducedX_ = NULL;}
420 if (reducedR_ != NULL) {HYPRE_IJVectorDestroy(reducedR_); reducedR_ = NULL;}
421 if (HYA21_ != NULL) {HYPRE_IJMatrixDestroy(HYA21_); HYA21_ = NULL;}
422 if (HYA12_ != NULL) {HYPRE_IJMatrixDestroy(HYA12_); HYA12_ = NULL;}
423 if (HYinvA22_ != NULL) {HYPRE_IJMatrixDestroy(HYinvA22_); HYinvA22_ = NULL;}
424
425 matrixVectorsCreated_ = 0;
426 systemAssembled_ = 0;
427 projectCurrSize_ = 0;
428
429 if ( colIndices_ != NULL )
430 {
431 for ( i = 0; i < localEndRow_-localStartRow_+1; i++ )
432 if ( colIndices_[i] != NULL ) delete [] colIndices_[i];
433 delete [] colIndices_;
434 colIndices_ = NULL;
435 }
436 if ( colValues_ != NULL )
437 {
438 for ( i = 0; i < localEndRow_-localStartRow_+1; i++ )
439 if ( colValues_[i] != NULL ) delete [] colValues_[i];
440 delete [] colValues_;
441 colValues_ = NULL;
442 }
443 if ( rowLengths_ != NULL )
444 {
445 delete [] rowLengths_;
446 rowLengths_ = NULL;
447 }
448 if ( rhsIDs_ != NULL ) delete [] rhsIDs_;
449 if ( storedIndices_ != NULL ) delete [] storedIndices_;
450 if ( auxStoredIndices_ != NULL ) delete [] auxStoredIndices_;
451 if ( mRHSNumGEqns_ > 0)
452 {
453 if (mRHSGEqnIDs_ != NULL) delete [] mRHSGEqnIDs_;
454 if (mRHSNEntries_ != NULL) delete [] mRHSNEntries_;
455 if (mRHSBCType_ != NULL) delete [] mRHSBCType_ ;
456 if (mRHSRowInds_ != NULL)
457 {
458 for (i = 0; i < mRHSNumGEqns_; i++) delete [] mRHSRowInds_[i];
459 delete [] mRHSRowInds_;
460 }
461 if (mRHSRowVals_ != NULL)
462 {
463 for (i = 0; i < mRHSNumGEqns_; i++) delete [] mRHSRowVals_[i];
464 delete [] mRHSRowVals_;
465 }
466 mRHSNumGEqns_ = 0;
467 mRHSGEqnIDs_ = NULL;
468 mRHSNEntries_ = NULL;
469 mRHSBCType_ = NULL;
470 mRHSRowInds_ = NULL;
471 mRHSRowVals_ = NULL;
472 }
473
474 //-------------------------------------------------------------------
475 // clean up direct matrix access variables
476 //-------------------------------------------------------------------
477
478 if ( mapFromSolnList_ != NULL )
479 {
480 delete [] mapFromSolnList_;
481 mapFromSolnList_ = NULL;
482 }
483 if ( mapFromSolnList2_ != NULL )
484 {
485 delete [] mapFromSolnList2_;
486 mapFromSolnList2_ = NULL;
487 }
488
489 //-------------------------------------------------------------------
490 // call solver destructors
491 //-------------------------------------------------------------------
492
493 if ( HYSolver_ != NULL )
494 {
495 if (HYSolverID_ == HYPCG) HYPRE_ParCSRPCGDestroy(HYSolver_);
496 if (HYSolverID_ == HYGMRES) HYPRE_ParCSRGMRESDestroy(HYSolver_);
497 if (HYSolverID_ == HYCGSTAB) HYPRE_ParCSRBiCGSTABDestroy(HYSolver_);
498 if (HYSolverID_ == HYCGSTABL) HYPRE_ParCSRBiCGSTABLDestroy(HYSolver_);
499 if (HYSolverID_ == HYAMG) HYPRE_BoomerAMGDestroy(HYSolver_);
500 if (HYSolverID_ == HYTFQMR) HYPRE_ParCSRTFQmrDestroy(HYSolver_);
501 HYSolver_ = NULL;
502 }
503 delete [] HYSolverName_;
504 HYSolverName_ = NULL;
505 #ifdef HAVE_AMGE
506 if ( HYSolverID_ == HYAMGE ) HYPRE_LSI_AMGeDestroy();
507 #endif
508
509 //-------------------------------------------------------------------
510 // call preconditioner destructors
511 //-------------------------------------------------------------------
512
513 if ( HYPrecon_ != NULL )
514 {
515 if ( HYPreconID_ == HYPILUT )
516 HYPRE_ParCSRPilutDestroy( HYPrecon_ );
517
518 else if ( HYPreconID_ == HYPARASAILS )
519 HYPRE_ParCSRParaSailsDestroy( HYPrecon_ );
520
521 else if ( HYPreconID_ == HYBOOMERAMG )
522 HYPRE_BoomerAMGDestroy( HYPrecon_ );
523
524 else if ( HYPreconID_ == HYDDILUT )
525 HYPRE_LSI_DDIlutDestroy( HYPrecon_ );
526
527 else if ( HYPreconID_ == HYSCHWARZ )
528 HYPRE_LSI_SchwarzDestroy( HYPrecon_ );
529
530 else if ( HYPreconID_ == HYPOLY )
531 HYPRE_LSI_PolyDestroy( HYPrecon_ );
532
533 else if ( HYPreconID_ == HYEUCLID )
534 HYPRE_EuclidDestroy( HYPrecon_ );
535
536 else if ( HYPreconID_ == HYBLOCK )
537 HYPRE_LSI_BlockPrecondDestroy( HYPrecon_ );
538
539 #ifdef HAVE_ML
540 else if ( HYPreconID_ == HYML )
541 HYPRE_LSI_MLDestroy( HYPrecon_ );
542 #endif
543
544 #ifdef HAVE_MLMAXWELL
545 else if ( HYPreconID_ == HYMLMAXWELL )
546 HYPRE_LSI_MLMaxwellDestroy( HYPrecon_ );
547 #endif
548 else if ( HYPreconID_ == HYMLI )
549 HYPRE_LSI_MLIDestroy( HYPrecon_ );
550
551 else if ( HYPreconID_ == HYAMS )
552 {
553 // Destroy G and coordinate vectors
554 // OLD WAY
555 if( amsG_ == NULL ) {
556 HYPRE_AMSFEIDestroy( HYPrecon_ );
557 }
558 HYPRE_AMSDestroy( HYPrecon_ );
559 }
560 #ifdef HAVE_SYSPDE
561 else if ( HYPreconID_ == HYSYSPDE )
562 HYPRE_ParCSRSysPDEDestroy( HYPrecon_ );
563 #endif
564 #ifdef HYPRE_USING_DSUPERLU
565 else if ( HYPreconID_ == HYDSLU )
566 HYPRE_LSI_DSuperLUDestroy(HYPrecon_);
567 #endif
568
569 HYPrecon_ = NULL;
570 }
571 delete [] HYPreconName_;
572 HYPreconName_ = NULL;
573
574 for (i = 0; i < euclidargc_*2; i++) delete [] euclidargv_[i];
575 delete [] euclidargv_;
576 euclidargv_ = NULL;
577
578 //-------------------------------------------------------------------
579 // clean up variable for various reduction
580 //-------------------------------------------------------------------
581
582 if ( constrList_ != NULL )
583 {
584 delete [] constrList_;
585 constrList_ = NULL;
586 }
587 if (selectedList_ != NULL)
588 {
589 delete [] selectedList_;
590 selectedList_ = NULL;
591 }
592 if (selectedListAux_ != NULL)
593 {
594 delete [] selectedListAux_;
595 selectedListAux_ = NULL;
596 }
597
598 //-------------------------------------------------------------------
599 // deallocate local storage for MLI
600 //-------------------------------------------------------------------
601
602 #ifdef HAVE_MLI
603 if ( feData_ != NULL )
604 {
605 if (haveFEData_ == 1) HYPRE_LSI_MLIFEDataDestroy(feData_);
606 else if (haveFEData_ == 2) HYPRE_LSI_MLISFEIDestroy(feData_);
607 feData_ = NULL;
608 }
609 if ( MLI_NodalCoord_ != NULL ) delete [] MLI_NodalCoord_;
610 if ( MLI_EqnNumbers_ != NULL ) delete [] MLI_EqnNumbers_;
611 #endif
612
613 if (maxwellANN_ != NULL)
614 {
615 HYPRE_ParCSRMatrixDestroy(maxwellANN_);
616 maxwellANN_ = NULL;
617 }
618 if (amsX_ != NULL) HYPRE_IJVectorDestroy(amsX_);
619 if (amsY_ != NULL) HYPRE_IJVectorDestroy(amsY_);
620 if (amsZ_ != NULL) HYPRE_IJVectorDestroy(amsZ_);
621 if (amsG_ != NULL) HYPRE_IJMatrixDestroy(amsG_);
622 if (amsD0_ != NULL) HYPRE_IJMatrixDestroy(amsD0_);
623 if (amsD1_ != NULL) HYPRE_IJMatrixDestroy(amsD1_);
624 // Users who copy this matrix in should be responsible for
625 // destroying this
626 //if (maxwellGEN_ != NULL)
627 //{
628 // HYPRE_ParCSRMatrixDestroy(maxwellGEN_);
629 // maxwellGEN_ = NULL;
630 //}
631 // This data seems to be freed by hypre_AMSFEIDestroy in ams.c
632 //if (AMSData_.EdgeNodeList_ != NULL) hypre_TFree(AMSData_.EdgeNodeList_, HYPRE_MEMORY_HOST);
633 if (AMSData_.NodeNumbers_ != NULL) hypre_TFree(AMSData_.NodeNumbers_, HYPRE_MEMORY_HOST);
634 if (AMSData_.NodalCoord_ != NULL) hypre_TFree(AMSData_.NodalCoord_, HYPRE_MEMORY_HOST);
635 if (FEI_mixedDiag_ != NULL) delete [] FEI_mixedDiag_;
636
637 //-------------------------------------------------------------------
638 // diagnostic message
639 //-------------------------------------------------------------------
640
641 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
642 printf("%4d : HYPRE_LSC::leaving destructor.\n",mypid_);
643 }
644
645 //***************************************************************************
646 // clone a copy of HYPRE_LinSysCore
647 //---------------------------------------------------------------------------
648
649 #ifndef NOFEI
clone()650 LinearSystemCore* HYPRE_LinSysCore::clone()
651 {
652 return(new HYPRE_LinSysCore(comm_));
653 }
654 #endif
655
656 //***************************************************************************
657 // passing a lookup table to this object
658 // (note : this is called in FEI_initComplete)
659 //---------------------------------------------------------------------------
660
setLookup(Lookup & lookup)661 int HYPRE_LinSysCore::setLookup(Lookup& lookup)
662 {
663 //-------------------------------------------------------------------
664 // diagnostic message
665 //-------------------------------------------------------------------
666
667 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
668 printf("%4d : HYPRE_LSC::entering setLookup.\n",mypid_);
669
670 //-------------------------------------------------------------------
671 // set the lookup object and initialize the MLI_FEData object
672 //-------------------------------------------------------------------
673
674 // RDF: The following line doesn't make sense and generates warnings with some compilers
675 //if (&lookup == NULL) return (0);
676 lookup_ = &lookup;
677 haveLookup_ = 1;
678
679 //-------------------------------------------------------------------
680 // diagnostic message
681 //-------------------------------------------------------------------
682
683 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
684 printf("%4d : HYPRE_LSC::leaving setLookup.\n",mypid_);
685 return (0);
686 }
687
688 //***************************************************************************
689 //This function is where we establish the structures/objects associated
690 //with the linear algebra library. i.e., do initial allocations, etc.
691 // Rows and columns are 1-based.
692 //---------------------------------------------------------------------------
693
createMatricesAndVectors(int numGlobalEqns,int firstLocalEqn,int numLocalEqns)694 int HYPRE_LinSysCore::createMatricesAndVectors(int numGlobalEqns,
695 int firstLocalEqn, int numLocalEqns)
696 {
697 int i;
698
699 //-------------------------------------------------------------------
700 // diagnostic message
701 //-------------------------------------------------------------------
702
703 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
704 {
705 printf("%4d : HYPRE_LSC::entering createMatricesAndVectors.\n",mypid_);
706 printf("%4d : HYPRE_LSC::startrow, endrow = %d %d\n",mypid_,
707 firstLocalEqn, firstLocalEqn+numLocalEqns-1);
708 }
709
710 //-------------------------------------------------------------------
711 // clean up previously allocated matrix
712 //-------------------------------------------------------------------
713
714 if ( rowLengths_ != NULL ) delete [] rowLengths_;
715 if ( colIndices_ != NULL )
716 {
717 int nrows = localEndRow_ - localStartRow_ + 1;
718 for ( i = 0; i < nrows; i++ )
719 if ( colIndices_[i] != NULL ) delete [] colIndices_[i];
720 delete [] colIndices_;
721 }
722 if ( colValues_ != NULL )
723 {
724 int nrows = localEndRow_ - localStartRow_ + 1;
725 for ( i = 0; i < nrows; i++ )
726 if ( colValues_[i] != NULL ) delete [] colValues_[i];
727 delete [] colValues_;
728 }
729 rowLengths_ = NULL;
730 colIndices_ = NULL;
731 colValues_ = NULL;
732
733 if ( mRHSNumGEqns_ > 0)
734 {
735 if (mRHSGEqnIDs_ != NULL) delete [] mRHSGEqnIDs_;
736 if (mRHSNEntries_ != NULL) delete [] mRHSNEntries_;
737 if (mRHSBCType_ != NULL) delete [] mRHSBCType_;
738 if (mRHSRowInds_ != NULL)
739 {
740 for (i = 0; i < mRHSNumGEqns_; i++) delete [] mRHSRowInds_[i];
741 delete [] mRHSRowInds_;
742 }
743 if (mRHSRowVals_ != NULL)
744 {
745 for (i = 0; i < mRHSNumGEqns_; i++) delete [] mRHSRowVals_[i];
746 delete [] mRHSRowVals_;
747 }
748 mRHSNumGEqns_ = 0;
749 mRHSGEqnIDs_ = NULL;
750 mRHSNEntries_ = NULL;
751 mRHSBCType_ = NULL;
752 mRHSRowInds_ = NULL;
753 mRHSRowVals_ = NULL;
754 }
755
756 //-------------------------------------------------------------------
757 // error checking
758 //-------------------------------------------------------------------
759
760 if ( ( firstLocalEqn <= 0 ) ||
761 ( firstLocalEqn+numLocalEqns-1) > numGlobalEqns)
762 {
763 printf("%4d : createMatricesVectors: invalid local equation nos.\n",
764 mypid_);
765 exit(1);
766 }
767
768 localStartRow_ = firstLocalEqn;
769 localEndRow_ = firstLocalEqn + numLocalEqns - 1;
770 numGlobalRows_ = numGlobalEqns;
771
772 //-------------------------------------------------------------------
773 // first clean up previous allocations
774 //-------------------------------------------------------------------
775
776 if ( matrixVectorsCreated_ )
777 {
778 if ( HYA_ != NULL ) {HYPRE_IJMatrixDestroy(HYA_); HYA_ = NULL;}
779 if ( HYx_ != NULL ) {HYPRE_IJVectorDestroy(HYx_); HYx_ = NULL;}
780 if ( HYr_ != NULL ) {HYPRE_IJVectorDestroy(HYr_); HYr_ = NULL;}
781 if ( HYbs_ != NULL )
782 {
783 for ( i = 0; i < numRHSs_; i++ )
784 if ( HYbs_[i] != NULL ) HYPRE_IJVectorDestroy(HYbs_[i]);
785 delete [] HYbs_;
786 HYbs_ = NULL;
787 }
788 if (reducedA_ != NULL) HYPRE_IJMatrixDestroy(reducedA_);
789 if (reducedB_ != NULL) HYPRE_IJVectorDestroy(reducedB_);
790 if (reducedX_ != NULL) HYPRE_IJVectorDestroy(reducedX_);
791 if (reducedR_ != NULL) HYPRE_IJVectorDestroy(reducedR_);
792 if (HYA21_ != NULL) HYPRE_IJMatrixDestroy(HYA21_);
793 if (HYA12_ != NULL) HYPRE_IJMatrixDestroy(HYA12_);
794 if (HYinvA22_ != NULL) HYPRE_IJMatrixDestroy(HYinvA22_);
795 reducedA_ = NULL;
796 reducedB_ = NULL;
797 reducedX_ = NULL;
798 reducedR_ = NULL;
799 HYA21_ = NULL;
800 HYA12_ = NULL;
801 HYinvA22_ = NULL;
802 A21NRows_ = A21NCols_ = reducedAStartRow_ = 0;
803 }
804
805 //-------------------------------------------------------------------
806 // instantiate the matrix (can also handle rectangular matrix)
807 //-------------------------------------------------------------------
808
809 if (localStartCol_ == -1)
810 HYPRE_IJMatrixCreate(comm_, localStartRow_-1,localEndRow_-1,
811 localStartRow_-1,localEndRow_-1, &HYA_);
812 else
813 HYPRE_IJMatrixCreate(comm_, localStartRow_-1,localEndRow_-1,
814 localStartCol_,localEndCol_, &HYA_);
815 HYPRE_IJMatrixSetObjectType(HYA_, HYPRE_PARCSR);
816
817 //-------------------------------------------------------------------
818 // instantiate the right hand vectors
819 //-------------------------------------------------------------------
820
821 HYbs_ = new HYPRE_IJVector[numRHSs_];
822 for ( i = 0; i < numRHSs_; i++ )
823 {
824 HYPRE_IJVectorCreate(comm_, localStartRow_-1, localEndRow_-1,
825 &(HYbs_[i]));
826 HYPRE_IJVectorSetObjectType(HYbs_[i], HYPRE_PARCSR);
827 HYPRE_IJVectorInitialize(HYbs_[i]);
828 HYPRE_IJVectorAssemble(HYbs_[i]);
829 }
830 HYb_ = HYbs_[0];
831
832 //-------------------------------------------------------------------
833 // instantiate the solution vector
834 //-------------------------------------------------------------------
835
836 if (localStartCol_ == -1)
837 HYPRE_IJVectorCreate(comm_,localStartRow_-1,localEndRow_-1,&HYx_);
838 else
839 HYPRE_IJVectorCreate(comm_,localStartCol_,localEndCol_,&HYx_);
840 HYPRE_IJVectorSetObjectType(HYx_, HYPRE_PARCSR);
841 HYPRE_IJVectorInitialize(HYx_);
842 HYPRE_IJVectorAssemble(HYx_);
843
844 //-------------------------------------------------------------------
845 // reset fedata
846 //-------------------------------------------------------------------
847
848 #ifdef HAVE_MLI
849 if ( feData_ != NULL )
850 {
851 if (haveFEData_ == 1) HYPRE_LSI_MLIFEDataDestroy(feData_);
852 else if (haveFEData_ == 2) HYPRE_LSI_MLISFEIDestroy(feData_);
853 feData_ = NULL;
854 if ( MLI_NodalCoord_ != NULL ) delete [] MLI_NodalCoord_;
855 if ( MLI_EqnNumbers_ != NULL ) delete [] MLI_EqnNumbers_;
856 MLI_NodalCoord_ = NULL;
857 MLI_EqnNumbers_ = NULL;
858 MLI_NumNodes_ = 0;
859 }
860 #endif
861
862 //-------------------------------------------------------------------
863 // for amge
864 //-------------------------------------------------------------------
865
866 #ifdef HAVE_AMGE
867 HYPRE_LSI_AMGeCreate();
868 HYPRE_LSI_AMGeSetNNodes(numGlobalRows_);
869 HYPRE_LSI_AMGeSetNElements(numGlobalRows_);
870 HYPRE_LSI_AMGeSetSystemSize(1);
871 #endif
872
873 //-------------------------------------------------------------------
874 // instantiate the residual vector
875 //-------------------------------------------------------------------
876
877 HYPRE_IJVectorCreate(comm_, localStartRow_-1, localEndRow_-1, &HYr_);
878 HYPRE_IJVectorSetObjectType(HYr_, HYPRE_PARCSR);
879 HYPRE_IJVectorInitialize(HYr_);
880 HYPRE_IJVectorAssemble(HYr_);
881 matrixVectorsCreated_ = 1;
882 schurReductionCreated_ = 0;
883 systemAssembled_ = 0;
884 normalEqnFlag_ &= 1;
885 if ( HYnormalA_ != NULL )
886 {
887 HYPRE_IJMatrixDestroy(HYnormalA_);
888 HYnormalA_ = NULL;
889 }
890 if ( HYnormalB_ != NULL )
891 {
892 HYPRE_IJVectorDestroy(HYnormalB_);
893 HYnormalB_ = NULL;
894 }
895
896 //-------------------------------------------------------------------
897 // diagnostic message
898 //-------------------------------------------------------------------
899
900 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
901 printf("%4d : HYPRE_LSC::leaving createMatricesAndVectors.\n",mypid_);
902 return (0);
903 }
904
905 //***************************************************************************
906 // set global and local number of equations
907 // (This is called in FEI_initComplete after setLookup)
908 //---------------------------------------------------------------------------
909
setGlobalOffsets(int leng,int * nodeOffsets,int * eqnOffsets,int * blkEqnOffsets)910 int HYPRE_LinSysCore::setGlobalOffsets(int leng, int* nodeOffsets,
911 int* eqnOffsets, int* blkEqnOffsets)
912 {
913 //-------------------------------------------------------------------
914 // diagnostic message
915 //-------------------------------------------------------------------
916
917 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
918 printf("%4d : HYPRE_LSC::entering setGlobalOffsets.\n",mypid_);
919
920 //-------------------------------------------------------------------
921 // set local range (incoming 0-based, locally 1-based)
922 //-------------------------------------------------------------------
923
924 (void) leng;
925 (void) nodeOffsets;
926 (void) blkEqnOffsets;
927 int firstLocalEqn = eqnOffsets[mypid_] + 1;
928 int numLocalEqns = eqnOffsets[mypid_+1] - firstLocalEqn + 1;
929 int numGlobalEqns = eqnOffsets[numProcs_];
930 createMatricesAndVectors(numGlobalEqns, firstLocalEqn, numLocalEqns);
931
932 //-------------------------------------------------------------------
933 // diagnostic message
934 //-------------------------------------------------------------------
935
936 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
937 {
938 printf("%4d : HYPRE_LSC::startrow, endrow = %d %d\n",mypid_,
939 localStartRow_, localEndRow_);
940 printf("%4d : HYPRE_LSC::leaving setGlobalOffsets.\n",mypid_);
941 }
942 return (0);
943 }
944
945 //***************************************************************************
946 // Grid related function : element node connectivities
947 //---------------------------------------------------------------------------
948
setConnectivities(GlobalID elemBlk,int nElems,int nNodesPerElem,const GlobalID * elemIDs,const int * const * connNodes)949 int HYPRE_LinSysCore::setConnectivities(GlobalID elemBlk, int nElems,
950 int nNodesPerElem, const GlobalID* elemIDs,
951 const int* const* connNodes)
952 {
953 #ifdef HAVE_MLI
954 (void) elemIDs;
955 (void) connNodes;
956 if ( HYPreconID_ == HYMLI && haveFEData_ == 2 )
957 {
958 if (feData_ == NULL) feData_ = (void *) HYPRE_LSI_MLISFEICreate(comm_);
959 HYPRE_LSI_MLISFEIAddNumElems(feData_,elemBlk,nElems,nNodesPerElem);
960 }
961 #else
962 (void) elemBlk;
963 (void) nElems;
964 (void) nNodesPerElem;
965 (void) elemIDs;
966 (void) connNodes;
967 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) > 3 )
968 printf("%4d : HYPRE_LSC::setConnectivities not implemented.\n",mypid_);
969 #endif
970 return (0);
971 }
972
973 //***************************************************************************
974 // Grid related function : element stiffness matrix loading
975 //---------------------------------------------------------------------------
976
setStiffnessMatrices(GlobalID elemBlk,int nElems,const GlobalID * elemIDs,const double * const * const * stiff,int nEqnsPerElem,const int * const * eqnIndices)977 int HYPRE_LinSysCore::setStiffnessMatrices(GlobalID elemBlk, int nElems,
978 const GlobalID* elemIDs,const double *const *const *stiff,
979 int nEqnsPerElem, const int *const * eqnIndices)
980 {
981 #ifdef HAVE_MLI
982 if ( HYPreconID_ == HYMLI && feData_ != NULL )
983 {
984 HYPRE_LSI_MLISFEILoadElemMatrices(feData_,elemBlk,nElems,(int*)elemIDs,
985 (double***)stiff,nEqnsPerElem,(int**)eqnIndices);
986 }
987 #else
988 (void) elemBlk;
989 (void) nElems;
990 (void) elemIDs;
991 (void) stiff;
992 (void) nEqnsPerElem;
993 (void) eqnIndices;
994 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) > 3 )
995 printf("%4d : HYPRE_LSC::setStiffnessMatrices not implemented.\n",
996 mypid_);
997 #endif
998 return (0);
999 }
1000
1001 //***************************************************************************
1002 // Grid related function : element load vector loading
1003 //---------------------------------------------------------------------------
1004
setLoadVectors(GlobalID elemBlock,int numElems,const GlobalID * elemIDs,const double * const * load,int numEqnsPerElem,const int * const * eqnIndices)1005 int HYPRE_LinSysCore::setLoadVectors(GlobalID elemBlock, int numElems,
1006 const GlobalID* elemIDs, const double *const *load,
1007 int numEqnsPerElem, const int *const * eqnIndices)
1008 {
1009 (void) elemBlock;
1010 (void) numElems;
1011 (void) elemIDs;
1012 (void) load;
1013 (void) numEqnsPerElem;
1014 (void) eqnIndices;
1015
1016 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) > 3 )
1017 printf("%4d : HYPRE_LSC::setLoadVectors not implemented.\n", mypid_);
1018 return (0);
1019 }
1020
1021 //***************************************************************************
1022 // Set the number of rows in the diagonal part and off diagonal part
1023 // of the matrix, using the structure of the matrix, stored in rows.
1024 // rows is an array that is 0-based. localStartRow and localEndRow are 1-based.
1025 //---------------------------------------------------------------------------
1026
allocateMatrix(int ** colIndices,int * rowLengths)1027 int HYPRE_LinSysCore::allocateMatrix(int **colIndices, int *rowLengths)
1028 {
1029 int i, j, nsize, rowLeng, maxSize, minSize, searchFlag, *indPtr, *indPtr2;
1030 double *vals;
1031
1032 //-------------------------------------------------------------------
1033 // diagnoistic message and error checking
1034 //-------------------------------------------------------------------
1035
1036 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1037 {
1038 printf("%4d : HYPRE_LSC::entering allocateMatrix.\n", mypid_);
1039 if ( localEndRow_ < localStartRow_ )
1040 {
1041 printf("allocateMatrix WARNING : createMatrixAndVectors should be\n");
1042 printf(" called before allocateMatrix.\n");
1043 }
1044 }
1045
1046 //-------------------------------------------------------------------
1047 // clean up previously allocated matrix
1048 //-------------------------------------------------------------------
1049
1050 if ( rowLengths_ != NULL ) delete [] rowLengths_;
1051 rowLengths_ = NULL;
1052 if ( colIndices_ != NULL )
1053 {
1054 for ( i = 0; i < localEndRow_-localStartRow_+1; i++ )
1055 if ( colIndices_[i] != NULL ) delete [] colIndices_[i];
1056 delete [] colIndices_;
1057 colIndices_ = NULL;
1058 }
1059 if ( colValues_ != NULL )
1060 {
1061 for ( i = 0; i < localEndRow_-localStartRow_+1; i++ )
1062 if ( colValues_[i] != NULL ) delete [] colValues_[i];
1063 delete [] colValues_;
1064 colValues_ = NULL;
1065 }
1066
1067 //-------------------------------------------------------------------
1068 // allocate and store the column index information
1069 //-------------------------------------------------------------------
1070
1071 nsize = localEndRow_ - localStartRow_ + 1;
1072 rowLengths_ = new int[nsize];
1073 colIndices_ = new int*[nsize];
1074 colValues_ = new double*[nsize];
1075 maxSize = 0;
1076 minSize = 1000000;
1077 for ( i = 0; i < nsize; i++ )
1078 {
1079 rowLeng = rowLengths_[i] = rowLengths[i];
1080 if ( rowLeng > 0 )
1081 {
1082 colIndices_[i] = new int[rowLeng];
1083 hypre_assert( colIndices_[i] != NULL );
1084 }
1085 else colIndices_[i] = NULL;
1086 indPtr = colIndices_[i];
1087 indPtr2 = colIndices[i];
1088 for ( j = 0; j < rowLeng; j++ ) indPtr[j] = indPtr2[j];
1089 searchFlag = 0;
1090 for ( j = 1; j < rowLeng; j++ )
1091 if ( indPtr[j] < indPtr[j-1]) {searchFlag = 1; break;}
1092 if ( searchFlag ) hypre_qsort0( indPtr, 0, rowLeng-1);
1093 maxSize = ( rowLeng > maxSize ) ? rowLeng : maxSize;
1094 minSize = ( rowLeng < minSize ) ? rowLeng : minSize;
1095 if ( rowLeng > 0 )
1096 {
1097 colValues_[i] = new double[rowLeng];
1098 hypre_assert( colValues_[i] != NULL );
1099 }
1100 vals = colValues_[i];
1101 for ( j = 0; j < rowLeng; j++ ) vals[j] = 0.0;
1102 }
1103 MPI_Allreduce(&maxSize, &pilutMaxNnzPerRow_,1,MPI_INT,MPI_MAX,comm_);
1104
1105 //-------------------------------------------------------------------
1106 // diagnoistic message
1107 //-------------------------------------------------------------------
1108
1109 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1110 {
1111 printf("%4d : allocateMatrix : max/min nnz/row = %d %d\n", mypid_,
1112 maxSize, minSize);
1113 printf("%4d : HYPRE_LSC::leaving allocateMatrix.\n", mypid_);
1114 }
1115 return (0);
1116 }
1117
1118 //***************************************************************************
1119 // to establish the structures/objects associated with the linear algebra
1120 // library. i.e., do initial allocations, etc.
1121 //---------------------------------------------------------------------------
1122
setMatrixStructure(int ** ptColIndices,int * ptRowLengths,int ** blkColIndices,int * blkRowLengths,int * ptRowsPerBlkRow)1123 int HYPRE_LinSysCore::setMatrixStructure(int** ptColIndices, int* ptRowLengths,
1124 int** blkColIndices, int* blkRowLengths,
1125 int* ptRowsPerBlkRow)
1126 {
1127 int i, j;
1128
1129 //-------------------------------------------------------------------
1130 // diagnoistic message
1131 //-------------------------------------------------------------------
1132
1133 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1134 {
1135 printf("%4d : HYPRE_LSC::entering setMatrixStructure.\n",mypid_);
1136 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 6 )
1137 {
1138 int nRows = localEndRow_ - localStartRow_ + 1;
1139 for (i = 0; i < nRows; i++)
1140 for (j = 0; j < ptRowLengths[i]; j++)
1141 printf(" %4d : row, col = %d %d\n",mypid_,
1142 localStartRow_+i, ptColIndices[i][j]+1);
1143 }
1144 }
1145 (void) blkColIndices;
1146 (void) blkRowLengths;
1147 (void) ptRowsPerBlkRow;
1148
1149 //-------------------------------------------------------------------
1150 // allocate local space for matrix
1151 //-------------------------------------------------------------------
1152
1153 int numLocalRows = localEndRow_ - localStartRow_ + 1;
1154 for ( i = 0; i < numLocalRows; i++ )
1155 for ( j = 0; j < ptRowLengths[i]; j++ ) ptColIndices[i][j]++;
1156
1157 allocateMatrix(ptColIndices, ptRowLengths);
1158
1159 for ( i = 0; i < numLocalRows; i++ )
1160 for ( j = 0; j < ptRowLengths[i]; j++ ) ptColIndices[i][j]--;
1161
1162 //-------------------------------------------------------------------
1163 // diagnoistic message
1164 //-------------------------------------------------------------------
1165
1166 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1167 printf("%4d : HYPRE_LSC::leaving setMatrixStructure.\n",mypid_);
1168 return (0);
1169 }
1170
1171 //***************************************************************************
1172 // set Lagrange multiplier equations
1173 //---------------------------------------------------------------------------
1174
setMultCREqns(int multCRSetID,int numCRs,int numNodesPerCR,int ** nodeNumbers,int ** eqnNumbers,int * fieldIDs,int * multiplierEqnNumbers)1175 int HYPRE_LinSysCore::setMultCREqns(int multCRSetID, int numCRs,
1176 int numNodesPerCR, int** nodeNumbers, int** eqnNumbers,
1177 int* fieldIDs, int* multiplierEqnNumbers)
1178 {
1179 (void) multCRSetID;
1180 (void) numCRs;
1181 (void) numNodesPerCR;
1182 (void) nodeNumbers;
1183 (void) eqnNumbers;
1184 (void) fieldIDs;
1185 (void) multiplierEqnNumbers;
1186
1187 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) > 3 )
1188 printf("%4d : HYPRE_LSC::setMultCREqns not implemented.\n",mypid_);
1189 return (0);
1190 }
1191
1192 //***************************************************************************
1193 // set penalty constraint equations
1194 //---------------------------------------------------------------------------
1195
setPenCREqns(int penCRSetID,int numCRs,int numNodesPerCR,int ** nodeNumbers,int ** eqnNumbers,int * fieldIDs)1196 int HYPRE_LinSysCore::setPenCREqns(int penCRSetID, int numCRs,
1197 int numNodesPerCR, int** nodeNumbers, int** eqnNumbers,
1198 int* fieldIDs)
1199 {
1200 (void) penCRSetID;
1201 (void) numCRs;
1202 (void) numNodesPerCR;
1203 (void) nodeNumbers;
1204 (void) eqnNumbers;
1205 (void) fieldIDs;
1206
1207 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) > 3 )
1208 printf("%4d : HYPRE_LSC::setPenCREqns not implemented.\n",mypid_);
1209 return (0);
1210 }
1211
1212 //***************************************************************************
1213 // This function is needed in order to construct a new problem with the
1214 // same sparsity pattern.
1215 //---------------------------------------------------------------------------
1216
resetMatrixAndVector(double setValue)1217 int HYPRE_LinSysCore::resetMatrixAndVector(double setValue)
1218 {
1219 int i, j, localNRows, *cols;
1220 double *vals;
1221
1222 //-------------------------------------------------------------------
1223 // diagnoistic message
1224 //-------------------------------------------------------------------
1225
1226 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1227 printf("%4d : HYPRE_LSC::entering resetMatrixAndVector.\n",mypid_);
1228 if ( setValue != 0.0 && mypid_ == 0 )
1229 {
1230 printf("resetMatrixAndVector ERROR : cannot take nonzeros.\n");
1231 exit(1);
1232 }
1233
1234 //-------------------------------------------------------------------
1235 // reset vector values
1236 //-------------------------------------------------------------------
1237
1238 localNRows = localEndRow_ - localStartRow_ + 1;
1239 cols = new int[localNRows];
1240 vals = new double[localNRows];
1241 for (i = 0; i < localNRows; i++)
1242 {
1243 cols[i] = localStartRow_ + i - 1;
1244 vals[i] = 0.0;
1245 }
1246 for (i = 0; i < numRHSs_; i++)
1247 HYPRE_IJVectorSetValues(HYbs_[i], localNRows, (const int *) cols,
1248 (const double *) vals);
1249
1250 delete [] cols;
1251 delete [] vals;
1252
1253 systemAssembled_ = 0;
1254 schurReductionCreated_ = 0;
1255 projectCurrSize_ = 0;
1256 normalEqnFlag_ &= 1;
1257 if ( HYnormalA_ != NULL )
1258 {
1259 HYPRE_IJMatrixDestroy(HYnormalA_);
1260 HYnormalA_ = NULL;
1261 }
1262 if ( HYnormalB_ != NULL )
1263 {
1264 HYPRE_IJVectorDestroy(HYnormalB_);
1265 HYnormalB_ = NULL;
1266 }
1267
1268 //-------------------------------------------------------------------
1269 // for now, since HYPRE does not yet support
1270 // re-initializing the matrix, restart the whole thing
1271 //-------------------------------------------------------------------
1272
1273 if ( HYA_ != NULL ) HYPRE_IJMatrixDestroy(HYA_);
1274 HYPRE_IJMatrixCreate(comm_, localStartRow_-1, localEndRow_-1,
1275 localStartRow_-1, localEndRow_-1, &HYA_);
1276 HYPRE_IJMatrixSetObjectType(HYA_, HYPRE_PARCSR);
1277
1278 //-------------------------------------------------------------------
1279 // clean the reduction stuff
1280 //-------------------------------------------------------------------
1281
1282 if (reducedA_ != NULL) {HYPRE_IJMatrixDestroy(reducedA_); reducedA_ = NULL;}
1283 if (reducedB_ != NULL) {HYPRE_IJVectorDestroy(reducedB_); reducedB_ = NULL;}
1284 if (reducedX_ != NULL) {HYPRE_IJVectorDestroy(reducedX_); reducedX_ = NULL;}
1285 if (reducedR_ != NULL) {HYPRE_IJVectorDestroy(reducedR_); reducedR_ = NULL;}
1286 if (HYA21_ != NULL) {HYPRE_IJMatrixDestroy(HYA21_); HYA21_ = NULL;}
1287 if (HYA12_ != NULL) {HYPRE_IJMatrixDestroy(HYA12_); HYA12_ = NULL;}
1288 if (HYinvA22_ != NULL) {HYPRE_IJMatrixDestroy(HYinvA22_); HYinvA22_ = NULL;}
1289 A21NRows_ = A21NCols_ = reducedAStartRow_ = 0;
1290
1291 //-------------------------------------------------------------------
1292 // allocate space for storing the matrix coefficient
1293 //-------------------------------------------------------------------
1294
1295 if ( colValues_ != NULL )
1296 {
1297 int nrows = localEndRow_ - localStartRow_ + 1;
1298 for ( i = 0; i < nrows; i++ )
1299 if ( colValues_[i] != NULL ) delete [] colValues_[i];
1300 delete [] colValues_;
1301 }
1302 colValues_ = NULL;
1303
1304 colValues_ = new double*[localNRows];
1305 for ( i = 0; i < localNRows; i++ )
1306 {
1307 if ( rowLengths_[i] > 0 ) colValues_[i] = new double[rowLengths_[i]];
1308 for ( j = 0; j < rowLengths_[i]; j++ ) colValues_[i][j] = 0.0;
1309 }
1310
1311 //-------------------------------------------------------------------
1312 // reset fedata
1313 //-------------------------------------------------------------------
1314
1315 #ifdef HAVE_MLI
1316 if ( feData_ != NULL )
1317 {
1318 if (haveFEData_ == 1) HYPRE_LSI_MLIFEDataDestroy(feData_);
1319 else if (haveFEData_ == 2) HYPRE_LSI_MLISFEIDestroy(feData_);
1320 feData_ = NULL;
1321 if ( MLI_NodalCoord_ != NULL ) delete [] MLI_NodalCoord_;
1322 if ( MLI_EqnNumbers_ != NULL ) delete [] MLI_EqnNumbers_;
1323 MLI_NodalCoord_ = NULL;
1324 MLI_EqnNumbers_ = NULL;
1325 MLI_NumNodes_ = 0;
1326 }
1327 #endif
1328
1329 //-------------------------------------------------------------------
1330 // diagnostic message
1331 //-------------------------------------------------------------------
1332
1333 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1334 printf("%4d : HYPRE_LSC::leaving resetMatrixAndVector.\n", mypid_);
1335 return (0);
1336 }
1337
1338 //***************************************************************************
1339 // new function to reset matrix independently
1340 //---------------------------------------------------------------------------
1341
resetMatrix(double setValue)1342 int HYPRE_LinSysCore::resetMatrix(double setValue)
1343 {
1344 int i, j, size;
1345
1346 //-------------------------------------------------------------------
1347 // diagnostic message
1348 //-------------------------------------------------------------------
1349
1350 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1351 printf("%4d : HYPRE_LSC::entering resetMatrix.\n",mypid_);
1352 if ( setValue != 0.0 && mypid_ == 0 )
1353 {
1354 printf("resetMatrix ERROR : cannot take nonzeros.\n");
1355 exit(1);
1356 }
1357
1358 //-------------------------------------------------------------------
1359 // clean the reduction stuff
1360 //-------------------------------------------------------------------
1361
1362 if (reducedA_ != NULL) {HYPRE_IJMatrixDestroy(reducedA_); reducedA_ = NULL;}
1363 if (reducedB_ != NULL) {HYPRE_IJVectorDestroy(reducedB_); reducedB_ = NULL;}
1364 if (reducedX_ != NULL) {HYPRE_IJVectorDestroy(reducedX_); reducedX_ = NULL;}
1365 if (reducedR_ != NULL) {HYPRE_IJVectorDestroy(reducedR_); reducedR_ = NULL;}
1366 if (HYA21_ != NULL) {HYPRE_IJMatrixDestroy(HYA21_); HYA21_ = NULL;}
1367 if (HYA12_ != NULL) {HYPRE_IJMatrixDestroy(HYA12_); HYA12_ = NULL;}
1368 if (HYinvA22_ != NULL) {HYPRE_IJMatrixDestroy(HYinvA22_); HYinvA22_ = NULL;}
1369 A21NRows_ = A21NCols_ = reducedAStartRow_ = 0;
1370
1371 //-------------------------------------------------------------------
1372 // for now, since HYPRE does not yet support
1373 // re-initializing the matrix, restart the whole thing
1374 //-------------------------------------------------------------------
1375
1376 if ( HYA_ != NULL ) HYPRE_IJMatrixDestroy(HYA_);
1377 size = localEndRow_ - localStartRow_ + 1;
1378 if (localStartCol_ == -1)
1379 HYPRE_IJMatrixCreate(comm_, localStartRow_-1, localEndRow_-1,
1380 localStartRow_-1, localEndRow_-1, &HYA_);
1381 else
1382 HYPRE_IJMatrixCreate(comm_, localStartRow_-1, localEndRow_-1,
1383 localStartCol_, localEndCol_, &HYA_);
1384 HYPRE_IJMatrixSetObjectType(HYA_, HYPRE_PARCSR);
1385
1386 //-------------------------------------------------------------------
1387 // allocate space for storing the matrix coefficient
1388 //-------------------------------------------------------------------
1389
1390 if ( colValues_ != NULL )
1391 {
1392 int nrows = localEndRow_ - localStartRow_ + 1;
1393 for ( i = 0; i < nrows; i++ )
1394 if ( colValues_[i] != NULL ) delete [] colValues_[i];
1395 delete [] colValues_;
1396 }
1397 colValues_ = NULL;
1398
1399 colValues_ = new double*[size];
1400 for ( i = 0; i < size; i++ )
1401 {
1402 if ( rowLengths_[i] > 0 ) colValues_[i] = new double[rowLengths_[i]];
1403 for ( j = 0; j < rowLengths_[i]; j++ ) colValues_[i][j] = 0.0;
1404 }
1405
1406 //-------------------------------------------------------------------
1407 // reset system flags
1408 //-------------------------------------------------------------------
1409
1410 systemAssembled_ = 0;
1411 schurReductionCreated_ = 0;
1412 projectCurrSize_ = 0;
1413 normalEqnFlag_ &= 5;
1414 if ( HYnormalA_ != NULL )
1415 {
1416 HYPRE_IJMatrixDestroy(HYnormalA_);
1417 HYnormalA_ = NULL;
1418 }
1419
1420 //-------------------------------------------------------------------
1421 // reset fedata
1422 //-------------------------------------------------------------------
1423
1424 #ifdef HAVE_MLI
1425 if ( feData_ != NULL )
1426 {
1427 if (haveFEData_ == 1) HYPRE_LSI_MLIFEDataDestroy(feData_);
1428 else if (haveFEData_ == 2) HYPRE_LSI_MLISFEIDestroy(feData_);
1429 feData_ = NULL;
1430 if ( MLI_NodalCoord_ != NULL ) delete [] MLI_NodalCoord_;
1431 if ( MLI_EqnNumbers_ != NULL ) delete [] MLI_EqnNumbers_;
1432 MLI_NodalCoord_ = NULL;
1433 MLI_EqnNumbers_ = NULL;
1434 MLI_NumNodes_ = 0;
1435 }
1436 #endif
1437
1438 //-------------------------------------------------------------------
1439 // diagnostic message
1440 //-------------------------------------------------------------------
1441
1442 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1443 printf("%4d : HYPRE_LSC::leaving resetMatrix.\n", mypid_);
1444 return (0);
1445 }
1446
1447 //***************************************************************************
1448 // new function to reset vectors independently
1449 //---------------------------------------------------------------------------
1450
resetRHSVector(double setValue)1451 int HYPRE_LinSysCore::resetRHSVector(double setValue)
1452 {
1453 int i, localNRows, *cols;
1454 double *vals;
1455
1456 //-------------------------------------------------------------------
1457 // diagnostic message
1458 //-------------------------------------------------------------------
1459
1460 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1461 printf("%4d : HYPRE_LSC::entering resetRHSVector.\n",mypid_);
1462
1463 //-------------------------------------------------------------------
1464 // reset right hand side vectors
1465 //-------------------------------------------------------------------
1466
1467 if ( HYbs_ != NULL )
1468 {
1469 localNRows = localEndRow_ - localStartRow_ + 1;
1470 cols = new int[localNRows];
1471 vals = new double[localNRows];
1472 for (i = 0; i < localNRows; i++)
1473 {
1474 cols[i] = localStartRow_ + i - 1;
1475 vals[i] = setValue;
1476 }
1477 for (i = 0; i < numRHSs_; i++)
1478 if ( HYbs_[i] != NULL )
1479 HYPRE_IJVectorSetValues(HYbs_[i], localNRows, (const int *) cols,
1480 (const double *) vals);
1481 delete [] cols;
1482 delete [] vals;
1483 }
1484 normalEqnFlag_ &= 3;
1485 if ( HYnormalB_ != NULL )
1486 {
1487 HYPRE_IJVectorDestroy(HYnormalB_);
1488 HYnormalB_ = NULL;
1489 }
1490
1491 //-------------------------------------------------------------------
1492 // diagnostic message
1493 //-------------------------------------------------------------------
1494
1495 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1496 printf("%4d : HYPRE_LSC::leaving resetRHSVector.\n",mypid_);
1497 return (0);
1498 }
1499
1500 //***************************************************************************
1501 // add nonzero entries into the matrix data structure (not in LSC but needed)
1502 //---------------------------------------------------------------------------
1503
sumIntoSystemMatrix(int row,int numValues,const double * values,const int * scatterIndices)1504 int HYPRE_LinSysCore::sumIntoSystemMatrix(int row, int numValues,
1505 const double* values, const int* scatterIndices)
1506 {
1507 int i, j, index, colIndex, localRow;
1508
1509 //-------------------------------------------------------------------
1510 // diagnostic message for high output level only
1511 //-------------------------------------------------------------------
1512
1513 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
1514 {
1515 printf("%4d : HYPRE_LSC::entering sumIntoSystemMatrix.\n",mypid_);
1516 printf("%4d : row number = %d.\n", mypid_, row);
1517 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 6 )
1518 for ( i = 0; i < numValues; i++ )
1519 printf(" %4d : row,col = %d %d, data = %e\n", mypid_,
1520 row+1, scatterIndices[i]+1, values[i]);
1521 }
1522
1523 //-------------------------------------------------------------------
1524 // error checking
1525 //-------------------------------------------------------------------
1526
1527 if ( systemAssembled_ == 1 )
1528 {
1529 printf("%4d : sumIntoSystemMatrix ERROR : matrix already assembled\n",
1530 mypid_);
1531 exit(1);
1532 }
1533 if ( row < localStartRow_ || row > localEndRow_ )
1534 {
1535 printf("%4d : sumIntoSystemMatrix ERROR : invalid row number %d.\n",
1536 mypid_,row);
1537 exit(1);
1538 }
1539 localRow = row - localStartRow_;
1540 if ( numValues > rowLengths_[localRow] )
1541 {
1542 printf("%4d : sumIntoSystemMatrix ERROR : row size too large.\n",mypid_);
1543 exit(1);
1544 }
1545
1546 //-------------------------------------------------------------------
1547 // load the local matrix
1548 //-------------------------------------------------------------------
1549
1550 for ( i = 0; i < numValues; i++ )
1551 {
1552 colIndex = scatterIndices[i];
1553 index = hypre_BinarySearch(colIndices_[localRow], colIndex,
1554 rowLengths_[localRow]);
1555 if ( index < 0 )
1556 {
1557 printf("%4d : sumIntoSystemMatrix ERROR - loading column",mypid_);
1558 printf(" that has not been declared before - %d.\n",colIndex);
1559 for ( j = 0; j < rowLengths_[localRow]; j++ )
1560 printf(" available column index = %d\n",
1561 colIndices_[localRow][j]);
1562 exit(1);
1563 }
1564 colValues_[localRow][index] += values[i];
1565 }
1566
1567 #ifdef HAVE_AMGE
1568 HYPRE_LSI_AMGePutRow(row,numValues,(double*) values,(int*)scatterIndices);
1569 #endif
1570
1571 //-------------------------------------------------------------------
1572 // diagnostic message
1573 //-------------------------------------------------------------------
1574
1575 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
1576 printf("%4d : HYPRE_LSC::leaving sumIntoSystemMatrix.\n",mypid_);
1577 return (0);
1578 }
1579
1580 //***************************************************************************
1581 // add nonzero entries into the matrix data structure
1582 //---------------------------------------------------------------------------
1583
sumIntoSystemMatrix(int numPtRows,const int * ptRows,int numPtCols,const int * ptCols,const double * const * values)1584 int HYPRE_LinSysCore::sumIntoSystemMatrix(int numPtRows, const int* ptRows,
1585 int numPtCols, const int* ptCols,
1586 const double* const* values)
1587 {
1588 int i, j, k, index, colIndex, localRow, orderFlag=0;
1589 int *indptr, rowLeng, useOld;
1590 double *valptr, *auxValues;
1591
1592 //-------------------------------------------------------------------
1593 // diagnostic message for high output level only
1594 //-------------------------------------------------------------------
1595
1596 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
1597 {
1598 printf("%4d : HYPRE_LSC::entering sumIntoSystemMatrix(2).\n",mypid_);
1599 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 6 )
1600 {
1601 for ( i = 0; i < numPtRows; i++ )
1602 {
1603 localRow = ptRows[i] - localStartRow_ + 1;
1604 for ( j = 0; j < numPtCols; j++ )
1605 printf(" %4d : row,col,val = %8d %8d %e\n",mypid_,
1606 ptRows[i]+1, ptCols[j]+1, values[i][j]);
1607 }
1608 }
1609 }
1610
1611 //-------------------------------------------------------------------
1612 // error checking
1613 //-------------------------------------------------------------------
1614
1615 if ( systemAssembled_ == 1 )
1616 {
1617 printf("sumIntoSystemMatrix ERROR : matrix already assembled\n");
1618 exit(1);
1619 }
1620 if (FEI_mixedDiagFlag_ && FEI_mixedDiag_ == NULL)
1621 {
1622 FEI_mixedDiag_ = new double[localEndRow_-localStartRow_+1];
1623 for ( i = 0; i < localEndRow_-localStartRow_+1; i++ )
1624 FEI_mixedDiag_[i] = 0.0;
1625 }
1626
1627 //-------------------------------------------------------------------
1628 // load the local matrix
1629 //-------------------------------------------------------------------
1630
1631 useOld = orderFlag = 0;
1632 if ( numPtCols == nStored_ && storedIndices_ != NULL )
1633 {
1634 for ( i = 0; i < numPtCols; i++ )
1635 if ( storedIndices_[i] != ptCols[i] ) break;
1636 if ( i == numPtCols ) useOld = 1;
1637 }
1638 if ( ! useOld )
1639 {
1640 for ( i = 1; i < numPtCols; i++ )
1641 if ( ptCols[i] < ptCols[i-1] ) { orderFlag = 1; break; }
1642 if ( orderFlag == 1 )
1643 {
1644 if ( numPtCols != nStored_ )
1645 {
1646 if ( storedIndices_ != NULL ) delete [] storedIndices_;
1647 if ( auxStoredIndices_ != NULL ) delete [] auxStoredIndices_;
1648 storedIndices_ = new int[numPtCols];
1649 auxStoredIndices_ = new int[numPtCols];
1650 nStored_ = numPtCols;
1651 }
1652 for ( i = 0; i < numPtCols; i++ )
1653 {
1654 storedIndices_[i] = ptCols[i];
1655 auxStoredIndices_[i] = i;
1656 }
1657 HYPRE_LSI_qsort1a(storedIndices_,auxStoredIndices_,0,numPtCols-1);
1658 for ( i = 0; i < numPtCols; i++ ) storedIndices_[i] = ptCols[i];
1659 }
1660 else
1661 {
1662 if ( storedIndices_ != NULL ) delete [] storedIndices_;
1663 if ( auxStoredIndices_ != NULL ) delete [] auxStoredIndices_;
1664 storedIndices_ = NULL;
1665 auxStoredIndices_ = NULL;
1666 nStored_ = 0;
1667 }
1668 }
1669 for ( i = 0; i < numPtRows; i++ )
1670 {
1671 localRow = ptRows[i] - localStartRow_ + 1;
1672 indptr = colIndices_[localRow];
1673 valptr = colValues_[localRow];
1674 rowLeng = rowLengths_[localRow];
1675 auxValues = (double *) values[i];
1676 index = 0;
1677 for ( j = 0; j < numPtCols; j++ )
1678 {
1679 if ( storedIndices_ )
1680 colIndex = storedIndices_[auxStoredIndices_[j]] + 1;
1681 else
1682 colIndex = ptCols[j] + 1;
1683
1684 if (FEI_mixedDiag_ && ptRows[i] == ptCols[j] && numPtRows > 1)
1685 FEI_mixedDiag_[ptCols[numPtCols-1]-localStartRow_+1] += auxValues[j];
1686
1687 while ( index < rowLeng && indptr[index] < colIndex ) index++;
1688 if ( index >= rowLeng )
1689 {
1690 printf("%4d : sumIntoSystemMatrix ERROR - loading column",mypid_);
1691 printf(" that has not been declared before - %d (row=%d).\n",
1692 colIndex, ptRows[i]+1);
1693 for ( k = 0; k < rowLeng; k++ )
1694 printf(" available column index = %d\n", indptr[k]);
1695 exit(1);
1696 }
1697 if ( auxStoredIndices_ )
1698 valptr[index] += auxValues[auxStoredIndices_[j]];
1699 else
1700 valptr[index] += auxValues[j];
1701 }
1702 }
1703 #ifdef HAVE_AMGE
1704 HYPRE_LSI_AMGePutRow(localRow,numPtCols,(double*) values[i],(int*)ptCols);
1705 #endif
1706
1707 //-------------------------------------------------------------------
1708 // diagnostic message
1709 //-------------------------------------------------------------------
1710
1711 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
1712 printf("%4d : HYPRE_LSC::leaving sumIntoSystemMatrix(2).\n",mypid_);
1713 return (0);
1714 }
1715
1716 //***************************************************************************
1717 // add nonzero entries into the matrix data structure
1718 //---------------------------------------------------------------------------
1719
sumIntoSystemMatrix(int numPtRows,const int * ptRows,int numPtCols,const int * ptCols,int numBlkRows,const int * blkRows,int numBlkCols,const int * blkCols,const double * const * values)1720 int HYPRE_LinSysCore::sumIntoSystemMatrix(int numPtRows, const int* ptRows,
1721 int numPtCols, const int* ptCols, int numBlkRows,
1722 const int* blkRows, int numBlkCols, const int* blkCols,
1723 const double* const* values)
1724 {
1725 (void) numBlkRows;
1726 (void) blkRows;
1727 (void) numBlkCols;
1728 (void) blkCols;
1729
1730 return(sumIntoSystemMatrix(numPtRows, ptRows, numPtCols, ptCols, values));
1731 }
1732
1733 //***************************************************************************
1734 // put nonzero entries into the matrix data structure
1735 //---------------------------------------------------------------------------
1736
putIntoSystemMatrix(int numPtRows,const int * ptRows,int numPtCols,const int * ptCols,const double * const * values)1737 int HYPRE_LinSysCore::putIntoSystemMatrix(int numPtRows, const int* ptRows,
1738 int numPtCols, const int* ptCols,
1739 const double* const* values)
1740 {
1741 int i, j, localRow, newLeng, *tempInd, colIndex, index, localNRows;
1742 int sortFlag;
1743 double *tempVal;
1744
1745 //-------------------------------------------------------------------
1746 // diagnostic message for high output level only
1747 //-------------------------------------------------------------------
1748
1749 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
1750 printf("%4d : HYPRE_LSC::entering putIntoSystemMatrix.\n",mypid_);
1751
1752 //-------------------------------------------------------------------
1753 // error checking
1754 //-------------------------------------------------------------------
1755
1756 if ( systemAssembled_ == 1 )
1757 {
1758 printf("putIntoSystemMatrix ERROR : matrix already assembled\n");
1759 exit(1);
1760 }
1761 if ( numPtRows <= 0 || numPtCols <= 0 )
1762 {
1763 printf("%4d : putIntoSystemMatrix ERROR : invalid numPt.\n",mypid_);
1764 return (-1);
1765 }
1766
1767 //-------------------------------------------------------------------
1768 // in case the matrix is loaded from scratch
1769 //-------------------------------------------------------------------
1770
1771 if ( rowLengths_ == NULL && colIndices_ == NULL )
1772 {
1773 localNRows = localEndRow_ - localStartRow_ + 1;
1774 if ( localNRows > 0 )
1775 {
1776 rowLengths_ = new int[localNRows];
1777 colIndices_ = new int*[localNRows];
1778 colValues_ = new double*[localNRows];
1779 }
1780 for ( i = 0; i < localNRows; i++ )
1781 {
1782 rowLengths_[i] = 0;
1783 colIndices_[i] = NULL;
1784 colValues_[i] = NULL;
1785 }
1786 }
1787
1788 //-------------------------------------------------------------------
1789 // first adjust memory allocation (conservative)
1790 //-------------------------------------------------------------------
1791
1792 for ( i = 0; i < numPtRows; i++ )
1793 {
1794 localRow = ptRows[i] - localStartRow_ + 1;
1795 if ( rowLengths_[localRow] > 0 )
1796 {
1797 newLeng = rowLengths_[localRow] + numPtCols;
1798 tempInd = new int[newLeng];
1799 tempVal = new double[newLeng];
1800 for ( j = 0; j < rowLengths_[localRow]; j++ )
1801 {
1802 tempVal[j] = colValues_[localRow][j];
1803 tempInd[j] = colIndices_[localRow][j];
1804 }
1805 delete [] colValues_[localRow];
1806 delete [] colIndices_[localRow];
1807 colValues_[localRow] = tempVal;
1808 colIndices_[localRow] = tempInd;
1809 }
1810 else
1811 {
1812 if ( colIndices_[localRow] != NULL ) delete [] colIndices_[localRow];
1813 if ( colValues_[localRow] != NULL ) delete [] colValues_[localRow];
1814 colIndices_[localRow] = new int[numPtCols];
1815 colValues_[localRow] = new double[numPtCols];
1816 }
1817 }
1818
1819 //-------------------------------------------------------------------
1820 // load the local matrix
1821 //-------------------------------------------------------------------
1822
1823 for ( i = 0; i < numPtRows; i++ )
1824 {
1825 localRow = ptRows[i] - localStartRow_ + 1;
1826 if ( rowLengths_[localRow] > 0 )
1827 {
1828 newLeng = rowLengths_[localRow];
1829 tempInd = colIndices_[localRow];
1830 tempVal = colValues_[localRow];
1831 for ( j = 0; j < numPtCols; j++ )
1832 {
1833 colIndex = ptCols[j] + 1;
1834 index = hypre_BinarySearch(tempInd, colIndex, newLeng);
1835 if ( index < 0 )
1836 {
1837 tempInd[rowLengths_[localRow]] = colIndex;
1838 tempVal[rowLengths_[localRow]++] = values[i][j];
1839 }
1840 else tempVal[index] = values[i][j];
1841 }
1842 newLeng = rowLengths_[localRow];
1843 hypre_qsort1( tempInd, tempVal, 0, newLeng-1 );
1844 }
1845 else
1846 {
1847 tempInd = colIndices_[localRow];
1848 tempVal = colValues_[localRow];
1849 for ( j = 0; j < numPtCols; j++ )
1850 {
1851 colIndex = ptCols[j] + 1;
1852 tempInd[j] = colIndex;
1853 tempVal[j] = values[i][j];
1854 }
1855 sortFlag = 0;
1856 for ( j = 1; j < numPtCols; j++ )
1857 if ( tempInd[j] < tempInd[j-1] ) sortFlag = 1;
1858 rowLengths_[localRow] = numPtCols;
1859 if ( sortFlag == 1 ) hypre_qsort1( tempInd, tempVal, 0, numPtCols-1 );
1860 }
1861 }
1862
1863 //-------------------------------------------------------------------
1864 // diagnostic message
1865 //-------------------------------------------------------------------
1866
1867 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
1868 printf("%4d : HYPRE_LSC::leaving putIntoSystemMatrix.\n",mypid_);
1869
1870 return (0);
1871 }
1872
1873 //***************************************************************************
1874 // get the length of the local row
1875 //---------------------------------------------------------------------------
1876
getMatrixRowLength(int row,int & length)1877 int HYPRE_LinSysCore::getMatrixRowLength(int row, int& length)
1878 {
1879 int *colInd, rowLeng;
1880 double *colVal;
1881 HYPRE_ParCSRMatrix A_csr;
1882
1883 if ((row+1) < localStartRow_ || (row+1) > localEndRow_) return (-1);
1884 if ( systemAssembled_ == 0 )
1885 {
1886 if ( rowLengths_ == NULL ) return (-1);
1887 length = rowLengths_[row+1];
1888 }
1889 else
1890 {
1891 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
1892 HYPRE_ParCSRMatrixGetRow(A_csr,row,&rowLeng,&colInd,&colVal);
1893 length = rowLeng;
1894 HYPRE_ParCSRMatrixRestoreRow(A_csr,row,&rowLeng,&colInd,&colVal);
1895 }
1896 return (0);
1897 }
1898
1899 //***************************************************************************
1900 // get the data of a local row
1901 //---------------------------------------------------------------------------
1902
getMatrixRow(int row,double * coefs,int * indices,int len,int & rowLength)1903 int HYPRE_LinSysCore::getMatrixRow(int row, double* coefs, int* indices,
1904 int len, int& rowLength)
1905 {
1906 int i, rowIndex, rowLeng, *colInd, minLeng;
1907 double *colVal;
1908 HYPRE_ParCSRMatrix A_csr;
1909
1910 if ( systemAssembled_ == 0 )
1911 {
1912 rowIndex = row + 1;
1913 if (rowIndex < localStartRow_ || rowIndex > localEndRow_) return (-1);
1914 if ( rowLengths_ == NULL || colIndices_ == NULL ) return (-1);
1915 rowLeng = rowLengths_[rowIndex];
1916 colInd = colIndices_[rowIndex];
1917 colVal = colValues_[rowIndex];
1918 minLeng = len;
1919 if ( minLeng > rowLeng ) minLeng = rowLeng;
1920 for( i = 0; i < minLeng; i++ )
1921 {
1922 coefs[i] = colVal[i];
1923 indices[i] = colInd[i];
1924 }
1925 rowLength = rowLeng;
1926 }
1927 else
1928 {
1929 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
1930 rowIndex = row + 1;
1931 if (rowIndex < localStartRow_ || rowIndex > localEndRow_) return (-1);
1932 rowIndex--;
1933 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowLeng,&colInd,&colVal);
1934 minLeng = len;
1935 if ( minLeng > rowLeng ) minLeng = rowLeng;
1936 for( i = 0; i < minLeng; i++ )
1937 {
1938 coefs[i] = colVal[i];
1939 indices[i] = colInd[i];
1940 }
1941 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowLeng,&colInd,&colVal);
1942 rowLength = rowLeng;
1943 }
1944 return (0);
1945 }
1946
1947 //***************************************************************************
1948 // input is 1-based, but HYPRE vectors are 0-based
1949 //---------------------------------------------------------------------------
1950
sumIntoRHSVector(int num,const double * values,const int * indices)1951 int HYPRE_LinSysCore::sumIntoRHSVector(int num, const double* values,
1952 const int* indices)
1953 {
1954 int i, *localInds;
1955
1956 //-------------------------------------------------------------------
1957 // diagnostic message
1958 //-------------------------------------------------------------------
1959
1960 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
1961 {
1962 printf("%4d : HYPRE_LSC::entering sumIntoRHSVector.\n", mypid_);
1963 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 6 )
1964 {
1965 for ( i = 0; i < num; i++ )
1966 printf("%d : sumIntoRHSVector - %d = %e.\n", mypid_, indices[i],
1967 values[i]);
1968 }
1969 }
1970
1971 //-------------------------------------------------------------------
1972 // change the incoming indices to 0-based before loading
1973 //-------------------------------------------------------------------
1974
1975 localInds = new int[num];
1976 for ( i = 0; i < num; i++ ) // change to 0-based
1977 {
1978 #if defined(NOFEI)
1979 if ( indices[i] >= localStartRow_ && indices[i] <= localEndRow_ )
1980 localInds[i] = indices[i] - 1;
1981 #else
1982 if ((indices[i]+1) >= localStartRow_ && (indices[i]+1) <= localEndRow_)
1983 localInds[i] = indices[i];
1984 #endif
1985 else
1986 {
1987 printf("%d : sumIntoRHSVector ERROR - index %d out of range.\n",
1988 mypid_, indices[i]);
1989 exit(1);
1990 }
1991 }
1992
1993 HYPRE_IJVectorAddToValues(HYb_, num, (const int *) localInds,
1994 (const double *) values);
1995 delete [] localInds;
1996
1997 //-------------------------------------------------------------------
1998 // diagnostic message
1999 //-------------------------------------------------------------------
2000
2001 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
2002 printf("%4d : HYPRE_LSC::leaving sumIntoRHSVector.\n", mypid_);
2003 return (0);
2004 }
2005
2006 //***************************************************************************
2007 // This function scatters (puts) values into the linear-system's
2008 // currently selected RHS vector.
2009 // num is how many values are being passed,
2010 // indices holds the global 0-based 'row-numbers' into which the values go,
2011 // and values holds the actual coefficients to be scattered.
2012 //---------------------------------------------------------------------------
2013
putIntoRHSVector(int num,const double * values,const int * indices)2014 int HYPRE_LinSysCore::putIntoRHSVector(int num, const double* values,
2015 const int* indices)
2016 {
2017 int i, index;
2018
2019 if ((numRHSs_ == 0) && (HYb_ == NULL)) return (0);
2020
2021 for ( i = 0; i < num; i++ )
2022 {
2023 index = indices[i];
2024 if (index < localStartRow_-1 || index >= localEndRow_) continue;
2025 HYPRE_IJVectorSetValues(HYb_, 1, (const int *) &index,
2026 (const double *) &(values[i]));
2027 }
2028 return (0);
2029 }
2030
2031 //***************************************************************************
2032 // This function gets values from the linear-system's
2033 // currently selected RHS vector.
2034 // num is how many values are being requested,
2035 // indices holds the requested global 0-based 'row-numbers',
2036 // and values will hold the returned coefficients.
2037 //---------------------------------------------------------------------------
2038
getFromRHSVector(int num,double * values,const int * indices)2039 int HYPRE_LinSysCore::getFromRHSVector(int num, double* values,
2040 const int* indices)
2041 {
2042 int i, index;
2043
2044 if ((numRHSs_ == 0) && (HYb_ == NULL)) return (0);
2045
2046 for ( i = 0; i < num; i++ )
2047 {
2048 index = indices[i];
2049 if (index < localStartRow_-1 || index >= localEndRow_) continue;
2050 HYPRE_IJVectorGetValues(HYb_,1,&index,&(values[i]));
2051 }
2052 return (0);
2053 }
2054
2055 //***************************************************************************
2056 // start assembling the matrix into its internal format
2057 //---------------------------------------------------------------------------
2058
matrixLoadComplete()2059 int HYPRE_LinSysCore::matrixLoadComplete()
2060 {
2061 int i, j, numLocalEqns, leng, eqnNum, nnz, *newColInd=NULL;
2062 int maxRowLeng, newLeng, rowSize, *colInd, nrows;
2063 double *newColVal=NULL, *colVal, value;
2064 char fname[40];
2065 FILE *fp;
2066 HYPRE_ParCSRMatrix A_csr;
2067 HYPRE_ParVector b_csr;
2068 HYPRE_SlideReduction *slideObj;
2069
2070 //-------------------------------------------------------------------
2071 // diagnostic message
2072 //-------------------------------------------------------------------
2073
2074 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
2075 printf("%4d : HYPRE_LSC::entering matrixLoadComplete.\n",mypid_);
2076
2077 //-------------------------------------------------------------------
2078 // Write MLI FEData information to a file
2079 //-------------------------------------------------------------------
2080
2081 #ifdef HAVE_MLI
2082 if ( haveFEData_ && feData_ != NULL )
2083 {
2084 char filename[100];
2085 if ( (HYOutputLevel_ & HYFEI_PRINTFEINFO) )
2086 {
2087 strcpy( filename, "fedata" );
2088 HYPRE_LSI_MLIFEDataWriteToFile( feData_, filename );
2089 }
2090 }
2091 #endif
2092
2093 if ( matrixPartition_ == 2 ) matrixPartition_ = 1;
2094
2095 //-------------------------------------------------------------------
2096 // if the matrix has not been assembled or it has been reset
2097 //-------------------------------------------------------------------
2098
2099 if ( systemAssembled_ != 1 )
2100 {
2101 //----------------------------------------------------------------
2102 // set system matrix initialization parameters
2103 //----------------------------------------------------------------
2104
2105 HYPRE_IJMatrixSetRowSizes(HYA_, rowLengths_);
2106 HYPRE_IJMatrixInitialize(HYA_);
2107
2108 //----------------------------------------------------------------
2109 // load the matrix stored locally to a HYPRE matrix
2110 //----------------------------------------------------------------
2111
2112 numLocalEqns = localEndRow_ - localStartRow_ + 1;
2113 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
2114 {
2115 printf("%4d : HYPRE_LSC::matrixLoadComplete - NEqns = %d.\n",
2116 mypid_, numLocalEqns);
2117 }
2118 maxRowLeng = 0;
2119 for ( i = 0; i < numLocalEqns; i++ )
2120 {
2121 leng = rowLengths_[i];
2122 if ( leng > maxRowLeng ) maxRowLeng = leng;
2123 }
2124 if ( maxRowLeng > 0 )
2125 {
2126 newColInd = new int[maxRowLeng];
2127 newColVal = new double[maxRowLeng];
2128 }
2129 nnz = 0;
2130 for ( i = 0; i < numLocalEqns; i++ )
2131 {
2132 eqnNum = localStartRow_ - 1 + i;
2133 leng = rowLengths_[i];
2134 newLeng = 0;
2135 for ( j = 0; j < leng; j++ )
2136 {
2137 if ( habs(colValues_[i][j]) >= truncThresh_ )
2138 {
2139 newColInd[newLeng] = colIndices_[i][j] - 1;
2140 newColVal[newLeng++] = colValues_[i][j];
2141 }
2142 }
2143 HYPRE_IJMatrixSetValues(HYA_, 1, &newLeng,(const int *) &eqnNum,
2144 (const int *) newColInd, (const double *) newColVal);
2145 delete [] colValues_[i];
2146 if ( memOptimizerFlag_ != 0 ) delete [] colIndices_[i];
2147 nnz += newLeng;
2148 }
2149 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
2150 {
2151 printf("%4d : HYPRE_LSC::matrixLoadComplete - nnz = %d.\n",
2152 mypid_, nnz);
2153 }
2154 delete [] colValues_;
2155 colValues_ = NULL;
2156 if ( memOptimizerFlag_ != 0 )
2157 {
2158 delete [] colIndices_;
2159 colIndices_ = NULL;
2160 }
2161 if ( maxRowLeng > 0 )
2162 {
2163 delete [] newColInd;
2164 delete [] newColVal;
2165 }
2166 HYPRE_IJMatrixAssemble(HYA_);
2167 systemAssembled_ = 1;
2168 projectCurrSize_ = 0;
2169 currA_ = HYA_;
2170 currB_ = HYb_;
2171 currX_ = HYx_;
2172 currR_ = HYr_;
2173 if (slideObj_ != NULL)
2174 {
2175 slideObj = (HYPRE_SlideReduction *) slideObj_;
2176 delete slideObj;
2177 }
2178 slideObj_ = NULL;
2179 }
2180
2181 //-------------------------------------------------------------------
2182 // diagnostics : print the matrix and rhs to files
2183 //-------------------------------------------------------------------
2184
2185 if ( (HYOutputLevel_ & HYFEI_PRINTMAT) &&
2186 (!(HYOutputLevel_ & HYFEI_PRINTREDMAT) ) )
2187 {
2188 if ( HYOutputLevel_ & HYFEI_PRINTPARCSRMAT )
2189 {
2190 printf("%4d : HYPRE_LSC::print the matrix/rhs to files(1)\n",mypid_);
2191 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
2192 sprintf(fname, "HYPRE_Mat");
2193 HYPRE_ParCSRMatrixPrint( A_csr, fname );
2194 HYPRE_IJVectorGetObject(HYb_, (void **) &b_csr);
2195 sprintf(fname, "HYPRE_RHS");
2196 HYPRE_ParVectorPrint( b_csr, fname );
2197 }
2198 else
2199 {
2200 printf("%4d : HYPRE_LSC::print the matrix/rhs to files(2)\n",mypid_);
2201 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
2202 sprintf(fname, "hypre_mat.out.%d",mypid_);
2203 fp = fopen(fname,"w");
2204 nrows = localEndRow_ - localStartRow_ + 1;
2205 nnz = 0;
2206 for ( i = localStartRow_-1; i <= localEndRow_-1; i++ )
2207 {
2208 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
2209 for ( j = 0; j < rowSize; j++ ) if ( colVal[j] != 0.0 ) nnz++;
2210 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
2211 }
2212 fprintf(fp, "%6d %7d \n", nrows, nnz);
2213 for ( i = localStartRow_-1; i <= localEndRow_-1; i++ )
2214 {
2215 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
2216 for (j = 0; j < rowSize; j++)
2217 {
2218 if ( colVal[j] != 0.0 )
2219 fprintf(fp, "%6d %6d %25.16e \n",i+1,colInd[j]+1,colVal[j]);
2220 }
2221 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
2222 }
2223 fclose(fp);
2224 sprintf(fname, "hypre_rhs.out.%d",mypid_);
2225 fp = fopen(fname,"w");
2226 fprintf(fp, "%6d \n", nrows);
2227 for ( i = localStartRow_-1; i <= localEndRow_-1; i++ )
2228 {
2229 HYPRE_IJVectorGetValues(HYb_, 1, &i, &value);
2230 fprintf(fp, "%6d %25.16e \n", i+1, value);
2231 }
2232 fclose(fp);
2233 MPI_Barrier(comm_);
2234 }
2235 if ( HYOutputLevel_ & HYFEI_STOPAFTERPRINT ) exit(1);
2236 }
2237 if (FEI_mixedDiagFlag_)
2238 {
2239 for ( i = 0; i < localEndRow_-localStartRow_+1; i++ )
2240 {
2241 FEI_mixedDiag_[i] *= 0.125;
2242 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
2243 printf("Mixed diag %5d = %e\n", i, FEI_mixedDiag_[i]);
2244 }
2245 }
2246
2247 //-------------------------------------------------------------------
2248 // diagnostic message
2249 //-------------------------------------------------------------------
2250
2251 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
2252 printf("%4d : HYPRE_LSC::leaving matrixLoadComplete.\n",mypid_);
2253 return (0);
2254 }
2255
2256 //***************************************************************************
2257 // put nodal information
2258 //---------------------------------------------------------------------------
2259
putNodalFieldData(int fieldID,int fieldSize,int * nodeNumbers,int numNodes,const double * data)2260 int HYPRE_LinSysCore::putNodalFieldData(int fieldID, int fieldSize,
2261 int* nodeNumbers, int numNodes, const double* data)
2262 {
2263 int i, j, **nodeFieldIDs, nodeFieldID, *procNRows, nRows, errCnt;
2264 int blockID, *blockIDs, *eqnNumbers, *iArray, newNumEdges;
2265 //int checkFieldSize;
2266 int *aleNodeNumbers, index, newNumNodes, numEdges;
2267 double *newData;
2268
2269 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
2270 {
2271 printf("%4d : HYPRE_LSC::entering putNodalFieldData.\n",mypid_);
2272 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 && mypid_ == 0 )
2273 {
2274 printf(" putNodalFieldData : fieldSize = %d\n", fieldSize);
2275 printf(" putNodalFieldData : fieldID = %d\n", fieldID);
2276 printf(" putNodalFieldData : numNodes = %d\n", numNodes);
2277 }
2278 }
2279
2280 //-------------------------------------------------------------------
2281 // This part is for loading the nodal coordinate information.
2282 // The node IDs in nodeNumbers are the one used in FEI (and thus
2283 // corresponds to the ones in the system matrix using lookup)
2284 //-------------------------------------------------------------------
2285
2286 if ( fieldID == -3 || fieldID == -25333 )
2287 {
2288 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
2289 {
2290 for ( i = 0; i < numNodes; i++ )
2291 for ( j = 0; j < fieldSize; j++ )
2292 printf("putNodalFieldData : %4d %2d = %e\n",i,j,
2293 data[i*fieldSize+j]);
2294 }
2295 if ( HYPreconID_ == HYMLI && lookup_ != NULL )
2296 {
2297 blockIDs = (int *) lookup_->getElemBlockIDs();
2298 blockID = blockIDs[0];
2299 nodeFieldIDs = (int **) lookup_->getFieldIDsTable(blockID);
2300 nodeFieldID = nodeFieldIDs[0][0];
2301 //checkFieldSize = lookup_->getFieldSize(nodeFieldID);
2302 //hypre_assert( checkFieldSize == fieldSize );
2303 eqnNumbers = new int[numNodes];
2304 newData = new double[numNodes*fieldSize];
2305 newNumNodes = 0;
2306 for ( i = 0; i < numNodes*fieldSize; i++ ) newData[i] = -99999.9;
2307 for ( i = 0; i < numNodes; i++ )
2308 {
2309 index = lookup_->getEqnNumber(nodeNumbers[i],nodeFieldID);
2310
2311 /* ======
2312 This should ultimately be taken out even for newer ale3d implementation
2313 =====*/
2314 if ( index >= localStartRow_-1 && index < localEndRow_)
2315 {
2316 if ( newData[newNumNodes*fieldSize] == -99999.9 )
2317 {
2318 for ( j = 0; j < fieldSize; j++ )
2319 newData[newNumNodes*fieldSize+j] = data[i*fieldSize+j];
2320 eqnNumbers[newNumNodes++] = index;
2321 }
2322 }
2323 }
2324 nRows = localEndRow_ - localStartRow_ + 1;
2325 if ( MLI_NodalCoord_ == NULL )
2326 {
2327 MLI_EqnNumbers_ = new int[nRows/fieldSize];
2328 for (i=0; i<nRows/fieldSize; i++)
2329 MLI_EqnNumbers_[i] = localStartRow_ - 1 + i * fieldSize;
2330 MLI_NodalCoord_ = new double[localEndRow_-localStartRow_+1];
2331 for (i=0; i<nRows; i++) MLI_NodalCoord_[i] = -99999.0;
2332 MLI_FieldSize_ = fieldSize;
2333 MLI_NumNodes_ = nRows / fieldSize;
2334 }
2335 for ( i = 0; i < newNumNodes; i++ )
2336 {
2337 index = eqnNumbers[i] - localStartRow_ + 1;
2338 for ( j = 0; j < fieldSize; j++ )
2339 MLI_NodalCoord_[index+j] = newData[i*fieldSize+j];
2340 }
2341 delete [] eqnNumbers;
2342 delete [] newData;
2343 errCnt = 0;
2344 for (i = 0; i < nRows; i++)
2345 if (MLI_NodalCoord_[i] == -99999.0) errCnt++;
2346 if (errCnt > 0)
2347 printf("putNodalFieldData ERROR:incomplete nodal coordinates (%d %d).\n",
2348 errCnt, nRows);
2349 }
2350 else
2351 {
2352 if (nodeNumbers != NULL && numNodes != 0)
2353 {
2354 printf("putNodalFieldData WARNING : \n");
2355 printf(" set nodeNumbers = NULL, set numNodes = 0.\n");
2356 }
2357 nRows = localEndRow_ - localStartRow_ + 1;
2358 MLI_NodalCoord_ = new double[nRows];
2359 for (i=0; i<nRows; i++) MLI_NodalCoord_[i] = data[i];
2360 }
2361 }
2362
2363 //-------------------------------------------------------------------
2364 // This part is for loading the edge to (hypre-compatible) node list
2365 // for AMS (the node list is ordered with the edge equation numbers
2366 // and the node numbers are true node equation numbers passed in by
2367 // the application which obtains the true node eqn numbers via the
2368 // nodal FEI) ===> EdgeNodeList
2369 //-------------------------------------------------------------------
2370
2371 if (fieldID == -4)
2372 {
2373 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5)
2374 {
2375 for (i = 0; i < numNodes; i++)
2376 for (j = 0; j < fieldSize; j++)
2377 printf("putNodalFieldData : %4d %2d = %e\n",i,j,
2378 data[i*fieldSize+j]);
2379 }
2380 if (lookup_ != NULL && fieldSize == 2 &&
2381 numNodes > 0)
2382 {
2383 blockIDs = (int *) lookup_->getElemBlockIDs();
2384 blockID = blockIDs[0];
2385 nodeFieldIDs = (int **) lookup_->getFieldIDsTable(blockID);
2386 nodeFieldID = nodeFieldIDs[0][0];
2387 numEdges = numNodes;
2388 eqnNumbers = new int[numEdges];
2389 iArray = new int[numEdges*fieldSize];
2390 newNumEdges = 0;
2391 for (i = 0; i < numEdges; i++)
2392 {
2393 index = lookup_->getEqnNumber(nodeNumbers[i],nodeFieldID);
2394 if (index >= localStartRow_-1 && index < localEndRow_)
2395 {
2396 for (j = 0; j < fieldSize; j++)
2397 iArray[newNumEdges*fieldSize+j] = (int) data[i*fieldSize+j];
2398 eqnNumbers[newNumEdges++] = index;
2399 }
2400 }
2401 nRows = localEndRow_ - localStartRow_ + 1;
2402 if (AMSData_.EdgeNodeList_ != NULL) hypre_TFree(AMSData_.EdgeNodeList_, HYPRE_MEMORY_HOST);
2403 AMSData_.EdgeNodeList_ = NULL;
2404 if (newNumEdges > 0)
2405 {
2406 AMSData_.numEdges_ = nRows;
2407 AMSData_.EdgeNodeList_ = hypre_CTAlloc(HYPRE_BigInt, nRows*fieldSize, HYPRE_MEMORY_HOST);
2408 for (i = 0; i < nRows*fieldSize; i++)
2409 AMSData_.EdgeNodeList_[i] = -99999;
2410 for (i = 0; i < newNumEdges; i++)
2411 {
2412 index = eqnNumbers[i] - localStartRow_ + 1;
2413 for (j = 0; j < fieldSize; j++ )
2414 AMSData_.EdgeNodeList_[index*fieldSize+j] = iArray[i*fieldSize+j];
2415 }
2416 errCnt = 0;
2417 for (i = 0; i < nRows*fieldSize; i++)
2418 if (AMSData_.EdgeNodeList_[i] == -99999) errCnt++;
2419 if (errCnt > 0)
2420 printf("putNodalFieldData ERROR:incomplete AMS edge vertex list\n");
2421 }
2422 delete [] eqnNumbers;
2423 delete [] iArray;
2424 }
2425 }
2426
2427 //-------------------------------------------------------------------
2428 // This part is for converting node numbers to equations as well as
2429 // for loading the nodal coordinate information
2430 // (stored in NodeNumbers, NodalCoord, numNodes, numLocalNodes)
2431 //-------------------------------------------------------------------
2432
2433 if (fieldID == -5)
2434 {
2435 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5)
2436 {
2437 for (i = 0; i < numNodes; i++)
2438 for (j = 0; j < fieldSize; j++)
2439 printf("putNodalFieldData : %4d %2d = %e\n",i,j,
2440 data[i*fieldSize+j]);
2441 }
2442 if (lookup_ != NULL && fieldSize == 3)
2443 {
2444 blockIDs = (int *) lookup_->getElemBlockIDs();
2445 blockID = blockIDs[0];
2446 nodeFieldIDs = (int **) lookup_->getFieldIDsTable(blockID);
2447 nodeFieldID = nodeFieldIDs[0][0];
2448 if (AMSData_.NodeNumbers_ != NULL) hypre_TFree(AMSData_.NodeNumbers_, HYPRE_MEMORY_HOST);
2449 if (AMSData_.NodalCoord_ != NULL) hypre_TFree(AMSData_.NodalCoord_, HYPRE_MEMORY_HOST);
2450 AMSData_.NodeNumbers_ = NULL;
2451 AMSData_.NodalCoord_ = NULL;
2452 AMSData_.numNodes_ = 0;
2453 if (numNodes > 0)
2454 {
2455 AMSData_.numNodes_ = numNodes;
2456 AMSData_.numLocalNodes_ = localEndRow_ - localStartRow_ + 1;
2457 AMSData_.NodeNumbers_ = hypre_CTAlloc(HYPRE_BigInt, numNodes, HYPRE_MEMORY_HOST);
2458 AMSData_.NodalCoord_ = hypre_CTAlloc(HYPRE_Real, fieldSize*numNodes, HYPRE_MEMORY_HOST);
2459 for (i = 0; i < numNodes; i++)
2460 {
2461 index = lookup_->getEqnNumber(nodeNumbers[i],nodeFieldID);
2462 AMSData_.NodeNumbers_[i] = index;
2463 for (j = 0; j < fieldSize; j++)
2464 AMSData_.NodalCoord_[i*fieldSize+j] = data[i*fieldSize+j];
2465 }
2466 }
2467 }
2468 }
2469
2470 //-------------------------------------------------------------------
2471 // this is needed to set up the correct node equation map
2472 // (the FEI remaps the node IDs in the incoming nodeNumbers array.
2473 // to revert to the original ALE3D node numbers, it is passed in
2474 // here as data)
2475 //-------------------------------------------------------------------
2476
2477 else if ( fieldID == -49773 )
2478 {
2479 if ( HYPreconID_ == HYMLI && lookup_ != NULL )
2480 {
2481 blockIDs = (int *) lookup_->getElemBlockIDs();
2482 blockID = blockIDs[0];
2483 nodeFieldIDs = (int **) lookup_->getFieldIDsTable(blockID);
2484 nodeFieldID = nodeFieldIDs[0][0];
2485 //checkFieldSize = lookup_->getFieldSize(nodeFieldID);
2486 hypre_assert( fieldSize == 1 );
2487 aleNodeNumbers = new int[numNodes];
2488 eqnNumbers = new int[numNodes];
2489 for ( i = 0; i < numNodes; i++ )
2490 {
2491 aleNodeNumbers[i] = (int) data[i];
2492 eqnNumbers[i] = lookup_->getEqnNumber(nodeNumbers[i],nodeFieldID);
2493 }
2494 procNRows = new int[numProcs_];
2495 for ( i = 0; i < numProcs_; i++ ) procNRows[i] = 0;
2496 procNRows[mypid_] = localEndRow_;
2497 iArray = procNRows;
2498 procNRows = new int[numProcs_+1];
2499 for ( i = 0; i <= numProcs_; i++ ) procNRows[i] = 0;
2500 MPI_Allreduce(iArray,&(procNRows[1]),numProcs_,MPI_INT,MPI_SUM,comm_);
2501 delete [] iArray;
2502 HYPRE_LSI_MLICreateNodeEqnMap(HYPrecon_, numNodes, aleNodeNumbers,
2503 eqnNumbers, procNRows);
2504 delete [] procNRows;
2505 delete [] eqnNumbers;
2506 delete [] aleNodeNumbers;
2507 }
2508 }
2509
2510 //-------------------------------------------------------------------
2511 // This part is for initializing the nodal coordinates vectors
2512 // for AMS in the New Fashion which allows multiply-connected edges
2513 //-------------------------------------------------------------------
2514
2515 if (fieldID == -100) {
2516
2517 errCnt = 0;
2518 //USAGE:
2519 //FieldID=-100: Initialize space for vectors X,Y,Z for AMS
2520 //double *data_dummy;
2521 //memAllocDoubles( 1, &data_dummy, NULL);
2522 //lsc->putNodalFieldData(-100,1,globalOffsets,numProcs+1,data_dummy);
2523 //MPI_Comm_rank(comm, &mypid_);
2524 //MPI_Comm_size(comm, &numProcs_);
2525 //localStartRowAMSV_ = nodeNumbers[mypid_] + 1;
2526 //localEndRowAMSV_ = nodeNumbers[mypid_+1];
2527 localStartRowAMSV_ = fieldSize;
2528 localEndRowAMSV_ = numNodes;
2529
2530 MPI_Comm_rank(comm_, &mypid_);
2531 MPI_Comm_size(comm_, &numProcs_);
2532
2533 errCnt += HYPRE_IJVectorCreate(comm_, localStartRowAMSV_, localEndRowAMSV_, &amsX_);
2534 errCnt += HYPRE_IJVectorCreate(comm_, localStartRowAMSV_, localEndRowAMSV_, &amsY_);
2535 errCnt += HYPRE_IJVectorCreate(comm_, localStartRowAMSV_, localEndRowAMSV_, &amsZ_);
2536 errCnt += HYPRE_IJVectorSetObjectType(amsX_, HYPRE_PARCSR);
2537 errCnt += HYPRE_IJVectorSetObjectType(amsY_, HYPRE_PARCSR);
2538 errCnt += HYPRE_IJVectorSetObjectType(amsZ_, HYPRE_PARCSR);
2539 errCnt += HYPRE_IJVectorInitialize(amsX_);
2540 errCnt += HYPRE_IJVectorInitialize(amsY_);
2541 errCnt += HYPRE_IJVectorInitialize(amsZ_);
2542 }
2543 if (fieldID == -101 || fieldID == -102 || fieldID == -103) {
2544
2545 errCnt = 0;
2546 //-------------------------------------------------------------------
2547 // This part is for loading the nodal coordinates vectors
2548 // for AMS in the New Fashion which allows multiply-connected edges
2549 //-------------------------------------------------------------------
2550
2551 //Usage
2552 //FieldID=-101,2,3 Put values for vectors X,Y,Z for AMS
2553 //lsc->putNodalFieldData(-101,1,nodeNumbers_dummy,numNodes,datax_dummy);
2554 //lsc->putNodalFieldData(-102,1,nodeNumbers_dummy,numNodes,datay_dummy);
2555 //lsc->putNodalFieldData(-103,1,nodeNumbers_dummy,numNodes,dataz_dummy);
2556
2557 switch( fieldID ) {
2558 case( -101 ) :
2559 //errCnt += HYPRE_IJVectorSetValues(amsX_, numNodes, localInds, data);
2560 errCnt += HYPRE_IJVectorSetValues(amsX_, numNodes, nodeNumbers, data);
2561 break;
2562 case( -102 ) :
2563 errCnt += HYPRE_IJVectorSetValues(amsY_, numNodes, nodeNumbers, data);
2564 break;
2565 case( -103 ) :
2566 errCnt += HYPRE_IJVectorSetValues(amsZ_, numNodes, nodeNumbers, data);
2567 break;
2568 default :
2569 printf("%d : PutNodalFieldData, FieldID=-101,-102,-103, ERROR - FieldID %d out of range.\n",
2570 mypid_,fieldID);
2571 exit(1);
2572 break;
2573 }
2574 }
2575 if (fieldID == -200) {
2576
2577 errCnt = 0;
2578 //USAGE
2579 //FieldID=-200: Initialize the Matrix G for AMS;
2580 //double *globalColOffsetsd;
2581 //memAllocDoubles( numProcs+1, &globalColOffsetsd, NULL);
2582 //for (int iproc=0;iproc<numProcs+1;iproc++){
2583 // globalColOffsetsd[iproc] = double( globalColOffsets[iproc] );
2584 //}
2585 //lsc->putNodalFieldData(-200,numProcs+1,globalRowOffsets,numProcs+1,globalColOffsetsd);
2586 //memFree(&globalRowOffsets,&globalColOffsets,&globalColOffsetsd,NULL);
2587
2588 if( amsG_ != NULL) errCnt += HYPRE_IJMatrixDestroy(amsG_);
2589 localStartRowAMSG_ = nodeNumbers[0];
2590 localEndRowAMSG_ = nodeNumbers[1];
2591 localStartColAMSG_ = int( data[0] );
2592 localEndColAMSG_ = int( data[1] );
2593
2594
2595 errCnt += HYPRE_IJMatrixCreate(comm_, localStartRowAMSG_, localEndRowAMSG_,
2596 localStartColAMSG_, localEndColAMSG_, &amsG_);
2597
2598
2599 errCnt += HYPRE_IJMatrixSetObjectType(amsG_, HYPRE_PARCSR);
2600
2601 }
2602
2603 if (fieldID == -201) {
2604 //USAGE
2605 //FieldID=-201: Set Row Sizes
2606 errCnt = 0;
2607 errCnt += HYPRE_IJMatrixSetRowSizes(amsG_, nodeNumbers);
2608 errCnt += HYPRE_IJMatrixInitialize(amsG_);
2609 }
2610
2611 if (fieldID == -202) {
2612 //USAGE
2613 //FieldID=-202: Put Values into the GMatrix
2614 //lsc->putNodalFieldData(-201,eindices[iedge],iindicesr,isize,icoefsr);
2615
2616 errCnt += HYPRE_IJMatrixSetValues(amsG_, 1, &numNodes, &fieldSize, nodeNumbers, data);
2617
2618 }
2619
2620 if (fieldID == -203) {
2621 //USAGE
2622 //FieldID=-203: Final Assembly of G Matrix
2623 //int *nodenumbers_dummy2;
2624 //double *data_dummy2;
2625 //memAllocInts( 1, &nodenumbers_dummy2, NULL);
2626 //memAllocDoubles( 1, &data_dummy2, NULL);
2627 //lsc->putNodalFieldData(-202,1,nodenumbers_dummy2,1,data_dummy2);
2628 //memFree( &nodenumbers_dummy2, &data_dummy2, NULL);
2629 errCnt = HYPRE_IJMatrixAssemble(amsG_);
2630 //this is also a good time to do this:
2631 errCnt += HYPRE_IJVectorAssemble(amsX_);
2632 errCnt += HYPRE_IJVectorAssemble(amsY_);
2633 errCnt += HYPRE_IJVectorAssemble(amsZ_);
2634
2635 }
2636
2637 // fieldID's -300, -301, -302, -303 are for calculating the D0 Matrix
2638 // This is purely for debugging.
2639 // The D0 Matrix maps nodes to edges
2640 if (fieldID == -300) {
2641
2642 errCnt = 0;
2643 //USAGE
2644 //FieldID=-300: Initialize the Matrix D0 for AMS;
2645 //double *globalColOffsetsd;
2646 //memAllocDoubles( numProcs+1, &globalColOffsetsd, NULL);
2647 //for (int iproc=0;iproc<numProcs+1;iproc++){
2648 // globalColOffsetsd[iproc] = double( globalColOffsets[iproc] );
2649 //}
2650 //lsc->putNodalFieldData(-200,numProcs+1,globalRowOffsets,numProcs+1,globalColOffsetsd);
2651 //memFree(&globalRowOffsets,&globalColOffsets,&globalColOffsetsd,NULL);
2652
2653 if( amsD0_ != NULL) errCnt += HYPRE_IJMatrixDestroy(amsD0_);
2654 int localStartRowAMSD_ = nodeNumbers[0];
2655 int localEndRowAMSD_ = nodeNumbers[1];
2656 int localStartColAMSD_ = int( data[0] );
2657 int localEndColAMSD_ = int( data[1] );
2658
2659
2660 errCnt += HYPRE_IJMatrixCreate(comm_, localStartRowAMSD_, localEndRowAMSD_,
2661 localStartColAMSD_, localEndColAMSD_, &amsD0_);
2662
2663
2664 errCnt += HYPRE_IJMatrixSetObjectType(amsD0_, HYPRE_PARCSR);
2665
2666 }
2667
2668 if (fieldID == -301) {
2669 //USAGE
2670 //FieldID=-301: Set Row Sizes for the D0 Matrix
2671 errCnt = 0;
2672 errCnt += HYPRE_IJMatrixSetRowSizes(amsD0_, nodeNumbers);
2673 errCnt += HYPRE_IJMatrixInitialize(amsD0_);
2674 }
2675
2676 if (fieldID == -302) {
2677 //USAGE
2678 //FieldID=-302: Put Values into the D0 Matrix
2679 //lsc->putNodalFieldData(-201,eindices[iedge],iindicesr,isize,icoefsr);
2680
2681 errCnt += HYPRE_IJMatrixSetValues(amsD0_, 1, &numNodes, &fieldSize, nodeNumbers, data);
2682
2683 }
2684
2685 if (fieldID == -303) {
2686 //USAGE
2687 //FieldID=-303: Final Assembly of D0 Matrix
2688 //int *nodenumbers_dummy2;
2689 //double *data_dummy2;
2690 //memAllocInts( 1, &nodenumbers_dummy2, NULL);
2691 //memAllocDoubles( 1, &data_dummy2, NULL);
2692 //lsc->putNodalFieldData(-202,1,nodenumbers_dummy2,1,data_dummy2);
2693 //memFree( &nodenumbers_dummy2, &data_dummy2, NULL);
2694 errCnt = HYPRE_IJMatrixAssemble(amsD0_);
2695
2696 bool printD0matrix = true;
2697 if(printD0matrix) {
2698 HYPRE_ParCSRMatrix D0_csr;
2699 HYPRE_IJMatrixGetObject(amsD0_, (void **) &D0_csr);
2700 HYPRE_ParCSRMatrixPrint(D0_csr, "D0.parmatrix");
2701 }
2702 }
2703
2704 // fieldID's -400, -401, -402, -403 are for calculating the D1 Matrix
2705 // This is purely for debugging.
2706 // The D1 Matrix maps edges to nodes
2707 if (fieldID == -400) {
2708
2709 errCnt = 0;
2710 //USAGE
2711 //FieldID=-400: Initialize the Matrix D1 for AMS;
2712 //double *globalColOffsetsd;
2713 //memAllocDoubles( numProcs+1, &globalColOffsetsd, NULL);
2714 //for (int iproc=0;iproc<numProcs+1;iproc++){
2715 // globalColOffsetsd[iproc] = double( globalColOffsets[iproc] );
2716 //}
2717 //lsc->putNodalFieldData(-200,numProcs+1,globalRowOffsets,numProcs+1,globalColOffsetsd);
2718 //memFree(&globalRowOffsets,&globalColOffsets,&globalColOffsetsd,NULL);
2719
2720 if( amsD1_ != NULL) errCnt += HYPRE_IJMatrixDestroy(amsD1_);
2721 int localStartRowAMSD_ = nodeNumbers[0];
2722 int localEndRowAMSD_ = nodeNumbers[1];
2723 int localStartColAMSD_ = int( data[0] );
2724 int localEndColAMSD_ = int( data[1] );
2725
2726
2727 errCnt += HYPRE_IJMatrixCreate(comm_, localStartRowAMSD_, localEndRowAMSD_,
2728 localStartColAMSD_, localEndColAMSD_, &amsD1_);
2729
2730
2731 errCnt += HYPRE_IJMatrixSetObjectType(amsD1_, HYPRE_PARCSR);
2732
2733 }
2734
2735 if (fieldID == -401) {
2736 //USAGE
2737 //FieldID=-401: Set Row Sizes for the D1 Matrix
2738 errCnt = 0;
2739 errCnt += HYPRE_IJMatrixSetRowSizes(amsD1_, nodeNumbers);
2740 errCnt += HYPRE_IJMatrixInitialize(amsD1_);
2741 }
2742
2743 if (fieldID == -402) {
2744 //USAGE
2745 //FieldID=-402: Put Values into the D1 Matrix
2746 //lsc->putNodalFieldData(-201,eindices[iedge],iindicesr,isize,icoefsr);
2747
2748 errCnt += HYPRE_IJMatrixSetValues(amsD1_, 1, &numNodes, &fieldSize, nodeNumbers, data);
2749
2750 }
2751
2752 if (fieldID == -403) {
2753 //USAGE
2754 //FieldID=-403: Final Assembly of D1 Matrix
2755 //int *nodenumbers_dummy2;
2756 //double *data_dummy2;
2757 //memAllocInts( 1, &nodenumbers_dummy2, NULL);
2758 //memAllocDoubles( 1, &data_dummy2, NULL);
2759 //lsc->putNodalFieldData(-202,1,nodenumbers_dummy2,1,data_dummy2);
2760 //memFree( &nodenumbers_dummy2, &data_dummy2, NULL);
2761 errCnt = HYPRE_IJMatrixAssemble(amsD1_);
2762
2763 bool printD1matrix = true;
2764 if(printD1matrix) {
2765 HYPRE_ParCSRMatrix D1_csr;
2766 HYPRE_IJMatrixGetObject(amsD1_, (void **) &D1_csr);
2767 HYPRE_ParCSRMatrixPrint(D1_csr, "D1.parmatrix");
2768 }
2769 }
2770
2771 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
2772 printf("%4d : HYPRE_LSC::leaving putNodalFieldData.\n",mypid_);
2773 return (0);
2774 }
2775
2776 #if 0
2777 /* ------------------------------------------------------ */
2778 /* This section has been replaced due to changes in ale3d */
2779 /* ------------------------------------------------------ */
2780 int HYPRE_LinSysCore::putNodalFieldData(int fieldID, int fieldSize,
2781 int* nodeNumbers, int numNodes, const double* data)
2782 {
2783 int i, **nodeFieldIDs, nodeFieldID, numFields, *procNRows;
2784 int blockID, *blockIDs, *eqnNumbers, *iArray, checkFieldSize;
2785 int *aleNodeNumbers;
2786
2787 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 )
2788 {
2789 printf("%4d : HYPRE_LSC::entering putNodalFieldData.\n",mypid_);
2790 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
2791 {
2792 printf(" putNodalFieldData : fieldSize = %d\n", fieldSize);
2793 printf(" putNodalFieldData : fieldID = %d\n", fieldID);
2794 printf(" putNodalFieldData : numNodes = %d\n", numNodes);
2795 }
2796 }
2797
2798 //-------------------------------------------------------------------
2799 // This part is for loading the nodal coordinate information.
2800 // The node IDs in nodeNumbers are the one used in FEI (and thus
2801 // corresponds to the ones in the system matrix using lookup)
2802 //-------------------------------------------------------------------
2803
2804 if ( fieldID == -3 || fieldID == -25333 )
2805 {
2806 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
2807 {
2808 for ( int i = 0; i < numNodes; i++ )
2809 for ( int j = 0; j < fieldSize; j++ )
2810 printf("putNodalFieldData : %4d %2d = %e\n",i,j,
2811 data[i*fieldSize+j]);
2812 }
2813 if ( HYPreconID_ == HYMLI && lookup_ != NULL )
2814 {
2815 blockIDs = (int *) lookup_->getElemBlockIDs();
2816 blockID = blockIDs[0];
2817 nodeFieldIDs = (int **) lookup_->getFieldIDsTable(blockID);
2818 nodeFieldID = nodeFieldIDs[0][0];
2819 checkFieldSize = lookup_->getFieldSize(nodeFieldID);
2820 //hypre_assert( checkFieldSize == fieldSize );
2821 eqnNumbers = new int[numNodes];
2822 for ( i = 0; i < numNodes; i++ )
2823 eqnNumbers[i] = lookup_->getEqnNumber(nodeNumbers[i],nodeFieldID);
2824 HYPRE_LSI_MLILoadNodalCoordinates(HYPrecon_, numNodes, checkFieldSize,
2825 eqnNumbers, fieldSize, (double *) data);
2826 delete [] eqnNumbers;
2827 }
2828 }
2829
2830 //-------------------------------------------------------------------
2831 // this is needed to set up the correct node equation map
2832 // (the FEI remaps the node IDs in the incoming nodeNumbers array.
2833 // to revert to the original ALE3D node numbers, it is passed in
2834 // here as data)
2835 //-------------------------------------------------------------------
2836
2837 else if ( fieldID == -49773 )
2838 {
2839 if ( HYPreconID_ == HYMLI && lookup_ != NULL )
2840 {
2841 blockIDs = (int *) lookup_->getElemBlockIDs();
2842 blockID = blockIDs[0];
2843 nodeFieldIDs = (int **) lookup_->getFieldIDsTable(blockID);
2844 nodeFieldID = nodeFieldIDs[0][0];
2845 checkFieldSize = lookup_->getFieldSize(nodeFieldID);
2846 hypre_assert( fieldSize == 1 );
2847 aleNodeNumbers = new int[numNodes];
2848 eqnNumbers = new int[numNodes];
2849 for ( i = 0; i < numNodes; i++ )
2850 {
2851 aleNodeNumbers[i] = (int) data[i];
2852 eqnNumbers[i] = lookup_->getEqnNumber(nodeNumbers[i],nodeFieldID);
2853 }
2854 procNRows = new int[numProcs_];
2855 for ( i = 0; i < numProcs_; i++ ) procNRows[i] = 0;
2856 procNRows[mypid_] = localEndRow_;
2857 iArray = procNRows;
2858 procNRows = new int[numProcs_+1];
2859 for ( i = 0; i <= numProcs_; i++ ) procNRows[i] = 0;
2860 MPI_Allreduce(iArray,&(procNRows[1]),numProcs_,MPI_INT,MPI_SUM,
2861 comm_);
2862 delete [] iArray;
2863 HYPRE_LSI_MLICreateNodeEqnMap(HYPrecon_, numNodes, aleNodeNumbers,
2864 eqnNumbers, procNRows);
2865 delete [] procNRows;
2866 delete [] eqnNumbers;
2867 delete [] aleNodeNumbers;
2868 }
2869 }
2870 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 )
2871 printf("%4d : HYPRE_LSC::leaving putNodalFieldData.\n",mypid_);
2872 return (0);
2873 }
2874 #endif
2875
2876 //***************************************************************************
2877 // This function must enforce an essential boundary condition on each local
2878 // equation in 'globalEqn'. This means, that the following modification
2879 // should be made to A and b, for each globalEqn:
2880 //
2881 // for(each local equation i){
2882 // for(each column j in row i) {
2883 // if (i==j) b[i] = gamma/alpha;
2884 // else b[j] -= (gamma/alpha)*A[j,i];
2885 // }
2886 // }
2887 // all of row 'globalEqn' and column 'globalEqn' in A should be zeroed,
2888 // except for 1.0 on the diagonal.
2889 //---------------------------------------------------------------------------
2890
enforceEssentialBC(int * globalEqn,double * alpha,double * gamma1,int leng)2891 int HYPRE_LinSysCore::enforceEssentialBC(int* globalEqn, double* alpha,
2892 double* gamma1, int leng)
2893 {
2894 int i, j, k, localEqnNum, colIndex, rowSize, *colInd, *iarray;
2895 int numLocalRows, eqnNum, rowSize2, *colInd2, numLabels, *labels;
2896 int **i2array, count;
2897 double rhs_term, val, *colVal2, *colVal, **d2array;
2898
2899 //-------------------------------------------------------------------
2900 // diagnostic message and error checking
2901 // (this function should be called before matrixLoadComplete)
2902 //-------------------------------------------------------------------
2903
2904 if ( (HYOutputLevel_ & HYFEI_IMPOSENOBC) != 0 ) return 0;
2905 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
2906 printf("%4d : HYPRE_LSC::entering enforceEssentialBC.\n",mypid_);
2907 if ( systemAssembled_ )
2908 {
2909 printf("enforceEssentialBC ERROR : system assembled already.\n");
2910 exit(1);
2911 }
2912
2913 //-------------------------------------------------------------------
2914 // if matrix partitioning is requested, do it.
2915 //-------------------------------------------------------------------
2916
2917 numLocalRows = localEndRow_ - localStartRow_ + 1;
2918 if ( matrixPartition_ == 1 && HYPreconID_ == HYMLI )
2919 {
2920 HYPRE_LSI_PartitionMatrix(numLocalRows,localStartRow_,rowLengths_,
2921 colIndices_,colValues_, &numLabels, &labels);
2922 HYPRE_LSI_MLILoadMaterialLabels(HYPrecon_, numLabels, labels);
2923 free( labels );
2924 matrixPartition_ = 2;
2925 }
2926
2927 //-------------------------------------------------------------------
2928 // examine each row individually
2929 //-------------------------------------------------------------------
2930
2931 //**/================================================================
2932 //**/ The following is for multiple right hand side (Mar 2009)
2933 if (mRHSFlag_ == 1 && currentRHS_ != 0 && mRHSNumGEqns_ > 0)
2934 {
2935 for( i = 0; i < leng; i++ )
2936 {
2937 for ( j = 0; j < mRHSNumGEqns_; j++ )
2938 if (mRHSGEqnIDs_[j] == globalEqn[i] && mRHSBCType_[j] == 1) break;
2939 if (j == mRHSNumGEqns_)
2940 {
2941 printf("%4d : HYPRE_LSC::enforceEssentialBC ERROR (1).\n",mypid_);
2942 return -1;
2943 }
2944 k = j;
2945 localEqnNum = globalEqn[i] + 1 - localStartRow_;
2946 if ( localEqnNum >= 0 && localEqnNum < numLocalRows )
2947 {
2948 for ( j = 0; j < mRHSNEntries_[k]; j++ )
2949 {
2950 rhs_term = gamma1[i] / alpha[i] * mRHSRowVals_[k][j];
2951 eqnNum = mRHSRowInds_[k][j] - 1;
2952 HYPRE_IJVectorGetValues(HYb_,1, &eqnNum, &val);
2953 val -= rhs_term;
2954 HYPRE_IJVectorSetValues(HYb_, 1, (const int *) &eqnNum,
2955 (const double *) &rhs_term);
2956 }
2957 }
2958 // Set rhs for boundary point
2959 rhs_term = gamma1[i] / alpha[i];
2960 eqnNum = globalEqn[i];
2961 HYPRE_IJVectorSetValues(HYb_, 1, (const int *) &eqnNum,
2962 (const double *) &rhs_term);
2963 }
2964 }
2965 //**/================================================================
2966 else
2967 {
2968 //**/=============================================================
2969 //**/ save the BC information
2970 if (mRHSFlag_ == 1)
2971 {
2972 if (mRHSNumGEqns_ == 0)
2973 {
2974 mRHSGEqnIDs_ = new int[leng];
2975 mRHSNEntries_ = new int[leng];
2976 mRHSBCType_ = new int[leng];
2977 mRHSRowInds_ = new int*[leng];
2978 mRHSRowVals_ = new double*[leng];
2979 for (i = 0; i < leng; i++) mRHSRowInds_[i] = NULL;
2980 for (i = 0; i < leng; i++) mRHSRowVals_[i] = NULL;
2981 }
2982 else
2983 {
2984 iarray = mRHSGEqnIDs_;
2985 mRHSGEqnIDs_ = new int[mRHSNumGEqns_+leng];
2986 for (i = 0; i < mRHSNumGEqns_; i++) mRHSGEqnIDs_[i] = iarray[i];
2987 iarray = mRHSNEntries_;
2988 mRHSNEntries_ = new int[mRHSNumGEqns_+leng];
2989 for (i = 0; i < mRHSNumGEqns_; i++) mRHSNEntries_[i] = iarray[i];
2990 iarray = mRHSBCType_;
2991 mRHSBCType_ = new int[mRHSNumGEqns_+leng];
2992 for (i = 0; i < mRHSNumGEqns_; i++) mRHSBCType_[i] = iarray[i];
2993 i2array = mRHSRowInds_;
2994 mRHSRowInds_ = new int*[mRHSNumGEqns_+leng];
2995 for (i = 0; i < mRHSNumGEqns_; i++) mRHSRowInds_[i] = i2array[i];
2996 d2array = mRHSRowVals_;
2997 for (i = 0; i < mRHSNumGEqns_; i++) mRHSRowInds_[i] = i2array[i];
2998 for (i = mRHSNumGEqns_; i < mRHSNumGEqns_+leng; i++)
2999 mRHSRowInds_[i] = NULL;
3000 mRHSRowVals_ = new double*[mRHSNumGEqns_+leng];
3001 for (i = 0; i < mRHSNumGEqns_; i++) mRHSRowVals_[i] = d2array[i];
3002 for (i = mRHSNumGEqns_; i < mRHSNumGEqns_+leng; i++)
3003 mRHSRowVals_[i] = NULL;
3004 }
3005 }
3006 //**/=============================================================
3007 for( i = 0; i < leng; i++ )
3008 {
3009 localEqnNum = globalEqn[i] + 1 - localStartRow_;
3010 if ( localEqnNum >= 0 && localEqnNum < numLocalRows )
3011 {
3012 rowSize = rowLengths_[localEqnNum];
3013 colInd = colIndices_[localEqnNum];
3014 colVal = colValues_[localEqnNum];
3015
3016 //===================================================
3017 // store the information for multiple right hand side
3018 if (mRHSFlag_ == 1)
3019 {
3020 count = 0;
3021 for ( j = 0; j < rowSize; j++ )
3022 {
3023 colIndex = colInd[j];
3024 if (colIndex >= localStartRow_ && colIndex <= localEndRow_)
3025 {
3026 if ( (colIndex-1) != globalEqn[i])
3027 {
3028 rowSize2 = rowLengths_[colIndex-localStartRow_];
3029 colInd2 = colIndices_[colIndex-localStartRow_];
3030 colVal2 = colValues_ [colIndex-localStartRow_];
3031 for( k = 0; k < rowSize2; k++ )
3032 {
3033 if (colInd2[k]-1 == globalEqn[i])
3034 count++;
3035 break;
3036 }
3037 }
3038 }
3039 }
3040 if (count > 0)
3041 {
3042 mRHSBCType_[mRHSNumGEqns_] = 1;
3043 mRHSGEqnIDs_[mRHSNumGEqns_] = globalEqn[i];
3044 mRHSNEntries_[mRHSNumGEqns_] = count;
3045 mRHSRowInds_[mRHSNumGEqns_] = new int[count];
3046 mRHSRowVals_[mRHSNumGEqns_] = new double[count];
3047 }
3048 count = 0;
3049 for ( j = 0; j < rowSize; j++ )
3050 {
3051 colIndex = colInd[j];
3052 if (colIndex >= localStartRow_ && colIndex <= localEndRow_)
3053 {
3054 if ( (colIndex-1) != globalEqn[i])
3055 {
3056 rowSize2 = rowLengths_[colIndex-localStartRow_];
3057 colInd2 = colIndices_[colIndex-localStartRow_];
3058 colVal2 = colValues_ [colIndex-localStartRow_];
3059 for( k = 0; k < rowSize2; k++ )
3060 {
3061 if ( colInd2[k]-1 == globalEqn[i] )
3062 {
3063 mRHSRowVals_[mRHSNumGEqns_][count] = colVal2[k];
3064 mRHSRowInds_[mRHSNumGEqns_][count] = colIndex;
3065 count++;
3066 break;
3067 }
3068 }
3069 }
3070 }
3071 }
3072 mRHSNumGEqns_++;
3073 }
3074 //===================================================
3075
3076 for ( j = 0; j < rowSize; j++ )
3077 {
3078 colIndex = colInd[j];
3079 if ( colIndex-1 == globalEqn[i] ) colVal[j] = 1.0;
3080 else colVal[j] = 0.0;
3081 if ( colIndex >= localStartRow_ && colIndex <= localEndRow_)
3082 {
3083 if ( (colIndex-1) != globalEqn[i])
3084 {
3085 rowSize2 = rowLengths_[colIndex-localStartRow_];
3086 colInd2 = colIndices_[colIndex-localStartRow_];
3087 colVal2 = colValues_ [colIndex-localStartRow_];
3088
3089 for( k = 0; k < rowSize2; k++ )
3090 {
3091 if ( colInd2[k]-1 == globalEqn[i] )
3092 {
3093 rhs_term = gamma1[i] / alpha[i] * colVal2[k];
3094 eqnNum = colIndex - 1;
3095 HYPRE_IJVectorGetValues(HYb_,1,&eqnNum, &val);
3096 val -= rhs_term;
3097 HYPRE_IJVectorSetValues(HYb_,1,(const int *) &eqnNum,
3098 (const double *) &val);
3099 colVal2[k] = 0.0;
3100 break;
3101 }
3102 }
3103 }
3104 }
3105 }// end for(j<rowSize) loop
3106
3107 // Set rhs for boundary point
3108 rhs_term = gamma1[i] / alpha[i];
3109 eqnNum = globalEqn[i];
3110 HYPRE_IJVectorSetValues(HYb_, 1, (const int *) &eqnNum,
3111 (const double *) &rhs_term);
3112
3113 }
3114 }
3115 }
3116
3117 //-------------------------------------------------------------------
3118 // set up the AMGe Dirichlet boundary conditions
3119 //-------------------------------------------------------------------
3120
3121 #ifdef HAVE_AMGE
3122 colInd = new int[leng];
3123 for( i = 0; i < leng; i++ ) colInd[i] = globalEqn[i];
3124 HYPRE_LSI_AMGeSetBoundary( leng, colInd );
3125 delete [] colInd;
3126 #endif
3127
3128 //-------------------------------------------------------------------
3129 // diagnostic message
3130 //-------------------------------------------------------------------
3131
3132 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
3133 printf("%4d : HYPRE_LSC::leaving enforceEssentialBC.\n",mypid_);
3134 return (0);
3135 }
3136
3137 //***************************************************************************
3138 // put in globalEqns should hold eqns that are owned locally, but which contain
3139 // column indices (the ones in colIndices) which are from remote equations
3140 // on which essential boundary-conditions need to be enforced.
3141 // This function will only make the modification if the above conditions
3142 // hold -- i.e., the equation is a locally-owned equation, and the column
3143 // index is NOT a locally owned equation.
3144 //---------------------------------------------------------------------------
3145
enforceRemoteEssBCs(int numEqns,int * globalEqns,int ** colIndices,int * colIndLen,double ** coefs)3146 int HYPRE_LinSysCore::enforceRemoteEssBCs(int numEqns, int* globalEqns,
3147 int** colIndices, int* colIndLen,
3148 double** coefs)
3149 {
3150 int i, j, k, numLocalRows, localEqnNum, rowLen, *colInd, eqnNum;
3151 int *iarray, **i2array, count;
3152 double bval, *colVal, rhs_term, **d2array;
3153
3154 //-------------------------------------------------------------------
3155 // diagnostic message and error checking
3156 // (this function should be called before matrixLoadComplete)
3157 //-------------------------------------------------------------------
3158
3159 if ( (HYOutputLevel_ & HYFEI_IMPOSENOBC) != 0 ) return 0;
3160 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
3161 printf("%4d : HYPRE_LSC::entering enforceRemoteEssBC.\n",mypid_);
3162 if ( systemAssembled_ )
3163 {
3164 printf("enforceRemoteEssBC ERROR : system assembled already.\n");
3165 exit(1);
3166 }
3167
3168 //-------------------------------------------------------------------
3169 // examine each row individually
3170 //-------------------------------------------------------------------
3171
3172 numLocalRows = localEndRow_ - localStartRow_ + 1;
3173
3174 //============================================================
3175 //**/ The following is for multiple right hand side (Mar 2009)
3176
3177 if (mRHSFlag_ == 1 && currentRHS_ != 0 && mRHSNumGEqns_ > 0)
3178 {
3179 for( i = 0; i < numEqns; i++ )
3180 {
3181 for ( j = 0; j < mRHSNumGEqns_; j++ )
3182 if (mRHSGEqnIDs_[j] == globalEqns[i] && mRHSBCType_[j] == 2) break;
3183 if (j == mRHSNumGEqns_)
3184 {
3185 printf("%4d : HYPRE_LSC::enforceRemoteEssBCs ERROR (1).\n",mypid_);
3186 return -1;
3187 }
3188 k = j;
3189 localEqnNum = globalEqns[i] + 1 - localStartRow_;
3190 if ( localEqnNum < 0 || localEqnNum >= numLocalRows )
3191 {
3192 continue;
3193 }
3194 eqnNum = globalEqns[i];
3195 rowLen = mRHSNEntries_[k];
3196 colInd = mRHSRowInds_[k];
3197 colVal = mRHSRowVals_[k];
3198 for ( j = 0; j < colIndLen[i]; j++)
3199 {
3200 for ( k = 0; k < rowLen; k++ )
3201 {
3202 if (colInd[k]-1 == colIndices[i][j])
3203 {
3204 rhs_term = colVal[k] * coefs[i][j];
3205 HYPRE_IJVectorGetValues(HYb_,1,&eqnNum,&bval);
3206 bval -= rhs_term;
3207 HYPRE_IJVectorSetValues(HYb_, 1, (const int *) &eqnNum,
3208 (const double *) &bval);
3209 }
3210 }
3211 }
3212 }
3213 }
3214 //**/================================================================
3215 else
3216 {
3217 //**/=============================================================
3218 //**/ save the BC information
3219 if (mRHSFlag_ == 1)
3220 {
3221 if (mRHSNumGEqns_ == 0)
3222 {
3223 mRHSGEqnIDs_ = new int[numEqns];
3224 mRHSNEntries_ = new int[numEqns];
3225 mRHSBCType_ = new int[numEqns];
3226 mRHSRowInds_ = new int*[numEqns];
3227 mRHSRowVals_ = new double*[numEqns];
3228 for (i = 0; i < numEqns; i++) mRHSRowInds_[i] = NULL;
3229 for (i = 0; i < numEqns; i++) mRHSRowVals_[i] = NULL;
3230 }
3231 else
3232 {
3233 iarray = mRHSGEqnIDs_;
3234 mRHSGEqnIDs_ = new int[mRHSNumGEqns_+numEqns];
3235 for (i = 0; i < mRHSNumGEqns_; i++) mRHSGEqnIDs_[i] = iarray[i];
3236 iarray = mRHSNEntries_;
3237 mRHSNEntries_ = new int[mRHSNumGEqns_+numEqns];
3238 for (i = 0; i < mRHSNumGEqns_; i++) mRHSNEntries_[i] = iarray[i];
3239 iarray = mRHSBCType_;
3240 mRHSBCType_ = new int[mRHSNumGEqns_+numEqns];
3241 for (i = 0; i < mRHSNumGEqns_; i++) mRHSBCType_[i] = iarray[i];
3242 i2array = mRHSRowInds_;
3243 mRHSRowInds_ = new int*[mRHSNumGEqns_+numEqns];
3244 for (i = 0; i < mRHSNumGEqns_; i++) mRHSRowInds_[i] = i2array[i];
3245 d2array = mRHSRowVals_;
3246 for (i = 0; i < mRHSNumGEqns_; i++) mRHSRowInds_[i] = i2array[i];
3247 for (i = mRHSNumGEqns_; i < mRHSNumGEqns_+numEqns; i++)
3248 mRHSRowInds_[i] = NULL;
3249 mRHSRowVals_ = new double*[mRHSNumGEqns_+numEqns];
3250 for (i = 0; i < mRHSNumGEqns_; i++) mRHSRowVals_[i] = d2array[i];
3251 for (i = mRHSNumGEqns_; i < mRHSNumGEqns_+numEqns; i++)
3252 mRHSRowVals_[i] = NULL;
3253 }
3254 }
3255 //**/=============================================================
3256
3257 for( i = 0; i < numEqns; i++ )
3258 {
3259 localEqnNum = globalEqns[i] + 1 - localStartRow_;
3260 if ( localEqnNum < 0 || localEqnNum >= numLocalRows )
3261 {
3262 continue;
3263 }
3264
3265 rowLen = rowLengths_[localEqnNum];
3266 colInd = colIndices_[localEqnNum];
3267 colVal = colValues_[localEqnNum];
3268 eqnNum = globalEqns[i];
3269
3270 //===================================================
3271 // store the information for multiple right hand side
3272 if (mRHSFlag_ == 1)
3273 {
3274 count = 0;
3275 for ( j = 0; j < colIndLen[i]; j++)
3276 {
3277 for ( k = 0; k < rowLen; k++ )
3278 if (colInd[k]-1 == colIndices[i][j]) count++;
3279 }
3280 if (count > 0)
3281 {
3282 mRHSGEqnIDs_[mRHSNumGEqns_] = globalEqns[i];
3283 mRHSBCType_[mRHSNumGEqns_] = 2;
3284 mRHSNEntries_[mRHSNumGEqns_] = count;
3285 mRHSRowInds_[mRHSNumGEqns_] = new int[count];
3286 mRHSRowVals_[mRHSNumGEqns_] = new double[count];
3287 }
3288 count = 0;
3289 for ( j = 0; j < colIndLen[i]; j++)
3290 {
3291 for ( k = 0; k < rowLen; k++ )
3292 {
3293 if (colInd[k]-1 == colIndices[i][j])
3294 {
3295 mRHSRowVals_[k][count] = colVal[k];
3296 mRHSRowInds_[k][count] = colInd[k];
3297 }
3298 }
3299 }
3300 }
3301 //===================================================
3302
3303 for ( j = 0; j < colIndLen[i]; j++)
3304 {
3305 for ( k = 0; k < rowLen; k++ )
3306 {
3307 if (colInd[k]-1 == colIndices[i][j])
3308 {
3309 rhs_term = colVal[k] * coefs[i][j];
3310 HYPRE_IJVectorGetValues(HYb_,1,&eqnNum,&bval);
3311 bval -= rhs_term;
3312 HYPRE_IJVectorSetValues(HYb_, 1, (const int *) &eqnNum,
3313 (const double *) &bval);
3314 colVal[k] = 0.0;
3315 }
3316 }
3317 }
3318 }
3319 }
3320
3321 //-------------------------------------------------------------------
3322 // diagnostic message
3323 //-------------------------------------------------------------------
3324
3325 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
3326 printf("%4d : HYPRE_LSC::leaving enforceRemoteEssBC.\n",mypid_);
3327 return (0);
3328 }
3329
3330 //***************************************************************************
3331 //This function must enforce a natural or mixed boundary condition on the
3332 //equations in 'globalEqn'. This means that the following modification should
3333 //be made to A and b:
3334 //
3335 //A[globalEqn,globalEqn] += alpha/beta;
3336 //b[globalEqn] += gamma/beta;
3337 //---------------------------------------------------------------------------
3338
enforceOtherBC(int * globalEqn,double * alpha,double * beta,double * gamma1,int leng)3339 int HYPRE_LinSysCore::enforceOtherBC(int* globalEqn, double* alpha,
3340 double* beta, double* gamma1, int leng)
3341 {
3342 int i, j, numLocalRows, localEqnNum, *colInd, rowSize, eqnNum;
3343 double val, *colVal, rhs_term;
3344
3345 //-------------------------------------------------------------------
3346 // diagnostic message and error checking
3347 // (this function should be called before matrixLoadComplete)
3348 //-------------------------------------------------------------------
3349
3350 if ( (HYOutputLevel_ & HYFEI_IMPOSENOBC) != 0 ) return 0;
3351 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
3352 printf("%4d : HYPRE_LSC::entering enforceOtherBC.\n",mypid_);
3353 if ( systemAssembled_ )
3354 {
3355 printf("enforceOtherBC ERROR : system assembled already.\n");
3356 exit(1);
3357 }
3358
3359 //-------------------------------------------------------------------
3360 // examine each row individually
3361 //-------------------------------------------------------------------
3362
3363 numLocalRows = localEndRow_ - localStartRow_ + 1;
3364
3365 //============================================================
3366 //**/ The following is for multiple right hand side (Mar 2009)
3367 if (mRHSFlag_ == 1 && currentRHS_ != 0)
3368 {
3369 for( i = 0; i < leng; i++ )
3370 {
3371 localEqnNum = globalEqn[i] + 1 - localStartRow_;
3372 if ( localEqnNum < 0 || localEqnNum >= numLocalRows )
3373 {
3374 continue;
3375 }
3376 eqnNum = globalEqn[i];
3377 rhs_term = gamma1[i] / beta[i];
3378 HYPRE_IJVectorGetValues(HYb_,1,&eqnNum,&val);
3379 val += rhs_term;
3380 HYPRE_IJVectorSetValues(HYb_, 1, (const int *) &eqnNum,
3381 (const double *) &val);
3382 }
3383 }
3384 //============================================================
3385 else
3386 {
3387 for( i = 0; i < leng; i++ )
3388 {
3389 localEqnNum = globalEqn[i] + 1 - localStartRow_;
3390 if ( localEqnNum < 0 || localEqnNum >= numLocalRows )
3391 {
3392 continue;
3393 }
3394
3395 rowSize = rowLengths_[localEqnNum];
3396 colVal = colValues_[localEqnNum];
3397 colInd = colIndices_[localEqnNum];
3398
3399 for ( j = 0; j < rowSize; j++)
3400 {
3401 if ((colInd[j]-1) == globalEqn[i])
3402 {
3403 colVal[j] += alpha[i]/beta[i];
3404 break;
3405 }
3406 }
3407
3408 //now make the rhs modification.
3409 // need to fetch matrix and put it back before assembled
3410
3411 eqnNum = globalEqn[i];
3412 rhs_term = gamma1[i] / beta[i];
3413 HYPRE_IJVectorGetValues(HYb_,1,&eqnNum,&val);
3414 val += rhs_term;
3415 HYPRE_IJVectorSetValues(HYb_, 1, (const int *) &eqnNum,
3416 (const double *) &val);
3417 }
3418 }
3419
3420 //-------------------------------------------------------------------
3421 // diagnostic message
3422 //-------------------------------------------------------------------
3423
3424 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
3425 printf("%4d : HYPRE_LSC::leaving enforceOtherBC.\n",mypid_);
3426 return (0);
3427 }
3428
3429 //***************************************************************************
3430 // put the pointer to the A matrix into the Data object
3431 //---------------------------------------------------------------------------
3432
3433 #ifndef NOFEI
getMatrixPtr(Data & data)3434 int HYPRE_LinSysCore::getMatrixPtr(Data& data)
3435 {
3436 (void) data;
3437 printf("%4d : HYPRE_LSC::getMatrixPtr ERROR - not implemented.\n",mypid_);
3438 exit(1);
3439 return (0);
3440 }
3441 #endif
3442
3443 //***************************************************************************
3444 // Overwrites the current internal matrix with a scaled copy of the
3445 // input argument.
3446 //---------------------------------------------------------------------------
3447
3448 #ifndef NOFEI
copyInMatrix(double scalar,const Data & data)3449 int HYPRE_LinSysCore::copyInMatrix(double scalar, const Data& data)
3450 {
3451 int i;
3452 char *name;
3453 HYPRE_FEI_AMSData *auxAMSData;
3454
3455 (void) scalar;
3456
3457 name = data.getTypeName();
3458 if (!strcmp(name, "ANN"))
3459 {
3460 maxwellANN_ = (HYPRE_ParCSRMatrix) data.getDataPtr();
3461 }
3462 else if (!strcmp(name, "GEN"))
3463 {
3464 maxwellGEN_ = (HYPRE_ParCSRMatrix) data.getDataPtr();
3465 }
3466 else if (!strcmp(name, "AMSBMATRIX"))
3467 {
3468 amsBetaPoisson_ = (HYPRE_ParCSRMatrix) data.getDataPtr();
3469 }
3470 else if (!strcmp(name, "AMSData"))
3471 {
3472 auxAMSData = (HYPRE_FEI_AMSData *) data.getDataPtr();
3473 if (AMSData_.NodeNumbers_ != NULL) hypre_TFree(AMSData_.NodeNumbers_, HYPRE_MEMORY_HOST);
3474 if (AMSData_.NodalCoord_ != NULL) hypre_TFree(AMSData_.NodalCoord_, HYPRE_MEMORY_HOST);
3475 AMSData_.NodeNumbers_ = NULL;
3476 AMSData_.NodalCoord_ = NULL;
3477 AMSData_.numNodes_ = auxAMSData->numNodes_;
3478 AMSData_.numLocalNodes_ = auxAMSData->numLocalNodes_;
3479 if (AMSData_.numNodes_ > 0)
3480 {
3481 AMSData_.NodeNumbers_ = hypre_CTAlloc(HYPRE_BigInt, AMSData_.numNodes_, HYPRE_MEMORY_HOST);
3482 AMSData_.NodalCoord_ = hypre_CTAlloc(HYPRE_Real, AMSData_.numNodes_*mlNumPDEs_, HYPRE_MEMORY_HOST);
3483 for (i = 0; i < AMSData_.numNodes_; i++)
3484 AMSData_.NodeNumbers_[i] = auxAMSData->NodeNumbers_[i];
3485 for (i = 0; i < AMSData_.numNodes_*mlNumPDEs_; i++)
3486 AMSData_.NodalCoord_[i] = auxAMSData->NodalCoord_[i];
3487 }
3488 }
3489 else
3490 {
3491 printf("%4d : HYPRE_LSC::copyInMatrix ERROR - invalid data.\n",mypid_);
3492 exit(1);
3493 }
3494 return (0);
3495 }
3496 #endif
3497
3498 //***************************************************************************
3499 //Passes out a scaled copy of the current internal matrix.
3500 //---------------------------------------------------------------------------
3501
3502 #ifndef NOFEI
copyOutMatrix(double scalar,Data & data)3503 int HYPRE_LinSysCore::copyOutMatrix(double scalar, Data& data)
3504 {
3505 char *name;
3506
3507 (void) scalar;
3508
3509 name = data.getTypeName();
3510
3511 if (!strcmp(name, "A"))
3512 {
3513 data.setDataPtr((void *) HYA_);
3514 }
3515 else if (!strcmp(name, "AMSData"))
3516 {
3517 data.setDataPtr((void *) &AMSData_);
3518 }
3519 else
3520 {
3521 printf("HYPRE_LSC::copyOutMatrix ERROR - invalid command.\n");
3522 exit(1);
3523 }
3524 return (0);
3525 }
3526 #endif
3527
3528 //***************************************************************************
3529 // add nonzero entries into the matrix data structure
3530 //---------------------------------------------------------------------------
3531
3532 #ifndef NOFEI
sumInMatrix(double scalar,const Data & data)3533 int HYPRE_LinSysCore::sumInMatrix(double scalar, const Data& data)
3534 {
3535 (void) scalar;
3536 (void) data;
3537 printf("%4d : HYPRE_LSC::sumInMatrix ERROR - not implemented.\n",mypid_);
3538 exit(1);
3539 return (0);
3540 }
3541 #endif
3542
3543 //***************************************************************************
3544 // get the data pointer for the right hand side
3545 //---------------------------------------------------------------------------
3546
3547 #ifndef NOFEI
getRHSVectorPtr(Data & data)3548 int HYPRE_LinSysCore::getRHSVectorPtr(Data& data)
3549 {
3550 //-------------------------------------------------------------------
3551 // diagnostic message
3552 //-------------------------------------------------------------------
3553
3554 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3555 printf("%4d : HYPRE_LSC::entering getRHSVectorPtr.\n",mypid_);
3556
3557 //-------------------------------------------------------------------
3558 // get the right hand side vector pointer
3559 //-------------------------------------------------------------------
3560
3561 data.setTypeName("IJ_Vector");
3562 data.setDataPtr((void*) HYb_);
3563
3564 //-------------------------------------------------------------------
3565 // diagnostic message
3566 //-------------------------------------------------------------------
3567
3568 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3569 printf("%4d : HYPRE_LSC::leaving getRHSVectorPtr.\n",mypid_);
3570 return (0);
3571 }
3572 #endif
3573
3574 //***************************************************************************
3575 // copy the content of the incoming vector to the right hand side vector
3576 //---------------------------------------------------------------------------
3577
3578 #ifndef NOFEI
copyInRHSVector(double scalar,const Data & data)3579 int HYPRE_LinSysCore::copyInRHSVector(double scalar, const Data& data)
3580 {
3581 //-------------------------------------------------------------------
3582 // diagnostic message
3583 //-------------------------------------------------------------------
3584
3585 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3586 printf("%4d : HYPRE_LSC::entering copyInRHSVector.\n",mypid_);
3587 if (strcmp("IJ_Vector", data.getTypeName()) &&
3588 strcmp("Sol_Vector", data.getTypeName()))
3589 {
3590 printf("copyInRHSVector: data's type string not compatible.\n");
3591 exit(1);
3592 }
3593
3594 //-------------------------------------------------------------------
3595 // copy the incoming vector to the internal right hand side
3596 //-------------------------------------------------------------------
3597
3598 HYPRE_IJVector inVec = (HYPRE_IJVector) data.getDataPtr();
3599 HYPRE_ParVector srcVec;
3600 HYPRE_ParVector destVec;
3601 HYPRE_IJVectorGetObject(inVec, (void **) &srcVec);
3602 if (!strcmp("Sol_Vector", data.getTypeName()))
3603 HYPRE_IJVectorGetObject(HYb_, (void **) &destVec);
3604 else
3605 HYPRE_IJVectorGetObject(HYx_, (void **) &destVec);
3606
3607 HYPRE_ParVectorCopy(srcVec, destVec);
3608
3609 if ( scalar != 1.0 ) HYPRE_ParVectorScale( scalar, destVec);
3610 // do not destroy the incoming vector
3611 //HYPRE_IJVectorDestroy(inVec);
3612
3613 //-------------------------------------------------------------------
3614 // diagnostic message
3615 //-------------------------------------------------------------------
3616
3617 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3618 printf("%4d : HYPRE_LSC::leaving copyInRHSVector.\n",mypid_);
3619 return (0);
3620 }
3621 #endif
3622
3623 //***************************************************************************
3624 // create an ParVector and copy the right hand side to it (scaled)
3625 //---------------------------------------------------------------------------
3626
3627 #ifndef NOFEI
copyOutRHSVector(double scalar,Data & data)3628 int HYPRE_LinSysCore::copyOutRHSVector(double scalar, Data& data)
3629 {
3630 //-------------------------------------------------------------------
3631 // diagnostic message
3632 //-------------------------------------------------------------------
3633
3634 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3635 printf("%4d : HYPRE_LSC::entering copyOutRHSVector.\n",mypid_);
3636
3637 //-------------------------------------------------------------------
3638 // extract the right hand side vector
3639 //-------------------------------------------------------------------
3640
3641 HYPRE_IJVector newVector;
3642 HYPRE_IJVectorCreate(comm_, localStartRow_-1, localEndRow_-1,
3643 &newVector);
3644 HYPRE_IJVectorSetObjectType(newVector, HYPRE_PARCSR);
3645 HYPRE_IJVectorInitialize(newVector);
3646 HYPRE_IJVectorAssemble(newVector);
3647
3648 HYPRE_ParVector Vec1;
3649 HYPRE_ParVector Vec2;
3650 HYPRE_IJVectorGetObject(HYb_, (void **) &Vec1);
3651 HYPRE_IJVectorGetObject(newVector, (void **) &Vec2);
3652
3653 HYPRE_ParVectorCopy( Vec1, Vec2);
3654 if ( scalar != 1.0 ) HYPRE_ParVectorScale( scalar, Vec2);
3655
3656 data.setTypeName("IJ_Vector");
3657 data.setDataPtr((void*) Vec2);
3658
3659 //-------------------------------------------------------------------
3660 // diagnostic message
3661 //-------------------------------------------------------------------
3662
3663 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3664 printf("%4d : HYPRE_LSC::leaving copyOutRHSVector.\n",mypid_);
3665 return (0);
3666 }
3667 #endif
3668
3669 //***************************************************************************
3670 // add the incoming ParCSR vector to the current right hand side (scaled)
3671 //---------------------------------------------------------------------------
3672
3673 #ifndef NOFEI
sumInRHSVector(double scalar,const Data & data)3674 int HYPRE_LinSysCore::sumInRHSVector(double scalar, const Data& data)
3675 {
3676 //-------------------------------------------------------------------
3677 // diagnostic message
3678 //-------------------------------------------------------------------
3679
3680 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3681 printf("%4d : HYPRE_LSC::entering sumInRHSVector.\n",mypid_);
3682 if (strcmp("IJ_Vector", data.getTypeName()))
3683 {
3684 printf("sumInRHSVector ERROR : data's type string not 'IJ_Vector'.\n");
3685 exit(1);
3686 }
3687
3688 //-------------------------------------------------------------------
3689 // add the incoming vector to the right hand side
3690 //-------------------------------------------------------------------
3691
3692 HYPRE_IJVector inVec = (HYPRE_IJVector) data.getDataPtr();
3693 HYPRE_ParVector xVec;
3694 HYPRE_ParVector yVec;
3695 HYPRE_IJVectorGetObject(inVec, (void **) &xVec);
3696 HYPRE_IJVectorGetObject(HYb_, (void **) &yVec);
3697 hypre_ParVectorAxpy(scalar,(hypre_ParVector*)xVec,(hypre_ParVector*)yVec);
3698
3699 //-------------------------------------------------------------------
3700 // diagnostic message
3701 //-------------------------------------------------------------------
3702
3703 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3704 printf("%4d : HYPRE_LSC::leaving sumInRHSVector.\n",mypid_);
3705 return (0);
3706 }
3707 #endif
3708
3709 //***************************************************************************
3710 // deallocate an incoming IJ matrix
3711 //---------------------------------------------------------------------------
3712
3713 #ifndef NOFEI
destroyMatrixData(Data & data)3714 int HYPRE_LinSysCore::destroyMatrixData(Data& data)
3715 {
3716 //-------------------------------------------------------------------
3717 // diagnostic message
3718 //-------------------------------------------------------------------
3719
3720 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3721 printf("%4d : HYPRE_LSC::entering destroyMatrixData.\n",mypid_);
3722 if (strcmp("IJ_Matrix", data.getTypeName()))
3723 {
3724 printf("destroyMatrixData ERROR : data doesn't contain a IJ_Matrix.\n");
3725 exit(1);
3726 }
3727
3728 //-------------------------------------------------------------------
3729 // destroy the incoming matrix data
3730 //-------------------------------------------------------------------
3731
3732 HYPRE_IJMatrix mat = (HYPRE_IJMatrix) data.getDataPtr();
3733 HYPRE_IJMatrixDestroy(mat);
3734
3735 //-------------------------------------------------------------------
3736 // diagnostic message
3737 //-------------------------------------------------------------------
3738
3739 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3740 printf("%4d : HYPRE_LSC::leaving destroyMatrixData.\n",mypid_);
3741 return (0);
3742 }
3743 #endif
3744
3745 //***************************************************************************
3746 // deallocate an incoming IJ vector
3747 //---------------------------------------------------------------------------
3748
3749 #ifndef NOFEI
destroyVectorData(Data & data)3750 int HYPRE_LinSysCore::destroyVectorData(Data& data)
3751 {
3752 //-------------------------------------------------------------------
3753 // diagnostic message
3754 //-------------------------------------------------------------------
3755
3756 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3757 printf("%4d : HYPRE_LSC::entering destroyVectorData.\n",mypid_);
3758 if (strcmp("IJ_Vector", data.getTypeName()))
3759 {
3760 printf("destroyVectorData ERROR : data doesn't contain a IJ_Vector.");
3761 exit(1);
3762 }
3763
3764 //-------------------------------------------------------------------
3765 // destroy the incoming vector data
3766 //-------------------------------------------------------------------
3767
3768 HYPRE_IJVector vec = (HYPRE_IJVector) data.getDataPtr();
3769 if ( vec != NULL ) HYPRE_IJVectorDestroy(vec);
3770
3771 //-------------------------------------------------------------------
3772 // diagnostic message
3773 //-------------------------------------------------------------------
3774
3775 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3776 printf("%4d : HYPRE_LSC::leaving destroyVectorData.\n",mypid_);
3777 return (0);
3778 }
3779 #endif
3780
3781 //***************************************************************************
3782 // set number of right hand side vectors
3783 //---------------------------------------------------------------------------
3784
setNumRHSVectors(int numRHSs,const int * rhsIDs)3785 int HYPRE_LinSysCore::setNumRHSVectors(int numRHSs, const int* rhsIDs)
3786 {
3787 //-------------------------------------------------------------------
3788 // diagnostic message
3789 //-------------------------------------------------------------------
3790
3791 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3792 {
3793 printf("%4d : HYPRE_LSC::entering setNumRHSVectors.\n",mypid_);
3794 printf("%4d : HYPRE_LSC::incoming numRHSs = %d\n",mypid_,numRHSs);
3795 for ( int i = 0; i < numRHSs_; i++ )
3796 printf("%4d : HYPRE_LSC::incoming RHSIDs = %d\n",mypid_,rhsIDs[i]);
3797 }
3798 if (numRHSs < 0)
3799 {
3800 printf("setNumRHSVectors ERROR : numRHSs < 0.\n");
3801 exit(1);
3802 }
3803
3804 //-------------------------------------------------------------------
3805 // first destroy the existing right hand side vectors
3806 //-------------------------------------------------------------------
3807
3808 if ( matrixVectorsCreated_ )
3809 {
3810 if ( HYbs_ != NULL )
3811 {
3812 for ( int i = 0; i < numRHSs_; i++ )
3813 if ( HYbs_[i] != NULL ) HYPRE_IJVectorDestroy(HYbs_[i]);
3814 delete [] HYbs_;
3815 HYbs_ = NULL;
3816 }
3817 }
3818 if (numRHSs == 0) return (0);
3819
3820 //-------------------------------------------------------------------
3821 // instantiate the right hand vectors
3822 //-------------------------------------------------------------------
3823
3824 if ( matrixVectorsCreated_ )
3825 {
3826 HYbs_ = new HYPRE_IJVector[numRHSs_];
3827 for ( int i = 0; i < numRHSs_; i++ )
3828 {
3829 HYPRE_IJVectorCreate(comm_, localStartRow_-1, localEndRow_-1,
3830 &(HYbs_[i]));
3831 HYPRE_IJVectorSetObjectType(HYbs_[i], HYPRE_PARCSR);
3832 HYPRE_IJVectorInitialize(HYbs_[i]);
3833 HYPRE_IJVectorAssemble(HYbs_[i]);
3834 }
3835 HYb_ = HYbs_[0];
3836 }
3837
3838 //-------------------------------------------------------------------
3839 // copy in the right hand side IDs
3840 //-------------------------------------------------------------------
3841
3842 delete [] rhsIDs_;
3843 numRHSs_ = numRHSs;
3844 rhsIDs_ = new int[numRHSs_];
3845
3846 for ( int i = 0; i < numRHSs; i++ ) rhsIDs_[i] = rhsIDs[i];
3847
3848 //-------------------------------------------------------------------
3849 // diagnostic message
3850 //-------------------------------------------------------------------
3851
3852 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3853 printf("%4d : HYPRE_LSC::leaving setNumRHSVectors.\n",mypid_);
3854 return (0);
3855 }
3856
3857 //***************************************************************************
3858 // select a right hand side vector
3859 //---------------------------------------------------------------------------
3860
setRHSID(int rhsID)3861 int HYPRE_LinSysCore::setRHSID(int rhsID)
3862 {
3863 //-------------------------------------------------------------------
3864 // diagnostic message
3865 //-------------------------------------------------------------------
3866
3867 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3868 printf("%4d : HYPRE_LSC::setRHSID = %d.\n",mypid_,rhsID);
3869
3870 //-------------------------------------------------------------------
3871 // set current right hand side vector ID
3872 //-------------------------------------------------------------------
3873
3874 for( int i = 0; i < numRHSs_; i++ )
3875 {
3876 if (rhsIDs_[i] == rhsID)
3877 {
3878 currentRHS_ = i;
3879 HYb_ = HYbs_[currentRHS_];
3880 currB_ = HYb_;
3881 return (0);
3882 }
3883 }
3884 printf("setRHSID ERROR : rhsID %d not found.\n", rhsID);
3885 exit(1);
3886 return (0);
3887 }
3888
3889 //***************************************************************************
3890 // used for initializing the initial guess
3891 //---------------------------------------------------------------------------
3892
putInitialGuess(const int * eqnNumbers,const double * values,int leng)3893 int HYPRE_LinSysCore::putInitialGuess(const int* eqnNumbers,
3894 const double* values, int leng)
3895 {
3896 int i, *localInds, *iarray, *iarray2;
3897
3898 //-------------------------------------------------------------------
3899 // diagnostic message
3900 //-------------------------------------------------------------------
3901
3902 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3903 printf("%4d : HYPRE_LSC::entering putInitalGuess.\n",mypid_);
3904
3905 //-------------------------------------------------------------------
3906 // this is to create a FEI to HYPRE equation node map
3907 //-------------------------------------------------------------------
3908
3909 if ( mapFromSolnFlag_ == 1 )
3910 {
3911 if ( (mapFromSolnLeng_+leng) >= mapFromSolnLengMax_ )
3912 {
3913 iarray = mapFromSolnList_;
3914 iarray2 = mapFromSolnList2_;
3915 mapFromSolnLengMax_ = mapFromSolnLengMax_ + 2 * leng;
3916 mapFromSolnList_ = new int[mapFromSolnLengMax_];
3917 mapFromSolnList2_ = new int[mapFromSolnLengMax_];
3918 for ( i = 0; i < mapFromSolnLeng_; i++ )
3919 {
3920 mapFromSolnList_[i] = iarray[i];
3921 mapFromSolnList2_[i] = iarray2[i];
3922 }
3923 if ( iarray != NULL ) delete [] iarray;
3924 if ( iarray2 != NULL ) delete [] iarray2;
3925 }
3926 }
3927
3928 localInds = new int[leng];
3929 for ( i = 0; i < leng; i++ ) // change to 0-based
3930 {
3931 if ((eqnNumbers[i]+1) >= localStartRow_ &&
3932 (eqnNumbers[i]+1) <= localEndRow_) localInds[i] = eqnNumbers[i];
3933 else
3934 {
3935 printf("%d : putInitialGuess ERROR - index %d out of range\n",
3936 mypid_, eqnNumbers[i]);
3937 exit(1);
3938 }
3939 if ( mapFromSolnFlag_ == 1 )
3940 {
3941 mapFromSolnList_[mapFromSolnLeng_] = eqnNumbers[i];
3942 mapFromSolnList2_[mapFromSolnLeng_++] = (int) values[i];
3943 }
3944 }
3945 HYPRE_IJVectorSetValues(HYx_, leng, (const int *) localInds,
3946 (const double *) values);
3947 delete [] localInds;
3948
3949 //-------------------------------------------------------------------
3950 // inject the initial guess into reduced systems, if any
3951 //-------------------------------------------------------------------
3952
3953 if ( schurReduction_ == 1 ) buildSchurInitialGuess();
3954
3955 //-------------------------------------------------------------------
3956 // diagnostic message
3957 //-------------------------------------------------------------------
3958
3959 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
3960 printf("%4d : HYPRE_LSC::leaving putInitalGuess.\n",mypid_);
3961 return (0);
3962 }
3963
3964 //***************************************************************************
3965 // This is a modified function for version 1.5
3966 // used for getting the solution out of the solver, and into the application
3967 //---------------------------------------------------------------------------
3968
getSolution(double * answers,int leng)3969 int HYPRE_LinSysCore::getSolution(double* answers,int leng)
3970 {
3971 int i, *equations;
3972
3973 //-------------------------------------------------------------------
3974 // diagnostic message
3975 //-------------------------------------------------------------------
3976
3977 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
3978 printf("%4d : HYPRE_LSC::entering getSolution.\n",mypid_);
3979 if (localStartCol_ == -1 && leng != (localEndRow_-localStartRow_+1))
3980 {
3981 printf("%4d : HYPRE_LSC ERROR : getSolution: leng != numLocalRows.\n",
3982 mypid_);
3983 exit(1);
3984 }
3985
3986 //-------------------------------------------------------------------
3987 // get the whole solution vector
3988 //-------------------------------------------------------------------
3989
3990 equations = new int[leng];
3991 if (localStartCol_ == -1)
3992 for ( i = 0; i < leng; i++ ) equations[i] = localStartRow_ + i - 1;
3993 else
3994 for ( i = 0; i < leng; i++ ) equations[i] = localStartCol_ + i;
3995
3996 HYPRE_IJVectorGetValues(HYx_,leng,equations,answers);
3997
3998 delete [] equations;
3999
4000 //-------------------------------------------------------------------
4001 // diagnostic message
4002 //-------------------------------------------------------------------
4003
4004 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
4005 printf("%4d : HYPRE_LSC::leaving getSolution.\n",mypid_);
4006 return (0);
4007 }
4008
4009 //***************************************************************************
4010 // used for getting the solution out of the solver, and into the application
4011 //---------------------------------------------------------------------------
4012
getSolnEntry(int eqnNumber,double & answer)4013 int HYPRE_LinSysCore::getSolnEntry(int eqnNumber, double& answer)
4014 {
4015 double val;
4016 int equation;
4017
4018 //-------------------------------------------------------------------
4019 // diagnostic message
4020 //-------------------------------------------------------------------
4021
4022 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
4023 printf("%4d : HYPRE_LSC::entering getSolnEntry.\n",mypid_);
4024
4025 //-------------------------------------------------------------------
4026 // get a single solution entry
4027 //-------------------------------------------------------------------
4028
4029 equation = eqnNumber; // incoming 0-based index
4030
4031 if (localStartCol_ == -1 && equation < localStartRow_-1 &&
4032 equation > localEndRow_ )
4033 {
4034 printf("%d : getSolnEntry ERROR - index out of range = %d.\n", mypid_,
4035 eqnNumber);
4036 exit(1);
4037 }
4038
4039 HYPRE_IJVectorGetValues(HYx_,1,&equation,&val);
4040 answer = val;
4041
4042 //-------------------------------------------------------------------
4043 // diagnostic message
4044 //-------------------------------------------------------------------
4045
4046 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 5 )
4047 printf("%4d : HYPRE_LSC::leaving getSolnEntry.\n",mypid_);
4048 return (0);
4049 }
4050
4051 //***************************************************************************
4052 // select which Krylov solver to use
4053 //---------------------------------------------------------------------------
4054
selectSolver(char * name)4055 void HYPRE_LinSysCore::selectSolver(char* name)
4056 {
4057 //-------------------------------------------------------------------
4058 // diagnostic message
4059 //-------------------------------------------------------------------
4060
4061 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
4062 {
4063 printf("%4d : HYPRE_LSC::entering selectSolver.\n",mypid_);
4064 printf("%4d : HYPRE_LSC::solver name = %s.\n",mypid_,name);
4065 }
4066
4067 //-------------------------------------------------------------------
4068 // if already been allocated, destroy it first
4069 //-------------------------------------------------------------------
4070
4071 if ( HYSolver_ != NULL )
4072 {
4073 if ( HYSolverID_ == HYPCG ) HYPRE_ParCSRPCGDestroy(HYSolver_);
4074 if ( HYSolverID_ == HYLSICG ) HYPRE_ParCSRLSICGDestroy(HYSolver_);
4075 if ( HYSolverID_ == HYHYBRID ) HYPRE_ParCSRHybridDestroy(HYSolver_);
4076 if ( HYSolverID_ == HYGMRES ) HYPRE_ParCSRGMRESDestroy(HYSolver_);
4077 if ( HYSolverID_ == HYFGMRES) HYPRE_ParCSRFGMRESDestroy(HYSolver_);
4078 if ( HYSolverID_ == HYCGSTAB) HYPRE_ParCSRBiCGSTABDestroy(HYSolver_);
4079 if ( HYSolverID_ == HYCGSTABL) HYPRE_ParCSRBiCGSTABLDestroy(HYSolver_);
4080 if ( HYSolverID_ == HYAMG) HYPRE_BoomerAMGDestroy(HYSolver_);
4081 if ( HYSolverID_ == HYTFQMR) HYPRE_ParCSRTFQmrDestroy(HYSolver_);
4082 if ( HYSolverID_ == HYBICGS) HYPRE_ParCSRBiCGSDestroy(HYSolver_);
4083 if ( HYSolverID_ == HYSYMQMR) HYPRE_ParCSRSymQMRDestroy(HYSolver_);
4084 }
4085
4086 //-------------------------------------------------------------------
4087 // check for the validity of the solver name
4088 //-------------------------------------------------------------------
4089
4090 if ( !strcmp(name, "cg" ) )
4091 {
4092 hypre_strcpy( HYSolverName_, name );
4093 HYSolverID_ = HYPCG;
4094 }
4095 else if ( !strcmp(name, "lsicg" ) )
4096 {
4097 hypre_strcpy( HYSolverName_, name );
4098 HYSolverID_ = HYLSICG;
4099 }
4100 else if ( !strcmp(name, "hybrid") )
4101 {
4102 hypre_strcpy( HYSolverName_, name );
4103 HYSolverID_ = HYHYBRID;
4104 }
4105 else if ( !strcmp(name, "gmres") )
4106 {
4107 hypre_strcpy( HYSolverName_, name );
4108 HYSolverID_ = HYGMRES;
4109 }
4110 else if ( !strcmp(name, "fgmres") )
4111 {
4112 hypre_strcpy( HYSolverName_, name );
4113 HYSolverID_ = HYFGMRES;
4114 }
4115 else if ( !strcmp(name, "bicgstab") )
4116 {
4117 hypre_strcpy( HYSolverName_, name );
4118 HYSolverID_ = HYCGSTAB;
4119 }
4120 else if ( !strcmp(name, "bicgstabl") )
4121 {
4122 hypre_strcpy( HYSolverName_, name );
4123 HYSolverID_ = HYCGSTABL;
4124 }
4125 else if ( !strcmp(name, "tfqmr") )
4126 {
4127 hypre_strcpy( HYSolverName_, name );
4128 HYSolverID_ = HYTFQMR;
4129 }
4130 else if ( !strcmp(name, "bicgs") )
4131 {
4132 hypre_strcpy( HYSolverName_, name );
4133 HYSolverID_ = HYBICGS;
4134 }
4135 else if ( !strcmp(name, "symqmr") )
4136 {
4137 hypre_strcpy( HYSolverName_, name );
4138 HYSolverID_ = HYSYMQMR;
4139 }
4140 else if ( !strcmp(name, "boomeramg") )
4141 {
4142 hypre_strcpy( HYSolverName_, name );
4143 HYSolverID_ = HYAMG;
4144 }
4145 else if ( !strcmp(name, "superlu") )
4146 {
4147 hypre_strcpy( HYSolverName_, name );
4148 HYSolverID_ = HYSUPERLU;
4149 }
4150 else if ( !strcmp(name, "superlux") )
4151 {
4152 hypre_strcpy( HYSolverName_, name );
4153 HYSolverID_ = HYSUPERLUX;
4154 }
4155 else if ( !strcmp(name, "dsuperlu") )
4156 {
4157 hypre_strcpy( HYSolverName_, name );
4158 #ifdef HYPRE_USING_DSUPERLU
4159 HYSolverID_ = HYDSUPERLU;
4160 #else
4161 printf("HYPRE_LinSysCore:: DSuperLU not available.\n");
4162 printf(" default solver to be GMRES.\n");
4163 HYSolverID_ = HYGMRES;
4164 #endif
4165 }
4166 else if ( !strcmp(name, "y12m") )
4167 {
4168 hypre_strcpy( HYSolverName_, name );
4169 HYSolverID_ = HYY12M;
4170 }
4171 else if ( !strcmp(name, "amge") )
4172 {
4173 hypre_strcpy( HYSolverName_, name );
4174 HYSolverID_ = HYAMGE;
4175 }
4176 else
4177 {
4178 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
4179 printf("HYPRE_LSC selectSolver : use default = gmres.\n");
4180 strcpy( HYSolverName_, "gmres" );
4181 HYSolverID_ = HYGMRES;
4182 }
4183
4184 //-------------------------------------------------------------------
4185 // instantiate solver
4186 //-------------------------------------------------------------------
4187
4188 switch ( HYSolverID_ )
4189 {
4190 case HYPCG :
4191 HYPRE_ParCSRPCGCreate(comm_, &HYSolver_);
4192 break;
4193 case HYLSICG :
4194 HYPRE_ParCSRLSICGCreate(comm_, &HYSolver_);
4195 break;
4196 case HYHYBRID :
4197 HYPRE_ParCSRHybridCreate(&HYSolver_);
4198 break;
4199 case HYGMRES :
4200 HYPRE_ParCSRGMRESCreate(comm_, &HYSolver_);
4201 break;
4202 case HYFGMRES :
4203 HYPRE_ParCSRFGMRESCreate(comm_, &HYSolver_);
4204 break;
4205 case HYCGSTAB :
4206 HYPRE_ParCSRBiCGSTABCreate(comm_, &HYSolver_);
4207 break;
4208 case HYCGSTABL :
4209 HYPRE_ParCSRBiCGSTABLCreate(comm_, &HYSolver_);
4210 break;
4211 case HYTFQMR :
4212 HYPRE_ParCSRTFQmrCreate(comm_, &HYSolver_);
4213 break;
4214 case HYBICGS :
4215 HYPRE_ParCSRBiCGSCreate(comm_, &HYSolver_);
4216 break;
4217 case HYSYMQMR :
4218 HYPRE_ParCSRSymQMRCreate(comm_, &HYSolver_);
4219 break;
4220 case HYAMG :
4221 HYPRE_BoomerAMGCreate( &HYSolver_);
4222 HYPRE_BoomerAMGSetCycleType(HYSolver_, 1);
4223 HYPRE_BoomerAMGSetMaxLevels(HYSolver_, 25);
4224 break;
4225 default:
4226 break;
4227 }
4228
4229 //-------------------------------------------------------------------
4230 // diagnostic message
4231 //-------------------------------------------------------------------
4232
4233 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
4234 printf("%4d : HYPRE_LSC::leaving selectSolver.\n",mypid_);
4235 return;
4236 }
4237
4238 //***************************************************************************
4239 // select which preconditioner to use
4240 //---------------------------------------------------------------------------
4241
selectPreconditioner(char * name)4242 void HYPRE_LinSysCore::selectPreconditioner(char *name)
4243 {
4244 //-------------------------------------------------------------------
4245 // diagnostic message
4246 //-------------------------------------------------------------------
4247
4248 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3)
4249 printf("%4d : HYPRE_LSC::entering selectPreconditioner = %s.\n",
4250 mypid_, name);
4251
4252 //-------------------------------------------------------------------
4253 // if already been allocated, destroy it first
4254 //-------------------------------------------------------------------
4255
4256 HYPreconSetup_ = 0;
4257 parasailsReuse_ = 0;
4258 if ( HYPrecon_ != NULL )
4259 {
4260 if (HYPreconID_ == HYPILUT)
4261 HYPRE_ParCSRPilutDestroy(HYPrecon_);
4262
4263 else if (HYPreconID_ == HYPARASAILS)
4264 HYPRE_ParCSRParaSailsDestroy(HYPrecon_);
4265
4266 else if (HYPreconID_ == HYBOOMERAMG)
4267 HYPRE_BoomerAMGDestroy(HYPrecon_);
4268
4269 else if (HYPreconID_ == HYDDILUT)
4270 HYPRE_LSI_DDIlutDestroy(HYPrecon_);
4271
4272 else if (HYPreconID_ == HYSCHWARZ)
4273 HYPRE_LSI_SchwarzDestroy(HYPrecon_);
4274
4275 else if (HYPreconID_ == HYDDICT)
4276 HYPRE_LSI_DDICTDestroy(HYPrecon_);
4277
4278 else if (HYPreconID_ == HYPOLY)
4279 HYPRE_LSI_PolyDestroy(HYPrecon_);
4280
4281 else if (HYPreconID_ == HYEUCLID)
4282 HYPRE_EuclidDestroy(HYPrecon_);
4283
4284 else if (HYPreconID_ == HYBLOCK)
4285 HYPRE_LSI_BlockPrecondDestroy(HYPrecon_);
4286
4287 #ifdef HAVE_ML
4288 else if (HYPreconID_ == HYML)
4289 HYPRE_LSI_MLDestroy(HYPrecon_);
4290 #endif
4291 #ifdef HAVE_MLI
4292 else if (HYPreconID_ == HYMLI)
4293 HYPRE_LSI_MLIDestroy(HYPrecon_);
4294 #endif
4295 else if (HYPreconID_ == HYUZAWA)
4296 HYPRE_LSI_UzawaDestroy(HYPrecon_);
4297 #ifdef HAVE_SYSPDE
4298 else if (HYPreconID_ == HYSYSPDE)
4299 HYPRE_ParCSRSysPDEDestroy(HYPrecon_);
4300 #endif
4301 #ifdef HYPRE_USING_DSUPERLU
4302 else if (HYPreconID_ == HYDSLU)
4303 HYPRE_LSI_DSuperLUDestroy(HYPrecon_);
4304 #endif
4305 }
4306
4307 //-------------------------------------------------------------------
4308 // check for the validity of the preconditioner name
4309 //-------------------------------------------------------------------
4310
4311 if (!strcmp(name, "identity"))
4312 {
4313 hypre_strcpy(HYPreconName_, name);
4314 HYPreconID_ = HYIDENTITY;
4315 }
4316 else if (!strcmp(name, "diagonal"))
4317 {
4318 hypre_strcpy(HYPreconName_, name);
4319 HYPreconID_ = HYDIAGONAL;
4320 }
4321 else if (!strcmp(name, "pilut"))
4322 {
4323 hypre_strcpy(HYPreconName_, name);
4324 HYPreconID_ = HYPILUT;
4325 }
4326 else if (!strcmp(name, "parasails"))
4327 {
4328 hypre_strcpy(HYPreconName_, name);
4329 HYPreconID_ = HYPARASAILS;
4330 }
4331 else if (!strcmp(name, "boomeramg"))
4332 {
4333 hypre_strcpy(HYPreconName_, name);
4334 HYPreconID_ = HYBOOMERAMG;
4335 }
4336 else if (!strcmp(name, "ddilut"))
4337 {
4338 hypre_strcpy(HYPreconName_, name);
4339 HYPreconID_ = HYDDILUT;
4340 }
4341 else if (!strcmp(name, "schwarz"))
4342 {
4343 hypre_strcpy(HYPreconName_, name);
4344 HYPreconID_ = HYSCHWARZ;
4345 }
4346 else if (!strcmp(name, "ddict"))
4347 {
4348 hypre_strcpy(HYPreconName_, name);
4349 HYPreconID_ = HYDDICT;
4350 }
4351 else if (!strcmp(name, "poly"))
4352 {
4353 hypre_strcpy(HYPreconName_, name);
4354 HYPreconID_ = HYPOLY;
4355 }
4356 else if (!strcmp(name, "euclid"))
4357 {
4358 hypre_strcpy(HYPreconName_, name);
4359 HYPreconID_ = HYEUCLID;
4360 }
4361 else if (!strcmp(name, "blockP"))
4362 {
4363 hypre_strcpy(HYPreconName_, name);
4364 HYPreconID_ = HYBLOCK;
4365 }
4366 else if (!strcmp(name, "ml"))
4367 {
4368 #ifdef HAVE_ML
4369 hypre_strcpy(HYPreconName_, name);
4370 HYPreconID_ = HYML;
4371 #else
4372 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3)
4373 {
4374 printf("selectPreconditioner - ML not available.\n");
4375 printf(" set default to diagonal.\n");
4376 }
4377 strcpy(HYPreconName_, "diagonal");
4378 HYPreconID_ = HYDIAGONAL;
4379 #endif
4380 }
4381 else if (!strcmp(name, "mlmaxwell"))
4382 {
4383 #ifdef HAVE_MLMAXWELL
4384 hypre_strcpy(HYPreconName_, name);
4385 HYPreconID_ = HYMLMAXWELL;
4386 #else
4387 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3)
4388 {
4389 printf("selectPreconditioner - MLMaxwell not available.\n");
4390 printf(" set default to diagonal.\n");
4391 }
4392 strcpy(HYPreconName_, "diagonal");
4393 HYPreconID_ = HYDIAGONAL;
4394 #endif
4395 }
4396 else if (!strcmp(name, "mli"))
4397 {
4398 #ifdef HAVE_MLI
4399 hypre_strcpy(HYPreconName_, name);
4400 HYPreconID_ = HYMLI;
4401 #else
4402 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3)
4403 {
4404 printf("selectPreconditioner - MLI not available.\n");
4405 printf(" set default to diagonal.\n");
4406 }
4407 strcpy(HYPreconName_, "diagonal");
4408 HYPreconID_ = HYDIAGONAL;
4409 #endif
4410 }
4411 else if (!strcmp(name, "ams"))
4412 {
4413 hypre_strcpy(HYPreconName_, name);
4414 HYPreconID_ = HYAMS;
4415 }
4416 else if (!strcmp(name, "uzawa"))
4417 {
4418 hypre_strcpy(HYPreconName_, name);
4419 HYPreconID_ = HYUZAWA;
4420 }
4421 #ifdef HAVE_SYSPDE
4422 else if (!strcmp(name, "syspde"))
4423 {
4424 hypre_strcpy(HYPreconName_, name);
4425 HYPreconID_ = HYSYSPDE;
4426 }
4427 #endif
4428 #ifdef HYPRE_USING_DSUPERLU
4429 else if (!strcmp(name, "dsuperlu"))
4430 {
4431 hypre_strcpy(HYPreconName_, name);
4432 HYPreconID_ = HYDSLU;
4433 }
4434 #endif
4435 else
4436 {
4437 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3)
4438 {
4439 printf("selectPreconditioner error : invalid option.\n");
4440 printf(" use default = diagonal.\n");
4441 }
4442 strcpy(HYPreconName_, "diagonal");
4443 HYPreconID_ = HYDIAGONAL;
4444 }
4445
4446 //-------------------------------------------------------------------
4447 // instantiate preconditioner
4448 //-------------------------------------------------------------------
4449
4450 switch (HYPreconID_)
4451 {
4452 case HYIDENTITY :
4453 HYPrecon_ = NULL;
4454 break;
4455
4456 case HYDIAGONAL :
4457 HYPrecon_ = NULL;
4458 break;
4459
4460 case HYPILUT :
4461 HYPRE_ParCSRPilutCreate(comm_, &HYPrecon_);
4462 HYPRE_ParCSRPilutSetMaxIter(HYPrecon_, 1);
4463 break;
4464
4465 case HYPARASAILS :
4466 HYPRE_ParCSRParaSailsCreate(comm_, &HYPrecon_);
4467 break;
4468
4469 case HYBOOMERAMG :
4470 HYPRE_BoomerAMGCreate(&HYPrecon_);
4471 HYPRE_BoomerAMGSetMaxIter(HYPrecon_, 1);
4472 HYPRE_BoomerAMGSetCycleType(HYPrecon_, 1);
4473 HYPRE_BoomerAMGSetMaxLevels(HYPrecon_, 25);
4474 HYPRE_BoomerAMGSetMeasureType(HYPrecon_, 0);
4475 break;
4476
4477 case HYDDILUT :
4478 HYPRE_LSI_DDIlutCreate(comm_, &HYPrecon_);
4479 break;
4480
4481 case HYSCHWARZ :
4482 HYPRE_LSI_SchwarzCreate(comm_, &HYPrecon_);
4483 break;
4484
4485 case HYDDICT :
4486 HYPRE_LSI_DDICTCreate(comm_, &HYPrecon_);
4487 break;
4488
4489 case HYPOLY :
4490 HYPRE_LSI_PolyCreate(comm_, &HYPrecon_);
4491 break;
4492
4493 case HYEUCLID :
4494 HYPRE_EuclidCreate(comm_, &HYPrecon_);
4495 break;
4496
4497 case HYBLOCK :
4498 HYPRE_LSI_BlockPrecondCreate(comm_, &HYPrecon_);
4499 break;
4500
4501 case HYML :
4502 #ifdef HAVE_ML
4503 HYPRE_LSI_MLCreate(comm_, &HYPrecon_);
4504 #else
4505 printf("HYPRE_LSC::selectPreconditioner - ML not supported.\n");
4506 #endif
4507 break;
4508 case HYMLI :
4509 #ifdef HAVE_MLI
4510 HYPRE_LSI_MLICreate(comm_, &HYPrecon_);
4511 #else
4512 printf("HYPRE_LSC::selectPreconditioner - MLI not supported.\n");
4513 #endif
4514 break;
4515 case HYMLMAXWELL :
4516 #ifdef HAVE_MLMAXWELL
4517 HYPRE_LSI_MLMaxwellCreate(comm_, &HYPrecon_);
4518 #else
4519 printf("HYPRE_LSC::selectPreconditioner-MLMaxwell unsupported.\n");
4520 #endif
4521 break;
4522 case HYAMS :
4523 HYPRE_AMSCreate(&HYPrecon_);
4524 break;
4525 case HYUZAWA :
4526 HYPRE_LSI_UzawaCreate(comm_, &HYPrecon_);
4527 break;
4528 case HYSYSPDE :
4529 #ifdef HAVE_SYSPDE
4530 HYPRE_ParCSRSysPDECreate(comm_, sysPDENVars_, &HYPrecon_);
4531 #else
4532 printf("HYPRE_LSC::selectPreconditioner-SYSPDE unsupported.\n");
4533 #endif
4534 break;
4535 case HYDSLU :
4536 #ifdef HYPRE_USING_DSUPERLU
4537 HYPRE_LSI_DSuperLUCreate(comm_, &HYPrecon_);
4538 #else
4539 printf("HYPRE_LSC::selectPreconditioner-DSUPERLU unsupported.\n");
4540 #endif
4541 break;
4542 }
4543
4544 //-------------------------------------------------------------------
4545 // diagnostic message
4546 //-------------------------------------------------------------------
4547
4548 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3)
4549 printf("%4d : HYPRE_LSC::leaving selectPreconditioner.\n",mypid_);
4550 }
4551
4552 //***************************************************************************
4553 // form the residual vector
4554 //---------------------------------------------------------------------------
4555
formResidual(double * values,int leng)4556 int HYPRE_LinSysCore::formResidual(double* values, int leng)
4557 {
4558 int i, index, nrows;
4559 HYPRE_ParCSRMatrix A_csr;
4560 HYPRE_ParVector x_csr;
4561 HYPRE_ParVector b_csr;
4562 HYPRE_ParVector r_csr;
4563
4564 //-------------------------------------------------------------------
4565 // diagnostic message and error checking
4566 //-------------------------------------------------------------------
4567
4568 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
4569 printf("%4d : HYPRE_LSC::entering formResidual.\n", mypid_);
4570
4571 nrows = localEndRow_ - localStartRow_ + 1;
4572 if (leng != nrows)
4573 {
4574 printf("%4d : HYPRE_LSC::formResidual ERROR - inleng != numLocalRows",
4575 mypid_);
4576 printf(" numLocalRows, inleng = %d %d",nrows,leng);
4577 return (0);
4578 }
4579 if ( ! systemAssembled_ )
4580 {
4581 printf("%4d : HYPRE_LSC formResidual ERROR : system not assembled.\n",
4582 mypid_);
4583 exit(1);
4584 }
4585
4586 //-------------------------------------------------------------------
4587 // fetch matrix and vector pointers
4588 //-------------------------------------------------------------------
4589
4590 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
4591 HYPRE_IJVectorGetObject(HYx_, (void **) &x_csr);
4592 HYPRE_IJVectorGetObject(HYb_, (void **) &b_csr);
4593 HYPRE_IJVectorGetObject(HYr_, (void **) &r_csr);
4594
4595 //-------------------------------------------------------------------
4596 // form residual vector
4597 //-------------------------------------------------------------------
4598
4599 HYPRE_ParVectorCopy( b_csr, r_csr );
4600 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
4601
4602 //-------------------------------------------------------------------
4603 // fetch residual vector
4604 //-------------------------------------------------------------------
4605
4606 for ( i = localStartRow_-1; i < localEndRow_; i++ )
4607 {
4608 index = i - localStartRow_ + 1;
4609 HYPRE_IJVectorGetValues(HYr_, 1, &i, &values[index]);
4610 }
4611
4612 //-------------------------------------------------------------------
4613 // diagnostic message
4614 //-------------------------------------------------------------------
4615
4616 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
4617 printf("%4d : HYPRE_LSC::leaving formResidual.\n", mypid_);
4618 return (0);
4619 }
4620
4621 //***************************************************************************
4622 // solve the linear system
4623 //---------------------------------------------------------------------------
4624
launchSolver(int & solveStatus,int & iterations)4625 int HYPRE_LinSysCore::launchSolver(int& solveStatus, int &iterations)
4626 {
4627 int i, j, numIterations=0, status, localNRows;
4628 int startRow, *procNRows, rowSize, *colInd, nnz, nrows;
4629 int slideCheck[2];
4630 #ifdef HAVE_MLI
4631 int *constrMap, *constrEqns, ncount, *iArray;
4632 double *tempNodalCoord;
4633 #endif
4634 int *numSweeps, *relaxType, reduceAFlag;
4635 int *matSizes, *rowInd, retFlag, tempIter, nTrials;
4636 double rnorm=0.0, ddata, *colVal, *relaxWt, *diagVals;
4637 double stime, etime, ptime, rtime1, rtime2, newnorm;
4638 double rnorm0, rnorm1, convRate, rateThresh;
4639 char fname[40], paramString[100];
4640 FILE *fp;
4641 HYPRE_IJMatrix TempA, IJI;
4642 HYPRE_IJVector TempX, TempB, TempR;
4643 #ifdef HAVE_MLI
4644 HYPRE_ParCSRMatrix perturb_csr;
4645 #endif
4646 HYPRE_ParCSRMatrix A_csr, I_csr, normalA_csr;
4647 HYPRE_ParVector x_csr, b_csr, r_csr;
4648 HYPRE_SlideReduction *slideObj = (HYPRE_SlideReduction *) slideObj_;
4649
4650 //-------------------------------------------------------------------
4651 // Clear Errors
4652 //-------------------------------------------------------------------
4653 HYPRE_ClearAllErrors();
4654
4655 //-------------------------------------------------------------------
4656 // diagnostic message
4657 //-------------------------------------------------------------------
4658
4659 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
4660 printf("%4d : HYPRE_LSC::entering launchSolver.\n", mypid_);
4661
4662 //-------------------------------------------------------------------
4663 // see if Schur or slide reduction is to be performed
4664 //-------------------------------------------------------------------
4665
4666 rnorm_ = 0.0;
4667 MPI_Barrier(comm_);
4668 rtime1 = LSC_Wtime();
4669 if ( schurReduction_ == 1 && schurReductionCreated_ == 0 )
4670 {
4671 buildSchurReducedSystem();
4672 schurReductionCreated_ = 1;
4673 }
4674 else if ( schurReduction_ == 1 ) buildSchurReducedRHS();
4675
4676 if ( schurReduction_ == 0 && slideReduction_ != 0 )
4677 {
4678 if ( constrList_ != NULL ) delete [] constrList_;
4679 constrList_ = NULL;
4680 if ( slideReduction_ == 1 ) buildSlideReducedSystem();
4681 else if ( slideReduction_ == 2 ) buildSlideReducedSystem2();
4682 else if ( slideReduction_ == 3 || slideReduction_ == 4 )
4683 {
4684 if (slideObj == NULL)
4685 {
4686 slideObj = new HYPRE_SlideReduction(comm_);
4687 slideObj_ = (void *) slideObj;
4688 }
4689 TempA = currA_;
4690 TempX = currX_;
4691 TempB = currB_;
4692 TempR = currR_;
4693 HYPRE_IJVectorGetLocalRange(HYb_,&slideCheck[0],&slideCheck[1]);
4694 // check to see if it has been reduced before
4695 // if so, need to redo B and X
4696 reduceAFlag = 1;
4697 if (currA_ != HYA_)
4698 {
4699 HYPRE_IJVectorDestroy(currB_);
4700 HYPRE_IJVectorDestroy(currX_);
4701 HYPRE_IJVectorDestroy(currR_);
4702 currB_ = HYb_;
4703 currX_ = HYx_;
4704 currR_ = HYr_;
4705 reduceAFlag = 0;
4706 }
4707
4708 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
4709 slideObj->setOutputLevel(1);
4710 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE2 )
4711 slideObj->setOutputLevel(2);
4712 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
4713 slideObj->setOutputLevel(3);
4714 if ( slideReductionMinNorm_ >= 0.0 )
4715 slideObj->setBlockMinNorm( slideReductionMinNorm_ );
4716 if ( slideReductionScaleMatrix_ == 1 )
4717 slideObj->setScaleMatrix();
4718 slideObj->setTruncationThreshold( truncThresh_ );
4719 if ( slideReduction_ == 4 ) slideObj->setUseSimpleScheme();
4720 slideObj->setup(currA_, currX_, currB_);
4721 if ( slideReductionScaleMatrix_ == 1 && HYPreconID_ == HYMLI )
4722 {
4723 diagVals = slideObj->getMatrixDiagonal();
4724 nrows = slideObj->getMatrixNumRows();
4725 HYPRE_LSI_MLILoadMatrixScalings(HYPrecon_, nrows, diagVals);
4726 }
4727 #ifdef HAVE_MLI
4728 if ( HYPreconID_ == HYMLI )
4729 {
4730 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
4731 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
4732 slideObj->getProcConstraintMap(&constrMap);
4733 HYPRE_LSI_MLIAdjustNodeEqnMap(HYPrecon_, procNRows, constrMap);
4734 j = constrMap[mypid_+1] - constrMap[mypid_];
4735 free(procNRows);
4736 slideObj->getSlaveEqnList(&constrEqns);
4737 slideObj->getPerturbationMatrix(&perturb_csr);
4738 HYPRE_LSI_MLIAdjustNullSpace(HYPrecon_,j,constrEqns,perturb_csr);
4739 }
4740 #endif
4741 if (reduceAFlag == 1)
4742 {
4743 slideObj->getReducedMatrix(&currA_);
4744 slideObj->getReducedAuxVector(&currR_);
4745 }
4746 slideObj->getReducedSolnVector(&currX_);
4747 slideObj->getReducedRHSVector(&currB_);
4748 if ( currA_ == NULL )
4749 {
4750 currA_ = TempA;
4751 currX_ = TempX;
4752 currB_ = TempB;
4753 currR_ = TempR;
4754 }
4755 }
4756 }
4757
4758 MPI_Barrier(comm_);
4759 rtime2 = LSC_Wtime();
4760
4761 //-------------------------------------------------------------------
4762 // if normal equation requested
4763 //-------------------------------------------------------------------
4764
4765 if ( (normalEqnFlag_ & 1) != 0 )
4766 {
4767 if ( (normalEqnFlag_ & 2) == 0 )
4768 {
4769 if ( HYnormalA_ != NULL ) HYPRE_IJMatrixDestroy(HYnormalA_);
4770 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
4771 HYPRE_IJMatrixCreate(comm_, localStartRow_-1,
4772 localEndRow_-1, localStartRow_-1, localEndRow_-1,&IJI);
4773 HYPRE_IJMatrixSetObjectType(IJI, HYPRE_PARCSR);
4774 localNRows = localEndRow_ - localStartRow_ + 1;
4775 matSizes = new int[localNRows];
4776 rowInd = new int[localNRows];
4777 colInd = new int[localNRows];
4778 colVal = new double[localNRows];
4779 for ( i = 0; i < localNRows; i++ )
4780 {
4781 matSizes[i] = 1;
4782 rowInd[i] = localStartRow_ - 1 + i;
4783 colInd[i] = rowInd[i];
4784 colVal[i] = 1.0;
4785 }
4786 HYPRE_IJMatrixSetRowSizes(IJI, matSizes);
4787 HYPRE_IJMatrixInitialize(IJI);
4788 HYPRE_IJMatrixSetValues(IJI, localNRows, matSizes,
4789 (const int *) rowInd, (const int *) colInd,
4790 (const double *) colVal);
4791 delete [] rowInd;
4792 delete [] colInd;
4793 delete [] colVal;
4794 HYPRE_IJMatrixAssemble(IJI);
4795 HYPRE_IJMatrixGetObject(IJI, (void **) &I_csr);
4796 hypre_BoomerAMGBuildCoarseOperator((hypre_ParCSRMatrix*) A_csr,
4797 (hypre_ParCSRMatrix*) I_csr, (hypre_ParCSRMatrix*) A_csr,
4798 (hypre_ParCSRMatrix**) &normalA_csr);
4799 HYPRE_IJMatrixDestroy( IJI );
4800 HYPRE_IJMatrixCreate(comm_, localStartRow_-1,
4801 localEndRow_-1, localStartRow_-1, localEndRow_-1,&HYnormalA_);
4802 HYPRE_IJMatrixSetObjectType(HYnormalA_, HYPRE_PARCSR);
4803 for ( i = localStartRow_-1; i < localEndRow_; i++ )
4804 {
4805 HYPRE_ParCSRMatrixGetRow(normalA_csr,i,&rowSize,NULL,NULL);
4806 matSizes[i-localStartRow_+1] = rowSize;
4807 HYPRE_ParCSRMatrixRestoreRow(normalA_csr,i,&rowSize,NULL,NULL);
4808 }
4809 HYPRE_IJMatrixSetRowSizes(HYnormalA_, matSizes);
4810 HYPRE_IJMatrixInitialize(HYnormalA_);
4811 for ( i = localStartRow_-1; i < localEndRow_; i++ )
4812 {
4813 HYPRE_ParCSRMatrixGetRow(normalA_csr,i,&rowSize,&colInd,&colVal);
4814 HYPRE_IJMatrixSetValues(HYnormalA_, 1, &rowSize,
4815 (const int *) &i, (const int *) colInd,
4816 (const double *) colVal);
4817 HYPRE_ParCSRMatrixRestoreRow(normalA_csr,i,&rowSize,&colInd,&colVal);
4818 }
4819 HYPRE_IJMatrixAssemble(HYnormalA_);
4820 delete [] matSizes;
4821 normalEqnFlag_ |= 2;
4822 }
4823 if ( (normalEqnFlag_ & 4) == 0 )
4824 {
4825 if ( HYnormalB_ != NULL ) HYPRE_IJVectorDestroy(HYnormalB_);
4826 HYPRE_IJVectorCreate(comm_, localStartRow_-1, localEndRow_-1,
4827 &HYnormalB_);
4828 HYPRE_IJVectorSetObjectType(HYnormalB_, HYPRE_PARCSR);
4829 HYPRE_IJVectorInitialize(HYnormalB_);
4830 HYPRE_IJVectorAssemble(HYnormalB_);
4831 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
4832 HYPRE_IJVectorGetObject(HYb_, (void **) &b_csr);
4833 HYPRE_IJVectorGetObject(HYnormalB_, (void **) &r_csr);
4834 HYPRE_ParCSRMatrixMatvecT( 1.0, A_csr, b_csr, 0.0, r_csr );
4835 normalEqnFlag_ |= 4;
4836 }
4837 }
4838
4839 //-------------------------------------------------------------------
4840 // fetch matrix and vector pointers
4841 //-------------------------------------------------------------------
4842
4843 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
4844 HYPRE_IJVectorGetObject(currX_, (void **) &x_csr);
4845 HYPRE_IJVectorGetObject(currB_, (void **) &b_csr);
4846 HYPRE_IJVectorGetObject(currR_, (void **) &r_csr);
4847 if ( A_csr == NULL || x_csr == NULL || b_csr == NULL || r_csr == NULL )
4848 {
4849 printf("%4d : HYPRE_LSC::launchSolver ERROR.\n",mypid_);
4850 printf(" csr pointers null \n");
4851 printf(" Did you forget to call matrixLoadComplete?\n");
4852 exit(1);
4853 }
4854 if ( (normalEqnFlag_ & 7) == 7 )
4855 {
4856 HYPRE_IJMatrixGetObject(HYnormalA_, (void **) &A_csr);
4857 HYPRE_IJVectorGetObject(HYnormalB_, (void **) &b_csr);
4858 }
4859
4860 //-------------------------------------------------------------------
4861 // diagnostics (print the reduced matrix to a file)
4862 //-------------------------------------------------------------------
4863
4864 if ( HYOutputLevel_ & HYFEI_PRINTREDMAT )
4865 {
4866 if ( HYOutputLevel_ & HYFEI_PRINTPARCSRMAT )
4867 {
4868 printf("%4d : HYPRE_LSC::print matrix/rhs to files(A)\n",mypid_);
4869 sprintf(fname, "HYPRE_Mat");
4870 HYPRE_ParCSRMatrixPrint( A_csr, fname);
4871 sprintf(fname, "HYPRE_RHS");
4872 HYPRE_ParVectorPrint( b_csr, fname);
4873 }
4874 else
4875 {
4876 printf("%4d : HYPRE_LSC::print matrix/rhs to files(B)\n",mypid_);
4877 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
4878 startRow = procNRows[mypid_];
4879 nrows = procNRows[mypid_+1] - startRow;
4880 free( procNRows );
4881
4882 sprintf(fname, "hypre_mat.out.%d", mypid_);
4883 fp = fopen( fname, "w");
4884 nnz = 0;
4885 for ( i = startRow; i < startRow+nrows; i++ )
4886 {
4887 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
4888 for ( j = 0; j < rowSize; j++ ) if ( colVal[j] != 0.0 ) nnz++;
4889 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
4890 }
4891 fprintf(fp, "%6d %7d \n", nrows, nnz);
4892 for ( i = startRow; i < startRow+nrows; i++ )
4893 {
4894 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
4895 for ( j = 0; j < rowSize; j++ )
4896 {
4897 if ( colVal[j] != 0.0 )
4898 fprintf(fp, "%6d %6d %25.8e\n",i+1,colInd[j]+1,colVal[j]);
4899 }
4900 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
4901 }
4902 fclose(fp);
4903
4904 sprintf(fname, "hypre_rhs.out.%d", mypid_);
4905 fp = fopen( fname, "w");
4906 fprintf(fp, "%6d \n", nrows);
4907 for ( i = startRow; i < startRow+nrows; i++ )
4908 {
4909 HYPRE_IJVectorGetValues(currB_, 1, &i, &ddata);
4910 fprintf(fp, "%6d %25.8e \n", i+1, ddata);
4911 }
4912 fclose(fp);
4913 MPI_Barrier(comm_);
4914 }
4915 if ( MLI_NumNodes_ > 0 )
4916 {
4917 fp = fopen("rbm","w");
4918 for (i = 0; i < MLI_NumNodes_; i++)
4919 for (j = 0; j < MLI_FieldSize_; j++)
4920 fprintf(fp,"%8d %25.16e\n", MLI_EqnNumbers_[i]+j+1,
4921 MLI_NodalCoord_[i*3+j]);
4922 fclose(fp);
4923 }
4924 if ( HYOutputLevel_ & HYFEI_STOPAFTERPRINT ) exit(1);
4925 }
4926
4927 #ifdef HAVE_AMGE
4928 if ( HYOutputLevel_ & HYFEI_PRINTFEINFO )
4929 {
4930 HYPRE_LSI_AMGeWriteToFile();
4931 }
4932 #endif
4933
4934 //-------------------------------------------------------------------
4935 // choose PCG, GMRES, ... or direct solver
4936 //-------------------------------------------------------------------
4937
4938 MPI_Barrier(comm_);
4939 status = 1;
4940 stime = LSC_Wtime();
4941 ptime = stime;
4942
4943 if ( projectionScheme_ == 1 )
4944 {
4945 computeAConjProjection(A_csr, x_csr, b_csr);
4946 }
4947 else if ( projectionScheme_ == 2 )
4948 {
4949 computeMinResProjection(A_csr, x_csr, b_csr);
4950 }
4951
4952 #ifdef HAVE_MLI
4953 if ( HYPreconID_ == HYMLI && feData_ != NULL )
4954 {
4955 if (haveFEData_ == 1) HYPRE_LSI_MLISetFEData( HYPrecon_, feData_ );
4956 if (haveFEData_ == 2) HYPRE_LSI_MLISetSFEI( HYPrecon_, feData_ );
4957 }
4958 if ( HYPreconID_ == HYMLI && MLI_EqnNumbers_ != NULL )
4959 {
4960 iArray = new int[MLI_NumNodes_];
4961 for (i = 0; i < MLI_NumNodes_; i++) iArray[i] = i;
4962 HYPRE_LSI_qsort1a(MLI_EqnNumbers_, iArray, 0, MLI_NumNodes_-1);
4963 tempNodalCoord = MLI_NodalCoord_;
4964 ncount = 1;
4965 for (i = 1; i < MLI_NumNodes_; i++)
4966 if (MLI_EqnNumbers_[i] != MLI_EqnNumbers_[ncount-1]) ncount++;
4967 MLI_NodalCoord_ = new double[ncount*MLI_FieldSize_];
4968 for (j = 0; j < MLI_FieldSize_; j++)
4969 MLI_NodalCoord_[j] = tempNodalCoord[iArray[0]*MLI_FieldSize_+j];
4970 ncount = 1;
4971 for (i = 1; i < MLI_NumNodes_; i++)
4972 {
4973 if (MLI_EqnNumbers_[i] != MLI_EqnNumbers_[ncount-1])
4974 {
4975 MLI_EqnNumbers_[ncount] = MLI_EqnNumbers_[i];
4976 for (j = 0; j < MLI_FieldSize_; j++)
4977 MLI_NodalCoord_[ncount*MLI_FieldSize_+j] =
4978 tempNodalCoord[iArray[i]*MLI_FieldSize_+j];
4979 ncount++;
4980 }
4981 }
4982 MLI_NumNodes_ = ncount;
4983 //hypre_assert((MLI_NumNodes_*MLI_FieldSize_)==(localEndRow_-localStartRow_+1));
4984 delete [] tempNodalCoord;
4985 delete [] iArray;
4986 for (i = 0; i < MLI_NumNodes_; i++)
4987 {
4988 if (MLI_NodalCoord_[i] == -99999.0)
4989 printf("%d : HYPRE launchSolver ERROR - coord %d not filled.\n",
4990 mypid_, i);
4991 }
4992 HYPRE_LSI_MLILoadNodalCoordinates(HYPrecon_, MLI_NumNodes_,
4993 MLI_FieldSize_, MLI_EqnNumbers_, MLI_FieldSize_,
4994 MLI_NodalCoord_, localEndRow_-localStartRow_+1);
4995 }
4996 #endif
4997 #if 0
4998 // replaced by better scheme, to be deleted later
4999 if ( HYPreconID_ == HYAMS && MLI_EqnNumbers_ != NULL )
5000 {
5001 HYPRE_LSI_BuildNodalCoordinates();
5002 }
5003 #endif
5004
5005 switch ( HYSolverID_ )
5006 {
5007 //----------------------------------------------------------------
5008 // choose PCG
5009 //----------------------------------------------------------------
5010
5011 case HYPCG :
5012 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5013 {
5014 printf("***************************************************\n");
5015 printf("* Preconditioned Conjugate Gradient solver \n");
5016 printf("* maximum no. of iterations = %d\n", maxIterations_);
5017 printf("* convergence tolerance = %e\n", tolerance_);
5018 printf("*--------------------------------------------------\n");
5019 }
5020 setupPCGPrecon();
5021 HYPRE_ParCSRPCGSetMaxIter(HYSolver_, maxIterations_);
5022 HYPRE_ParCSRPCGSetRelChange(HYSolver_, 0);
5023 HYPRE_ParCSRPCGSetTwoNorm(HYSolver_, 1);
5024 HYPRE_PCGSetRecomputeResidual(HYSolver_, pcgRecomputeRes_);
5025 if ( normAbsRel_ == 0 )
5026 {
5027 HYPRE_PCGSetTol(HYSolver_, tolerance_);
5028 }
5029 else
5030 {
5031 HYPRE_PCGSetTol(HYSolver_, 0.0);
5032 HYPRE_PCGSetAbsoluteTol(HYSolver_, tolerance_);
5033 }
5034 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5035 {
5036 if ( mypid_ == 0 )
5037 printf("***************************************************\n");
5038 HYPRE_ParCSRPCGSetPrintLevel(HYSolver_, 1);
5039 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 )
5040 HYPRE_ParCSRPCGSetPrintLevel(HYSolver_, 2);
5041 }
5042 retFlag = HYPRE_ParCSRPCGSetup(HYSolver_, A_csr, b_csr, x_csr);
5043 if ( retFlag != 0 )
5044 {
5045 printf("HYPRE_LSC::launchSolver ERROR : in PCG setup.\n");
5046 return retFlag;
5047 }
5048 // if variable mli preconditioner (SA and GSA)
5049 if ( MLI_Hybrid_GSA_ && HYPreconID_ == HYMLI )
5050 {
5051 if ( normAbsRel_ == 0 )
5052 {
5053 HYPRE_ParVectorCopy( b_csr, r_csr );
5054 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5055 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm0);
5056 rnorm0 = sqrt( rnorm0 );
5057 }
5058 else rnorm0 = 1.0;
5059 HYPRE_ParCSRPCGSetMaxIter(HYSolver_, MLI_Hybrid_MaxIter_);
5060 rateThresh = 1.0;
5061 for ( i = 0; i < MLI_Hybrid_MaxIter_; i++ )
5062 rateThresh *= MLI_Hybrid_ConvRate_;
5063 }
5064 MPI_Barrier( comm_ );
5065 ptime = LSC_Wtime();
5066 retFlag = HYPRE_ParCSRPCGSolve(HYSolver_, A_csr, b_csr, x_csr);
5067 HYPRE_ParCSRPCGGetNumIterations(HYSolver_, &numIterations);
5068 HYPRE_ParVectorCopy( b_csr, r_csr );
5069 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5070 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5071 rnorm = sqrt( rnorm );
5072 // if variable mli preconditioner (SA and GSA)
5073 if ( MLI_Hybrid_GSA_ && HYPreconID_ == HYMLI )
5074 {
5075 nTrials = 1;
5076 if (rnorm/rnorm0 >= tolerance_)
5077 {
5078 HYPRE_ParCSRPCGSolve(HYSolver_, A_csr, b_csr, x_csr);
5079 HYPRE_ParCSRPCGGetNumIterations(HYSolver_, &tempIter);
5080 numIterations += tempIter;
5081 HYPRE_ParVectorCopy( b_csr, r_csr );
5082 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5083 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm1);
5084 rnorm1 = sqrt( rnorm1 );
5085 convRate = rnorm1 / rnorm;
5086 rnorm = rnorm1;
5087 }
5088 while ((rnorm/rnorm0)>=tolerance_ && nTrials<MLI_Hybrid_NTrials_)
5089 {
5090 nTrials++;
5091 if ( convRate > rateThresh )
5092 {
5093 if ( MLI_Hybrid_NSIncr_ > 1 )
5094 sprintf(paramString, "MLI incrNullSpaceDim %d",
5095 MLI_Hybrid_NSIncr_);
5096 else
5097 sprintf(paramString, "MLI incrNullSpaceDim 2");
5098 HYPRE_LSI_MLISetParams(HYPrecon_, paramString);
5099 HYPRE_ParCSRPCGSetup(HYSolver_, A_csr, b_csr, x_csr);
5100 }
5101 HYPRE_ParCSRPCGSolve(HYSolver_, A_csr, b_csr, x_csr);
5102 HYPRE_ParCSRPCGGetNumIterations(HYSolver_, &tempIter);
5103 numIterations += tempIter;
5104 HYPRE_ParVectorCopy( b_csr, r_csr );
5105 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5106 rnorm1 = rnorm;
5107 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5108 rnorm = sqrt( rnorm );
5109 convRate = rnorm / rnorm1;
5110 }
5111 if (rnorm/rnorm0 < tolerance_) retFlag = 0;
5112 else if (numIterations < maxIterations_)
5113 {
5114 HYPRE_ParCSRPCGSetMaxIter(HYSolver_,maxIterations_-numIterations);
5115 retFlag = HYPRE_ParCSRPCGSolve(HYSolver_, A_csr, b_csr, x_csr);
5116 HYPRE_ParCSRPCGGetNumIterations(HYSolver_, &tempIter);
5117 numIterations += tempIter;
5118 }
5119 else retFlag = 1;
5120 }
5121 if ( retFlag != 0 )
5122 {
5123 printf("HYPRE_LSC::launchSolver ERROR : in PCG solve.\n");
5124 return retFlag;
5125 }
5126 switch ( projectionScheme_ )
5127 {
5128 case 1 : addToAConjProjectionSpace(currX_,currB_); break;
5129 case 2 : addToMinResProjectionSpace(currX_,currB_); break;
5130 }
5131 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5132 break;
5133
5134 //----------------------------------------------------------------
5135 // choose LSICG
5136 //----------------------------------------------------------------
5137
5138 case HYLSICG :
5139 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5140 {
5141 printf("***************************************************\n");
5142 printf("* Conjugate Gradient solver \n");
5143 printf("* maximum no. of iterations = %d\n", maxIterations_);
5144 printf("* convergence tolerance = %e\n", tolerance_);
5145 printf("*--------------------------------------------------\n");
5146 }
5147 setupLSICGPrecon();
5148 HYPRE_ParCSRLSICGSetMaxIter(HYSolver_, maxIterations_);
5149 HYPRE_ParCSRLSICGSetTol(HYSolver_, tolerance_);
5150 if (normAbsRel_ == 0) HYPRE_ParCSRLSICGSetStopCrit(HYSolver_,0);
5151 else HYPRE_ParCSRLSICGSetStopCrit(HYSolver_,1);
5152 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5153 {
5154 if ( mypid_ == 0 )
5155 printf("***************************************************\n");
5156 HYPRE_ParCSRLSICGSetLogging(HYSolver_, 1);
5157 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 )
5158 HYPRE_ParCSRLSICGSetLogging(HYSolver_, 2);
5159 }
5160 retFlag = HYPRE_ParCSRLSICGSetup(HYSolver_, A_csr, b_csr, x_csr);
5161 if ( retFlag != 0 )
5162 {
5163 printf("HYPRE_LSC::launchSolver ERROR : in LSICG setup.\n");
5164 return retFlag;
5165 }
5166 MPI_Barrier( comm_ );
5167 ptime = LSC_Wtime();
5168 retFlag = HYPRE_ParCSRLSICGSolve(HYSolver_, A_csr, b_csr, x_csr);
5169 if ( retFlag != 0 )
5170 {
5171 printf("HYPRE_LSC::launchSolver ERROR : in LSICG solve.\n");
5172 return retFlag;
5173 }
5174 HYPRE_ParCSRLSICGGetNumIterations(HYSolver_, &numIterations);
5175 HYPRE_ParVectorCopy( b_csr, r_csr );
5176 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5177 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5178 rnorm = sqrt( rnorm );
5179 switch ( projectionScheme_ )
5180 {
5181 case 1 : addToAConjProjectionSpace(currX_,currB_); break;
5182 case 2 : addToMinResProjectionSpace(currX_,currB_); break;
5183 }
5184 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5185 break;
5186
5187 //----------------------------------------------------------------
5188 // choose hybrid method : CG with diagonal/BoomerAMG preconditioner
5189 //----------------------------------------------------------------
5190
5191 case HYHYBRID :
5192 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5193 {
5194 printf("***************************************************\n");
5195 printf("* PCG with hybrid diagonal/BoomerAMG preconditioner\n");
5196 printf("* maximum no. of iterations = %d\n", maxIterations_);
5197 printf("* convergence tolerance = %e\n", tolerance_);
5198 printf("*--------------------------------------------------\n");
5199 }
5200 HYPRE_ParCSRHybridSetPCGMaxIter(HYSolver_, maxIterations_);
5201 HYPRE_ParCSRHybridSetTol(HYSolver_, tolerance_);
5202 HYPRE_ParCSRHybridSetRelChange(HYSolver_, 0);
5203 HYPRE_ParCSRHybridSetTwoNorm(HYSolver_, 1);
5204 HYPRE_ParCSRHybridSetConvergenceTol(HYSolver_, 0.9);
5205 HYPRE_ParCSRHybridSetDSCGMaxIter(HYSolver_, 20);
5206 if ( HYOutputLevel_ & HYFEI_AMGDEBUG )
5207 HYPRE_ParCSRHybridSetPrintLevel(HYSolver_, 32);
5208 HYPRE_ParCSRHybridSetCoarsenType(HYSolver_, amgCoarsenType_);
5209 HYPRE_ParCSRHybridSetMeasureType(HYSolver_, amgMeasureType_);
5210 HYPRE_ParCSRHybridSetStrongThreshold(HYSolver_,amgStrongThreshold_);
5211 numSweeps = hypre_CTAlloc(int,4,HYPRE_MEMORY_HOST);
5212 for ( i = 0; i < 4; i++ ) numSweeps[i] = amgNumSweeps_[i];
5213 HYPRE_ParCSRHybridSetNumGridSweeps(HYSolver_, numSweeps);
5214 relaxType = hypre_CTAlloc(int,4,HYPRE_MEMORY_HOST);
5215 for ( i = 0; i < 4; i++ ) relaxType[i] = amgRelaxType_[i];
5216 HYPRE_ParCSRHybridSetGridRelaxType(HYSolver_, relaxType);
5217 relaxWt = hypre_CTAlloc(double,25,HYPRE_MEMORY_HOST);
5218 for ( i = 0; i < 25; i++ ) relaxWt[i] = amgRelaxWeight_[i];
5219 HYPRE_ParCSRHybridSetRelaxWeight(HYSolver_, relaxWt);
5220 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5221 {
5222 if ( mypid_ == 0 )
5223 printf("***************************************************\n");
5224 HYPRE_ParCSRHybridSetLogging(HYSolver_, 1);
5225 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 )
5226 HYPRE_ParCSRHybridSetLogging(HYSolver_, 2);
5227 }
5228 retFlag = HYPRE_ParCSRHybridSetup(HYSolver_, A_csr, b_csr, x_csr);
5229 if ( retFlag != 0 )
5230 {
5231 printf("HYPRE_LSC::launchSolver ERROR : in Hybrid setup.\n");
5232 return retFlag;
5233 }
5234 MPI_Barrier( comm_ );
5235 ptime = LSC_Wtime();
5236 retFlag = HYPRE_ParCSRHybridSolve(HYSolver_, A_csr, b_csr, x_csr);
5237 if ( retFlag != 0 )
5238 {
5239 printf("HYPRE_LSC::launchSolver ERROR : in Hybrid solve.\n");
5240 return retFlag;
5241 }
5242 HYPRE_ParCSRHybridGetNumIterations(HYSolver_, &numIterations);
5243 HYPRE_ParVectorCopy( b_csr, r_csr );
5244 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5245 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5246 rnorm = sqrt( rnorm );
5247 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5248 break;
5249
5250 //----------------------------------------------------------------
5251 // choose GMRES
5252 //----------------------------------------------------------------
5253
5254 case HYGMRES :
5255 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5256 {
5257 printf("***************************************************\n");
5258 printf("* Generalized Minimal Residual (GMRES) solver \n");
5259 printf("* restart size = %d\n", gmresDim_);
5260 printf("* maximum no. of iterations = %d\n", maxIterations_);
5261 printf("* convergence tolerance = %e\n", tolerance_);
5262 printf("*--------------------------------------------------\n");
5263 }
5264 setupGMRESPrecon();
5265 HYPRE_ParCSRGMRESSetKDim(HYSolver_, gmresDim_);
5266 HYPRE_ParCSRGMRESSetMaxIter(HYSolver_, maxIterations_);
5267 if ( normAbsRel_ == 0 )
5268 {
5269 HYPRE_GMRESSetTol(HYSolver_, tolerance_);
5270 }
5271 else
5272 {
5273 HYPRE_GMRESSetTol(HYSolver_, 0.0);
5274 HYPRE_GMRESSetAbsoluteTol(HYSolver_, tolerance_);
5275 }
5276 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5277 {
5278 HYPRE_ParCSRGMRESSetPrintLevel(HYSolver_, 1);
5279 if ( mypid_ == 0 )
5280 printf("***************************************************\n");
5281 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 )
5282 HYPRE_ParCSRGMRESSetPrintLevel(HYSolver_, 2);
5283 }
5284 retFlag = HYPRE_ParCSRGMRESSetup(HYSolver_, A_csr, b_csr, x_csr);
5285 if ( retFlag != 0 )
5286 {
5287 printf("HYPRE_LSC::launchSolver ERROR : in GMRES setup.\n");
5288 return retFlag;
5289 }
5290 // if variable mli preconditioner (SA and GSA)
5291 if ( MLI_Hybrid_GSA_ && HYPreconID_ == HYMLI )
5292 {
5293 if ( normAbsRel_ == 0 )
5294 {
5295 HYPRE_ParVectorCopy( b_csr, r_csr );
5296 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5297 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm0);
5298 rnorm0 = sqrt( rnorm0 );
5299 }
5300 else rnorm0 = 1.0;
5301 HYPRE_ParCSRGMRESSetMaxIter(HYSolver_, MLI_Hybrid_MaxIter_);
5302 rateThresh = 1.0;
5303 for ( i = 0; i < MLI_Hybrid_MaxIter_; i++ )
5304 rateThresh *= MLI_Hybrid_ConvRate_;
5305 }
5306 MPI_Barrier( comm_ );
5307 ptime = LSC_Wtime();
5308 retFlag = HYPRE_ParCSRGMRESSolve(HYSolver_, A_csr, b_csr, x_csr);
5309 HYPRE_ParCSRGMRESGetNumIterations(HYSolver_, &numIterations);
5310 HYPRE_ParVectorCopy( b_csr, r_csr );
5311 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5312 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5313 rnorm = sqrt( rnorm );
5314 // if variable mli preconditioner (SA and GSA)
5315 if ( MLI_Hybrid_GSA_ && HYPreconID_ == HYMLI )
5316 {
5317 nTrials = 0;
5318 //if (rnorm/rnorm0 >= tolerance_)
5319 //{
5320 // HYPRE_ParCSRGMRESSolve(HYSolver_, A_csr, b_csr, x_csr);
5321 // HYPRE_ParCSRGMRESGetNumIterations(HYSolver_, &tempIter);
5322 // numIterations += tempIter;
5323 // HYPRE_ParVectorCopy( b_csr, r_csr );
5324 // HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5325 // HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm1);
5326 // rnorm1 = sqrt( rnorm1 );
5327 // convRate = rnorm1 / rnorm;
5328 // rnorm = rnorm1;
5329 //}
5330 convRate = rnorm / rnorm0;
5331 while ((rnorm/rnorm0)>=tolerance_ && nTrials<MLI_Hybrid_NTrials_)
5332 {
5333 nTrials++;
5334 if ( convRate > rateThresh )
5335 {
5336 if ( MLI_Hybrid_NSIncr_ > 1 )
5337 sprintf(paramString, "MLI incrNullSpaceDim %d",
5338 MLI_Hybrid_NSIncr_);
5339 else
5340 sprintf(paramString, "MLI incrNullSpaceDim 2");
5341 HYPRE_LSI_MLISetParams(HYPrecon_, paramString);
5342 HYPRE_ParCSRGMRESSetup(HYSolver_, A_csr, b_csr, x_csr);
5343 }
5344 HYPRE_ParCSRGMRESSolve(HYSolver_, A_csr, b_csr, x_csr);
5345 HYPRE_ParCSRGMRESGetNumIterations(HYSolver_, &tempIter);
5346 numIterations += tempIter;
5347 HYPRE_ParVectorCopy( b_csr, r_csr );
5348 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5349 rnorm1 = rnorm;
5350 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5351 rnorm = sqrt( rnorm );
5352 convRate = rnorm / rnorm1;
5353 }
5354 if (rnorm/rnorm0 < tolerance_) retFlag = 0;
5355 else if (numIterations < maxIterations_)
5356 {
5357 HYPRE_ParCSRGMRESSetMaxIter(HYSolver_,
5358 maxIterations_-numIterations);
5359 retFlag = HYPRE_ParCSRGMRESSolve(HYSolver_,A_csr,b_csr,x_csr);
5360 HYPRE_ParCSRGMRESGetNumIterations(HYSolver_, &tempIter);
5361 numIterations += tempIter;
5362 }
5363 else retFlag = 1;
5364 }
5365 if ( retFlag != 0 )
5366 {
5367 printf("HYPRE_LSC::launchSolver ERROR : in GMRES solve.\n");
5368 return retFlag;
5369 }
5370 switch ( projectionScheme_ )
5371 {
5372 case 1 : addToAConjProjectionSpace(currX_,currB_); break;
5373 case 2 : addToMinResProjectionSpace(currX_,currB_); break;
5374 }
5375 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5376 break;
5377
5378 //----------------------------------------------------------------
5379 // choose flexible GMRES
5380 //----------------------------------------------------------------
5381
5382 case HYFGMRES :
5383 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5384 {
5385 printf("***************************************************\n");
5386 printf("* Flexible GMRES solver \n");
5387 printf("* restart size = %d\n", gmresDim_);
5388 printf("* maximum no. of iterations = %d\n", maxIterations_);
5389 printf("* convergence tolerance = %e\n", tolerance_);
5390 printf("*--------------------------------------------------\n");
5391 }
5392 setupFGMRESPrecon();
5393 HYPRE_ParCSRFGMRESSetKDim(HYSolver_, gmresDim_);
5394 HYPRE_ParCSRFGMRESSetMaxIter(HYSolver_, maxIterations_);
5395 HYPRE_ParCSRFGMRESSetTol(HYSolver_, tolerance_);
5396 if ( normAbsRel_ == 0 ) HYPRE_ParCSRFGMRESSetStopCrit(HYSolver_,0);
5397 else HYPRE_ParCSRFGMRESSetStopCrit(HYSolver_,1);
5398 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5399 {
5400 if ( mypid_ == 0 )
5401 printf("***************************************************\n");
5402 HYPRE_ParCSRFGMRESSetLogging(HYSolver_, 1);
5403 }
5404
5405 retFlag = HYPRE_ParCSRFGMRESSetup(HYSolver_, A_csr, b_csr, x_csr);
5406 if ( retFlag != 0 )
5407 {
5408 printf("HYPRE_LSC::launchSolver ERROR : in FGMRES setup.\n");
5409 return retFlag;
5410 }
5411
5412 if ( fgmresUpdateTol_ && HYPreconID_ == HYBLOCK )
5413 HYPRE_ParCSRFGMRESUpdatePrecondTolerance(HYSolver_,
5414 HYPRE_LSI_BlockPrecondSetA11Tolerance);
5415
5416 MPI_Barrier( comm_ );
5417 ptime = LSC_Wtime();
5418 retFlag = HYPRE_ParCSRFGMRESSolve(HYSolver_, A_csr, b_csr, x_csr);
5419 if ( retFlag != 0 )
5420 {
5421 printf("HYPRE_LSC::launchSolver ERROR : in FGMRES solve.\n");
5422 return retFlag;
5423 }
5424 HYPRE_ParCSRFGMRESGetNumIterations(HYSolver_, &numIterations);
5425 HYPRE_ParVectorCopy( b_csr, r_csr );
5426 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5427 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5428 rnorm = sqrt( rnorm );
5429 switch ( projectionScheme_ )
5430 {
5431 case 1 : addToAConjProjectionSpace(currX_,currB_); break;
5432 case 2 : addToMinResProjectionSpace(currX_,currB_); break;
5433 }
5434 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5435 break;
5436
5437 //----------------------------------------------------------------
5438 // choose BiCGSTAB
5439 //----------------------------------------------------------------
5440
5441 case HYCGSTAB :
5442 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5443 {
5444 printf("***************************************************\n");
5445 printf("* BiCGSTAB solver \n");
5446 printf("* maximum no. of iterations = %d\n", maxIterations_);
5447 printf("* convergence tolerance = %e\n", tolerance_);
5448 printf("*--------------------------------------------------\n");
5449 }
5450 setupBiCGSTABPrecon();
5451 HYPRE_ParCSRBiCGSTABSetMaxIter(HYSolver_, maxIterations_);
5452 HYPRE_ParCSRBiCGSTABSetTol(HYSolver_, tolerance_);
5453 if ( normAbsRel_ == 0 ) HYPRE_ParCSRBiCGSTABSetStopCrit(HYSolver_,0);
5454 else HYPRE_ParCSRBiCGSTABSetStopCrit(HYSolver_,1);
5455 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5456 {
5457 HYPRE_ParCSRBiCGSTABSetPrintLevel(HYSolver_, 1);
5458 if ( mypid_ == 0 )
5459 printf("***************************************************\n");
5460 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 )
5461 HYPRE_ParCSRBiCGSTABSetPrintLevel(HYSolver_, 2);
5462 }
5463 retFlag = HYPRE_ParCSRBiCGSTABSetup(HYSolver_, A_csr, b_csr, x_csr);
5464 if ( retFlag != 0 )
5465 {
5466 printf("HYPRE_LSC::launchSolver ERROR : in BiCGSTAB setup.\n");
5467 return retFlag;
5468 }
5469 MPI_Barrier( comm_ );
5470 ptime = LSC_Wtime();
5471 retFlag = HYPRE_ParCSRBiCGSTABSolve(HYSolver_, A_csr, b_csr, x_csr);
5472 if ( retFlag != 0 )
5473 {
5474 printf("HYPRE_LSC::launchSolver ERROR : in BiCGSTAB solve.\n");
5475 return retFlag;
5476 }
5477 HYPRE_ParCSRBiCGSTABGetNumIterations(HYSolver_, &numIterations);
5478 HYPRE_ParVectorCopy( b_csr, r_csr );
5479 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5480 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5481 rnorm = sqrt( rnorm );
5482 switch ( projectionScheme_ )
5483 {
5484 case 1 : addToAConjProjectionSpace(currX_,currB_); break;
5485 case 2 : addToMinResProjectionSpace(currX_,currB_); break;
5486 }
5487 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5488 break;
5489
5490 //----------------------------------------------------------------
5491 // choose BiCGSTABL
5492 //----------------------------------------------------------------
5493
5494 case HYCGSTABL :
5495 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5496 {
5497 printf("***************************************************\n");
5498 printf("* BiCGSTAB(2) solver \n");
5499 printf("* maximum no. of iterations = %d\n", maxIterations_);
5500 printf("* convergence tolerance = %e\n", tolerance_);
5501 printf("*--------------------------------------------------\n");
5502 }
5503 setupBiCGSTABLPrecon();
5504 HYPRE_ParCSRBiCGSTABLSetMaxIter(HYSolver_, maxIterations_);
5505 HYPRE_ParCSRBiCGSTABLSetTol(HYSolver_, tolerance_);
5506 if (normAbsRel_ == 0) HYPRE_ParCSRBiCGSTABLSetStopCrit(HYSolver_,0);
5507 else HYPRE_ParCSRBiCGSTABLSetStopCrit(HYSolver_,1);
5508 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5509 {
5510 if ( mypid_ == 0 )
5511 printf("***************************************************\n");
5512 HYPRE_ParCSRBiCGSTABLSetLogging(HYSolver_, 1);
5513 }
5514 retFlag = HYPRE_ParCSRBiCGSTABLSetup(HYSolver_, A_csr, b_csr, x_csr);
5515 if ( retFlag != 0 )
5516 {
5517 printf("HYPRE_LSC::launchSolver ERROR : in BiCGSTABL setup.\n");
5518 return retFlag;
5519 }
5520 MPI_Barrier( comm_ );
5521 ptime = LSC_Wtime();
5522 retFlag = HYPRE_ParCSRBiCGSTABLSolve(HYSolver_, A_csr, b_csr, x_csr);
5523 if ( retFlag != 0 )
5524 {
5525 printf("HYPRE_LSC::launchSolver ERROR : in BiCGSTABL solve.\n");
5526 return retFlag;
5527 }
5528 HYPRE_ParCSRBiCGSTABLGetNumIterations(HYSolver_, &numIterations);
5529 HYPRE_ParVectorCopy( b_csr, r_csr );
5530 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5531 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5532 rnorm = sqrt( rnorm );
5533 switch ( projectionScheme_ )
5534 {
5535 case 1 : addToAConjProjectionSpace(currX_,currB_); break;
5536 case 2 : addToMinResProjectionSpace(currX_,currB_); break;
5537 }
5538 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5539 break;
5540
5541 //----------------------------------------------------------------
5542 // choose TFQMR
5543 //----------------------------------------------------------------
5544
5545 case HYTFQMR :
5546 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5547 {
5548 printf("***************************************************\n");
5549 printf("* TFQMR solver \n");
5550 printf("* maximum no. of iterations = %d\n", maxIterations_);
5551 printf("* convergence tolerance = %e\n", tolerance_);
5552 printf("*--------------------------------------------------\n");
5553 }
5554 setupTFQmrPrecon();
5555 HYPRE_ParCSRTFQmrSetMaxIter(HYSolver_, maxIterations_);
5556 HYPRE_ParCSRTFQmrSetTol(HYSolver_, tolerance_);
5557 if ( normAbsRel_ == 0 ) HYPRE_ParCSRTFQmrSetStopCrit(HYSolver_,0);
5558 else HYPRE_ParCSRTFQmrSetStopCrit(HYSolver_,1);
5559 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5560 {
5561 if ( mypid_ == 0 )
5562 printf("***************************************************\n");
5563 HYPRE_ParCSRTFQmrSetLogging(HYSolver_, 1);
5564 }
5565 retFlag = HYPRE_ParCSRTFQmrSetup(HYSolver_, A_csr, b_csr, x_csr);
5566 if ( retFlag != 0 )
5567 {
5568 printf("HYPRE_LSC::launchSolver ERROR : in TFQMR setup.\n");
5569 return retFlag;
5570 }
5571 MPI_Barrier( comm_ );
5572 ptime = LSC_Wtime();
5573 retFlag = HYPRE_ParCSRTFQmrSolve(HYSolver_, A_csr, b_csr, x_csr);
5574 if ( retFlag != 0 )
5575 {
5576 printf("HYPRE_LSC::launchSolver ERROR : in TFQMR solve.\n");
5577 return retFlag;
5578 }
5579 HYPRE_ParCSRTFQmrGetNumIterations(HYSolver_, &numIterations);
5580 HYPRE_ParVectorCopy( b_csr, r_csr );
5581 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5582 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5583 rnorm = sqrt( rnorm );
5584 switch ( projectionScheme_ )
5585 {
5586 case 1 : addToAConjProjectionSpace(currX_,currB_); break;
5587 case 2 : addToMinResProjectionSpace(currX_,currB_); break;
5588 }
5589 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5590 break;
5591
5592 //----------------------------------------------------------------
5593 // choose BiCGS
5594 //----------------------------------------------------------------
5595
5596 case HYBICGS :
5597 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5598 {
5599 printf("***************************************************\n");
5600 printf("* BiCGS solver \n");
5601 printf("* maximum no. of iterations = %d\n", maxIterations_);
5602 printf("* convergence tolerance = %e\n", tolerance_);
5603 printf("*--------------------------------------------------\n");
5604 }
5605 setupBiCGSPrecon();
5606 HYPRE_ParCSRBiCGSSetMaxIter(HYSolver_, maxIterations_);
5607 HYPRE_ParCSRBiCGSSetTol(HYSolver_, tolerance_);
5608 if ( normAbsRel_ == 0 ) HYPRE_ParCSRBiCGSSetStopCrit(HYSolver_,0);
5609 else HYPRE_ParCSRBiCGSSetStopCrit(HYSolver_,1);
5610 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5611 {
5612 if ( mypid_ == 0 )
5613 printf("***************************************************\n");
5614 HYPRE_ParCSRBiCGSSetLogging(HYSolver_, 1);
5615 }
5616 retFlag = HYPRE_ParCSRBiCGSSetup(HYSolver_, A_csr, b_csr, x_csr);
5617 if ( retFlag != 0 )
5618 {
5619 printf("HYPRE_LSC::launchSolver ERROR : in CGS setup.\n");
5620 return retFlag;
5621 }
5622 MPI_Barrier( comm_ );
5623 ptime = LSC_Wtime();
5624 retFlag = HYPRE_ParCSRBiCGSSolve(HYSolver_, A_csr, b_csr, x_csr);
5625 if ( retFlag != 0 )
5626 {
5627 printf("HYPRE_LSC::launchSolver ERROR : in CGS solve.\n");
5628 return retFlag;
5629 }
5630 HYPRE_ParCSRBiCGSGetNumIterations(HYSolver_, &numIterations);
5631 HYPRE_ParVectorCopy( b_csr, r_csr );
5632 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5633 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5634 rnorm = sqrt( rnorm );
5635 switch ( projectionScheme_ )
5636 {
5637 case 1 : addToAConjProjectionSpace(currX_,currB_); break;
5638 case 2 : addToMinResProjectionSpace(currX_,currB_); break;
5639 }
5640 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5641 break;
5642
5643 //----------------------------------------------------------------
5644 // choose Symmetric QMR
5645 //----------------------------------------------------------------
5646
5647 case HYSYMQMR :
5648 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5649 {
5650 printf("***************************************************\n");
5651 printf("* SymQMR solver (for symmetric matrices and precond) \n");
5652 printf("* maximum no. of iterations = %d\n", maxIterations_);
5653 printf("* convergence tolerance = %e\n", tolerance_);
5654 printf("*--------------------------------------------------\n");
5655 }
5656 setupSymQMRPrecon();
5657 HYPRE_ParCSRSymQMRSetMaxIter(HYSolver_, maxIterations_);
5658 HYPRE_ParCSRSymQMRSetTol(HYSolver_, tolerance_);
5659 if ( normAbsRel_ == 0 ) HYPRE_ParCSRSymQMRSetStopCrit(HYSolver_,0);
5660 else HYPRE_ParCSRSymQMRSetStopCrit(HYSolver_,1);
5661 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5662 {
5663 if ( mypid_ == 0 )
5664 printf("***************************************************\n");
5665 HYPRE_ParCSRSymQMRSetLogging(HYSolver_, 1);
5666 }
5667 retFlag = HYPRE_ParCSRSymQMRSetup(HYSolver_, A_csr, b_csr, x_csr);
5668 if ( retFlag != 0 )
5669 {
5670 printf("HYPRE_LSC::launchSolver ERROR : in SymQMR setup.\n");
5671 return retFlag;
5672 }
5673 MPI_Barrier( comm_ );
5674 ptime = LSC_Wtime();
5675 retFlag = HYPRE_ParCSRSymQMRSolve(HYSolver_, A_csr, b_csr, x_csr);
5676 if ( retFlag != 0 )
5677 {
5678 printf("HYPRE_LSC::launchSolver ERROR : in SymQMR solve.\n");
5679 return retFlag;
5680 }
5681 HYPRE_ParCSRSymQMRGetNumIterations(HYSolver_, &numIterations);
5682 HYPRE_ParVectorCopy( b_csr, r_csr );
5683 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5684 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5685 rnorm = sqrt( rnorm );
5686 switch ( projectionScheme_ )
5687 {
5688 case 1 : addToAConjProjectionSpace(currX_,currB_); break;
5689 case 2 : addToMinResProjectionSpace(currX_,currB_); break;
5690 }
5691 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5692 break;
5693
5694 //----------------------------------------------------------------
5695 // choose Boomeramg
5696 //----------------------------------------------------------------
5697
5698 case HYAMG :
5699 solveUsingBoomeramg(status);
5700 HYPRE_BoomerAMGGetNumIterations(HYSolver_, &numIterations);
5701 HYPRE_ParVectorCopy( b_csr, r_csr );
5702 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5703 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5704 rnorm = sqrt( rnorm );
5705 if ( numIterations >= maxIterations_ ) status = 1; else status = 0;
5706 ptime = stime;
5707 //printf("Boomeramg solver - return status = %d\n",status);
5708 break;
5709
5710 //----------------------------------------------------------------
5711 // choose SuperLU (single processor)
5712 //----------------------------------------------------------------
5713
5714 case HYSUPERLU :
5715 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5716 {
5717 printf("***************************************************\n");
5718 printf("* SuperLU (sequential) solver \n");
5719 printf("*--------------------------------------------------\n");
5720 }
5721 rnorm = solveUsingSuperLU(status);
5722 #ifndef NOFEI
5723 if ( status == 1 ) status = 0;
5724 #endif
5725 numIterations = 1;
5726 ptime = stime;
5727 //printf("SuperLU solver - return status = %d\n",status);
5728 break;
5729
5730 //----------------------------------------------------------------
5731 // choose SuperLU (single processor)
5732 //----------------------------------------------------------------
5733
5734 case HYSUPERLUX :
5735 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5736 {
5737 printf("***************************************************\n");
5738 printf("* SuperLU (sequential) solver with refinement \n");
5739 printf("*--------------------------------------------------\n");
5740 }
5741 rnorm = solveUsingSuperLUX(status);
5742 #ifndef NOFEI
5743 if ( status == 1 ) status = 0;
5744 #endif
5745 numIterations = 1;
5746 //printf("SuperLUX solver - return status = %d\n",status);
5747 break;
5748
5749 //----------------------------------------------------------------
5750 // choose distributed SuperLU
5751 //----------------------------------------------------------------
5752
5753 case HYDSUPERLU :
5754 #ifdef HYPRE_USING_DSUPERLU
5755 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5756 {
5757 printf("***************************************************\n");
5758 printf("* distributed SuperLU solver \n");
5759 printf("*--------------------------------------------------\n");
5760 }
5761 rnorm = solveUsingDSuperLU(status);
5762 #ifndef NOFEI
5763 if ( status == 1 ) status = 0;
5764 #endif
5765 numIterations = 1;
5766 #else
5767 printf("distributed SuperLU not available.\n");
5768 exit(1);
5769 #endif
5770 break;
5771
5772 //----------------------------------------------------------------
5773 // choose Y12M (single processor)
5774 //----------------------------------------------------------------
5775
5776 case HYY12M :
5777 #ifdef Y12M
5778 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5779 {
5780 printf("***************************************************\n");
5781 printf("* Y12M (sequential) solver\n");
5782 printf("*--------------------------------------------------\n");
5783 }
5784 solveUsingY12M(status);
5785 #ifndef NOFEI
5786 if ( status == 1 ) status = 0;
5787 #endif
5788 numIterations = 1;
5789 ptime = stime;
5790 //printf("Y12M solver - return status = %d\n",status);
5791 break;
5792 #else
5793 printf("HYPRE_LSC : Y12M not available. \n");
5794 exit(1);
5795 break;
5796 #endif
5797
5798 //----------------------------------------------------------------
5799 // choose AMGE (single processor)
5800 //----------------------------------------------------------------
5801
5802 case HYAMGE :
5803 #ifdef HAVE_AMGE
5804 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5805 {
5806 printf("***************************************************\n");
5807 printf("* AMGe (sequential) solver\n");
5808 printf("*--------------------------------------------------\n");
5809 }
5810 solveUsingAMGe(numIterations);
5811 if ( numIterations >= maxIterations_ ) status = 1;
5812 ptime = stime;
5813 #else
5814 printf("AMGe not supported.\n");
5815 #endif
5816 break;
5817 }
5818
5819 //-------------------------------------------------------------------
5820 // recover solution for reduced system
5821 //-------------------------------------------------------------------
5822
5823 if ( slideReduction_ == 1 )
5824 {
5825 newnorm = rnorm;
5826 rnorm = buildSlideReducedSoln();
5827 }
5828 else if ( slideReduction_ == 2 )
5829 {
5830 newnorm = rnorm;
5831 rnorm = buildSlideReducedSoln2();
5832 }
5833 else if ( slideReduction_ == 3 || slideReduction_ == 4 )
5834 {
5835 newnorm = rnorm;
5836 currA_ = TempA;
5837 currX_ = TempX;
5838 currB_ = TempB;
5839 currR_ = TempR;
5840 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
5841 HYPRE_IJVectorGetObject(currX_, (void **) &x_csr);
5842 HYPRE_IJVectorGetObject(currB_, (void **) &b_csr);
5843 HYPRE_IJVectorGetObject(currR_, (void **) &r_csr);
5844 if ( slideReduction_ == 3 )
5845 slideObj->buildReducedSolnVector(currX_, currB_);
5846 else slideObj->buildModifiedSolnVector(currX_);
5847 HYPRE_ParVectorCopy( b_csr, r_csr );
5848 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5849 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5850 rnorm = sqrt( rnorm );
5851 }
5852 else if ( schurReduction_ == 1 )
5853 {
5854 newnorm = rnorm;
5855 rnorm = buildSchurReducedSoln();
5856 }
5857 if ( (normalEqnFlag_ & 7) == 7 )
5858 {
5859 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
5860 HYPRE_IJVectorGetObject(currX_, (void **) &x_csr);
5861 HYPRE_IJVectorGetObject(currB_, (void **) &b_csr);
5862 HYPRE_IJVectorGetObject(currR_, (void **) &r_csr);
5863 HYPRE_ParVectorCopy( b_csr, r_csr );
5864 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5865 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5866 rnorm = sqrt( rnorm );
5867 }
5868
5869 //-------------------------------------------------------------------
5870 // register solver return information and print timing information
5871 //-------------------------------------------------------------------
5872
5873 solveStatus = status;
5874 iterations = numIterations;
5875 rnorm_ = rnorm;
5876
5877 MPI_Barrier(comm_);
5878 etime = LSC_Wtime();
5879 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
5880 {
5881 printf("***************************************************\n");
5882 printf("* Solver Statistics *\n");
5883 printf("*-------------------------------------------------*\n");
5884 if ( slideReduction_ || schurReduction_ )
5885 printf("** HYPRE matrix reduction time = %e\n",rtime2-rtime1);
5886 printf("** HYPRE preconditioner setup time = %e\n", ptime - stime);
5887 printf("** HYPRE solution time = %e\n", etime - ptime);
5888 printf("** HYPRE total time = %e\n", etime - stime);
5889 printf("** HYPRE number of iterations = %d\n", numIterations);
5890 if ( slideReduction_ || schurReduction_ )
5891 printf("** HYPRE reduced residual norm = %e\n", newnorm);
5892 printf("** HYPRE final residual norm = %e\n", rnorm);
5893 printf("***************************************************\n");
5894 }
5895
5896 //-------------------------------------------------------------------
5897 // write solution to an output file
5898 //-------------------------------------------------------------------
5899
5900 if ( HYOutputLevel_ & HYFEI_PRINTSOL )
5901 {
5902 nrows = localEndRow_ - localStartRow_ + 1;
5903 startRow = localStartRow_ - 1;
5904 sprintf(fname, "hypre_sol.out.%d", mypid_);
5905 fp = fopen( fname, "w");
5906 fprintf(fp, "%6d \n", nrows);
5907 for ( i = startRow; i < startRow+nrows; i++ )
5908 {
5909 HYPRE_IJVectorGetValues(currX_, 1, &i, &ddata);
5910 fprintf(fp, "%6d %25.16e \n", i+1, ddata);
5911 }
5912 fclose(fp);
5913 MPI_Barrier(comm_);
5914 }
5915
5916 //-------------------------------------------------------------------
5917 // diagnostic message
5918 //-------------------------------------------------------------------
5919
5920 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
5921 printf("%4d : HYPRE_LSC::leaving launchSolver.\n", mypid_);
5922 return (0);
5923 }
5924
5925 //***************************************************************************
5926 // this function extracts the matrix in a CSR format
5927 //---------------------------------------------------------------------------
5928
writeSystem(const char * name)5929 int HYPRE_LinSysCore::writeSystem(const char *name)
5930 {
5931 (void) name;
5932 printf("HYPRE_LinsysCore : writeSystem not implemented.\n");
5933 return (0);
5934 }
5935
5936 //***************************************************************************
5937 // this function computes matrix vector product
5938 //---------------------------------------------------------------------------
5939
HYPRE_LSC_Matvec(void * x,void * y)5940 int HYPRE_LinSysCore::HYPRE_LSC_Matvec(void *x, void *y)
5941 {
5942 HYPRE_ParCSRMatrix A_csr;
5943 HYPRE_ParVector x_csr = (HYPRE_ParVector) x;
5944 HYPRE_ParVector y_csr = (HYPRE_ParVector) y;
5945 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
5946 HYPRE_ParCSRMatrixMatvec(1.0, A_csr, x_csr, 0.0, y_csr);
5947 return (0);
5948 }
5949
5950 //***************************************************************************
5951 // this function computes vector multiply and add
5952 //---------------------------------------------------------------------------
5953
HYPRE_LSC_Axpby(double a,void * x,double b,void * y)5954 int HYPRE_LinSysCore::HYPRE_LSC_Axpby(double a, void *x, double b, void *y)
5955 {
5956 HYPRE_ParVector x_csr = (HYPRE_ParVector) x;
5957 HYPRE_ParVector y_csr = (HYPRE_ParVector) y;
5958 if ( b != 1.0 ) HYPRE_ParVectorScale( b, y_csr);
5959 hypre_ParVectorAxpy(a, (hypre_ParVector*) x_csr, (hypre_ParVector*) y_csr);
5960 return (0);
5961 }
5962
5963 //***************************************************************************
5964 // this function fetches the right hand side vector
5965 //---------------------------------------------------------------------------
5966
HYPRE_LSC_GetRHSVector()5967 void *HYPRE_LinSysCore::HYPRE_LSC_GetRHSVector()
5968 {
5969 HYPRE_ParVector b_csr;
5970 HYPRE_IJVectorGetObject(HYb_, (void **) &b_csr);
5971 return (void *) b_csr;
5972 }
5973
5974 //***************************************************************************
5975 // this function fetches the solution vector
5976 //---------------------------------------------------------------------------
5977
HYPRE_LSC_GetSolVector()5978 void *HYPRE_LinSysCore::HYPRE_LSC_GetSolVector()
5979 {
5980 HYPRE_ParVector x_csr;
5981 HYPRE_IJVectorGetObject(HYx_, (void **) &x_csr);
5982 return (void *) x_csr;
5983 }
5984
5985 //***************************************************************************
5986 // this function fetches the matrix
5987 //---------------------------------------------------------------------------
5988
HYPRE_LSC_GetMatrix()5989 void *HYPRE_LinSysCore::HYPRE_LSC_GetMatrix()
5990 {
5991 HYPRE_ParCSRMatrix A_csr;
5992 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
5993 return (void *) A_csr;
5994 }
5995
5996 //***************************************************************************
5997 // this function set column ranges
5998 //---------------------------------------------------------------------------
5999
HYPRE_LSC_SetColMap(int start,int end)6000 void *HYPRE_LinSysCore::HYPRE_LSC_SetColMap(int start, int end)
6001 {
6002 localStartCol_ = start;
6003 localEndCol_ = end;
6004 return (void *) NULL;
6005 }
6006
6007 //***************************************************************************
6008 // this function set column ranges
6009 //---------------------------------------------------------------------------
6010
HYPRE_LSC_MatMatMult(void * inMat)6011 void *HYPRE_LinSysCore::HYPRE_LSC_MatMatMult(void *inMat)
6012 {
6013 HYPRE_ParCSRMatrix A_csr;
6014 hypre_ParCSRMatrix *B_csr, *C_csr;
6015 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
6016 B_csr = (hypre_ParCSRMatrix *) inMat;
6017 C_csr = hypre_ParMatmul((hypre_ParCSRMatrix *)A_csr,B_csr);
6018 return (void *) C_csr;
6019 }
6020
6021 //***************************************************************************
6022 // this function returns the residual norm
6023 //---------------------------------------------------------------------------
6024
HYPRE_LSC_GetRNorm()6025 double HYPRE_LinSysCore::HYPRE_LSC_GetRNorm()
6026 {
6027 return rnorm_;
6028 }
6029
6030