FEI Version of the Day
Loading...
Searching...
No Matches
fei_AztecDMSR_Matrix.cpp
1/*
2// @HEADER
3// ************************************************************************
4// FEI: Finite Element Interface to Linear Solvers
5// Copyright (2005) Sandia Corporation.
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8// U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Alan Williams (william@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41*/
42
43
44#include <fei_trilinos_macros.hpp>
45#include <fei_iostream.hpp>
46
47#ifdef HAVE_FEI_AZTECOO
48
49#include <assert.h>
50#include <stdlib.h>
51#include <math.h>
52
53#include <fei_mpi.h>
54
55#ifndef FEI_SER
56
57#define AZTEC_MPI AZTEC_MPI
58#define AZ_MPI AZ_MPI
59#ifndef MPI
60#define MPI MPI
61#endif
62
63#endif
64
65#include <az_aztec.h>
66
67#include <fei_ArrayUtils.hpp>
68
69#include <fei_Aztec_Map.hpp>
70#include <fei_Aztec_LSVector.hpp>
71#include <fei_AztecDMSR_Matrix.hpp>
72
73#define ADMSR_LOCAL_ROW_ALLOC_LEN(localRow) 1+bindx[localRow+1]-bindx[localRow]
74
75#define ADMSR_LOCAL_ROW_LEN(localRow) 1+rowLengths_[localRow]
76
77namespace fei_trilinos {
78
79 //==============================================================================
80 AztecDMSR_Matrix::AztecDMSR_Matrix(fei::SharedPtr<Aztec_Map> map)
81 : isFilled_(false),
82 isAllocated_(false),
83 localOffset_(map->localOffset()),
84 localSize_(map->localSize()),
85 amap_(map),
86 Amat_(NULL),
87 arraysAllocated_(false),
88 val(NULL),
89 bindx(NULL),
90 rowLengths_(NULL),
91 nnzeros_(0),
92 N_update_(map->localSize()),
93 tmp_array_(0),
94 tmp_array_len_(0),
95 dtmp_array_(0),
96 dtmp_array_len_(0),
97 azTransformed_(false)
98 {
99 if (N_update_ > 0) {
100 rowLengths_ = new int[N_update_];
101 for(int i=0; i<N_update_; i++) {
102 rowLengths_[i] = 0;
103 }
104 }
105
106 Amat_ = AZ_matrix_create(N_update_);
107 }
108
109 //==============================================================================
110 AztecDMSR_Matrix::AztecDMSR_Matrix(const AztecDMSR_Matrix& src)
111 : isFilled_(src.isFilled_),
112 isAllocated_(src.isAllocated_),
113 localOffset_(src.localOffset_),
114 localSize_(src.localSize_),
115 amap_(src.amap_),
116 Amat_(NULL),
117 arraysAllocated_(src.arraysAllocated_),
118 val(NULL),
119 bindx(NULL),
120 rowLengths_(NULL),
121 nnzeros_(src.nnzeros_),
122 N_update_(src.N_update_),
123 tmp_array_(0),
124 tmp_array_len_(0),
125 dtmp_array_(0),
126 dtmp_array_len_(0),
127 azTransformed_(src.azTransformed_)
128 {
129 expand_array(tmp_array_, tmp_array_len_, src.tmp_array_len_);
130 expand_array(dtmp_array_, dtmp_array_len_, src.dtmp_array_len_);
131
132 if (N_update_ > 0) {
133 rowLengths_ = new int[N_update_];
134 for(int i=0; i<N_update_; i++) {
135 rowLengths_[i] = src.rowLengths_[i];
136 }
137 }
138
139 Amat_ = AZ_matrix_create(N_update_);
140
141 if (isAllocated_ && nnzeros_ > 0) {
142 val = new double[nnzeros_+1];
143 bindx = new int[nnzeros_+1];
144
145 for(int i=0; i<nnzeros_+1; ++i) {
146 val[i] = src.val[i];
147 bindx[i] = src.bindx[i];
148 }
149
150 if (isFilled_) {
151 AZ_set_MSR(Amat_, bindx, val, amap_->data_org, 0, NULL, AZ_LOCAL);
152 }
153 }
154 }
155
156 //==============================================================================
157 AztecDMSR_Matrix::~AztecDMSR_Matrix()
158 {
159 if (arraysAllocated_) {
160 delete [] val;
161 val = NULL;
162 delete [] bindx;
163 bindx = NULL;
164 arraysAllocated_ = false;
165 }
166
167 if (dtmp_array_) delete [] dtmp_array_; dtmp_array_ = NULL;
168
169 if (N_update_ > 0) {
170 delete [] rowLengths_;
171 rowLengths_ = NULL;
172 N_update_ = 0;
173 }
174
175 if (azTransformed_) {
176 azTransformed_ = false;
177 }
178
179 AZ_matrix_destroy(&Amat_);
180 Amat_ = NULL;
181
182 isFilled_ = false;
183 isAllocated_ = false;
184 arraysAllocated_ = false;
185
186 if (tmp_array_len_ > 0) {
187 delete [] tmp_array_;
188 tmp_array_ = 0;
189 tmp_array_len_ = 0;
190 }
191 }
192
193 //==============================================================================
194 void AztecDMSR_Matrix::expand_array(int*& array, int& arraylen, int newlen)
195 {
196 if (arraylen < newlen) {
197 delete [] array;
198 array = new int[newlen];
199 for(int i=0; i<newlen; ++i) array[i] = -999;
200 arraylen = newlen;
201 }
202 }
203
204 //==============================================================================
205 void AztecDMSR_Matrix::expand_array(double*& array, int& arraylen, int newlen)
206 {
207 if (arraylen < newlen) {
208 delete [] array;
209 array = new double[newlen];
210 for(int i=0; i<newlen; ++i) array[i] = -999.9;
211 arraylen = newlen;
212 }
213 }
214
215 //==============================================================================
216 void AztecDMSR_Matrix::scale(double scalar)
217 {
218 if (scalar == 1.0) return;
219
220 if (val != NULL) {
221 for(int i=0; i<nnzeros_+1; ++i) {
222 val[i] *= scalar;
223 }
224 }
225 }
226
227 //==============================================================================
228 void AztecDMSR_Matrix::matvec(const Aztec_LSVector& x,
229 Aztec_LSVector& y) const
230 {
231 //
232 // This function forms the product y = Ax
233 //
234
235 assert(isFilled());
236
237 int *proc_config = amap_->getProcConfig();
238 // int *idummy = 0, one=1;
239 double *b = (double*)x.startPointer();
240 double *c = (double*)y.startPointer();
241
242 AZ_MSR_matvec_mult(b, c, Amat_, proc_config);
243 //val,idummy,bindx,idummy,idummy,idummy,b,c,one, data_org_);
244
245 return;
246 }
247
248 //==============================================================================
249 void AztecDMSR_Matrix::put(double s)
250 {
251 if(isAllocated()){
252 for(int i=0; i<bindx[N_update_]; i++) val[i] = s;
253 }
254 else {
255 fei::console_out() << "AztecDMSR_Matrix::put - ERROR, can't do put until allocated"
256 << FEI_ENDL;
257 }
258 return;
259 }
260
261 //==============================================================================
262 int AztecDMSR_Matrix::rowLength(int row) const
263 {
264 int localRow;
265
266 int thisRow = row;
267
268 if(amap_->inUpdate(thisRow,localRow)){
269 return(ADMSR_LOCAL_ROW_ALLOC_LEN(localRow));
270 }
271 else {
272 fei::console_out() << "AztecDMSR_Matrix::rowLength: ERROR row " << row
273 << " not in local update set." << FEI_ENDL;
274 abort();
275 return(-1);
276 }
277 }
278
280 int AztecDMSR_Matrix::setDiagEntry(int row, double value)
281 {
282 int thisRow = row;
283 int localRow = -1;
284 if(!amap_->inUpdate(thisRow,localRow)){
285 fei::console_out() << "AztecDMSR_Matrix::setDiagEntry: ERROR - row " << row
286 << " not in local update set." << FEI_ENDL;
287 abort(); return(-1);
288 }
289
290 val[localRow] = value;
291 return(0);
292 }
293
295 double AztecDMSR_Matrix::getDiagEntry(int row) const
296 {
297 int thisRow = row;
298 int localRow = -1;
299 if(!amap_->inUpdate(thisRow,localRow)){
300 fei::console_out() << "AztecDMSR_Matrix::getDiagEntry: ERROR - row " << row
301 << " not in local update set." << FEI_ENDL;
302 abort(); return val[0];
303 }
304
305 return(val[localRow]);
306 }
307
309 int AztecDMSR_Matrix::getOffDiagRowPointers(int row, int*& colIndices,
310 double*& coefs,
311 int& offDiagRowLength)
312 {
313 int thisRow = row;
314 int localRow = -1;
315 if(!amap_->inUpdate(thisRow,localRow)){
316 fei::console_out() << "AztecDMSR_Matrix::getOffDiagRowPointers: ERROR - row " << row
317 << " not in local update set." << FEI_ENDL;
318 abort(); return(-1);
319 }
320
321 offDiagRowLength = rowLengths_[localRow];
322 if (isFilled_) offDiagRowLength = ADMSR_LOCAL_ROW_ALLOC_LEN(localRow)-1;
323 int start = bindx[localRow];
324 colIndices = &(bindx[start]);
325 coefs = &(val[start]);
326
327 return(0);
328 }
329
331 void AztecDMSR_Matrix::getRow(int row,
332 int &length,
333 double *coefs,
334 int *colInd) const
335 {
336 //
337 //Note to myself:
338 //getRow, putRow and sumIntoRow are incredibly ugly, and need re-written.
339 //
340
341 bool foundDiagonal = false;
342 double dtmp;
343 int j, localRow, itmp;
344
345 int thisRow = row;
346
347 if(!amap_->inUpdate(thisRow,localRow)){
348 fei::console_out() << "AztecDMSR_Matrix::getRow: ERROR - row " << row
349 << " not in local update set." << FEI_ENDL;
350 length = 0;
351 return;
352 }
353
354 int start = bindx[localRow];
355 int end = bindx[localRow+1]-1;
356
357 j = 0;
358 for(int i=start; i<=end; i++){
359 coefs[j] = val[i];
360 colInd[j] = amap_->getTransformedEqn(bindx[i]);
361
362 if (colInd[j]==row) {
363 //we're at the diagonal element, so put it in.
364 dtmp = coefs[j];
365 itmp = colInd[j];
366 coefs[j] = val[localRow];
367 colInd[j] = row;
368 j++;
369 coefs[j] = dtmp;
370 colInd[j] = itmp;
371 foundDiagonal = true;
372 }
373 j++;
374 }
375
376 if(!foundDiagonal){ // still need to put in the diagonal
377 coefs[j] = val[localRow];
378 colInd[j] = row;
379 }
380
381 length = j+1;
382 return;
383 }
384
386 int AztecDMSR_Matrix::putRow(int row, int len, const double *coefs,
387 const int *colInd)
388 {
389 //
390 //Note to myself:
391 //getRow, putRow and sumIntoRow are incredibly ugly, and need re-written.
392 //
393
394 int j, localRow, globalColIndex;
395
396 int thisRow = row;
397
398 if (!amap_->inUpdate(thisRow,localRow)){
399 fei::console_out() << "AztecDMSR_Matrix::putRow: ERROR row " << row
400 << " not in local update set." << FEI_ENDL;
401 return(-1);
402 }
403
404 int offDiagRowLen = rowLengths_[localRow];
405 int rowAllocLen = ADMSR_LOCAL_ROW_ALLOC_LEN(localRow);
406 if (isFilled_) offDiagRowLen = rowAllocLen - 1;
407 int offDiagRowAllocLen = rowAllocLen - 1;
408
409 if (len > rowAllocLen) {
410 fei::console_out() << "AztecDMSR_Matrix::putRow. too many coefs, row " << row << FEI_ENDL;
411 return(-1);
412 }
413
414 int jLimit = bindx[localRow+1] - 1;
415 int jStart = bindx[localRow];
416 int* colInds = &(bindx[jStart]);
417 int colStart = *colInds;
418 double* rowCoefs = &(val[jStart]);
419
420 for(int i=0; i<len; i++){
421 if (colInd[i] == row){ //it's on the diagonal
422 val[localRow] = coefs[i];
423 continue;
424 }
425
426 //look along the row until we find the right col index
427 //or an empty spot
428 j=jStart;
429 if (isFilled()){
430 int colIndex = colStart;
431 int col = colInd[i];
432 while (j <= jLimit) {
433 globalColIndex = amap_->getTransformedEqn(colIndex);
434
435 if (globalColIndex == col) break;
436
437 colIndex = bindx[++j];
438 }
439
440 //now put the coefficient in if we haven't gone too far
441 //along the row
442 if (j <= jLimit) {
443 val[j] = coefs[i];
444 }
445 else {
446 fei::console_out() << "AztecDMSR_Matrix::putRow: ERROR didn't "
447 << " find col. index " << colInd[i] << " in "
448 << "row " << row << FEI_ENDL;
449 return(-1);
450 }
451 }
452 else { // !isFilled()
453
454 //first, look for colInd[i] in the row
455 int col = colInd[i];
456 int insertPoint = -1;
457 int index = fei::binarySearch<int>(col, colInds, offDiagRowLen, insertPoint);
458
459 if (index >= offDiagRowAllocLen){ //bad index
460 fei::console_out() << "AztecDMSR_Matrix::putRow, ERROR: "
461 << "row " << row << ", colInd["<<i<<"] " << colInd[i]
462 << ", index = " << index << FEI_ENDL;
463 return(-1);
464 }
465
466 if (index >= 0) {
467 rowCoefs[index] = coefs[i];
468 }
469 else {
470 int tmp = offDiagRowLen;
471 int err = insert(col, insertPoint, colInds,
472 tmp, offDiagRowAllocLen);
473 err += insert(coefs[i], insertPoint, rowCoefs,
474 offDiagRowLen, offDiagRowAllocLen);
475 if (err != 0) {
476 fei::console_out() << "AztecDMSR_Matrix::putRow ERROR: failed to add "
477 << "value for index " << col << " to row " << row << FEI_ENDL;
478 return(-1);
479 }
480 rowLengths_[localRow]++;
481 }
482 }
483 }
484
485 return(0);
486 }
487
489 int AztecDMSR_Matrix::sumIntoRow(int numRows, const int* rows,
490 int numCols, const int *colInd,
491 const double* const* coefs)
492 {
493 if (numRows == 0 || numCols == 0) return(0);
494
495 if (!isFilled_) {
496 int err = 0;
497 for(int i=0; i<numRows; ++i) {
498 err = sumIntoRow(rows[i], numCols, coefs[i], colInd);
499 if (err != 0) return(err);
500 }
501
502 return(0);
503 }
504
505 //Now for the harder (but more important) case where isFilled_ == true.
506
507 //first compute max-row-length:
508 int maxRowLen = 0;
509 for(int i=0; i<numRows; ++i) {
510 int row = rows[i];
511 int localRow;
512 if (!amap_->inUpdate(row, localRow)) {
513 fei::console_out() << "AztecDMSR_Matrix::sumIntoRow: ERROR row " << row
514 << " not in local update set [" << amap_->getUpdate()[0] << " ... "
515 << amap_->getUpdate()[N_update_-1] << "]." << FEI_ENDL;
516 return(-1);
517 }
518
519 int rowlen = bindx[localRow+1]-bindx[localRow];
520 if (maxRowLen < rowlen) maxRowLen = rowlen;
521 }
522
523 if (maxRowLen+2*numCols > tmp_array_len_) {
524 expand_array(tmp_array_, tmp_array_len_, maxRowLen+2*numCols);
525 }
526
527 int* incols = &tmp_array_[maxRowLen];
528 int* indirect = incols+numCols;
529
530 for(int jj=0; jj<numCols; ++jj) {
531 incols[jj] = colInd[jj];
532 indirect[jj] = jj;
533 }
534
535 fei::insertion_sort_with_companions<int>(numCols, incols, indirect);
536
537 int row, localRow;
538
539 for(int i=0; i<numRows; ++i) {
540 row = rows[i];
541 if (!amap_->inUpdate(row, localRow)) {
542 fei::console_out() << "AztecDMSR_Matrix::sumIntoRow: ERROR row " << row
543 << " not in local update set [" << amap_->getUpdate()[0] << " ... "
544 << amap_->getUpdate()[N_update_-1] << "]." << FEI_ENDL;
545 return(-1);
546 }
547
548 int jStart = bindx[localRow];
549 double* rowCoefs = val+jStart;
550 int* rowColInds = bindx+jStart;
551 int rowLen= bindx[localRow+1]-jStart;
552
553 for(int jj=0; jj<rowLen; ++jj) {
554 tmp_array_[jj] = amap_->getTransformedEqn(rowColInds[jj]);
555 }
556
557 const double* coefs_i = coefs[i];
558
559 int inoffset = 0;
560 int incol = incols[inoffset];
561 while (incol == row) {
562 val[localRow] += coefs_i[indirect[inoffset++]];
563 if (inoffset >= numCols) break;
564 incol = incols[inoffset];
565 }
566
567 if (inoffset >= numCols) continue;
568
569 //rowOffset is the offset into the row at which incol appears.
570
571 int rowOffset = fei::binarySearch<int>(incol, tmp_array_, rowLen);
572 if (rowOffset < 0) {
573 fei::console_out() << "AztecDMSR_Matrix::sumIntoRow, ERROR: "
574 << "row " << row << ", col not found: "
575 << incol << FEI_ENDL;
576 return(-1);
577 }
578
579 rowCoefs[rowOffset++] += coefs_i[indirect[inoffset++]];
580
581 //check whether incols has a repeated column-index
582 if (inoffset>0 && incols[inoffset] == incols[inoffset-1]) --rowOffset;
583
584 while(inoffset < numCols) {
585 incol = incols[inoffset];
586
587 if (incol == row) {
588 val[localRow] += coefs_i[indirect[inoffset++]];
589 continue;
590 }
591
592 while(tmp_array_[rowOffset] != incol) {
593 ++rowOffset;
594 if (rowOffset >= rowLen) {
595 fei::console_out() << "AztecDMSR_Matrix::sumIntoRow, ERROR, col "
596 << incol << " not found in row " << row << FEI_ENDL;
597 return(-1);
598 }
599 }
600
601 rowCoefs[rowOffset++] += coefs_i[indirect[inoffset++]];
602 if (inoffset>0 && incols[inoffset] == incols[inoffset-1]) --rowOffset;
603 }
604 }
605
606 return(0);
607 }
608
610 int AztecDMSR_Matrix::sumIntoRow(int row, int len, const double *coefs,
611 const int *colInd)
612 {
613 //
614 //Note to myself:
615 //getRow, putRow and sumIntoRow are incredibly ugly, and need re-written.
616 //
617
618 int localRow, thisRow = row ;
619
620 if (!amap_->inUpdate(thisRow,localRow)) {
621 fei::console_out() << "AztecDMSR_Matrix::sumIntoRow: ERROR row " << row
622 << " not in local update set." << FEI_ENDL;
623 return(-1);
624 }
625
626 int jLimit = bindx[localRow+1] - 1;
627 int jStart = bindx[localRow];
628 int jLen = jLimit-jStart+1;
629 int* colInds = &(bindx[jStart]);
630 double* rowCoefs = &(val[jStart]);
631
632 if (isFilled_) {
633 if (jLen+len > tmp_array_len_) {
634 expand_array(tmp_array_, tmp_array_len_, jLen+len);
635 expand_array(dtmp_array_, dtmp_array_len_, len);
636 }
637 for(int jj=0; jj<jLen; ++jj) {
638 tmp_array_[jj] = amap_->getTransformedEqn(colInds[jj]);
639 }
640
641 int* incols = &tmp_array_[jLen];
642 int doffs = 0;
643 for(int jj=0; jj<len; ++jj) {
644 int col = colInd[jj];
645 if (col == row) {
646 val[localRow] += coefs[jj];
647 }
648 else {
649 incols[doffs] = col;
650 dtmp_array_[doffs++] = coefs[jj];
651 }
652 }
653 fei::insertion_sort_with_companions<double>(doffs, incols, dtmp_array_);
654
655 int ioffset = 0;
656 int offset = fei::binarySearch<int>(incols[ioffset], tmp_array_, jLen);
657 if (offset < 0) {
658 fei::console_out() << "AztecDMSR_Matrix::sumIntoRow, ERROR: "
659 << "row " << row << ", col not found: "
660 << colInd[ioffset] << FEI_ENDL;
661 return(-1);
662 }
663
664 rowCoefs[offset++] += dtmp_array_[ioffset++];
665 if (incols[ioffset] == tmp_array_[offset-1]) --offset;
666
667 while(ioffset < doffs) {
668 int incol = incols[ioffset];
669
670 while(tmp_array_[offset] != incol) {
671 ++offset;
672 if (offset >= jLen) {
673 fei::console_out() << "AztecDMSR_Matrix::sumIntoRow, ERROR, col "
674 << incols[ioffset] << " not found in row " << row << FEI_ENDL;
675 return(-1);
676 }
677 }
678 rowCoefs[offset++] += dtmp_array_[ioffset++];
679 if (incols[ioffset] == tmp_array_[offset-1]) --offset;
680 }
681
682 return(0);
683 }
684
685 //if we get to here, then we know that isFilled_ is false...
686
687 int rowAllocLen = ADMSR_LOCAL_ROW_ALLOC_LEN(localRow);
688 int offDiagRowLen = isFilled_ ? rowAllocLen - 1 : rowLengths_[localRow];
689 int offDiagRowAllocLen = rowAllocLen - 1;
690
691 for(int i=0; i<len; i++){
692 if (colInd[i] == row){ //it's on the diagonal
693 val[localRow] += coefs[i];
694 continue;
695 }
696
697 //find the right col index in the row, or an empty spot
698 int col = colInd[i];
699 int insertPoint = -1;
700 int index = fei::binarySearch<int>(col, colInds, offDiagRowLen, insertPoint);
701
702 if (index >= 0) {
703 rowCoefs[index] += coefs[i];
704 }
705 else {
706 int tmp = offDiagRowLen;
707 int err = insert(col, insertPoint, colInds,
708 tmp, offDiagRowAllocLen);
709 err += insert(coefs[i], insertPoint, rowCoefs,
710 offDiagRowLen, offDiagRowAllocLen);
711 if (err != 0) {
712 fei::console_out() << "AztecDMSR_Matrix::sumIntoRow ERROR: failed to add "
713 << "value for index " << col << " to row " << row << FEI_ENDL;
714 return(-1);
715 }
716 rowLengths_[localRow]++;
717 }
718 }
719
720 return(0);
721 }
722
723 //==============================================================================
724 int AztecDMSR_Matrix::addScaledMatrix(double scalar,
725 const AztecDMSR_Matrix& source)
726 {
727 if (N_update_ != source.N_update_ ||
728 nnzeros_ != source.nnzeros_ ||
729 isFilled_ != source.isFilled_) {
730 fei::console_out() << "AztecDMSR_Matrix::addScaledMatrix ERROR, not compatible"
731 << FEI_ENDL;
732 return(-1);
733 }
734
735 const double* src_val = source.val;
736 int i;
737 for(i=0; i<N_update_; ++i) {
738 val[i] += scalar*src_val[i];
739 }
740
741 //val[N_update_] is not used.
742
743 if (scalar == 1.0) {
744 for(i=N_update_+1; i<nnzeros_+1; ++i) {
745 val[i] += src_val[i];
746 }
747 }
748 else {
749 for(i=N_update_+1; i<nnzeros_+1; ++i) {
750 val[i] += scalar*src_val[i];
751 }
752 }
753
754 return(0);
755 }
756
757 //==============================================================================
758 int AztecDMSR_Matrix::insert(int item, int offset, int* list,
759 int& len, int allocLen)
760 {
761 if (len >= allocLen) return(-1);
762
763 for(int i=len; i>offset; i--) list[i] = list[i-1];
764
765 list[offset] = item;
766 len++;
767
768 return(0);
769 }
770
771 //==============================================================================
772 int AztecDMSR_Matrix::insert(double item, int offset, double* list,
773 int& len, int allocLen)
774 {
775 if (len >= allocLen) return(-1);
776
777 for(int i=len; i>offset; i--) list[i] = list[i-1];
778
779 list[offset] = item;
780 len++;
781
782 return(0);
783 }
784
785 //==============================================================================
786 void AztecDMSR_Matrix::getDiagonal(Aztec_LSVector& diagVector) const {
787
791 // have each processor form its piece of the diagonal
792 double *pdv = (double*)diagVector.startPointer();
793
794 for(int i=0; i<N_update_; i++){
795 pdv[i] = val[i];
796 }
797 }
798
800 void AztecDMSR_Matrix::allocate(int *rowLengths)
801 {
802 //
803 //We assume that 'rowLengths' is of length N_update_.
804 //
805 //rowLengths contains the length of each local row, *NOT* including the
806 //coefficient on the diagonal.
807 //
808 int i;
809
810 //first, count how many non-zeros there are in the local submatrix
811
812 nnzeros_ = 0;
813 for(i=0; i<N_update_; i++){
814 if (rowLengths[i] < 0) {
815 messageAbort("allocate: negative row length");
816 }
817 nnzeros_ += rowLengths[i] + 1;
818 }
819
820 if (bindx != NULL) {
821 delete [] bindx;
822 }
823 bindx = new int[nnzeros_+1];
824
825 if (val != NULL) {
826 delete [] val;
827 }
828 val = new double[nnzeros_+1];
829
830 arraysAllocated_ = true;
831
832 for(i=0; i<nnzeros_+1; i++){
833 val[i] = 0.0;
834 bindx[i] = -1;
835 }
836
837 bindx[0] = N_update_+1;
838
839 for(i=0; i<N_update_; i++){
840 bindx[i+1] = bindx[i] + rowLengths[i];
841 if (bindx[i+1] < 0) {
842 messageAbort("allocate: bindx row length negative.");
843 }
844 }
845
846 //val[N_update_] not used by aztec but we'll initialize it anyway...
847 val[N_update_] = 0.0;
848
849 AZ_set_MSR(Amat_, bindx, val,amap_->data_org, 0, NULL, AZ_LOCAL);
850
851 setAllocated(true);
852 return;
853 }
854
856 void AztecDMSR_Matrix::allocate(int *rowLengths,
857 const int* const* colIndices)
858 {
859 allocate(rowLengths);
860 int col;
861
862 int offset = N_update_+1;
863 for(int i=0; i<N_update_; ++i) {
864 const int* row_cols = colIndices[i];
865 int row_len = rowLengths[i];
866 rowLengths_[i] = row_len-1;
867
868 int prev_col = -999;
869 int coffset = 0;
870 for(int j=0; j<row_len; ++j) {
871 col = row_cols[coffset++];
872
873 if (col == localOffset_+i) {
874 col = row_cols[coffset++];
875 }
876
877 if (col <= prev_col) {
878 messageAbort("allocate: column-indices not sorted.");
879 }
880
881 prev_col = col;
882
883 bindx[offset++] = col;
884 }
885 }
886 }
887
889 double AztecDMSR_Matrix::rowMax(int row) const {
890 int localRow;
891 double max = 0.0;
892
893 if(!amap_->inUpdate(row,localRow)){
894 fei::console_out() << "AztecDMSR_Matrix::rowMax: ERROR row " << row
895 << " not in local update set." << FEI_ENDL;
896 return(-1.0);
897 }
898
899 max = fabs(val[localRow]);
900
901 for(int i=bindx[localRow]; i<bindx[localRow+1]; i++)
902 if(fabs(val[i])>max)max = fabs(val[i]);
903
904 return(max);
905 }
906
908 void AztecDMSR_Matrix::fillComplete() {
909 /*
910 This is where we call the Aztec function AZ_transform, which calculates
911 communication parameters and re-orders the equations for use as a
912 global distributed matrix.
913 */
914 if (isFilled_ || azTransformed_) {
915 isFilled_ = true;
916 azTransformed_ = true;
917 return;
918 }
919
920 int *proc_config = amap_->getProcConfig();
921 int *dummy = 0;
922
923 //before we turn Aztec loose on the matrix, lets do a quick check on the
924 //indices to try to make sure none of them are garbage...
925 int globalSize = amap_->globalSize();
926 for(int i=N_update_+1; i<nnzeros_+1; i++) {
927 if (bindx[i] < 0 || bindx[i] >= globalSize) {
928 fei::console_out() << "AztecDMSR_Matrix: ERROR, bindx["<<i<<"]: " << bindx[i]
929 << ", globalSize: " << globalSize << FEI_ENDL;
930#ifndef FEI_SER
931 MPI_Comm thisComm = amap_->getCommunicator();
932 MPI_Abort(thisComm, -1);
933#endif
934 }
935 }
936
937 AZ_transform(proc_config, &amap_->external, bindx, val,
938 amap_->getUpdate(), &amap_->update_index,
939 &amap_->extern_index, &amap_->data_org, N_update_,
940 dummy, dummy, dummy, &dummy, AZ_MSR_MATRIX);
941
942 //AZ_transform allocates these arrays:
943 // amap_->external
944 // amap_->update_index
945 // amap_->extern_index
946 // amap_->data_org
947 //
948 //On return from AZ_transform, the array update_index contains a mapping
949 //to the local re-ordering of the indices of the update array. Now we will fill
950 //the orderingUpdate array with the reverse of that mapping. i.e., a record
951 //of how to get back to the original ordering of the update indices.
952
953 AZ_set_MSR(Amat_, bindx, val, amap_->data_org, 0, NULL, AZ_LOCAL);
954
955 amap_->orderingUpdate.resize(N_update_);
956 for(int ii=0; ii<N_update_; ii++) {
957 amap_->orderingUpdate[amap_->update_index[ii]] = ii;
958 }
959
960 amap_->az_transformed = true;
961 azTransformed_ = true;
962
963 setFilled(true);
964 return;
965 }
966
967 //==============================================================================
968 void AztecDMSR_Matrix::copyStructure(AztecDMSR_Matrix& source)
969 {
970 //
971 //This function copies the structure (essentially just the bindx and
972 //rowLengths_ arrays) and other relevant variables from the 'source' matrix.
973 //The result is that 'this' matrix is laid out the same as 'source'.
974 //
975 nnzeros_ = source.nnzeros_;
976
977 if (arraysAllocated_) {
978 delete [] val;
979 delete [] bindx;
980 arraysAllocated_ = false;
981 }
982
983 val = new double[nnzeros_+1];
984 bindx = new int[nnzeros_+1];
985
986 int i;
987 for(i=0; i<nnzeros_+1; i++) {
988 val[i] = 0.0;
989 bindx[i] = source.bindx[i];
990 }
991
992 for(i=0; i<N_update_; ++i) rowLengths_[i] = source.rowLengths_[i];
993
994 amap_ = source.amap_;
995
996 AZ_set_MSR(Amat_, bindx, val, amap_->data_org, 0, NULL, AZ_LOCAL);
997
998 isFilled_ = source.isFilled_;
999 azTransformed_ = source.azTransformed_;
1000
1001 arraysAllocated_ = true;
1002 setAllocated(true);
1003 }
1004
1005 //=============================================================================
1006 bool AztecDMSR_Matrix::readFromFile(const char *filename)
1007 {
1008 /*
1009 This function reads the matrix data from a matrix-market format data file,
1010 which looks like:
1011
1012 n
1013 i j val
1014 .
1015 .
1016
1017 Important note: we are going to assume that the index-base of the indices in
1018 the file are 0.
1019 */
1020 int i, j, dummy;
1021 double value;
1022 char line[128];
1023
1024 FILE *mfp = fopen(filename,"r");
1025
1026 if(!mfp){
1027 fei::console_out() << "AztecDMSR_Matrix::readFromFile - couldn't open matrix file."
1028 << FEI_ENDL;
1029 return(false);
1030 }
1031
1032 if (strstr(filename, ".mtx") == NULL) {
1033 fei::console_out() << "AztecDMSR_Matrix::readFromFile: filename doesn't contain "
1034 << "'.mtx'. File should be a MatrixMarket file." << FEI_ENDL;
1035 return(false);
1036 }
1037
1038 do {
1039 fgets(line,128,mfp);
1040 } while(strchr(line,'%'));
1041
1042 while(!feof(mfp)){
1043 do {
1044 fgets(line,128,mfp);
1045 } while(strchr(line,'%'));
1046 if(feof(mfp)){
1047 fclose(mfp);
1048 return(true);
1049 }
1050 sscanf(line,"%d %d %le",&i,&j,&value);
1051
1052 if(amap_->inUpdate(i, dummy)) {
1053 if (putRow(i, 1, &value, &j) != 0) return(false);
1054 }
1055 }
1056
1057 fclose(mfp);
1058 return(true);
1059 }
1060
1062 bool AztecDMSR_Matrix::writeToFile(const char *fileName) const
1063 {
1064 /* Write the matrix into the file "fileName", using the format:
1065 n
1066 i j val
1067 .
1068 .
1069 */
1070
1071 int numProcs = amap_->getProcConfig()[AZ_N_procs];
1072 int thisProc = amap_->getProcConfig()[AZ_node];
1073 int masterRank = 0;
1074
1075
1076 int localNNZ = nnzeros_;
1077 int globalNNZ = localNNZ;
1078#ifndef FEI_SER
1079 MPI_Comm thisComm = amap_->getCommunicator();
1080 MPI_Allreduce(&localNNZ, &globalNNZ, 1, MPI_INT, MPI_SUM, thisComm);
1081#endif
1082
1083 for(int p=0; p<numProcs; p++){
1084
1085 //A barrier inside the loop so each processor waits its turn.
1086#ifndef FEI_SER
1087 MPI_Barrier(thisComm);
1088#endif
1089
1090 if (p == thisProc) {
1091 FILE *file = NULL;
1092
1093 if (masterRank == thisProc) {
1094 //This is the master processor, open a new file.
1095 file = fopen(fileName,"w");
1096
1097 //Write the matrix dimensions n and n (rows==cols) into the file,
1098 //along with the global number-of-nonzeros globalNNZ.
1099
1100 int n = amap_->globalSize();
1101 fprintf(file,"%d %d %d\n",n, n, globalNNZ);
1102 }
1103 else {
1104 //This is a non-master node, open the file for appending to.
1105 file = fopen(fileName,"a");
1106 }
1107
1108 //Now loop over the local portion of the matrix.
1109 for(int i=0; i<N_update_; i++){
1110 int row = localOffset_+i;
1111
1112 int localRow = -1;
1113 if (!amap_->inUpdate(row, localRow)) return(false);
1114
1115 int offDiagRowLen = ADMSR_LOCAL_ROW_LEN(localRow) - 1;
1116 if (isFilled_) offDiagRowLen = ADMSR_LOCAL_ROW_ALLOC_LEN(localRow) - 1;
1117 int* colInds = &(bindx[bindx[localRow]]);
1118 double* coefs = &(val[bindx[localRow]]);
1119
1120 bool wroteDiagonal = false;
1121 for(int j=0; j<offDiagRowLen; j++) {
1122 int col = colInds[j];
1123 int globalCol = col;
1124 if (isFilled()) {
1125 globalCol = amap_->getTransformedEqn(col);
1126 }
1127
1128 if (globalCol >= row && !wroteDiagonal) {
1129 fprintf(file,"%d %d %20.13e\n", row, row,
1130 val[localRow]);
1131 wroteDiagonal = true;
1132 }
1133
1134 fprintf(file,"%d %d %20.13e\n", row, globalCol, coefs[j]);
1135 }
1136 if (!wroteDiagonal) {
1137 fprintf(file,"%d %d %20.13e\n", row, row,
1138 val[localRow]);
1139 wroteDiagonal = true;
1140 }
1141 }
1142 fclose(file);
1143 }
1144 }
1145
1146 return(true);
1147 }
1148
1149 //==============================================================================
1150 void AztecDMSR_Matrix::messageAbort(const char* mesg) {
1151 fei::console_out() << "AztecDMSR_Matrix: ERROR: " << mesg << " Aborting." << FEI_ENDL;
1152 abort();
1153 }
1154
1155}//namespace fei_trilinos
1156
1157#endif
1158//HAVE_FEI_AZTECOO
1159
std::ostream & console_out()
int numProcs(MPI_Comm comm)