1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 //***************************************************************************
9 // system includes
10 //---------------------------------------------------------------------------
11
12 #include <stdlib.h>
13 #include <string.h>
14 #include <stdio.h>
15
16 //***************************************************************************
17 // HYPRE includes
18 //---------------------------------------------------------------------------
19
20 #include "HYPRE.h"
21 #include "utilities/_hypre_utilities.h"
22 #include "IJ_mv/HYPRE_IJ_mv.h"
23 #include "parcsr_mv/HYPRE_parcsr_mv.h"
24 #include "parcsr_ls/HYPRE_parcsr_ls.h"
25 #include "parcsr_mv/_hypre_parcsr_mv.h"
26 #include "HYPRE_LinSysCore.h"
27 #include "HYPRE_LSI_mli.h"
28
29 //***************************************************************************
30 // local defines and external functions
31 //---------------------------------------------------------------------------
32
33 #define habs(x) (((x) > 0.0) ? x : -(x))
34
35 extern "C"
36 {
37 int hypre_BoomerAMGBuildCoarseOperator(hypre_ParCSRMatrix*,
38 hypre_ParCSRMatrix*,
39 hypre_ParCSRMatrix*,
40 hypre_ParCSRMatrix**);
41 void hypre_qsort0(int *, int, int);
42 void hypre_qsort1(int *, double *, int, int);
43 int HYPRE_LSI_Search(int*, int, int);
44 }
45
46 //*****************************************************************************
47 // Given the matrix (A) within the object, compute the reduced system and put
48 // it in place. Additional information given are :
49 //
50 // Additional assumptions are :
51 //
52 // - a given slave equation and the corresponding constraint equation
53 // reside in the same processor
54 // - constraint equations are given at the end of the local matrix
55 // (hence given by EndRow_-nConstr to EndRow_)
56 // - each processor gets a contiguous block of equations, and processor
57 // i+1 has equation numbers higher than those of processor i
58 //-----------------------------------------------------------------------------
59
buildSlideReducedSystem()60 void HYPRE_LinSysCore::buildSlideReducedSystem()
61 {
62 int i, j, StartRow, EndRow, rowSize, *colInd, globalNConstr;
63 int nRows, *ProcNRows, *tempList, globalNRows, ncnt;
64 int globalNSelected, *globalSelectedList, *globalSelectedListAux;
65 int *ProcNConstr, isAConstr, ncnt2;
66 double *colVal;
67 HYPRE_ParCSRMatrix A_csr, RAP_csr;
68
69 //******************************************************************
70 // initial set up
71 //------------------------------------------------------------------
72
73 if ( mypid_ == 0 && (HYOutputLevel_ & HYFEI_SLIDEREDUCE1) )
74 printf("%4d : SlideReduction begins....\n",mypid_);
75 StartRow = localStartRow_ - 1;
76 EndRow = localEndRow_ - 1;
77 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
78 printf("%4d : SlideReduction - StartRow/EndRow = %d %d\n",mypid_,
79 StartRow,EndRow);
80
81 //******************************************************************
82 // construct local and global information about where the constraints
83 // are (this is given by user or searched within this code)
84 //------------------------------------------------------------------
85
86 //------------------------------------------------------------------
87 // get the CSR matrix for A
88 //------------------------------------------------------------------
89
90 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
91
92 //------------------------------------------------------------------
93 // search the entire local matrix to find where the constraint
94 // equations are, if not already given
95 //------------------------------------------------------------------
96
97 MPI_Allreduce(&nConstraints_,&globalNConstr,1,MPI_INT,MPI_SUM,comm_);
98 if ( globalNConstr == 0 )
99 {
100 for ( i = EndRow; i >= StartRow; i-- )
101 {
102 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
103 isAConstr = 1;
104 for (j = 0; j < rowSize; j++)
105 if ( colInd[j] == i && colVal[j] != 0.0 ) {isAConstr = 0; break;}
106 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
107 if ( isAConstr ) nConstraints_++;
108 else break;
109 }
110 }
111 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
112 printf("%4d : SlideReduction - no. constr = %d\n",mypid_,nConstraints_);
113
114 MPI_Allreduce(&nConstraints_, &globalNConstr, 1, MPI_INT, MPI_SUM, comm_);
115 if ( globalNConstr == 0 ) return;
116
117 //------------------------------------------------------------------
118 // get information about nRows from all processors
119 //------------------------------------------------------------------
120
121 nRows = localEndRow_ - localStartRow_ + 1;
122 ProcNRows = new int[numProcs_];
123 tempList = new int[numProcs_];
124 for ( i = 0; i < numProcs_; i++ ) tempList[i] = 0;
125 tempList[mypid_] = nRows;
126 MPI_Allreduce(tempList, ProcNRows, numProcs_, MPI_INT, MPI_SUM, comm_);
127 delete [] tempList;
128 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
129 printf("%4d : SlideReduction - localNRows = %d\n", mypid_, nRows);
130
131 //------------------------------------------------------------------
132 // compute the base NRows on each processor
133 // (This is needed later on for column index conversion)
134 //------------------------------------------------------------------
135
136 globalNRows = 0;
137 ncnt = 0;
138 for ( i = 0; i < numProcs_; i++ )
139 {
140 globalNRows += ProcNRows[i];
141 ncnt2 = ProcNRows[i];
142 ProcNRows[i] = ncnt;
143 ncnt += ncnt2;
144 }
145
146 //------------------------------------------------------------------
147 // compose a global array marking where the constraint equations are
148 //------------------------------------------------------------------
149
150 globalNConstr = 0;
151 tempList = new int[numProcs_];
152 ProcNConstr = new int[numProcs_];
153 for ( i = 0; i < numProcs_; i++ ) tempList[i] = 0;
154 tempList[mypid_] = nConstraints_;
155 MPI_Allreduce(tempList,ProcNConstr,numProcs_,MPI_INT,MPI_SUM,comm_);
156 delete [] tempList;
157
158 //------------------------------------------------------------------
159 // compute the base nConstraints on each processor
160 // (This is needed later on for column index conversion)
161 //------------------------------------------------------------------
162
163 ncnt = 0;
164 for ( i = 0; i < numProcs_; i++ )
165 {
166 globalNConstr += ProcNConstr[i];
167 ncnt2 = ProcNConstr[i];
168 ProcNConstr[i] = ncnt;
169 ncnt += ncnt2;
170 }
171
172 //******************************************************************
173 // compose the local and global selected node lists
174 //------------------------------------------------------------------
175
176 //------------------------------------------------------------------
177 // allocate array for storing indices of selected nodes
178 //------------------------------------------------------------------
179
180 globalNSelected = globalNConstr;
181 if (globalNSelected > 0)
182 {
183 globalSelectedList = new int[globalNSelected];
184 globalSelectedListAux = new int[globalNSelected];
185 }
186 else globalSelectedList = globalSelectedListAux = NULL;
187 if ( selectedList_ != NULL ) delete [] selectedList_;
188 if ( selectedListAux_ != NULL ) delete [] selectedListAux_;
189 if ( nConstraints_ > 0 )
190 {
191 selectedList_ = new int[nConstraints_];
192 selectedListAux_ = new int[nConstraints_];
193 }
194 else selectedList_ = selectedListAux_ = NULL;
195
196 //------------------------------------------------------------------
197 // call the three parts
198 //------------------------------------------------------------------
199
200 buildSlideReducedSystemPartA(ProcNRows,ProcNConstr,globalNRows,
201 globalNConstr,globalSelectedList,
202 globalSelectedListAux);
203 buildSlideReducedSystemPartB(ProcNRows,ProcNConstr,globalNRows,
204 globalNConstr,globalSelectedList,
205 globalSelectedListAux, &RAP_csr);
206 buildSlideReducedSystemPartC(ProcNRows,ProcNConstr,globalNRows,
207 globalNConstr,globalSelectedList,
208 globalSelectedListAux, RAP_csr);
209
210 //------------------------------------------------------------------
211 // initialize global variables and clean up
212 //------------------------------------------------------------------
213
214 currA_ = reducedA_;
215 currB_ = reducedB_;
216 currR_ = reducedR_;
217 currX_ = reducedX_;
218 delete [] globalSelectedList;
219 delete [] globalSelectedListAux;
220 delete [] ProcNRows;
221 delete [] ProcNConstr;
222 HYPRE_ParCSRMatrixDestroy(RAP_csr);
223 //if ( HYA_ != NULL ) {HYPRE_IJMatrixDestroy(HYA_); HYA_ = NULL;}
224 if ( colIndices_ != NULL )
225 {
226 for ( i = 0; i < localEndRow_-localStartRow_+1; i++ )
227 if ( colIndices_[i] != NULL ) delete [] colIndices_[i];
228 delete [] colIndices_;
229 colIndices_ = NULL;
230 }
231 if ( colValues_ != NULL )
232 {
233 for ( j = 0; j < localEndRow_-localStartRow_+1; j++ )
234 if ( colValues_[j] != NULL ) delete [] colValues_[j];
235 delete [] colValues_;
236 colValues_ = NULL;
237 if ( rowLengths_ != NULL )
238 {
239 delete [] rowLengths_;
240 rowLengths_ = NULL;
241 }
242 }
243 }
244
245 //*****************************************************************************
246 // Part A of buildSlideReducedSystem : generate a selected equation list
247 //-----------------------------------------------------------------------------
248
buildSlideReducedSystemPartA(int * ProcNRows,int * ProcNConstr,int globalNRows,int globalNConstr,int * globalSelectedList,int * globalSelectedListAux)249 void HYPRE_LinSysCore::buildSlideReducedSystemPartA(int *ProcNRows,
250 int *ProcNConstr, int globalNRows,
251 int globalNConstr, int *globalSelectedList,
252 int *globalSelectedListAux)
253 {
254 int i, ncnt2, StartRow, EndRow, ncnt;;
255 int nSlaves, *constrListAux, colIndex, searchIndex, procIndex, ubound;
256 int j, k, ierr, rowSize, *colInd, *colInd2, rowIndex, nSelected;
257 int *recvCntArray, *displArray, *selectedList, *selectedListAux;
258 int rowSize2;
259 double *colVal, searchValue, *dble_array, *colVal2;
260 HYPRE_ParCSRMatrix A_csr;
261
262 //------------------------------------------------------------------
263 // get matrix information
264 //------------------------------------------------------------------
265
266 StartRow = localStartRow_ - 1;
267 EndRow = localEndRow_ - 1;
268 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
269 nSelected = nConstraints_;
270 selectedList = selectedList_;
271 selectedListAux = selectedListAux_;
272
273 //------------------------------------------------------------------
274 // compose candidate slave list
275 //------------------------------------------------------------------
276
277 nSlaves = 0;
278 constrListAux = NULL;
279 if ( nConstraints_ > 0 && constrList_ == NULL )
280 {
281 constrList_ = new int[EndRow-nConstraints_-StartRow+1];
282 constrListAux = new int[EndRow-nConstraints_-StartRow+1];
283
284 //------------------------------------------------------------------
285 // candidates are those with 1 link to the constraint list
286 //------------------------------------------------------------------
287
288 for ( i = StartRow; i <= EndRow-nConstraints_; i++ )
289 {
290 ierr = HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
291 hypre_assert(!ierr);
292 ncnt = 0;
293 for (j = 0; j < rowSize; j++)
294 {
295 colIndex = colInd[j];
296 for (procIndex=0; procIndex < numProcs_; procIndex++ )
297 if ( colIndex < ProcNRows[procIndex] ) break;
298 if ( procIndex == numProcs_ )
299 ubound = globalNRows -
300 (globalNConstr-ProcNConstr[procIndex-1]);
301 else
302 ubound = ProcNRows[procIndex] - (ProcNConstr[procIndex] -
303 ProcNConstr[procIndex-1]);
304
305 //Note : include structural zeros by not checking for nonzero
306 //if ( colIndex >= ubound && colVal[j] != 0.0 )
307 if ( colIndex >= ubound )
308 {
309 ncnt++;
310 searchIndex = colIndex;
311 }
312 if ( ncnt > 1 ) break;
313 }
314 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
315 if ( j == rowSize && ncnt == 1 )
316 {
317 constrListAux[nSlaves] = searchIndex;
318 constrList_[nSlaves++] = i;
319 }
320 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE2 )
321 {
322 if ( j == rowSize && ncnt == 1 )
323 printf("%4d : SlideReductionA - slave candidate %d = %d(%d)\n",
324 mypid_, nSlaves-1, i, constrListAux[nSlaves-1]);
325 }
326 }
327 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
328 {
329 printf("%4d : SlideReductionA - nSlave Candidate, nConstr = %d %d\n",
330 mypid_,nSlaves, nConstraints_);
331 }
332 }
333 else if ( constrList_ != NULL )
334 {
335 if ( mypid_ == 0 && ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 ) )
336 printf("%4d : SlideReductionA WARNING - constraint list not empty\n",
337 mypid_);
338 }
339
340 //---------------------------------------------------------------------
341 // search the constraint equations for the selected nodes
342 // (search for candidates column index with maximum magnitude)
343 //---------------------------------------------------------------------
344
345 nSelected = 0;
346 rowIndex = -1;
347 searchIndex = 0;
348
349 for ( i = EndRow-nConstraints_+1; i <= EndRow; i++ )
350 {
351 ierr = HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
352 hypre_assert(!ierr);
353 searchIndex = -1;
354 searchValue = -1.0E10;
355 for (j = 0; j < rowSize; j++)
356 {
357 if (colVal[j] != 0.0 && colInd[j] >= StartRow
358 && colInd[j] <= (EndRow-nConstraints_))
359 {
360 colIndex = hypre_BinarySearch(constrList_,colInd[j],nSlaves);
361 if ( colIndex >= 0 && constrListAux[colIndex] != -1)
362 {
363 if ( habs(colVal[j]) > searchValue )
364 {
365 if (i != constrListAux[colIndex])
366 {
367 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
368 {
369 printf("%4d : SlideReductionA WARNING - slave %d",
370 mypid_, colInd[j]);
371 printf(" candidate does not have constr %d\n", i);
372 }
373 }
374 searchValue = habs(colVal[j]);
375 searchIndex = colInd[j];
376 }
377 }
378 }
379 }
380 if ( searchIndex >= 0 )
381 {
382 selectedList[nSelected++] = searchIndex;
383 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE2 )
384 {
385 printf("%4d : SlideReductionA - constraint %4d <=> slave %d \n",
386 mypid_,i,searchIndex);
387 }
388 }
389 else
390 {
391 // get ready for error processing
392
393 colInd2 = new int[rowSize];
394 colVal2 = new double[rowSize];
395 for ( j = 0; j < rowSize; j++ )
396 {
397 colInd2[j] = colInd[j];
398 colVal2[j] = colVal[j];
399 }
400 rowIndex = i;
401 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
402 colInd = colInd2;
403 colVal = colVal2;
404 break;
405 }
406 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
407 }
408
409 //---------------------------------------------------------------------
410 // error processing
411 //---------------------------------------------------------------------
412
413 if ( searchIndex < 0 ) searchIndex = 1; else searchIndex = 0;
414 MPI_Allreduce(&searchIndex, &ncnt,1,MPI_INT,MPI_MAX,comm_);
415
416 if ( ncnt > 0 )
417 {
418 ncnt2 = 0;
419 while ( ncnt2 < numProcs_ )
420 {
421 if ( ncnt2 == mypid_ && rowIndex >= 0 )
422 {
423 printf("%4d : SlideReductionA ERROR - constraint number",mypid_);
424 printf(" cannot be found for row %d\n", rowIndex);
425 for (j = 0; j < rowSize; j++)
426 {
427 printf("ROW %4d COL = %d VAL = %e\n",rowIndex,colInd[j],
428 colVal[j]);
429 if (colVal[j] != 0.0 && colInd[j] >= StartRow
430 && colInd[j] <= (EndRow-nConstraints_))
431 {
432 colIndex = colInd[j];
433 HYPRE_ParCSRMatrixGetRow(A_csr,colIndex,&rowSize2,&colInd2,
434 &colVal2);
435 printf(" row %4d (%d) : \n",colIndex, rowSize2);
436 for (k = 0; k < rowSize2; k++)
437 printf(" row %4d col = %d val = %e\n",colIndex,
438 colInd2[k],colVal2[k]);
439 HYPRE_ParCSRMatrixRestoreRow(A_csr,colIndex,&rowSize2,
440 &colInd2,&colVal2);
441 }
442 }
443 printf("===================================================\n");
444 }
445 ncnt2++;
446 MPI_Barrier(comm_);
447 }
448 MPI_Finalize();
449 exit(1);
450 }
451 delete [] constrListAux;
452
453 //------------------------------------------------------------------
454 // sort the local selected node list and its auxiliary list, then
455 // form a global list of selected nodes on each processor
456 // form the corresponding auxiliary list for later pruning
457 //------------------------------------------------------------------
458
459 dble_array = new double[nSelected];
460 for ( i = 0; i < nSelected; i++ ) dble_array[i] = (double) i;
461 if ( nSelected > 1 ) hypre_qsort1(selectedList, dble_array, 0, nSelected-1);
462 for (i = 1; i < nSelected; i++)
463 {
464 if ( selectedList[i] == selectedList[i-1] )
465 {
466 printf("%4d : SlideReductionA ERROR - repeated selected nodes %d \n",
467 mypid_, selectedList[i]);
468 exit(1);
469 }
470 }
471 for (i = 0; i < nSelected; i++) selectedListAux[i] = (int) dble_array[i];
472 delete [] dble_array;
473
474 recvCntArray = new int[numProcs_];
475 displArray = new int[numProcs_];
476 MPI_Allgather(&nSelected, 1, MPI_INT, recvCntArray, 1,MPI_INT, comm_);
477 displArray[0] = 0;
478 for ( i = 1; i < numProcs_; i++ )
479 displArray[i] = displArray[i-1] + recvCntArray[i-1];
480 for ( i = 0; i < nSelected; i++ )
481 selectedListAux[i] += displArray[mypid_];
482 MPI_Allgatherv(selectedList, nSelected, MPI_INT, globalSelectedList,
483 recvCntArray, displArray, MPI_INT, comm_);
484 MPI_Allgatherv(selectedListAux, nSelected, MPI_INT, globalSelectedListAux,
485 recvCntArray, displArray, MPI_INT, comm_);
486 for ( i = 0; i < nSelected; i++ )
487 selectedListAux[i] -= displArray[mypid_];
488 delete [] recvCntArray;
489 delete [] displArray;
490
491 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE2 )
492 {
493 for ( i = 0; i < nSelected; i++ )
494 printf("%4d : SlideReductionA - selectedList %d = %d(%d)\n",mypid_,
495 i,selectedList[i],selectedListAux[i]);
496 }
497 }
498
499 //*****************************************************************************
500 // Part B of buildSlideReducedSystem : create submatrices
501 //-----------------------------------------------------------------------------
502
buildSlideReducedSystemPartB(int * ProcNRows,int * ProcNConstr,int globalNRows,int globalNConstr,int * globalSelectedList,int * globalSelectedListAux,HYPRE_ParCSRMatrix * rap_csr)503 void HYPRE_LinSysCore::buildSlideReducedSystemPartB(int *ProcNRows,
504 int *ProcNConstr, int globalNRows,
505 int globalNConstr, int *globalSelectedList,
506 int *globalSelectedListAux,
507 HYPRE_ParCSRMatrix *rap_csr)
508 {
509 int A21NRows, A21GlobalNRows, A21NCols, A21GlobalNCols, A21StartRow;
510 int i, j, nRows, ncnt, StartRow, A21StartCol;
511 int ierr, rowCount, maxRowSize, newEndRow, *A21MatSize, EndRow;
512 int rowIndex, rowSize, *colInd, rowSize2, colIndex, searchIndex;
513 int nnzA21, *newColInd, diagCount, newRowSize, procIndex;
514 int *recvCntArray, *displArray, *invA22MatSize, one=1;
515 int invA22NRows, invA22GlobalNRows, invA22NCols, invA22GlobalNCols;
516 int nSelected, *selectedList, *selectedListAux, globalNSelected;
517 double *colVal, *newColVal, *diagonal, *extDiagonal;
518 HYPRE_ParCSRMatrix A_csr, A21_csr, invA22_csr, RAP_csr;
519 HYPRE_IJMatrix A21, invA22;
520
521 //******************************************************************
522 // get matrix and constraint information
523 //------------------------------------------------------------------
524
525 nRows = localEndRow_ - localStartRow_ + 1;
526 StartRow = localStartRow_ - 1;
527 EndRow = localEndRow_ - 1;
528 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
529 globalNSelected = globalNConstr;
530 nSelected = nConstraints_;
531 selectedList = selectedList_;
532 selectedListAux = selectedListAux_;
533
534 //------------------------------------------------------------------
535 // calculate the dimension of A21
536 //------------------------------------------------------------------
537
538 A21NRows = 2 * nConstraints_;
539 A21NCols = nRows - 2 * nConstraints_;
540 A21GlobalNRows = 2 * globalNConstr;
541 A21GlobalNCols = globalNRows - 2 * globalNConstr;
542 A21StartRow = 2 * ProcNConstr[mypid_];
543 A21StartCol = ProcNRows[mypid_] - 2 * ProcNConstr[mypid_];
544
545 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
546 {
547 printf("%4d : SlideReductionB - A21StartRow = %d\n", mypid_,
548 A21StartRow);
549 printf("%4d : SlideReductionB - A21GlobalDim = %d %d\n", mypid_,
550 A21GlobalNRows, A21GlobalNCols);
551 printf("%4d : SlideReductionB - A21LocalDim = %d %d\n",mypid_,
552 A21NRows, A21NCols);
553 }
554
555 //------------------------------------------------------------------
556 // create a matrix context for A21
557 //------------------------------------------------------------------
558
559 ierr = HYPRE_IJMatrixCreate(comm_, A21StartRow, A21StartRow+A21NRows-1,
560 A21StartCol, A21StartCol+A21NCols-1, &A21);
561 ierr += HYPRE_IJMatrixSetObjectType(A21, HYPRE_PARCSR);
562 hypre_assert(!ierr);
563
564 //------------------------------------------------------------------
565 // compute the number of nonzeros in the first nConstraint row of A21
566 // (which consists of the rows in selectedList), the nnz will
567 // be reduced by excluding the constraint and selected slave columns
568 //------------------------------------------------------------------
569
570 rowCount = 0;
571 maxRowSize = 0;
572 newEndRow = EndRow - nConstraints_;
573 A21MatSize = new int[A21NRows];
574
575 for ( i = 0; i < nSelected; i++ )
576 {
577 for ( j = 0; j < nSelected; j++ )
578 {
579 if ( selectedListAux[j] == i )
580 {
581 rowIndex = selectedList[j];
582 break;
583 }
584 }
585 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
586 rowSize2 = 0;
587 for (j = 0; j < rowSize; j++)
588 {
589 colIndex = colInd[j];
590 if ( colVal[j] != 0.0 )
591 {
592 searchIndex = hypre_BinarySearch(globalSelectedList,colIndex,
593 globalNSelected);
594 if (searchIndex < 0 &&
595 (colIndex <= newEndRow || colIndex >= localEndRow_))
596 rowSize2++;
597 }
598 }
599 A21MatSize[rowCount] = rowSize2;
600 maxRowSize = ( rowSize2 > maxRowSize ) ? rowSize2 : maxRowSize;
601 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
602 rowCount++;
603 }
604
605 //------------------------------------------------------------------
606 // compute the number of nonzeros in the second nConstraint row of A21
607 // (which consists of the rows in constraint equations), the nnz will
608 // be reduced by excluding the selected slave columns only (since the
609 // entries corresponding to the constraint columns are 0, and since
610 // the selected matrix is a diagonal matrix, there is no need to
611 // search for slave equations in the off-processor list)
612 //------------------------------------------------------------------
613
614 rowCount = nSelected;
615 for ( i = EndRow-nConstraints_+1; i <= EndRow; i++ )
616 {
617 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
618 rowSize2 = 0;
619 for (j = 0; j < rowSize; j++)
620 {
621 if ( colVal[j] != 0.0 )
622 {
623 colIndex = colInd[j];
624 searchIndex = hypre_BinarySearch(globalSelectedList,colIndex,
625 globalNSelected);
626 if ( searchIndex < 0 ) rowSize2++;
627 }
628 }
629 A21MatSize[rowCount] = rowSize2;
630 maxRowSize = ( rowSize2 > maxRowSize ) ? rowSize2 : maxRowSize;
631 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
632 rowCount++;
633 }
634 nnzA21 = 0;
635 for ( i = 0; i < 2*nConstraints_; i++ ) nnzA21 += A21MatSize[i];
636
637 //------------------------------------------------------------------
638 // after fetching the row sizes, set up A21 with such sizes
639 //------------------------------------------------------------------
640
641 ierr = HYPRE_IJMatrixSetRowSizes(A21, A21MatSize);
642 ierr += HYPRE_IJMatrixInitialize(A21);
643 hypre_assert(!ierr);
644 delete [] A21MatSize;
645
646 //------------------------------------------------------------------
647 // next load the first nConstraint row to A21 extracted from A
648 // (at the same time, the D block is saved for future use)
649 //------------------------------------------------------------------
650
651 rowCount = A21StartRow;
652 if ( nConstraints_ > 0 ) diagonal = new double[nConstraints_];
653 else diagonal = NULL;
654 newColInd = new int[maxRowSize+1];
655 newColVal = new double[maxRowSize+1];
656
657 diagCount = 0;
658 for ( i = 0; i < nSelected; i++ )
659 {
660 for ( j = 0; j < nSelected; j++ )
661 {
662 if ( selectedListAux[j] == i )
663 {
664 rowIndex = selectedList[j];
665 break;
666 }
667 }
668 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
669 newRowSize = 0;
670 for (j = 0; j < rowSize; j++)
671 {
672 if ( colVal[j] != 0.0 )
673 {
674 colIndex = colInd[j];
675 if (colIndex <= newEndRow || colIndex >= localEndRow_)
676 {
677 searchIndex = HYPRE_LSI_Search(globalSelectedList,colIndex,
678 globalNSelected);
679 if ( searchIndex < 0 )
680 {
681 searchIndex = - searchIndex - 1;
682 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
683 if ( ProcNRows[procIndex] > colIndex ) break;
684 procIndex--;
685 colIndex = colInd[j]-ProcNConstr[procIndex]-searchIndex;
686 newColInd[newRowSize] = colIndex;
687 newColVal[newRowSize++] = colVal[j];
688 if ( colIndex < 0 || colIndex >= A21GlobalNCols )
689 {
690 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
691 {
692 printf("%4d : SlideReductionB WARNING - A21 ",mypid_);
693 printf("out of range (%d,%d (%d))\n", rowCount,
694 colIndex, A21GlobalNCols);
695 }
696 }
697 if ( newRowSize > maxRowSize+1 )
698 {
699 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
700 {
701 printf("%4d : SlideReductionB WARNING - ",mypid_);
702 printf("passing array boundary(1).\n");
703 }
704 }
705 }
706 }
707 else if ( colIndex > newEndRow && colIndex <= EndRow )
708 {
709 if ( colVal[j] != 0.0 ) diagonal[diagCount++] = colVal[j];
710 if ( habs(colVal[j]) < 1.0E-8 )
711 {
712 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
713 {
714 printf("%4d : SlideReductionB WARNING - large entry ",
715 mypid_);
716 printf("in invA22\n");
717 }
718 }
719 }
720 }
721 }
722
723 HYPRE_IJMatrixSetValues(A21, 1, &newRowSize, (const int *) &rowCount,
724 (const int *) newColInd, (const double *) newColVal);
725 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
726 if ( diagCount != (i+1) )
727 {
728 printf("%4d : SlideReductionB ERROR (3) - %d %d.\n", mypid_,
729 diagCount,i+1);
730 exit(1);
731 }
732 rowCount++;
733 }
734
735 //------------------------------------------------------------------
736 // send the diagonal to each processor that needs them
737 //------------------------------------------------------------------
738
739 recvCntArray = new int[numProcs_];
740 displArray = new int[numProcs_];
741 MPI_Allgather(&diagCount, 1, MPI_INT, recvCntArray, 1, MPI_INT, comm_);
742 displArray[0] = 0;
743 for ( i = 1; i < numProcs_; i++ )
744 displArray[i] = displArray[i-1] + recvCntArray[i-1];
745 ncnt = displArray[numProcs_-1] + recvCntArray[numProcs_-1];
746 if ( ncnt > 0 ) extDiagonal = new double[ncnt];
747 else extDiagonal = NULL;
748 MPI_Allgatherv(diagonal, diagCount, MPI_DOUBLE, extDiagonal,
749 recvCntArray, displArray, MPI_DOUBLE, comm_);
750 diagCount = ncnt;
751 delete [] recvCntArray;
752 delete [] displArray;
753 if ( diagonal != NULL ) delete [] diagonal;
754
755 //------------------------------------------------------------------
756 // next load the second nConstraint rows to A21 extracted from A
757 //------------------------------------------------------------------
758
759 for ( i = EndRow-nConstraints_+1; i <= EndRow; i++ )
760 {
761 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
762 newRowSize = 0;
763 for (j = 0; j < rowSize; j++)
764 {
765 colIndex = colInd[j];
766 searchIndex = HYPRE_LSI_Search(globalSelectedList,colIndex,
767 globalNSelected);
768 if ( searchIndex < 0 && colVal[j] != 0.0 )
769 {
770 searchIndex = - searchIndex - 1;
771 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
772 if ( ProcNRows[procIndex] > colIndex ) break;
773 procIndex--;
774 colIndex = colInd[j] - ProcNConstr[procIndex] - searchIndex;
775 newColInd[newRowSize] = colIndex;
776 newColVal[newRowSize++] = colVal[j];
777 if ( colIndex < 0 || colIndex >= A21GlobalNCols )
778 {
779 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
780 printf("%4d : SlideReductionB WARNING - A21(%d,%d(%d))\n",
781 mypid_, rowCount, colIndex, A21GlobalNCols);
782 }
783 if ( newRowSize > maxRowSize+1 )
784 {
785 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
786 {
787 printf("%4d : SlideReductionB WARNING - ",mypid_);
788 printf("passing array boundary(2).\n");
789 }
790 }
791 }
792 }
793 HYPRE_IJMatrixSetValues(A21, 1, &newRowSize, (const int *) &rowCount,
794 (const int *) newColInd, (const double *) newColVal);
795 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
796 rowCount++;
797 }
798 delete [] newColInd;
799 delete [] newColVal;
800
801 //------------------------------------------------------------------
802 // finally assemble the matrix and sanitize
803 //------------------------------------------------------------------
804
805 HYPRE_IJMatrixAssemble(A21);
806 HYPRE_IJMatrixGetObject(A21, (void **) &A21_csr);
807 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) A21_csr);
808
809 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
810 {
811 ncnt = 0;
812 MPI_Barrier(comm_);
813 while ( ncnt < numProcs_ )
814 {
815 if ( mypid_ == ncnt )
816 {
817 printf("====================================================\n");
818 printf("%4d : SlideReductionB - matrix A21 assembled %d.\n",
819 mypid_,A21StartRow);
820 fflush(stdout);
821 for ( i = A21StartRow; i < A21StartRow+2*nConstraints_; i++ )
822 {
823 HYPRE_ParCSRMatrixGetRow(A21_csr,i,&rowSize,&colInd,&colVal);
824 printf("A21 ROW = %6d (%d)\n", i, rowSize);
825 for ( j = 0; j < rowSize; j++ )
826 printf(" col = %6d, val = %e \n", colInd[j], colVal[j]);
827 HYPRE_ParCSRMatrixRestoreRow(A21_csr, i, &rowSize,
828 &colInd, &colVal);
829 }
830 printf("====================================================\n");
831 }
832 ncnt++;
833 MPI_Barrier(comm_);
834 }
835 }
836
837 //******************************************************************
838 // construct invA22
839 //------------------------------------------------------------------
840
841 //------------------------------------------------------------------
842 // calculate the dimension of invA22
843 //------------------------------------------------------------------
844
845 invA22NRows = A21NRows;
846 invA22NCols = invA22NRows;
847 invA22GlobalNRows = A21GlobalNRows;
848 invA22GlobalNCols = invA22GlobalNRows;
849 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
850 {
851 printf("%4d : SlideReductionB - A22GlobalDim = %d %d\n", mypid_,
852 invA22GlobalNRows, invA22GlobalNCols);
853 printf("%4d : SlideReductionB - A22LocalDim = %d %d\n", mypid_,
854 invA22NRows, invA22NCols);
855 }
856
857 //------------------------------------------------------------------
858 // create a matrix context for A22
859 //------------------------------------------------------------------
860
861 ierr = HYPRE_IJMatrixCreate(comm_, A21StartRow, A21StartRow+invA22NRows-1,
862 A21StartRow, A21StartRow+invA22NCols-1, &invA22);
863 ierr += HYPRE_IJMatrixSetObjectType(invA22, HYPRE_PARCSR);
864 hypre_assert(!ierr);
865
866 //------------------------------------------------------------------
867 // compute the no. of nonzeros in the first nConstraint row of invA22
868 //------------------------------------------------------------------
869
870 maxRowSize = 0;
871 invA22MatSize = new int[invA22NRows];
872 for ( i = 0; i < nConstraints_; i++ ) invA22MatSize[i] = 1;
873
874 //------------------------------------------------------------------
875 // compute the number of nonzeros in the second nConstraints row of
876 // invA22 (consisting of [D and A22 block])
877 //------------------------------------------------------------------
878
879 for ( i = 0; i < nSelected; i++ )
880 {
881 for ( j = 0; j < nSelected; j++ )
882 {
883 if ( selectedListAux[j] == i )
884 {
885 rowIndex = selectedList[j];
886 break;
887 }
888 }
889 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
890 rowSize2 = 1;
891 for (j = 0; j < rowSize; j++)
892 {
893 colIndex = colInd[j];
894 if ( colVal[j] != 0.0 )
895 {
896 if ( colIndex >= StartRow && colIndex <= newEndRow )
897 {
898 searchIndex = hypre_BinarySearch(selectedList, colIndex,
899 nSelected);
900 if ( searchIndex >= 0 ) rowSize2++;
901 }
902 else if ( colIndex < StartRow || colIndex > EndRow )
903 {
904 searchIndex = hypre_BinarySearch(globalSelectedList,colIndex,
905 globalNSelected);
906 if ( searchIndex >= 0 ) rowSize2++;
907 }
908 }
909 }
910 invA22MatSize[nConstraints_+i] = rowSize2;
911 maxRowSize = ( rowSize2 > maxRowSize ) ? rowSize2 : maxRowSize;
912 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
913 }
914
915 //------------------------------------------------------------------
916 // after fetching the row sizes, set up invA22 with such sizes
917 //------------------------------------------------------------------
918
919 ierr = HYPRE_IJMatrixSetRowSizes(invA22, invA22MatSize);
920 ierr += HYPRE_IJMatrixInitialize(invA22);
921 hypre_assert(!ierr);
922 delete [] invA22MatSize;
923
924 //------------------------------------------------------------------
925 // next load the first nConstraints_ row to invA22 extracted from A
926 // (that is, the D block)
927 //------------------------------------------------------------------
928
929 maxRowSize++;
930 newColInd = new int[maxRowSize];
931 newColVal = new double[maxRowSize];
932
933 for ( i = 0; i < diagCount; i++ )
934 {
935 extDiagonal[i] = 1.0 / extDiagonal[i];
936 }
937 for ( i = 0; i < nConstraints_; i++ )
938 {
939 newColInd[0] = A21StartRow + nConstraints_ + i;
940 rowIndex = A21StartRow + i;
941 if ( newColInd[0] < 0 || newColInd[0] >= invA22GlobalNCols )
942 {
943 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
944 printf("%4d : SlideReductionB WARNING - A22 (%d, %d (%d))\n",
945 mypid_, rowIndex, newColInd[0], invA22GlobalNCols);
946 }
947 newColVal[0] = extDiagonal[A21StartRow/2+i];
948 ierr = HYPRE_IJMatrixSetValues(invA22, 1, &one, (const int *) &rowIndex,
949 (const int *) newColInd, (const double *) newColVal);
950 hypre_assert(!ierr);
951 }
952
953 //------------------------------------------------------------------
954 // next load the second nConstraints_ rows to A22 extracted from A
955 //------------------------------------------------------------------
956
957 for ( i = 0; i < nSelected; i++ )
958 {
959 for ( j = 0; j < nSelected; j++ )
960 {
961 if ( selectedListAux[j] == i )
962 {
963 rowIndex = selectedList[j];
964 break;
965 }
966 }
967 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
968 newRowSize = 1;
969 newColInd[0] = A21StartRow + i;
970 newColVal[0] = extDiagonal[A21StartRow/2+i];
971 for (j = 0; j < rowSize; j++)
972 {
973 colIndex = colInd[j];
974 if ( colVal[j] != 0.0 )
975 {
976 searchIndex = hypre_BinarySearch(globalSelectedList,colIndex,
977 globalNSelected);
978 if ( searchIndex >= 0 )
979 {
980 searchIndex = globalSelectedListAux[searchIndex];
981 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
982 if ( ProcNRows[procIndex] > colIndex ) break;
983 if ( procIndex == numProcs_ )
984 {
985 newColInd[newRowSize] = searchIndex + globalNConstr;
986 }
987 else
988 {
989 newColInd[newRowSize] = searchIndex +
990 ProcNConstr[procIndex];
991 }
992 if ( newColInd[newRowSize] < 0 ||
993 newColInd[newRowSize] >= invA22GlobalNCols )
994 {
995 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
996 printf("%4d : SlideReductionB WARNING - A22(%d,%d,%d)\n",
997 mypid_,rowCount,newColInd[newRowSize],
998 invA22GlobalNCols);
999 }
1000 newColVal[newRowSize++] = - extDiagonal[A21StartRow/2+i] *
1001 colVal[j] * extDiagonal[searchIndex];
1002 if ( newRowSize > maxRowSize )
1003 {
1004 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1005 {
1006 printf("%4d : SlideReductionB WARNING - ",mypid_);
1007 printf("passing array boundary(3).\n");
1008 }
1009 }
1010 }
1011 }
1012 }
1013 rowCount = A21StartRow + nConstraints_ + i;
1014 ierr = HYPRE_IJMatrixSetValues(invA22, 1, &newRowSize,
1015 (const int *) &rowCount, (const int *) newColInd,
1016 (const double *) newColVal);
1017 hypre_assert(!ierr);
1018 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
1019 }
1020 delete [] newColInd;
1021 delete [] newColVal;
1022 delete [] extDiagonal;
1023
1024 //------------------------------------------------------------------
1025 // finally assemble the matrix and sanitize
1026 //------------------------------------------------------------------
1027
1028 HYPRE_IJMatrixAssemble(invA22);
1029 HYPRE_IJMatrixGetObject(invA22, (void **) &invA22_csr);
1030 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) invA22_csr);
1031
1032 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
1033 {
1034 ncnt = 0;
1035 MPI_Barrier(comm_);
1036 while ( ncnt < numProcs_ )
1037 {
1038 if ( mypid_ == ncnt )
1039 {
1040 printf("====================================================\n");
1041 printf("%4d : SlideReductionB - invA22 \n", mypid_);
1042 for ( i = A21StartRow; i < A21StartRow+2*nConstraints_; i++ )
1043 {
1044 HYPRE_ParCSRMatrixGetRow(invA22_csr,i,&rowSize,&colInd,
1045 &colVal);
1046 printf("invA22 ROW = %6d (%d)\n", i, rowSize);
1047 for ( j = 0; j < rowSize; j++ )
1048 printf(" col = %6d, val = %e \n", colInd[j], colVal[j]);
1049 HYPRE_ParCSRMatrixRestoreRow(invA22_csr,i,&rowSize,&colInd,
1050 &colVal);
1051 }
1052 printf("====================================================\n");
1053 }
1054 MPI_Barrier(comm_);
1055 ncnt++;
1056 }
1057 }
1058
1059 //******************************************************************
1060 // perform the triple matrix product
1061 //------------------------------------------------------------------
1062
1063 HYPRE_IJMatrixGetObject(A21, (void **) &A21_csr);
1064 HYPRE_IJMatrixGetObject(invA22, (void **) &invA22_csr);
1065 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1066 printf("%4d : SlideReductionB - Triple matrix product starts\n",mypid_);
1067
1068 hypre_BoomerAMGBuildCoarseOperator( (hypre_ParCSRMatrix *) A21_csr,
1069 (hypre_ParCSRMatrix *) invA22_csr,
1070 (hypre_ParCSRMatrix *) A21_csr,
1071 (hypre_ParCSRMatrix **) &RAP_csr);
1072
1073 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1074 printf("%4d : SlideReductionB - Triple matrix product ends\n", mypid_);
1075
1076 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
1077 {
1078 MPI_Barrier(comm_);
1079 ncnt = 0;
1080 while ( ncnt < numProcs_ )
1081 {
1082 if ( mypid_ == ncnt )
1083 {
1084 for ( i = A21StartRow; i < A21StartRow+A21NCols; i++ ) {
1085 HYPRE_ParCSRMatrixGetRow(RAP_csr,i,&rowSize,&colInd, &colVal);
1086 printf("RAP ROW = %6d (%d)\n", i, rowSize);
1087 for ( j = 0; j < rowSize; j++ )
1088 printf(" col = %6d, val = %e \n", colInd[j], colVal[j]);
1089 HYPRE_ParCSRMatrixRestoreRow(RAP_csr,i,&rowSize,&colInd,
1090 &colVal);
1091 }
1092 }
1093 MPI_Barrier(comm_);
1094 ncnt++;
1095 }
1096 }
1097
1098 //******************************************************************
1099 // set global objects and checking
1100 //------------------------------------------------------------------
1101
1102 HYA21_ = A21;
1103 HYinvA22_ = invA22;
1104 (*rap_csr) = RAP_csr;
1105
1106 MPI_Allreduce(&nnzA21,&ncnt,1,MPI_INT,MPI_SUM,comm_);
1107 if ( mypid_ == 0 && ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 ))
1108 printf(" SlideReductionB : NNZ of A21 = %d\n", ncnt);
1109 }
1110
1111 //*****************************************************************************
1112 // Part C of buildSlideReducedSystem : create subvectors and wrap up
1113 //-----------------------------------------------------------------------------
1114
buildSlideReducedSystemPartC(int * ProcNRows,int * ProcNConstr,int globalNRows,int globalNConstr,int * globalSelectedList,int * globalSelectedListAux,HYPRE_ParCSRMatrix RAP_csr)1115 void HYPRE_LinSysCore::buildSlideReducedSystemPartC(int *ProcNRows,
1116 int *ProcNConstr, int globalNRows,
1117 int globalNConstr, int *globalSelectedList,
1118 int *globalSelectedListAux,
1119 HYPRE_ParCSRMatrix RAP_csr)
1120 {
1121 int i, j, nRows, StartRow, EndRow;
1122 int newNRows, *reducedAMatSize, reducedAStartRow;
1123 int rowCount, rowIndex, newRowSize, rowSize, rowSize2, *newColInd;
1124 int *colInd, *colInd2, colIndex, searchIndex, ubound, ncnt, ierr;
1125 int A21NRows, A21NCols, A21GlobalNRows, A21GlobalNCols, A21StartRow;
1126 int A12NRows, A12NCols, A12GlobalNRows, A12GlobalNCols, A12StartRow;
1127 int *A12MatSize, newEndRow, globalNSelected, procIndex, nnzA12;
1128 int nSelected, *selectedList, *selectedListAux, A21StartCol;
1129 double *colVal, *colVal2, *newColVal, ddata;
1130 HYPRE_IJMatrix A12, reducedA;
1131 HYPRE_ParCSRMatrix A_csr, A12_csr, reducedA_csr, invA22_csr;
1132 HYPRE_IJVector f2, f2hat;
1133 HYPRE_ParVector f2_csr, f2hat_csr, reducedB_csr;
1134
1135 //------------------------------------------------------------------
1136 // get matrix and constraint information
1137 //------------------------------------------------------------------
1138
1139 nRows = localEndRow_ - localStartRow_ + 1;
1140 StartRow = localStartRow_ - 1;
1141 EndRow = localEndRow_ - 1;
1142 newEndRow = EndRow - nConstraints_;
1143 newEndRow = EndRow - nConstraints_;
1144 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
1145 globalNSelected = globalNConstr;
1146 nSelected = nConstraints_;
1147 selectedList = selectedList_;
1148 selectedListAux = selectedListAux_;
1149
1150 //------------------------------------------------------------------
1151 // first calculate the dimension of the reduced matrix
1152 //------------------------------------------------------------------
1153
1154 newNRows = nRows - 2 * nConstraints_;
1155 A21StartCol = ProcNRows[mypid_] - 2 * ProcNConstr[mypid_];
1156 ierr = HYPRE_IJMatrixCreate(comm_, A21StartCol, A21StartCol+newNRows-1,
1157 A21StartCol, A21StartCol+newNRows-1, &reducedA);
1158 ierr += HYPRE_IJMatrixSetObjectType(reducedA, HYPRE_PARCSR);
1159 hypre_assert(!ierr);
1160
1161 //------------------------------------------------------------------
1162 // set up reducedA with proper sizes
1163 //------------------------------------------------------------------
1164
1165 reducedAMatSize = new int[newNRows];
1166 reducedAStartRow = ProcNRows[mypid_] - 2 * ProcNConstr[mypid_];
1167 rowCount = reducedAStartRow;
1168 rowIndex = 0;
1169
1170 for ( i = StartRow; i <= newEndRow; i++ )
1171 {
1172 searchIndex = hypre_BinarySearch(selectedList, i, nSelected);
1173 if ( searchIndex < 0 )
1174 {
1175 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
1176 ierr = HYPRE_ParCSRMatrixGetRow(RAP_csr,rowCount,&rowSize2,
1177 &colInd2, &colVal2);
1178 hypre_assert( !ierr );
1179 newRowSize = rowSize + rowSize2;
1180 newColInd = new int[newRowSize];
1181 for (j = 0; j < rowSize; j++) newColInd[j] = colInd[j];
1182 for (j = 0; j < rowSize2; j++) newColInd[rowSize+j] = colInd2[j];
1183 hypre_qsort0(newColInd, 0, newRowSize-1);
1184 ncnt = 0;
1185 for ( j = 1; j < newRowSize; j++ )
1186 {
1187 if ( newColInd[j] != newColInd[ncnt] )
1188 {
1189 ncnt++;
1190 newColInd[ncnt] = newColInd[j];
1191 }
1192 }
1193 if ( newRowSize > 0 ) ncnt++;
1194 reducedAMatSize[rowIndex++] = ncnt;
1195
1196 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
1197 ierr = HYPRE_ParCSRMatrixRestoreRow(RAP_csr,rowCount,&rowSize2,
1198 &colInd2,&colVal2);
1199 delete [] newColInd;
1200 hypre_assert( !ierr );
1201 rowCount++;
1202 }
1203 }
1204
1205 //------------------------------------------------------------------
1206 // create a matrix context for reducedA
1207 //------------------------------------------------------------------
1208
1209 ierr = HYPRE_IJMatrixSetRowSizes(reducedA, reducedAMatSize);
1210 ierr += HYPRE_IJMatrixInitialize(reducedA);
1211 hypre_assert(!ierr);
1212 delete [] reducedAMatSize;
1213
1214 //------------------------------------------------------------------
1215 // load the reducedA matrix
1216 //------------------------------------------------------------------
1217
1218 rowCount = reducedAStartRow;
1219 for ( i = StartRow; i <= newEndRow; i++ )
1220 {
1221 searchIndex = hypre_BinarySearch(selectedList, i, nSelected);
1222 if ( searchIndex < 0 )
1223 {
1224 HYPRE_ParCSRMatrixGetRow(A_csr, i, &rowSize, &colInd, &colVal);
1225 HYPRE_ParCSRMatrixGetRow(RAP_csr,rowCount,&rowSize2,&colInd2,
1226 &colVal2);
1227 newRowSize = rowSize + rowSize2;
1228 newColInd = new int[newRowSize];
1229 newColVal = new double[newRowSize];
1230 ncnt = 0;
1231
1232 for ( j = 0; j < rowSize; j++ )
1233 {
1234 colIndex = colInd[j];
1235 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
1236 if ( ProcNRows[procIndex] > colIndex ) break;
1237 if ( procIndex == numProcs_ )
1238 ubound = globalNRows-(globalNConstr-ProcNConstr[numProcs_-1]);
1239 else
1240 ubound = ProcNRows[procIndex] -
1241 (ProcNConstr[procIndex]-ProcNConstr[procIndex-1]);
1242 procIndex--;
1243 if ( colIndex < ubound )
1244 {
1245 searchIndex = HYPRE_LSI_Search(globalSelectedList,colIndex,
1246 globalNSelected);
1247 if ( searchIndex < 0 )
1248 {
1249 searchIndex = - searchIndex - 1;
1250 newColInd[ncnt] = colIndex - ProcNConstr[procIndex] -
1251 searchIndex;
1252 newColVal[ncnt++] = colVal[j];
1253 }
1254 }
1255 }
1256 for ( j = 0; j < rowSize2; j++ )
1257 {
1258 newColInd[ncnt+j] = colInd2[j];
1259 newColVal[ncnt+j] = - colVal2[j];
1260 }
1261 newRowSize = ncnt + rowSize2;
1262 hypre_qsort1(newColInd, newColVal, 0, newRowSize-1);
1263 ncnt = 0;
1264 for ( j = 0; j < newRowSize; j++ )
1265 {
1266 if ( j != ncnt && newColInd[j] == newColInd[ncnt] )
1267 newColVal[ncnt] += newColVal[j];
1268 else if ( newColInd[j] != newColInd[ncnt] )
1269 {
1270 ncnt++;
1271 newColVal[ncnt] = newColVal[j];
1272 newColInd[ncnt] = newColInd[j];
1273 }
1274 }
1275 newRowSize = ncnt + 1;
1276 ierr = HYPRE_IJMatrixSetValues(reducedA, 1, &newRowSize,
1277 (const int *) &rowCount, (const int *) newColInd,
1278 (const double *) newColVal);
1279 hypre_assert(!ierr);
1280 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
1281 HYPRE_ParCSRMatrixRestoreRow(RAP_csr,rowCount,&rowSize2,&colInd2,
1282 &colVal2);
1283 rowCount++;
1284 delete [] newColInd;
1285 delete [] newColVal;
1286 }
1287 }
1288
1289 //------------------------------------------------------------------
1290 // assemble the reduced matrix
1291 //------------------------------------------------------------------
1292
1293 HYPRE_IJMatrixAssemble(reducedA);
1294 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1295 printf("%4d : SlideReductionC - reducedAStartRow = %d\n", mypid_,
1296 reducedAStartRow);
1297
1298 HYPRE_IJMatrixGetObject(reducedA, (void **) &reducedA_csr);
1299
1300 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
1301 {
1302 MPI_Barrier(comm_);
1303 ncnt = 0;
1304 while ( ncnt < numProcs_ )
1305 {
1306 if ( mypid_ == ncnt )
1307 {
1308 printf("====================================================\n");
1309 for ( i = reducedAStartRow;
1310 i < reducedAStartRow+nRows-2*nConstraints_; i++ )
1311 {
1312 printf("%d : reducedA ROW %d\n", mypid_, i);
1313 ierr = HYPRE_ParCSRMatrixGetRow(reducedA_csr,i,&rowSize,
1314 &colInd, &colVal);
1315 //hypre_qsort1(colInd, colVal, 0, rowSize-1);
1316 for ( j = 0; j < rowSize; j++ )
1317 if ( colVal[j] != 0.0 )
1318 printf("%4d %4d %20.13e\n", i, colInd[j], colVal[j]);
1319 HYPRE_ParCSRMatrixRestoreRow(reducedA_csr,i,&rowSize,&colInd,
1320 &colVal);
1321 }
1322 printf("====================================================\n");
1323 }
1324 MPI_Barrier(comm_);
1325 ncnt++;
1326 }
1327 }
1328
1329 // *****************************************************************
1330 // form modified right hand side (f1 = f1 - A12*invA22*f2)
1331 // *****************************************************************
1332
1333 // *****************************************************************
1334 // form f2hat = invA22 * f2
1335 //------------------------------------------------------------------
1336
1337 A21NRows = 2 * nConstraints_;
1338 A21NCols = nRows - 2 * nConstraints_;
1339 A21GlobalNRows = 2 * globalNConstr;
1340 A21GlobalNCols = globalNRows - 2 * globalNConstr;
1341 A21StartRow = 2 * ProcNConstr[mypid_];
1342
1343 HYPRE_IJVectorCreate(comm_, A21StartRow, A21StartRow+A21NRows-1, &f2);
1344 HYPRE_IJVectorSetObjectType(f2, HYPRE_PARCSR);
1345 ierr += HYPRE_IJVectorInitialize(f2);
1346 ierr += HYPRE_IJVectorAssemble(f2);
1347 hypre_assert(!ierr);
1348
1349 HYPRE_IJVectorCreate(comm_, A21StartRow, A21StartRow+A21NRows-1, &f2hat);
1350 HYPRE_IJVectorSetObjectType(f2hat, HYPRE_PARCSR);
1351 ierr += HYPRE_IJVectorInitialize(f2hat);
1352 ierr += HYPRE_IJVectorAssemble(f2hat);
1353 hypre_assert(!ierr);
1354
1355 colInd = new int[nSelected*2];
1356 colVal = new double[nSelected*2];
1357
1358 for ( i = 0; i < nSelected; i++ )
1359 {
1360 for ( j = 0; j < nSelected; j++ )
1361 {
1362 if ( selectedListAux[j] == i )
1363 {
1364 colInd[i] = selectedList[j];
1365 break;
1366 }
1367 }
1368 if ( colInd[i] < 0 )
1369 {
1370 printf("%4d : SlideReductionC ERROR - out of range %d\n", mypid_,
1371 colInd[i]);
1372 exit(1);
1373 }
1374 }
1375 for ( i = 0; i < nSelected; i++ )
1376 {
1377 colInd[nSelected+i] = EndRow - nConstraints_ + i + 1;
1378 }
1379 HYPRE_IJVectorGetValues(HYb_, 2*nSelected, colInd, colVal);
1380 for ( i = 0; i < nSelected*2; i++ ) colInd[i] = A21StartRow + i;
1381 ierr = HYPRE_IJVectorSetValues(f2,2*nSelected,(const int *) colInd,
1382 (const double *) colVal);
1383 hypre_assert( !ierr );
1384 HYPRE_IJVectorGetObject(f2, (void **) &f2_csr);
1385 HYPRE_IJVectorGetObject(f2hat, (void **) &f2hat_csr);
1386 HYPRE_IJMatrixGetObject(HYinvA22_, (void **) &invA22_csr);
1387 HYPRE_ParCSRMatrixMatvec( 1.0, invA22_csr, f2_csr, 0.0, f2hat_csr );
1388 delete [] colVal;
1389 delete [] colInd;
1390 HYPRE_IJVectorDestroy(f2);
1391
1392 // *****************************************************************
1393 // set up A12 with proper sizes before forming f2til = A12 * f2hat
1394 //------------------------------------------------------------------
1395
1396 //------------------------------------------------------------------
1397 // calculate the dimension of A12
1398 //------------------------------------------------------------------
1399
1400 A12NRows = A21NCols;
1401 A12NCols = A21NRows;
1402 A12GlobalNRows = A21GlobalNCols;
1403 A12GlobalNCols = A21GlobalNRows;
1404 A12MatSize = new int[A12NRows];
1405 A12StartRow = ProcNRows[mypid_] - 2 * ProcNConstr[mypid_];
1406 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1407 {
1408 printf("%4d : SlideReductionC - A12GlobalDim = %d %d\n", mypid_,
1409 A12GlobalNRows, A12GlobalNCols);
1410 printf("%4d : SlideReductionC - A12LocalDim = %d %d\n", mypid_,
1411 A12NRows, A12NCols);
1412 }
1413
1414 //------------------------------------------------------------------
1415 // create a matrix context for A12
1416 //------------------------------------------------------------------
1417
1418 ierr = HYPRE_IJMatrixCreate(comm_, A21StartCol, A21StartCol+A12NRows-1,
1419 A21StartRow, A21StartRow+A12NCols-1, &A12);
1420 ierr += HYPRE_IJMatrixSetObjectType(A12, HYPRE_PARCSR);
1421 hypre_assert(!ierr);
1422
1423 //------------------------------------------------------------------
1424 // compute the number of nonzeros in each row of A12
1425 // (which consists of the rows in selectedList and the constraints)
1426 //------------------------------------------------------------------
1427
1428 rowCount = A12StartRow;
1429 rowIndex = 0;
1430
1431 for ( i = StartRow; i <= newEndRow; i++ )
1432 {
1433 searchIndex = hypre_BinarySearch(selectedList, i, nSelected);
1434 if ( searchIndex < 0 )
1435 {
1436 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
1437 newRowSize = 0;
1438 for (j = 0; j < rowSize; j++)
1439 {
1440 if ( colVal[j] != 0.0 )
1441 {
1442 colIndex = colInd[j];
1443 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
1444 if ( ProcNRows[procIndex] > colIndex ) break;
1445 if ( procIndex == numProcs_ )
1446 ubound = globalNRows -
1447 (globalNConstr - ProcNConstr[numProcs_-1]);
1448 else
1449 ubound = ProcNRows[procIndex] -
1450 (ProcNConstr[procIndex]-ProcNConstr[procIndex-1]);
1451 procIndex--;
1452 if ( colIndex >= ubound ) newRowSize++;
1453 else
1454 {
1455 if (hypre_BinarySearch(globalSelectedList,colIndex,
1456 globalNSelected) >= 0)
1457 newRowSize++;
1458 }
1459 }
1460 }
1461 A12MatSize[rowIndex++] = newRowSize;
1462 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
1463 rowCount++;
1464 }
1465 }
1466
1467 //------------------------------------------------------------------
1468 // after fetching the row sizes, set up A12 with such sizes
1469 //------------------------------------------------------------------
1470
1471 nnzA12 = 0;
1472 for ( i = 0; i < A12NRows; i++ ) nnzA12 += A12MatSize[i];
1473 ierr = HYPRE_IJMatrixSetRowSizes(A12, A12MatSize);
1474 ierr += HYPRE_IJMatrixInitialize(A12);
1475 hypre_assert(!ierr);
1476 delete [] A12MatSize;
1477
1478 //------------------------------------------------------------------
1479 // load the A12 matrix
1480 //------------------------------------------------------------------
1481
1482 rowCount = A12StartRow;
1483 for ( i = StartRow; i <= newEndRow; i++ )
1484 {
1485 searchIndex = hypre_BinarySearch(selectedList, i, nSelected);
1486 if ( searchIndex < 0 )
1487 {
1488 HYPRE_ParCSRMatrixGetRow(A_csr, i, &rowSize, &colInd, &colVal);
1489 newRowSize = 0;
1490 newColInd = new int[rowSize];
1491 newColVal = new double[rowSize];
1492 for (j = 0; j < rowSize; j++)
1493 {
1494 colIndex = colInd[j];
1495 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
1496 if ( ProcNRows[procIndex] > colIndex ) break;
1497 if ( procIndex == numProcs_ )
1498 ubound = globalNRows-(globalNConstr-ProcNConstr[numProcs_-1]);
1499 else
1500 ubound = ProcNRows[procIndex] -
1501 (ProcNConstr[procIndex]-ProcNConstr[procIndex-1]);
1502 procIndex--;
1503 if ( colIndex >= ubound ) {
1504 if ( procIndex != numProcs_ - 1 )
1505 {
1506 newColInd[newRowSize] = colInd[j] - ubound +
1507 ProcNConstr[procIndex] +
1508 ProcNConstr[procIndex+1];
1509 }
1510 else
1511 {
1512 newColInd[newRowSize] = colInd[j] - ubound +
1513 ProcNConstr[procIndex] +
1514 globalNConstr;
1515 }
1516 if ( newColInd[newRowSize] < 0 ||
1517 newColInd[newRowSize] >= A12GlobalNCols )
1518 {
1519 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1520 {
1521 printf("%4d : SlideReductionC WARNING - A12 col index ",
1522 mypid_);
1523 printf("out of range %d %d(%d)\n", i,
1524 newColInd[newRowSize], A12GlobalNCols);
1525 }
1526 }
1527 newColVal[newRowSize++] = colVal[j];
1528 } else
1529 {
1530 searchIndex = HYPRE_LSI_Search(globalSelectedList,colInd[j],
1531 globalNSelected);
1532 if ( searchIndex >= 0)
1533 {
1534 searchIndex = globalSelectedListAux[searchIndex];
1535 newColInd[newRowSize] = searchIndex +
1536 ProcNConstr[procIndex];
1537 if ( newColInd[newRowSize] < 0 ||
1538 newColInd[newRowSize] >= A12GlobalNCols )
1539 {
1540 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1541 {
1542 printf("%4d : SlideReductionC WARNING - \n",mypid_);
1543 printf(" A12(%d,%d,%d))\n", i,
1544 newColInd[newRowSize], A12GlobalNCols);
1545 }
1546 }
1547 newColVal[newRowSize++] = colVal[j];
1548 }
1549 }
1550 }
1551 ierr = HYPRE_IJMatrixSetValues(A12, 1, &newRowSize,
1552 (const int *) &rowCount, (const int *) newColInd,
1553 (const double *) newColVal);
1554 hypre_assert(!ierr);
1555 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
1556
1557 rowCount++;
1558 delete [] newColInd;
1559 delete [] newColVal;
1560 }
1561 }
1562 MPI_Allreduce(&nnzA12,&ncnt,1,MPI_INT,MPI_SUM,comm_);
1563 if ( mypid_ == 0 && ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 ))
1564 printf(" SlideReductionC : NNZ of A12 = %d\n", ncnt);
1565
1566 //------------------------------------------------------------------
1567 // assemble the A12 matrix
1568 //------------------------------------------------------------------
1569
1570 ierr = HYPRE_IJMatrixAssemble(A12);
1571 hypre_assert( !ierr );
1572 HYPRE_IJMatrixGetObject(A12, (void **) &A12_csr);
1573
1574 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
1575 {
1576 MPI_Barrier(comm_);
1577 ncnt = 0;
1578 while ( ncnt < numProcs_ )
1579 {
1580 if ( mypid_ == ncnt )
1581 {
1582 printf("====================================================\n");
1583 for (i=A12StartRow;i<A12StartRow+A12NRows;i++)
1584 {
1585 printf("%d : A12 ROW %d\n", mypid_, i);
1586 HYPRE_ParCSRMatrixGetRow(A12_csr,i,&rowSize,&colInd,&colVal);
1587 //hypre_qsort1(colInd, colVal, 0, rowSize-1);
1588 for ( j = 0; j < rowSize; j++ )
1589 if ( colVal[j] != 0.0 )
1590 printf(" A12 %d %d %20.13e\n", i, colInd[j], colVal[j]);
1591 HYPRE_ParCSRMatrixRestoreRow(A12_csr,i,&rowSize,&colInd,
1592 &colVal);
1593 }
1594 printf("====================================================\n");
1595 }
1596 MPI_Barrier(comm_);
1597 ncnt++;
1598 }
1599 }
1600
1601 //------------------------------------------------------------------
1602 // form reducedB_ = A12 * f2hat
1603 //------------------------------------------------------------------
1604
1605 ierr = HYPRE_IJVectorCreate(comm_, reducedAStartRow,
1606 reducedAStartRow+newNRows-1, &reducedB_);
1607 ierr += HYPRE_IJVectorSetObjectType(reducedB_, HYPRE_PARCSR);
1608 ierr += HYPRE_IJVectorInitialize(reducedB_);
1609 ierr += HYPRE_IJVectorAssemble(reducedB_);
1610 hypre_assert( !ierr );
1611
1612 HYPRE_IJVectorGetObject(reducedB_, (void **) &reducedB_csr);
1613 HYPRE_ParCSRMatrixMatvec( -1.0, A12_csr, f2hat_csr, 0.0, reducedB_csr );
1614 HYPRE_IJMatrixDestroy(A12);
1615 HYPRE_IJVectorDestroy(f2hat);
1616
1617 //------------------------------------------------------------------
1618 // finally form reducedB = f1 - f2til
1619 //------------------------------------------------------------------
1620
1621 rowCount = reducedAStartRow;
1622 for ( i = StartRow; i <= newEndRow; i++ )
1623 {
1624 if ( hypre_BinarySearch(selectedList, i, nSelected) < 0 )
1625 {
1626 HYPRE_IJVectorGetValues(HYb_, 1, &i, &ddata);
1627 HYPRE_IJVectorAddToValues(reducedB_, 1, (const int *) &rowCount,
1628 (const double *) &ddata);
1629 rowCount++;
1630 }
1631 }
1632
1633 //******************************************************************
1634 // set up the system with the new matrix
1635 //------------------------------------------------------------------
1636
1637 reducedA_ = reducedA;
1638 ierr = HYPRE_IJVectorCreate(comm_, reducedAStartRow,
1639 reducedAStartRow+newNRows-1, &reducedX_);
1640 ierr = HYPRE_IJVectorSetObjectType(reducedX_, HYPRE_PARCSR);
1641 ierr = HYPRE_IJVectorInitialize(reducedX_);
1642 ierr = HYPRE_IJVectorAssemble(reducedX_);
1643 hypre_assert(!ierr);
1644
1645 ierr = HYPRE_IJVectorCreate(comm_, reducedAStartRow,
1646 reducedAStartRow+newNRows-1, &reducedR_);
1647 ierr = HYPRE_IJVectorSetObjectType(reducedR_, HYPRE_PARCSR);
1648 ierr = HYPRE_IJVectorInitialize(reducedR_);
1649 ierr = HYPRE_IJVectorAssemble(reducedR_);
1650 hypre_assert(!ierr);
1651 }
1652
1653 //*****************************************************************************
1654 // Given the matrix (A) within the object, compute the reduced system and put
1655 // it in place. Additional information given are :
1656 //
1657 // Additional assumptions are :
1658 //
1659 // - a given slave equation and the corresponding constraint equation
1660 // reside in the same processor
1661 // - constraint equations are given at the end of the local matrix
1662 // (hence given by EndRow_-nConstr to EndRow_)
1663 // - each processor gets a contiguous block of equations, and processor
1664 // i+1 has equation numbers higher than those of processor i
1665 //-----------------------------------------------------------------------------
1666 // This version replaces the selected slave equation with an identity row
1667 //-----------------------------------------------------------------------------
1668
buildSlideReducedSystem2()1669 void HYPRE_LinSysCore::buildSlideReducedSystem2()
1670 {
1671 int i, j, nRows, globalNRows, colIndex;
1672 int globalNConstr, globalNSelected, *globalSelectedList;
1673 int *globalSelectedListAux;
1674 int nSelected, *tempList, reducedAStartRow;
1675 int searchIndex, procIndex, A21StartRow, A12StartRow, A12NRows;
1676 int rowSize, *colInd, A21NRows, A21GlobalNRows;
1677 int A21NCols, A21GlobalNCols, rowCount, maxRowSize, newEndRow;
1678 int A12NCols, A12GlobalNCols;
1679 int *A21MatSize, rowIndex, *A12MatSize, A12GlobalNRows;
1680 int *newColInd, diagCount, newRowSize, ierr;
1681 int invA22NRows, invA22GlobalNRows, invA22NCols, invA22GlobalNCols;
1682 int *invA22MatSize, newNRows, newColIndex;
1683 int *colInd2, ncnt, ubound, one=1;
1684 int rowSize2, *recvCntArray, *displArray, ncnt2, isAConstr;
1685 int StartRow, EndRow, *reducedAMatSize;
1686 int *ProcNRows, *ProcNConstr, nnzA21, nnzA12;
1687 int A21StartCol;
1688 double *colVal, *colVal2, *newColVal, *diagonal;
1689 double *extDiagonal, ddata;
1690 HYPRE_IJMatrix A12, A21, invA22, reducedA;
1691 HYPRE_ParCSRMatrix A_csr, A12_csr, A21_csr, invA22_csr, RAP_csr;
1692 HYPRE_ParCSRMatrix reducedA_csr;
1693 HYPRE_IJVector f2, f2hat;
1694 HYPRE_ParVector f2_csr, f2hat_csr, reducedB_csr;
1695
1696 //******************************************************************
1697 // fetch local matrix information
1698 //------------------------------------------------------------------
1699
1700 if ( mypid_ == 0 && (HYOutputLevel_ & HYFEI_SLIDEREDUCE1) )
1701 printf("%4d : SlideReduction2 begins....\n",mypid_);
1702 StartRow = localStartRow_ - 1;
1703 EndRow = localEndRow_ - 1;
1704 nRows = localEndRow_ - localStartRow_ + 1;
1705 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1706 printf("%4d : SlideReduction2 : StartRow/EndRow = %d %d\n",mypid_,
1707 StartRow,EndRow);
1708
1709 //******************************************************************
1710 // construct local and global information about where the constraints
1711 // are (this is given by user or searched within this code)
1712 //------------------------------------------------------------------
1713
1714 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
1715
1716 //------------------------------------------------------------------
1717 // search the entire local matrix to find where the constraint
1718 // equations are, if not already given
1719 // (The constraint equations are assumed to be at the end of the
1720 // matrix) ==> nConstraints, globalNConstr
1721 //------------------------------------------------------------------
1722
1723 MPI_Allreduce(&nConstraints_,&globalNConstr,1,MPI_INT,MPI_SUM,comm_);
1724 if ( globalNConstr == 0 )
1725 {
1726 for ( i = EndRow; i >= StartRow; i-- )
1727 {
1728 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
1729 isAConstr = 1;
1730 for (j = 0; j < rowSize; j++)
1731 if (colInd[j] == i && colVal[j] != 0.0) {isAConstr = 0; break;}
1732 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
1733 if ( isAConstr ) nConstraints_++;
1734 else break;
1735 }
1736 }
1737 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1738 printf("%4d : SlideReduction2 : no. constr = %d\n", mypid_,
1739 nConstraints_);
1740
1741 MPI_Allreduce(&nConstraints_, &globalNConstr, 1, MPI_INT, MPI_SUM, comm_);
1742 if ( globalNConstr == 0 ) return;
1743
1744 //------------------------------------------------------------------
1745 // get information about nRows from all processors, and then
1746 // compute the base NRows on each processor
1747 // (This is needed later on for column index conversion)
1748 // ==> ProcNRows, globalNRows
1749 //------------------------------------------------------------------
1750
1751 ProcNRows = new int[numProcs_];
1752 tempList = new int[numProcs_];
1753 for ( i = 0; i < numProcs_; i++ ) tempList[i] = 0;
1754 tempList[mypid_] = nRows;
1755 MPI_Allreduce(tempList, ProcNRows, numProcs_, MPI_INT, MPI_SUM, comm_);
1756 delete [] tempList;
1757 globalNRows = 0;
1758 ncnt = 0;
1759 for ( i = 0; i < numProcs_; i++ )
1760 {
1761 globalNRows += ProcNRows[i];
1762 ncnt2 = ProcNRows[i];
1763 ProcNRows[i] = ncnt;
1764 ncnt += ncnt2;
1765 }
1766 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1767 printf("%4d : SlideReduction2 : localNRows = %d\n", mypid_, nRows);
1768
1769 //------------------------------------------------------------------
1770 // compose a global array marking where the constraint equations are,
1771 // then compute the base nConstraints on each processor
1772 // (This is needed later on for column index conversion)
1773 // ==> ProcNConstr, globalNConstr
1774 //------------------------------------------------------------------
1775
1776 globalNConstr = 0;
1777 tempList = new int[numProcs_];
1778 ProcNConstr = new int[numProcs_];
1779 for ( i = 0; i < numProcs_; i++ ) tempList[i] = 0;
1780 tempList[mypid_] = nConstraints_;
1781 MPI_Allreduce(tempList,ProcNConstr,numProcs_,MPI_INT,MPI_SUM,comm_);
1782 delete [] tempList;
1783 ncnt = 0;
1784 for ( i = 0; i < numProcs_; i++ )
1785 {
1786 globalNConstr += ProcNConstr[i];
1787 ncnt2 = ProcNConstr[i];
1788 ProcNConstr[i] = ncnt;
1789 ncnt += ncnt2;
1790 }
1791 #ifdef HAVE_MLI
1792 if ( HYPreconID_ == HYMLI )
1793 HYPRE_LSI_MLIAdjustNodeEqnMap(HYPrecon_, ProcNRows, ProcNConstr);
1794 #endif
1795
1796 //******************************************************************
1797 // compose the local and global selected node lists
1798 //------------------------------------------------------------------
1799
1800 if ( selectedList_ != NULL ) delete [] selectedList_;
1801 if ( selectedListAux_ != NULL ) delete [] selectedListAux_;
1802 nSelected = nConstraints_;
1803 if ( nConstraints_ > 0 )
1804 {
1805 selectedList_ = new int[nConstraints_];
1806 selectedListAux_ = new int[nConstraints_];
1807 }
1808 else selectedList_ = selectedListAux_ = NULL;
1809 globalNSelected = globalNConstr;
1810 if (globalNSelected > 0)
1811 {
1812 globalSelectedList = new int[globalNSelected];
1813 globalSelectedListAux = new int[globalNSelected];
1814 }
1815 else globalSelectedList = globalSelectedListAux = NULL;
1816
1817 buildSlideReducedSystemPartA(ProcNRows,ProcNConstr,globalNRows,
1818 globalNSelected,globalSelectedList,
1819 globalSelectedListAux);
1820
1821 //******************************************************************
1822 // construct A21
1823 //------------------------------------------------------------------
1824
1825 //------------------------------------------------------------------
1826 // calculate the dimension of A21
1827 //------------------------------------------------------------------
1828
1829 A21NRows = 2 * nConstraints_;
1830 A21NCols = nRows - nConstraints_;
1831 A21GlobalNRows = 2 * globalNConstr;
1832 A21GlobalNCols = globalNRows - globalNConstr;
1833 A21StartRow = 2 * ProcNConstr[mypid_];
1834 A21StartCol = ProcNRows[mypid_] - ProcNConstr[mypid_];
1835
1836 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
1837 {
1838 printf("%4d : SlideReduction2 : A21StartRow = %d\n", mypid_,
1839 A21StartRow);
1840 printf("%4d : SlideReduction2 : A21GlobalDim = %d %d\n", mypid_,
1841 A21GlobalNRows, A21GlobalNCols);
1842 printf("%4d : SlideReduction2 : A21LocalDim = %d %d\n",mypid_,
1843 A21NRows, A21NCols);
1844 }
1845
1846 //------------------------------------------------------------------
1847 // create a matrix context for A21
1848 //------------------------------------------------------------------
1849
1850 ierr = HYPRE_IJMatrixCreate(comm_, A21StartRow, A21StartRow+A21NRows-1,
1851 A21StartCol, A21StartCol+A21NCols-1, &A21);
1852 ierr += HYPRE_IJMatrixSetObjectType(A21, HYPRE_PARCSR);
1853 hypre_assert(!ierr);
1854
1855 //------------------------------------------------------------------
1856 // compute the number of nonzeros in the first nConstraint row of A21
1857 // (which consists of the rows in selectedList), the nnz will
1858 // be reduced by excluding the constraint and selected slave columns
1859 //------------------------------------------------------------------
1860
1861 rowCount = 0;
1862 maxRowSize = 0;
1863 newEndRow = EndRow - nConstraints_;
1864 A21MatSize = new int[A21NRows];
1865
1866 for ( i = 0; i < nSelected; i++ )
1867 {
1868 for ( j = 0; j < nSelected; j++ )
1869 {
1870 if ( selectedListAux_[j] == i )
1871 {
1872 rowIndex = selectedList_[j];
1873 break;
1874 }
1875 }
1876 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
1877
1878 newRowSize = 0;
1879 for (j = 0; j < rowSize; j++)
1880 {
1881 colIndex = colInd[j];
1882 if ( colVal[j] != 0.0 )
1883 {
1884 if (colIndex <= newEndRow || colIndex >= localEndRow_)
1885 {
1886 if (colIndex >= StartRow && colIndex <= newEndRow )
1887 searchIndex = hypre_BinarySearch(selectedList_,colIndex,
1888 nSelected);
1889 else
1890 searchIndex = hypre_BinarySearch(globalSelectedList,
1891 colIndex, globalNSelected);
1892 if (searchIndex < 0 ) newRowSize++;
1893 }
1894 }
1895 }
1896 A21MatSize[rowCount] = newRowSize;
1897 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
1898 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
1899 rowCount++;
1900 }
1901
1902 //------------------------------------------------------------------
1903 // compute the number of nonzeros in the second nConstraint row of A21
1904 // (which consists of the rows in constraint equations), the nnz will
1905 // be reduced by excluding the selected slave columns only (since the
1906 // entries corresponding to the constraint columns are 0, and since
1907 // the selected matrix is a diagonal matrix, there is no need to
1908 // search for slave equations in the off-processor list)
1909 //------------------------------------------------------------------
1910
1911 rowCount = nSelected;
1912 for ( i = EndRow-nConstraints_+1; i <= EndRow; i++ )
1913 {
1914 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
1915 newRowSize = 0;
1916 for (j = 0; j < rowSize; j++)
1917 {
1918 if ( colVal[j] != 0.0 )
1919 {
1920 colIndex = colInd[j];
1921 if (colIndex <= newEndRow || colIndex >= localEndRow_)
1922 {
1923 if (colIndex >= StartRow && colIndex <= newEndRow )
1924 searchIndex = hypre_BinarySearch(selectedList_,colIndex,
1925 nSelected);
1926 else
1927 searchIndex = hypre_BinarySearch(globalSelectedList,
1928 colIndex, globalNSelected);
1929 if ( searchIndex < 0 ) newRowSize++;
1930 }
1931 }
1932 }
1933 A21MatSize[rowCount] = newRowSize;
1934 maxRowSize = ( newRowSize > maxRowSize ) ? newRowSize : maxRowSize;
1935 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
1936 rowCount++;
1937 }
1938 nnzA21 = 0;
1939 for ( i = 0; i < 2*nConstraints_; i++ ) nnzA21 += A21MatSize[i];
1940
1941 //------------------------------------------------------------------
1942 // after fetching the row sizes, set up A21 with such sizes
1943 //------------------------------------------------------------------
1944
1945 ierr = HYPRE_IJMatrixSetRowSizes(A21, A21MatSize);
1946 ierr += HYPRE_IJMatrixInitialize(A21);
1947 hypre_assert(!ierr);
1948 delete [] A21MatSize;
1949
1950 //------------------------------------------------------------------
1951 // next load the first nConstraint row to A21 extracted from A
1952 // (at the same time, the D block is saved for future use)
1953 //------------------------------------------------------------------
1954
1955 rowCount = A21StartRow;
1956 if ( nConstraints_ > 0 ) diagonal = new double[nConstraints_];
1957 else diagonal = NULL;
1958 newColInd = new int[maxRowSize+1];
1959 newColVal = new double[maxRowSize+1];
1960
1961 diagCount = 0;
1962 for ( i = 0; i < nSelected; i++ )
1963 {
1964 for ( j = 0; j < nSelected; j++ )
1965 {
1966 if ( selectedListAux_[j] == i )
1967 {
1968 rowIndex = selectedList_[j];
1969 break;
1970 }
1971 }
1972 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
1973 newRowSize = 0;
1974 for (j = 0; j < rowSize; j++)
1975 {
1976 if ( colVal[j] != 0.0 )
1977 {
1978 colIndex = colInd[j];
1979
1980 if (colIndex <= newEndRow || colIndex >= localEndRow_)
1981 {
1982 searchIndex = hypre_BinarySearch(globalSelectedList,
1983 colIndex, globalNSelected);
1984 if ( searchIndex < 0 )
1985 {
1986 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
1987 if ( ProcNRows[procIndex] > colIndex ) break;
1988 procIndex--;
1989 newColIndex = colInd[j] - ProcNConstr[procIndex];
1990 newColInd[newRowSize] = newColIndex;
1991 newColVal[newRowSize++] = colVal[j];
1992 if ( newColIndex < 0 || newColIndex >= A21GlobalNCols )
1993 {
1994 if (HYOutputLevel_ & HYFEI_SLIDEREDUCE1)
1995 {
1996 printf("%4d : SlideReduction2 WARNING - ",mypid_);
1997 printf(" A21(%d,%d(%d))\n", rowCount,
1998 colIndex, A21GlobalNCols);
1999 }
2000 }
2001 if ( newRowSize > maxRowSize+1 )
2002 {
2003 if (HYOutputLevel_ & HYFEI_SLIDEREDUCE1)
2004 {
2005 printf("%4d : SlideReduction2 : WARNING - ",mypid_);
2006 printf("cross array boundary(1).\n");
2007 }
2008 }
2009 }
2010 }
2011
2012 //---------------------------------------------------------
2013 // slave equations should only have one nonzeros
2014 // corresponding to the D in A22
2015 //---------------------------------------------------------
2016
2017 else if ( colIndex > newEndRow && colIndex <= EndRow )
2018 {
2019 if ( colVal[j] != 0.0 ) diagonal[diagCount++] = colVal[j];
2020 if ( habs(colVal[j]) < 1.0E-8 )
2021 {
2022 if (HYOutputLevel_ & HYFEI_SLIDEREDUCE1)
2023 {
2024 printf("%4d : SlideReduction2 WARNING : large ",mypid_);
2025 printf("entry in invA22\n");
2026 }
2027 }
2028 }
2029 }
2030 }
2031
2032 HYPRE_IJMatrixSetValues(A21, 1, &newRowSize, (const int *) &rowCount,
2033 (const int *) newColInd, (const double *) newColVal);
2034 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2035 if ( diagCount != (i+1) )
2036 {
2037 printf("%4d : SlideReduction2 ERROR (3) : %d %d.\n", mypid_,
2038 diagCount,i+1);
2039 exit(1);
2040 }
2041 rowCount++;
2042 }
2043
2044 //------------------------------------------------------------------
2045 // send the diagonal to each processor that needs them
2046 //------------------------------------------------------------------
2047
2048 recvCntArray = new int[numProcs_];
2049 displArray = new int[numProcs_];
2050 MPI_Allgather(&diagCount, 1, MPI_INT, recvCntArray, 1, MPI_INT, comm_);
2051 displArray[0] = 0;
2052 for ( i = 1; i < numProcs_; i++ )
2053 displArray[i] = displArray[i-1] + recvCntArray[i-1];
2054 ncnt = displArray[numProcs_-1] + recvCntArray[numProcs_-1];
2055 if ( ncnt > 0 ) extDiagonal = new double[ncnt];
2056 else extDiagonal = NULL;
2057 MPI_Allgatherv(diagonal, diagCount, MPI_DOUBLE, extDiagonal,
2058 recvCntArray, displArray, MPI_DOUBLE, comm_);
2059 diagCount = ncnt;
2060 delete [] recvCntArray;
2061 delete [] displArray;
2062 if ( diagonal != NULL ) delete [] diagonal;
2063
2064 //------------------------------------------------------------------
2065 // next load the second nConstraint rows to A21 extracted from A
2066 // (assume the constraint-constraint block is 0 )
2067 //------------------------------------------------------------------
2068
2069 for ( i = EndRow-nConstraints_+1; i <= EndRow; i++ )
2070 {
2071 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
2072 newRowSize = 0;
2073 for (j = 0; j < rowSize; j++)
2074 {
2075 colIndex = colInd[j];
2076 if ( colVal[j] != 0.0 )
2077 {
2078 if (colIndex <= newEndRow || colIndex >= localEndRow_)
2079 {
2080 if (colIndex >= StartRow && colIndex <= newEndRow )
2081 searchIndex = hypre_BinarySearch(selectedList_,colIndex,
2082 nSelected);
2083 else
2084 searchIndex = hypre_BinarySearch(globalSelectedList,
2085 colIndex, globalNSelected);
2086
2087 if ( searchIndex < 0 )
2088 {
2089 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
2090 if ( ProcNRows[procIndex] > colIndex ) break;
2091 procIndex--;
2092 newColIndex = colInd[j] - ProcNConstr[procIndex];
2093 newColInd[newRowSize] = newColIndex;
2094 newColVal[newRowSize++] = colVal[j];
2095 if ( newColIndex < 0 || newColIndex >= A21GlobalNCols )
2096 {
2097 if (HYOutputLevel_ & HYFEI_SLIDEREDUCE1)
2098 {
2099 printf("%4d : SlideReduction2 WARNING - ",mypid_);
2100 printf(" A21(%d,%d(%d))\n", rowCount,
2101 colIndex, A21GlobalNCols);
2102 }
2103 }
2104 if ( newRowSize > maxRowSize+1 )
2105 {
2106 if (HYOutputLevel_ & HYFEI_SLIDEREDUCE1)
2107 {
2108 printf("%4d : SlideReductionWARNING : ",mypid_);
2109 printf("crossing array boundary(2).\n");
2110 }
2111 }
2112 }
2113 }
2114 }
2115 }
2116 if ( newRowSize == 0 && (HYOutputLevel_ & HYFEI_SLIDEREDUCE1))
2117 printf("%4d : SlideReduction2 WARNING : loading all 0 to A21\n",
2118 mypid_);
2119 HYPRE_IJMatrixSetValues(A21, 1, &newRowSize, (const int *) &rowCount,
2120 (const int *) newColInd, (const double *) newColVal);
2121 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
2122 rowCount++;
2123 }
2124 delete [] newColInd;
2125 delete [] newColVal;
2126
2127 //------------------------------------------------------------------
2128 // finally assemble the matrix and sanitize
2129 //------------------------------------------------------------------
2130
2131 HYPRE_IJMatrixAssemble(A21);
2132
2133 HYPRE_IJMatrixGetObject(A21, (void **) &A21_csr);
2134 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) A21_csr);
2135
2136 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
2137 {
2138 ncnt = 0;
2139 MPI_Barrier(comm_);
2140 while ( ncnt < numProcs_ )
2141 {
2142 if ( mypid_ == ncnt )
2143 {
2144 printf("====================================================\n");
2145 printf("%4d : SlideReduction2 : matrix A21 assembled %d.\n",
2146 mypid_,A21StartRow);
2147 fflush(stdout);
2148 for ( i = A21StartRow; i < A21StartRow+2*nConstraints_; i++ )
2149 {
2150 HYPRE_ParCSRMatrixGetRow(A21_csr,i,&rowSize,&colInd,&colVal);
2151 printf("A21 ROW = %6d (%d)\n", i, rowSize);
2152 for ( j = 0; j < rowSize; j++ )
2153 printf(" col = %6d, val = %e \n", colInd[j], colVal[j]);
2154 HYPRE_ParCSRMatrixRestoreRow(A21_csr,i,&rowSize,&colInd,
2155 &colVal);
2156 }
2157 printf("====================================================\n");
2158 }
2159 ncnt++;
2160 MPI_Barrier(comm_);
2161 }
2162 }
2163
2164 //******************************************************************
2165 // construct invA22
2166 //------------------------------------------------------------------
2167
2168 //------------------------------------------------------------------
2169 // calculate the dimension of invA22
2170 //------------------------------------------------------------------
2171
2172 invA22NRows = A21NRows;
2173 invA22NCols = invA22NRows;
2174 invA22GlobalNRows = A21GlobalNRows;
2175 invA22GlobalNCols = invA22GlobalNRows;
2176 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2177 {
2178 printf("%4d : SlideReduction2 - A22GlobalDim = %d %d\n", mypid_,
2179 invA22GlobalNRows, invA22GlobalNCols);
2180 printf("%4d : SlideReduction2 - A22LocalDim = %d %d\n", mypid_,
2181 invA22NRows, invA22NCols);
2182 }
2183
2184 //------------------------------------------------------------------
2185 // create a matrix context for A22
2186 //------------------------------------------------------------------
2187
2188 ierr = HYPRE_IJMatrixCreate(comm_, A21StartRow, A21StartRow+invA22NRows-1,
2189 A21StartRow, A21StartRow+invA22NCols-1, &invA22);
2190 ierr += HYPRE_IJMatrixSetObjectType(invA22, HYPRE_PARCSR);
2191 hypre_assert(!ierr);
2192
2193 //------------------------------------------------------------------
2194 // compute the no. of nonzeros in the first nConstraint row of invA22
2195 //------------------------------------------------------------------
2196
2197 maxRowSize = 0;
2198 invA22MatSize = new int[invA22NRows];
2199 for ( i = 0; i < nConstraints_; i++ ) invA22MatSize[i] = 1;
2200
2201 //------------------------------------------------------------------
2202 // compute the number of nonzeros in the second nConstraints row of
2203 // invA22 (consisting of [D and A22 block])
2204 //------------------------------------------------------------------
2205
2206 for ( i = 0; i < nSelected; i++ )
2207 {
2208 for ( j = 0; j < nSelected; j++ )
2209 {
2210 if ( selectedListAux_[j] == i )
2211 {
2212 rowIndex = selectedList_[j];
2213 break;
2214 }
2215 }
2216 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2217 rowSize2 = 1;
2218 for (j = 0; j < rowSize; j++)
2219 {
2220 colIndex = colInd[j];
2221 if ( colVal[j] != 0.0 )
2222 {
2223 if ( colIndex >= StartRow && colIndex <= newEndRow )
2224 {
2225 searchIndex = hypre_BinarySearch(selectedList_, colIndex,
2226 nSelected);
2227 if ( searchIndex >= 0 ) rowSize2++;
2228 }
2229 else if ( colIndex < StartRow || colIndex > EndRow )
2230 {
2231 searchIndex = hypre_BinarySearch(globalSelectedList,colIndex,
2232 globalNSelected);
2233 if ( searchIndex >= 0 ) rowSize2++;
2234 }
2235 }
2236 }
2237 invA22MatSize[nConstraints_+i] = rowSize2;
2238 maxRowSize = ( rowSize2 > maxRowSize ) ? rowSize2 : maxRowSize;
2239 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2240 }
2241
2242 //------------------------------------------------------------------
2243 // after fetching the row sizes, set up invA22 with such sizes
2244 //------------------------------------------------------------------
2245
2246 ierr = HYPRE_IJMatrixSetRowSizes(invA22, invA22MatSize);
2247 ierr += HYPRE_IJMatrixInitialize(invA22);
2248 hypre_assert(!ierr);
2249 delete [] invA22MatSize;
2250
2251 //------------------------------------------------------------------
2252 // next load the first nConstraints_ row to invA22 extracted from A
2253 // (that is, the D block)
2254 //------------------------------------------------------------------
2255
2256 maxRowSize++;
2257 newColInd = new int[maxRowSize];
2258 newColVal = new double[maxRowSize];
2259
2260 for ( i = 0; i < diagCount; i++ ) extDiagonal[i] = 1.0 / extDiagonal[i];
2261 for ( i = 0; i < nConstraints_; i++ )
2262 {
2263 newColInd[0] = A21StartRow + nConstraints_ + i;
2264 rowIndex = A21StartRow + i;
2265 if ( newColInd[0] < 0 || newColInd[0] >= invA22GlobalNCols )
2266 {
2267 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2268 printf("%4d : SlideReduction2 WARNING - A22(%d,%d(%d))\n",
2269 mypid_, rowIndex, newColInd[0], invA22GlobalNCols);
2270 }
2271 newColVal[0] = extDiagonal[A21StartRow/2+i];
2272 ierr = HYPRE_IJMatrixSetValues(invA22, 1, &one, (const int *) &rowIndex,
2273 (const int *) newColInd, (const double *) newColVal);
2274 hypre_assert(!ierr);
2275 }
2276
2277 //------------------------------------------------------------------
2278 // next load the second nConstraints_ rows to A22 extracted from A
2279 //------------------------------------------------------------------
2280
2281 for ( i = 0; i < nSelected; i++ )
2282 {
2283 for ( j = 0; j < nSelected; j++ )
2284 {
2285 if ( selectedListAux_[j] == i )
2286 {
2287 rowIndex = selectedList_[j];
2288 break;
2289 }
2290 }
2291 HYPRE_ParCSRMatrixGetRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2292 newRowSize = 1;
2293 newColInd[0] = A21StartRow + i;
2294 newColVal[0] = extDiagonal[A21StartRow/2+i];
2295 for (j = 0; j < rowSize; j++)
2296 {
2297 colIndex = colInd[j];
2298 if ( colVal[j] != 0.0 )
2299 {
2300 searchIndex = hypre_BinarySearch(globalSelectedList,
2301 colIndex,globalNSelected);
2302 if ( searchIndex >= 0 )
2303 {
2304 searchIndex = globalSelectedListAux[searchIndex];
2305 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
2306 if ( ProcNRows[procIndex] > colIndex ) break;
2307 if ( procIndex == numProcs_ )
2308 newColInd[newRowSize] = searchIndex + globalNConstr;
2309 else
2310 newColInd[newRowSize] = searchIndex+ProcNConstr[procIndex];
2311 if ( newColInd[newRowSize] < 0 ||
2312 newColInd[newRowSize] >= invA22GlobalNCols )
2313 {
2314 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2315 printf("%4d : SlideReduction2 WARNING - A22(%d,%d,%d)\n",
2316 mypid_, rowCount, newColInd[newRowSize],
2317 invA22GlobalNCols);
2318 }
2319 newColVal[newRowSize++] = - extDiagonal[A21StartRow/2+i] *
2320 colVal[j] * extDiagonal[searchIndex];
2321 if ( newRowSize > maxRowSize )
2322 {
2323 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2324 {
2325 printf("%4d : SlideReduction2 WARNING - ",mypid_);
2326 printf("passing array boundary(3).\n");
2327 }
2328 }
2329 }
2330 }
2331 }
2332 rowCount = A21StartRow + nConstraints_ + i;
2333 ierr = HYPRE_IJMatrixSetValues(invA22, 1, &newRowSize,
2334 (const int *) &rowCount, (const int *) newColInd,
2335 (const double *) newColVal);
2336 hypre_assert(!ierr);
2337 HYPRE_ParCSRMatrixRestoreRow(A_csr,rowIndex,&rowSize,&colInd,&colVal);
2338 }
2339 delete [] newColInd;
2340 delete [] newColVal;
2341 delete [] extDiagonal;
2342
2343 //------------------------------------------------------------------
2344 // finally assemble the matrix and sanitize
2345 //------------------------------------------------------------------
2346
2347 HYPRE_IJMatrixAssemble(invA22);
2348 HYPRE_IJMatrixGetObject(invA22, (void **) &invA22_csr);
2349 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) invA22_csr);
2350
2351 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
2352 {
2353 ncnt = 0;
2354 MPI_Barrier(comm_);
2355 while ( ncnt < numProcs_ )
2356 {
2357 if ( mypid_ == ncnt )
2358 {
2359 printf("====================================================\n");
2360 printf("%4d : SlideReduction - invA22 \n", mypid_);
2361 for ( i = A21StartRow; i < A21StartRow+2*nConstraints_; i++ )
2362 {
2363 HYPRE_ParCSRMatrixGetRow(invA22_csr,i,&rowSize,&colInd,
2364 &colVal);
2365 printf("invA22 ROW = %6d (%d)\n", i, rowSize);
2366 for ( j = 0; j < rowSize; j++ )
2367 printf(" col = %6d, val = %e \n", colInd[j], colVal[j]);
2368 HYPRE_ParCSRMatrixRestoreRow(invA22_csr,i,&rowSize,&colInd,
2369 &colVal);
2370 }
2371 printf("====================================================\n");
2372 }
2373 MPI_Barrier(comm_);
2374 ncnt++;
2375 }
2376 }
2377
2378 //******************************************************************
2379 // perform the triple matrix product
2380 //------------------------------------------------------------------
2381
2382 HYPRE_IJMatrixGetObject(A21, (void **) &A21_csr);
2383 HYPRE_IJMatrixGetObject(invA22, (void **) &invA22_csr);
2384 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2385 printf("%4d : SlideReduction2 - Triple matrix product starts\n",mypid_);
2386
2387 hypre_BoomerAMGBuildCoarseOperator( (hypre_ParCSRMatrix *) A21_csr,
2388 (hypre_ParCSRMatrix *) invA22_csr,
2389 (hypre_ParCSRMatrix *) A21_csr,
2390 (hypre_ParCSRMatrix **) &RAP_csr);
2391
2392 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2393 printf("%4d : SlideReduction2 - Triple matrix product ends\n", mypid_);
2394
2395 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
2396 {
2397 MPI_Barrier(comm_);
2398 ncnt = 0;
2399 while ( ncnt < numProcs_ )
2400 {
2401 if ( mypid_ == ncnt )
2402 {
2403 for ( i = A21StartRow; i < A21StartRow+A21NCols; i++ )
2404 {
2405 HYPRE_ParCSRMatrixGetRow(RAP_csr,i,&rowSize,&colInd, &colVal);
2406 printf("RAP ROW = %6d (%d)\n", i, rowSize);
2407 for ( j = 0; j < rowSize; j++ )
2408 printf(" col = %6d, val = %e \n", colInd[j], colVal[j]);
2409 HYPRE_ParCSRMatrixRestoreRow(RAP_csr,i,&rowSize,&colInd,
2410 &colVal);
2411 }
2412 }
2413 MPI_Barrier(comm_);
2414 ncnt++;
2415 }
2416 }
2417
2418 //******************************************************************
2419 // finally formed the Schur complement reduced system by
2420 // extracting the A11 part of A and subtracting the RAP
2421 //------------------------------------------------------------------
2422
2423 //------------------------------------------------------------------
2424 // first calculate the dimension of the reduced matrix
2425 //------------------------------------------------------------------
2426
2427 newNRows = nRows - nConstraints_;
2428 ierr = HYPRE_IJMatrixCreate(comm_, A21StartCol, A21StartCol+newNRows-1,
2429 A21StartCol, A21StartCol+newNRows-1, &reducedA);
2430 ierr += HYPRE_IJMatrixSetObjectType(reducedA, HYPRE_PARCSR);
2431 hypre_assert(!ierr);
2432
2433 //------------------------------------------------------------------
2434 // set up reducedA with proper sizes
2435 //------------------------------------------------------------------
2436
2437 reducedAMatSize = new int[newNRows];
2438 reducedAStartRow = ProcNRows[mypid_] - ProcNConstr[mypid_];
2439 rowCount = reducedAStartRow;
2440 rowIndex = 0;
2441
2442 for ( i = StartRow; i <= newEndRow; i++ )
2443 {
2444 searchIndex = hypre_BinarySearch(selectedList_, i, nSelected);
2445 if ( searchIndex < 0 )
2446 {
2447 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
2448 ierr = HYPRE_ParCSRMatrixGetRow(RAP_csr,rowCount,&rowSize2,
2449 &colInd2, &colVal2);
2450 hypre_assert( !ierr );
2451 newRowSize = rowSize + rowSize2;
2452 newColInd = new int[newRowSize];
2453 for (j = 0; j < rowSize; j++) newColInd[j] = colInd[j];
2454 for (j = 0; j < rowSize2; j++) newColInd[rowSize+j] = colInd2[j];
2455 hypre_qsort0(newColInd, 0, newRowSize-1);
2456 ncnt = 0;
2457 for ( j = 0; j < newRowSize; j++ )
2458 {
2459 if ( newColInd[j] != newColInd[ncnt] )
2460 {
2461 ncnt++;
2462 newColInd[ncnt] = newColInd[j];
2463 }
2464 }
2465 reducedAMatSize[rowIndex++] = ncnt;
2466
2467 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
2468 ierr = HYPRE_ParCSRMatrixRestoreRow(RAP_csr,rowCount,&rowSize2,
2469 &colInd2,&colVal2);
2470 delete [] newColInd;
2471 hypre_assert( !ierr );
2472 rowCount++;
2473 }
2474 else
2475 {
2476 reducedAMatSize[rowIndex++] = 1;
2477 rowCount++;
2478 }
2479 }
2480
2481 //------------------------------------------------------------------
2482 // create a matrix context for reducedA
2483 //------------------------------------------------------------------
2484
2485 ierr = HYPRE_IJMatrixSetRowSizes(reducedA, reducedAMatSize);
2486 ierr += HYPRE_IJMatrixInitialize(reducedA);
2487 hypre_assert(!ierr);
2488 delete [] reducedAMatSize;
2489
2490 //------------------------------------------------------------------
2491 // load the reducedA matrix
2492 //------------------------------------------------------------------
2493
2494 rowCount = reducedAStartRow;
2495 for ( i = StartRow; i <= newEndRow; i++ )
2496 {
2497 searchIndex = hypre_BinarySearch(selectedList_, i, nSelected);
2498 if ( searchIndex < 0 )
2499 {
2500 HYPRE_ParCSRMatrixGetRow(A_csr, i, &rowSize, &colInd, &colVal);
2501 HYPRE_ParCSRMatrixGetRow(RAP_csr,rowCount,&rowSize2,&colInd2,
2502 &colVal2);
2503 newRowSize = rowSize + rowSize2;
2504 newColInd = new int[newRowSize];
2505 newColVal = new double[newRowSize];
2506 ncnt = 0;
2507
2508 for ( j = 0; j < rowSize; j++ )
2509 {
2510 colIndex = colInd[j];
2511 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
2512 if ( ProcNRows[procIndex] > colIndex ) break;
2513 if ( procIndex == numProcs_ )
2514 ubound = globalNRows-(globalNConstr-ProcNConstr[numProcs_-1]);
2515 else
2516 ubound = ProcNRows[procIndex] -
2517 (ProcNConstr[procIndex]-ProcNConstr[procIndex-1]);
2518 procIndex--;
2519 if ( colIndex < ubound )
2520 {
2521 if ( colIndex >= StartRow && colIndex <= EndRow )
2522 searchIndex = HYPRE_LSI_Search(selectedList_,colIndex,
2523 nSelected);
2524 else
2525 searchIndex = HYPRE_LSI_Search(globalSelectedList,colIndex,
2526 globalNSelected);
2527
2528 if ( searchIndex < 0 )
2529 {
2530 newColInd[ncnt] = colIndex - ProcNConstr[procIndex];
2531 newColVal[ncnt++] = colVal[j];
2532 }
2533 }
2534 }
2535 for ( j = 0; j < rowSize2; j++ )
2536 {
2537 newColInd[ncnt+j] = colInd2[j];
2538 newColVal[ncnt+j] = - colVal2[j];
2539 }
2540 newRowSize = ncnt + rowSize2;
2541 hypre_qsort1(newColInd, newColVal, 0, newRowSize-1);
2542 ncnt = 0;
2543 for ( j = 0; j < newRowSize; j++ )
2544 {
2545 if ( j != ncnt && newColInd[j] == newColInd[ncnt] )
2546 newColVal[ncnt] += newColVal[j];
2547 else if ( newColInd[j] != newColInd[ncnt] )
2548 {
2549 ncnt++;
2550 newColVal[ncnt] = newColVal[j];
2551 newColInd[ncnt] = newColInd[j];
2552 }
2553 }
2554 newRowSize = ncnt + 1;
2555 ierr = HYPRE_IJMatrixSetValues(reducedA, 1, &newRowSize,
2556 (const int *) &rowCount, (const int *) newColInd,
2557 (const double *) newColVal);
2558 hypre_assert(!ierr);
2559 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
2560 HYPRE_ParCSRMatrixRestoreRow(RAP_csr,rowCount,&rowSize2,&colInd2,
2561 &colVal2);
2562 rowCount++;
2563 delete [] newColInd;
2564 delete [] newColVal;
2565 }
2566 else
2567 {
2568 newRowSize = 1;
2569 newColInd = new int[newRowSize];
2570 newColVal = new double[newRowSize];
2571 newColInd[0] = rowCount;
2572 newColVal[0] = 1.0;
2573 ierr = HYPRE_IJMatrixSetValues(reducedA, 1, &newRowSize,
2574 (const int *) &rowCount, (const int *) newColInd,
2575 (const double *) newColVal);
2576 hypre_assert(!ierr);
2577 rowCount++;
2578 delete [] newColInd;
2579 delete [] newColVal;
2580 }
2581 }
2582
2583 //------------------------------------------------------------------
2584 // assemble the reduced matrix
2585 //------------------------------------------------------------------
2586
2587 HYPRE_IJMatrixAssemble(reducedA);
2588 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2589 printf("%4d : SlideReduction2 - reducedA - StartRow = %d\n",
2590 mypid_, reducedAStartRow);
2591
2592 HYPRE_IJMatrixGetObject(reducedA, (void **) &reducedA_csr);
2593
2594 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
2595 {
2596 MPI_Barrier(comm_);
2597 ncnt = 0;
2598 while ( ncnt < numProcs_ )
2599 {
2600 if ( mypid_ == ncnt )
2601 {
2602 printf("====================================================\n");
2603 for ( i = reducedAStartRow;
2604 i < reducedAStartRow+nRows-2*nConstraints_; i++ )
2605 {
2606 printf("%d : reducedA ROW %d\n", mypid_, i);
2607 ierr = HYPRE_ParCSRMatrixGetRow(reducedA_csr,i,&rowSize,
2608 &colInd,&colVal);
2609 //hypre_qsort1(colInd, colVal, 0, rowSize-1);
2610 for ( j = 0; j < rowSize; j++ )
2611 if ( colVal[j] != 0.0 )
2612 printf("%4d %4d %20.13e\n", i+1, colInd[j]+1, colVal[j]);
2613 HYPRE_ParCSRMatrixRestoreRow(reducedA_csr,i,&rowSize,&colInd,
2614 &colVal);
2615 }
2616 printf("====================================================\n");
2617 }
2618 MPI_Barrier(comm_);
2619 ncnt++;
2620 }
2621 }
2622
2623 // *****************************************************************
2624 // form modified right hand side (f1 = f1 - A12*invA22*f2)
2625 // *****************************************************************
2626
2627 // *****************************************************************
2628 // form f2hat = invA22 * f2
2629 //------------------------------------------------------------------
2630
2631 HYPRE_IJVectorCreate(comm_, A21StartRow, A21StartRow+A21NRows-1, &f2);
2632 HYPRE_IJVectorSetObjectType(f2, HYPRE_PARCSR);
2633 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2634 printf("%4d : SlideReduction2 - A21 dims = %d %d %d\n", mypid_,
2635 A21StartRow, A21NRows, A21GlobalNRows);
2636 ierr += HYPRE_IJVectorInitialize(f2);
2637 ierr += HYPRE_IJVectorAssemble(f2);
2638 hypre_assert(!ierr);
2639
2640 HYPRE_IJVectorCreate(comm_, A21StartRow, A21StartRow+A21NRows-1, &f2hat);
2641 HYPRE_IJVectorSetObjectType(f2hat, HYPRE_PARCSR);
2642 ierr += HYPRE_IJVectorInitialize(f2hat);
2643 ierr += HYPRE_IJVectorAssemble(f2hat);
2644 hypre_assert(!ierr);
2645
2646 colInd = new int[nSelected*2];
2647 colVal = new double[nSelected*2];
2648
2649 for ( i = 0; i < nSelected; i++ )
2650 {
2651 for ( j = 0; j < nSelected; j++ )
2652 {
2653 if ( selectedListAux_[j] == i )
2654 {
2655 colInd[i] = selectedList_[j];
2656 break;
2657 }
2658 }
2659 if ( colInd[i] < 0 )
2660 {
2661 printf("%4d : SlideReduction2 ERROR - out of range %d\n", mypid_,
2662 colInd[i]);
2663 exit(1);
2664 }
2665 }
2666 for ( i = 0; i < nSelected; i++ )
2667 {
2668 colInd[nSelected+i] = EndRow - nConstraints_ + i + 1;
2669 }
2670 HYPRE_IJVectorGetValues(HYb_, 2*nSelected, colInd, colVal);
2671 for ( i = 0; i < nSelected*2; i++ ) colInd[i] = A21StartRow + i;
2672 ierr = HYPRE_IJVectorSetValues(f2, 2*nSelected, (const int *) colInd,
2673 (const double *) colVal);
2674 hypre_assert( !ierr );
2675 HYPRE_IJVectorGetObject(f2, (void **) &f2_csr);
2676 HYPRE_IJVectorGetObject(f2hat, (void **) &f2hat_csr);
2677 HYPRE_ParCSRMatrixMatvec( 1.0, invA22_csr, f2_csr, 0.0, f2hat_csr );
2678 delete [] colVal;
2679 delete [] colInd;
2680 HYPRE_IJVectorDestroy(f2);
2681
2682 // *****************************************************************
2683 // set up A12 with proper sizes before forming f2til = A12 * f2hat
2684 //------------------------------------------------------------------
2685
2686 //------------------------------------------------------------------
2687 // calculate the dimension of A12
2688 //------------------------------------------------------------------
2689
2690 A12NRows = A21NCols;
2691 A12NCols = A21NRows;
2692 A12GlobalNRows = A21GlobalNCols;
2693 A12GlobalNCols = A21GlobalNRows;
2694 A12MatSize = new int[A12NRows];
2695 A12StartRow = ProcNRows[mypid_] - ProcNConstr[mypid_];
2696 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2697 {
2698 printf("%4d : SlideReduction2 - A12GlobalDim = %d %d\n", mypid_,
2699 A12GlobalNRows, A12GlobalNCols);
2700 printf("%4d : SlideReduction2 - A12LocalDim = %d %d\n", mypid_,
2701 A12NRows, A12NCols);
2702 }
2703
2704 //------------------------------------------------------------------
2705 // create a matrix context for A12
2706 //------------------------------------------------------------------
2707
2708 ierr = HYPRE_IJMatrixCreate(comm_, A21StartCol, A21StartCol+A12NRows-1,
2709 A21StartRow, A21StartRow+A12NCols-1, &A12);
2710 ierr += HYPRE_IJMatrixSetObjectType(A12, HYPRE_PARCSR);
2711 hypre_assert(!ierr);
2712
2713 //------------------------------------------------------------------
2714 // compute the number of nonzeros in each row of A12
2715 // (which consists of the rows in selectedList and the constraints)
2716 //------------------------------------------------------------------
2717
2718 rowCount = A12StartRow;
2719 rowIndex = 0;
2720 nnzA12 = 0;
2721
2722 for ( i = StartRow; i <= newEndRow; i++ )
2723 {
2724 searchIndex = hypre_BinarySearch(selectedList_, i, nSelected);
2725 if ( searchIndex < 0 )
2726 {
2727 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
2728 newRowSize = 0;
2729 for (j = 0; j < rowSize; j++)
2730 {
2731 if ( colVal[j] != 0.0 )
2732 {
2733 colIndex = colInd[j];
2734 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
2735 if ( ProcNRows[procIndex] > colIndex ) break;
2736 if ( procIndex == numProcs_ )
2737 ubound = globalNRows -
2738 (globalNConstr-ProcNConstr[numProcs_-1]);
2739 else
2740 ubound = ProcNRows[procIndex] -
2741 (ProcNConstr[procIndex]-ProcNConstr[procIndex-1]);
2742 procIndex--;
2743 if ( colIndex >= ubound ) newRowSize++;
2744 else if (colIndex >= StartRow && colIndex <= EndRow)
2745 {
2746 if (hypre_BinarySearch(selectedList_,colIndex,nSelected)>=0)
2747 newRowSize++;
2748 }
2749 else
2750 {
2751 if (hypre_BinarySearch(globalSelectedList,colIndex,
2752 globalNSelected) >= 0)
2753 newRowSize++;
2754 }
2755 }
2756 }
2757 A12MatSize[rowIndex++] = newRowSize;
2758 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
2759 rowCount++;
2760 }
2761 else
2762 {
2763 A12MatSize[rowIndex++] = 1;
2764 rowCount++;
2765 nnzA12--;
2766 }
2767 }
2768
2769 //------------------------------------------------------------------
2770 // after fetching the row sizes, set up A12 with such sizes
2771 //------------------------------------------------------------------
2772
2773 for ( i = 0; i < A12NRows; i++ ) nnzA12 += A12MatSize[i];
2774 ierr = HYPRE_IJMatrixSetRowSizes(A12, A12MatSize);
2775 ierr += HYPRE_IJMatrixInitialize(A12);
2776 hypre_assert(!ierr);
2777 delete [] A12MatSize;
2778
2779 //------------------------------------------------------------------
2780 // load the A12 matrix
2781 //------------------------------------------------------------------
2782
2783 rowCount = A12StartRow;
2784 for ( i = StartRow; i <= newEndRow; i++ )
2785 {
2786 searchIndex = hypre_BinarySearch(selectedList_, i, nSelected);
2787 if ( searchIndex < 0 )
2788 {
2789 HYPRE_ParCSRMatrixGetRow(A_csr, i, &rowSize, &colInd, &colVal);
2790 newRowSize = 0;
2791 newColInd = new int[rowSize];
2792 newColVal = new double[rowSize];
2793 for (j = 0; j < rowSize; j++)
2794 {
2795 colIndex = colInd[j];
2796 if ( colVal[j] != 0.0 )
2797 {
2798 for ( procIndex = 0; procIndex < numProcs_; procIndex++ )
2799 if ( ProcNRows[procIndex] > colIndex ) break;
2800 if ( procIndex == numProcs_ )
2801 ubound = globalNRows -
2802 (globalNConstr - ProcNConstr[numProcs_-1]);
2803 else
2804 ubound = ProcNRows[procIndex] -
2805 (ProcNConstr[procIndex]-ProcNConstr[procIndex-1]);
2806 procIndex--;
2807 if ( colIndex >= ubound ) {
2808 if ( procIndex != numProcs_ - 1 )
2809 {
2810 newColInd[newRowSize] = colInd[j] - ubound +
2811 ProcNConstr[procIndex] +
2812 ProcNConstr[procIndex+1];
2813 }
2814 else
2815 {
2816 newColInd[newRowSize] = colInd[j] - ubound +
2817 ProcNConstr[procIndex] +
2818 globalNConstr;
2819 }
2820 if ( newColInd[newRowSize] < 0 ||
2821 newColInd[newRowSize] >= A12GlobalNCols )
2822 {
2823 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2824 {
2825 printf("%4d : SlideReduction WARNING - A12 col index",
2826 mypid_);
2827 printf(" out of range %d %d(%d)\n", i,
2828 newColInd[newRowSize], A12GlobalNCols);
2829 }
2830 }
2831 newColVal[newRowSize++] = colVal[j];
2832 }
2833 else
2834 {
2835 if ( colInd[j] >= StartRow && colInd[j] <= EndRow )
2836 {
2837 searchIndex = HYPRE_LSI_Search(selectedList_,colInd[j],
2838 nSelected);
2839 if ( searchIndex >= 0 )
2840 searchIndex = selectedListAux_[searchIndex] +
2841 ProcNConstr[mypid_];
2842 } else {
2843 searchIndex = HYPRE_LSI_Search(globalSelectedList,
2844 colInd[j], globalNSelected);
2845 if ( searchIndex >= 0 )
2846 searchIndex = globalSelectedListAux[searchIndex];
2847 }
2848 if ( searchIndex >= 0)
2849 {
2850 newColInd[newRowSize] = searchIndex +
2851 ProcNConstr[procIndex];
2852 if ( newColInd[newRowSize] < 0 ||
2853 newColInd[newRowSize] >= A12GlobalNCols )
2854 {
2855 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2856 {
2857 printf("%4d : SlideReduction WARNING - A12 col ",
2858 mypid_);
2859 printf(" index out of range %d %d(%d)\n",i,
2860 newColInd[newRowSize], A12GlobalNCols);
2861 }
2862 }
2863 newColVal[newRowSize++] = colVal[j];
2864 }
2865 }
2866 }
2867 }
2868 ierr = HYPRE_IJMatrixSetValues(A12, 1, &newRowSize,
2869 (const int *) &rowCount, (const int *) newColInd,
2870 (const double *) newColVal);
2871 hypre_assert(!ierr);
2872 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
2873 rowCount++;
2874 delete [] newColInd;
2875 delete [] newColVal;
2876 }
2877 else
2878 {
2879 newRowSize = 1;
2880 newColInd = new int[newRowSize];
2881 newColVal = new double[newRowSize];
2882 newColInd[0] = A21StartRow;
2883 newColVal[0] = 0.0;
2884 ierr = HYPRE_IJMatrixSetValues(A12, 1, &newRowSize,
2885 (const int *) &rowCount, (const int *) newColInd,
2886 (const double *) newColVal);
2887 hypre_assert(!ierr);
2888 rowCount++;
2889 delete [] newColInd;
2890 delete [] newColVal;
2891 }
2892 }
2893
2894 //------------------------------------------------------------------
2895 // assemble the A12 matrix
2896 //------------------------------------------------------------------
2897
2898 ierr = HYPRE_IJMatrixAssemble(A12);
2899 hypre_assert( !ierr );
2900 HYPRE_IJMatrixGetObject(A12, (void **) &A12_csr);
2901
2902 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE3 )
2903 {
2904 MPI_Barrier(comm_);
2905 ncnt = 0;
2906 while ( ncnt < numProcs_ )
2907 {
2908 if ( mypid_ == ncnt )
2909 {
2910 printf("====================================================\n");
2911 for (i=A12StartRow;i<A12StartRow+A12NRows;i++)
2912 {
2913 printf("%d : A12 ROW %d\n", mypid_, i+1);
2914 HYPRE_ParCSRMatrixGetRow(A12_csr,i,&rowSize,&colInd,&colVal);
2915 //hypre_qsort1(colInd, colVal, 0, rowSize-1);
2916 for ( j = 0; j < rowSize; j++ )
2917 if ( colVal[j] != 0.0 )
2918 printf(" A12 %d %d %20.13e\n",i+1,colInd[j]+1,colVal[j]);
2919 HYPRE_ParCSRMatrixRestoreRow(A12_csr,i,&rowSize,&colInd,
2920 &colVal);
2921 }
2922 printf("====================================================\n");
2923 }
2924 MPI_Barrier(comm_);
2925 ncnt++;
2926 }
2927 }
2928
2929 //------------------------------------------------------------------
2930 // form reducedB_ = A12 * f2hat
2931 //------------------------------------------------------------------
2932
2933 if ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 )
2934 printf("%4d : SlideReduction2 - form reduced right hand side\n",mypid_);
2935 ierr = HYPRE_IJVectorCreate(comm_, reducedAStartRow,
2936 reducedAStartRow+newNRows-1, &reducedB_);
2937 ierr += HYPRE_IJVectorSetObjectType(reducedB_, HYPRE_PARCSR);
2938 ierr += HYPRE_IJVectorInitialize(reducedB_);
2939 ierr += HYPRE_IJVectorAssemble(reducedB_);
2940 hypre_assert( !ierr );
2941 HYPRE_IJVectorGetObject(reducedB_, (void **) &reducedB_csr);
2942
2943 HYPRE_ParCSRMatrixMatvec( -1.0, A12_csr, f2hat_csr, 0.0, reducedB_csr );
2944 HYPRE_IJMatrixDestroy(A12);
2945 HYPRE_IJVectorDestroy(f2hat);
2946
2947 //------------------------------------------------------------------
2948 // finally form reducedB = f1 - f2til
2949 //------------------------------------------------------------------
2950
2951 rowCount = reducedAStartRow;
2952 for ( i = StartRow; i <= newEndRow; i++ )
2953 {
2954 if ( hypre_BinarySearch(selectedList_, i, nSelected) < 0 )
2955 {
2956 HYPRE_IJVectorGetValues(HYb_, 1, &i, &ddata);
2957 HYPRE_IJVectorAddToValues(reducedB_, 1, (const int *) &rowCount,
2958 (const double *) &ddata);
2959 rowCount++;
2960 }
2961 else
2962 {
2963 ddata = 0.0;
2964 HYPRE_IJVectorSetValues(reducedB_, 1, (const int *) &rowCount,
2965 (const double *) &ddata);
2966 rowCount++;
2967 }
2968 }
2969
2970 //******************************************************************
2971 // set up the system with the new matrix
2972 //------------------------------------------------------------------
2973
2974 reducedA_ = reducedA;
2975 ierr = HYPRE_IJVectorCreate(comm_, reducedAStartRow,
2976 reducedAStartRow+newNRows-1, &reducedX_);
2977 ierr = HYPRE_IJVectorSetObjectType(reducedX_, HYPRE_PARCSR);
2978 ierr = HYPRE_IJVectorInitialize(reducedX_);
2979 ierr = HYPRE_IJVectorAssemble(reducedX_);
2980 hypre_assert(!ierr);
2981
2982 ierr = HYPRE_IJVectorCreate(comm_, reducedAStartRow,
2983 reducedAStartRow+newNRows-1, &reducedR_);
2984 ierr = HYPRE_IJVectorSetObjectType(reducedR_, HYPRE_PARCSR);
2985 ierr = HYPRE_IJVectorInitialize(reducedR_);
2986 ierr = HYPRE_IJVectorAssemble(reducedR_);
2987 hypre_assert(!ierr);
2988
2989 currA_ = reducedA_;
2990 currB_ = reducedB_;
2991 currR_ = reducedR_;
2992 currX_ = reducedX_;
2993
2994 //******************************************************************
2995 // save A21 and invA22 for solution recovery
2996 //------------------------------------------------------------------
2997
2998 HYA21_ = A21;
2999 HYinvA22_ = invA22;
3000
3001 //------------------------------------------------------------------
3002 // final clean up
3003 //------------------------------------------------------------------
3004
3005 if (globalSelectedList != NULL) delete [] globalSelectedList;
3006 if (globalSelectedListAux != NULL) delete [] globalSelectedListAux;
3007 if (ProcNRows != NULL) delete [] ProcNRows;
3008 if (ProcNConstr != NULL) delete [] ProcNConstr;
3009
3010 HYPRE_ParCSRMatrixDestroy(RAP_csr);
3011 if ( colIndices_ != NULL )
3012 {
3013 for ( i = 0; i < localEndRow_-localStartRow_+1; i++ )
3014 if ( colIndices_[i] != NULL ) delete [] colIndices_[i];
3015 delete [] colIndices_;
3016 colIndices_ = NULL;
3017 }
3018 if ( colValues_ != NULL )
3019 {
3020 for ( j = 0; j < localEndRow_-localStartRow_+1; j++ )
3021 if ( colValues_[j] != NULL ) delete [] colValues_[j];
3022 delete [] colValues_;
3023 colValues_ = NULL;
3024 if ( rowLengths_ != NULL )
3025 {
3026 delete [] rowLengths_;
3027 rowLengths_ = NULL;
3028 }
3029 }
3030
3031 //------------------------------------------------------------------
3032 // checking
3033 //------------------------------------------------------------------
3034
3035 MPI_Allreduce(&nnzA12,&ncnt,1,MPI_INT,MPI_SUM,comm_);
3036 if ( mypid_ == 0 && ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1) )
3037 printf(" SlideReduction2 - NNZ of A12 = %d\n", ncnt);
3038 MPI_Allreduce(&nnzA21,&ncnt,1,MPI_INT,MPI_SUM,comm_);
3039 if ( mypid_ == 0 && ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1) )
3040 printf(" SlideReduction2 - NNZ of A21 = %d\n", ncnt);
3041 }
3042
3043 //*****************************************************************************
3044 // Given the matrix (A) within the object, compute the reduced system and put
3045 // it in place. Additional information given are :
3046 //-----------------------------------------------------------------------------
3047
buildSlideReducedSoln()3048 double HYPRE_LinSysCore::buildSlideReducedSoln()
3049 {
3050 int i, j, *int_array, *gint_array, x2NRows, x2GlobalNRows;
3051 int ierr, rowNum, startRow, startRow2, index, localNRows;
3052 double ddata, rnorm;
3053 HYPRE_ParCSRMatrix A_csr, A21_csr, A22_csr;
3054 HYPRE_ParVector x_csr, x2_csr, r_csr, b_csr;
3055 HYPRE_IJVector R1, x2;
3056
3057 if ( HYA21_ == NULL || HYinvA22_ == NULL )
3058 {
3059 printf("buildSlideReducedSoln WARNING : A21 or A22 absent.\n");
3060 return (0.0);
3061 }
3062 else
3063 {
3064 //------------------------------------------------------------------
3065 // compute A21 * sol
3066 //------------------------------------------------------------------
3067
3068 x2NRows = 2 * nConstraints_;
3069 int_array = new int[numProcs_];
3070 gint_array = new int[numProcs_];
3071 for ( i = 0; i < numProcs_; i++ ) int_array[i] = 0;
3072 int_array[mypid_] = x2NRows;
3073 MPI_Allreduce(int_array,gint_array,numProcs_,MPI_INT,MPI_SUM,comm_);
3074 x2GlobalNRows = 0;
3075 for ( i = 0; i < numProcs_; i++ ) x2GlobalNRows += gint_array[i];
3076 rowNum = 0;
3077 for ( i = 0; i < mypid_; i++ ) rowNum += gint_array[i];
3078 startRow = rowNum;
3079 startRow2 = localStartRow_ - 1 - rowNum;
3080 delete [] int_array;
3081 delete [] gint_array;
3082 ierr = HYPRE_IJVectorCreate(comm_, startRow, startRow+x2NRows-1, &R1);
3083 ierr = HYPRE_IJVectorSetObjectType(R1, HYPRE_PARCSR);
3084 ierr = HYPRE_IJVectorInitialize(R1);
3085 ierr = HYPRE_IJVectorAssemble(R1);
3086 hypre_assert(!ierr);
3087 HYPRE_IJMatrixGetObject(HYA21_, (void **) &A21_csr);
3088 HYPRE_IJVectorGetObject(currX_, (void **) &x_csr);
3089 HYPRE_IJVectorGetObject(R1, (void **) &r_csr);
3090 HYPRE_ParCSRMatrixMatvec( -1.0, A21_csr, x_csr, 0.0, r_csr );
3091
3092 //------------------------------------------------------------------
3093 // f2 - A21 * sol
3094 //------------------------------------------------------------------
3095
3096 for ( i = 0; i < nConstraints_; i++ )
3097 {
3098 for ( j = 0; j < nConstraints_; j++ )
3099 {
3100 if ( selectedListAux_[j] == i )
3101 {
3102 index = selectedList_[j];
3103 break;
3104 }
3105 }
3106 HYPRE_IJVectorGetValues(HYb_, 1, &index, &ddata);
3107 HYPRE_IJVectorAddToValues(R1, 1, (const int *) &rowNum,
3108 (const double *) &ddata);
3109 rowNum++;
3110 }
3111 for ( i = localEndRow_-nConstraints_; i < localEndRow_; i++ )
3112 {
3113 HYPRE_IJVectorGetValues(HYb_, 1, &i, &ddata);
3114 HYPRE_IJVectorAddToValues(R1, 1, (const int *) &rowNum,
3115 (const double *) &ddata);
3116 rowNum++;
3117 }
3118
3119 //-------------------------------------------------------------
3120 // inv(A22) * (f2 - A21 * sol)
3121 //-------------------------------------------------------------
3122
3123 ierr = HYPRE_IJVectorCreate(comm_, startRow, startRow+x2NRows-1, &x2);
3124 ierr = HYPRE_IJVectorSetObjectType(x2, HYPRE_PARCSR);
3125 ierr = HYPRE_IJVectorInitialize(x2);
3126 ierr = HYPRE_IJVectorAssemble(x2);
3127 hypre_assert(!ierr);
3128 HYPRE_IJMatrixGetObject(HYinvA22_, (void **) &A22_csr);
3129 HYPRE_IJVectorGetObject(R1, (void **) &r_csr);
3130 HYPRE_IJVectorGetObject(x2, (void **) &x2_csr);
3131 HYPRE_ParCSRMatrixMatvec( 1.0, A22_csr, r_csr, 0.0, x2_csr );
3132
3133 //-------------------------------------------------------------
3134 // inject final solution to the solution vector
3135 //-------------------------------------------------------------
3136
3137 localNRows = localEndRow_ - localStartRow_ + 1 - 2 * nConstraints_;
3138 rowNum = localStartRow_ - 1;
3139 for ( i = startRow2; i < startRow2+localNRows; i++ )
3140 {
3141 HYPRE_IJVectorGetValues(reducedX_, 1, &i, &ddata);
3142 while (HYPRE_LSI_Search(selectedList_,rowNum,nConstraints_)>=0)
3143 rowNum++;
3144 HYPRE_IJVectorSetValues(HYx_, 1, (const int *) &rowNum,
3145 (const double *) &ddata);
3146 rowNum++;
3147 }
3148 for ( i = 0; i < nConstraints_; i++ )
3149 {
3150 for ( j = 0; j < nConstraints_; j++ )
3151 {
3152 if ( selectedListAux_[j] == i )
3153 {
3154 index = selectedList_[j];
3155 break;
3156 }
3157 }
3158 j = i + startRow;
3159 HYPRE_IJVectorGetValues(x2, 1, &j, &ddata);
3160 HYPRE_IJVectorSetValues(HYx_, 1, (const int *) &index,
3161 (const double *) &ddata);
3162 }
3163 for ( i = nConstraints_; i < 2*nConstraints_; i++ )
3164 {
3165 j = startRow + i;
3166 HYPRE_IJVectorGetValues(x2, 1, &j, &ddata);
3167 index = localEndRow_ - 2 * nConstraints_ + i;
3168 HYPRE_IJVectorSetValues(HYx_, 1, (const int *) &index,
3169 (const double *) &ddata);
3170 }
3171
3172 //-------------------------------------------------------------
3173 // residual norm check
3174 //-------------------------------------------------------------
3175
3176 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
3177 HYPRE_IJVectorGetObject(HYx_, (void **) &x_csr);
3178 HYPRE_IJVectorGetObject(HYb_, (void **) &b_csr);
3179 HYPRE_IJVectorGetObject(HYr_, (void **) &r_csr);
3180 HYPRE_ParVectorCopy( b_csr, r_csr );
3181 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
3182 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
3183 rnorm = sqrt( rnorm );
3184 if ( mypid_ == 0 && ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 ))
3185 printf("buildSlideReducedSoln::final residual norm = %e\n", rnorm);
3186 }
3187 currX_ = HYx_;
3188
3189 //****************************************************************
3190 // clean up
3191 //----------------------------------------------------------------
3192
3193 HYPRE_IJVectorDestroy(R1);
3194 HYPRE_IJVectorDestroy(x2);
3195 return rnorm;
3196 }
3197
3198 //*****************************************************************************
3199 // Given the matrix (A) within the object, compute the reduced system and put
3200 // it in place. Additional information given are :
3201 //-----------------------------------------------------------------------------
3202
buildSlideReducedSoln2()3203 double HYPRE_LinSysCore::buildSlideReducedSoln2()
3204 {
3205 int i, j, *int_array, *gint_array, x2NRows, x2GlobalNRows;
3206 int ierr, rowNum, startRow, startRow2, index, localNRows;
3207 int index2;
3208 double ddata, rnorm;
3209 HYPRE_ParCSRMatrix A_csr, A21_csr, A22_csr;
3210 HYPRE_ParVector x_csr, x2_csr, r_csr, b_csr;
3211 HYPRE_IJVector R1, x2;
3212
3213 if ( HYA21_ == NULL || HYinvA22_ == NULL )
3214 {
3215 printf("buildSlideReducedSoln2 WARNING : A21 or A22 absent.\n");
3216 return (0.0);
3217 }
3218 else
3219 {
3220 //------------------------------------------------------------------
3221 // compute A21 * sol
3222 //------------------------------------------------------------------
3223
3224 x2NRows = 2 * nConstraints_;
3225 int_array = new int[numProcs_];
3226 gint_array = new int[numProcs_];
3227 for ( i = 0; i < numProcs_; i++ ) int_array[i] = 0;
3228 int_array[mypid_] = x2NRows;
3229 MPI_Allreduce(int_array,gint_array,numProcs_,MPI_INT,MPI_SUM,comm_);
3230 x2GlobalNRows = 0;
3231 for ( i = 0; i < numProcs_; i++ ) x2GlobalNRows += gint_array[i];
3232 rowNum = 0;
3233 for ( i = 0; i < mypid_; i++ ) rowNum += gint_array[i];
3234 startRow = rowNum;
3235 startRow2 = localStartRow_ - 1 - rowNum/2;
3236 delete [] int_array;
3237 delete [] gint_array;
3238 ierr = HYPRE_IJVectorCreate(comm_, startRow, startRow+x2NRows-1, &R1);
3239 ierr = HYPRE_IJVectorSetObjectType(R1, HYPRE_PARCSR);
3240 ierr = HYPRE_IJVectorInitialize(R1);
3241 ierr = HYPRE_IJVectorAssemble(R1);
3242 hypre_assert(!ierr);
3243 HYPRE_IJMatrixGetObject(HYA21_, (void **) &A21_csr);
3244 HYPRE_IJVectorGetObject(currX_, (void **) &x_csr);
3245 HYPRE_IJVectorGetObject(R1, (void **) &r_csr);
3246 HYPRE_ParCSRMatrixMatvec( -1.0, A21_csr, x_csr, 0.0, r_csr );
3247
3248 //------------------------------------------------------------------
3249 // f2 - A21 * sol
3250 //------------------------------------------------------------------
3251
3252 for ( i = 0; i < nConstraints_; i++ )
3253 {
3254 for ( j = 0; j < nConstraints_; j++ )
3255 {
3256 if ( selectedListAux_[j] == i )
3257 {
3258 index = selectedList_[j];
3259 break;
3260 }
3261 }
3262 HYPRE_IJVectorGetValues(HYb_, 1, &index, &ddata);
3263 HYPRE_IJVectorAddToValues(R1, 1, (const int *) &rowNum,
3264 (const double *) &ddata);
3265 rowNum++;
3266 }
3267 for ( i = localEndRow_-nConstraints_; i < localEndRow_; i++ )
3268 {
3269 HYPRE_IJVectorGetValues(HYb_, 1, &i, &ddata);
3270 HYPRE_IJVectorAddToValues(R1, 1, (const int *) &rowNum,
3271 (const double *) &ddata);
3272 rowNum++;
3273 }
3274
3275 //-------------------------------------------------------------
3276 // inv(A22) * (f2 - A21 * sol)
3277 //-------------------------------------------------------------
3278
3279 ierr = HYPRE_IJVectorCreate(comm_, startRow, startRow+x2NRows-1, &x2);
3280 ierr = HYPRE_IJVectorSetObjectType(x2, HYPRE_PARCSR);
3281 ierr = HYPRE_IJVectorInitialize(x2);
3282 ierr = HYPRE_IJVectorAssemble(x2);
3283 hypre_assert(!ierr);
3284 HYPRE_IJMatrixGetObject(HYinvA22_, (void **) &A22_csr );
3285 HYPRE_IJVectorGetObject(R1, (void **) &r_csr );
3286 HYPRE_IJVectorGetObject(x2, (void **) &x2_csr );
3287 HYPRE_ParCSRMatrixMatvec( 1.0, A22_csr, r_csr, 0.0, x2_csr );
3288
3289 //-------------------------------------------------------------
3290 // inject final solution to the solution vector
3291 //-------------------------------------------------------------
3292
3293 localNRows = localEndRow_ - localStartRow_ + 1 - nConstraints_;
3294 for ( i = 0; i < localNRows; i++ )
3295 {
3296 index = startRow2 + i;
3297 HYPRE_IJVectorGetValues(reducedX_, 1, &index, &ddata);
3298 index2 = localStartRow_ - 1 + i;
3299 HYPRE_IJVectorSetValues(HYx_, 1, (const int *) &index2,
3300 (const double *) &ddata);
3301 }
3302 for ( i = 0; i < nConstraints_; i++ )
3303 {
3304 for ( j = 0; j < nConstraints_; j++ )
3305 {
3306 if ( selectedListAux_[j] == i )
3307 {
3308 index = selectedList_[j];
3309 break;
3310 }
3311 }
3312 j = i + startRow;
3313 HYPRE_IJVectorGetValues(x2, 1, &j, &ddata);
3314 HYPRE_IJVectorSetValues(HYx_, 1, (const int *) &index,
3315 (const double *) &ddata);
3316 }
3317 for ( i = nConstraints_; i < 2*nConstraints_; i++ )
3318 {
3319 j = startRow + i;
3320 HYPRE_IJVectorGetValues(x2, 1, &j, &ddata);
3321 index = localEndRow_ - 2 * nConstraints_ + i;
3322 HYPRE_IJVectorSetValues(HYx_, 1, (const int *) &index,
3323 (const double *) &ddata);
3324 }
3325
3326 //-------------------------------------------------------------
3327 // residual norm check
3328 //-------------------------------------------------------------
3329
3330 HYPRE_IJMatrixGetObject(HYA_, (void **) &A_csr);
3331 HYPRE_IJVectorGetObject(HYx_, (void **) &x_csr);
3332 HYPRE_IJVectorGetObject(HYb_, (void **) &b_csr);
3333 HYPRE_IJVectorGetObject(HYr_, (void **) &r_csr);
3334 HYPRE_ParVectorCopy( b_csr, r_csr );
3335 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
3336 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
3337 rnorm = sqrt( rnorm );
3338 if ( mypid_ == 0 && ( HYOutputLevel_ & HYFEI_SLIDEREDUCE1 ))
3339 printf("buildSlideReducedSoln::final residual norm = %e\n", rnorm);
3340 }
3341 currX_ = HYx_;
3342
3343 //****************************************************************
3344 // clean up
3345 //----------------------------------------------------------------
3346
3347 HYPRE_IJVectorDestroy(R1);
3348 HYPRE_IJVectorDestroy(x2);
3349 return rnorm;
3350 }
3351
3352