1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
#ifndef PARSPARSEMATRIX_H
#define PARSPARSEMATRIX_H
#include "mpi.h"
#include "MPI_Wrappers.h"
#include "SparseMatrix.h"
#include "DiagonalMatrix.h"
#include <algorithm>
namespace ATC_matrix {
/**
* @class ParSparseMatrix
* @brief Parallelized version of SparseMatrix class.
*
* ParSparseMatrix<double>::MultMv is used in LinearSolver, which is then
* used in NonLinearSolver, PoissonSolver, and SchrodingerSolver. These
* parallelized solvers are used in the following locations:
*
* - LinearSolver
* - ExtrinsicModelDriftDiffusion.cpp (lines 511 and 522)
* - AtomicRegulator.cpp (line 926)
* - TransferLibrary.cpp (lines 72 and 260)
* - PoissonSolver
* - ExtrinsicModelDriftDiffusion.cpp (line 232)
* - SchrodingerSolver
* - ExtrinsicModelDriftDiffusion.cpp (line 251)
* - SliceSchrodingerSolver
* - ExtrinsicModelDriftDiffusion.cpp (line 246)
*/
template <typename T>
class ParSparseMatrix : public SparseMatrix<T>
{
public:
ParSparseMatrix(MPI_Comm comm, INDEX rows = 0, INDEX cols = 0)
: SparseMatrix<T>(rows, cols), _comm(comm){}
ParSparseMatrix(MPI_Comm comm, const SparseMatrix<T> &c)
: SparseMatrix<T>(c), _comm(comm){}
ParSparseMatrix(MPI_Comm comm, INDEX* rows, INDEX* cols, T* vals,
INDEX size, INDEX nRows, INDEX nCols, INDEX nRowsCRS)
: SparseMatrix<T>(rows, cols, vals, size, nRows, nCols, nRowsCRS)
,_comm(comm){}
ParSparseMatrix(MPI_Comm comm)
: SparseMatrix<T>(), _comm(comm){}
virtual void operator=(const SparseMatrix<T> &source)
{
copy(source);
}
template<typename U>
friend void ParMultAB(MPI_Comm comm, const SparseMatrix<U>& A,
const Matrix<U>& B, DenseMatrix<U>& C);
private:
MPI_Comm _comm;
};
template <>
class ParSparseMatrix<double> : public SparseMatrix<double>
{
public:
// All the same constructors as for SparseMatrix
ParSparseMatrix(MPI_Comm comm, INDEX rows = 0, INDEX cols=0);
ParSparseMatrix(MPI_Comm comm, const SparseMatrix<double> &c);
ParSparseMatrix(MPI_Comm comm, INDEX* rows, INDEX* cols, double* vals, INDEX size,
INDEX nRows, INDEX nCols, INDEX nRowsCRS);
// Parallel sparse matrix multiplication functions
void MultMv(const Vector<double>& v, DenseVector<double>& c) const;
DenseVector<double> transMat(const Vector<double>& v) const;
void MultAB(const Matrix<double>& B, DenseMatrix<double>& C) const;
DenseMatrix<double> transMat(const DenseMatrix<double>& B) const;
DenseMatrix<double> transMat(const SparseMatrix<double>& B) const;
virtual void operator=(const SparseMatrix<double> &source);
template<typename U>
friend void ParMultAB(MPI_Comm comm, const SparseMatrix<U>& A, const Matrix<U>& B, DenseMatrix<U>& C);
private:
void partition(ParSparseMatrix<double>& A_local) const;
void finalize();
MPI_Comm _comm;
};
// The SparseMatrix versions of these functions will call the correct
// MultMv/MultAB:
// DenseVector<double> operator*(const ParSparseMatrix<double> &A, const Vector<double> &v);
// DenseVector<double> operator*(const Vector<double> &v, const ParSparseMatrix<double> &A);
// DenseMatrix<double> operator*(const ParSparseMatrix<double> &A, const Matrix<double> &B);
template<typename T>
void ParMultAB(MPI_Comm comm, const SparseMatrix<T>& A, const Matrix<T>& B, DenseMatrix<T>& C)
{
SparseMatrix<T>::compress(A);
INDEX M = A.nRows(), N = B.nCols();
if (!C.is_size(M, N))
{
C.resize(M, N);
C.zero();
}
// Temporarily put fields into a ParSparseMatrix for distributed multiplication
ParSparseMatrix<T> Ap(comm);
Ap._nRows = A._nRows;
Ap._nCols = A._nCols;
Ap._size = A._size;
Ap._nRowsCRS = A._nRowsCRS;
Ap._val = A._val;
Ap._ja = A._ja;
Ap._ia = A._ia;
Ap.hasTemplate_ = A.hasTemplate_;
// MultAB calls compress(), but we hope that does nothing because we just
// compressed A. If it did something, it might mess up other members
// (e.g. _tri).
Ap.MultAB(B, C);
// We're not changing the matrix's values, so we can justify calling A const.
SparseMatrix<T> &Avar = const_cast<SparseMatrix<T> &>(A);
Avar._nRows = Ap._nRows;
Avar._nCols = Ap._nCols;
Avar._size = Ap._size;
Avar._nRowsCRS = Ap._nRowsCRS;
Avar._val = Ap._val;
Avar._ja = Ap._ja;
Avar._ia = Ap._ia;
Avar.hasTemplate_ = Ap.hasTemplate_;
// Avoid catastrophe
Ap._val = NULL;
Ap._ja = NULL;
Ap._ia = NULL;
}
/*SparseMatrix<double> operator*(const ParSparseMatrix<double> &A, const SparseMatrix<double> &B);
SparseMatrix<double> operator*(const ParSparseMatrix<double> &A, const DiagonalMatrix<double> &B);
*/
} // end namespace
#endif