00001
00002
00003
00004
00005
00006 #ifndef RANDOMSPARSEMATRIX_BUILDER_IMPL_HPP
00007 #define RANDOMSPARSEMATRIX_BUILDER_IMPL_HPP
00008
00009 #include "PlayaRandomSparseMatrixBuilderDecl.hpp"
00010 #include "PlayaIncrementallyConfigurableMatrixFactory.hpp"
00011 #include "PlayaLoadableMatrix.hpp"
00012
00013
00014 using namespace Playa;
00015 using namespace Teuchos;
00016
00017
00018 namespace Playa
00019 {
00020
00021 template <class Scalar>
00022 inline RandomSparseMatrixBuilder<Scalar>
00023 ::RandomSparseMatrixBuilder(int nLocalRows, int nLocalCols,
00024 double onProcDensity,
00025 double offProcDensity,
00026 const VectorType<double>& type)
00027 : OperatorBuilder<double>(nLocalRows, nLocalCols, type), op_()
00028 {
00029 initOp(onProcDensity, offProcDensity);
00030 }
00031
00032
00033 template <class Scalar>
00034 inline RandomSparseMatrixBuilder<Scalar>
00035 ::RandomSparseMatrixBuilder(const VectorSpace<Scalar>& d,
00036 const VectorSpace<Scalar>& r,
00037 double onProcDensity,
00038 double offProcDensity,
00039 const VectorType<double>& type)
00040 : OperatorBuilder<double>(d, r, type), op_()
00041 {
00042 initOp(onProcDensity, offProcDensity);
00043 }
00044
00045
00046 template <class Scalar>
00047 inline void RandomSparseMatrixBuilder<Scalar>
00048 ::initOp(double onProcDensity,
00049 double offProcDensity)
00050 {
00051 int rank = MPIComm::world().getRank();
00052 int nProc = MPIComm::world().getNProc();
00053
00054 RCP<MatrixFactory<double> > mFact
00055 = this->vecType().createMatrixFactory(this->domain(), this->range());
00056
00057 int colDimension = this->domain().dim();
00058 int rowDimension = this->range().dim();
00059 int numLocalCols = colDimension / nProc;
00060 int numLocalRows = rowDimension / nProc;
00061 int lowestLocalRow = numLocalRows * rank;
00062
00063 int lowestLocalCol = numLocalCols * rank;
00064 int highestLocalCol = numLocalCols * (rank+1) - 1;
00065
00066
00067 IncrementallyConfigurableMatrixFactory* icmf
00068 = dynamic_cast<IncrementallyConfigurableMatrixFactory*>(mFact.get());
00069 Array<Array<int> > colIndices(numLocalRows);
00070 for (int i=0; i<numLocalRows; i++)
00071 {
00072 int row = lowestLocalRow + i;
00073
00074 Array<int>& cols = colIndices[i];
00075
00076 while (cols.size() == 0)
00077 {
00078 for (int j=0; j<colDimension; j++)
00079 {
00080 double acceptProb;
00081 if (j >= lowestLocalCol && j <= highestLocalCol)
00082 {
00083 acceptProb = onProcDensity;
00084 }
00085 else
00086 {
00087 acceptProb = offProcDensity;
00088 }
00089 double p = 0.5*(ScalarTraits<double>::random() + 1.0);
00090
00091 if (p < acceptProb)
00092 {
00093 cols.append(j);
00094 }
00095 }
00096 if (cols.size()>0)
00097 {
00098 icmf->initializeNonzerosInRow(row, colIndices[i].size(),
00099 &(colIndices[i][0]));
00100 }
00101 }
00102
00103 }
00104 icmf->finalize();
00105
00106 op_ = mFact->createMatrix();
00107
00108 RCP<LoadableMatrix<double> > mat = op_.matrix();
00109
00110
00111 for (int i=0; i<numLocalRows; i++)
00112 {
00113 int row = lowestLocalRow + i;
00114 const Array<int>& cols = colIndices[i];
00115 Array<Scalar> colVals(cols.size());
00116 for (int j=0; j<cols.size(); j++)
00117 {
00118 colVals[j] = ScalarTraits<Scalar>::random();
00119 }
00120 if (cols.size() > 0)
00121 {
00122 mat->addToRow(row, colIndices[i].size(),
00123 &(colIndices[i][0]), &(colVals[0]));
00124 }
00125 }
00126 }
00127 }
00128
00129 #endif