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 // Date : Apr 26, 2002 (This version works sequentially for up to 10000 elems)
10 //***************************************************************************
11 // system includes
12 //---------------------------------------------------------------------------
13
14 #include <stdlib.h>
15 #include <string.h>
16 #include <stdio.h>
17 #include <math.h>
18
19 #define HYPRE_SLIDEMAX 100
20 #define HYPRE_BITMASK2 3
21
22 //***************************************************************************
23 // local includes
24 //---------------------------------------------------------------------------
25
26 #include "HYPRE.h"
27 #include "HYPRE_SlideReduction.h"
28 #include "HYPRE_LSI_mli.h"
29 #include "parcsr_mv/_hypre_parcsr_mv.h"
30 #include "parcsr_ls/_hypre_parcsr_ls.h"
31 #include "seq_mv/seq_mv.h"
32 #include "HYPRE_FEI.h"
33
34 //***************************************************************************
35 // local defines and external functions
36 //---------------------------------------------------------------------------
37
38 #define habs(x) (((x) > 0.0) ? x : -(x))
39
40 extern "C"
41 {
42 // int hypre_BoomerAMGBuildCoarseOperator(hypre_ParCSRMatrix*,
43 // hypre_ParCSRMatrix*, hypre_ParCSRMatrix*, hypre_ParCSRMatrix**);
44 //void hypre_qsort0(int *, int, int);
45 //void hypre_qsort1(int *, double *, int, int);
46 //int HYPRE_LSI_Search(int*, int, int);
47 //int HYPRE_LSI_qsort1a(int *, int *, int, int);
48 //int HYPRE_LSI_MatrixInverse(double **, int, double ***);
49 }
50
51 //***************************************************************************
52 // Constructor
53 //---------------------------------------------------------------------------
54
HYPRE_SlideReduction(MPI_Comm comm)55 HYPRE_SlideReduction::HYPRE_SlideReduction(MPI_Comm comm)
56 {
57 Amat_ = NULL;
58 A21mat_ = NULL;
59 invA22mat_ = NULL;
60 reducedAmat_ = NULL;
61 reducedBvec_ = NULL;
62 reducedXvec_ = NULL;
63 reducedRvec_ = NULL;
64 mpiComm_ = comm;
65 outputLevel_ = 0;
66 procNConstr_ = NULL;
67 slaveEqnList_ = NULL;
68 slaveEqnListAux_ = NULL;
69 gSlaveEqnList_ = NULL;
70 gSlaveEqnListAux_ = NULL;
71 constrBlkInfo_ = NULL;
72 constrBlkSizes_ = NULL;
73 eqnStatuses_ = NULL;
74 blockMinNorm_ = 1.0e-4;
75 hypreRAP_ = NULL;
76 truncTol_ = 1.0e-20;
77 scaleMatrixFlag_ = 0;
78 ADiagISqrts_ = NULL;
79 useSimpleScheme_ = 0;
80 }
81
82 //***************************************************************************
83 // destructor
84 //---------------------------------------------------------------------------
85
~HYPRE_SlideReduction()86 HYPRE_SlideReduction::~HYPRE_SlideReduction()
87 {
88 Amat_ = NULL;
89 mpiComm_ = 0;
90 if ( procNConstr_ != NULL ) delete [] procNConstr_;
91 if ( slaveEqnList_ != NULL ) delete [] slaveEqnList_;
92 if ( slaveEqnListAux_ != NULL ) delete [] slaveEqnListAux_;
93 if ( eqnStatuses_ != NULL ) delete [] eqnStatuses_;
94 if ( gSlaveEqnList_ != NULL ) delete [] gSlaveEqnList_;
95 if ( gSlaveEqnListAux_ != NULL ) delete [] gSlaveEqnListAux_;
96 if ( constrBlkInfo_ != NULL ) delete [] constrBlkInfo_;
97 if ( constrBlkSizes_ != NULL ) delete [] constrBlkSizes_;
98 if ( A21mat_ != NULL ) HYPRE_IJMatrixDestroy(A21mat_);
99 if ( invA22mat_ != NULL ) HYPRE_IJMatrixDestroy(invA22mat_);
100 if ( reducedAmat_ != NULL ) HYPRE_IJMatrixDestroy(reducedAmat_);
101 if ( reducedBvec_ != NULL ) HYPRE_IJVectorDestroy(reducedBvec_);
102 if ( reducedXvec_ != NULL ) HYPRE_IJVectorDestroy(reducedXvec_);
103 if ( reducedRvec_ != NULL ) HYPRE_IJVectorDestroy(reducedRvec_);
104 if ( hypreRAP_ != NULL ) HYPRE_ParCSRMatrixDestroy(hypreRAP_);
105 if ( ADiagISqrts_ != NULL ) delete [] ADiagISqrts_;
106 }
107
108 //***************************************************************************
109 // set output level
110 //---------------------------------------------------------------------------
111
setOutputLevel(int level)112 int HYPRE_SlideReduction::setOutputLevel( int level )
113 {
114 if ( level == 1 ) outputLevel_ |= 1;
115 if ( level == 2 ) outputLevel_ |= 2;
116 if ( level == 3 ) outputLevel_ |= 4;
117 return 0;
118 }
119
120 //***************************************************************************
121 // set use simple scheme
122 //---------------------------------------------------------------------------
123
setUseSimpleScheme()124 int HYPRE_SlideReduction::setUseSimpleScheme()
125 {
126 useSimpleScheme_ = 1;
127 return 0;
128 }
129
130 //***************************************************************************
131 // set truncation threshold for matrix
132 //---------------------------------------------------------------------------
133
setTruncationThreshold(double trunc)134 int HYPRE_SlideReduction::setTruncationThreshold(double trunc)
135 {
136 truncTol_ = trunc;
137 return 0;
138 }
139
140 //***************************************************************************
141 // enable scaling the matrix
142 //---------------------------------------------------------------------------
143
setScaleMatrix()144 int HYPRE_SlideReduction::setScaleMatrix()
145 {
146 scaleMatrixFlag_ = 1;
147 return 0;
148 }
149
150 //***************************************************************************
151 // set the minimum norm for stable blocks
152 //---------------------------------------------------------------------------
153
setBlockMinNorm(double norm)154 int HYPRE_SlideReduction::setBlockMinNorm(double norm)
155 {
156 blockMinNorm_ = norm;
157 return 0;
158 }
159
160 //***************************************************************************
161 // get matrix number of rows
162 //---------------------------------------------------------------------------
163
getMatrixNumRows()164 int HYPRE_SlideReduction::getMatrixNumRows()
165 {
166 int mypid, nprocs, *procNRows, localNRows, nConstraints;
167 HYPRE_ParCSRMatrix A_csr;
168
169 MPI_Comm_rank( mpiComm_, &mypid );
170 MPI_Comm_size( mpiComm_, &nprocs );
171 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
172 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
173 localNRows = procNRows[mypid+1] - procNRows[mypid];
174 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
175 hypre_TFree( procNRows,HYPRE_MEMORY_HOST );
176 return (localNRows-nConstraints);
177 }
178
179 //***************************************************************************
180 // get matrix diagonal
181 //---------------------------------------------------------------------------
182
getMatrixDiagonal()183 double *HYPRE_SlideReduction::getMatrixDiagonal()
184 {
185 return ADiagISqrts_;
186 }
187
188 //***************************************************************************
189 // get reduced matrix
190 //---------------------------------------------------------------------------
191
getReducedMatrix(HYPRE_IJMatrix * mat)192 int HYPRE_SlideReduction::getReducedMatrix(HYPRE_IJMatrix *mat)
193 {
194 (*mat) = reducedAmat_;
195 return 0;
196 }
197
198 //***************************************************************************
199 // get reduced rhs
200 //---------------------------------------------------------------------------
201
getReducedRHSVector(HYPRE_IJVector * rhs)202 int HYPRE_SlideReduction::getReducedRHSVector(HYPRE_IJVector *rhs)
203 {
204 (*rhs) = reducedBvec_;
205 return 0;
206 }
207
208 //***************************************************************************
209 // get reduced solution vector
210 //---------------------------------------------------------------------------
211
getReducedSolnVector(HYPRE_IJVector * sol)212 int HYPRE_SlideReduction::getReducedSolnVector(HYPRE_IJVector *sol)
213 {
214 (*sol) = reducedXvec_;
215 return 0;
216 }
217
218 //***************************************************************************
219 // get auxiliary (temporary) vector
220 //---------------------------------------------------------------------------
221
getReducedAuxVector(HYPRE_IJVector * auxV)222 int HYPRE_SlideReduction::getReducedAuxVector(HYPRE_IJVector *auxV )
223 {
224 (*auxV) = reducedRvec_;
225 return 0;
226 }
227
228 //***************************************************************************
229 // get processor to constraint map
230 //---------------------------------------------------------------------------
231
getProcConstraintMap(int ** map)232 int HYPRE_SlideReduction::getProcConstraintMap(int **map)
233 {
234 (*map) = procNConstr_;
235 return 0;
236 }
237
238 //***************************************************************************
239 // get slave equation list
240 //---------------------------------------------------------------------------
241
getSlaveEqnList(int ** slist)242 int HYPRE_SlideReduction::getSlaveEqnList(int **slist)
243 {
244 (*slist) = slaveEqnList_;
245 return 0;
246 }
247
248 //***************************************************************************
249 // get perturbation matrix (reduced = sub(A) - perturb(A))
250 // (for oorrecting the null space)
251 //---------------------------------------------------------------------------
252
getPerturbationMatrix(HYPRE_ParCSRMatrix * matrix)253 int HYPRE_SlideReduction::getPerturbationMatrix(HYPRE_ParCSRMatrix *matrix)
254 {
255 (*matrix) = hypreRAP_;
256 hypreRAP_ = NULL;
257 return 0;
258 }
259
260 //***************************************************************************
261 // Given the matrix (A) within the object, compute the reduced system and put
262 // it in place. Additional information given are :
263 //
264 // Additional assumptions are :
265 //
266 // - a given slave equation and the corresponding constraint equation
267 // reside in the same processor
268 // - constraint equations are given at the end of the local matrix
269 // (hence given by endRow-nConstr to endRow)
270 // - each processor gets a contiguous block of equations, and processor
271 // i+1 has equation numbers higher than those of processor i
272 //---------------------------------------------------------------------------
273
setup(HYPRE_IJMatrix A,HYPRE_IJVector x,HYPRE_IJVector b)274 int HYPRE_SlideReduction::setup(HYPRE_IJMatrix A, HYPRE_IJVector x,
275 HYPRE_IJVector b)
276 {
277 int mypid, nprocs, ierr, maxBSize=HYPRE_SLIDEMAX, bSize=2;
278 int *procNRows, nrows1, nrows2, reduceAFlag;
279 HYPRE_ParCSRMatrix A_csr;
280 HYPRE_ParVector b_csr;
281
282 //------------------------------------------------------------------
283 // initial set up
284 //------------------------------------------------------------------
285
286 MPI_Comm_rank( mpiComm_, &mypid );
287 MPI_Comm_size( mpiComm_, &nprocs );
288 if ( mypid == 0 && (outputLevel_ & HYPRE_BITMASK2) >= 1 )
289 printf("%4d : HYPRE_SlideReduction begins....\n", mypid);
290
291
292 //------------------------------------------------------------------
293 // check matrix and vector compatibility
294 //------------------------------------------------------------------
295
296 reduceAFlag = 1;
297 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
298 HYPRE_ParCSRMatrixGetRowPartitioning(A_csr, &procNRows);
299 nrows1 = procNRows[nprocs] - procNRows[0];
300 free(procNRows);
301 HYPRE_IJMatrixGetObject(A, (void **) &A_csr);
302 HYPRE_ParCSRMatrixGetRowPartitioning(A_csr, &procNRows);
303 nrows2 = procNRows[nprocs] - procNRows[0];
304 free(procNRows);
305 if (nrows1 != nrows2) reduceAFlag = 0;
306 if (reduceAFlag == 0)
307 {
308 HYPRE_IJVectorGetObject(b, (void **) &b_csr);
309 procNRows = hypre_ParVectorPartitioning((hypre_ParVector *) b_csr);
310 nrows2 = procNRows[nprocs] - procNRows[0];
311 if (nrows1 != nrows2)
312 {
313 if (mypid == 0)
314 printf("HYPRE_SlideReduction ERROR - A,b dim mismatch (reuse)!\n");
315 exit(1);
316 }
317 }
318
319 //------------------------------------------------------------------
320 // clean up first
321 //------------------------------------------------------------------
322
323 if (reduceAFlag == 1)
324 {
325 Amat_ = A;
326 if ( procNConstr_ != NULL ) delete [] procNConstr_;
327 if ( slaveEqnList_ != NULL ) delete [] slaveEqnList_;
328 if ( slaveEqnListAux_ != NULL ) delete [] slaveEqnListAux_;
329 if ( gSlaveEqnList_ != NULL ) delete [] gSlaveEqnList_;
330 if ( gSlaveEqnListAux_ != NULL ) delete [] gSlaveEqnListAux_;
331 if ( constrBlkInfo_ != NULL ) delete [] constrBlkInfo_;
332 if ( constrBlkSizes_ != NULL ) delete [] constrBlkSizes_;
333 if ( eqnStatuses_ != NULL ) delete [] eqnStatuses_;
334 if ( invA22mat_ != NULL ) HYPRE_IJMatrixDestroy(invA22mat_);
335 if ( A21mat_ != NULL ) HYPRE_IJMatrixDestroy(A21mat_);
336 if ( reducedAmat_ != NULL ) HYPRE_IJMatrixDestroy(reducedAmat_);
337 if ( reducedBvec_ != NULL ) HYPRE_IJVectorDestroy(reducedBvec_);
338 if ( reducedXvec_ != NULL ) HYPRE_IJVectorDestroy(reducedXvec_);
339 if ( reducedRvec_ != NULL ) HYPRE_IJVectorDestroy(reducedRvec_);
340 procNConstr_ = NULL;
341 slaveEqnList_ = NULL;
342 slaveEqnListAux_ = NULL;
343 gSlaveEqnList_ = NULL;
344 gSlaveEqnListAux_ = NULL;
345 eqnStatuses_ = NULL;
346 constrBlkInfo_ = NULL;
347 constrBlkSizes_ = NULL;
348 reducedAmat_ = NULL;
349 invA22mat_ = NULL;
350 A21mat_ = NULL;
351 reducedBvec_ = NULL;
352 reducedXvec_ = NULL;
353 reducedRvec_ = NULL;
354 }
355 else
356 {
357 if ( reducedBvec_ != NULL ) HYPRE_IJVectorDestroy(reducedBvec_);
358 if ( reducedXvec_ != NULL ) HYPRE_IJVectorDestroy(reducedXvec_);
359 if ( reducedRvec_ != NULL ) HYPRE_IJVectorDestroy(reducedRvec_);
360 reducedBvec_ = NULL;
361 reducedXvec_ = NULL;
362 reducedRvec_ = NULL;
363 }
364
365 //------------------------------------------------------------------
366 // find the number of constraints in the local processor
367 //------------------------------------------------------------------
368
369 if (reduceAFlag == 1)
370 {
371 if ( findConstraints() == 0 ) return 0;
372 }
373
374 //------------------------------------------------------------------
375 // see if we can find a set of slave nodes for the constraint eqns
376 // If not, search for block size of 2 or higher.
377 //------------------------------------------------------------------
378
379 if (reduceAFlag == 1)
380 {
381 if ( useSimpleScheme_ == 0 )
382 {
383 ierr = findSlaveEqns1();
384 while (ierr < 0 && bSize <= maxBSize)
385 ierr = findSlaveEqnsBlock(bSize++);
386 if ( ierr < 0 )
387 {
388 printf("%4d : HYPRE_SlideReduction ERROR - fail !\n", mypid);
389 exit(1);
390 }
391 composeGlobalList();
392 }
393 }
394
395 //------------------------------------------------------------------
396 // build the reduced matrix
397 //------------------------------------------------------------------
398
399 if (reduceAFlag == 1)
400 {
401 if (useSimpleScheme_ == 0) buildReducedMatrix();
402 else buildSubMatrices();
403 }
404
405 //------------------------------------------------------------------
406 // build the reduced right hand side vector
407 //------------------------------------------------------------------
408
409 if (useSimpleScheme_ == 0) buildReducedRHSVector(b);
410 else buildModifiedRHSVector(x,b);
411
412 //------------------------------------------------------------------
413 // if scale matrix is request, scale matrix and vector
414 //------------------------------------------------------------------
415
416 if ( scaleMatrixFlag_ == 1 )
417 {
418 if (reduceAFlag == 1) scaleMatrixVector();
419 else
420 {
421 if (mypid == 0)
422 printf("HYPRE_SlideReduction ERROR - reuse & scale don't match!\n");
423 exit(1);
424 }
425 }
426
427 //------------------------------------------------------------------
428 // clean up and return
429 //------------------------------------------------------------------
430
431 if ( mypid == 0 && ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
432 printf("%4d : HYPRE_SlideReduction ends.\n", mypid);
433 return 0;
434 }
435
436 //***************************************************************************
437 // search for local constraints (end of the matrix block)
438 //---------------------------------------------------------------------------
439
findConstraints()440 int HYPRE_SlideReduction::findConstraints()
441 {
442 int mypid, nprocs, *procNRows, startRow, endRow;
443 int nConstraints, irow, ncnt, isAConstr, jcol, rowSize, *colInd;
444 int *iTempList, ip, globalNConstr;
445 double *colVal;
446 HYPRE_ParCSRMatrix A_csr;
447
448 //------------------------------------------------------------------
449 // get matrix information
450 //------------------------------------------------------------------
451
452 MPI_Comm_rank(mpiComm_, &mypid);
453 MPI_Comm_size(mpiComm_, &nprocs);
454 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
455 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
456 startRow = procNRows[mypid];
457 endRow = procNRows[mypid+1] - 1;
458 free( procNRows );
459
460 //------------------------------------------------------------------
461 // search for number of local constraints
462 // (==> nConstraints)
463 //------------------------------------------------------------------
464
465 //#define PRINTC
466 #ifdef PRINTC
467 int localNRows = endRow - startRow + 1;
468 char filename[100];
469 FILE *fp;
470 sprintf( filename, "Constr.%d", localNRows);
471 fp = fopen( filename, "w" );
472 #endif
473 nConstraints = 0;
474 for ( irow = endRow; irow >= startRow; irow-- )
475 {
476 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
477 isAConstr = 1;
478 for ( jcol = 0; jcol < rowSize; jcol++ )
479 {
480 if ( colInd[jcol] == irow && colVal[jcol] != 0.0 )
481 {
482 isAConstr = 0;
483 break;
484 }
485 }
486 #ifdef PRINTC
487 if ( isAConstr )
488 {
489 for ( jcol = 0; jcol < rowSize; jcol++ )
490 fprintf(fp,"%8d %8d %e\n",nConstraints+1,colInd[jcol]+1,
491 colVal[jcol]);
492 }
493 #endif
494 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
495 if ( isAConstr ) nConstraints++;
496 else break;
497 }
498 #ifdef PRINTC
499 fclose(fp);
500 #endif
501 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
502 printf("%4d : findConstraints - number of constraints = %d\n",
503 mypid, nConstraints);
504
505 //------------------------------------------------------------------
506 // compute the base nConstraints on each processor
507 // (==> globalNConstr, procNConstr)
508 //------------------------------------------------------------------
509
510 iTempList = new int[nprocs];
511 if ( procNConstr_ != NULL ) delete [] procNConstr_;
512 procNConstr_ = new int[nprocs+1];
513 for ( ip = 0; ip < nprocs; ip++ ) iTempList[ip] = 0;
514 iTempList[mypid] = nConstraints;
515 MPI_Allreduce(iTempList,procNConstr_,nprocs,MPI_INT,MPI_SUM,mpiComm_);
516 delete [] iTempList;
517 globalNConstr = 0;
518 ncnt = 0;
519 for ( ip = 0; ip < nprocs; ip++ )
520 {
521 ncnt = procNConstr_[ip];
522 procNConstr_[ip] = globalNConstr;
523 globalNConstr += ncnt;
524 }
525 procNConstr_[nprocs] = globalNConstr;
526 if ( slaveEqnList_ != NULL ) delete [] slaveEqnList_;
527 if ( nConstraints > 0 ) slaveEqnList_ = new int[nConstraints];
528 else slaveEqnList_ = NULL;
529 for ( irow = 0; irow < nConstraints; irow++ ) slaveEqnList_[irow] = -1;
530 if ( constrBlkInfo_ != NULL ) delete [] constrBlkInfo_;
531 if ( nConstraints > 0 ) constrBlkInfo_ = new int[nConstraints];
532 else constrBlkInfo_ = NULL;
533 for ( irow = 0; irow < nConstraints; irow++ ) constrBlkInfo_[irow] = -1;
534 if ( constrBlkSizes_ != NULL ) delete [] constrBlkSizes_;
535 if ( nConstraints > 0 ) constrBlkSizes_ = new int[nConstraints];
536 else constrBlkSizes_ = NULL;
537 for ( irow = 0; irow < nConstraints; irow++ ) constrBlkSizes_[irow] = 0;
538 if ( nConstraints > 0 )
539 {
540 eqnStatuses_ = new int[endRow-nConstraints-startRow+1];
541 for (irow = 0; irow < endRow-nConstraints-startRow+1; irow++ )
542 eqnStatuses_[irow] = 0;
543 }
544 else eqnStatuses_ = NULL;
545 return globalNConstr;
546 }
547
548 //***************************************************************************
549 // search for a slave equation list (block size = 1)
550 //---------------------------------------------------------------------------
551
findSlaveEqns1()552 int HYPRE_SlideReduction::findSlaveEqns1()
553 {
554 int mypid, nprocs, *procNRows, startRow, endRow;
555 int nConstraints, irow, jcol, rowSize, ncnt, *colInd, index;
556 int nCandidates, *candidateList;
557 int *constrListAux, colIndex, searchIndex, procIndex, uBound;
558 int nSum, newEndRow;
559 double *colVal, searchValue;
560 HYPRE_ParCSRMatrix A_csr;
561
562 //------------------------------------------------------------------
563 // get matrix information
564 //------------------------------------------------------------------
565
566 MPI_Comm_rank(mpiComm_, &mypid);
567 MPI_Comm_size(mpiComm_, &nprocs);
568 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
569 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
570 startRow = procNRows[mypid];
571 endRow = procNRows[mypid+1] - 1;
572 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
573 newEndRow = endRow - nConstraints;
574
575 //------------------------------------------------------------------
576 // compose candidate slave list (slaves in candidateList, corresponding
577 // constraint equation in constrListAux)
578 //------------------------------------------------------------------
579
580 nCandidates = 0;
581 candidateList = NULL;
582 constrListAux = NULL;
583 if ( nConstraints > 0 )
584 {
585 candidateList = new int[newEndRow-startRow+1];
586 constrListAux = new int[newEndRow-startRow+1];
587
588 //------------------------------------------------------------------
589 // candidates are those with 1 link to the constraint list
590 //------------------------------------------------------------------
591
592 for ( irow = startRow; irow <= endRow-nConstraints; irow++ )
593 {
594 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
595 ncnt = 0;
596 constrListAux[irow-startRow] = -1;
597 for ( jcol = 0; jcol < rowSize; jcol++ )
598 {
599 colIndex = colInd[jcol];
600 for ( procIndex = 1; procIndex <= nprocs; procIndex++ )
601 if ( colIndex < procNRows[procIndex] ) break;
602 uBound = procNRows[procIndex] - (procNConstr_[procIndex] -
603 procNConstr_[procIndex-1]);
604 if ( colIndex >= uBound && procIndex == (mypid+1) )
605 {
606 ncnt++;
607 searchIndex = colIndex;
608 }
609 else if ( colIndex >= uBound && procIndex != (mypid+1) ) ncnt = 2;
610 if ( ncnt > 1 ) break;
611 }
612 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
613 if (ncnt == 1 && searchIndex > newEndRow && searchIndex <= endRow)
614 {
615 constrListAux[nCandidates] = searchIndex;
616 candidateList[nCandidates++] = irow;
617 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 3 )
618 printf("%4d : findSlaveEqns1 - candidate %d = %d(%d)\n",
619 mypid, nCandidates-1, irow, searchIndex);
620 }
621 }
622 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
623 printf("%4d : findSlaveEqns1 - nCandidates, nConstr = %d %d\n",
624 mypid, nCandidates, nConstraints);
625 }
626
627 //---------------------------------------------------------------------
628 // search the constraint equations for the selected slave equations
629 // (search for candidates column index with maximum magnitude)
630 // ==> slaveEqnList_
631 //---------------------------------------------------------------------
632
633 searchIndex = 0;
634 for ( irow = endRow-nConstraints+1; irow <= endRow; irow++ )
635 {
636 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
637 searchIndex = -1;
638 searchValue = 1.0E-6;
639 for ( jcol = 0; jcol < rowSize; jcol++ )
640 {
641 if (colVal[jcol] != 0.0 && colInd[jcol] >= startRow &&
642 colInd[jcol] <= (endRow-nConstraints) &&
643 eqnStatuses_[colInd[jcol]-startRow] == 0)
644 {
645 colIndex = hypre_BinarySearch(candidateList, colInd[jcol],
646 nCandidates);
647 if ( colIndex >= 0 && habs(colVal[jcol]) > searchValue )
648 {
649 if ( irow != constrListAux[colIndex] ) break;
650 searchValue = habs(colVal[jcol]);
651 searchIndex = colInd[jcol];
652 }
653 }
654 }
655 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
656 if ( searchIndex >= 0 )
657 {
658 index = irow - endRow + nConstraints - 1;
659 slaveEqnList_[index] = searchIndex;
660 constrBlkInfo_[index] = index;
661 constrBlkSizes_[index] = 1;
662 eqnStatuses_[searchIndex-startRow] = 1;
663 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 2 )
664 printf("%4d : findSlaveEqns1 - constr %7d <=> slave %d\n",
665 mypid, irow, searchIndex);
666 }
667 else
668 {
669 slaveEqnList_[irow-endRow+nConstraints-1] = -1;
670 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 2 )
671 {
672 printf("%4d : findSlaveEqns1 - constraint %4d fails",mypid,irow);
673 printf(" to find a slave.\n");
674 }
675 }
676 }
677 if ( nConstraints > 0 )
678 {
679 delete [] constrListAux;
680 delete [] candidateList;
681 }
682 free( procNRows );
683
684 //---------------------------------------------------------------------
685 // if not all constraint-slave pairs can be found, return -1
686 //---------------------------------------------------------------------
687
688 ncnt = 0;
689 for ( irow = 0; irow < nConstraints; irow++ )
690 if ( slaveEqnList_[irow] == -1 ) ncnt++;
691 MPI_Allreduce(&ncnt, &nSum, 1, MPI_INT, MPI_SUM, mpiComm_);
692 if ( nSum > 0 )
693 {
694 if ( mypid == 0 && ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
695 {
696 printf("%4d : findSlaveEqns1 fails - total number of unsatisfied",
697 mypid);
698 printf(" constraints = %d \n", nSum);
699 }
700 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
701 {
702 for ( irow = 0; irow < nConstraints; irow++ )
703 if ( slaveEqnList_[irow] == -1 )
704 {
705 printf("%4d : findSlaveEqns1 - unsatisfied constraint",mypid);
706 printf(" equation = %d\n", irow+endRow-nConstraints+1);
707 }
708 }
709 return -1;
710 } else return 0;
711 }
712
713 //***************************************************************************
714 // search for a slave equation list (block size of up to blkSizeMax)
715 //---------------------------------------------------------------------------
716
findSlaveEqnsBlock(int blkSize)717 int HYPRE_SlideReduction::findSlaveEqnsBlock(int blkSize)
718 {
719 int mypid, nprocs, *procNRows, startRow, endRow, localNRows;
720 int nConstraints, irow, jcol, rowSize, ncnt, *colInd;
721 int nCandidates, *candidateList, **constrListAuxs, colIndex, rowSize2;
722 int ic, ii, jj, searchIndex, searchInd2, newEndRow;
723 int blkSizeMax=HYPRE_SLIDEMAX, *colInd2, searchInd3;
724 int constrIndex, uBound, lBound, nSum, isACandidate, newBlkSize, *colTmp;
725 int constrIndex2, searchBlkSize, newIndex, oldIndex, irowLocal, ip;
726 int *blkInfo, blkInfoCnt;
727 double *colVal, searchValue, retVal, *colVal2;
728 HYPRE_ParCSRMatrix A_csr;
729
730 //------------------------------------------------------------------
731 // if block size too large, return error
732 //------------------------------------------------------------------
733
734 if ( blkSize > blkSizeMax ) return -1;
735
736 //------------------------------------------------------------------
737 // get matrix information
738 //------------------------------------------------------------------
739
740 MPI_Comm_rank(mpiComm_, &mypid);
741 MPI_Comm_size(mpiComm_, &nprocs);
742 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
743 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
744 startRow = procNRows[mypid];
745 endRow = procNRows[mypid+1] - 1;
746 localNRows = endRow - startRow + 1;
747 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
748 newEndRow = endRow - nConstraints;
749 if ( mypid == 0 && ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
750 printf("%4d : findSlaveEqnsBlock - size = %d\n", mypid, blkSize);
751
752 //------------------------------------------------------------------
753 // compose candidate slave list (slaves in candidateList, corresponding
754 // constraint equation in constrListAuxs)
755 //------------------------------------------------------------------
756
757 nCandidates = 0;
758 if ( nConstraints > 0 )
759 {
760 candidateList = new int[localNRows-nConstraints];
761 constrListAuxs = new int*[localNRows-nConstraints];
762 for ( ic = 0; ic < localNRows-nConstraints; ic++ )
763 {
764 constrListAuxs[ic] = new int[blkSize];
765 for (jcol = 0; jcol < blkSize; jcol++) constrListAuxs[ic][jcol] = -1;
766 }
767
768 //---------------------------------------------------------------
769 // candidates are those with <input> links to the constraint list
770 //---------------------------------------------------------------
771
772 for ( irow = startRow; irow <= endRow-nConstraints; irow++ )
773 {
774 if ( eqnStatuses_[irow-startRow] == 1 ) continue;
775 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
776 ncnt = 0;
777 for ( jcol = 0; jcol < rowSize; jcol++ )
778 {
779 colIndex = colInd[jcol];
780 for ( ip = 0; ip < nprocs; ip++ )
781 {
782 uBound = procNRows[ip+1];
783 lBound = uBound - (procNConstr_[ip+1] - procNConstr_[ip]);
784 if ( colIndex >= lBound && colIndex < uBound && ip == mypid )
785 {
786 ncnt++;
787 if ( ncnt <= blkSize )
788 constrListAuxs[nCandidates][ncnt-1] = colIndex;
789 }
790 else if (colIndex >= lBound && colIndex < uBound && ip != mypid)
791 {
792 ncnt = blkSize + 1;
793 break;
794 }
795 }
796 if ( ncnt > blkSize ) break;
797 }
798 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
799 if ( ncnt >= 1 && ncnt <= blkSize )
800 {
801 isACandidate = 1;
802 for ( ic = 0; ic < ncnt; ic++ )
803 {
804 if ( constrListAuxs[nCandidates][ic] <= newEndRow ||
805 constrListAuxs[nCandidates][ic] > endRow )
806 {
807 isACandidate = 0;
808 break;
809 }
810 }
811 if ( isACandidate )
812 {
813 candidateList[nCandidates++] = irow;
814 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 3 )
815 printf("%4d : findSlaveEqnsBlock - candidate %d = %d\n",
816 mypid, nCandidates-1, irow);
817 }
818 }
819 }
820 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
821 printf("%4d : findSlaveEqnsBlock - nCandidates, nConstr = %d %d\n",
822 mypid, nCandidates, nConstraints);
823 }
824
825 //---------------------------------------------------------------------
826 // revise candidates taking into consideration nonsymmetry
827 //---------------------------------------------------------------------
828
829 int *tempSlaveList, *tempSlaveListAux;
830 if ( nConstraints > 0 ) tempSlaveList = new int[nConstraints];
831 if ( nConstraints > 0 ) tempSlaveListAux = new int[nConstraints];
832 for (irow = 0; irow < nConstraints; irow++)
833 {
834 tempSlaveList[irow] = slaveEqnList_[irow];
835 tempSlaveListAux[irow] = irow;
836 }
837 HYPRE_LSI_qsort1a(tempSlaveList, tempSlaveListAux, 0, nConstraints-1);
838
839 /* for each of the candidates, examine all associated constraints dof */
840
841 for ( irow = 0; irow < nCandidates; irow++ )
842 {
843 for ( ic = 0; ic < blkSize; ic++ )
844 {
845 constrIndex = constrListAuxs[irow][ic];
846 /* if valid constraint number */
847 if ( constrIndex >= 0 )
848 {
849 /* get the constraint row */
850 HYPRE_ParCSRMatrixGetRow(A_csr,constrIndex,&rowSize,&colInd,NULL);
851
852 /* for each nonzero entry of the constraint row */
853 /* - see if the column number is an already selected slave */
854 /* - if so, find the corresponding constraint no. of that slave */
855 /* - add that constraint to my list */
856 for ( jcol = 0; jcol < rowSize; jcol++ )
857 {
858 colIndex = colInd[jcol];
859 searchIndex = hypre_BinarySearch(tempSlaveList,colIndex,
860 nConstraints);
861 if ( searchIndex >= 0 )
862 {
863 searchInd2 = tempSlaveListAux[searchIndex] + newEndRow + 1;
864 for ( ip = 0; ip < blkSize; ip++ )
865 if ( constrListAuxs[irow][ip] == searchInd2 ||
866 constrListAuxs[irow][ip] == -1 ) break;
867 if ( ip == blkSize ) constrListAuxs[irow][0] = -5;
868 else if ( constrListAuxs[irow][ip] == -1 )
869 {
870 constrListAuxs[irow][ip] = searchInd2;
871 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 2 )
872 printf("Slave candidate %d adds new constr %d\n",
873 candidateList[irow], searchInd2);
874 }
875 }
876 }
877 HYPRE_ParCSRMatrixRestoreRow(A_csr,constrIndex,&rowSize,&colInd,
878 NULL);
879 }
880 }
881 }
882
883 /* delete candidates that gives larger than expected blocksize */
884
885 ncnt = 0;
886 for ( irow = 0; irow < nCandidates; irow++ )
887 {
888 if ( constrListAuxs[irow][0] != -5 )
889 {
890 if ( irow != ncnt )
891 {
892 if ( constrListAuxs[ncnt] != NULL ) delete [] constrListAuxs[ncnt];
893 constrListAuxs[ncnt] = constrListAuxs[irow];
894 constrListAuxs[irow] = NULL;
895 candidateList[ncnt++] = candidateList[irow];
896 }
897 else ncnt++;
898 }
899 }
900 nCandidates = ncnt;
901 if ( nConstraints > 0 ) delete [] tempSlaveList;
902 if ( nConstraints > 0 ) delete [] tempSlaveListAux;
903
904 //---------------------------------------------------------------------
905 // search the constraint equations for the selected slave equations
906 // (search for candidates column index with maximum magnitude)
907 // ==> slaveEqnList_
908 //---------------------------------------------------------------------
909
910 searchIndex = 0;
911
912 blkInfo = new int[blkSize+HYPRE_SLIDEMAX];
913 for ( irow = newEndRow+1; irow <= endRow; irow++ )
914 {
915 /* -- if slave variable has not been picked for constraint irow -- */
916
917 irowLocal = irow - endRow + nConstraints - 1;
918
919 if ( slaveEqnList_[irowLocal] == -1 )
920 {
921 /* -- get the constraint row, and search for nonzero entries -- */
922
923 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize2,&colInd2,&colVal2);
924 rowSize = rowSize2;
925 colInd = new int[rowSize];
926 colVal = new double[rowSize];
927 for ( jcol = 0; jcol < rowSize; jcol++ )
928 {
929 colInd[jcol] = colInd2[jcol];
930 colVal[jcol] = colVal2[jcol];
931 }
932 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize2,&colInd2,&colVal2);
933 searchIndex = -1;
934 searchValue = blockMinNorm_;
935 for ( jcol = 0; jcol < rowSize; jcol++ )
936 {
937 colIndex = colInd[jcol];
938 if (colVal[jcol] != 0.0 && colIndex >= startRow
939 && colIndex <= newEndRow)
940 {
941 /* -- see if the nonzero entry is a potential candidate -- */
942
943 searchInd2 = hypre_BinarySearch(candidateList, colIndex,
944 nCandidates);
945
946 /* -- if the nonzero entry is a potential candidate, see -- */
947 /* -- if the use of this as slave will give block size -- */
948 /* -- that is too large. -- */
949
950 if (searchInd2 >= 0 && eqnStatuses_[colIndex-startRow] != 1)
951 {
952 newBlkSize = 1;
953 blkInfoCnt = 0;
954 for ( ic = 0; ic < blkSize; ic++ )
955 {
956 constrIndex = constrListAuxs[searchInd2][ic];
957 if ( constrIndex != -1 )
958 {
959 constrIndex2 = constrIndex - endRow + nConstraints - 1;
960 if ( constrIndex != irow &&
961 slaveEqnList_[constrIndex2] != -1)
962 {
963 for ( ip = 0; ip < blkInfoCnt; ip++ )
964 if ( blkInfo[ip] == constrBlkInfo_[constrIndex2] )
965 break;
966 if ( ip == blkInfoCnt )
967 {
968 newBlkSize += constrBlkSizes_[constrIndex2];
969 blkInfo[blkInfoCnt]=constrBlkInfo_[constrIndex2];
970 blkInfoCnt++;
971 }
972 }
973 /*
974 else if (constrIndex != irow ) newBlkSize++;
975 */
976 }
977 }
978 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 2 )
979 {
980 printf("%4d : constraint %d - candidate %d (%d) ", mypid,
981 irow, searchInd2, candidateList[searchInd2]);
982 printf("gives blksize = %d\n", newBlkSize);
983 }
984 /*
985 if (newBlkSize > 1 && newBlkSize <= blkSize)
986 */
987 if (newBlkSize <= blkSize)
988 {
989 retVal = matrixCondEst(irow,colIndex,blkInfo,blkInfoCnt);
990 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 2 )
991 printf("%4d : pivot = %e (%e) : %d\n", mypid, retVal,
992 searchValue,newBlkSize);
993 if ( retVal > searchValue )
994 {
995 searchValue = habs(colVal[jcol]);
996 searchIndex = colIndex;
997 searchBlkSize = newBlkSize;
998 }
999 }
1000 }
1001 }
1002 }
1003 delete [] colInd;
1004 delete [] colVal;
1005 if ( searchIndex >= 0 && searchValue > blockMinNorm_ )
1006 {
1007 searchInd2 = hypre_BinarySearch(candidateList,searchIndex,
1008 nCandidates);
1009 newIndex = -9;
1010 for ( ic = 0; ic < blkSize; ic++ )
1011 {
1012 constrIndex = constrListAuxs[searchInd2][ic];
1013 if ( constrIndex != -1 )
1014 {
1015 constrIndex2 = constrIndex - endRow + nConstraints - 1;
1016 if (constrIndex != irow && slaveEqnList_[constrIndex2] != -1)
1017 {
1018 if (newIndex == -9) newIndex=constrBlkInfo_[constrIndex2];
1019 oldIndex = constrBlkInfo_[constrIndex2];
1020 for ( ii = 0; ii < nConstraints; ii++ )
1021 {
1022 if ( constrBlkInfo_[ii] == oldIndex )
1023 {
1024 constrBlkInfo_[ii] = newIndex;
1025 constrBlkSizes_[ii] = searchBlkSize;
1026 }
1027 }
1028 }
1029 }
1030 }
1031 if (newIndex == -9) newIndex = irowLocal;
1032 constrBlkInfo_[irowLocal] = newIndex;
1033 constrBlkSizes_[irowLocal] = searchBlkSize;
1034 slaveEqnList_[irowLocal] = searchIndex;
1035 searchInd2 = hypre_BinarySearch(candidateList, searchIndex,
1036 nCandidates);
1037 eqnStatuses_[searchIndex-startRow] = 1;
1038
1039 /* update the constrListAux - first get selected slave row */
1040
1041 for ( ii = 0; ii < blkSize; ii++ )
1042 {
1043 constrIndex2 = constrListAuxs[searchInd2][ii];
1044 if ( constrIndex2 != -1 )
1045 {
1046 HYPRE_ParCSRMatrixGetRow(A_csr,constrIndex2,&rowSize2,
1047 &colInd2,&colVal2);
1048 for ( jj = 0; jj < rowSize2; jj++ )
1049 {
1050 searchInd3 = hypre_BinarySearch(candidateList,
1051 colInd2[jj],nCandidates);
1052 if ( searchInd3 >= 0 )
1053 {
1054 for ( ip = 0; ip < blkSize; ip++ )
1055 {
1056 if ( constrListAuxs[searchInd3][ip] == irow ||
1057 constrListAuxs[searchInd3][ip] == -1 ) break;
1058 }
1059 if ( ip == blkSize )
1060 {
1061 constrListAuxs[searchInd3][0] = -5;
1062 eqnStatuses_[colInd2[jj]-startRow] = 1;
1063 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 3 )
1064 printf("*Slave candidate %d disabled.\n",
1065 candidateList[searchInd3]);
1066 }
1067 else if ( constrListAuxs[searchInd3][ip] == -1 )
1068 {
1069 constrListAuxs[searchInd3][ip] = irow;
1070 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 3 )
1071 printf("*Slave candidate %d adds new constr %d\n",
1072 candidateList[searchInd3], irow);
1073 }
1074 }
1075 }
1076 HYPRE_ParCSRMatrixRestoreRow(A_csr,constrIndex2,&rowSize2,
1077 &colInd2,&colVal2);
1078 }
1079 }
1080 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 2 )
1081 printf("%4d : findSlaveEqnsBlock - constr %d <=> slave %d (%d)\n",
1082 mypid, irow, searchIndex, newIndex);
1083 }
1084 else
1085 {
1086 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 2 )
1087 {
1088 if ( searchIndex < 0 && searchValue > blockMinNorm_ )
1089 {
1090 printf("%4d : findSlaveEqnsBlock - constraint %4d fails (0)",
1091 mypid, irow);
1092 printf(" to find a slave.\n");
1093 }
1094 else if ( searchIndex >= 0 && searchValue <= blockMinNorm_ )
1095 {
1096 printf("%4d : findSlaveEqnsBlock - constraint %4d fails (1)",
1097 mypid, irow);
1098 printf(" to find a slave.\n");
1099 }
1100 else
1101 {
1102 printf("%4d : findSlaveEqnsBlock - constraint %4d fails (2)",
1103 mypid, irow);
1104 printf(" to find a slave.\n");
1105 }
1106 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 3 )
1107 {
1108 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
1109 colTmp = new int[rowSize];
1110 rowSize2 = rowSize;
1111 for ( ii = 0; ii < rowSize; ii++ ) colTmp[ii] = colInd[ii];
1112 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,
1113 &colVal);
1114 for ( jcol = 0; jcol < rowSize2; jcol++ )
1115 {
1116 colIndex = colTmp[jcol];
1117 printf("%4d : row %d has col %d (%d,%d) (%d,%d)\n",mypid,
1118 irow,colIndex,jcol,rowSize,procNRows[mypid],
1119 procNRows[mypid+1]);
1120 if ( colIndex >= procNRows[mypid] &&
1121 colIndex < procNRows[mypid+1])
1122 {
1123 HYPRE_ParCSRMatrixGetRow(A_csr,colIndex,&rowSize,
1124 &colInd,NULL);
1125 for ( ii = 0; ii < rowSize; ii++ )
1126 printf("%4d : col %d has col %d (%d,%d)\n",mypid,
1127 colIndex,colInd[ii],ii,rowSize);
1128 HYPRE_ParCSRMatrixRestoreRow(A_csr,colIndex,&rowSize,
1129 &colInd,NULL);
1130 }
1131 }
1132 delete [] colTmp;
1133 }
1134 }
1135 }
1136 }
1137 }
1138 delete [] blkInfo;
1139 if ( nConstraints > 0 )
1140 {
1141 for ( ic = 0; ic < localNRows-nConstraints; ic++ )
1142 if ( constrListAuxs[ic] != NULL ) delete [] constrListAuxs[ic];
1143 delete [] constrListAuxs;
1144 delete [] candidateList;
1145 }
1146 free( procNRows );
1147
1148 #if 0
1149 int is, *iArray1, *iArray2;
1150 if ( constrBlkInfo_ != NULL )
1151 {
1152 iArray1 = new int[nConstraints];
1153 iArray2 = new int[nConstraints];
1154 for ( is = 0; is < nConstraints; is++ )
1155 {
1156 iArray1[is] = constrBlkInfo_[is];
1157 iArray2[is] = constrBlkSizes_[is];
1158 }
1159 HYPRE_LSI_qsort1a(iArray1, iArray2, 0, nConstraints-1);
1160 ip = -1; ncnt = 0;
1161 for ( is = 0; is < nConstraints; is++ )
1162 {
1163 if ( iArray1[is] != ip )
1164 {
1165 iArray1[ncnt] = iArray1[is];
1166 iArray2[ncnt] = iArray2[is];
1167 ncnt++;
1168 ip = iArray1[is];
1169 }
1170 }
1171 HYPRE_LSI_qsort1a(iArray2, iArray1, 0, ncnt-1);
1172 ip = 1;
1173 for ( is = 1; is < ncnt; is++ )
1174 {
1175 if ( iArray2[is] == iArray2[is-1] ) ip++;
1176 else
1177 {
1178 printf("%4d : number of blocks with blksize %6d = %d\n",
1179 mypid, iArray2[is-1], ip);
1180 ip = 1;
1181 }
1182 }
1183 printf("%4d : number of blocks with blksize %6d = %d\n",
1184 mypid, iArray2[ncnt-1], ip);
1185 delete [] iArray1;
1186 delete [] iArray2;
1187 }
1188 #endif
1189
1190 //---------------------------------------------------------------------
1191 // if not all constraint-slave pairs can be found, return -1
1192 //---------------------------------------------------------------------
1193
1194 ncnt = 0;
1195 for ( irow = 0; irow < nConstraints; irow++ )
1196 if ( slaveEqnList_[irow] == -1 ) ncnt++;
1197 MPI_Allreduce(&ncnt, &nSum, 1, MPI_INT, MPI_SUM, mpiComm_);
1198 if ( nSum > 0 )
1199 {
1200 if ( mypid == 0 && ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
1201 {
1202 printf("%4d : findSlaveEqnsBlock fails - total number of unsatisfied",
1203 mypid);
1204 printf(" constraints = %d \n", nSum);
1205 }
1206 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
1207 {
1208 for ( irow = 0; irow < nConstraints; irow++ )
1209 if ( slaveEqnList_[irow] == -1 )
1210 {
1211 printf("%4d : findSlaveEqnsBlock - unsatisfied constraint",mypid);
1212 printf(" equation = %d\n", irow+endRow-nConstraints+1);
1213 }
1214 }
1215 return -1;
1216 }
1217 else return 0;
1218 }
1219
1220 //***************************************************************************
1221 // compose global slave equation list
1222 //---------------------------------------------------------------------------
1223
composeGlobalList()1224 int HYPRE_SlideReduction::composeGlobalList()
1225 {
1226 int mypid, nprocs, nConstraints, is, ip, *recvCntArray, *displArray;
1227 int globalNConstr, ierr, ncnt, *iArray1, *iArray2;
1228
1229 //------------------------------------------------------------------
1230 // fetch machine and constraint parameters
1231 //------------------------------------------------------------------
1232
1233 MPI_Comm_rank( mpiComm_, &mypid );
1234 MPI_Comm_size( mpiComm_, &nprocs );
1235 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
1236 globalNConstr = procNConstr_[nprocs];
1237
1238 //------------------------------------------------------------------
1239 // sort the local selected node list and its auxiliary list, then
1240 // form a global list of slave nodes on each processor
1241 // form the corresponding auxiliary list for later pruning
1242 // ==> slaveEqnListAux_, gSlaveEqnList_, gSlaveEqnListAux_
1243 //------------------------------------------------------------------
1244
1245 if ( slaveEqnListAux_ != NULL ) delete [] slaveEqnListAux_;
1246 if ( gSlaveEqnList_ != NULL ) delete [] gSlaveEqnList_;
1247 if ( gSlaveEqnListAux_ != NULL ) delete [] gSlaveEqnListAux_;
1248 slaveEqnListAux_ = NULL;
1249 if ( nConstraints > 0 )
1250 {
1251 slaveEqnListAux_ = new int[nConstraints];
1252 for ( is = 0; is < nConstraints; is++ ) slaveEqnListAux_[is] = is;
1253 HYPRE_LSI_qsort1a(slaveEqnList_, slaveEqnListAux_, 0, nConstraints-1);
1254 ierr = 0;
1255 for ( is = 1; is < nConstraints; is++ )
1256 {
1257 if ( slaveEqnList_[is] == slaveEqnList_[is-1] )
1258 {
1259 printf("%4d : HYPRE_SlideReduction ERROR - repeated slave",mypid);
1260 printf(" equation %d\n", slaveEqnList_[is]);
1261 ierr = 1;
1262 break;
1263 }
1264 }
1265 if ( ierr )
1266 {
1267 for ( is = 0; is < nConstraints; is++ )
1268 {
1269 printf("%4d : HYPRE_SlideReduction slave %d = %d \n",mypid,is,
1270 slaveEqnList_[is]);
1271 }
1272 exit(1);
1273 }
1274 }
1275 gSlaveEqnList_ = new int[globalNConstr];
1276 gSlaveEqnListAux_ = new int[globalNConstr];
1277
1278 //------------------------------------------------------------------
1279 // compose global slave equation list
1280 //------------------------------------------------------------------
1281
1282 recvCntArray = new int[nprocs];
1283 displArray = new int[nprocs];
1284 MPI_Allgather(&nConstraints,1,MPI_INT,recvCntArray,1,MPI_INT,mpiComm_);
1285 displArray[0] = 0;
1286 for ( ip = 1; ip < nprocs; ip++ )
1287 displArray[ip] = displArray[ip-1] + recvCntArray[ip-1];
1288 for ( ip = 0; ip < nConstraints; ip++ )
1289 slaveEqnListAux_[ip] += displArray[mypid];
1290 MPI_Allgatherv(slaveEqnList_, nConstraints, MPI_INT, gSlaveEqnList_,
1291 recvCntArray, displArray, MPI_INT, mpiComm_);
1292 MPI_Allgatherv(slaveEqnListAux_, nConstraints, MPI_INT, gSlaveEqnListAux_,
1293 recvCntArray, displArray, MPI_INT, mpiComm_);
1294 for ( is = 0; is < nConstraints; is++ )
1295 slaveEqnListAux_[is] -= displArray[mypid];
1296 delete [] recvCntArray;
1297 delete [] displArray;
1298
1299 if ( constrBlkInfo_ != NULL && ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
1300 {
1301 iArray1 = new int[nConstraints];
1302 iArray2 = new int[nConstraints];
1303 for ( is = 0; is < nConstraints; is++ )
1304 {
1305 iArray1[is] = constrBlkInfo_[is];
1306 iArray2[is] = constrBlkSizes_[is];
1307 }
1308 HYPRE_LSI_qsort1a(iArray1, iArray2, 0, nConstraints-1);
1309 ip = -1; ncnt = 0;
1310 for ( is = 0; is < nConstraints; is++ )
1311 {
1312 if ( iArray1[is] != ip )
1313 {
1314 iArray1[ncnt] = iArray1[is];
1315 iArray2[ncnt] = iArray2[is];
1316 ncnt++;
1317 ip = iArray1[is];
1318 }
1319 }
1320 HYPRE_LSI_qsort1a(iArray2, iArray1, 0, ncnt-1);
1321 ip = 1;
1322 for ( is = 1; is < ncnt; is++ )
1323 {
1324 if ( iArray2[is] == iArray2[is-1] ) ip++;
1325 else
1326 {
1327 printf("%4d : number of blocks with blksize %6d = %d\n",
1328 mypid, iArray2[is-1], ip);
1329 ip = 1;
1330 }
1331 }
1332 printf("%4d : number of blocks with blksize %6d = %d\n",
1333 mypid, iArray2[ncnt-1], ip);
1334 delete [] iArray1;
1335 delete [] iArray2;
1336 }
1337
1338 if ( ( outputLevel_ & HYPRE_BITMASK2 ) > 1 )
1339 for ( is = 0; is < nConstraints; is++ )
1340 printf("%4d : HYPRE_SlideReduction - slaveEqnList %d = %d(%d)\n",
1341 mypid, is, slaveEqnList_[is], slaveEqnListAux_[is]);
1342
1343 return 0;
1344 }
1345
1346 //****************************************************************************
1347 // build the submatrix matrix
1348 //----------------------------------------------------------------------------
1349
buildSubMatrices()1350 int HYPRE_SlideReduction::buildSubMatrices()
1351 {
1352 int mypid, nprocs, *procNRows, startRow, endRow, localNRows;
1353 int globalNConstr, globalNRows, nConstraints, ncnt;
1354 int A21NRows, A21NCols, A21GlobalNRows, A21GlobalNCols;
1355 int A21StartRow, A21StartCol, ierr, maxRowSize, *A21MatSize;
1356 int newEndRow, irow, rowIndex, newRowSize, nnzA21, rowCount;
1357 int *colInd, *newColInd, rowSize, *reducedAMatSize, uBound;
1358 int reducedANRows, reducedANCols, reducedAStartRow, reducedAStartCol;
1359 int reducedAGlobalNRows, reducedAGlobalNCols, jcol;
1360 int procIndex, colIndex, totalNNZ, newColIndex;
1361 double *colVal, *newColVal;
1362 HYPRE_ParCSRMatrix A_csr, A21_csr, reducedA_csr;
1363
1364 //------------------------------------------------------------------
1365 // get matrix information
1366 //------------------------------------------------------------------
1367
1368 MPI_Comm_rank(mpiComm_, &mypid);
1369 MPI_Comm_size(mpiComm_, &nprocs);
1370 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
1371 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
1372 startRow = procNRows[mypid];
1373 endRow = procNRows[mypid+1] - 1;
1374 localNRows = endRow - startRow + 1;
1375 globalNConstr = procNConstr_[nprocs];
1376 globalNRows = procNRows[nprocs];
1377 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
1378
1379 //******************************************************************
1380 // extract A21 from A
1381 //------------------------------------------------------------------
1382 // calculate the dimension of A21
1383 //------------------------------------------------------------------
1384
1385 A21NRows = nConstraints;
1386 A21NCols = localNRows - nConstraints;
1387 A21GlobalNRows = globalNConstr;
1388 A21GlobalNCols = globalNRows - globalNConstr;
1389 A21StartRow = procNConstr_[mypid];
1390 A21StartCol = procNRows[mypid] - procNConstr_[mypid];
1391
1392 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
1393 {
1394 printf("%4d : buildA21Mat(2) - A21StartRow = %d\n", mypid, A21StartRow);
1395 printf("%4d : buildA21Mat(2) - A21GlobalDim = %d %d\n", mypid,
1396 A21GlobalNRows, A21GlobalNCols);
1397 printf("%4d : buildA21Mat(2) - A21LocalDim = %d %d\n",mypid,
1398 A21NRows, A21NCols);
1399 }
1400
1401 //------------------------------------------------------------------
1402 // create a matrix context for A21
1403 //------------------------------------------------------------------
1404
1405 ierr = HYPRE_IJMatrixCreate(mpiComm_,A21StartRow,A21StartRow+A21NRows-1,
1406 A21StartCol,A21StartCol+A21NCols-1,&A21mat_);
1407 ierr += HYPRE_IJMatrixSetObjectType(A21mat_, HYPRE_PARCSR);
1408 hypre_assert(!ierr);
1409
1410 //------------------------------------------------------------------
1411 // compute the number of nonzeros in the nConstraint row of A21
1412 //------------------------------------------------------------------
1413
1414 rowCount = maxRowSize = 0;
1415 newEndRow = endRow - nConstraints;
1416 if ( A21NRows > 0 ) A21MatSize = new int[A21NRows];
1417 else A21MatSize = NULL;
1418 for ( irow = endRow-nConstraints+1; irow <= endRow; irow++ )
1419 {
1420 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
1421 newRowSize = 0;
1422 for ( jcol = 0; jcol < rowSize; jcol++ )
1423 {
1424 colIndex = colInd[jcol];
1425 if (colVal[jcol] != 0.0 &&
1426 (colIndex <= newEndRow || colIndex > endRow)) newRowSize++;
1427 }
1428 A21MatSize[irow-newEndRow-1] = newRowSize;
1429 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
1430 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
1431 }
1432 nnzA21 = 0;
1433 for ( irow = 0; irow < nConstraints; irow++ ) nnzA21 += A21MatSize[irow];
1434 MPI_Allreduce(&nnzA21,&ncnt,1,MPI_INT,MPI_SUM,mpiComm_);
1435 if ( mypid == 0 && ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
1436 printf(" 0 : buildSubMatrices : NNZ of A21 = %d\n", ncnt);
1437
1438 //------------------------------------------------------------------
1439 // after fetching the row sizes, set up A21 with such sizes
1440 //------------------------------------------------------------------
1441
1442 ierr = HYPRE_IJMatrixSetRowSizes(A21mat_, A21MatSize);
1443 ierr += HYPRE_IJMatrixInitialize(A21mat_);
1444 hypre_assert(!ierr);
1445 if ( A21NRows > 0 ) delete [] A21MatSize;
1446
1447 //------------------------------------------------------------------
1448 // next load the first nConstraint row to A21 extracted from A
1449 // (at the same time, the D block is saved for future use)
1450 //------------------------------------------------------------------
1451
1452 rowCount = A21StartRow;
1453 newColInd = new int[maxRowSize+1];
1454 newColVal = new double[maxRowSize+1];
1455
1456 for ( irow = endRow-nConstraints+1; irow <= endRow; irow++ )
1457 {
1458 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
1459 newRowSize = 0;
1460 for ( jcol = 0; jcol < rowSize; jcol++ )
1461 {
1462 colIndex = colInd[jcol];
1463 if (colVal[jcol] != 0.0 &&
1464 (colIndex <= newEndRow || colIndex > endRow))
1465 {
1466 for ( procIndex = 0; procIndex < nprocs; procIndex++ )
1467 if ( procNRows[procIndex] > colIndex ) break;
1468 procIndex--;
1469 newColIndex = colInd[jcol] - procNConstr_[procIndex];
1470 newColInd[newRowSize] = newColIndex;
1471 newColVal[newRowSize++] = colVal[jcol];
1472 }
1473 }
1474 HYPRE_IJMatrixSetValues(A21mat_, 1, &newRowSize, (const int *) &rowCount,
1475 (const int *) newColInd, (const double *) newColVal);
1476 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
1477 rowCount++;
1478 }
1479 delete [] newColInd;
1480 delete [] newColVal;
1481
1482 HYPRE_IJMatrixAssemble(A21mat_);
1483 HYPRE_IJMatrixGetObject(A21mat_, (void **) &A21_csr);
1484 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) A21_csr);
1485
1486 //******************************************************************
1487 // extract submatrix of A not corresponding to constraints
1488 //------------------------------------------------------------------
1489
1490 reducedANRows = localNRows - nConstraints;
1491 reducedANCols = reducedANRows;
1492 reducedAStartRow = procNRows[mypid] - procNConstr_[mypid];
1493 reducedAStartCol = reducedAStartRow;
1494 reducedAGlobalNRows = globalNRows - globalNConstr;
1495 reducedAGlobalNCols = reducedAGlobalNRows;
1496
1497 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
1498 {
1499 printf("%4d : buildReducedMatrix - reduceAGlobalDim = %d %d\n", mypid,
1500 reducedAGlobalNRows, reducedAGlobalNCols);
1501 printf("%4d : buildReducedMatrix - reducedALocalDim = %d %d\n", mypid,
1502 reducedANRows, reducedANCols);
1503 }
1504
1505 //------------------------------------------------------------------
1506 // create a matrix context for reducedA
1507 //------------------------------------------------------------------
1508
1509 ierr = HYPRE_IJMatrixCreate(mpiComm_,reducedAStartRow,
1510 reducedAStartRow+reducedANRows-1, reducedAStartCol,
1511 reducedAStartCol+reducedANCols-1,&reducedAmat_);
1512 ierr += HYPRE_IJMatrixSetObjectType(reducedAmat_, HYPRE_PARCSR);
1513 hypre_assert(!ierr);
1514
1515 //------------------------------------------------------------------
1516 // compute row sizes for reducedA
1517 //------------------------------------------------------------------
1518
1519 reducedAMatSize = new int[reducedANRows];
1520 rowCount = maxRowSize = totalNNZ = 0;
1521 for ( irow = startRow; irow <= newEndRow; irow++ )
1522 {
1523 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
1524 newRowSize = 0;
1525 for ( jcol = 0; jcol < rowSize; jcol++ )
1526 {
1527 colIndex = colInd[jcol];
1528 for ( procIndex = 0; procIndex < nprocs; procIndex++ )
1529 if ( procNRows[procIndex] > colIndex ) break;
1530 uBound = procNRows[procIndex] -
1531 (procNConstr_[procIndex]-procNConstr_[procIndex-1]);
1532 procIndex--;
1533 if (colIndex < uBound) newRowSize++;
1534 }
1535 rowIndex = reducedAStartRow + rowCount;
1536 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
1537 reducedAMatSize[rowCount++] = newRowSize;
1538 totalNNZ += newRowSize;
1539 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
1540 }
1541 ierr = HYPRE_IJMatrixSetRowSizes(reducedAmat_, reducedAMatSize);
1542 ierr += HYPRE_IJMatrixInitialize(reducedAmat_);
1543 hypre_assert(!ierr);
1544 delete [] reducedAMatSize;
1545
1546 //------------------------------------------------------------------
1547 // load the reducedA matrix
1548 //------------------------------------------------------------------
1549
1550 rowCount = 0;
1551 newColInd = new int[maxRowSize+1];
1552 newColVal = new double[maxRowSize+1];
1553 for ( irow = startRow; irow <= newEndRow; irow++ )
1554 {
1555 HYPRE_ParCSRMatrixGetRow(A_csr, irow, &rowSize, &colInd, &colVal);
1556 newRowSize = 0;
1557 for (jcol = 0; jcol < rowSize; jcol++)
1558 {
1559 colIndex = colInd[jcol];
1560 for ( procIndex = 0; procIndex < nprocs; procIndex++ )
1561 if ( procNRows[procIndex] > colIndex ) break;
1562 uBound = procNRows[procIndex] -
1563 (procNConstr_[procIndex]-procNConstr_[procIndex-1]);
1564 procIndex--;
1565 if ( colIndex < uBound )
1566 {
1567 newColInd[newRowSize] = colIndex - procNConstr_[procIndex];
1568 newColVal[newRowSize++] = colVal[jcol];
1569 }
1570 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
1571 }
1572 rowIndex = reducedAStartRow + rowCount;
1573 ierr = HYPRE_IJMatrixSetValues(reducedAmat_, 1, &newRowSize,
1574 (const int *) &rowIndex, (const int *) newColInd,
1575 (const double *) newColVal);
1576 hypre_assert(!ierr);
1577 rowCount++;
1578 }
1579 delete [] newColInd;
1580 delete [] newColVal;
1581
1582 free( procNRows );
1583
1584 //------------------------------------------------------------------
1585 // assemble the reduced matrix
1586 //------------------------------------------------------------------
1587
1588 HYPRE_IJMatrixAssemble(reducedAmat_);
1589 HYPRE_IJMatrixGetObject(reducedAmat_, (void **) &reducedA_csr);
1590
1591 return 0;
1592 }
1593
1594 //****************************************************************************
1595 // build reduced rhs vector (form red_f1 = f1 - A12*x2)
1596 //----------------------------------------------------------------------------
1597
buildModifiedRHSVector(HYPRE_IJVector x,HYPRE_IJVector b)1598 int HYPRE_SlideReduction::buildModifiedRHSVector(HYPRE_IJVector x,
1599 HYPRE_IJVector b)
1600 {
1601 int nprocs, mypid, *procNRows, startRow, endRow, localNRows;
1602 int nConstraints, redBNRows, redBStart, ierr, x2NRows, x2Start;
1603 int vecIndex, irow;
1604 double *b_data, *rb_data, *x_data, *x2_data;
1605 HYPRE_ParCSRMatrix A_csr, A21_csr;
1606 HYPRE_IJVector x2;
1607 HYPRE_ParVector b_csr, rb_csr, x_csr, x2_csr;
1608 hypre_Vector *b_local, *rb_local, *x_local, *x2_local;
1609
1610 //------------------------------------------------------------------
1611 // sanitize
1612 //------------------------------------------------------------------
1613
1614 if (reducedBvec_ != NULL ) HYPRE_IJVectorDestroy(reducedBvec_);
1615 if (reducedXvec_ != NULL ) HYPRE_IJVectorDestroy(reducedXvec_);
1616 if (reducedRvec_ != NULL ) HYPRE_IJVectorDestroy(reducedRvec_);
1617 reducedBvec_ = NULL;
1618 reducedXvec_ = NULL;
1619 reducedRvec_ = NULL;
1620
1621 //------------------------------------------------------------------
1622 // get machine and matrix information
1623 //------------------------------------------------------------------
1624
1625 if (reducedAmat_ == NULL) return 0;
1626 MPI_Comm_rank( mpiComm_, &mypid );
1627 MPI_Comm_size( mpiComm_, &nprocs );
1628 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
1629 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
1630 startRow = procNRows[mypid];
1631 endRow = procNRows[mypid+1] - 1;
1632 localNRows = endRow - startRow + 1;
1633 if (procNConstr_ == NULL || procNConstr_[nprocs] == 0)
1634 {
1635 printf("%4d : buildModifiedRHSVector WARNING - no local data.\n",mypid);
1636 free(procNRows);
1637 return 1;
1638 }
1639
1640 //------------------------------------------------------------------
1641 // create reducedB
1642 //------------------------------------------------------------------
1643
1644 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
1645 redBNRows = localNRows - nConstraints;
1646 redBStart = procNRows[mypid] - procNConstr_[mypid];
1647 ierr = HYPRE_IJVectorCreate(mpiComm_, redBStart,
1648 redBStart+redBNRows-1, &reducedBvec_);
1649 ierr += HYPRE_IJVectorSetObjectType(reducedBvec_, HYPRE_PARCSR);
1650 ierr += HYPRE_IJVectorInitialize(reducedBvec_);
1651 ierr += HYPRE_IJVectorAssemble(reducedBvec_);
1652 hypre_assert( !ierr );
1653 HYPRE_IJVectorGetObject(reducedBvec_, (void **) &rb_csr);
1654 HYPRE_IJVectorGetObject(b, (void **) &b_csr);
1655 b_local = hypre_ParVectorLocalVector((hypre_ParVector *) b_csr);
1656 b_data = (double *) hypre_VectorData(b_local);
1657 rb_local = hypre_ParVectorLocalVector((hypre_ParVector *) rb_csr);
1658 rb_data = (double *) hypre_VectorData(rb_local);
1659 for ( irow = 0; irow < localNRows-nConstraints; irow++ )
1660 rb_data[irow] = b_data[irow];
1661
1662 //------------------------------------------------------------------
1663 // create x_2
1664 //------------------------------------------------------------------
1665
1666 x2NRows = nConstraints;
1667 x2Start = procNConstr_[mypid];
1668 HYPRE_IJVectorCreate(mpiComm_, x2Start, x2Start+x2NRows-1, &x2);
1669 HYPRE_IJVectorSetObjectType(x2, HYPRE_PARCSR);
1670 ierr = HYPRE_IJVectorInitialize(x2);
1671 ierr += HYPRE_IJVectorAssemble(x2);
1672 hypre_assert(!ierr);
1673 HYPRE_IJVectorGetObject(x2, (void **) &x2_csr);
1674 HYPRE_IJVectorGetObject(x, (void **) &x_csr);
1675 x_local = hypre_ParVectorLocalVector((hypre_ParVector *) x_csr);
1676 x_data = (double *) hypre_VectorData(x_local);
1677 x2_local = hypre_ParVectorLocalVector((hypre_ParVector *) x2_csr);
1678 x2_data = (double *) hypre_VectorData(x2_local);
1679 for ( irow = 0; irow < nConstraints; irow++ )
1680 {
1681 vecIndex = localNRows - nConstraints + irow;
1682 x2_data[irow] = x_data[vecIndex];
1683 }
1684
1685 //------------------------------------------------------------------
1686 // form reducedB = A21^T * u_2
1687 //------------------------------------------------------------------
1688
1689 HYPRE_IJMatrixGetObject(A21mat_, (void **) &A21_csr);
1690 HYPRE_ParCSRMatrixMatvecT(-1.0, A21_csr, x2_csr, 1.0, rb_csr);
1691 HYPRE_IJVectorDestroy(x2);
1692
1693 //------------------------------------------------------------------
1694 // create a few more vectors
1695 //------------------------------------------------------------------
1696
1697 ierr = HYPRE_IJVectorCreate(mpiComm_, redBStart,
1698 redBStart+redBNRows-1, &reducedXvec_);
1699 ierr += HYPRE_IJVectorSetObjectType(reducedXvec_, HYPRE_PARCSR);
1700 ierr += HYPRE_IJVectorInitialize(reducedXvec_);
1701 ierr += HYPRE_IJVectorAssemble(reducedXvec_);
1702 hypre_assert( !ierr );
1703
1704 ierr = HYPRE_IJVectorCreate(mpiComm_, redBStart,
1705 redBStart+redBNRows-1, &reducedRvec_);
1706 ierr += HYPRE_IJVectorSetObjectType(reducedRvec_, HYPRE_PARCSR);
1707 ierr += HYPRE_IJVectorInitialize(reducedRvec_);
1708 ierr += HYPRE_IJVectorAssemble(reducedRvec_);
1709 hypre_assert( !ierr );
1710 free( procNRows );
1711
1712 return 0;
1713 }
1714
1715 //*****************************************************************************
1716 // given the solution vector, copy the actual solution
1717 //-----------------------------------------------------------------------------
1718
buildModifiedSolnVector(HYPRE_IJVector x)1719 int HYPRE_SlideReduction::buildModifiedSolnVector(HYPRE_IJVector x)
1720 {
1721 int mypid, nprocs, *procNRows, startRow, endRow, localNRows;
1722 int nConstraints, irow;
1723 double *x_data, *rx_data;
1724 HYPRE_ParCSRMatrix A_csr;
1725 HYPRE_ParVector x_csr, rx_csr;
1726 hypre_Vector *x_local, *rx_local;
1727
1728 //------------------------------------------------------------------
1729 // get machine and matrix information
1730 //------------------------------------------------------------------
1731
1732 if ( reducedXvec_ == NULL ) return -1;
1733 MPI_Comm_rank( mpiComm_, &mypid );
1734 MPI_Comm_size( mpiComm_, &nprocs );
1735 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
1736 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
1737 startRow = procNRows[mypid];
1738 endRow = procNRows[mypid+1] - 1;
1739 localNRows = endRow - startRow + 1;
1740 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
1741 free( procNRows );
1742 if (( outputLevel_ & HYPRE_BITMASK2 ) >= 1 &&
1743 (procNConstr_==NULL || procNConstr_[nprocs]==0))
1744 {
1745 printf("%4d : buildModifiedSolnVector WARNING - no local entry.\n",
1746 mypid);
1747 return 1;
1748 }
1749
1750 //------------------------------------------------------------------
1751 // compute b2 - A21 * sol (v1 = b2 + v1)
1752 //------------------------------------------------------------------
1753
1754 HYPRE_IJVectorGetObject(x, (void **) &x_csr);
1755 x_local = hypre_ParVectorLocalVector((hypre_ParVector *) x_csr);
1756 x_data = (double *) hypre_VectorData(x_local);
1757 HYPRE_IJVectorGetObject(reducedXvec_, (void **) &rx_csr);
1758 rx_local = hypre_ParVectorLocalVector((hypre_ParVector *) rx_csr);
1759 rx_data = (double *) hypre_VectorData(rx_local);
1760 for ( irow = 0; irow < localNRows-nConstraints; irow++ )
1761 x_data[irow] = rx_data[irow];
1762
1763 return 0;
1764 }
1765
1766 //****************************************************************************
1767 // build the reduced matrix
1768 //----------------------------------------------------------------------------
1769
buildReducedMatrix()1770 int HYPRE_SlideReduction::buildReducedMatrix()
1771 {
1772 int mypid, nprocs, *procNRows, startRow, endRow, localNRows;
1773 int globalNConstr, globalNRows, nConstraints, newEndRow;
1774 int reducedANRows, reducedANCols, reducedAStartRow, reducedAStartCol;
1775 int reducedAGlobalNRows, reducedAGlobalNCols, ncnt, irow, jcol;
1776 int rowSize, *colInd, *reducedAMatSize, rowCount, maxRowSize;
1777 int rowSize2, *colInd2, newRowSize, rowIndex, searchIndex, uBound;
1778 int procIndex, colIndex, ierr, *newColInd, totalNNZ;
1779 double *colVal, *colVal2, *newColVal;
1780 char fname[40];
1781 HYPRE_ParCSRMatrix A_csr, A21_csr, invA22_csr, RAP_csr, reducedA_csr;
1782
1783 //------------------------------------------------------------------
1784 // first compute A21 and invA22
1785 //------------------------------------------------------------------
1786
1787 buildA21Mat();
1788 buildInvA22Mat();
1789
1790 //------------------------------------------------------------------
1791 // get matrix information
1792 //------------------------------------------------------------------
1793
1794 MPI_Comm_rank(mpiComm_, &mypid);
1795 MPI_Comm_size(mpiComm_, &nprocs);
1796 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
1797 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
1798 startRow = procNRows[mypid];
1799 endRow = procNRows[mypid+1] - 1;
1800 localNRows = endRow - startRow + 1;
1801 globalNConstr = procNConstr_[nprocs];
1802 globalNRows = procNRows[nprocs];
1803 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
1804 newEndRow = endRow - nConstraints;
1805
1806 reducedANRows = localNRows - nConstraints;
1807 reducedANCols = reducedANRows;
1808 reducedAStartRow = procNRows[mypid] - procNConstr_[mypid];
1809 reducedAStartCol = reducedAStartRow;
1810 reducedAGlobalNRows = globalNRows - globalNConstr;
1811 reducedAGlobalNCols = reducedAGlobalNRows;
1812
1813 //------------------------------------------------------------------
1814 // perform the triple matrix product A12 * invA22 * A21
1815 //------------------------------------------------------------------
1816
1817 HYPRE_IJMatrixGetObject(A21mat_, (void **) &A21_csr);
1818 HYPRE_IJMatrixGetObject(invA22mat_, (void **) &invA22_csr);
1819 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
1820 printf("%4d : buildReducedMatrix - Triple matrix product starts\n",
1821 mypid);
1822
1823 hypre_BoomerAMGBuildCoarseOperator((hypre_ParCSRMatrix *) A21_csr,
1824 (hypre_ParCSRMatrix *) invA22_csr,
1825 (hypre_ParCSRMatrix *) A21_csr,
1826 (hypre_ParCSRMatrix **) &RAP_csr);
1827
1828 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
1829 printf("%4d : buildReducedMatrix - Triple matrix product ends\n",
1830 mypid);
1831
1832 if ( outputLevel_ >= 4 )
1833 {
1834 sprintf(fname, "rap.%d", mypid);
1835 FILE *fp = fopen(fname, "w");
1836
1837 if ( mypid == 0 )
1838 {
1839 printf("====================================================\n");
1840 printf("%4d : Printing RAP matrix... \n", mypid);
1841 fflush(stdout);
1842 }
1843 for (irow=reducedAStartRow; irow<reducedAStartRow+reducedANRows;irow++)
1844 {
1845 HYPRE_ParCSRMatrixGetRow(RAP_csr,irow,&rowSize,&colInd,&colVal);
1846 for ( jcol = 0; jcol < rowSize; jcol++ )
1847 if ( colVal[jcol] != 0.0 )
1848 fprintf(fp,"%6d %6d %25.16e \n",irow+1,colInd[jcol]+1,
1849 colVal[jcol]);
1850 HYPRE_ParCSRMatrixRestoreRow(RAP_csr,irow,&rowSize,&colInd,&colVal);
1851 }
1852 fclose(fp);
1853 if ( mypid == 0 )
1854 printf("====================================================\n");
1855 }
1856
1857 //******************************************************************
1858 // form reduceA = A11 - A12 * invA22 * A21
1859 //------------------------------------------------------------------
1860
1861 //------------------------------------------------------------------
1862 // compute row sizes
1863 //------------------------------------------------------------------
1864
1865 reducedAMatSize = new int[reducedANRows];
1866
1867 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
1868 {
1869 printf("%4d : buildReducedMatrix - reduceAGlobalDim = %d %d\n", mypid,
1870 reducedAGlobalNRows, reducedAGlobalNCols);
1871 printf("%4d : buildReducedMatrix - reducedALocalDim = %d %d\n", mypid,
1872 reducedANRows, reducedANCols);
1873 }
1874
1875 //------------------------------------------------------------------
1876 // create a matrix context for reducedA
1877 //------------------------------------------------------------------
1878
1879 ierr = HYPRE_IJMatrixCreate(mpiComm_,reducedAStartRow,
1880 reducedAStartRow+reducedANRows-1, reducedAStartCol,
1881 reducedAStartCol+reducedANCols-1,&reducedAmat_);
1882 ierr += HYPRE_IJMatrixSetObjectType(reducedAmat_, HYPRE_PARCSR);
1883 hypre_assert(!ierr);
1884
1885 //------------------------------------------------------------------
1886 // compute row sizes for reducedA
1887 //------------------------------------------------------------------
1888
1889 rowCount = maxRowSize = totalNNZ = 0;
1890 for ( irow = startRow; irow <= newEndRow; irow++ )
1891 {
1892 searchIndex = hypre_BinarySearch(slaveEqnList_, irow, nConstraints);
1893 if ( searchIndex >= 0 ) reducedAMatSize[rowCount++] = 1;
1894 else
1895 {
1896 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
1897 rowIndex = reducedAStartRow + rowCount;
1898 ierr = HYPRE_ParCSRMatrixGetRow(RAP_csr,rowIndex,&rowSize2,
1899 &colInd2, &colVal2);
1900 hypre_assert( !ierr );
1901 newRowSize = rowSize + rowSize2;
1902 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
1903 newColInd = new int[newRowSize];
1904 for (jcol = 0; jcol < rowSize; jcol++)
1905 newColInd[jcol] = colInd[jcol];
1906 for (jcol = 0; jcol < rowSize2; jcol++)
1907 newColInd[rowSize+jcol] = colInd2[jcol];
1908 hypre_qsort0(newColInd, 0, newRowSize-1);
1909 ncnt = 0;
1910 for ( jcol = 1; jcol < newRowSize; jcol++ )
1911 if (newColInd[jcol] != newColInd[ncnt])
1912 newColInd[++ncnt] = newColInd[jcol];
1913 if ( newRowSize > 0 ) ncnt++;
1914 reducedAMatSize[rowCount++] = ncnt;
1915 totalNNZ += ncnt;
1916 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
1917 ierr = HYPRE_ParCSRMatrixRestoreRow(RAP_csr,rowIndex,&rowSize2,
1918 &colInd2,&colVal2);
1919 delete [] newColInd;
1920 hypre_assert( !ierr );
1921 }
1922 }
1923 ierr = HYPRE_IJMatrixSetRowSizes(reducedAmat_, reducedAMatSize);
1924 ierr += HYPRE_IJMatrixInitialize(reducedAmat_);
1925 hypre_assert(!ierr);
1926 delete [] reducedAMatSize;
1927
1928 int totalNNZA = 0;
1929 for ( irow = startRow; irow <= endRow; irow++ )
1930 {
1931 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,NULL,NULL);
1932 totalNNZA += rowSize;
1933 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,NULL,NULL);
1934 }
1935 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
1936 {
1937 printf("%4d : buildReducedMatrix - NNZ of reducedA = %d %d %e\n", mypid,
1938 totalNNZ, totalNNZA, 1.0*totalNNZ/totalNNZA);
1939 }
1940
1941 //------------------------------------------------------------------
1942 // load the reducedA matrix
1943 //------------------------------------------------------------------
1944
1945 rowCount = 0;
1946 newColInd = new int[maxRowSize+1];
1947 newColVal = new double[maxRowSize+1];
1948 for ( irow = startRow; irow <= newEndRow; irow++ )
1949 {
1950 searchIndex = hypre_BinarySearch(slaveEqnList_, irow, nConstraints);
1951 rowIndex = reducedAStartRow + rowCount;
1952 if ( searchIndex >= 0 )
1953 {
1954 newRowSize = 1;
1955 newColInd[0] = reducedAStartRow + rowCount;
1956 newColVal[0] = 1.0;
1957 }
1958 else
1959 {
1960 HYPRE_ParCSRMatrixGetRow(A_csr, irow, &rowSize, &colInd, &colVal);
1961 HYPRE_ParCSRMatrixGetRow(RAP_csr,rowIndex,&rowSize2,&colInd2,
1962 &colVal2);
1963 newRowSize = rowSize + rowSize2;
1964 ncnt = 0;
1965 for ( jcol = 0; jcol < rowSize; jcol++ )
1966 {
1967 colIndex = colInd[jcol];
1968 for ( procIndex = 0; procIndex < nprocs; procIndex++ )
1969 if ( procNRows[procIndex] > colIndex ) break;
1970 uBound = procNRows[procIndex] -
1971 (procNConstr_[procIndex]-procNConstr_[procIndex-1]);
1972 procIndex--;
1973 if ( colIndex < uBound )
1974 {
1975 searchIndex = hypre_BinarySearch(gSlaveEqnList_, colIndex,
1976 globalNConstr);
1977 if ( searchIndex < 0 )
1978 {
1979 newColInd[ncnt] = colIndex - procNConstr_[procIndex];
1980 newColVal[ncnt++] = colVal[jcol];
1981 }
1982 }
1983 }
1984 for ( jcol = 0; jcol < rowSize2; jcol++ )
1985 {
1986 newColInd[ncnt+jcol] = colInd2[jcol];
1987 newColVal[ncnt+jcol] = - colVal2[jcol];
1988 }
1989 newRowSize = ncnt + rowSize2;
1990 hypre_qsort1(newColInd, newColVal, 0, newRowSize-1);
1991 ncnt = 0;
1992 for ( jcol = 0; jcol < newRowSize; jcol++ )
1993 {
1994 if ( jcol != ncnt && newColInd[jcol] == newColInd[ncnt] )
1995 newColVal[ncnt] += newColVal[jcol];
1996 else if ( newColInd[jcol] != newColInd[ncnt] )
1997 {
1998 ncnt++;
1999 newColVal[ncnt] = newColVal[jcol];
2000 newColInd[ncnt] = newColInd[jcol];
2001 }
2002 }
2003 newRowSize = ncnt + 1;
2004 ncnt = 0;
2005 for ( jcol = 0; jcol < newRowSize; jcol++ )
2006 {
2007 if ( habs(newColVal[jcol]) >= truncTol_ )
2008 {
2009 newColInd[ncnt] = newColInd[jcol];
2010 newColVal[ncnt++] = newColVal[jcol];
2011 }
2012 }
2013 newRowSize = ncnt;
2014 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
2015 HYPRE_ParCSRMatrixRestoreRow(RAP_csr,rowIndex,&rowSize2,&colInd2,
2016 &colVal2);
2017 }
2018 ierr = HYPRE_IJMatrixSetValues(reducedAmat_, 1, &newRowSize,
2019 (const int *) &rowIndex, (const int *) newColInd,
2020 (const double *) newColVal);
2021 hypre_assert(!ierr);
2022 rowCount++;
2023 }
2024 delete [] newColInd;
2025 delete [] newColVal;
2026 hypreRAP_ = RAP_csr;
2027 free( procNRows );
2028
2029 //------------------------------------------------------------------
2030 // assemble the reduced matrix
2031 //------------------------------------------------------------------
2032
2033 HYPRE_IJMatrixAssemble(reducedAmat_);
2034 HYPRE_IJMatrixGetObject(reducedAmat_, (void **) &reducedA_csr);
2035
2036 if ( outputLevel_ >= 5 )
2037 {
2038 sprintf(fname, "reducedA.%d", mypid);
2039 FILE *fp = fopen(fname, "w");
2040
2041 if ( mypid == 0 )
2042 {
2043 printf("====================================================\n");
2044 printf("%4d : Printing reducedA matrix... \n", mypid);
2045 fflush(stdout);
2046 }
2047 for ( irow = reducedAStartRow;
2048 irow < reducedAStartRow+localNRows-nConstraints; irow++ )
2049 {
2050 //printf("%d : reducedA ROW %d\n", mypid, irow);
2051 ierr = HYPRE_ParCSRMatrixGetRow(reducedA_csr,irow,&rowSize,
2052 &colInd, &colVal);
2053 //hypre_qsort1(colInd, colVal, 0, rowSize-1);
2054 for ( jcol = 0; jcol < rowSize; jcol++ )
2055 if ( colVal[jcol] != 0.0 )
2056 fprintf(fp,"%6d %6d %25.16e \n",irow+1,colInd[jcol]+1,
2057 colVal[jcol]);
2058 HYPRE_ParCSRMatrixRestoreRow(reducedA_csr,irow,&rowSize,&colInd,
2059 &colVal);
2060 }
2061 fclose(fp);
2062 if ( mypid == 0 )
2063 printf("====================================================\n");
2064 }
2065 return 0;
2066 }
2067
2068 //****************************************************************************
2069 // build reduced rhs vector (form red_f1 = f1 - A12*invA22*f2)
2070 //----------------------------------------------------------------------------
2071
buildReducedRHSVector(HYPRE_IJVector b)2072 int HYPRE_SlideReduction::buildReducedRHSVector(HYPRE_IJVector b)
2073 {
2074 int nprocs, mypid, *procNRows, startRow, endRow, localNRows;
2075 int nConstraints, newEndRow, f2LocalLength, f2Start, ierr;
2076 int irow, jcol, vecIndex, redBLocalLength, redBStart, rowIndex;
2077 double *b_data, *f2_data, ddata;
2078 HYPRE_ParCSRMatrix A_csr, A21_csr, invA22_csr;
2079 HYPRE_IJVector f2, f2hat;
2080 HYPRE_ParVector b_csr, f2_csr, f2hat_csr, rb_csr;
2081 hypre_Vector *b_local, *f2_local;
2082
2083 //------------------------------------------------------------------
2084 // get machine and matrix information
2085 //------------------------------------------------------------------
2086
2087 if ( reducedAmat_ == NULL ) return 0;
2088 MPI_Comm_rank( mpiComm_, &mypid );
2089 MPI_Comm_size( mpiComm_, &nprocs );
2090 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
2091 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
2092 startRow = procNRows[mypid];
2093 endRow = procNRows[mypid+1] - 1;
2094 localNRows = endRow - startRow + 1;
2095 if ( procNConstr_ == NULL || procNConstr_[nprocs] == 0 )
2096 {
2097 printf("%4d : buildReducedRHSVector WARNING - no local entries.\n",mypid);
2098 free(procNRows);
2099 return 1;
2100 }
2101
2102 //------------------------------------------------------------------
2103 // form f2hat = invA22 * f2
2104 //------------------------------------------------------------------
2105
2106 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
2107 newEndRow = endRow - nConstraints;
2108 f2LocalLength = 2 * nConstraints;
2109 f2Start = 2 * procNConstr_[mypid];
2110
2111 HYPRE_IJVectorCreate(mpiComm_, f2Start, f2Start+f2LocalLength-1, &f2);
2112 HYPRE_IJVectorSetObjectType(f2, HYPRE_PARCSR);
2113 ierr = HYPRE_IJVectorInitialize(f2);
2114 ierr += HYPRE_IJVectorAssemble(f2);
2115 hypre_assert(!ierr);
2116 HYPRE_IJVectorGetObject(f2, (void **) &f2_csr);
2117
2118 HYPRE_IJVectorCreate(mpiComm_, f2Start, f2Start+f2LocalLength-1, &f2hat);
2119 HYPRE_IJVectorSetObjectType(f2hat, HYPRE_PARCSR);
2120 ierr = HYPRE_IJVectorInitialize(f2hat);
2121 ierr += HYPRE_IJVectorAssemble(f2hat);
2122 hypre_assert(!ierr);
2123 HYPRE_IJVectorGetObject(f2hat, (void **) &f2hat_csr);
2124
2125 HYPRE_IJVectorGetObject(b, (void **) &b_csr);
2126 b_local = hypre_ParVectorLocalVector((hypre_ParVector *) b_csr);
2127 b_data = (double *) hypre_VectorData(b_local);
2128 f2_local = hypre_ParVectorLocalVector((hypre_ParVector *) f2_csr);
2129 f2_data = (double *) hypre_VectorData(f2_local);
2130
2131 for ( irow = 0; irow < nConstraints; irow++ )
2132 {
2133 vecIndex = -1;
2134 for ( jcol = 0; jcol < nConstraints; jcol++ )
2135 {
2136 if ( slaveEqnListAux_[jcol] == irow )
2137 {
2138 vecIndex = slaveEqnList_[jcol];
2139 break;
2140 }
2141 }
2142 hypre_assert( vecIndex >= startRow );
2143 hypre_assert( vecIndex <= endRow );
2144 f2_data[irow] = b_data[vecIndex-startRow];
2145 }
2146 for ( irow = 0; irow < nConstraints; irow++ )
2147 {
2148 vecIndex = localNRows - nConstraints + irow;
2149 f2_data[irow+nConstraints] = b_data[vecIndex];
2150 }
2151
2152 HYPRE_IJMatrixGetObject(invA22mat_, (void **) &invA22_csr);
2153 HYPRE_ParCSRMatrixMatvec( 1.0, invA22_csr, f2_csr, 0.0, f2hat_csr );
2154 HYPRE_IJVectorDestroy(f2);
2155
2156 //------------------------------------------------------------------
2157 // form reducedB = A21^T * f2hat
2158 //------------------------------------------------------------------
2159
2160 redBLocalLength = localNRows - nConstraints;
2161 redBStart = procNRows[mypid] - procNConstr_[mypid];
2162
2163 ierr = HYPRE_IJVectorCreate(mpiComm_, redBStart,
2164 redBStart+redBLocalLength-1, &reducedBvec_);
2165 ierr += HYPRE_IJVectorSetObjectType(reducedBvec_, HYPRE_PARCSR);
2166 ierr += HYPRE_IJVectorInitialize(reducedBvec_);
2167 ierr += HYPRE_IJVectorAssemble(reducedBvec_);
2168 hypre_assert( !ierr );
2169
2170 HYPRE_IJVectorGetObject(reducedBvec_, (void **) &rb_csr);
2171 HYPRE_IJMatrixGetObject(A21mat_, (void **) &A21_csr);
2172 HYPRE_ParCSRMatrixMatvecT(-1.0, A21_csr, f2hat_csr, 0.0, rb_csr);
2173 HYPRE_IJVectorDestroy(f2hat);
2174
2175 //------------------------------------------------------------------
2176 // finally form reducedB = f1 - f2til
2177 //------------------------------------------------------------------
2178
2179 rowIndex = redBStart;
2180 for ( irow = startRow; irow <= newEndRow; irow++ )
2181 {
2182 if ( hypre_BinarySearch(slaveEqnList_, irow, nConstraints) < 0 )
2183 {
2184 ddata = b_data[irow-startRow];
2185 HYPRE_IJVectorAddToValues(reducedBvec_, 1, (const int *) &rowIndex,
2186 (const double *) &ddata);
2187 }
2188 else
2189 {
2190 ddata = 0.0;
2191 HYPRE_IJVectorSetValues(reducedBvec_, 1, (const int *) &rowIndex,
2192 (const double *) &ddata);
2193 }
2194 rowIndex++;
2195 }
2196 HYPRE_IJVectorGetObject(reducedBvec_, (void **) &rb_csr);
2197
2198 //------------------------------------------------------------------
2199 // create a few more vectors
2200 //------------------------------------------------------------------
2201
2202 ierr = HYPRE_IJVectorCreate(mpiComm_, redBStart,
2203 redBStart+redBLocalLength-1, &reducedXvec_);
2204 ierr += HYPRE_IJVectorSetObjectType(reducedXvec_, HYPRE_PARCSR);
2205 ierr += HYPRE_IJVectorInitialize(reducedXvec_);
2206 ierr += HYPRE_IJVectorAssemble(reducedXvec_);
2207 hypre_assert( !ierr );
2208
2209 ierr = HYPRE_IJVectorCreate(mpiComm_, redBStart,
2210 redBStart+redBLocalLength-1, &reducedRvec_);
2211 ierr += HYPRE_IJVectorSetObjectType(reducedRvec_, HYPRE_PARCSR);
2212 ierr += HYPRE_IJVectorInitialize(reducedRvec_);
2213 ierr += HYPRE_IJVectorAssemble(reducedRvec_);
2214 hypre_assert( !ierr );
2215 free( procNRows );
2216
2217 return 0;
2218 }
2219
2220 //*****************************************************************************
2221 // given the submatrices and the solution vector, restore the actual solution
2222 // A21 x_1 + A22 x_2 = b2
2223 // x_2 = invA22 * ( b2 - A21 x_1 )
2224 //-----------------------------------------------------------------------------
2225
buildReducedSolnVector(HYPRE_IJVector x,HYPRE_IJVector b)2226 int HYPRE_SlideReduction::buildReducedSolnVector(HYPRE_IJVector x,
2227 HYPRE_IJVector b)
2228 {
2229 int mypid, nprocs, *procNRows, startRow, endRow, localNRows;
2230 int nConstraints, newEndRow, vecStart, vecLocalLength, ierr;
2231 int irow, jcol, rowIndex, searchIndex, length;
2232 double *b_data, *v1_data, *rx_data, *x_data, *x2_data;
2233 HYPRE_ParCSRMatrix A_csr, A21_csr, invA22_csr;
2234 HYPRE_ParVector x_csr, x2_csr, v1_csr, b_csr, rx_csr;
2235 HYPRE_IJVector v1, x2;
2236 hypre_Vector *b_local, *v1_local, *rx_local, *x_local, *x2_local;
2237
2238 //------------------------------------------------------------------
2239 // get machine and matrix information
2240 //------------------------------------------------------------------
2241
2242 if ( reducedAmat_ == NULL ) return 0;
2243 MPI_Comm_rank( mpiComm_, &mypid );
2244 MPI_Comm_size( mpiComm_, &nprocs );
2245 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
2246 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
2247 startRow = procNRows[mypid];
2248 endRow = procNRows[mypid+1] - 1;
2249 localNRows = endRow - startRow + 1;
2250 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
2251 newEndRow = endRow - nConstraints;
2252 if (( outputLevel_ & HYPRE_BITMASK2 ) >= 1 &&
2253 (procNConstr_==NULL || procNConstr_[nprocs]==0))
2254 {
2255 printf("%4d : buildReducedSolnVector WARNING - no local entry.\n",mypid);
2256 return 1;
2257 }
2258
2259 //------------------------------------------------------------------
2260 // compute v1 = - A21 * sol
2261 //------------------------------------------------------------------
2262
2263 vecStart = 2 * procNConstr_[mypid];
2264 vecLocalLength = 2 * nConstraints;
2265 ierr = HYPRE_IJVectorCreate(mpiComm_, vecStart,
2266 vecStart+vecLocalLength-1, &v1);
2267 ierr += HYPRE_IJVectorSetObjectType(v1, HYPRE_PARCSR);
2268 ierr += HYPRE_IJVectorInitialize(v1);
2269 ierr += HYPRE_IJVectorAssemble(v1);
2270 hypre_assert(!ierr);
2271 HYPRE_IJVectorGetObject(v1, (void **) &v1_csr);
2272 HYPRE_IJMatrixGetObject(A21mat_, (void **) &A21_csr);
2273 HYPRE_IJVectorGetObject(reducedXvec_, (void **) &rx_csr);
2274 if ( scaleMatrixFlag_ == 1 && ADiagISqrts_ != NULL )
2275 {
2276 rx_local = hypre_ParVectorLocalVector((hypre_ParVector *) rx_csr);
2277 rx_data = (double *) hypre_VectorData(rx_local);
2278 length = hypre_VectorSize(rx_local);
2279 for ( irow = 0; irow < length; irow++ )
2280 rx_data[irow] *= ADiagISqrts_[irow];
2281 }
2282 HYPRE_ParCSRMatrixMatvec( -1.0, A21_csr, rx_csr, 0.0, v1_csr );
2283
2284 //------------------------------------------------------------------
2285 // compute b2 - A21 * sol (v1 = b2 + v1)
2286 //------------------------------------------------------------------
2287
2288 HYPRE_IJVectorGetObject(b, (void **) &b_csr);
2289 b_local = hypre_ParVectorLocalVector((hypre_ParVector *) b_csr);
2290 b_data = (double *) hypre_VectorData(b_local);
2291 v1_local = hypre_ParVectorLocalVector((hypre_ParVector *) v1_csr);
2292 v1_data = (double *) hypre_VectorData(v1_local);
2293
2294 rowIndex = 0;
2295 for ( irow = 0; irow < nConstraints; irow++ )
2296 {
2297 searchIndex = -1;
2298 for ( jcol = 0; jcol < nConstraints; jcol++ )
2299 {
2300 if ( slaveEqnListAux_[jcol] == irow )
2301 {
2302 searchIndex = slaveEqnList_[jcol];
2303 break;
2304 }
2305 }
2306 hypre_assert( searchIndex >= startRow );
2307 hypre_assert( searchIndex <= newEndRow );
2308 v1_data[rowIndex++] += b_data[searchIndex-startRow];
2309 }
2310 for ( irow = endRow-nConstraints+1; irow <= endRow; irow++ )
2311 v1_data[rowIndex++] += b_data[irow-startRow];
2312
2313 //-------------------------------------------------------------
2314 // compute inv(A22) * (f2 - A21 * sol) --> x2 = invA22 * v1
2315 //-------------------------------------------------------------
2316
2317 ierr = HYPRE_IJVectorCreate(mpiComm_, vecStart,
2318 vecStart+vecLocalLength-1, &x2);
2319 ierr += HYPRE_IJVectorSetObjectType(x2, HYPRE_PARCSR);
2320 ierr += HYPRE_IJVectorInitialize(x2);
2321 ierr += HYPRE_IJVectorAssemble(x2);
2322 hypre_assert(!ierr);
2323 HYPRE_IJVectorGetObject(x2, (void **) &x2_csr );
2324 HYPRE_IJMatrixGetObject(invA22mat_, (void **) &invA22_csr );
2325 HYPRE_ParCSRMatrixMatvec(1.0, invA22_csr, v1_csr, 0.0, x2_csr);
2326
2327 #if 0
2328 FILE *fp = fopen("rhs.m", "w");
2329 for ( irow = 0; irow < 2*nConstraints; irow++ )
2330 fprintf(fp, " %6d %25.16e \n", irow+1, v1_data[irow]);
2331 fclose(fp);
2332 #endif
2333
2334 //-------------------------------------------------------------
2335 // inject final solution to the solution vector x
2336 //-------------------------------------------------------------
2337
2338 HYPRE_IJVectorGetObject(x, (void **) &x_csr );
2339 rx_local = hypre_ParVectorLocalVector((hypre_ParVector *) rx_csr);
2340 rx_data = (double *) hypre_VectorData(rx_local);
2341 x_local = hypre_ParVectorLocalVector((hypre_ParVector *) x_csr);
2342 x_data = (double *) hypre_VectorData(x_local);
2343 x2_local = hypre_ParVectorLocalVector((hypre_ParVector *) x2_csr);
2344 x2_data = (double *) hypre_VectorData(x2_local);
2345
2346 for ( irow = 0; irow < localNRows-nConstraints; irow++ )
2347 x_data[irow] = rx_data[irow];
2348
2349 for ( irow = 0; irow < nConstraints; irow++ )
2350 {
2351 for ( jcol = 0; jcol < nConstraints; jcol++ )
2352 {
2353 if ( slaveEqnListAux_[jcol] == irow )
2354 {
2355 searchIndex = slaveEqnList_[jcol];
2356 break;
2357 }
2358 }
2359 x_data[searchIndex-startRow] = x2_data[irow];
2360 }
2361 for ( irow = nConstraints; irow < 2*nConstraints; irow++ )
2362 x_data[localNRows-2*nConstraints+irow] = x2_data[irow];
2363
2364 //------------------------------------------------------------------
2365 // compute true residual
2366 //------------------------------------------------------------------
2367
2368 #if 0
2369 double rnorm;
2370 HYPRE_IJVector R;
2371 HYPRE_ParVector R_csr;
2372
2373 ierr = HYPRE_IJVectorCreate(mpiComm_, procNRows[mypid],
2374 procNRows[mypid+1]-1, &R);
2375 ierr += HYPRE_IJVectorSetObjectType(R, HYPRE_PARCSR);
2376 ierr += HYPRE_IJVectorInitialize(R);
2377 ierr += HYPRE_IJVectorAssemble(R);
2378 hypre_assert(!ierr);
2379 HYPRE_IJVectorGetObject(R, (void **) &R_csr);
2380 HYPRE_ParVectorCopy( b_csr, R_csr );
2381 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, R_csr );
2382 HYPRE_ParVectorInnerProd( R_csr, R_csr, &rnorm);
2383 hypre_Vector *R_local=hypre_ParVectorLocalVector((hypre_ParVector*) R_csr);
2384 double *R_data = (double *) hypre_VectorData(R_local);
2385 HYPRE_ParVectorInnerProd( R_csr, R_csr, &rnorm);
2386 double rnorm2 = 0.0;
2387 for ( irow = 0; irow < nConstraints; irow++ )
2388 {
2389 searchIndex = -1;
2390 for ( jcol = 0; jcol < nConstraints; jcol++ )
2391 {
2392 if ( slaveEqnListAux_[jcol] == irow )
2393 {
2394 searchIndex = slaveEqnList_[jcol];
2395 break;
2396 }
2397 }
2398 rnorm2 += (R_data[searchIndex-startRow] * R_data[searchIndex-startRow]);
2399 }
2400 for ( irow = endRow-nConstraints+1; irow <= endRow; irow++ )
2401 rnorm2 += (R_data[irow-startRow] * R_data[irow-startRow]);
2402 HYPRE_IJVectorDestroy(R);
2403 if ( mypid == 0 )
2404 printf("HYPRE_SlideRedction norm check = %e %e %e\n", sqrt(rnorm),
2405 sqrt(rnorm-rnorm2), sqrt(rnorm2));
2406 #endif
2407
2408 //----------------------------------------------------------------
2409 // clean up
2410 //----------------------------------------------------------------
2411
2412 HYPRE_IJVectorDestroy(v1);
2413 HYPRE_IJVectorDestroy(x2);
2414 free( procNRows );
2415 return 0;
2416 }
2417
2418 //****************************************************************************
2419 // build A21 matrix
2420 //----------------------------------------------------------------------------
2421
buildA21Mat()2422 int HYPRE_SlideReduction::buildA21Mat()
2423 {
2424 int mypid, nprocs, *procNRows, startRow, endRow, localNRows;
2425 int globalNConstr, globalNRows, nConstraints, A21NRows, A21NCols;
2426 int A21GlobalNRows, A21GlobalNCols, A21StartRow, A21StartCol, ierr;
2427 int rowCount, maxRowSize, newEndRow, *A21MatSize, irow, is, rowIndex;
2428 int rowSize, *colInd, newRowSize, jcol, colIndex, searchIndex;
2429 int nnzA21, *newColInd, procIndex, ncnt, newColIndex;
2430 double *colVal, *newColVal;
2431 char fname[40];
2432 HYPRE_ParCSRMatrix A_csr, A21_csr;
2433
2434 //------------------------------------------------------------------
2435 // get matrix information
2436 //------------------------------------------------------------------
2437
2438 MPI_Comm_rank(mpiComm_, &mypid);
2439 MPI_Comm_size(mpiComm_, &nprocs);
2440 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
2441 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
2442 startRow = procNRows[mypid];
2443 endRow = procNRows[mypid+1] - 1;
2444 localNRows = endRow - startRow + 1;
2445 globalNConstr = procNConstr_[nprocs];
2446 globalNRows = procNRows[nprocs];
2447 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
2448
2449 //******************************************************************
2450 // extract A21 from A
2451 //------------------------------------------------------------------
2452 // calculate the dimension of A21
2453 //------------------------------------------------------------------
2454
2455 A21NRows = 2 * nConstraints;
2456 A21NCols = localNRows - nConstraints;
2457 A21GlobalNRows = 2 * globalNConstr;
2458 A21GlobalNCols = globalNRows - globalNConstr;
2459 A21StartRow = 2 * procNConstr_[mypid];
2460 A21StartCol = procNRows[mypid] - procNConstr_[mypid];
2461
2462 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
2463 {
2464 printf("%4d : buildA21Mat - A21StartRow = %d\n", mypid, A21StartRow);
2465 printf("%4d : buildA21Mat - A21GlobalDim = %d %d\n", mypid,
2466 A21GlobalNRows, A21GlobalNCols);
2467 printf("%4d : buildA21Mat - A21LocalDim = %d %d\n",mypid,
2468 A21NRows, A21NCols);
2469 }
2470
2471 //------------------------------------------------------------------
2472 // create a matrix context for A21
2473 //------------------------------------------------------------------
2474
2475 ierr = HYPRE_IJMatrixCreate(mpiComm_,A21StartRow,A21StartRow+A21NRows-1,
2476 A21StartCol,A21StartCol+A21NCols-1,&A21mat_);
2477 ierr += HYPRE_IJMatrixSetObjectType(A21mat_, HYPRE_PARCSR);
2478 hypre_assert(!ierr);
2479
2480 //------------------------------------------------------------------
2481 // compute the number of nonzeros in the first nConstraint row of A21
2482 // (which consists of the rows in selectedList), the nnz will
2483 // be reduced by excluding the constraint and selected slave columns
2484 //------------------------------------------------------------------
2485
2486 rowCount = maxRowSize = 0;
2487 newEndRow = endRow - nConstraints;
2488 if ( A21NRows > 0 ) A21MatSize = new int[A21NRows];
2489 else A21MatSize = NULL;
2490 for ( irow = 0; irow < nConstraints; irow++ )
2491 {
2492 for ( is = 0; is < nConstraints; is++ )
2493 {
2494 if ( slaveEqnListAux_[is] == irow )
2495 {
2496 rowIndex = slaveEqnList_[is];
2497 break;
2498 }
2499 }
2500 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2501 newRowSize = 0;
2502 for ( jcol = 0; jcol < rowSize; jcol++ )
2503 {
2504 colIndex = colInd[jcol];
2505 if ( colVal[jcol] != 0.0 )
2506 {
2507 if ( colIndex <= newEndRow || colIndex > endRow )
2508 {
2509 searchIndex = hypre_BinarySearch(gSlaveEqnList_,colIndex,
2510 globalNConstr);
2511 if ( searchIndex < 0 ) newRowSize++;
2512 }
2513 }
2514 }
2515 A21MatSize[irow] = newRowSize;
2516 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
2517 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2518 }
2519
2520 //------------------------------------------------------------------
2521 // compute the number of nonzeros in the second nConstraint row of A21
2522 // (which consists of the rows in constraint equations)
2523 //------------------------------------------------------------------
2524
2525 rowCount = nConstraints;
2526 for ( irow = endRow-nConstraints+1; irow <= endRow; irow++ )
2527 {
2528 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
2529 newRowSize = 0;
2530 for ( jcol = 0; jcol < rowSize; jcol++ )
2531 {
2532 if ( colVal[jcol] != 0.0 )
2533 {
2534 colIndex = colInd[jcol];
2535 if ( colIndex <= newEndRow || colIndex > endRow )
2536 {
2537 searchIndex = hypre_BinarySearch(gSlaveEqnList_,colIndex,
2538 globalNConstr);
2539 if ( searchIndex < 0 ) newRowSize++;
2540 }
2541 }
2542 }
2543 A21MatSize[rowCount] = newRowSize;
2544 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
2545 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
2546 rowCount++;
2547 }
2548 nnzA21 = 0;
2549 for ( irow = 0; irow < 2*nConstraints; irow++ ) nnzA21 += A21MatSize[irow];
2550 MPI_Allreduce(&nnzA21,&ncnt,1,MPI_INT,MPI_SUM,mpiComm_);
2551 if ( mypid == 0 && ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
2552 printf(" 0 : buildA21Mat : NNZ of A21 = %d\n", ncnt);
2553
2554 //------------------------------------------------------------------
2555 // after fetching the row sizes, set up A21 with such sizes
2556 //------------------------------------------------------------------
2557
2558 ierr = HYPRE_IJMatrixSetRowSizes(A21mat_, A21MatSize);
2559 ierr += HYPRE_IJMatrixInitialize(A21mat_);
2560 hypre_assert(!ierr);
2561 if ( A21NRows > 0 ) delete [] A21MatSize;
2562
2563 //------------------------------------------------------------------
2564 // next load the first nConstraint row to A21 extracted from A
2565 // (at the same time, the D block is saved for future use)
2566 //------------------------------------------------------------------
2567
2568 rowCount = A21StartRow;
2569 newColInd = new int[maxRowSize+1];
2570 newColVal = new double[maxRowSize+1];
2571
2572 for ( irow = 0; irow < nConstraints; irow++ )
2573 {
2574 for ( is = 0; is < nConstraints; is++ )
2575 {
2576 if ( slaveEqnListAux_[is] == irow )
2577 {
2578 rowIndex = slaveEqnList_[is];
2579 break;
2580 }
2581 }
2582 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2583 newRowSize = 0;
2584 for ( jcol = 0; jcol < rowSize; jcol++ )
2585 {
2586 if ( colVal[jcol] != 0.0 )
2587 {
2588 colIndex = colInd[jcol];
2589 if ( colIndex <= newEndRow || colIndex > endRow )
2590 {
2591 searchIndex = HYPRE_LSI_Search(gSlaveEqnList_,colIndex,
2592 globalNConstr);
2593 if ( searchIndex < 0 )
2594 {
2595 for ( procIndex = 0; procIndex < nprocs; procIndex++ )
2596 if ( procNRows[procIndex] > colIndex ) break;
2597 procIndex--;
2598 newColIndex = colIndex - procNConstr_[procIndex];
2599 newColInd[newRowSize] = newColIndex;
2600 newColVal[newRowSize++] = colVal[jcol];
2601 if ( newColIndex < 0 || newColIndex >= A21GlobalNCols )
2602 {
2603 printf("%4d : buildA21Mat ERROR - ",mypid);
2604 printf(" out of range (%d,%d (%d))\n", rowCount,
2605 colIndex, A21GlobalNCols);
2606 for ( is = 0; is < rowSize; is++ )
2607 printf("%4d : row %7d has col = %7d\n",mypid,rowIndex,
2608 colInd[is]);
2609 exit(1);
2610 }
2611 if ( newRowSize > maxRowSize+1 )
2612 {
2613 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 2 )
2614 {
2615 printf("%4d : buildA21Mat WARNING - ",mypid);
2616 printf("passing array boundary(1).\n");
2617 }
2618 }
2619 }
2620 }
2621 }
2622 }
2623 HYPRE_IJMatrixSetValues(A21mat_, 1, &newRowSize, (const int *) &rowCount,
2624 (const int *) newColInd, (const double *) newColVal);
2625 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2626 rowCount++;
2627 }
2628
2629 //------------------------------------------------------------------
2630 // next load the second nConstraint rows to A21 extracted from A
2631 //------------------------------------------------------------------
2632
2633 for ( irow = endRow-nConstraints+1; irow <= endRow; irow++ )
2634 {
2635 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
2636 newRowSize = 0;
2637 for ( jcol = 0; jcol < rowSize; jcol++ )
2638 {
2639 colIndex = colInd[jcol];
2640 if (colVal[jcol] != 0.0 &&
2641 (colIndex <= newEndRow || colIndex > endRow))
2642 {
2643 searchIndex = hypre_BinarySearch(gSlaveEqnList_, colIndex,
2644 globalNConstr);
2645 if ( searchIndex < 0 )
2646 {
2647 for ( procIndex = 0; procIndex < nprocs; procIndex++ )
2648 if ( procNRows[procIndex] > colIndex ) break;
2649 procIndex--;
2650 newColIndex = colInd[jcol] - procNConstr_[procIndex];
2651 newColInd[newRowSize] = newColIndex;
2652 newColVal[newRowSize++] = colVal[jcol];
2653 }
2654 }
2655 }
2656 HYPRE_IJMatrixSetValues(A21mat_, 1, &newRowSize, (const int *) &rowCount,
2657 (const int *) newColInd, (const double *) newColVal);
2658 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
2659 rowCount++;
2660 }
2661 delete [] newColInd;
2662 delete [] newColVal;
2663 free( procNRows );
2664
2665 //------------------------------------------------------------------
2666 // finally assemble the matrix and sanitize
2667 //------------------------------------------------------------------
2668
2669 HYPRE_IJMatrixAssemble(A21mat_);
2670 HYPRE_IJMatrixGetObject(A21mat_, (void **) &A21_csr);
2671 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) A21_csr);
2672
2673 if ( outputLevel_ >= 5 )
2674 {
2675 sprintf(fname, "A21.%d", mypid);
2676 FILE *fp = fopen(fname, "w");
2677
2678 if ( mypid == 0 )
2679 {
2680 printf("====================================================\n");
2681 printf("%4d : Printing A21 matrix... \n", mypid);
2682 fflush(stdout);
2683 }
2684 for (irow = A21StartRow;irow < A21StartRow+2*nConstraints;irow++)
2685 {
2686 HYPRE_ParCSRMatrixGetRow(A21_csr,irow,&rowSize,&colInd,&colVal);
2687 for ( jcol = 0; jcol < rowSize; jcol++ )
2688 if ( colVal[jcol] != 0.0 )
2689 fprintf(fp, "%6d %6d %25.16e \n",irow+1,colInd[jcol]+1,
2690 colVal[jcol]);
2691 HYPRE_ParCSRMatrixRestoreRow(A21_csr, irow, &rowSize, &colInd,
2692 &colVal);
2693 }
2694 fclose(fp);
2695 if ( mypid == 0 )
2696 printf("====================================================\n");
2697 }
2698 return 0;
2699 }
2700
2701 //****************************************************************************
2702 // build invA22 matrix
2703 // - given A22 = | B C |, compute | 0 C^{-T} |
2704 // | C^T 0 | | C^{-1} -C^{-1} B C^{-T} |
2705 //----------------------------------------------------------------------------
2706
buildInvA22Mat()2707 int HYPRE_SlideReduction::buildInvA22Mat()
2708 {
2709 int mypid, nprocs, *procNRows, endRow;
2710 int is, globalNConstr, nConstraints, irow, jcol, rowSize;
2711 int *colInd, newEndRow, rowIndex, colIndex, searchIndex;
2712 int ig, ir, ic, ierr, index, offset, newRowSize, *newColInd;
2713 int maxBlkSize=HYPRE_SLIDEMAX, procIndex;
2714 int nGroups, *groupIDs, **groupRowNums, *groupSizes, *iTempList;
2715 int rowCount, *colInd2, rowSize2;
2716 double *colVal, *colVal2, **Imat, **Imat2, *newColVal;
2717 char fname[40];
2718 HYPRE_ParCSRMatrix A_csr, invA22_csr;
2719
2720 //------------------------------------------------------------------
2721 // get matrix information
2722 //------------------------------------------------------------------
2723
2724 MPI_Comm_rank(mpiComm_, &mypid);
2725 MPI_Comm_size(mpiComm_, &nprocs);
2726 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
2727 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
2728 endRow = procNRows[mypid+1] - 1;
2729 globalNConstr = procNConstr_[nprocs];
2730 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
2731 newEndRow = endRow - nConstraints;
2732
2733 //------------------------------------------------------------------
2734 // construct the group information
2735 //------------------------------------------------------------------
2736
2737 nGroups = 0;
2738 if ( nConstraints > 0 )
2739 {
2740 iTempList = new int[nConstraints];
2741 for ( irow = 0; irow < nConstraints; irow++ )
2742 iTempList[irow] = constrBlkInfo_[irow];
2743 hypre_qsort0( iTempList, 0, nConstraints-1 );
2744 nGroups = 1;
2745 for ( irow = 1; irow < nConstraints; irow++ )
2746 if ( iTempList[irow] != iTempList[irow-1] ) nGroups++;
2747 groupIDs = new int[nGroups];
2748 groupSizes = new int[nGroups];
2749 groupIDs[0] = iTempList[0];
2750 groupSizes[0] = 1;
2751 nGroups = 1;
2752 for ( irow = 1; irow < nConstraints; irow++ )
2753 {
2754 if ( iTempList[irow] != iTempList[irow-1] )
2755 {
2756 groupSizes[nGroups] = 1;
2757 groupIDs[nGroups++] = iTempList[irow];
2758 }
2759 else groupSizes[nGroups-1]++;
2760 }
2761 groupRowNums = new int*[nGroups];
2762 for ( ig = 0; ig < nGroups; ig++ )
2763 {
2764 if ( groupSizes[ig] > maxBlkSize )
2765 {
2766 printf("%4d : buildInvA22 ERROR - block Size %d >= %d\n", mypid,
2767 groupSizes[ig], maxBlkSize);
2768 printf("%4d : buildInvA22 ERROR - group ID = %d\n", mypid,
2769 groupIDs[ig]);
2770 exit(1);
2771 }
2772 }
2773 for ( ig = 0; ig < nGroups; ig++ )
2774 {
2775 groupRowNums[ig] = new int[groupSizes[ig]];
2776 groupSizes[ig] = 0;
2777 }
2778 for ( irow = 0; irow < nConstraints; irow++ )
2779 {
2780 index = constrBlkInfo_[irow];
2781 searchIndex = hypre_BinarySearch(groupIDs, index, nGroups);
2782 groupRowNums[searchIndex][groupSizes[searchIndex]++] = irow;
2783 }
2784 delete [] iTempList;
2785 }
2786
2787 //------------------------------------------------------------------
2788 // first extract the (2,1) block of A22
2789 // ( constraints-to-local slaves )
2790 //------------------------------------------------------------------
2791
2792 int *CT_JA = NULL, CTRowSize, CTOffset;
2793 double *CT_AA = NULL;
2794 if ( nConstraints > 0 )
2795 {
2796 CT_JA = new int[nConstraints*maxBlkSize];
2797 CT_AA = new double[nConstraints*maxBlkSize];
2798 for ( irow = 0; irow < nConstraints*maxBlkSize; irow++ ) CT_JA[irow] = -1;
2799 }
2800 #if 0
2801 FILE *fp = fopen("CT.m","w");
2802 #endif
2803 for ( irow = 0; irow < nConstraints; irow++ )
2804 {
2805 rowIndex = newEndRow + 1 + irow;
2806 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2807 CTOffset = maxBlkSize * irow;
2808 CTRowSize = 0;
2809 for ( jcol = 0; jcol < rowSize; jcol++ )
2810 {
2811 colIndex = colInd[jcol];
2812 searchIndex = hypre_BinarySearch(slaveEqnList_,colIndex,nConstraints);
2813 if ( searchIndex >= 0 )
2814 {
2815 CT_JA[CTOffset+CTRowSize] = slaveEqnListAux_[searchIndex];
2816 CT_AA[CTOffset+CTRowSize] = colVal[jcol];
2817 CTRowSize++;
2818 #if 0
2819 fprintf(fp,"%d %d %25.16e\n",irow+1,CT_JA[CTOffset+CTRowSize-1]+1,colVal[jcol]);
2820 #endif
2821 }
2822 }
2823 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2824 }
2825 #if 0
2826 fclose(fp);
2827 #endif
2828
2829 //------------------------------------------------------------------
2830 // invert the (2,1) block of A22
2831 //------------------------------------------------------------------
2832
2833 #if 0
2834 FILE *fp2 = fopen("invCT.m","w");
2835 #endif
2836 Imat = hypre_TAlloc(double*, maxBlkSize , HYPRE_MEMORY_HOST);
2837 for ( ir = 0; ir < maxBlkSize; ir++ )
2838 Imat[ir] = hypre_TAlloc(double, maxBlkSize , HYPRE_MEMORY_HOST);
2839
2840 for ( ig = 0; ig < nGroups; ig++ )
2841 {
2842 for ( ir = 0; ir < groupSizes[ig]; ir++ )
2843 for ( ic = 0; ic < groupSizes[ig]; ic++ ) Imat[ir][ic] = 0.0;
2844 for ( ir = 0; ir < groupSizes[ig]; ir++ )
2845 {
2846 rowIndex = groupRowNums[ig][ir];
2847 offset = rowIndex * maxBlkSize;
2848 for ( ic = 0; ic < maxBlkSize; ic++ )
2849 {
2850 colIndex = CT_JA[offset+ic];
2851 if ( colIndex != -1 )
2852 {
2853 for ( is = 0; is < groupSizes[ig]; is++ )
2854 if ( colIndex == groupRowNums[ig][is] ) break;
2855 Imat[ir][is] = CT_AA[offset+ic];
2856 }
2857 }
2858 }
2859 ierr = HYPRE_LSI_MatrixInverse((double**) Imat, groupSizes[ig], &Imat2);
2860 if ( ierr )
2861 {
2862 printf("Failed Block %d has indices (%d) : ", ig, groupSizes[ig]);
2863 for ( ir = 0; ir < groupSizes[ig]; ir++ )
2864 printf(" %d ", groupRowNums[ig][ir]);
2865 printf("\n");
2866 for ( ir = 0; ir < groupSizes[ig]; ir++ )
2867 {
2868 for ( ic = 0; ic < groupSizes[ig]; ic++ )
2869 printf(" %e ", Imat[ir][ic]);
2870 printf("\n");
2871 }
2872 printf("\n");
2873 for ( ir = 0; ir < groupSizes[ig]; ir++ )
2874 {
2875 for ( ic = 0; ic < groupSizes[ig]; ic++ )
2876 printf(" %e ", Imat2[ir][ic]);
2877 printf("\n");
2878 }
2879 printf("\n");
2880 for ( ir = 0; ir < groupSizes[ig]; ir++ )
2881 {
2882 rowIndex = groupRowNums[ig][ir];
2883 offset = rowIndex * maxBlkSize;
2884 for ( ic = 0; ic < maxBlkSize; ic++ )
2885 {
2886 colIndex = CT_JA[offset+ic];
2887 if ( colIndex != -1 )
2888 {
2889 printf(" rowIndex,colIndex,val = %d %d %e\n",rowIndex,
2890 colIndex,CT_AA[offset+ic]);
2891 }
2892 }
2893 }
2894 }
2895 hypre_assert( !ierr );
2896 for ( ir = 0; ir < groupSizes[ig]; ir++ )
2897 {
2898 rowIndex = groupRowNums[ig][ir];
2899 offset = rowIndex * maxBlkSize;
2900 for (ic = 0; ic < maxBlkSize; ic++) CT_JA[offset+ic] = -1;
2901 for ( ic = 0; ic < groupSizes[ig]; ic++ )
2902 {
2903 if ( Imat2[ir][ic] != 0.0 )
2904 {
2905 CT_JA[offset+ic] = groupRowNums[ig][ic];
2906 CT_AA[offset+ic] = Imat2[ir][ic];
2907 #if 0
2908 fprintf(fp2,"%d %d %25.16e\n",rowIndex+1,CT_JA[offset+ic]+1,CT_AA[offset+ic]);
2909 #endif
2910 }
2911 }
2912 free( Imat2[ir] );
2913 }
2914 free( Imat2 );
2915 }
2916 #if 0
2917 fclose(fp2);
2918 #endif
2919 for ( ir = 0; ir < maxBlkSize; ir++ ) free( Imat[ir] );
2920 free( Imat );
2921
2922 //------------------------------------------------------------------
2923 // form ParCSRMatrix of the (2,1) block of A22
2924 //------------------------------------------------------------------
2925
2926 int *hypreCTMatSize, maxRowSize;
2927 hypre_ParCSRMatrix *hypreCT;
2928 HYPRE_IJMatrix IJCT;
2929
2930 ierr = HYPRE_IJMatrixCreate(mpiComm_, procNConstr_[mypid],
2931 procNConstr_[mypid]+nConstraints-1, procNConstr_[mypid],
2932 procNConstr_[mypid]+nConstraints-1, &IJCT);
2933 ierr += HYPRE_IJMatrixSetObjectType(IJCT, HYPRE_PARCSR);
2934 hypre_assert(!ierr);
2935 if ( nConstraints > 0 ) hypreCTMatSize = new int[nConstraints];
2936 else hypreCTMatSize = NULL;
2937 maxRowSize = 0;
2938 for ( irow = 0; irow < nConstraints; irow++ )
2939 {
2940 newRowSize = 0;
2941 offset = irow * maxBlkSize;
2942 for ( ic = 0; ic < maxBlkSize; ic++ )
2943 if ( CT_JA[offset+ic] != -1 ) newRowSize++;
2944 hypreCTMatSize[irow] = newRowSize;
2945 maxRowSize = (newRowSize > maxRowSize) ? newRowSize : maxRowSize;
2946 }
2947 ierr = HYPRE_IJMatrixSetRowSizes(IJCT, hypreCTMatSize);
2948 ierr = HYPRE_IJMatrixInitialize(IJCT);
2949 hypre_assert(!ierr);
2950 if ( nConstraints > 0 ) delete [] hypreCTMatSize;
2951 newColInd = new int[maxRowSize+1];
2952 newColVal = new double[maxRowSize+1];
2953 for ( irow = 0; irow < nConstraints; irow++ )
2954 {
2955 rowIndex = procNConstr_[mypid] + irow;
2956 offset = irow * maxBlkSize;
2957 newRowSize = 0;
2958 for ( ic = 0; ic < maxBlkSize; ic++ )
2959 {
2960 if ( CT_JA[offset+ic] != -1 )
2961 {
2962 newColInd[newRowSize] = CT_JA[offset+ic] + procNConstr_[mypid];
2963 newColVal[newRowSize++] = CT_AA[offset+ic];
2964 }
2965 }
2966 ierr = HYPRE_IJMatrixSetValues(IJCT, 1, &newRowSize,
2967 (const int *) &rowIndex, (const int *) newColInd,
2968 (const double *) newColVal);
2969 hypre_assert( !ierr );
2970 }
2971 delete [] newColInd;
2972 delete [] newColVal;
2973 if ( nConstraints > 0 )
2974 {
2975 delete [] CT_JA;
2976 delete [] CT_AA;
2977 }
2978 HYPRE_IJMatrixAssemble(IJCT);
2979 HYPRE_IJMatrixGetObject(IJCT, (void **) &hypreCT);
2980 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) hypreCT);
2981
2982 //------------------------------------------------------------------
2983 // next extract the (1,2) block of A22
2984 // ( local slaves-to-constraints )
2985 //------------------------------------------------------------------
2986
2987 int *C_JA = NULL, CRowSize, COffset;
2988 double *C_AA = NULL;
2989 if ( nConstraints > 0 )
2990 {
2991 C_JA = new int[nConstraints*maxBlkSize];
2992 C_AA = new double[nConstraints*maxBlkSize];
2993 for ( irow = 0; irow < nConstraints*maxBlkSize; irow++ ) C_JA[irow] = -1;
2994 }
2995 for ( irow = 0; irow < nConstraints; irow++ )
2996 {
2997 for ( is = 0; is < nConstraints; is++ )
2998 {
2999 if ( slaveEqnListAux_[is] == irow )
3000 {
3001 rowIndex = slaveEqnList_[is];
3002 break;
3003 }
3004 }
3005 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
3006 CRowSize = 0;
3007 COffset = maxBlkSize * irow;
3008 for ( jcol = 0; jcol < rowSize; jcol++ )
3009 {
3010 colIndex = colInd[jcol];
3011 if ( colIndex > newEndRow && colIndex <= endRow )
3012 {
3013 C_JA[COffset+CRowSize] = colIndex - newEndRow - 1;
3014 C_AA[COffset+CRowSize] = colVal[jcol];
3015 CRowSize++;
3016 }
3017 }
3018 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
3019 }
3020
3021 //------------------------------------------------------------------
3022 // invert the (2,1) block of A22
3023 //------------------------------------------------------------------
3024
3025 Imat = hypre_TAlloc(double*, maxBlkSize , HYPRE_MEMORY_HOST);
3026 for ( ir = 0; ir < maxBlkSize; ir++ )
3027 Imat[ir] = hypre_TAlloc(double, maxBlkSize , HYPRE_MEMORY_HOST);
3028
3029 for ( ig = 0; ig < nGroups; ig++ )
3030 {
3031 for ( ir = 0; ir < groupSizes[ig]; ir++ )
3032 for ( ic = 0; ic < groupSizes[ig]; ic++ ) Imat[ir][ic] = 0.0;
3033 for ( ir = 0; ir < groupSizes[ig]; ir++ )
3034 {
3035 rowIndex = groupRowNums[ig][ir];
3036 offset = rowIndex * maxBlkSize;
3037 for ( ic = 0; ic < maxBlkSize; ic++ )
3038 {
3039 colIndex = C_JA[offset+ic];
3040 if ( colIndex != -1 )
3041 {
3042 for ( is = 0; is < groupSizes[ig]; is++ )
3043 if ( colIndex == groupRowNums[ig][is] ) break;
3044 Imat[ir][is] = C_AA[offset+ic];
3045 }
3046 }
3047 }
3048 ierr = HYPRE_LSI_MatrixInverse((double**) Imat, groupSizes[ig], &Imat2);
3049 hypre_assert( !ierr );
3050 for ( ir = 0; ir < groupSizes[ig]; ir++ )
3051 {
3052 rowIndex = groupRowNums[ig][ir];
3053 offset = rowIndex * maxBlkSize;
3054 for (ic = 0; ic < maxBlkSize; ic++) C_JA[offset+ic] = -1;
3055 for ( ic = 0; ic < groupSizes[ig]; ic++ )
3056 {
3057 if ( Imat2[ir][ic] != 0.0 )
3058 {
3059 C_JA[offset+ic] = groupRowNums[ig][ic];
3060 C_AA[offset+ic] = Imat2[ir][ic];
3061 }
3062 }
3063 free( Imat2[ir] );
3064 }
3065 free( Imat2 );
3066 }
3067 for ( ir = 0; ir < maxBlkSize; ir++ ) free( Imat[ir] );
3068 free( Imat );
3069
3070 //------------------------------------------------------------------
3071 // form ParCSRMatrix of the (1,2) block of A22
3072 //------------------------------------------------------------------
3073
3074 int *hypreCMatSize;
3075 hypre_ParCSRMatrix *hypreC;
3076 HYPRE_IJMatrix IJC;
3077
3078 ierr = HYPRE_IJMatrixCreate(mpiComm_, procNConstr_[mypid],
3079 procNConstr_[mypid]+nConstraints-1, procNConstr_[mypid],
3080 procNConstr_[mypid]+nConstraints-1, &IJC);
3081 ierr += HYPRE_IJMatrixSetObjectType(IJC, HYPRE_PARCSR);
3082 hypre_assert(!ierr);
3083 if ( nConstraints > 0 ) hypreCMatSize = new int[nConstraints];
3084 else hypreCMatSize = NULL;
3085 maxRowSize = 0;
3086 for ( irow = 0; irow < nConstraints; irow++ )
3087 {
3088 newRowSize = 0;
3089 offset = irow * maxBlkSize;
3090 for ( ic = 0; ic < maxBlkSize; ic++ )
3091 if ( C_JA[offset+ic] != -1 ) newRowSize++;
3092 hypreCMatSize[irow] = newRowSize;
3093 maxRowSize = (newRowSize > maxRowSize) ? newRowSize : maxRowSize;
3094 }
3095 ierr = HYPRE_IJMatrixSetRowSizes(IJC, hypreCMatSize);
3096 ierr = HYPRE_IJMatrixInitialize(IJC);
3097 hypre_assert(!ierr);
3098 if ( nConstraints > 0 ) delete [] hypreCMatSize;
3099 newColInd = new int[maxRowSize+1];
3100 newColVal = new double[maxRowSize+1];
3101 for ( irow = 0; irow < nConstraints; irow++ )
3102 {
3103 rowIndex = procNConstr_[mypid] + irow;
3104 offset = irow * maxBlkSize;
3105 newRowSize = 0;
3106 for ( ic = 0; ic < maxBlkSize; ic++ )
3107 {
3108 if ( C_JA[offset+ic] != -1 )
3109 {
3110 newColInd[newRowSize] = C_JA[offset+ic] + procNConstr_[mypid];
3111 newColVal[newRowSize++] = C_AA[offset+ic];
3112 }
3113 }
3114 ierr = HYPRE_IJMatrixSetValues(IJC, 1, &newRowSize,
3115 (const int *) &rowIndex, (const int *) newColInd,
3116 (const double *) newColVal);
3117 hypre_assert( !ierr );
3118 }
3119 delete [] newColInd;
3120 delete [] newColVal;
3121 if ( nConstraints > 0 )
3122 {
3123 delete [] C_JA;
3124 delete [] C_AA;
3125 }
3126 HYPRE_IJMatrixAssemble(IJC);
3127 HYPRE_IJMatrixGetObject(IJC, (void **) &hypreC);
3128 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) hypreC);
3129 if ( nConstraints > 0 )
3130 {
3131 delete [] groupIDs;
3132 delete [] groupSizes;
3133 for ( ig = 0; ig < nGroups; ig++ ) delete [] groupRowNums[ig];
3134 delete [] groupRowNums;
3135 }
3136
3137 //------------------------------------------------------------------
3138 // form ParCSRMatrix of the (2,2) block of the invA22 matrix
3139 //------------------------------------------------------------------
3140
3141 int *hypreBMatSize=NULL;
3142 hypre_ParCSRMatrix *hypreB;
3143 HYPRE_IJMatrix IJB;
3144
3145 ierr = HYPRE_IJMatrixCreate(mpiComm_, procNConstr_[mypid],
3146 procNConstr_[mypid]+nConstraints-1, procNConstr_[mypid],
3147 procNConstr_[mypid]+nConstraints-1, &IJB);
3148 ierr = HYPRE_IJMatrixSetObjectType(IJB, HYPRE_PARCSR);
3149 hypre_assert(!ierr);
3150 if ( nConstraints > 0 ) hypreBMatSize = new int[nConstraints];
3151 maxRowSize = 0;
3152 for ( irow = 0; irow < nConstraints; irow++ )
3153 {
3154 for ( is = 0; is < nConstraints; is++ )
3155 {
3156 if ( slaveEqnListAux_[is] == irow )
3157 {
3158 rowIndex = slaveEqnList_[is];
3159 break;
3160 }
3161 }
3162 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
3163 newRowSize = 0;
3164 for ( jcol = 0; jcol < rowSize; jcol++ )
3165 {
3166 colIndex = colInd[jcol];
3167 searchIndex = hypre_BinarySearch(gSlaveEqnList_, colIndex,
3168 globalNConstr);
3169 if ( searchIndex >= 0 ) newRowSize++;
3170 }
3171 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
3172 hypreBMatSize[irow] = newRowSize;
3173 maxRowSize = (newRowSize > maxRowSize) ? newRowSize : maxRowSize;
3174 }
3175 ierr = HYPRE_IJMatrixSetRowSizes(IJB, hypreBMatSize);
3176 ierr = HYPRE_IJMatrixInitialize(IJB);
3177 hypre_assert(!ierr);
3178 if ( nConstraints > 0 ) delete [] hypreBMatSize;
3179
3180 newColInd = new int[maxRowSize+1];
3181 newColVal = new double[maxRowSize+1];
3182
3183 for ( irow = 0; irow < nConstraints; irow++ )
3184 {
3185 for ( is = 0; is < nConstraints; is++ )
3186 {
3187 if ( slaveEqnListAux_[is] == irow )
3188 {
3189 rowIndex = slaveEqnList_[is];
3190 break;
3191 }
3192 }
3193 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
3194 newRowSize = 0;
3195 for ( jcol = 0; jcol < rowSize; jcol++ )
3196 {
3197 colIndex = colInd[jcol];
3198 searchIndex = hypre_BinarySearch(gSlaveEqnList_, colIndex,
3199 globalNConstr);
3200 if ( searchIndex >= 0 )
3201 {
3202 newColInd[newRowSize] = gSlaveEqnListAux_[searchIndex];
3203 newColVal[newRowSize++] = - colVal[jcol];
3204 }
3205 }
3206 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
3207 rowIndex = procNConstr_[mypid] + irow;
3208 ierr = HYPRE_IJMatrixSetValues(IJB, 1, &newRowSize,
3209 (const int *) &rowIndex, (const int *) newColInd,
3210 (const double *) newColVal);
3211 hypre_assert( !ierr );
3212 }
3213 HYPRE_IJMatrixAssemble(IJB);
3214 HYPRE_IJMatrixGetObject(IJB, (void **) &hypreB);
3215 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) hypreB);
3216 delete [] newColInd;
3217 delete [] newColVal;
3218
3219 //------------------------------------------------------------------
3220 // perform triple matrix product - C^{-1} B C^{-T}
3221 //------------------------------------------------------------------
3222
3223 HYPRE_ParCSRMatrix hypreCBC;
3224
3225 #if 0
3226 strcpy( fname, "hypreCT" );
3227 HYPRE_ParCSRMatrixPrint((HYPRE_ParCSRMatrix) hypreCT, fname);
3228 strcpy( fname, "hypreB" );
3229 HYPRE_ParCSRMatrixPrint((HYPRE_ParCSRMatrix) hypreB, fname);
3230 #endif
3231
3232 hypre_BoomerAMGBuildCoarseOperator(hypreCT, hypreB, hypreCT,
3233 (hypre_ParCSRMatrix **) &hypreCBC);
3234 #if 0
3235 for ( irow = 0; irow < nConstraints; irow++ )
3236 {
3237 HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCT,irow,&rowSize,
3238 &colInd,&colVal);
3239 for (jcol = 0; jcol < rowSize; jcol++ )
3240 printf("CT : %5d %5d %25.16e\n",irow+1,colInd[jcol]+1,colVal[jcol]);
3241 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCT,irow,&rowSize,
3242 &colInd,&colVal);
3243 }
3244 for ( irow = 0; irow < nConstraints; irow++ )
3245 {
3246 HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreB,irow,&rowSize,
3247 &colInd,&colVal);
3248 for (jcol = 0; jcol < rowSize; jcol++ )
3249 printf("B : %5d %5d %25.16e\n",irow+1,colInd[jcol]+1,colVal[jcol]);
3250 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreB,irow,&rowSize,
3251 &colInd,&colVal);
3252 }
3253 for ( irow = 0; irow < nConstraints; irow++ )
3254 {
3255 HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
3256 &colInd,&colVal);
3257 for (jcol = 0; jcol < rowSize; jcol++ )
3258 printf("CBC : %5d %5d %25.16e\n",irow+1,colInd[jcol]+1,colVal[jcol]);
3259 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
3260 &colInd,&colVal);
3261 }
3262 #endif
3263
3264 HYPRE_IJMatrixDestroy( IJB );
3265
3266 //------------------------------------------------------------------
3267 // calculate the dimension of invA22
3268 //------------------------------------------------------------------
3269
3270 int invA22NRows = 2 * nConstraints;
3271 int invA22NCols = invA22NRows;
3272 int invA22StartRow = 2 * procNConstr_[mypid];
3273 int invA22StartCol = invA22StartRow;
3274 int *invA22MatSize=NULL;
3275
3276 //------------------------------------------------------------------
3277 // create a matrix context for A22
3278 //------------------------------------------------------------------
3279
3280 ierr = HYPRE_IJMatrixCreate(mpiComm_, invA22StartRow,
3281 invA22StartRow+invA22NRows-1, invA22StartCol,
3282 invA22StartCol+invA22NCols-1, &invA22mat_);
3283 ierr += HYPRE_IJMatrixSetObjectType(invA22mat_, HYPRE_PARCSR);
3284 hypre_assert(!ierr);
3285
3286 //------------------------------------------------------------------
3287 // compute the no. of nonzeros in the first nConstraint row of invA22
3288 //------------------------------------------------------------------
3289
3290 maxRowSize = 0;
3291 if ( invA22NRows > 0 ) invA22MatSize = new int[invA22NRows];
3292 for ( irow = 0; irow < nConstraints; irow++ )
3293 {
3294 rowIndex = procNConstr_[mypid] + irow;
3295 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCT,rowIndex,
3296 &rowSize,NULL,NULL);
3297 hypre_assert( !ierr );
3298 invA22MatSize[irow] = rowSize;
3299 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCT,rowIndex,
3300 &rowSize,NULL,NULL);
3301 maxRowSize = ( rowSize > maxRowSize ) ? rowSize : maxRowSize;
3302 }
3303
3304 //------------------------------------------------------------------
3305 // compute the number of nonzeros in the second nConstraints row of
3306 // invA22 (consisting of [D and A22 block])
3307 //------------------------------------------------------------------
3308
3309 #if 0
3310 for ( irow = 0; irow < nConstraints; irow++ )
3311 {
3312 HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
3313 &colInd,&colVal);
3314 for (jcol = 0; jcol < rowSize; jcol++ )
3315 printf("CBC1 : %5d %5d %25.16e\n",irow+1,colInd[jcol]+1,colVal[jcol]);
3316 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
3317 &colInd,&colVal);
3318 }
3319 #endif
3320 for ( irow = 0; irow < nConstraints; irow++ )
3321 {
3322 rowIndex = procNConstr_[mypid] + irow;
3323 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreC,rowIndex,
3324 &rowSize,&colInd,&colVal);
3325 hypre_assert( !ierr );
3326 newRowSize = rowSize;
3327 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreC,rowIndex,
3328 &rowSize,&colInd,&colVal);
3329 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCBC,rowIndex,
3330 &rowSize,&colInd,&colVal);
3331 hypre_assert( !ierr );
3332 newRowSize += rowSize;
3333 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCBC,rowIndex,
3334 &rowSize,&colInd,&colVal);
3335 invA22MatSize[nConstraints+irow] = newRowSize;
3336 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
3337 }
3338
3339 //------------------------------------------------------------------
3340 // after fetching the row sizes, set up invA22 with such sizes
3341 //------------------------------------------------------------------
3342
3343 #if 0
3344 for ( irow = 0; irow < nConstraints; irow++ )
3345 {
3346 HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
3347 &colInd,&colVal);
3348 for (jcol = 0; jcol < rowSize; jcol++ )
3349 printf("CBC2 : %5d %5d %25.16e\n",irow+1,colInd[jcol]+1,colVal[jcol]);
3350 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
3351 &colInd,&colVal);
3352 }
3353 #endif
3354 ierr = HYPRE_IJMatrixSetRowSizes(invA22mat_, invA22MatSize);
3355 ierr += HYPRE_IJMatrixInitialize(invA22mat_);
3356 hypre_assert(!ierr);
3357 if ( invA22NRows > 0 ) delete [] invA22MatSize;
3358
3359 //------------------------------------------------------------------
3360 // next load the first nConstraints row to invA22 extracted from A
3361 // (that is, the D block)
3362 //------------------------------------------------------------------
3363
3364 newColInd = new int[maxRowSize+1];
3365 newColVal = new double[maxRowSize+1];
3366 for ( irow = 0; irow < nConstraints; irow++ )
3367 {
3368 rowIndex = procNConstr_[mypid] + irow;
3369 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCT,rowIndex,
3370 &rowSize,&colInd,&colVal);
3371 hypre_assert( !ierr );
3372 newRowSize = 0;
3373 for ( jcol = 0; jcol < rowSize; jcol++ )
3374 {
3375 newColInd[newRowSize] = colInd[jcol] + procNConstr_[mypid] +
3376 nConstraints;
3377 newColVal[newRowSize++] = colVal[jcol];
3378 }
3379 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCT,rowIndex,
3380 &rowSize,&colInd,&colVal);
3381 rowCount = invA22StartRow + irow;
3382 ierr = HYPRE_IJMatrixSetValues(invA22mat_, 1, &rowSize,
3383 (const int *) &rowCount, (const int *) newColInd,
3384 (const double *) newColVal);
3385 hypre_assert(!ierr);
3386 }
3387
3388 //------------------------------------------------------------------
3389 // next load the second nConstraints rows to A22 extracted from A
3390 //------------------------------------------------------------------
3391
3392 for ( irow = 0; irow < nConstraints; irow++ )
3393 {
3394 rowIndex = procNConstr_[mypid] + irow;
3395 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreC,rowIndex,
3396 &rowSize,&colInd,&colVal);
3397 hypre_assert( !ierr );
3398 newRowSize = 0;
3399 for ( jcol = 0; jcol < rowSize; jcol++ )
3400 {
3401 newColInd[newRowSize] = colInd[jcol] + procNConstr_[mypid];
3402 newColVal[newRowSize++] = colVal[jcol];
3403 }
3404 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreC,rowIndex,
3405 &rowSize,&colInd,&colVal);
3406 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCBC,rowIndex,
3407 &rowSize2,&colInd2,&colVal2);
3408 hypre_assert( !ierr );
3409 for ( jcol = 0; jcol < rowSize2; jcol++ )
3410 {
3411 colIndex = colInd2[jcol];
3412 for ( procIndex = 0; procIndex <= nprocs; procIndex++ )
3413 if ( procNConstr_[procIndex] > colIndex ) break;
3414 newColInd[newRowSize] = colIndex + procNConstr_[procIndex];
3415 newColVal[newRowSize++] = colVal2[jcol];
3416 }
3417 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCBC,rowIndex,
3418 &rowSize2,&colInd2,&colVal2);
3419 rowCount = invA22StartRow + nConstraints + irow;
3420 ierr = HYPRE_IJMatrixSetValues(invA22mat_, 1, &newRowSize,
3421 (const int *) &rowCount, (const int *) newColInd,
3422 (const double *) newColVal);
3423 hypre_assert(!ierr);
3424 }
3425 delete [] newColInd;
3426 delete [] newColVal;
3427 free( procNRows );
3428 HYPRE_IJMatrixDestroy( IJC );
3429 HYPRE_IJMatrixDestroy( IJCT );
3430 HYPRE_ParCSRMatrixDestroy( (HYPRE_ParCSRMatrix) hypreCBC );
3431
3432 //------------------------------------------------------------------
3433 // finally assemble the matrix and sanitize
3434 //------------------------------------------------------------------
3435
3436 HYPRE_IJMatrixAssemble(invA22mat_);
3437 HYPRE_IJMatrixGetObject(invA22mat_, (void **) &invA22_csr);
3438 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) invA22_csr);
3439
3440 if ( outputLevel_ >= 5 )
3441 {
3442 sprintf( fname, "invA.%d", mypid );
3443 FILE *fp = fopen( fname, "w");
3444
3445 if ( mypid == 0 )
3446 {
3447 printf("====================================================\n");
3448 printf("%4d : Printing invA22 matrix... \n", mypid);
3449 fflush(stdout);
3450 }
3451 for (irow=invA22StartRow; irow < invA22StartRow+invA22NRows;irow++)
3452 {
3453 HYPRE_ParCSRMatrixGetRow(invA22_csr,irow,&rowSize,&colInd,&colVal);
3454 for ( jcol = 0; jcol < rowSize; jcol++ )
3455 if ( colVal[jcol] != 0.0 )
3456 fprintf(fp,"%6d %6d %25.16e \n",irow+1,colInd[jcol]+1,
3457 colVal[jcol]);
3458 HYPRE_ParCSRMatrixRestoreRow(invA22_csr,irow,&rowSize,&colInd,
3459 &colVal);
3460 }
3461 fclose(fp);
3462 if ( mypid == 0 )
3463 printf("====================================================\n");
3464 }
3465 return 0;
3466 }
3467
3468 //***************************************************************************
3469 // scale matrix
3470 //---------------------------------------------------------------------------
3471
scaleMatrixVector()3472 int HYPRE_SlideReduction::scaleMatrixVector()
3473 {
3474 int *partition, startRow, localNRows, index, offset;
3475 int irow, jcol, iP, rowSize, *colInd, *rowLengs;
3476 int mypid, nprocs;
3477 int nSends, *sendStarts, *sendMap, *offdMap, ierr;
3478 int nRecvs, *recvStarts, pstart, pend, maxRowLeng;
3479 int *ADiagI, *ADiagJ, *AOffdI, *AOffdJ, rowInd;
3480 double *ADiagA, *AOffdA, *bData, *b2Data;
3481 double *scaleVec, *extScaleVec, *colVal, *sBuffer;
3482 HYPRE_IJMatrix newA;
3483 HYPRE_IJVector newB;
3484 hypre_ParCSRMatrix *A_csr;
3485 hypre_CSRMatrix *ADiag, *AOffd;
3486 hypre_ParVector *b_csr, *b2_csr;
3487 hypre_ParCSRCommPkg *commPkg;
3488 hypre_ParCSRCommHandle *commHandle;
3489
3490 //-----------------------------------------------------------------------
3491 // fetch matrix and parameters
3492 //-----------------------------------------------------------------------
3493
3494 MPI_Comm_rank( mpiComm_, &mypid );
3495 MPI_Comm_size( mpiComm_, &nprocs );
3496 HYPRE_IJMatrixGetObject(reducedAmat_, (void **) &A_csr);
3497 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) A_csr);
3498 HYPRE_ParCSRMatrixGetRowPartitioning((HYPRE_ParCSRMatrix) A_csr,&partition);
3499 startRow = partition[mypid];
3500 localNRows = partition[mypid+1] - startRow;
3501 free( partition );
3502 ADiag = hypre_ParCSRMatrixDiag(A_csr);
3503 ADiagI = hypre_CSRMatrixI(ADiag);
3504 ADiagJ = hypre_CSRMatrixJ(ADiag);
3505 ADiagA = hypre_CSRMatrixData(ADiag);
3506 AOffd = hypre_ParCSRMatrixOffd(A_csr);
3507 AOffdI = hypre_CSRMatrixI(AOffd);
3508 AOffdJ = hypre_CSRMatrixJ(AOffd);
3509 AOffdA = hypre_CSRMatrixData(AOffd);
3510 HYPRE_IJVectorGetObject(reducedBvec_, (void **) &b_csr);
3511 bData = hypre_VectorData(hypre_ParVectorLocalVector(b_csr));
3512
3513 offdMap = hypre_ParCSRMatrixColMapOffd(A_csr);
3514 commPkg = hypre_ParCSRMatrixCommPkg((hypre_ParCSRMatrix *) A_csr);
3515 nSends = hypre_ParCSRCommPkgNumSends(commPkg);
3516 nRecvs = hypre_ParCSRCommPkgNumRecvs(commPkg);
3517 recvStarts = hypre_ParCSRCommPkgRecvVecStarts(commPkg);
3518 sendStarts = hypre_ParCSRCommPkgSendMapStarts(commPkg);
3519 sendMap = hypre_ParCSRCommPkgSendMapElmts(commPkg);
3520
3521 //-----------------------------------------------------------------------
3522 // fetch diagonal of A
3523 //-----------------------------------------------------------------------
3524
3525 scaleVec = new double[localNRows];
3526 rowLengs = new int[localNRows];
3527 extScaleVec = NULL;
3528 if ( nRecvs > 0 ) extScaleVec = new double[recvStarts[nRecvs]];
3529
3530 maxRowLeng = 0;
3531 for ( irow = 0; irow < localNRows; irow++ )
3532 {
3533 scaleVec[irow] = 0.0;
3534 rowLengs[irow] = ADiagI[irow+1] - ADiagI[irow] +
3535 AOffdI[irow+1] - AOffdI[irow];
3536 if ( rowLengs[irow] > maxRowLeng ) maxRowLeng = rowLengs[irow];
3537 for ( jcol = ADiagI[irow]; jcol < ADiagI[irow+1]; jcol++ )
3538 if ( ADiagJ[jcol] == irow ) scaleVec[irow] = ADiagA[jcol];
3539 }
3540 for ( irow = 0; irow < localNRows; irow++ )
3541 {
3542 if ( habs( scaleVec[irow] ) == 0.0 )
3543 {
3544 printf("%d : scaleMatrixVector - diag %d = %e <= 0 \n",mypid,irow,
3545 scaleVec[irow]);
3546 exit(1);
3547 }
3548 scaleVec[irow] = 1.0/sqrt(scaleVec[irow]);
3549 }
3550
3551 //-----------------------------------------------------------------------
3552 // exchange diagonal of A
3553 //-----------------------------------------------------------------------
3554
3555 if ( nSends > 0 )
3556 {
3557 sBuffer = new double[sendStarts[nSends]];
3558 offset = 0;
3559 for ( iP = 0; iP < nSends; iP++ )
3560 {
3561 pstart = sendStarts[iP];
3562 pend = sendStarts[iP+1];
3563 for ( jcol = pstart; jcol < pend; jcol++ )
3564 {
3565 index = sendMap[jcol];
3566 sBuffer[offset++] = scaleVec[index];
3567 }
3568 }
3569 }
3570 else sBuffer = NULL;
3571
3572 commHandle = hypre_ParCSRCommHandleCreate(1,commPkg,sBuffer,extScaleVec);
3573 hypre_ParCSRCommHandleDestroy(commHandle);
3574
3575 if ( nSends > 0 ) delete [] sBuffer;
3576
3577 //-----------------------------------------------------------------------
3578 // construct new matrix
3579 //-----------------------------------------------------------------------
3580
3581 HYPRE_IJMatrixCreate(mpiComm_, startRow, startRow+localNRows-1,
3582 startRow, startRow+localNRows-1, &newA);
3583 HYPRE_IJMatrixSetObjectType(newA, HYPRE_PARCSR);
3584 HYPRE_IJMatrixSetRowSizes(newA, rowLengs);
3585 HYPRE_IJMatrixInitialize(newA);
3586 delete [] rowLengs;
3587 colInd = new int[maxRowLeng];
3588 colVal = new double[maxRowLeng];
3589 for ( irow = 0; irow < localNRows; irow++ )
3590 {
3591 rowSize = 0;
3592 for ( jcol = ADiagI[irow]; jcol < ADiagI[irow+1]; jcol++ )
3593 {
3594 index = ADiagJ[jcol];
3595 colInd[rowSize] = index + startRow;
3596 colVal[rowSize++] = scaleVec[irow]*scaleVec[index]*ADiagA[jcol];
3597 }
3598 for ( jcol = AOffdI[irow]; jcol < AOffdI[irow+1]; jcol++ )
3599 {
3600 index = AOffdJ[jcol];
3601 colInd[rowSize] = offdMap[index];
3602 colVal[rowSize++] = scaleVec[irow]*extScaleVec[index]*AOffdA[jcol];
3603 }
3604 rowInd = irow + startRow;
3605 HYPRE_IJMatrixSetValues(newA, 1, &rowSize, (const int *) &rowInd,
3606 (const int *) colInd, (const double *) colVal);
3607 }
3608 HYPRE_IJMatrixAssemble(newA);
3609 delete [] colInd;
3610 delete [] colVal;
3611 delete [] extScaleVec;
3612
3613 //-----------------------------------------------------------------------
3614 // construct new vector
3615 //-----------------------------------------------------------------------
3616
3617 ierr = HYPRE_IJVectorCreate(mpiComm_,startRow,startRow+localNRows-1,&newB);
3618 ierr += HYPRE_IJVectorSetObjectType(newB, HYPRE_PARCSR);
3619 ierr += HYPRE_IJVectorInitialize(newB);
3620 ierr += HYPRE_IJVectorAssemble(newB);
3621 ierr += HYPRE_IJVectorGetObject(newB, (void **) &b2_csr);
3622 b2Data = hypre_VectorData(hypre_ParVectorLocalVector(b2_csr));
3623 hypre_assert( !ierr );
3624 for ( irow = 0; irow < localNRows; irow++ )
3625 b2Data[irow] = bData[irow] * scaleVec[irow];
3626
3627 ADiagISqrts_ = scaleVec;
3628 reducedAmat_ = newA;
3629 reducedBvec_ = newB;
3630 return 0;
3631 }
3632
3633 //****************************************************************************
3634 // estimate conditioning of a small block
3635 //----------------------------------------------------------------------------
3636
matrixCondEst(int globalRowID,int globalColID,int * blkInfo,int blkCnt)3637 double HYPRE_SlideReduction::matrixCondEst(int globalRowID, int globalColID,
3638 int *blkInfo,int blkCnt)
3639 {
3640 int mypid, nprocs, *procNRows, endRow, nConstraints;
3641 int localBlkCnt, *localBlkInfo, irow, jcol, matDim, searchIndex;
3642 int *rowIndices, rowSize, *colInd, rowIndex, index2, ierr;
3643 int *localSlaveEqns, *localSlaveAuxs;
3644 double *colVal, **matrix, **matrix2, retVal, value;
3645 HYPRE_ParCSRMatrix A_csr;
3646
3647 //------------------------------------------------------------------
3648 // get matrix information
3649 //------------------------------------------------------------------
3650
3651 MPI_Comm_rank(mpiComm_, &mypid);
3652 MPI_Comm_size(mpiComm_, &nprocs);
3653 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
3654 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
3655 endRow = procNRows[mypid+1] - 1;
3656 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
3657 free( procNRows );
3658
3659 //------------------------------------------------------------------
3660 // collect all row indices
3661 //------------------------------------------------------------------
3662
3663 localBlkCnt = blkCnt;
3664 localBlkInfo = new int[blkCnt];
3665 for (irow = 0; irow < blkCnt; irow++) localBlkInfo[irow] = blkInfo[irow];
3666
3667 hypre_qsort0(localBlkInfo, 0, localBlkCnt-1);
3668 matDim = 1;
3669 for ( irow = 0; irow < nConstraints; irow++ )
3670 {
3671 searchIndex = hypre_BinarySearch(localBlkInfo, constrBlkInfo_[irow],
3672 localBlkCnt);
3673 if ( searchIndex >= 0 ) matDim++;
3674 }
3675 rowIndices = new int[matDim];
3676 matDim = 0;
3677 rowIndices[matDim++] = globalRowID;
3678 for ( irow = 0; irow < nConstraints; irow++ )
3679 {
3680 searchIndex = hypre_BinarySearch(localBlkInfo, constrBlkInfo_[irow],
3681 localBlkCnt);
3682 if ( searchIndex >= 0 )
3683 rowIndices[matDim++] = endRow - nConstraints + irow + 1;
3684 }
3685 hypre_qsort0(rowIndices, 0, matDim-1);
3686 matrix = hypre_TAlloc(double*, matDim , HYPRE_MEMORY_HOST);
3687 localSlaveEqns = new int[nConstraints];
3688 localSlaveAuxs = new int[nConstraints];
3689 for ( irow = 0; irow < nConstraints; irow++ )
3690 localSlaveEqns[irow] = slaveEqnList_[irow];
3691 localSlaveEqns[globalRowID-(endRow+1-nConstraints)] = globalColID;
3692 for ( irow = 0; irow < nConstraints; irow++ )
3693 localSlaveAuxs[irow] = irow;
3694 HYPRE_LSI_qsort1a(localSlaveEqns, localSlaveAuxs, 0, nConstraints-1);
3695
3696 for ( irow = 0; irow < matDim; irow++ )
3697 {
3698 matrix[irow] = hypre_TAlloc(double, matDim , HYPRE_MEMORY_HOST);
3699 for ( jcol = 0; jcol < matDim; jcol++ ) matrix[irow][jcol] = 0.0;
3700 }
3701 for ( irow = 0; irow < matDim; irow++ )
3702 {
3703 rowIndex = rowIndices[irow];
3704 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
3705 for ( jcol = 0; jcol < rowSize; jcol++ )
3706 {
3707 searchIndex = hypre_BinarySearch(localSlaveEqns,colInd[jcol],
3708 nConstraints);
3709 if ( searchIndex >= 0 )
3710 {
3711 index2 = localSlaveAuxs[searchIndex] + endRow - nConstraints + 1;
3712 searchIndex = hypre_BinarySearch(rowIndices,index2,matDim);
3713 if ( searchIndex >= 0 ) matrix[irow][searchIndex] = colVal[jcol];
3714 }
3715 }
3716 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
3717 }
3718 #if 0
3719 if ( matDim <= 4 )
3720 for ( irow = 0; irow < matDim; irow++ )
3721 {
3722 for ( jcol = 0; jcol < matDim; jcol++ )
3723 printf(" %e ", matrix[irow][jcol]);
3724 printf("\n");
3725 }
3726 #endif
3727 ierr = HYPRE_LSI_MatrixInverse((double**) matrix,matDim,&matrix2);
3728 if ( ierr ) retVal = 1.0e-10;
3729 else
3730 {
3731 retVal = 0.0;
3732 for ( irow = 0; irow < matDim; irow++ )
3733 {
3734 for ( jcol = 0; jcol < matDim; jcol++ )
3735 {
3736 value = habs(matrix2[irow][jcol]);
3737 retVal = ( value > retVal ) ? value : retVal;
3738 }
3739 }
3740 retVal = 1.0 / retVal;
3741 for ( irow = 0; irow < matDim; irow++ ) free(matrix2[irow]);
3742 free( matrix2 );
3743 }
3744 for ( irow = 0; irow < matDim; irow++ ) free(matrix[irow]);
3745 free( matrix );
3746 delete [] localBlkInfo;
3747 delete [] rowIndices;
3748 delete [] localSlaveEqns;
3749 delete [] localSlaveAuxs;
3750 return retVal;
3751 }
3752
3753 //***************************************************************************
3754 //***************************************************************************
3755 // Obsolete, but keep it for now
3756 //***************************************************************************
3757 // search for a slave equation list (block size = 2)
3758 //---------------------------------------------------------------------------
3759
findSlaveEqns2(int ** couplings)3760 int HYPRE_SlideReduction::findSlaveEqns2(int **couplings)
3761 {
3762 int mypid, nprocs, *procNRows, startRow, endRow;
3763 int nConstraints, irow, jcol, rowSize, ncnt, *colInd;
3764 int nCandidates, *candidateList;
3765 int *constrListAux, colIndex, searchIndex, newEndRow;
3766 int *constrListAux2;
3767 int constrIndex, uBound, lBound, nSum, nPairs, index;
3768 double *colVal, searchValue;
3769 HYPRE_ParCSRMatrix A_csr;
3770
3771 //------------------------------------------------------------------
3772 // get matrix information
3773 //------------------------------------------------------------------
3774
3775 MPI_Comm_rank(mpiComm_, &mypid);
3776 MPI_Comm_size(mpiComm_, &nprocs);
3777 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
3778 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
3779 startRow = procNRows[mypid];
3780 endRow = procNRows[mypid+1] - 1;
3781 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
3782 newEndRow = endRow - nConstraints;
3783 nPairs = 0;
3784 for ( irow = 0; irow < nConstraints; irow++ )
3785 if ( slaveEqnList_[irow] == -1 ) nPairs++;
3786 (*couplings) = new int[nPairs*2+1];
3787 (*couplings)[0] = nPairs;
3788
3789 //------------------------------------------------------------------
3790 // compose candidate slave list (slaves in candidateList, corresponding
3791 // constraint equation in constrListAux and constrListAux2)
3792 //------------------------------------------------------------------
3793
3794 nCandidates = 0;
3795 candidateList = NULL;
3796 constrListAux = NULL;
3797 constrListAux2 = NULL;
3798 if ( nConstraints > 0 )
3799 {
3800 candidateList = new int[endRow-nConstraints-startRow+1];
3801 constrListAux = new int[endRow-nConstraints-startRow+1];
3802 constrListAux2 = new int[endRow-nConstraints-startRow+1];
3803
3804 //---------------------------------------------------------------
3805 // candidates are those with 2 links to the constraint list
3806 //---------------------------------------------------------------
3807
3808 uBound = procNRows[mypid+1];
3809 lBound = uBound - nConstraints;
3810
3811 for ( irow = startRow; irow <= endRow-nConstraints; irow++ )
3812 {
3813 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
3814 ncnt = 0;
3815 constrListAux[nCandidates] = -1;
3816 constrListAux2[nCandidates] = -1;
3817 for ( jcol = 0; jcol < rowSize; jcol++ )
3818 {
3819 colIndex = colInd[jcol];
3820 if ( colIndex >= lBound && colIndex < uBound )
3821 {
3822 ncnt++;
3823
3824 if ( ncnt == 1 && constrListAux[nCandidates] == -1 )
3825 constrListAux[nCandidates] = colIndex;
3826 else if ( ncnt == 2 && constrListAux2[nCandidates] == -1 )
3827 constrListAux2[nCandidates] = colIndex;
3828 }
3829 if ( ncnt > 2 ) break;
3830 }
3831 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
3832 if ( ncnt == 2 )
3833 {
3834 if ( constrListAux[nCandidates] > newEndRow &&
3835 constrListAux[nCandidates] <= endRow &&
3836 constrListAux2[nCandidates] > newEndRow &&
3837 constrListAux2[nCandidates] <= endRow )
3838 {
3839 candidateList[nCandidates++] = irow;
3840 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
3841 printf("%4d : findSlaveEqns2 - candidate %d = %d\n",
3842 mypid, nCandidates-1, irow);
3843 }
3844 }
3845 }
3846 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
3847 printf("%4d : findSlaveEqns2 - nCandidates, nConstr = %d %d\n",
3848 mypid, nCandidates, nConstraints);
3849 }
3850
3851 //---------------------------------------------------------------------
3852 // search the constraint equations for the selected slave equations
3853 // (search for candidates column index with maximum magnitude)
3854 // ==> slaveEqnList_
3855 //---------------------------------------------------------------------
3856
3857 nPairs = 0;
3858 searchIndex = 0;
3859
3860 for ( irow = endRow-nConstraints+1; irow <= endRow; irow++ )
3861 {
3862 if ( slaveEqnList_[irow-endRow+nConstraints-1] == -1 )
3863 {
3864 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
3865 searchIndex = -1;
3866 searchValue = -1.0E10;
3867 for ( jcol = 0; jcol < rowSize; jcol++ )
3868 {
3869 if (colVal[jcol] != 0.0 && colInd[jcol] >= startRow
3870 && colInd[jcol] <= (endRow-nConstraints))
3871 {
3872 colIndex = hypre_BinarySearch(candidateList, colInd[jcol],
3873 nCandidates);
3874
3875 /* -- if the column corresponds to a candidate, then -- */
3876 /* -- see if that candidate has a constraint connection -- */
3877 /* -- that has been satisfied (slaveEqnList_[irow]) -- */
3878
3879 if ( colIndex >= 0 )
3880 {
3881 constrIndex = constrListAux[colIndex];
3882 if ( constrIndex == irow )
3883 constrIndex = constrListAux2[colIndex];
3884 if (slaveEqnList_[constrIndex-endRow+nConstraints-1] != -1)
3885 {
3886 if ( habs(colVal[jcol]) > searchValue )
3887 {
3888 searchValue = habs(colVal[jcol]);
3889 searchIndex = colInd[jcol];
3890 }
3891 }
3892 }
3893 }
3894 }
3895 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
3896 if ( searchIndex >= 0 )
3897 {
3898 slaveEqnList_[irow-endRow+nConstraints-1] = searchIndex;
3899 index = hypre_BinarySearch(candidateList,searchIndex,nCandidates);
3900 (*couplings)[nPairs*2+1] = constrListAux[index];
3901 (*couplings)[nPairs*2+2] = constrListAux2[index];
3902 nPairs++;
3903 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
3904 printf("%4d : findSlaveEqns2 - constr %d <=> slave %d\n",
3905 mypid, irow, searchIndex);
3906 }
3907 else
3908 {
3909 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
3910 {
3911 printf("%4d : findSlaveEqns2 - constraint %4d fails", mypid,
3912 irow);
3913 printf(" to find a slave.\n");
3914 }
3915 break;
3916 }
3917 }
3918 }
3919 if ( nConstraints > 0 )
3920 {
3921 delete [] constrListAux;
3922 delete [] constrListAux2;
3923 delete [] candidateList;
3924 }
3925 free( procNRows );
3926
3927 //---------------------------------------------------------------------
3928 // if not all constraint-slave pairs can be found, return -1
3929 //---------------------------------------------------------------------
3930
3931 ncnt = 0;
3932 for ( irow = 0; irow < nConstraints; irow++ )
3933 if ( slaveEqnList_[irow] == -1 ) ncnt++;
3934 MPI_Allreduce(&ncnt, &nSum, 1, MPI_INT, MPI_SUM, mpiComm_);
3935 if ( nSum > 0 )
3936 {
3937 if ( mypid == 0 && ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
3938 {
3939 printf("%4d : findSlaveEqns2 fails - total number of unsatisfied",
3940 mypid);
3941 printf(" constraints = %d \n", nSum);
3942 }
3943 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
3944 {
3945 for ( irow = 0; irow < nConstraints; irow++ )
3946 if ( slaveEqnList_[irow] == -1 )
3947 {
3948 printf("%4d : findSlaveEqns2 - unsatisfied constraint",mypid);
3949 printf(" equation = %d\n", irow+endRow-nConstraints+1);
3950 }
3951 }
3952 return -1;
3953 }
3954 else return 0;
3955 }
3956
3957 //****************************************************************************
3958 // build reduced matrix
3959 //----------------------------------------------------------------------------
3960
buildReducedMatrix2()3961 int HYPRE_SlideReduction::buildReducedMatrix2()
3962 {
3963 int mypid, nprocs, *procNRows, startRow, endRow, localNRows;
3964 int globalNConstr, globalNRows, nConstraints, A21NRows, A21NCols;
3965 int A21GlobalNRows, A21GlobalNCols, A21StartRow, A21StartCol, ierr;
3966 int rowCount, maxRowSize, newEndRow, *A21MatSize, irow, is, rowIndex;
3967 int rowSize, *colInd, newRowSize, jcol, colIndex, searchIndex;
3968 int nnzA21, *newColInd, procIndex, ncnt, uBound, *colInd2, rowSize2;
3969 int newColIndex, *rowTags=NULL;
3970 double *colVal, *newColVal, mat2X2[4], denom, *colVal2;
3971 char fname[40];
3972 HYPRE_ParCSRMatrix A_csr, A21_csr, invA22_csr, RAP_csr, reducedA_csr;
3973
3974 //------------------------------------------------------------------
3975 // get matrix information
3976 //------------------------------------------------------------------
3977
3978 MPI_Comm_rank(mpiComm_, &mypid);
3979 MPI_Comm_size(mpiComm_, &nprocs);
3980 HYPRE_IJMatrixGetObject(Amat_, (void **) &A_csr);
3981 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &procNRows );
3982 startRow = procNRows[mypid];
3983 endRow = procNRows[mypid+1] - 1;
3984 localNRows = endRow - startRow + 1;
3985 globalNConstr = procNConstr_[nprocs];
3986 globalNRows = procNRows[nprocs];
3987 nConstraints = procNConstr_[mypid+1] - procNConstr_[mypid];
3988
3989 //******************************************************************
3990 // extract A21 from A
3991 //------------------------------------------------------------------
3992 // calculate the dimension of A21
3993 //------------------------------------------------------------------
3994
3995 A21NRows = 2 * nConstraints;
3996 A21NCols = localNRows - nConstraints;
3997 A21GlobalNRows = 2 * globalNConstr;
3998 A21GlobalNCols = globalNRows - globalNConstr;
3999 A21StartRow = 2 * procNConstr_[mypid];
4000 A21StartCol = procNRows[mypid] - procNConstr_[mypid];
4001
4002 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
4003 {
4004 printf("%4d : buildReducedMatrix - A21StartRow = %d\n", mypid,
4005 A21StartRow);
4006 printf("%4d : buildReducedMatrix - A21GlobalDim = %d %d\n", mypid,
4007 A21GlobalNRows, A21GlobalNCols);
4008 printf("%4d : buildReducedMatrix - A21LocalDim = %d %d\n",mypid,
4009 A21NRows, A21NCols);
4010 }
4011
4012 //------------------------------------------------------------------
4013 // create a matrix context for A21
4014 //------------------------------------------------------------------
4015
4016 ierr = HYPRE_IJMatrixCreate(mpiComm_,A21StartRow,A21StartRow+A21NRows-1,
4017 A21StartCol,A21StartCol+A21NCols-1,&A21mat_);
4018 ierr += HYPRE_IJMatrixSetObjectType(A21mat_, HYPRE_PARCSR);
4019 hypre_assert(!ierr);
4020
4021 //------------------------------------------------------------------
4022 // compute the number of nonzeros in the first nConstraint row of A21
4023 // (which consists of the rows in selectedList), the nnz will
4024 // be reduced by excluding the constraint and selected slave columns
4025 //------------------------------------------------------------------
4026
4027 rowCount = maxRowSize = 0;
4028 newEndRow = endRow - nConstraints;
4029 A21MatSize = new int[A21NRows];
4030 for ( irow = 0; irow < nConstraints; irow++ )
4031 {
4032 for ( is = 0; is < nConstraints; is++ )
4033 {
4034 if ( slaveEqnListAux_[is] == irow )
4035 {
4036 rowIndex = slaveEqnList_[is];
4037 break;
4038 }
4039 }
4040 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4041 newRowSize = 0;
4042 for ( jcol = 0; jcol < rowSize; jcol++ )
4043 {
4044 colIndex = colInd[jcol];
4045 if ( colVal[jcol] != 0.0 )
4046 {
4047 if ( colIndex <= newEndRow || colIndex > endRow )
4048 {
4049 searchIndex = hypre_BinarySearch(gSlaveEqnList_,colIndex,
4050 globalNConstr);
4051 if ( searchIndex < 0 ) newRowSize++;
4052 }
4053 }
4054 }
4055 A21MatSize[irow] = newRowSize;
4056 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
4057 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4058 }
4059
4060 //------------------------------------------------------------------
4061 // compute the number of nonzeros in the second nConstraint row of A21
4062 // (which consists of the rows in constraint equations)
4063 //------------------------------------------------------------------
4064
4065 rowCount = nConstraints;
4066 for ( irow = endRow-nConstraints+1; irow <= endRow; irow++ )
4067 {
4068 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
4069 newRowSize = 0;
4070 for ( jcol = 0; jcol < rowSize; jcol++ )
4071 {
4072 if ( colVal[jcol] != 0.0 )
4073 {
4074 colIndex = colInd[jcol];
4075 if ( colIndex <= newEndRow || colIndex > endRow )
4076 {
4077 searchIndex = hypre_BinarySearch(gSlaveEqnList_,colIndex,
4078 globalNConstr);
4079 if ( searchIndex < 0 ) newRowSize++;
4080 }
4081 }
4082 }
4083 A21MatSize[rowCount] = newRowSize;
4084 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
4085 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
4086 rowCount++;
4087 }
4088 nnzA21 = 0;
4089 for ( irow = 0; irow < 2*nConstraints; irow++ ) nnzA21 += A21MatSize[irow];
4090 MPI_Allreduce(&nnzA21,&ncnt,1,MPI_INT,MPI_SUM,mpiComm_);
4091 if ( mypid == 0 && ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
4092 printf(" 0 : buildReducedMatrix : NNZ of A21 = %d\n", ncnt);
4093
4094 //------------------------------------------------------------------
4095 // after fetching the row sizes, set up A21 with such sizes
4096 //------------------------------------------------------------------
4097
4098 ierr = HYPRE_IJMatrixSetRowSizes(A21mat_, A21MatSize);
4099 ierr += HYPRE_IJMatrixInitialize(A21mat_);
4100 hypre_assert(!ierr);
4101 delete [] A21MatSize;
4102
4103 //------------------------------------------------------------------
4104 // next load the first nConstraint row to A21 extracted from A
4105 // (at the same time, the D block is saved for future use)
4106 //------------------------------------------------------------------
4107
4108 rowCount = A21StartRow;
4109 newColInd = new int[maxRowSize+1];
4110 newColVal = new double[maxRowSize+1];
4111
4112 for ( irow = 0; irow < nConstraints; irow++ )
4113 {
4114 for ( is = 0; is < nConstraints; is++ )
4115 {
4116 if ( slaveEqnListAux_[is] == irow )
4117 {
4118 rowIndex = slaveEqnList_[is];
4119 break;
4120 }
4121 }
4122 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4123 newRowSize = 0;
4124 for ( jcol = 0; jcol < rowSize; jcol++ )
4125 {
4126 if ( colVal[jcol] != 0.0 )
4127 {
4128 colIndex = colInd[jcol];
4129 if ( colIndex <= newEndRow || colIndex > endRow )
4130 {
4131 searchIndex = HYPRE_LSI_Search(gSlaveEqnList_,colIndex,
4132 globalNConstr);
4133 if ( searchIndex < 0 )
4134 {
4135 for ( procIndex = 0; procIndex < nprocs; procIndex++ )
4136 if ( procNRows[procIndex] > colIndex ) break;
4137 procIndex--;
4138 newColIndex = colIndex - procNConstr_[procIndex];
4139 newColInd[newRowSize] = newColIndex;
4140 newColVal[newRowSize++] = colVal[jcol];
4141 if ( newColIndex < 0 || newColIndex >= A21GlobalNCols )
4142 {
4143 printf("%4d : buildReducedMatrix ERROR - A21",mypid);
4144 printf(" out of range (%d,%d (%d))\n", rowCount,
4145 colIndex, A21GlobalNCols);
4146 for ( is = 0; is < rowSize; is++ )
4147 printf("%4d : row %7d has col = %7d\n",mypid,rowIndex,
4148 colInd[is]);
4149 exit(1);
4150 }
4151 if ( newRowSize > maxRowSize+1 )
4152 {
4153 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 2 )
4154 {
4155 printf("%4d : buildReducedMatrix WARNING - ",mypid);
4156 printf("passing array boundary(1).\n");
4157 }
4158 }
4159 }
4160 }
4161 }
4162 }
4163 HYPRE_IJMatrixSetValues(A21mat_,1,&newRowSize,(const int *) &rowCount,
4164 (const int *) newColInd, (const double *) newColVal);
4165 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4166 rowCount++;
4167 }
4168
4169 //------------------------------------------------------------------
4170 // next load the second nConstraint rows to A21 extracted from A
4171 //------------------------------------------------------------------
4172
4173 for ( irow = endRow-nConstraints+1; irow <= endRow; irow++ )
4174 {
4175 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,&colVal);
4176 newRowSize = 0;
4177 for ( jcol = 0; jcol < rowSize; jcol++ )
4178 {
4179 colIndex = colInd[jcol];
4180 if (colVal[jcol] != 0.0 &&
4181 (colIndex <= newEndRow || colIndex > endRow))
4182 {
4183 searchIndex = hypre_BinarySearch(gSlaveEqnList_, colIndex,
4184 globalNConstr);
4185 if ( searchIndex < 0 )
4186 {
4187 for ( procIndex = 0; procIndex < nprocs; procIndex++ )
4188 if ( procNRows[procIndex] > colIndex ) break;
4189 procIndex--;
4190 newColIndex = colInd[jcol] - procNConstr_[procIndex];
4191 newColInd[newRowSize] = newColIndex;
4192 newColVal[newRowSize++] = colVal[jcol];
4193 }
4194 }
4195 }
4196 HYPRE_IJMatrixSetValues(A21mat_,1,&newRowSize,(const int *) &rowCount,
4197 (const int *) newColInd, (const double *) newColVal);
4198 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
4199 rowCount++;
4200 }
4201 delete [] newColInd;
4202 delete [] newColVal;
4203
4204 //------------------------------------------------------------------
4205 // finally assemble the matrix and sanitize
4206 //------------------------------------------------------------------
4207
4208 HYPRE_IJMatrixAssemble(A21mat_);
4209 HYPRE_IJMatrixGetObject(A21mat_, (void **) &A21_csr);
4210 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) A21_csr);
4211
4212 if ( outputLevel_ >= 5 )
4213 {
4214 sprintf(fname, "A21.%d", mypid);
4215 FILE *fp = fopen(fname, "w");
4216
4217 if ( mypid == 0 )
4218 {
4219 printf("====================================================\n");
4220 printf("%4d : Printing A21 matrix... \n", mypid);
4221 fflush(stdout);
4222 }
4223 for (irow = A21StartRow;irow < A21StartRow+2*nConstraints;irow++)
4224 {
4225 HYPRE_ParCSRMatrixGetRow(A21_csr,irow,&rowSize,&colInd,&colVal);
4226 for ( jcol = 0; jcol < rowSize; jcol++ )
4227 if ( colVal[jcol] != 0.0 )
4228 fprintf(fp,"%6d %6d %25.16e \n",irow+1,colInd[jcol]+1,
4229 colVal[jcol]);
4230 HYPRE_ParCSRMatrixRestoreRow(A21_csr, irow, &rowSize,
4231 &colInd, &colVal);
4232 }
4233 if ( mypid == 0 )
4234 printf("====================================================\n");
4235 fclose(fp);
4236 }
4237
4238 //******************************************************************
4239 // construct invA22
4240 // - given A22 = | B C |, compute | 0 C^{-T} |
4241 // | C^T 0 | | C^{-1} -C^{-1} B C^{-T} |
4242 //------------------------------------------------------------------
4243
4244 //------------------------------------------------------------------
4245 // first extract the (2,1) block of A22
4246 // ( constraints-to-local slaves )
4247 //------------------------------------------------------------------
4248
4249 int *CT_JA = NULL;
4250 double *CT_AA = NULL;
4251 if ( nConstraints > 0 )
4252 {
4253 CT_JA = new int[nConstraints*2];
4254 CT_AA = new double[nConstraints*2];
4255 }
4256 for ( irow = 0; irow < nConstraints; irow++ )
4257 {
4258 rowIndex = newEndRow + 1 + irow;
4259 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4260 CT_JA[irow*2] = CT_JA[irow*2+1] = -1;
4261 CT_AA[irow*2] = CT_AA[irow*2+1] = 0.0;
4262 for ( jcol = 0; jcol < rowSize; jcol++ )
4263 {
4264 colIndex = colInd[jcol];
4265 searchIndex = hypre_BinarySearch(slaveEqnList_,colIndex,nConstraints);
4266 if ( searchIndex >= 0 )
4267 {
4268 if ( CT_JA[irow*2] == -1 )
4269 {
4270 CT_JA[irow*2] = slaveEqnListAux_[searchIndex];
4271 CT_AA[irow*2] = colVal[jcol];
4272 }
4273 else
4274 {
4275 CT_JA[irow*2+1] = slaveEqnListAux_[searchIndex];
4276 CT_AA[irow*2+1] = colVal[jcol];
4277 }
4278 }
4279 }
4280 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4281 }
4282
4283 //------------------------------------------------------------------
4284 // invert the (2,1) block of A22
4285 //------------------------------------------------------------------
4286
4287 if ( nConstraints > 0 ) rowTags = new int[nConstraints];
4288 for ( irow = 0; irow < nConstraints; irow++ ) rowTags[irow] = -1;
4289
4290 for ( irow = 0; irow < nConstraints; irow++ )
4291 {
4292 if ( rowTags[irow] == -1 )
4293 {
4294 if ( CT_JA[irow*2+1] == -1 )
4295 CT_AA[irow*2] = 1.0 / CT_AA[irow*2];
4296 else
4297 {
4298 if ( CT_JA[2*irow] == irow )
4299 {
4300 mat2X2[0] = CT_AA[2*irow];
4301 mat2X2[2] = CT_AA[2*irow+1];
4302 rowIndex = CT_JA[2*irow+1];
4303 }
4304 else
4305 {
4306 mat2X2[0] = CT_AA[2*irow+1];
4307 mat2X2[2] = CT_AA[2*irow];
4308 rowIndex = CT_JA[2*irow];
4309 }
4310 if ( rowTags[rowIndex] != -1 )
4311 CT_AA[rowIndex*2] = 1.0 / CT_AA[rowIndex*2];
4312 if ( CT_JA[2*rowIndex] == rowIndex )
4313 {
4314 mat2X2[3] = CT_AA[2*rowIndex];
4315 mat2X2[1] = CT_AA[2*rowIndex+1];
4316 }
4317 else
4318 {
4319 mat2X2[3] = CT_AA[2*rowIndex+1];
4320 mat2X2[1] = CT_AA[2*rowIndex];
4321 }
4322 rowTags[rowIndex] = 0;
4323 denom = mat2X2[0] * mat2X2[3] - mat2X2[1] * mat2X2[2];
4324 denom = 1.0 / denom;
4325 CT_JA[irow*2] = irow;
4326 CT_AA[irow*2] = mat2X2[3] * denom;
4327 CT_JA[irow*2+1] = rowIndex;
4328 CT_AA[irow*2+1] = - mat2X2[2] * denom;
4329 CT_JA[rowIndex*2] = rowIndex;
4330 CT_AA[rowIndex*2] = mat2X2[0] * denom;
4331 CT_JA[rowIndex*2+1] = irow;
4332 CT_AA[rowIndex*2+1] = - mat2X2[1] * denom;
4333 }
4334 rowTags[irow] = 0;
4335 }
4336 }
4337 if ( nConstraints > 0 ) delete [] rowTags;
4338
4339 //------------------------------------------------------------------
4340 // form ParCSRMatrix of the (2,1) block of A22
4341 //------------------------------------------------------------------
4342
4343 int *hypreCTMatSize;
4344 hypre_ParCSRMatrix *hypreCT;
4345 HYPRE_IJMatrix IJCT;
4346
4347 ierr = HYPRE_IJMatrixCreate(mpiComm_, procNConstr_[mypid],
4348 procNConstr_[mypid]+nConstraints-1, procNConstr_[mypid],
4349 procNConstr_[mypid]+nConstraints-1, &IJCT);
4350 ierr += HYPRE_IJMatrixSetObjectType(IJCT, HYPRE_PARCSR);
4351 hypre_assert(!ierr);
4352 if ( nConstraints > 0 ) hypreCTMatSize = new int[nConstraints];
4353 maxRowSize = 0;
4354 for ( irow = 0; irow < nConstraints; irow++ )
4355 {
4356 hypreCTMatSize[irow] = 1;
4357 if ( CT_JA[irow*2+1] != -1 && CT_AA[irow*2+1] != 0.0 )
4358 hypreCTMatSize[irow]++;
4359 }
4360 ierr = HYPRE_IJMatrixSetRowSizes(IJCT, hypreCTMatSize);
4361 ierr = HYPRE_IJMatrixInitialize(IJCT);
4362 hypre_assert(!ierr);
4363 if ( nConstraints > 0 ) delete [] hypreCTMatSize;
4364 newColInd = new int[2];
4365 newColVal = new double[2];
4366 for ( irow = 0; irow < nConstraints; irow++ )
4367 {
4368 rowIndex = procNConstr_[mypid] + irow;
4369 newColInd[0] = CT_JA[irow*2] + procNConstr_[mypid];
4370 newColVal[0] = CT_AA[irow*2];
4371 newRowSize = 1;
4372 if ( CT_JA[irow*2+1] != -1 && CT_AA[irow*2+1] != 0.0 )
4373 {
4374 newColInd[1] = CT_JA[irow*2+1] + procNConstr_[mypid];
4375 newColVal[1] = CT_AA[irow*2+1];
4376 newRowSize++;
4377 }
4378 ierr = HYPRE_IJMatrixSetValues(IJCT, 1, &newRowSize,
4379 (const int *) &rowIndex, (const int *) newColInd,
4380 (const double *) newColVal);
4381 hypre_assert( !ierr );
4382 }
4383 delete [] newColInd;
4384 delete [] newColVal;
4385 if ( nConstraints > 0 )
4386 {
4387 delete [] CT_JA;
4388 delete [] CT_AA;
4389 }
4390 HYPRE_IJMatrixAssemble(IJCT);
4391 HYPRE_IJMatrixGetObject(IJCT, (void **) &hypreCT);
4392 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) hypreCT);
4393
4394 //------------------------------------------------------------------
4395 // next extract the (1,2) block of A22
4396 // ( local slaves-to-constraints )
4397 //------------------------------------------------------------------
4398
4399 int *C_JA = NULL;
4400 double *C_AA = NULL;
4401 if ( nConstraints > 0 )
4402 {
4403 C_JA = new int[nConstraints*2];
4404 C_AA = new double[nConstraints*2];
4405 }
4406 for ( irow = 0; irow < nConstraints; irow++ )
4407 {
4408 for ( is = 0; is < nConstraints; is++ )
4409 {
4410 if ( slaveEqnListAux_[is] == irow )
4411 {
4412 rowIndex = slaveEqnList_[is];
4413 break;
4414 }
4415 }
4416 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4417 C_JA[irow*2] = C_JA[irow*2+1] = -1;
4418 C_AA[irow*2] = C_AA[irow*2+1] = 0.0;
4419 for ( jcol = 0; jcol < rowSize; jcol++ )
4420 {
4421 colIndex = colInd[jcol];
4422 if ( colIndex > newEndRow && colIndex <= endRow )
4423 {
4424 if ( C_JA[irow*2] == -1 )
4425 {
4426 C_JA[irow*2] = colIndex - newEndRow - 1;
4427 C_AA[irow*2] = colVal[jcol];
4428 }
4429 else
4430 {
4431 C_JA[irow*2+1] = colIndex - newEndRow - 1;
4432 C_AA[irow*2+1] = colVal[jcol];
4433 }
4434 }
4435 }
4436 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4437 }
4438
4439 //------------------------------------------------------------------
4440 // invert the (1,2) block of A22
4441 //------------------------------------------------------------------
4442
4443 if ( nConstraints > 0 ) rowTags = new int[nConstraints];
4444 for ( irow = 0; irow < nConstraints; irow++ ) rowTags[irow] = -1;
4445
4446 for ( irow = 0; irow < nConstraints; irow++ )
4447 {
4448 if ( rowTags[irow] == -1 )
4449 {
4450 if ( C_JA[irow*2+1] == -1 )
4451 C_AA[irow*2] = 1.0 / C_AA[irow*2];
4452 else
4453 {
4454 if ( C_JA[2*irow] == irow )
4455 {
4456 mat2X2[0] = C_AA[2*irow];
4457 mat2X2[2] = C_AA[2*irow+1];
4458 rowIndex = C_JA[2*irow+1];
4459 }
4460 else
4461 {
4462 mat2X2[0] = C_AA[2*irow+1];
4463 mat2X2[2] = C_AA[2*irow];
4464 rowIndex = C_JA[2*irow];
4465 }
4466 if ( rowTags[rowIndex] != -1 )
4467 C_AA[rowIndex*2] = 1.0 / C_AA[rowIndex*2];
4468 if ( C_JA[2*rowIndex] == rowIndex )
4469 {
4470 mat2X2[3] = C_AA[2*rowIndex];
4471 mat2X2[1] = C_AA[2*rowIndex+1];
4472 }
4473 else
4474 {
4475 mat2X2[3] = C_AA[2*rowIndex+1];
4476 mat2X2[1] = C_AA[2*rowIndex];
4477 }
4478 rowTags[rowIndex] = 0;
4479 denom = mat2X2[0] * mat2X2[3] - mat2X2[1] * mat2X2[3];
4480 denom = 1.0 / denom;
4481 C_JA[irow*2] = irow;
4482 C_AA[irow*2] = mat2X2[3] * denom;
4483 C_JA[irow*2+1] = rowIndex;
4484 C_AA[irow*2+1] = - mat2X2[2] * denom;
4485 C_JA[rowIndex*2] = rowIndex;
4486 C_AA[rowIndex*2] = mat2X2[0] * denom;
4487 C_JA[rowIndex*2+1] = irow;
4488 C_AA[rowIndex*2+1] = - mat2X2[1] * denom;
4489 }
4490 rowTags[irow] = 0;
4491 }
4492 }
4493 if ( nConstraints > 0 ) delete [] rowTags;
4494
4495 //------------------------------------------------------------------
4496 // form ParCSRMatrix of the (1,2) block of A22
4497 //------------------------------------------------------------------
4498
4499 int *hypreCMatSize;
4500 hypre_ParCSRMatrix *hypreC;
4501 HYPRE_IJMatrix IJC;
4502
4503 ierr = HYPRE_IJMatrixCreate(mpiComm_, procNConstr_[mypid],
4504 procNConstr_[mypid]+nConstraints-1, procNConstr_[mypid],
4505 procNConstr_[mypid]+nConstraints-1, &IJC);
4506 ierr += HYPRE_IJMatrixSetObjectType(IJC, HYPRE_PARCSR);
4507 hypre_assert(!ierr);
4508 if ( nConstraints > 0 ) hypreCMatSize = new int[nConstraints];
4509 maxRowSize = 0;
4510 for ( irow = 0; irow < nConstraints; irow++ )
4511 {
4512 hypreCMatSize[irow] = 1;
4513 if ( C_JA[irow*2+1] != -1 && C_AA[irow*2+1] != 0.0 )
4514 hypreCMatSize[irow]++;
4515 }
4516 ierr = HYPRE_IJMatrixSetRowSizes(IJC, hypreCMatSize);
4517 ierr = HYPRE_IJMatrixInitialize(IJC);
4518 hypre_assert(!ierr);
4519 if ( nConstraints > 0 ) delete [] hypreCMatSize;
4520 newColInd = new int[2];
4521 newColVal = new double[2];
4522 for ( irow = 0; irow < nConstraints; irow++ )
4523 {
4524 rowIndex = procNConstr_[mypid] + irow;
4525 newColInd[0] = C_JA[irow*2] + procNConstr_[mypid];
4526 newColVal[0] = C_AA[irow*2];
4527 newRowSize = 1;
4528 if ( C_JA[irow*2+1] != -1 && C_AA[irow*2+1] != 0.0 )
4529 {
4530 newColInd[1] = C_JA[irow*2+1] + procNConstr_[mypid];
4531 newColVal[1] = C_AA[irow*2+1];
4532 newRowSize++;
4533 }
4534 ierr = HYPRE_IJMatrixSetValues(IJC, 1, &newRowSize,
4535 (const int *) &rowIndex, (const int *) newColInd,
4536 (const double *) newColVal);
4537 hypre_assert( !ierr );
4538 }
4539 delete [] newColInd;
4540 delete [] newColVal;
4541 if ( nConstraints > 0 )
4542 {
4543 delete [] C_JA;
4544 delete [] C_AA;
4545 }
4546 HYPRE_IJMatrixAssemble(IJC);
4547 HYPRE_IJMatrixGetObject(IJC, (void **) &hypreC);
4548 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) hypreC);
4549
4550 //------------------------------------------------------------------
4551 // form ParCSRMatrix of the (2,2) block of the invA22 matrix
4552 //------------------------------------------------------------------
4553
4554 int *hypreBMatSize;
4555 hypre_ParCSRMatrix *hypreB;
4556 HYPRE_IJMatrix IJB;
4557
4558 ierr = HYPRE_IJMatrixCreate(mpiComm_, procNConstr_[mypid],
4559 procNConstr_[mypid]+nConstraints-1, procNConstr_[mypid],
4560 procNConstr_[mypid]+nConstraints-1, &IJB);
4561 ierr = HYPRE_IJMatrixSetObjectType(IJB, HYPRE_PARCSR);
4562 hypre_assert(!ierr);
4563 if ( nConstraints > 0 ) hypreBMatSize = new int[nConstraints];
4564 maxRowSize = 0;
4565 for ( irow = 0; irow < nConstraints; irow++ )
4566 {
4567 for ( is = 0; is < nConstraints; is++ )
4568 {
4569 if ( slaveEqnListAux_[is] == irow )
4570 {
4571 rowIndex = slaveEqnList_[is];
4572 break;
4573 }
4574 }
4575 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4576 newRowSize = 0;
4577 for ( jcol = 0; jcol < rowSize; jcol++ )
4578 {
4579 colIndex = colInd[jcol];
4580 searchIndex = hypre_BinarySearch(gSlaveEqnList_, colIndex,
4581 globalNConstr);
4582 if ( searchIndex >= 0 ) newRowSize++;
4583 }
4584 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4585 hypreBMatSize[irow] = newRowSize;
4586 maxRowSize = (newRowSize > maxRowSize) ? newRowSize : maxRowSize;
4587 }
4588 ierr = HYPRE_IJMatrixSetRowSizes(IJB, hypreBMatSize);
4589 ierr = HYPRE_IJMatrixInitialize(IJB);
4590 hypre_assert(!ierr);
4591 if ( nConstraints > 0 ) delete [] hypreBMatSize;
4592
4593 if ( maxRowSize > 0 )
4594 {
4595 newColInd = new int[maxRowSize];
4596 newColVal = new double[maxRowSize];
4597 }
4598 for ( irow = 0; irow < nConstraints; irow++ )
4599 {
4600 for ( is = 0; is < nConstraints; is++ )
4601 {
4602 if ( slaveEqnListAux_[is] == irow )
4603 {
4604 rowIndex = slaveEqnList_[is];
4605 break;
4606 }
4607 }
4608 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4609 newRowSize = 0;
4610 for ( jcol = 0; jcol < rowSize; jcol++ )
4611 {
4612 colIndex = colInd[jcol];
4613 searchIndex = hypre_BinarySearch(gSlaveEqnList_, colIndex,
4614 globalNConstr);
4615 if ( searchIndex >= 0 )
4616 {
4617 newColInd[newRowSize] = gSlaveEqnListAux_[searchIndex];
4618 newColVal[newRowSize++] = - colVal[jcol];
4619 }
4620 }
4621 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
4622 rowIndex = procNConstr_[mypid] + irow;
4623 ierr = HYPRE_IJMatrixSetValues(IJB, 1, &newRowSize,
4624 (const int *) &rowIndex, (const int *) newColInd,
4625 (const double *) newColVal);
4626 hypre_assert( !ierr );
4627 }
4628 HYPRE_IJMatrixAssemble(IJB);
4629 HYPRE_IJMatrixGetObject(IJB, (void **) &hypreB);
4630 if ( maxRowSize > 0 )
4631 {
4632 delete [] newColInd;
4633 delete [] newColVal;
4634 }
4635
4636 //------------------------------------------------------------------
4637 // perform triple matrix product - C^{-1} B C^{-T}
4638 //------------------------------------------------------------------
4639
4640 HYPRE_ParCSRMatrix hypreCBC;
4641
4642 strcpy( fname, "hypreCT" );
4643 HYPRE_ParCSRMatrixPrint((HYPRE_ParCSRMatrix) hypreCT, fname);
4644 strcpy( fname, "hypreB" );
4645 HYPRE_ParCSRMatrixPrint((HYPRE_ParCSRMatrix) hypreB, fname);
4646 hypre_BoomerAMGBuildCoarseOperator(hypreCT, hypreB, hypreCT,
4647 (hypre_ParCSRMatrix **) &hypreCBC);
4648 #if 0
4649 for ( irow = 0; irow < nConstraints; irow++ )
4650 {
4651 HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCT,irow,&rowSize,
4652 &colInd,&colVal);
4653 for (jcol = 0; jcol < rowSize; jcol++ )
4654 printf("CT : %5d %5d %25.16e\n",irow+1,colInd[jcol]+1,colVal[jcol]);
4655 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCT,irow,&rowSize,
4656 &colInd,&colVal);
4657 }
4658 for ( irow = 0; irow < nConstraints; irow++ )
4659 {
4660 HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreB,irow,&rowSize,
4661 &colInd,&colVal);
4662 for (jcol = 0; jcol < rowSize; jcol++ )
4663 printf("B : %5d %5d %25.16e\n",irow+1,colInd[jcol]+1,colVal[jcol]);
4664 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreB,irow,&rowSize,
4665 &colInd,&colVal);
4666 }
4667 for ( irow = 0; irow < nConstraints; irow++ )
4668 {
4669 HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
4670 &colInd,&colVal);
4671 for (jcol = 0; jcol < rowSize; jcol++ )
4672 printf("CBC : %5d %5d %25.16e\n",irow+1,colInd[jcol]+1,colVal[jcol]);
4673 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
4674 &colInd,&colVal);
4675 }
4676 #endif
4677
4678 HYPRE_IJMatrixDestroy( IJB );
4679
4680 //------------------------------------------------------------------
4681 // calculate the dimension of invA22
4682 //------------------------------------------------------------------
4683
4684 int invA22NRows = A21NRows;
4685 int invA22NCols = invA22NRows;
4686 int invA22StartRow = A21StartRow;
4687 int invA22StartCol = invA22StartRow;
4688 int *invA22MatSize;
4689
4690 //------------------------------------------------------------------
4691 // create a matrix context for A22
4692 //------------------------------------------------------------------
4693
4694 ierr = HYPRE_IJMatrixCreate(mpiComm_, invA22StartRow,
4695 invA22StartRow+invA22NRows-1, invA22StartCol,
4696 invA22StartCol+invA22NCols-1, &invA22mat_);
4697 ierr += HYPRE_IJMatrixSetObjectType(invA22mat_, HYPRE_PARCSR);
4698 hypre_assert(!ierr);
4699
4700 //------------------------------------------------------------------
4701 // compute the no. of nonzeros in the first nConstraint row of invA22
4702 //------------------------------------------------------------------
4703
4704 maxRowSize = 0;
4705 invA22MatSize = new int[invA22NRows];
4706 for ( irow = 0; irow < nConstraints; irow++ )
4707 {
4708 rowIndex = procNConstr_[mypid] + irow;
4709 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCT,rowIndex,
4710 &rowSize,NULL,NULL);
4711 hypre_assert( !ierr );
4712 invA22MatSize[irow] = rowSize;
4713 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCT,rowIndex,
4714 &rowSize,NULL,NULL);
4715 maxRowSize = ( rowSize > maxRowSize ) ? rowSize : maxRowSize;
4716 }
4717
4718 //------------------------------------------------------------------
4719 // compute the number of nonzeros in the second nConstraints row of
4720 // invA22 (consisting of [D and A22 block])
4721 //------------------------------------------------------------------
4722
4723 #if 0
4724 for ( irow = 0; irow < nConstraints; irow++ )
4725 {
4726 HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
4727 &colInd,&colVal);
4728 for (jcol = 0; jcol < rowSize; jcol++ )
4729 printf("CBC1 : %5d %5d %25.16e\n",irow+1,colInd[jcol]+1,colVal[jcol]);
4730 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
4731 &colInd,&colVal);
4732 }
4733 #endif
4734 for ( irow = 0; irow < nConstraints; irow++ )
4735 {
4736 rowIndex = procNConstr_[mypid] + irow;
4737 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreC,rowIndex,
4738 &rowSize,&colInd,&colVal);
4739 hypre_assert( !ierr );
4740 newRowSize = rowSize;
4741 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreC,rowIndex,
4742 &rowSize,&colInd,&colVal);
4743 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCBC,rowIndex,
4744 &rowSize,&colInd,&colVal);
4745 hypre_assert( !ierr );
4746 newRowSize += rowSize;
4747 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCBC,rowIndex,
4748 &rowSize,&colInd,&colVal);
4749 invA22MatSize[nConstraints+irow] = newRowSize;
4750 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
4751 }
4752
4753 //------------------------------------------------------------------
4754 // after fetching the row sizes, set up invA22 with such sizes
4755 //------------------------------------------------------------------
4756
4757 #if 0
4758 for ( irow = 0; irow < nConstraints; irow++ )
4759 {
4760 HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
4761 &colInd,&colVal);
4762 for (jcol = 0; jcol < rowSize; jcol++ )
4763 printf("CBC2 : %5d %5d %25.16e\n",irow+1,colInd[jcol]+1,colVal[jcol]);
4764 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCBC,irow,&rowSize,
4765 &colInd,&colVal);
4766 }
4767 #endif
4768 ierr = HYPRE_IJMatrixSetRowSizes(invA22mat_, invA22MatSize);
4769 ierr += HYPRE_IJMatrixInitialize(invA22mat_);
4770 hypre_assert(!ierr);
4771 delete [] invA22MatSize;
4772
4773 //------------------------------------------------------------------
4774 // next load the first nConstraints row to invA22 extracted from A
4775 // (that is, the D block)
4776 //------------------------------------------------------------------
4777
4778 newColInd = new int[maxRowSize+1];
4779 newColVal = new double[maxRowSize+1];
4780 for ( irow = 0; irow < nConstraints; irow++ )
4781 {
4782 rowIndex = procNConstr_[mypid] + irow;
4783 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCT,rowIndex,
4784 &rowSize,&colInd,&colVal);
4785 hypre_assert( !ierr );
4786 newRowSize = 0;
4787 for ( jcol = 0; jcol < rowSize; jcol++ )
4788 {
4789 newColInd[newRowSize] = colInd[jcol] + procNConstr_[mypid] +
4790 nConstraints;
4791 newColVal[newRowSize++] = colVal[jcol];
4792 }
4793 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCT,rowIndex,
4794 &rowSize,&colInd,&colVal);
4795 rowCount = invA22StartRow + irow;
4796 ierr = HYPRE_IJMatrixSetValues(invA22mat_, 1, &rowSize,
4797 (const int *) &rowCount, (const int *) newColInd,
4798 (const double *) newColVal);
4799 hypre_assert(!ierr);
4800
4801 }
4802
4803 //------------------------------------------------------------------
4804 // next load the second nConstraints rows to A22 extracted from A
4805 //------------------------------------------------------------------
4806
4807 for ( irow = 0; irow < nConstraints; irow++ )
4808 {
4809 rowIndex = procNConstr_[mypid] + irow;
4810 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreC,rowIndex,
4811 &rowSize,&colInd,&colVal);
4812 hypre_assert( !ierr );
4813 newRowSize = 0;
4814 for ( jcol = 0; jcol < rowSize; jcol++ )
4815 {
4816 newColInd[newRowSize] = colInd[jcol] + procNConstr_[mypid];
4817 newColVal[newRowSize++] = colVal[jcol];
4818 }
4819 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreC,rowIndex,
4820 &rowSize,&colInd,&colVal);
4821 ierr = HYPRE_ParCSRMatrixGetRow((HYPRE_ParCSRMatrix) hypreCBC,rowIndex,
4822 &rowSize2,&colInd2,&colVal2);
4823 hypre_assert( !ierr );
4824 for ( jcol = 0; jcol < rowSize2; jcol++ )
4825 {
4826 newColInd[newRowSize] = colInd2[jcol] + procNConstr_[mypid] +
4827 nConstraints;
4828 newColVal[newRowSize++] = colVal2[jcol];
4829 }
4830 HYPRE_ParCSRMatrixRestoreRow((HYPRE_ParCSRMatrix) hypreCBC,rowIndex,
4831 &rowSize2,&colInd2,&colVal2);
4832 rowCount = invA22StartRow + nConstraints + irow;
4833 ierr = HYPRE_IJMatrixSetValues(invA22mat_, 1, &newRowSize,
4834 (const int *) &rowCount, (const int *) newColInd,
4835 (const double *) newColVal);
4836 hypre_assert(!ierr);
4837 }
4838 delete [] newColInd;
4839 delete [] newColVal;
4840 HYPRE_IJMatrixDestroy( IJC );
4841 HYPRE_IJMatrixDestroy( IJCT );
4842 HYPRE_ParCSRMatrixDestroy( (HYPRE_ParCSRMatrix) hypreCBC );
4843
4844 //------------------------------------------------------------------
4845 // finally assemble the matrix and sanitize
4846 //------------------------------------------------------------------
4847
4848 HYPRE_IJMatrixAssemble(invA22mat_);
4849 HYPRE_IJMatrixGetObject(invA22mat_, (void **) &invA22_csr);
4850 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) invA22_csr);
4851
4852 if ( outputLevel_ >= 5 )
4853 {
4854 sprintf(fname, "invA22.%d", mypid);
4855 FILE *fp = fopen(fname, "w");
4856
4857 if ( mypid == ncnt )
4858 {
4859 printf("====================================================\n");
4860 printf("%4d : Printing invA22 matrix... \n", mypid);
4861 fflush(stdout);
4862 }
4863 for (irow=invA22StartRow; irow < invA22StartRow+invA22NRows;irow++)
4864 {
4865 HYPRE_ParCSRMatrixGetRow(invA22_csr,irow,&rowSize,&colInd,&colVal);
4866 for ( jcol = 0; jcol < rowSize; jcol++ )
4867 if ( colVal[jcol] != 0.0 )
4868 fprintf(fp,"%6d %6d %25.16e \n",irow+1,colInd[jcol]+1,
4869 colVal[jcol]);
4870 HYPRE_ParCSRMatrixRestoreRow(invA22_csr,irow,&rowSize,&colInd,
4871 &colVal);
4872 }
4873 if ( mypid == ncnt )
4874 printf("====================================================\n");
4875 fclose(fp);
4876 }
4877
4878 //******************************************************************
4879 // perform the triple matrix product A12 * invA22 * A21
4880 //------------------------------------------------------------------
4881
4882 HYPRE_IJMatrixGetObject(A21mat_, (void **) &A21_csr);
4883 HYPRE_IJMatrixGetObject(invA22mat_, (void **) &invA22_csr);
4884 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
4885 printf("%4d : buildReducedMatrix - Triple matrix product starts\n",
4886 mypid);
4887
4888 hypre_BoomerAMGBuildCoarseOperator((hypre_ParCSRMatrix *) A21_csr,
4889 (hypre_ParCSRMatrix *) invA22_csr,
4890 (hypre_ParCSRMatrix *) A21_csr,
4891 (hypre_ParCSRMatrix **) &RAP_csr);
4892
4893 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
4894 printf("%4d : buildReducedMatrix - Triple matrix product ends\n",
4895 mypid);
4896
4897 if ( outputLevel_ >= 4 )
4898 {
4899 sprintf(fname, "rap.%d", mypid);
4900 FILE *fp = fopen(fname, "w");
4901
4902 if ( mypid == 0 )
4903 {
4904 printf("====================================================\n");
4905 printf("%4d : Printing RAP matrix... \n", mypid);
4906 fflush(stdout);
4907 }
4908 for ( irow = A21StartRow; irow < A21StartRow+A21NCols; irow++ )
4909 {
4910 HYPRE_ParCSRMatrixGetRow(RAP_csr,irow,&rowSize,&colInd,&colVal);
4911 for ( jcol = 0; jcol < rowSize; jcol++ )
4912 if ( colVal[jcol] != 0.0 )
4913 printf("%6d %6d %25.16e \n",irow+1,colInd[jcol]+1,
4914 colVal[jcol]);
4915 HYPRE_ParCSRMatrixRestoreRow(RAP_csr,irow,&rowSize,&colInd,
4916 &colVal);
4917 }
4918 fclose(fp);
4919 if ( mypid == 0 )
4920 printf("====================================================\n");
4921 }
4922
4923 //******************************************************************
4924 // finally form reduceA = A11 - A12 * invA22 * A21
4925 //------------------------------------------------------------------
4926
4927 //------------------------------------------------------------------
4928 // compute row sizes
4929 //------------------------------------------------------------------
4930
4931 int reducedANRows = localNRows - nConstraints;
4932 int reducedANCols = reducedANRows;
4933 int reducedAStartRow = procNRows[mypid] - procNConstr_[mypid];
4934 int reducedAStartCol = reducedAStartRow;
4935 int reducedAGlobalNRows = globalNRows - globalNConstr;
4936 int reducedAGlobalNCols = reducedAGlobalNRows;
4937 int *reducedAMatSize = new int[reducedANRows];
4938
4939 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 1 )
4940 {
4941 printf("%4d : buildReducedMatrix - reduceAGlobalDim = %d %d\n", mypid,
4942 reducedAGlobalNRows, reducedAGlobalNCols);
4943 printf("%4d : buildReducedMatrix - reducedALocalDim = %d %d\n", mypid,
4944 reducedANRows, reducedANCols);
4945 }
4946
4947 //------------------------------------------------------------------
4948 // create a matrix context for reducedA
4949 //------------------------------------------------------------------
4950
4951 ierr = HYPRE_IJMatrixCreate(mpiComm_,reducedAStartRow,
4952 reducedAStartRow+reducedANRows-1, reducedAStartCol,
4953 reducedAStartCol+reducedANCols-1,&reducedAmat_);
4954 ierr += HYPRE_IJMatrixSetObjectType(reducedAmat_, HYPRE_PARCSR);
4955 hypre_assert(!ierr);
4956
4957 //------------------------------------------------------------------
4958 // compute row sizes for reducedA
4959 //------------------------------------------------------------------
4960
4961 rowCount = maxRowSize = 0;
4962 for ( irow = startRow; irow <= newEndRow; irow++ )
4963 {
4964 searchIndex = hypre_BinarySearch(slaveEqnList_, irow, nConstraints);
4965 if ( searchIndex >= 0 ) reducedAMatSize[rowCount++] = 1;
4966 else
4967 {
4968 HYPRE_ParCSRMatrixGetRow(A_csr,irow,&rowSize,&colInd,NULL);
4969 rowIndex = reducedAStartRow + rowCount;
4970 ierr = HYPRE_ParCSRMatrixGetRow(RAP_csr,rowIndex,&rowSize2,
4971 &colInd2, NULL);
4972 hypre_assert( !ierr );
4973 newRowSize = rowSize + rowSize2;
4974 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
4975 newColInd = new int[newRowSize];
4976 for (jcol = 0; jcol < rowSize; jcol++) newColInd[jcol] = colInd[jcol];
4977 for (jcol = 0; jcol < rowSize2; jcol++)
4978 newColInd[rowSize+jcol] = colInd2[jcol];
4979 hypre_qsort0(newColInd, 0, newRowSize-1);
4980 ncnt = 0;
4981 for ( jcol = 1; jcol < newRowSize; jcol++ )
4982 if (newColInd[jcol] != newColInd[ncnt])
4983 newColInd[++ncnt] = newColInd[jcol];
4984 if ( newRowSize > 0 ) ncnt++;
4985 reducedAMatSize[rowIndex++] = ncnt;
4986 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,NULL);
4987 ierr = HYPRE_ParCSRMatrixRestoreRow(RAP_csr,rowIndex,&rowSize2,
4988 &colInd2,NULL);
4989 delete [] newColInd;
4990 hypre_assert( !ierr );
4991 rowCount++;
4992 }
4993 }
4994 ierr = HYPRE_IJMatrixSetRowSizes(reducedAmat_, reducedAMatSize);
4995 ierr += HYPRE_IJMatrixInitialize(reducedAmat_);
4996 hypre_assert(!ierr);
4997 delete [] reducedAMatSize;
4998
4999 //------------------------------------------------------------------
5000 // load the reducedA matrix
5001 //------------------------------------------------------------------
5002
5003 rowCount = 0;
5004 newColInd = new int[maxRowSize+1];
5005 newColVal = new double[maxRowSize+1];
5006 for ( irow = startRow; irow <= newEndRow; irow++ )
5007 {
5008 searchIndex = hypre_BinarySearch(slaveEqnList_, irow, nConstraints);
5009 rowIndex = reducedAStartRow + rowCount;
5010 if ( searchIndex >= 0 )
5011 {
5012 newRowSize = 1;
5013 newColInd[0] = reducedAStartRow + rowCount;
5014 newColVal[0] = 1.0;
5015 }
5016 else
5017 {
5018 HYPRE_ParCSRMatrixGetRow(A_csr, irow, &rowSize, &colInd, &colVal);
5019 HYPRE_ParCSRMatrixGetRow(RAP_csr,rowIndex,&rowSize2,&colInd2,
5020 &colVal2);
5021 newRowSize = rowSize + rowSize2;
5022 ncnt = 0;
5023 for ( jcol = 0; jcol < rowSize; jcol++ )
5024 {
5025 colIndex = colInd[jcol];
5026 for ( procIndex = 0; procIndex < nprocs; procIndex++ )
5027 if ( procNRows[procIndex] > colIndex ) break;
5028 uBound = procNRows[procIndex] -
5029 (procNConstr_[procIndex]-procNConstr_[procIndex-1]);
5030 procIndex--;
5031 if ( colIndex < uBound )
5032 {
5033 searchIndex = hypre_BinarySearch(gSlaveEqnList_, colIndex,
5034 globalNConstr);
5035 if ( searchIndex < 0 )
5036 {
5037 newColInd[ncnt] = colIndex - procNConstr_[procIndex];
5038 newColVal[ncnt++] = colVal[jcol];
5039 }
5040 }
5041 }
5042 for ( jcol = 0; jcol < rowSize2; jcol++ )
5043 {
5044 newColInd[ncnt+jcol] = colInd2[jcol];
5045 newColVal[ncnt+jcol] = - colVal2[jcol];
5046 }
5047 newRowSize = ncnt + rowSize2;
5048 hypre_qsort1(newColInd, newColVal, 0, newRowSize-1);
5049 ncnt = 0;
5050 for ( jcol = 0; jcol < newRowSize; jcol++ )
5051 {
5052 if ( jcol != ncnt && newColInd[jcol] == newColInd[ncnt] )
5053 newColVal[ncnt] += newColVal[jcol];
5054 else if ( newColInd[jcol] != newColInd[ncnt] )
5055 {
5056 ncnt++;
5057 newColVal[ncnt] = newColVal[jcol];
5058 newColInd[ncnt] = newColInd[jcol];
5059 }
5060 }
5061 newRowSize = ncnt + 1;
5062 HYPRE_ParCSRMatrixRestoreRow(A_csr,irow,&rowSize,&colInd,&colVal);
5063 HYPRE_ParCSRMatrixRestoreRow(RAP_csr,rowIndex,&rowSize2,&colInd2,
5064 &colVal2);
5065 }
5066 ierr = HYPRE_IJMatrixSetValues(reducedAmat_, 1, &newRowSize,
5067 (const int *) &rowIndex, (const int *) newColInd,
5068 (const double *) newColVal);
5069 hypre_assert(!ierr);
5070 rowCount++;
5071 }
5072 delete [] newColInd;
5073 delete [] newColVal;
5074
5075 //------------------------------------------------------------------
5076 // assemble the reduced matrix
5077 //------------------------------------------------------------------
5078
5079 HYPRE_IJMatrixAssemble(reducedAmat_);
5080 HYPRE_IJMatrixGetObject(reducedAmat_, (void **) &reducedA_csr);
5081
5082 if ( ( outputLevel_ & HYPRE_BITMASK2 ) >= 5 )
5083 {
5084 MPI_Barrier(mpiComm_);
5085 ncnt = 0;
5086 while ( ncnt < nprocs )
5087 {
5088 if ( mypid == ncnt )
5089 {
5090 printf("====================================================\n");
5091 printf("%4d : Printing reducedA matrix... \n", mypid);
5092 fflush(stdout);
5093 for ( irow = reducedAStartRow;
5094 irow < reducedAStartRow+localNRows-nConstraints; irow++ )
5095 {
5096 //printf("%d : reducedA ROW %d\n", mypid, irow);
5097 ierr = HYPRE_ParCSRMatrixGetRow(reducedA_csr,irow,&rowSize,
5098 &colInd, &colVal);
5099 //hypre_qsort1(colInd, colVal, 0, rowSize-1);
5100 for ( jcol = 0; jcol < rowSize; jcol++ )
5101 if ( colVal[jcol] != 0.0 )
5102 printf("%6d %6d %25.8e \n",irow+1,colInd[jcol]+1,
5103 colVal[jcol]);
5104 HYPRE_ParCSRMatrixRestoreRow(reducedA_csr,irow,&rowSize,
5105 &colInd, &colVal);
5106 }
5107 printf("====================================================\n");
5108 }
5109 MPI_Barrier(mpiComm_);
5110 ncnt++;
5111 }
5112 }
5113
5114 //------------------------------------------------------------------
5115 // store away matrix and clean up
5116 //------------------------------------------------------------------
5117
5118 free( procNRows );
5119 HYPRE_ParCSRMatrixDestroy(RAP_csr);
5120 return 0;
5121 }
5122
5123