FEI Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
fei_AztecDVBR_Matrix.cpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ************************************************************************
4// FEI: Finite Element Interface to Linear Solvers
5// Copyright (2005) Sandia Corporation.
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8// U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Alan Williams (william@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41*/
42
43
45#include <fei_iostream.hpp>
46
47#ifdef HAVE_FEI_AZTECOO
48
49#include <assert.h>
50#include <stdlib.h>
51#include <stdio.h>
52#include <cstring>
53
54#include <fei_mpi.h>
55
56#ifndef FEI_SER
57
58#define AZTEC_MPI AZTEC_MPI
59#define AZ_MPI AZ_MPI
60#ifndef MPI
61#define MPI MPI
62#endif
63
64#endif
65
66#include <az_aztec.h>
67#include <fei_Aztec_Map.hpp>
71
72namespace fei_trilinos {
73
74 //==============================================================================
76 : amap_(map),
77 Amat_(NULL),
78 N_update_(map->getNumLocalBlocks()),
79 external_(NULL),
80 extern_index_(NULL),
81 update_index_(NULL),
82 data_org_(NULL),
83 orderingUpdate_(NULL),
84 isLoaded_(false),
85 isAllocated_(false),
86 localNNZ_(0),
87 nnzPerRow_(NULL),
88 numRemoteBlocks_(0),
89 remoteInds_(NULL),
90 remoteBlockSizes_(NULL)
91 {
92 nnzPerRow_ = new int[N_update_];
93
94 for(int i=0; i<N_update_; i++) {
95 nnzPerRow_[i] = 0;
96 }
97
98 Amat_ = AZ_matrix_create(N_update_);
99 Amat_->matrix_type = AZ_VBR_MATRIX;
100 Amat_->matvec =
101 (void(*)(double*,double*,AZ_MATRIX_STRUCT*,int*))AZ_VBR_matvec_mult;
102
103 //now we can allocate and fill the rpntr array.
104 Amat_->rpntr = NULL;
105 calcRpntr();
106 }
107
108 //==============================================================================
109 AztecDVBR_Matrix::AztecDVBR_Matrix(const AztecDVBR_Matrix& src)
110 : amap_(src.amap_),
111 Amat_(NULL),
112 N_update_(src.N_update_),
113 external_(NULL),
114 extern_index_(NULL),
115 update_index_(NULL),
116 data_org_(NULL),
117 orderingUpdate_(NULL),
118 isLoaded_(src.isLoaded_),
119 isAllocated_(src.isAllocated_),
120 localNNZ_(src.localNNZ_),
121 nnzPerRow_(NULL),
122 numRemoteBlocks_(0),
123 remoteInds_(NULL),
124 remoteBlockSizes_(NULL)
125 {
126 //
127 //This copy constructor just takes a reference to src's amap_ (above), but
128 //'deep'-copies everything else. i.e., all arrays are allocated here and copied
129 //from src, etc.
130 //
131 //When this constructor completes, this matrix should be in the same state as
132 //the 'src' matrix. i.e., if src is already allocated, then this matrix will
133 //have its structure allocated. If src is already loaded (AZ_transform'd) then
134 //this matrix will have all arrays resulting from AZ_transform allocated too.
135 //The only thing this matrix won't get is the coefficient data from src.
136 //
137 nnzPerRow_ = new int[N_update_];
138
139 int i;
140 for(i=0; i<N_update_; i++) {
141 nnzPerRow_[i] = src.nnzPerRow_[i];
142 }
143
144 Amat_ = AZ_matrix_create(N_update_);
145 Amat_->matrix_type = AZ_VBR_MATRIX;
146 Amat_->matvec =
147 (void(*)(double*,double*,AZ_MATRIX_STRUCT*,int*))AZ_VBR_matvec_mult;
148
149 Amat_->rpntr = NULL;
150 calcRpntr();
151
152 if (isAllocated_) {
153 Amat_->bpntr = new int[N_update_+1];
154 for(i=0; i<N_update_+1; i++) Amat_->bpntr[i] = src.Amat_->bpntr[i];
155
156 int totalNNZBlks = Amat_->bpntr[N_update_];
157 Amat_->bindx = new int[totalNNZBlks];
158 for(i=0; i<totalNNZBlks; i++) Amat_->bindx[i] = src.Amat_->bindx[i];
159
160 Amat_->indx = new int[totalNNZBlks+1];
161 for(i=0; i<totalNNZBlks+1; i++) Amat_->indx[i] = src.Amat_->indx[i];
162
163 Amat_->val = new double[localNNZ_];
164 for(i=0; i<localNNZ_; i++) Amat_->val[i] = 0.0;
165 }
166
167 if (isLoaded_) {
168 int dataOrgLength = src.data_org_[AZ_total_send] + AZ_send_list;
169 data_org_ = (int*) AZ_allocate(dataOrgLength * sizeof(int));
170 for(i=0; i<dataOrgLength; i++) data_org_[i] = src.data_org_[i];
171
172 Amat_->data_org = data_org_;
173
174 int extLength = src.data_org_[AZ_N_ext_blk];
175 external_ = (int*) AZ_allocate(extLength * sizeof(int));
176 extern_index_ = (int*) AZ_allocate(extLength * sizeof(int));
177 for(i=0; i<extLength; i++) {
178 external_[i] = src.external_[i];
179 extern_index_[i] = src.extern_index_[i];
180 }
181
182 update_index_ = (int*) AZ_allocate(N_update_ * sizeof(int));
183 orderingUpdate_ = new int[N_update_];
184 for(i=0; i<N_update_; i++) {
185 update_index_[i] = src.update_index_[i];
186 orderingUpdate_[i] = src.orderingUpdate_[i];
187 }
188
189 int cpntrLength = N_update_ + src.data_org_[AZ_N_ext_blk] + 1;
190 Amat_->cpntr = (int*) AZ_allocate(cpntrLength * sizeof(int));
191 for(i=0; i<cpntrLength; i++) Amat_->cpntr[i] = src.Amat_->cpntr[i];
192 }
193 }
194
195 //==============================================================================
196 AztecDVBR_Matrix::~AztecDVBR_Matrix(){
197
198 if (isAllocated()) {
199 delete [] Amat_->val;
200 delete [] Amat_->bindx;
201 delete [] Amat_->rpntr;
202 delete [] Amat_->bpntr;
203 delete [] Amat_->indx;
204
205 delete [] remoteInds_;
206 delete [] remoteBlockSizes_;
207
208 setAllocated(false);
209 }
210
211 if (isLoaded()) {
212 free(Amat_->cpntr);
213 free(external_);
214 free(extern_index_);
215 free(update_index_);
216 free(data_org_);
217 delete [] orderingUpdate_;
218
219 setLoaded(false);
220 }
221
222 delete [] nnzPerRow_;
223 localNNZ_ = 0;
224
225 AZ_matrix_destroy(&Amat_);
226 Amat_ = NULL;
227 }
228
229 //==============================================================================
230 int AztecDVBR_Matrix::getNumBlocksPerRow(int blkRow, int& nnzBlksPerRow) const {
231 //
232 //On return, nnzBlksPerRow will be the number of nonzero blocks in row blkRow.
233 //
234 if (!isAllocated()) return(1);
235
236 int index;
237
238 if (!inUpdate(blkRow, index)) {
239 fei::console_out() << "AztecDVBR_Matrix::getNumBlocksPerRow: ERROR: blkRow "
240 << blkRow << " not in local update list." << FEI_ENDL;
241 return(1);
242 }
243
244 nnzBlksPerRow = Amat_->bpntr[index+1] - Amat_->bpntr[index];
245
246 return(0);
247 }
248
249 //==============================================================================
250 int AztecDVBR_Matrix::getNumBlocksPerRow(int* nnzBlksPerRow) const {
251 //
252 //nnzBlksPerRow must be allocated by the calling code.
253 //
254 //nnzBlksPerRow is a list of length number-of-local-block-rows, and
255 //nnzBlksPerRow[i] gives the number of nonzeros blocks in row i.
256 //
257 if (!isAllocated()) return(1);
258
259 for(int i=0; i<amap_->getNumLocalBlocks(); i++) {
260 nnzBlksPerRow[i] = Amat_->bpntr[i+1] - Amat_->bpntr[i];
261 }
262
263 return(0);
264 }
265
266 //==============================================================================
267 int AztecDVBR_Matrix::getNumNonzerosPerRow(int blkRow, int& nnzPerRow) const {
268 //
269 //This function finds nnzPerRow, the number of nonzero *point* entries for
270 //row 'blkRow'.
271 //
272 if (!isAllocated()) return(1);
273
274 int index;
275
276 if (!inUpdate(blkRow, index)) {
277 fei::console_out() << "AztecDVBR_Matrix::getNumNonzerosPerRow: ERROR: blkRow "
278 << blkRow << " not in local update list." << FEI_ENDL;
279 return(1);
280 }
281
282 nnzPerRow = nnzPerRow_[index];
283
284 return(0);
285 }
286
287 //==============================================================================
288 int AztecDVBR_Matrix::getNumNonzerosPerRow(int* nnzPerRow) const {
289 //
290 //nnzPerRow must be allocated by the calling code,
291 //length number-of-local-block-rows.
292 //
293 //This function fills nnzPerRow so that nnzPerRow[i] gives the
294 //number of nonzero *point* entries for row i.
295 //
296 if (!isAllocated()) return(1);
297
298 for(int i=0; i<amap_->getNumLocalBlocks(); i++) {
299 nnzPerRow[i] = nnzPerRow_[i];
300 }
301
302 return(0);
303 }
304
305 //==============================================================================
306 int AztecDVBR_Matrix::getBlockSize(int blkRow, int blkCol,
307 int& ptRows, int& ptCols) {
308 int index;
309
310 ptRows = 0;
311 ptCols = 0;
312
313 if (!inUpdate(blkRow, index)) {
314 fei::console_out() << "AztecDVBR_Matrix::getBlockSize: ERROR: blkRow "
315 << blkRow << " not in local update list." << FEI_ENDL;
316 return(1);
317 }
318
319 ptRows = Amat_->rpntr[index+1] - Amat_->rpntr[index];
320
321 int local = inUpdate(blkCol, index);
322
323 if (local) {
324 ptCols = Amat_->rpntr[index+1] - Amat_->rpntr[index];
325 }
326 else {
327 index = AZ_find_index(blkCol, remoteInds_, numRemoteBlocks_);
328
329 if (index < 0) return(1);
330 else ptCols = remoteBlockSizes_[index];
331 }
332
333 return(0);
334 }
335
336 //==============================================================================
337 void AztecDVBR_Matrix::matvec(const Aztec_LSVector& x, Aztec_LSVector& y) const {
338
339 // AztecDVBR_Matrix::matvec --- form y = Ax
340
341 assert(isLoaded());
342
343 int *proc_config = amap_->getProcConfig();
344 double *b = (double*)x.startPointer();
345 double *c = (double*)y.startPointer();
346
347 AZ_VBR_matvec_mult(b, c, Amat_, proc_config);
348
349 return;
350 }
351
352 //==============================================================================
353 void AztecDVBR_Matrix::put(double s){
354
355 if (!isAllocated()) return;
356
357 for(int i=0; i<localNNZ_; i++) {
358 Amat_->val[i] = s;
359 }
360
361 return;
362 }
363
364 //==============================================================================
365 int AztecDVBR_Matrix::getBlockRow(int blkRow,
366 double* val,
367 int* blkColInds,
368 int numNzBlks) const {
369
370 if (!isAllocated()) return(1);
371
372 int index;
373
374 if (!inUpdate(blkRow, index)) {
375 fei::console_out() << "AztecDVBR_Matrix::getBlockRow: ERROR: blkRow "
376 << blkRow << " not in local update list." << FEI_ENDL;
377 return(1);
378 }
379
380 //for each block, we need to find its block column index
381 //in the bindx array, then go to that same position in the indx
382 //array to find out how many point-entries are in that block.
383 //We can then use the indx entry to go to the val array and get
384 //the data.
385
386 int nnzBlks = 0, nnzPts = 0;
387 int err = getNumBlocksPerRow(blkRow, nnzBlks);
388 if (err) return(err);
389 err = getNumNonzerosPerRow(blkRow, nnzPts);
390 if (err) return(err);
391
392 if (numNzBlks != nnzBlks) return(1);
393
394 int offset = 0;
395 int blkCounter = 0;
396 const int* blkUpdate = amap_->getBlockUpdate();
397 for(int indb = Amat_->bpntr[index]; indb<Amat_->bpntr[index+1]; indb++) {
398
399 int numEntries = Amat_->indx[indb+1] - Amat_->indx[indb];
400 int valOffset = Amat_->indx[indb];
401
402 if (isLoaded()) {
403 int ind = Amat_->bindx[indb];
404 if (ind < N_update_) {
405 blkColInds[blkCounter++] = blkUpdate[orderingUpdate_[ind]];
406 }
407 else {
408 blkColInds[blkCounter++] = external_[ind-N_update_];
409 }
410 }
411 else {
412 blkColInds[blkCounter++] = Amat_->bindx[indb];
413 }
414
415 //ok, now we're ready to get the stuff.
416 for(int i=0; i<numEntries; i++) {
417 val[offset + i] = Amat_->val[valOffset + i];
418 }
419
420 offset += numEntries;
421 }
422
423 return(0);
424 }
425
426 //==============================================================================
427 int AztecDVBR_Matrix::putBlockRow(int blkRow,
428 double* val,
429 int* blkColInds,
430 int numNzBlks) const {
431
432 if (!isAllocated()) return(1);
433
434 int index;
435
436 if (!inUpdate(blkRow, index)) {
437 fei::console_out() << "AztecDVBR_Matrix::putBlockRow: ERROR: blkRow "
438 << blkRow << " not in local update list." << FEI_ENDL;
439 return(1);
440 }
441
442 //for each incoming block, we need to find its block column index
443 //in the bindx array, then go to that same position in the indx
444 //array to find out how many (point) entries are in that block.
445 //We can then use the indx entry to go to the val array and store
446 //the data.
447
448 int offset = 0;
449 for(int blk = 0; blk<numNzBlks; blk++) {
450 int indb = getBindxOffset(blkColInds[blk],
451 Amat_->bpntr[index], Amat_->bpntr[index+1]-1);
452
453 if (indb < 0) messageAbort("putBlockRow: blk col not found in row.");
454
455 int numEntries = Amat_->indx[indb+1] - Amat_->indx[indb];
456 int valOffset = Amat_->indx[indb];
457
458 //ok, now we're ready to store the stuff.
459 for(int i=0; i<numEntries; i++) {
460 Amat_->val[valOffset + i] = val[offset + i];
461 }
462
463 offset += numEntries;
464 }
465
466 return(0);
467 }
468
469 //==============================================================================
470 int AztecDVBR_Matrix::sumIntoBlockRow(int blkRow,
471 double* val,
472 int* blkColInds,
473 int numNzBlks) const
474 {
475 //
476 //This function is the same as putBlockRow, except the values
477 //are summed into any existing values rather than overwriting
478 //them.
479 //
480 if (!isAllocated()) return(1);
481
482 int index;
483
484 if (!inUpdate(blkRow, index)) {
485 fei::console_out() << "AztecDVBR_Matrix::sumIntoBlockRow: ERROR: blkRow "
486 << blkRow << " not in local update list." << FEI_ENDL;
487 return(1);
488 }
489
490 //for each incoming block, we need to find its block column index
491 //in the bindx array, then go to that same position in the indx
492 //array to find out how many (point) entries are in that block.
493 //We can then use the indx entry to go to the val array and store
494 //the data.
495
496 int offset = 0;
497
498 for(int blk = 0; blk<numNzBlks; blk++) {
499 int indb = getBindxOffset(blkColInds[blk],
500 Amat_->bpntr[index], Amat_->bpntr[index+1]-1);
501
502 if (indb < 0) {
503 fei::console_out() << "AztecDVBR_Matrix::sumIntoBlockRow: blk col "
504 << blkColInds[blk] << " not found in row " << blkRow << FEI_ENDL;
505 abort();
506 }
507
508 int numEntries = Amat_->indx[indb+1] - Amat_->indx[indb];
509 int valOffset = Amat_->indx[indb];
510
511 //ok, now we're ready to store the stuff.
512 for(int i=0; i<numEntries; i++) {
513 Amat_->val[valOffset + i] += val[offset + i];
514 }
515
516 offset += numEntries;
517 }
518
519 return(0);
520 }
521
522 //==============================================================================
523 void AztecDVBR_Matrix::allocate(int* numNzBlks, int* blkColInds) {
524 //
525 // This function builds the structure of the matrix. i.e., does the
526 // memory allocation, etc.
527 //
528
529 //calculate the bpntr array, which holds info about the number of
530 //nonzero blocks per row.
531 calcBpntr(numNzBlks);
532
533 //we can now get the total number of nonzero blocks from the last
534 //entry in bpntr.
535 int totalNumNzBlks = Amat_->bpntr[N_update_];
536
537 //now we can set the bindx array, which holds block column indices.
538 setBindx(totalNumNzBlks, blkColInds);
539
540 //and now we're ready to allocate and fill the indx array, which
541 //holds info on the number of point entries in each nonzero block.
542 calcIndx(totalNumNzBlks);
543
544 //the last thing we need to do is allocate and initialize the val array.
545 Amat_->val = new double[localNNZ_];
546
547 for(int i=0; i<localNNZ_; i++) {
548 Amat_->val[i] = 0.0;
549 }
550
551 setAllocated(true);
552 return;
553 }
554
555 //==============================================================================
556 void AztecDVBR_Matrix::loadComplete() {
557 //
558 // This is where we call the Aztec function AZ_transform, which calculates
559 // communication parameters and re-orders the equations for use as a
560 // global distributed matrix.
561 //
562 MPI_Comm thisComm = amap_->getCommunicator();
563
564 // Sync processors.
565
566 MPI_Barrier(thisComm);
567
568#ifndef FEI_SER
569 int thisProc = 0;
570 MPI_Comm_rank(thisComm, &thisProc);
571#endif
572
573 AZ_transform(amap_->getProcConfig(), &external_, Amat_->bindx, Amat_->val,
574 amap_->getBlockUpdate(), &update_index_, &extern_index_, &data_org_,
575 N_update_, Amat_->indx, Amat_->bpntr, Amat_->rpntr,
576 &(Amat_->cpntr), AZ_VBR_MATRIX);
577
578 data_org_[AZ_internal_use] = 1;
579
580 Amat_->data_org = data_org_;
581
582 //On return from AZ_transform, the array update_index_ contains a mapping
583 //to the local re-ordering of the indices of the update_ array. Now we will
584 //fill the orderingUpdate array with the reverse of that mapping. i.e., a
585 //record of how to get back to the original ordering of the update indices.
586
587 orderingUpdate_ = new int[N_update_];
588 for(int ii=0; ii<N_update_; ii++)
589 orderingUpdate_[update_index_[ii]] = ii;
590
591 // Sync processors.
592#ifndef FEI_SER
593 MPI_Barrier(thisComm);
594#endif
595
596 setLoaded(true);
597
598 return;
599 }
600
601 //==============================================================================
602 bool AztecDVBR_Matrix::readFromFile(const char *filename){
603 //
604 //readFromFile should be able to be called after the matrix is constructed,
605 //and before allocate has been called. i.e., example usage should include:
606 //
607 // AztecDVBR_Matrix A(map, update);
608 // A.readFromFile(fileName);
609 // A.matvec(b, c);
610 //
611 //i.e., readFromFile can take the place of the allocate and loadComplete
612 //calls.
613 //
614 FILE *infile = NULL;
615 MPI_Comm thisComm = amap_->getCommunicator();
616
617 MPI_Barrier(thisComm);
618
619 infile = fopen(filename, "r");
620 if (!infile) messageAbort("readFromFile: couldn't open file.");
621
622 int* num_nz_blocks = NULL;
623 int* blk_col_inds = NULL;
624
625 readAllocateInfo(infile, num_nz_blocks, blk_col_inds);
626
627 allocate(num_nz_blocks, blk_col_inds);
628
629 delete [] num_nz_blocks;
630 delete [] blk_col_inds;
631
632 fclose(infile);
633 infile = fopen(filename, "r");
634
635 readMatrixData(infile);
636
637 fclose(infile);
638
639 loadComplete();
640
641 return(true);
642 }
643
644 //==============================================================================
645 void AztecDVBR_Matrix::readAllocateInfo(FILE* infile,
646 int*& num_nz_blocks,
647 int*& blk_col_inds) {
648 //
649 //This function will read through infile and construct the lists
650 //num_nz_blocks (which is the number of nonzero blocks per row) and
651 //blk_col_inds (which is the block-column indices of those blocks).
652 //
653 //It is assumed that these two lists are empty when this function is
654 //called.
655
656 int i;
657
658 if (num_nz_blocks) delete [] num_nz_blocks;
659 if (blk_col_inds) delete [] blk_col_inds;
660
661 num_nz_blocks = new int[N_update_];
662
663 //we'll use a 2-D array for constructing the set of block column indices,
664 //because we need to keep them grouped by rows, and we aren't guaranteed
665 //that they'll be grouped by rows in the file.
666 int totalNumBlks = 0;
667 int** blkColInds = new int*[N_update_];
668
669 for(i=0; i<N_update_; i++) {
670 num_nz_blocks[i] = 0;
671 blkColInds[i] = NULL;
672 }
673
674 int blkRows, blkCols, rows, cols;
675 char line[256];
676
677 do {
678 fgets(line,256,infile);
679 } while(strchr(line,'%'));
680 sscanf(line,"%d %d %d %d",&blkRows, &blkCols, &rows, &cols);
681
682 if ((blkRows != blkCols) || (rows != cols))
683 messageAbort("readAllocateInfo: non-square matrix not allowed.");
684
685 int br, bc, pr, pc, index;
686
687 while (!feof(infile)) {
688 do {
689 fgets(line,256,infile);
690 } while(strchr(line,'%'));
691
692 if(feof(infile))break;
693
694 sscanf(line, "%d %d %d %d", &br, &bc, &pr, &pc);
695
696 if (inUpdate(br, index)) {
697 if ((bc < 0) || bc >= blkCols) {
698 char mesg[80];
699 sprintf(mesg,"readAllocateInfo: blkCols %d, 0-based col ind %d",
700 blkCols, bc);
701 fclose(infile);
702 messageAbort(mesg);
703 }
704 insertList(bc, blkColInds[index], num_nz_blocks[index]);
705 totalNumBlks++;
706 }
707 }
708
709 //so we've read the whole file, now flatten the 2-D list blkColInds
710 //into the required 1-D list blk_col_inds.
711 blk_col_inds = new int[totalNumBlks];
712
713 int offset = 0;
714 for(i=0; i<N_update_; i++) {
715 for(int j=0; j<num_nz_blocks[i]; j++) {
716 blk_col_inds[offset++] = blkColInds[i][j];
717 }
718
719 delete [] blkColInds[i];
720 }
721
722 delete [] blkColInds;
723 }
724
725 //==============================================================================
726 void AztecDVBR_Matrix::readMatrixData(FILE* infile) {
727
728 int blkRows, blkCols, rows, cols;
729 int br, bc, pr, pc, nnz, index;
730 double* blockValues = NULL;
731 char line[256];
732
733 do {
734 fgets(line,256,infile);
735 } while(strchr(line,'%'));
736 sscanf(line,"%d %d %d %d",&blkRows, &blkCols, &rows, &cols);
737
738 while (!feof(infile)) {
739 do {
740 fgets(line,256,infile);
741 } while(strchr(line,'%'));
742
743 if(feof(infile))break;
744
745 sscanf(line, "%d %d %d %d %d", &br, &bc, &pr, &pc, &nnz);
746
747 if (inUpdate(br, index)){
748 //br (block-row) is in the local row-space, so let's
749 //plug this block into our matrix data structure.
750
751 blockValues = new double[nnz];
752 getValuesFromString(line, std::strlen(line)+1, blockValues, nnz);
753
754 putBlockRow(br, blockValues, &bc, 1);
755
756 delete [] blockValues;
757 }
758 }
759 }
760
762 void AztecDVBR_Matrix::getValuesFromString(char *line, int len, double *values,
763 int lenValues){
764 (void)len;
765 int i;
766
767 //first, we know that 'line' contains 5 integers at the beginning, separated
768 //by spaces, which we want to jump over. So we need the offset of the 5th
769 //space in 'line'.
770
771 char *offset = &line[0];
772 for(i=0; i<5; i++){
773 offset = strchr(offset, ' ');
774 offset++;
775 }
776
777 //now we're ready to pick out the numbers to put into the values array.
778 for(i=0; i<lenValues; i++){
779 sscanf(offset,"%le",&(values[i]));
780 offset = strchr(offset, ' ');
781 offset++;
782 }
783
784 return;
785 }
786
788 bool AztecDVBR_Matrix::writeToFile(const char *fileName) const {
789
790 int thisProc = amap_->getProcConfig()[AZ_node];
791 int numProcs = amap_->getProcConfig()[AZ_N_procs];
792 MPI_Comm thisComm = amap_->getCommunicator();
793
794 int numGlobalBlocks = amap_->getNumGlobalBlocks();
795 int numGlobalPtEqns = amap_->globalSize();
796
797 int numLocalBlocks = N_update_;
798
799 FILE* file = NULL;
800
801 for(int p=0; p<numProcs; p++) {
802 MPI_Barrier(thisComm);
803
804 if (p == thisProc) {
805 if (thisProc == 0) {
806 //open the file for writing and write the first line, which is
807 //num-blk-rows num-block-cols num-pt-rows num-pt-cols
808
809 file = fopen(fileName, "w");
810 fprintf(file, "%d %d %d %d\n", numGlobalBlocks, numGlobalBlocks,
811 numGlobalPtEqns, numGlobalPtEqns);
812 }
813 else {
814 //open the file for appending
815 file = fopen(fileName, "a");
816 }
817
818 //now loop over the local portion of the matrix, writing it to file,
819 //one nonzer-block per line. Each line of the file will contain these
820 //numbers separated by spaces:
821 //blk-row-index
822 //blk-col-index
823 //blk-row-size (number of pt-rows in this block-row)
824 //nnz (number of nonzero point-entries in this block-entry)
825 //nonzero1 nonzero2 ... nonzero<nnz> (the nonzeros for this block)
826
827 for(int brow=0; brow<numLocalBlocks; brow++) {
828 int bcolind1 = Amat_->bpntr[brow];
829 int bcolind2 = Amat_->bpntr[brow+1];
830
831 for(int ind=bcolind1; ind<bcolind2; ind++) {
832 int nnzPts = Amat_->indx[ind+1] - Amat_->indx[ind];
833
834 if (nnzPts <= 0) continue;
835
836 int blkRowSize = Amat_->rpntr[brow+1]-Amat_->rpntr[brow];
837
838 int globCol = -1;
839 int lookup = Amat_->bindx[ind];
840 if (isLoaded()) {
841 if (lookup < N_update_) {
842 globCol = amap_->getBlockUpdate()[orderingUpdate_[lookup]];
843 }
844 else {
845 globCol = external_[lookup-N_update_];
846 }
847 }
848 else {
849 globCol = lookup;
850 }
851
852 int globalRow = amap_->getBlockUpdate()[brow];
853 if (isLoaded()) globalRow = amap_->getBlockUpdate()[orderingUpdate_[brow]];
854
855 fprintf(file, "%d %d %d %d ", globalRow, globCol,
856 blkRowSize, nnzPts);
857
858 int offset = Amat_->indx[ind];
859 for(int i=0; i<nnzPts; i++) {
860 fprintf(file, "%20.13e ", Amat_->val[offset+i]);
861 }
862 fprintf(file, "\n");
863 }
864 }
865
866 fclose(file);
867 }
868 }
869
870 return(true);
871 }
872
873 //==============================================================================
874 int AztecDVBR_Matrix::inUpdate(int globalIndex, int& localIndex) const {
875 //
876 // This function determines whether globalIndex is in the local update set,
877 // and if it is, returns in localIndex the local index for it. If update_index_
878 // has already been allocated and set (by AZ_transform) then localIndex is
879 // taken from there.
880 // If globalIndex is not in the update set, inUpdate returns 0.
881 //
882 localIndex = AZ_find_index(globalIndex, amap_->getBlockUpdate(), N_update_);
883
884 if(localIndex==-1)return(0);
885
886 if(isLoaded_){
887 localIndex = update_index_[localIndex];
888 }
889
890 return(1);
891 }
892
893 //==============================================================================
894 void AztecDVBR_Matrix::calcRpntr() {
895 //
896 //This function will use information from the Aztec_BlockMap 'amap_'
897 //to set the Amat_->rpntr array.
898 //
899 //rpntr[0..M] (where M = number-of-blocks)
900 //rpntr[0] = 0
901 //rpntr[k+1] - rpntr[k] = size of block k
902 //
903 const int* blkSizes = amap_->getBlockSizes();
904
905 Amat_->rpntr = new int[N_update_+1];
906
907 Amat_->rpntr[0] = 0;
908
909 for(int i=0; i<N_update_; i++) {
910 Amat_->rpntr[i+1] = Amat_->rpntr[i] + blkSizes[i];
911 if (blkSizes[i] < 0)
912 messageAbort("allocate: negative block size.");
913 }
914 }
915
916 //==============================================================================
917 void AztecDVBR_Matrix::calcBpntr(int* numNzBlks) {
918 //
919 //This function will set the Amat_->bpntr array.
920 //
921 //bpntr[0..M] (where M = number-of-blocks)
922 //bpntr[0] = 0
923 //bpntr[k+1]-bpntr[k] = number of nonzero blocks in row k
924 //
925 Amat_->bpntr = new int[N_update_+1];
926
927 Amat_->bpntr[0] = 0;
928
929 for(int i=0; i<N_update_; i++) {
930 Amat_->bpntr[i+1] = Amat_->bpntr[i] + numNzBlks[i];
931 }
932 }
933
934 //==============================================================================
935 void AztecDVBR_Matrix::setBindx(int nnzBlks, int* blkColInds) {
936 //
937 //This function simply allocates and fills the Amat_->bindx array.
938 //
939 Amat_->bindx = new int[nnzBlks];
940
941 for(int i=0; i<nnzBlks; i++) {
942 Amat_->bindx[i] = blkColInds[i];
943 if (blkColInds[i] < 0)
944 messageAbort("setBindx: negative block col index.");
945 }
946 }
947
948 //==============================================================================
949 void AztecDVBR_Matrix::messageAbort(const char* mesg) const {
950 fei::console_out() << "AztecDVBR_Matrix: ERROR: " << mesg << " Aborting." << FEI_ENDL;
951 abort();
952 }
953
954 //==============================================================================
955 void AztecDVBR_Matrix::calcIndx(int nnzBlks) {
956 //
957 //This function allocates and fills the Amat_->indx array, which holds info
958 //on the number of entries in each nonzero block.
959 //
960 //indx[0..bpntr[M]], (where M = number of local block rows)
961 //indx[0] = 0
962 //indx[k+1]-indx[k] = number of entries in nonzero block k
963 //
964
965 Amat_->indx = new int[nnzBlks+1];
966
967 //we need to obtain block sizes for all local nonzero blocks. rpntr
968 //gives us the sizes for the blocks with column indices in the local
969 //update set, but we'll have to do some message passing to obtain the
970 //sizes of blocks with column indices in other procs' update sets.
971
972 int numProcs = amap_->getProcConfig()[AZ_N_procs];
973
974 if (numProcs > 1) {
975 //form a list of the column indices that are not local.
976 calcRemoteInds(remoteInds_, numRemoteBlocks_);
977
978 //now get sizes of blocks that correspond to remote rows.
979 remoteBlockSizes_ = new int[numRemoteBlocks_];
980 getRemoteBlkSizes(remoteBlockSizes_, remoteInds_, numRemoteBlocks_);
981 }
982
983 //now we're ready to set the block sizes in Amat_->indx.
984 int index;
985
986 Amat_->indx[0] = 0;
987
988 for(int i=0; i<amap_->getNumLocalBlocks(); i++) {
989 int rowBlkSize = Amat_->rpntr[i+1] - Amat_->rpntr[i];
990
991 int colStart = Amat_->bpntr[i];
992 int colEnd = Amat_->bpntr[i+1] - 1;
993
994 for(int j=colStart; j<=colEnd; j++) {
995 if (inUpdate(Amat_->bindx[j], index)) {
996 int colBlkSize = Amat_->rpntr[index+1] - Amat_->rpntr[index];
997
998 Amat_->indx[j+1] = Amat_->indx[j] + rowBlkSize*colBlkSize;
999 }
1000 else { //it's a remoteIndex
1001 if (numProcs == 1) {
1002 char mesg[80];
1003 sprintf(mesg,"calcIndx: blk col index %d not in update set.",
1004 Amat_->bindx[j]);
1005 messageAbort(mesg);
1006 }
1007
1008 index = AZ_find_index(Amat_->bindx[j], remoteInds_,
1009 numRemoteBlocks_);
1010 if (index >= 0) {
1011 Amat_->indx[j+1] = Amat_->indx[j] +
1012 rowBlkSize*remoteBlockSizes_[index];
1013 }
1014 else { //if it wasn't in update or remoteInds, then panic!
1015 messageAbort("calcIndx: block column index not found.");
1016 }
1017 }
1018 } // end for j loop
1019
1020 nnzPerRow_[i] = Amat_->indx[colEnd+1] - Amat_->indx[colStart];
1021 } // end for i loop
1022
1023 localNNZ_ = Amat_->indx[nnzBlks];
1024 }
1025
1026 //==============================================================================
1027 int AztecDVBR_Matrix::getBindxOffset(int blkInd,
1028 int bpntrStart, int bpntrEnd) const {
1029 //
1030 //This function returns the index of blkInd in the bindx array,
1031 //searching positions bindx[bpntrStart..bpntrEnd].
1032 //
1033 //If blkInd is not found, -1 is returned.
1034 //
1035 for(int i=bpntrStart; i<=bpntrEnd; i++) {
1036 int ind = Amat_->bindx[i];
1037 int globalCol = -1;
1038 if (isLoaded()) {
1039 if (ind < N_update_) {
1040 globalCol = amap_->getBlockUpdate()[orderingUpdate_[ind]];
1041 }
1042 else {
1043 globalCol = external_[ind-N_update_];
1044 }
1045 }
1046 else globalCol = ind;
1047
1048 if (globalCol == blkInd) return(i);
1049 }
1050
1051 return(-1);
1052 }
1053
1054 //==============================================================================
1055 void AztecDVBR_Matrix::calcRemoteInds(int*& remoteInds, int& len) {
1056 //
1057 //Form a list of the block column indices that are not in the local
1058 //update set.
1059 //
1060 int nnzBlks = Amat_->bpntr[amap_->getNumLocalBlocks()];
1061 int local;
1062
1063 for(int i=0; i<nnzBlks; i++) {
1064 if (!inUpdate(Amat_->bindx[i], local)) {
1065 insertList(Amat_->bindx[i], remoteInds, len);
1066 }
1067 }
1068 }
1069
1070 //==============================================================================
1071 void AztecDVBR_Matrix::getRemoteBlkSizes(int* remoteBlkSizes,
1072 int* remoteInds,
1073 int len)
1074 {
1075 //
1076 //remoteInds is a sorted list of indices that correspond to rows
1077 //in remote processors' update lists. This function will spread the
1078 //indices to all processors so that they can provide the blk sizes,
1079 //then spread that information back to all processors.
1080 //
1081#ifdef FEI_SER
1082 return;
1083#else
1084 int numProcs = amap_->getProcConfig()[AZ_N_procs];
1085 int thisProc = amap_->getProcConfig()[AZ_node];
1086 MPI_Comm comm = amap_->getCommunicator();
1087
1088 int* lengths = new int[numProcs];
1089 lengths[0] = 0;
1090
1091 //gather up the lengths of the lists that each proc will be sending.
1092 MPI_Allgather(&len, 1, MPI_INT, lengths, 1, MPI_INT, comm);
1093
1094 //now form a list of the offset at which each proc's contribution will
1095 //be placed in the all-gathered list.
1096 int* offsets = new int[numProcs];
1097
1098 offsets[0] = 0;
1099 int totalLength = lengths[0];
1100 for(int i=1; i<numProcs; i++) {
1101 offsets[i] = offsets[i-1] + lengths[i-1];
1102 totalLength += lengths[i];
1103 }
1104
1105 //now we can allocate the list to recv into.
1106 int* recvBuf = new int[totalLength];
1107
1108 //now we're ready to do the gather.
1109 MPI_Allgatherv(remoteInds, len, MPI_INT, recvBuf, lengths, offsets,
1110 MPI_INT, comm);
1111
1112 //now we'll run through the list and put block sizes into a list of
1113 //the same length as the total recvBuf list.
1114 int* blkSizes = new int[totalLength];
1115 int index;
1116
1117 for(int j=0; j<totalLength; j++) {
1118 if (inUpdate(recvBuf[j], index)) {
1119 blkSizes[j] = Amat_->rpntr[index+1]-Amat_->rpntr[index];
1120 }
1121 else blkSizes[j] = 0;
1122 }
1123
1124 //now we'll reduce this info back onto all processors. We'll use MPI_SUM.
1125 //Since the sizes we did NOT supply hold a 0, and each spot in the list
1126 //should only have a nonzero size from 1 processor, the result will be
1127 //that each spot in the result list has the correct value.
1128 int* recvSizes = new int[totalLength];
1129
1130 MPI_Allreduce(blkSizes, recvSizes, totalLength, MPI_INT, MPI_SUM, comm);
1131
1132 //and finally, we just need to run our section of the list of recv'd sizes,
1133 //and transfer them into the remoteBlkSizes list.
1134 int offset = offsets[thisProc];
1135 for(int k=0; k<len; k++) {
1136 remoteBlkSizes[k] = recvSizes[offset + k];
1137 if (recvSizes[offset+k] <= 0)
1138 messageAbort("getRemoteBlkSizes: recvd a size <= 0.");
1139 }
1140
1141 delete [] lengths;
1142 delete [] offsets;
1143 delete [] recvBuf;
1144 delete [] blkSizes;
1145 delete [] recvSizes;
1146#endif
1147 }
1148
1149 //==============================================================================
1150 void AztecDVBR_Matrix::insertList(int item, int*& list, int& len) {
1151 //
1152 //insert 'item' in 'list', if it's not already in there,
1153 //and update the list's length, 'len'.
1154 //
1155 //We want to keep the list ordered, so we'll insert item in
1156 //the list after the biggest existing entry that's smaller, and
1157 //before the smallest existing entry that's bigger.
1158 //
1159
1160 if (len <= 0) {
1161 list = new int[1];
1162 list[0] = item;
1163 len = 1;
1164 return;
1165 }
1166
1167 int index = AZ_find_index(item, list, len);
1168
1169 if (index >= 0) return;
1170
1171 int* newList = new int[len+1];
1172
1173 //bring over the contents of the old list, putting in the new
1174 //one at the appropriate point.
1175 int inserted = 0;
1176 for(int i=0; i<len; i++) {
1177 if (!inserted) {
1178 if (list[i] < item) newList[i] = list[i];
1179 else {
1180 newList[i] = item;
1181 inserted = 1;
1182 }
1183 }
1184 else newList[i] = list[i-1];
1185 }
1186
1187 //now put in the last list entry
1188 if (inserted) newList[len] = list[len-1];
1189 else newList[len] = item;
1190
1191 //delete the old memory and reset the pointer.
1192 if (len > 0) delete [] list;
1193 list = newList;
1194
1195 //update the length.
1196 len++;
1197 }
1198
1199}//namespace fei_trilinos
1200
1201#endif
1202//HAVE_FEI_AZTECOO
1203
AztecDVBR_Matrix(fei::SharedPtr< Aztec_BlockMap > map)
#define FEI_ENDL
#define MPI_Barrier(a)
Definition fei_mpi.h:61
#define MPI_Comm
Definition fei_mpi.h:56
std::ostream & console_out()
int numProcs(MPI_Comm comm)