FEI Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fei_AztecDMSR_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 
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 
77 namespace 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 
#define MPI_Abort(a, b)
Definition: fei_mpi.h:59
#define MPI_Comm
Definition: fei_mpi.h:56
#define MPI_Barrier(a)
Definition: fei_mpi.h:61
#define FEI_ENDL
std::ostream & console_out()
int n
int numProcs(MPI_Comm comm)