FEI Version of the Day
Loading...
Searching...
No Matches
fei_Matrix_Impl.hpp
1/*--------------------------------------------------------------------*/
2/* Copyright 2005 Sandia Corporation. */
3/* Under the terms of Contract DE-AC04-94AL85000, there is a */
4/* non-exclusive license for use of this work by or on behalf */
5/* of the U.S. Government. Export of this program may require */
6/* a license from the United States Government. */
7/*--------------------------------------------------------------------*/
8
9#ifndef _fei_Matrix_Impl_hpp_
10#define _fei_Matrix_Impl_hpp_
11
12#include <fei_fwd.hpp>
13#include <fei_mpi.h>
14#include <fei_defs.h>
15#include "fei_iostream.hpp"
16#include "fei_fstream.hpp"
17#include "fei_sstream.hpp"
18#include <fei_ArrayUtils.hpp>
19#include <fei_MatrixTraits.hpp>
20#include <fei_MatrixTraits_LinProbMgr.hpp>
21#include <fei_MatrixTraits_LinSysCore.hpp>
22#include <fei_MatrixTraits_FEData.hpp>
23#include <fei_MatrixTraits_FillableMat.hpp>
24#include <fei_FillableMat.hpp>
25
26#include <snl_fei_FEMatrixTraits.hpp>
27#include <snl_fei_FEMatrixTraits_FED.hpp>
28#include <snl_fei_BlockMatrixTraits.hpp>
29#include <fei_Matrix.hpp>
30#include <fei_MatrixGraph.hpp>
31#include <fei_Matrix_core.hpp>
32#include <snl_fei_Utils.hpp>
33
34#undef fei_file
35#define fei_file "fei_Matrix_Impl.hpp"
36#include <fei_ErrMacros.hpp>
37
38namespace fei {
39
52 template<typename T>
53 class Matrix_Impl : public fei::Matrix, public fei::Matrix_core {
54 public:
58 int numLocalEqns,
59 bool zeroSharedRows=true);
60
62 virtual ~Matrix_Impl();
63
67 const char* typeName()
68 {
69 if (haveBlockMatrix()) {
71 }
72 else if (haveFEMatrix()) {
74 }
75 else {
77 }
78 }
79
82 int parameters(const fei::ParameterSet& paramset);
83
87 fei::SharedPtr<T> getMatrix() { return( matrix_ ); }
88 const fei::SharedPtr<T> getMatrix() const { return( matrix_ ); }
89
91 {return( Matrix_core::getMatrixGraph() ); }
92
95
98 int getGlobalNumRows() const
99 {
100 if ((int)(globalOffsets().size()) < numProcs()+1) return(-1);
101 return(globalOffsets()[numProcs()]);
102 }
103
106 int getLocalNumRows() const
107 {
108 return(lastLocalOffset() - firstLocalOffset() + 1);
109 }
110
112 int putScalar(double scalar);
113
119 int getRowLength(int row, int& length) const;
120
130 int copyOutRow(int row, int len, double* coefs, int* indices) const;
131
143 int sumIn(int numRows, const int* rows,
144 int numCols, const int* cols,
145 const double* const* values,
146 int format=0);
147
159 int copyIn(int numRows, const int* rows,
160 int numCols, const int* cols,
161 const double* const* values,
162 int format=0);
163
179 int sumInFieldData(int fieldID,
180 int idType,
181 int rowID,
182 int colID,
183 const double* const* data,
184 int format=0);
185
203 int sumInFieldData(int fieldID,
204 int idType,
205 int rowID,
206 int colID,
207 const double* data,
208 int format=0);
209
219 int sumIn(int blockID, int connectivityID,
220 const double* const* values,
221 int format=0);
222
227 int globalAssemble();
228
231 int multiply(fei::Vector* x,
232 fei::Vector* y);
233
234 void setCommSizes();
235
241 int gatherFromOverlap(bool accumulate = true);
242
244 int writeToFile(const char* filename,
245 bool matrixMarketFormat=true);
246
249 int writeToStream(FEI_OSTREAM& ostrm,
250 bool matrixMarketFormat=true);
251
253 {
254 return(haveBlockMatrix());
255 }
256
258 int giveToUnderlyingMatrix(int numRows, const int* rows,
259 int numCols, const int* cols,
260 const double* const* values,
261 bool sumInto,
262 int format);
263
265 int giveToUnderlyingBlockMatrix(int row,
266 int rowDim,
267 int numCols,
268 const int* cols,
269 const int* LDAs,
270 const int* colDims,
271 const double* const* values,
272 bool sumInto);
273
274 void markState();
275
276 bool changedSinceMark();
277
278 double* getBeginPointer()
279 {
280 return fei::MatrixTraits<T>::getBeginPointer(matrix_.get());
281 }
282
283 int getOffset(int row, int col)
284 {
286 fei::SharedPtr<fei::VectorSpace> rowspace = mgraph->getRowSpace();
287 fei::SharedPtr<fei::VectorSpace> colspace = mgraph->getColSpace();
288 int row_index, col_index;
289 int nodeType = 0;//fix this!!! hard-coded 0
290 rowspace->getGlobalIndex(nodeType, row, row_index);
291 colspace->getGlobalIndex(nodeType, col, col_index);
292
293 return fei::MatrixTraits<T>::getOffset(matrix_.get(),row_index,col_index);
294 }
295
296 private:
297 int giveToMatrix(int numRows, const int* rows,
298 int numCols, const int* cols,
299 const double* const* values,
300 bool sumInto,
301 int format);
302
303 int giveToBlockMatrix(int numRows, const int* rows,
304 int numCols, const int* cols,
305 const double* const* values,
306 bool sumInto);
307
308 fei::SharedPtr<T> matrix_;
309 bool globalAssembleCalled_;
310 bool changedSinceMark_;
311 std::string dbgprefix_;
312 };//class Matrix_Impl
313}//namespace fei
314
315//----------------------------------------------------------------------------
316template<typename T>
318{
319 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
320 FEI_OSTREAM& os = *output_stream_;
321 os << dbgprefix_<<"globalAssemble"<<FEI_ENDL;
322 }
323
324 globalAssembleCalled_ = true;
325
326 if (haveBlockMatrix()) {
327 return( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
328 }
329 return( fei::MatrixTraits<T>::globalAssemble(matrix_.get()) );
330}
331
332//----------------------------------------------------------------------------
333template<typename T>
335{
336 Matrix_core::setMatrixGraph(matrixGraph);
337}
338
339//----------------------------------------------------------------------------
340template<typename T>
342{
343if(changedSinceMark_)
344 changedSinceMark_ = false;
345}
346
347//----------------------------------------------------------------------------
348template<typename T>
350{
351 return(changedSinceMark_);
352}
353
354//----------------------------------------------------------------------------
355template<typename T>
356int fei::Matrix_Impl<T>::giveToUnderlyingMatrix(int numRows, const int* rows,
357 int numCols, const int* cols,
358 const double* const* values,
359 bool sumInto,
360 int format)
361{
362 if (format != FEI_DENSE_ROW) {
363 ERReturn(-1);
364 }
365
366 if (output_level_ > fei::BRIEF_LOGS && output_stream_ != 0) {
367 FEI_OSTREAM& os = *output_stream_;
368 for(int i=0; i<numRows; ++i) {
369 os << dbgprefix_<<"giveToUnderlyingMatrix ";
370 for(int j=0; j<numCols; ++j) {
371 os << "("<<rows[i]<<","<<cols[j]<<","<<values[i][j]<<") ";
372 }
373 os << FEI_ENDL;
374 }
375 }
376
377 int err = fei::MatrixTraits<T>::putValuesIn(matrix_.get(), numRows, rows,
378 numCols, cols, values, sumInto);
379 if (err != 0) {
380 return(err);
381 }
382
383 changedSinceMark_ = true;
384 return(0);
385}
386
387//----------------------------------------------------------------------------
388template<typename T>
390 int rowDim,
391 int numCols,
392 const int* cols,
393 const int* LDAs,
394 const int* colDims,
395 const double* const* values,
396 bool sumInto)
397{
398 if (sumInto) {
399 if ( snl_fei::BlockMatrixTraits<T>::sumIn(matrix_.get(),
400 row, rowDim, numCols, cols,
401 LDAs, colDims, values) != 0) {
402 ERReturn(-1);
403 }
404 }
405 else {
406 if ( snl_fei::BlockMatrixTraits<T>::copyIn(matrix_.get(),
407 row, rowDim, numCols, cols,
408 LDAs, colDims, values) != 0) {
409 ERReturn(-1);
410 }
411 }
412
413 changedSinceMark_ = true;
414
415 return(0);
416}
417
418#include <fei_macros.hpp>
419#include <fei_chk_mpi.hpp>
420
421#include <fei_TemplateUtils.hpp>
422
423#include <fei_Pattern.hpp>
424#include <fei_VectorSpace.hpp>
425
426#include <fei_ConnectivityBlock.hpp>
427
428#include <snl_fei_PointBlockMap.hpp>
429
430#include <fei_Record.hpp>
431
432#include <snl_fei_Utils.hpp>
433
434//----------------------------------------------------------------------------
435template<typename T>
438 int numLocalEqns,
439 bool zeroSharedRows)
440 : Matrix_core(matrixGraph, numLocalEqns),
441 matrix_(matrix),
442 globalAssembleCalled_(false),
443 changedSinceMark_(true),
444 dbgprefix_("MatImpl: ")
445{
446 if (strcmp(snl_fei::FEMatrixTraits<T>::typeName(), "unsupported")) {
447 setFEMatrix(true);
448 }
449 else {
450 setFEMatrix(false);
451 }
452
453 if (strcmp(snl_fei::BlockMatrixTraits<T>::typeName(), "unsupported")) {
454 setBlockMatrix(true);
455 }
456 else {
457 setBlockMatrix(false);
458 }
459
460 if (zeroSharedRows && matrixGraph->getGlobalNumSlaveConstraints() == 0) {
461 std::vector<double> zeros;
462 fei::SharedPtr<fei::SparseRowGraph> srg = matrixGraph->getRemotelyOwnedGraphRows();
463 if (srg.get() != NULL) {
464 for(size_t row=0; row<srg->rowNumbers.size(); ++row) {
465 int rowLength = srg->rowOffsets[row+1] - srg->rowOffsets[row];
466 if (rowLength == 0) continue;
467 zeros.resize(rowLength, 0.0);
468 const double* zerosPtr = &zeros[0];
469 const int* cols = &srg->packedColumnIndices[srg->rowOffsets[row]];
470 sumIn(1, &srg->rowNumbers[row], rowLength, cols, &zerosPtr);
471 }
472 setCommSizes();
473 std::map<int,FillableMat*>& remote = getRemotelyOwnedMatrices();
474 for(std::map<int,FillableMat*>::iterator iter=remote.begin(), end=remote.end(); iter!=end; ++iter)
475 {
476 iter->second->clear();
477 }
478 }
479 }
480}
481
482//----------------------------------------------------------------------------
483template<typename T>
487
488//----------------------------------------------------------------------------
489template<typename T>
491{
492 Matrix_core::parameters(paramset);
493 return 0;
494}
495
496//----------------------------------------------------------------------------
497template<typename T>
498int fei::Matrix_Impl<T>::getRowLength(int row, int& length) const
499{
500 if (haveBlockMatrix()) {
501 return( snl_fei::BlockMatrixTraits<T>::getPointRowLength(matrix_.get(), row, length) );
502 }
503 else {
504 int code = fei::MatrixTraits<T>::getRowLength(matrix_.get(), row, length);
505 if (code < 0) {
506 //matrix row probably not local. maybe it's shared.
507 int proc = getOwnerProc(row);
508 const FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
509 if (remote_mat->hasRow(row)) {
510 const CSVec* row_entries = remote_mat->getRow(row);
511 length = row_entries->size();
512 }
513 else {
514 length = 0;
515 }
516 }
517 }
518 return 0;
519}
520
521//----------------------------------------------------------------------------
522template<typename T>
524{
525 if (output_level_ > fei::BRIEF_LOGS && output_stream_ != 0) {
526 FEI_OSTREAM& os = *output_stream_;
527 os << dbgprefix_<<"putScalar("<<scalar<<")"<<FEI_ENDL;
528 }
529
530 if (haveFEMatrix()) {
531 if (scalar != 0.0) return(-1);
532 CHK_ERR( snl_fei::FEMatrixTraits<T>::reset(matrix_.get()) );
533 }
534 else if (haveBlockMatrix()) {
535 if (globalAssembleCalled_ == true) {
536 CHK_ERR( snl_fei::BlockMatrixTraits<T>::putScalar(matrix_.get(), scalar) );
537 }
538 }
539 else {
540 CHK_ERR( fei::MatrixTraits<T>::setValues(matrix_.get(), scalar) );
541 }
542
543 putScalar_remotelyOwned(scalar);
544
545 changedSinceMark_ = true;
546
547 return(0);
548}
549
550//----------------------------------------------------------------------------
551template<typename T>
553 double* coefs, int* indices) const
554{
555 if (len < 1) {
556 return 0;
557 }
558
559 if (haveBlockMatrix()) {
560 int dummy;
561 return( snl_fei::BlockMatrixTraits<T>::copyOutPointRow(matrix_.get(), firstLocalOffset(),
562 row, len,
563 coefs, indices, dummy));
564 }
565 else {
566 int code = fei::MatrixTraits<T>::copyOutRow(matrix_.get(), row, len,
567 coefs, indices);
568 if (code < 0) {
569 int proc = getOwnerProc(row);
570 const FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
571 if (remote_mat->hasRow(row)) {
572 const CSVec* row_entries = remote_mat->getRow(row);
573 const std::vector<int>& row_indices = row_entries->indices();
574 const std::vector<double>& row_coefs = row_entries->coefs();
575 for(size_t i=0; i<row_indices.size(); ++i) {
576 indices[i] = row_indices[i];
577 coefs[i] = row_coefs[i];
578 }
579 }
580 }
581 }
582 return 0;
583}
584
585//----------------------------------------------------------------------------
586template<typename T>
587int fei::Matrix_Impl<T>::sumIn(int numRows, const int* rows,
588 int numCols, const int* cols,
589 const double* const* values,
590 int format)
591{
592 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
593 FEI_OSTREAM& os = *output_stream_;
594 os << dbgprefix_<<"sumIn"<<FEI_ENDL;
595 if (output_level_ >= fei::FULL_LOGS) {
596 for(int i=0; i<numRows; ++i) {
597 for(int j=0; j<numCols; ++j) {
598 os << "("<<rows[i]<<","<<cols[j]<<","<<values[i][j]<<") ";
599 }
600 os << FEI_ENDL;
601 }
602 }
603 }
604
605 return( giveToMatrix( numRows, rows, numCols, cols, values, true, format) );
606}
607
608//----------------------------------------------------------------------------
609template<typename T>
610int fei::Matrix_Impl<T>::copyIn(int numRows, const int* rows,
611 int numCols, const int* cols,
612 const double* const* values,
613 int format)
614{
615 if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
616 FEI_OSTREAM& os = *output_stream_;
617 os << "copyIn"<<FEI_ENDL;
618 if (output_level_ >= fei::FULL_LOGS) {
619 for(int i=0; i<numRows; ++i) {
620 for(int j=0; j<numCols; ++j) {
621 os << "("<<rows[i]<<","<<cols[j]<<","<<values[i][j]<<") ";
622 }
623 os << FEI_ENDL;
624 }
625 }
626 }
627 return( giveToMatrix( numRows, rows, numCols, cols, values, false, format) );
628}
629
630//----------------------------------------------------------------------------
631template<typename T>
633 int idType,
634 int rowID,
635 int colID,
636 const double* const* data,
637 int format)
638{
639 fei::SharedPtr<fei::VectorSpace> vspace = vecSpace();
640
641 int fieldSize = vspace->getFieldSize(fieldID);
642
643 if (fieldSize <= 0) ERReturn(-1);
644
645 work_indices_.resize(fieldSize*2);
646 int* indicesPtr = &work_indices_[0];
647 int i;
648
649 CHK_ERR( vspace->getGlobalIndices(1, &rowID, idType, fieldID, indicesPtr));
650 for(i=1; i<fieldSize; ++i) {
651 indicesPtr[i] = indicesPtr[i-1]+1;
652 }
653
654 CHK_ERR( vspace->getGlobalIndices(1, &colID, idType, fieldID,
655 &(indicesPtr[fieldSize])) );
656 for(i=fieldSize+1; i<2*fieldSize; ++i) {
657 indicesPtr[i] = indicesPtr[i-1]+1;
658 }
659
660 CHK_ERR( sumIn(fieldSize, indicesPtr, fieldSize, &(indicesPtr[fieldSize]),
661 data, format) );
662
663 return(0);
664}
665
666//----------------------------------------------------------------------------
667template<typename T>
669 int idType,
670 int rowID,
671 int colID,
672 const double* data,
673 int format)
674{
675 fei::SharedPtr<fei::VectorSpace> vspace = vecSpace();
676
677 int fieldSize = vspace->getFieldSize(fieldID);
678
679 if (fieldSize <= 0) ERReturn(-1);
680
681 work_data2D_.resize(fieldSize);
682
683 const double** data2DPtr = &work_data2D_[0];
684 for(int i=0; i<fieldSize; ++i) {
685 data2DPtr[i] = &(data[i*fieldSize]);
686 }
687
688 CHK_ERR( sumInFieldData(fieldID, idType, rowID, colID, data2DPtr, format) );
689
690 return(0);
691}
692
693//----------------------------------------------------------------------------
694template<typename T>
695int fei::Matrix_Impl<T>::sumIn(int blockID, int connectivityID,
696 const double* const* values, int format)
697{
698 if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
699 FEI_OSTREAM& os = *output_stream_;
700 os << dbgprefix_<<"sumIn blkID=" << blockID
701 << ", connID=" << connectivityID << FEI_ENDL;
702 }
703
704 fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
705
706 if (mgraph.get() == NULL) ERReturn(-1);
707
708 const fei::ConnectivityBlock* cblock = mgraph->getConnectivityBlock(blockID);
709 if (cblock==NULL) {
710 FEI_OSTRINGSTREAM osstr;
711 osstr << "fei::Matrix_Impl::sumIn ERROR, unable to "
712 << "look up connectivity-block with ID "<<blockID;
713 throw std::runtime_error(osstr.str());
714 }
715
716 bool symmetric = cblock->isSymmetric();
717 const fei::Pattern* pattern = cblock->getRowPattern();
718 const fei::Pattern* colpattern = symmetric ? NULL : cblock->getColPattern();
719 const int* rowConn = cblock->getRowConnectivity(connectivityID);
720
721 fei::SharedPtr<fei::VectorSpace> rspace = mgraph->getRowSpace();
722 rspace->getGlobalIndicesL(pattern, rowConn, work_indices2_);
723
724 int numRowIndices = work_indices2_.size();
725 int* rowIndices = &work_indices2_[0];
726
727 if (haveFEMatrix() || haveBlockMatrix()) {
728 FieldDofMap<int>& fdofmap = rspace->getFieldDofMap();
729 const std::map<int,int>& connIDs = cblock->getConnectivityIDs();
730 std::map<int,int>::const_iterator
731 iter = connIDs.find(connectivityID);
732 if (iter == connIDs.end()) ERReturn(-1);
733 int connOffset = iter->second;
734 int numIDs = pattern->getNumIDs();
735 const int* indicesPerID = pattern->getNumIndicesPerID();
736 const int* fieldsPerID = pattern->getNumFieldsPerID();
737 const int* fieldIDs = pattern->getFieldIDs();
738
739 int numDofs = 0;
740 for(int ii=0; ii<numIDs; ++ii) {
741 numDofs += indicesPerID[ii];
742 }
743
744 const int* numIndicesPerID = pattern->getNumIndicesPerID();
745 work_indices_.resize(numIDs+numDofs);
746 int i, *nodeNumbers = &work_indices_[0];
747 int* dof_ids = nodeNumbers+numIDs;
748
749 int nodeType = 0;
750 snl_fei::RecordCollection* records = NULL;
751 rspace->getRecordCollection(nodeType, records);
752 int foffset = 0;
753 int doffset = 0;
754 for(i=0; i<numIDs; ++i) {
755 fei::Record<int>* rec = records->getRecordWithLocalID(rowConn[i]);
756 nodeNumbers[i] = rec->getNumber();
757 for(int ii=0; ii<fieldsPerID[i]; ++ii) {
758 int fieldSize = rspace->getFieldSize(fieldIDs[foffset]);
759 int dof_id = fdofmap.get_dof_id(fieldIDs[foffset++], 0);
760 for(int j=0; j<fieldSize; ++j) {
761 dof_ids[doffset++] = dof_id + j;
762 }
763 }
764 }
765
766 if (haveFEMatrix()) {
767 CHK_ERR( snl_fei::FEMatrixTraits<T>::sumInElemMatrix(matrix_.get(), blockID, connOffset,
768 numIDs, nodeNumbers,
769 numIndicesPerID, dof_ids, values) );
770 changedSinceMark_ = true;
771 }
772
773 if (haveBlockMatrix()) {
774 if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL &&
775 format != FEI_BLOCK_DIAGONAL_ROW) {
776 fei::console_out() << "fei::Matrix_Impl::sumIn ERROR, for block-matrix, format must"
777 << " be FEI_DENSE_ROW or FEI_DENSE_COL."<<FEI_ENDL;
778 ERReturn(-1);
779 }
780
781 int numPtIndices = pattern->getNumIndices();
782 int* ptIndices = &work_indices2_[0];
783
784 int numPtColIndices = symmetric ? numPtIndices : colpattern->getNumIndices();
785
786 int len = numPtIndices*numPtColIndices;
787
788 if (format == FEI_BLOCK_DIAGONAL_ROW) {
789 len = 0;
790 for(i=0; i<numIDs; ++i) {
791 len += numIndicesPerID[i]*numIndicesPerID[i];
792 }
793 }
794
795 double* ccvalues = new double[len];
796 //Ouch! Data copy! Remember to optimize this later...
797 if (format == FEI_BLOCK_DIAGONAL_ROW) {
798 snl_fei::copy2DBlockDiagToColumnContig(numIDs, numIndicesPerID,
799 values, format, ccvalues);
800 }
801 else {
802 snl_fei::copy2DToColumnContig(numPtIndices, numPtColIndices, values, format,
803 ccvalues);
804 }
805
806 int ptRowOffset = 0;
807 for(i=0; i<numIDs; ++i) {
808
809 fei::Record<int>* record = records->getRecordWithLocalID(rowConn[i]);
810 if (record->getOwnerProc() == localProc()) {
811
812 int numColIDs = numIDs;
813 int* colNodeNums = nodeNumbers;
814 const int* colDims = numIndicesPerID;
815 int LDA = numPtColIndices;
816 if (format == FEI_BLOCK_DIAGONAL_ROW) {
817 numColIDs = 1;
818 colNodeNums = &(nodeNumbers[i]);
819 colDims = &(numIndicesPerID[i]);
820 LDA = numIndicesPerID[i];
821 }
822
823 CHK_ERR( snl_fei::BlockMatrixTraits<T>::sumIn(matrix_.get(),
824 nodeNumbers[i],
825 numIndicesPerID[i],
826 numColIDs, colNodeNums,
827 colDims, LDA,
828 &(ccvalues[ptRowOffset])) );
829 changedSinceMark_ = true;
830 }
831 else {
832 if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
833 FEI_OSTREAM& os = *output_stream_;
834 for(int ii=0; ii<numIndicesPerID[i]; ++ii) {
835 os << "# remote pt-row " << ptIndices[ptRowOffset]+ii <<" ";
836 for(int jj=0; jj<numPtIndices; ++jj) {
837 os << "("<<ptIndices[jj]<<","<<values[ptRowOffset+ii][jj]<<") ";
838 }
839 os << FEI_ENDL;
840 }
841 }
842
843 for(int ii=0; ii<numIndicesPerID[i]; ++ii) {
844 int p=eqnComm_->getOwnerProc(ptIndices[ptRowOffset]+ii);
845 FillableMat* remote_mat = getRemotelyOwnedMatrix(p);
846 remote_mat->sumInRow(ptIndices[ptRowOffset]+ii,
847 ptIndices,
848 values[ptRowOffset+ii],
849 numPtIndices);
850 changedSinceMark_ = true;
851 }
852 }
853
854 ptRowOffset += numIndicesPerID[i];
855 }
856
857 delete [] ccvalues;
858 }
859
860 return(0);
861 }
862
863 int numColIndices = symmetric ? numRowIndices : colpattern->getNumIndices();
864 int* colIndices = rowIndices;
865 const int* colConn = NULL;
866
867 if (!symmetric) {
868 colConn = cblock->getColConnectivity(connectivityID);
869 mgraph->getColSpace()->getGlobalIndicesL(colpattern,
870 colConn, work_indices_);
871 colIndices = &work_indices_[0];
872 }
873
874 if (symmetric) {
875 if (format == FEI_DENSE_ROW || format == FEI_DENSE_COL) {
876 CHK_ERR( sumIn(numRowIndices, rowIndices, numRowIndices, rowIndices,
877 values, format) );
878 }
879 else if (format == FEI_BLOCK_DIAGONAL_ROW) {
880 int totalnumfields = pattern->getTotalNumFields();
881 const int* fieldIDs = pattern->getFieldIDs();
882 fei::SharedPtr<fei::VectorSpace> rowspace = mgraph->getRowSpace();
883 fei::VectorSpace* rowspaceptr = rowspace.get();
884 int ioffset = 0;
885 for(int i=0; i<totalnumfields; ++i) {
886 int fieldsize = rowspaceptr->getFieldSize(fieldIDs[i]);
887 if (ioffset+fieldsize > numRowIndices) {
888 FEI_OSTRINGSTREAM osstr;
889 osstr<<"snl_fei::sumIn, format=FEI_BLOCK_DIAGONAL_ROW, block-sizes"
890 << " not consistent with total num-indices."<<FEI_ENDL;
891 throw std::runtime_error(osstr.str());
892 }
893
894 CHK_ERR( sumIn(fieldsize, &(rowIndices[ioffset]),
895 fieldsize, &(rowIndices[ioffset]),
896 &(values[ioffset]), FEI_DENSE_ROW) );
897 ioffset += fieldsize;
898 changedSinceMark_ = true;
899 }
900 }
901 else {
902 FEI_OSTRINGSTREAM osstr;
903 osstr << "fei::Matrix_Impl::sumIn, format="<<format<<" not supported."
904 << FEI_ENDL;
905 throw std::runtime_error(osstr.str());
906 }
907 }
908 else {
909 if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
910 FEI_OSTRINGSTREAM osstr;
911 osstr << "fei::Matrix_Impl::sumIn, format="<<format<<" not valid with"
912 << " un-symmetric matrix contributions."<<FEI_ENDL;
913 throw std::runtime_error(osstr.str());
914 }
915
916 CHK_ERR( sumIn(numRowIndices, rowIndices, numColIndices, colIndices,
917 values, format) );
918 changedSinceMark_ = true;
919 }
920
921 return(0);
922}
923
924//----------------------------------------------------------------------------
925template<typename T>
927 fei::Vector* y)
928{
929 return( fei::MatrixTraits<T>::matvec(matrix_.get(), x, y) );
930}
931
932//----------------------------------------------------------------------------
933template<typename T>
935{
936 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
937 (*output_stream_) << dbgprefix_<<"setCommSizes"<<FEI_ENDL;
938 }
939
940 Matrix_core::setCommSizes();
941}
942
943//----------------------------------------------------------------------------
944template<typename T>
946{
947 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
948 (*output_stream_) << dbgprefix_<<"gatherFromOverlap"<<FEI_ENDL;
949 }
950
951 CHK_ERR( Matrix_core::gatherFromOverlap(accumulate) );
952
953 return(0);
954}
955
956//----------------------------------------------------------------------------
957template<typename T>
958int fei::Matrix_Impl<T>::giveToMatrix(int numRows, const int* rows,
959 int numCols, const int* cols,
960 const double* const* values,
961 bool sumInto,
962 int format)
963{
964 if (numRows == 0 || numCols == 0) {
965 return(0);
966 }
967
968 if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
969 return(-1);
970 }
971
972 const double** myvalues = const_cast<const double**>(values);
973 if (format != FEI_DENSE_ROW) {
974 copyTransposeToWorkArrays(numRows, numCols, values,
975 work_data1D_, work_data2D_);
976 myvalues = &work_data2D_[0];
977 }
978
979 if ((int)work_ints_.size() < numRows) {
980 work_ints_.resize(numRows);
981 }
982
983 if (haveBlockMatrix()) {
984 return( giveToBlockMatrix(numRows, rows, numCols, cols,
985 myvalues, sumInto) );
986 }
987
988 int i;
989 int numRemote = 0;
990 int* workIntPtr = &work_ints_[0];
991 for(i=0; i<numRows; ++i) {
992 int row = rows[i];
993 if (row < firstLocalOffset() || row > lastLocalOffset()) {
994 ++numRemote;
995 workIntPtr[i] = 1;
996 }
997 else {
998 workIntPtr[i] = 0;
999 }
1000 }
1001
1002 if (numRemote < 1) {
1003 int err = giveToUnderlyingMatrix(numRows, rows, numCols, cols, myvalues,
1004 sumInto, FEI_DENSE_ROW);
1005 changedSinceMark_ = true;
1006 if (err != 0) {
1007 FEI_OSTRINGSTREAM osstr;
1008 osstr << "fei::Matrix_Impl::giveToMatrix ERROR: err="<<err
1009 << " returned from giveToUnderlyingMatrix.";
1010 throw std::runtime_error(osstr.str());
1011 }
1012 return(0);
1013 }
1014
1015 for(i=0; i<numRows; ++i) {
1016 int row = rows[i];
1017 const double*const rowvalues = myvalues[i];
1018
1019 if (workIntPtr[i] > 0) {
1020 int proc = eqnComm_->getOwnerProc(row);
1021 FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
1022
1023 if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
1024 FEI_OSTREAM& os = *output_stream_;
1025 os << dbgprefix_<<" remote_mat["<<proc<<"]: ";
1026 for(int jj=0; jj<numCols; ++jj) {
1027 os << "("<<row<<","<<cols[jj]<<","<<rowvalues[jj]<<") ";
1028 }
1029 os << FEI_ENDL;
1030 }
1031
1032 if (sumInto) {
1033 remote_mat->sumInRow(row, cols, rowvalues, numCols);
1034 }
1035 else {
1036 remote_mat->putRow(row, cols, rowvalues, numCols);
1037 }
1038 }
1039 else {
1040 CHK_ERR( giveToUnderlyingMatrix(1, &row, numCols, cols, &rowvalues,
1041 sumInto, 0) );
1042 }
1043 }
1044
1045 changedSinceMark_ = true;
1046
1047 return(0);
1048}
1049
1050//----------------------------------------------------------------------------
1051template<typename T>
1052int fei::Matrix_Impl<T>::giveToBlockMatrix(int numRows, const int* rows,
1053 int numCols, const int* cols,
1054 const double* const* values,
1055 bool sumInto)
1056{
1057 if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
1058 FEI_OSTREAM& os = *output_stream_;
1059 os << "# giveToBlockMatrix numRows: " << numRows
1060 << ", numCols: " << numCols << FEI_ENDL;
1061 }
1062
1063 if (numRows == 0 || numCols == 0) {
1064 return(0);
1065 }
1066
1067 snl_fei::PointBlockMap* pointBlockMap = vecSpace()->getPointBlockMap();
1068
1069 if (sumInto && numProcs() == 1) {
1070 std::vector<int> temp((numRows+numCols)*2);
1071 int* tempPtr = &temp[0];
1072 int* blkRows = tempPtr;
1073 int* blkRowOffsets = tempPtr+numRows;
1074 int* blkCols = blkRowOffsets+numRows;
1075 int* blkColOffsets = blkCols+numCols;
1076
1077 CHK_ERR( convertPtToBlk(numRows, rows, numCols, cols,
1078 blkRows, blkRowOffsets,
1079 blkCols, blkColOffsets) );
1080
1081 std::vector<int> blockRows, blockRowSizes;
1082 std::vector<int> blockCols, blockColSizes;
1083 for(int i=0; i<numRows; ++i) fei::sortedListInsert(blkRows[i], blockRows);
1084 for(int i=0; i<numCols; ++i) fei::sortedListInsert(blkCols[i], blockCols);
1085
1086 int rowSizeTotal = 0, colSizeTotal = 0;
1087
1088 for(size_t i=0; i<blockRows.size(); ++i) {
1089 int size = pointBlockMap->getBlkEqnSize(blockRows[i]);
1090 blockRowSizes.push_back(size);
1091 rowSizeTotal += size;
1092 }
1093 for(size_t i=0; i<blockCols.size(); ++i) {
1094 int size = pointBlockMap->getBlkEqnSize(blockCols[i]);
1095 blockColSizes.push_back(size);
1096 colSizeTotal += size;
1097 }
1098 std::vector<double> coefs_1d(rowSizeTotal*colSizeTotal, 0.0);
1099 double* coefs1dPtr = &coefs_1d[0];
1100 std::vector<double*> coefs_2d(blockRows.size()*blockCols.size());
1101 double** coefs2dPtr = &coefs_2d[0];
1102
1103 int blkCounter = 0;
1104 int offset = 0;
1105 for(size_t i=0; i<blockRows.size(); ++i) {
1106 for(size_t j=0; j<blockCols.size(); ++j) {
1107 coefs2dPtr[blkCounter++] = &(coefs1dPtr[offset]);
1108 offset += blockRowSizes[i]*blockColSizes[j];
1109 }
1110 }
1111
1112 for(int i=0; i<numRows; ++i) {
1113 int rowind = fei::binarySearch(blkRows[i], blockRows);
1114 int rowsize = blockRowSizes[rowind];
1115
1116 for(int j=0; j<numCols; ++j) {
1117 int colind = fei::binarySearch(blkCols[j], blockCols);
1118 int pos = blkColOffsets[j]*rowsize + blkRowOffsets[i];
1119
1120 coefs2dPtr[rowind*blockCols.size()+colind][pos] += values[i][j];
1121 }
1122 }
1123
1124 for(size_t i=0; i<blockRows.size(); ++i) {
1125 CHK_ERR( giveToUnderlyingBlockMatrix(blockRows[i],
1126 blockRowSizes[i],
1127 blockCols.size(),
1128 &blockCols[0],
1129 &blockColSizes[0],
1130 &blockColSizes[0],
1131 &coefs2dPtr[i*blockCols.size()],
1132 true) );
1133 }
1134 changedSinceMark_ = true;
1135
1136 return(0);
1137 }
1138
1139 int maxBlkEqnSize = pointBlockMap->getMaxBlkEqnSize();
1140 int coefBlkLen = maxBlkEqnSize*maxBlkEqnSize*2;
1141
1142 for(int i=0; i<numRows; ++i) {
1143 int row = rows[i];
1144
1145 if (row < firstLocalOffset() || row > lastLocalOffset()) {
1146 int proc = eqnComm_->getOwnerProc(row);
1147 FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
1148 if (sumInto) {
1149 remote_mat->sumInRow(row, cols, values[i], numCols);
1150 }
1151 else {
1152 remote_mat->putRow(row, cols, values[i], numCols);
1153 }
1154 changedSinceMark_ = true;
1155 continue;
1156 }
1157
1158 int blockRow = pointBlockMap->eqnToBlkEqn(row);
1159 int blockRowSize = pointBlockMap->getBlkEqnSize(blockRow);
1160 int blockRowOffset = pointBlockMap->getBlkEqnOffset(blockRow, row);
1161
1162 int blockRowLength = 0;
1163 CHK_ERR( snl_fei::BlockMatrixTraits<T>::getRowLength(matrix_.get(), blockRow,
1164 blockRowLength) );
1165
1166 std::vector<int> blkCols(blockRowLength);
1167 int* blkCols_ptr = &blkCols[0];
1168 std::vector<int> blkColDims(blockRowLength);
1169 int* blkColDims_ptr = &blkColDims[0];
1170 std::vector<double> coefs_1D(blockRowLength*coefBlkLen);
1171 double* coefs_1D_ptr = &coefs_1D[0];
1172 int coefsLen = coefs_1D.size();
1173 std::vector<double*> coefs_2D(blockRowLength);
1174 double** coefs_2D_ptr = &coefs_2D[0];
1175
1176 std::vector<int> LDAs(blockRowLength);
1177 std::fill(LDAs.begin(), LDAs.end(), blockRowSize);
1178 std::fill(coefs_1D.begin(), coefs_1D.end(), 0.0);
1179
1180 int checkRowLen = 0;
1181 CHK_ERR( snl_fei::BlockMatrixTraits<T>::copyOutRow(matrix_.get(),
1182 blockRow, blockRowLength,
1183 blockRowSize,
1184 blkCols_ptr,
1185 blkColDims_ptr,
1186 coefs_1D_ptr,
1187 coefsLen,
1188 checkRowLen) );
1189 int coefs_1D_offset = 0;
1190 for(int j=0; j<checkRowLen; ++j) {
1191 coefs_2D_ptr[j] = &(coefs_1D_ptr[coefs_1D_offset]);
1192 coefs_1D_offset += blockRowSize*blkColDims_ptr[j];
1193 }
1194
1195 for(int j=0; j<numCols; ++j) {
1196 int blockCol, blkOffset;
1197 CHK_ERR( pointBlockMap->getPtEqnInfo(cols[j], blockCol, blkOffset) );
1198
1199 for(int jj=0; jj<blockRowLength; ++jj) {
1200
1201 if (blockCol == blkCols_ptr[jj]) {
1202 if (sumInto) {
1203 coefs_2D_ptr[jj][blkOffset*blockRowSize+blockRowOffset] += values[i][j];
1204 }
1205 else {
1206 coefs_2D_ptr[jj][blkOffset*blockRowSize+blockRowOffset] = values[i][j];
1207 }
1208
1209 break;
1210 }
1211 }
1212 }
1213
1214 //Now put the block-row back into the matrix
1215 CHK_ERR( giveToUnderlyingBlockMatrix(blockRow, blockRowSize,
1216 blockRowLength, blkCols_ptr,
1217 &LDAs[0],
1218 blkColDims_ptr,
1219 coefs_2D_ptr,
1220 false) );
1221
1222 changedSinceMark_ = true;
1223 }
1224
1225 return(0);
1226}
1227
1228//----------------------------------------------------------------------------
1229template<typename T>
1230int fei::Matrix_Impl<T>::writeToFile(const char* filename,
1231 bool matrixMarketFormat)
1232{
1233 fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
1234 if (mgraph.get() == NULL) {
1235 return(-1);
1236 }
1237
1238 if (haveFEMatrix()) {
1239 return(0);//temporary no-op...
1240 }
1241
1242 if (haveBlockMatrix()) {
1243 CHK_ERR( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
1244 }
1245
1246 static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
1247
1248 fei::SharedPtr<fei::VectorSpace> vspace = mgraph->getRowSpace();
1249
1250 int globalNNZ = 0;
1251 int globalNumRows = vspace->getGlobalNumIndices();
1252 int localNumRows = vspace->getNumIndices_Owned();
1253
1254 fei::SharedPtr<fei::VectorSpace> cspace = mgraph->getColSpace();
1255 int globalNumCols = globalNumRows;
1256 if (cspace.get() != NULL) {
1257 globalNumCols = cspace->getGlobalNumIndices();
1258 }
1259
1260 std::vector<int> indices_owned;
1261 int localNNZ = 0;
1262 CHK_ERR( vspace->getIndices_Owned(indices_owned) );
1263 int* rowsPtr = &indices_owned[0];
1264 for(int i=0; i<localNumRows; ++i) {
1265 int len;
1266 CHK_ERR( getRowLength(rowsPtr[i], len) );
1267 localNNZ += len;
1268 }
1269
1270 CHK_MPI( GlobalSum(getCommunicator(), localNNZ, globalNNZ) );
1271
1272 for(int p=0; p<numProcs(); ++p) {
1273 fei::Barrier(getCommunicator());
1274 if (p != localProc()) continue;
1275
1276 FEI_OFSTREAM* outFile = NULL;
1277 if (p==0) {
1278 outFile = new FEI_OFSTREAM(filename, IOS_OUT);
1279 FEI_OFSTREAM& ofs = *outFile;
1280 if (matrixMarketFormat) {
1281 ofs << mmbanner << FEI_ENDL;
1282 ofs <<globalNumRows << " " <<globalNumCols << " " <<globalNNZ <<FEI_ENDL;
1283 }
1284 else {
1285 ofs <<globalNumRows << " " <<globalNumCols <<FEI_ENDL;
1286 }
1287 }
1288 else outFile = new FEI_OFSTREAM(filename, IOS_APP);
1289
1290 outFile->setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
1291 outFile->precision(13);
1292 FEI_OFSTREAM& ofs = *outFile;
1293
1294 int rowLength;
1295
1296 for(int i=firstLocalOffset(); i<=lastLocalOffset(); ++i) {
1297 CHK_ERR( getRowLength(i, rowLength) );
1298
1299 work_indices_.resize(rowLength);
1300 work_data1D_.resize(rowLength);
1301
1302 int* indPtr = &work_indices_[0];
1303 double* coefPtr = &work_data1D_[0];
1304
1305 CHK_ERR( copyOutRow(i, rowLength, coefPtr, indPtr) );
1306
1307 fei::insertion_sort_with_companions<double>(rowLength,
1308 indPtr, coefPtr);
1309
1310 for(int j=0; j<rowLength; ++j) {
1311 if (matrixMarketFormat) {
1312 ofs << i+1 <<" "<<indPtr[j]+1<<" "<<coefPtr[j]<<FEI_ENDL;
1313 }
1314 else {
1315 ofs << i <<" "<<indPtr[j]<<" "<<coefPtr[j]<<FEI_ENDL;
1316 }
1317 }
1318 }
1319
1320 delete outFile;
1321 }
1322
1323 return(0);
1324}
1325
1326//----------------------------------------------------------------------------
1327template<typename T>
1329 bool matrixMarketFormat)
1330{
1331 fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
1332 if (mgraph.get() == NULL) {
1333 return(-1);
1334 }
1335
1336 if (haveFEMatrix()) {
1337 return(0);//temporary no-op...
1338 }
1339
1340 if (haveBlockMatrix()) {
1341 CHK_ERR( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
1342 }
1343
1344 static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
1345
1346 fei::SharedPtr<fei::VectorSpace> vspace = mgraph->getRowSpace();
1347
1348 int globalNNZ = 0;
1349 int localNumRows = vspace->getNumIndices_Owned();
1350 std::vector<int> indices_owned;
1351 int localNNZ = 0;
1352 CHK_ERR( vspace->getIndices_Owned(indices_owned));
1353 int* rowsPtr = &indices_owned[0];
1354 for(int i=0; i<localNumRows; ++i) {
1355 int len;
1356 CHK_ERR( getRowLength(rowsPtr[i], len) );
1357 localNNZ += len;
1358 }
1359
1360 CHK_MPI( fei::GlobalSum(getCommunicator(), localNNZ, globalNNZ) );
1361
1362 IOS_FMTFLAGS oldf = ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
1363
1364 for(int p=0; p<numProcs(); ++p) {
1365 fei::Barrier(getCommunicator());
1366 if (p != localProc()) continue;
1367
1368 if (p==0) {
1369 int globalSize = globalOffsets()[numProcs()]-1;
1370 if (matrixMarketFormat) {
1371 ostrm << mmbanner << FEI_ENDL;
1372 ostrm << globalSize << " " << globalSize << " " << globalNNZ << FEI_ENDL;
1373 }
1374 else {
1375 ostrm << globalSize << " " << globalSize << FEI_ENDL;
1376 }
1377 }
1378
1379 int rowLength;
1380
1381 for(int i=firstLocalOffset(); i<=lastLocalOffset(); ++i) {
1382 CHK_ERR( getRowLength(i, rowLength) );
1383
1384 work_indices_.resize(rowLength);
1385 work_data1D_.resize(rowLength);
1386
1387 int* indPtr = &work_indices_[0];
1388 double* coefPtr = &work_data1D_[0];
1389
1390 CHK_ERR( copyOutRow(i, rowLength, coefPtr, indPtr) );
1391
1392 for(int j=0; j<rowLength; ++j) {
1393 ostrm << i << " " << indPtr[j] << " " << coefPtr[j] << FEI_ENDL;
1394 }
1395 }
1396 }
1397
1398 ostrm.setf(oldf, IOS_FLOATFIELD);
1399
1400 return(0);
1401}
1402
1403#endif // _fei_Matrix_Impl_hpp_
1404
const fei::Pattern * getColPattern() const
const int * getColConnectivity(int ID) const
const int * getRowConnectivity(int ID) const
const std::map< int, int > & getConnectivityIDs() const
const fei::Pattern * getRowPattern() const
int multiply(fei::Vector *x, fei::Vector *y)
int writeToStream(FEI_OSTREAM &ostrm, bool matrixMarketFormat=true)
int gatherFromOverlap(bool accumulate=true)
int getLocalNumRows() const
int giveToUnderlyingMatrix(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, bool sumInto, int format)
Matrix_Impl(fei::SharedPtr< T > matrix, fei::SharedPtr< fei::MatrixGraph > matrixGraph, int numLocalEqns, bool zeroSharedRows=true)
int sumInFieldData(int fieldID, int idType, int rowID, int colID, const double *const *data, int format=0)
int putScalar(double scalar)
int getRowLength(int row, int &length) const
int getGlobalNumRows() const
int parameters(const fei::ParameterSet &paramset)
int copyOutRow(int row, int len, double *coefs, int *indices) const
int copyIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)
void setMatrixGraph(fei::SharedPtr< fei::MatrixGraph > matrixGraph)
int sumIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)
fei::SharedPtr< fei::MatrixGraph > getMatrixGraph() const
int writeToFile(const char *filename, bool matrixMarketFormat=true)
fei::SharedPtr< T > getMatrix()
int giveToUnderlyingBlockMatrix(int row, int rowDim, int numCols, const int *cols, const int *LDAs, const int *colDims, const double *const *values, bool sumInto)
const char * typeName()
const int * getNumFieldsPerID() const
int getNumIDs() const
const int * getNumIndicesPerID() const
int getNumIndices() const
int getTotalNumFields() const
const int * getFieldIDs() const
GlobalIDType getNumber() const
int getOwnerProc() const
unsigned getFieldSize(int fieldID)
int getPtEqnInfo(int ptEqn, int &blkEqn, int &blkOffset)
int getBlkEqnOffset(int blkEqn, int ptEqn)
int localProc(MPI_Comm comm)
int binarySearch(const T &item, const T *list, int len)
std::ostream & console_out()
int sortedListInsert(const T &item, std::vector< T > &list)
int numProcs(MPI_Comm comm)
int GlobalSum(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
static int getRowLength(T *mat, int row, int &length)
static int putValuesIn(T *mat, int numRows, const int *rows, int numCols, const int *cols, const double *const *values, bool sum_into)
static int copyOutRow(T *mat, int row, int len, double *coefs, int *indices)