My Project
dgbmatrix-_dgbmatrix.hpp
1 //=============================================================================
3 inline dgbmatrix& dgbmatrix::operator=(const _dgbmatrix& mat)
4 {
5 #ifdef CPPL_VERBOSE
6  std::cerr << "# [MARK] dgbmatrix::operator=(const _dgbmatrix&)"
7  << std::endl;
8 #endif//CPPL_VERBOSE
9 
10  shallow_copy(mat);
11  return *this;
12 }
13 
17 
18 //=============================================================================
21 inline dgbmatrix& dgbmatrix::operator+=(const _dgbmatrix& mat)
22 {
23 #ifdef CPPL_VERBOSE
24  std::cerr << "# [MARK] dgbmatrix::operator+=(const _dgbmatrix&)"
25  << std::endl;
26 #endif//CPPL_VERBOSE
27 
28 #ifdef CPPL_DEBUG
29  if(N!=mat.N || M!=mat.M){
30  std::cerr << "[ERROR] dgbmatrix::operator+=(_dgbmatrix&)" << std::endl
31  << "These two matrises can not make a summation." << std::endl
32  << "Your input was"
33  << "(" << M <<"x"<< N <<","<< KL <<":"<< KU << ") "<< "+="
34  << "("<< mat.M <<"x"<< mat.N <<","<< mat.KL <<":"<< mat.KU <<") "
35  << std::endl;
36  exit(1);
37  }
38 #endif//CPPL_DEBUG
39 
40  if(KL>=mat.KL && KU>=mat.KU){
41  for(long i=0; i<M; i++){
42  for(long j=max(0,i-mat.KL); j<min(N,i+mat.KU+1); j++){
43  operator()(i,j)+=mat(i,j);
44  }
45  }
46 
47  mat.destroy();
48  return *this;
49  }
50  else{
51  dgbmatrix newmat(M,N,max(KL,mat.KL),max(KU,mat.KU));
52  newmat.zero();
53  for(long i=0; i<M; i++){
54  for(long j=max(0,i-KL); j<min(N,i+KU+1); j++){
55  newmat(i,j)+=operator()(i,j);
56  }
57  for(long j=max(0,i-mat.KL); j<min(mat.N,i+mat.KU+1); j++){
58  newmat(i,j)+=mat(i,j);
59  }
60  }
61 
62  swap(*this,newmat);
63  mat.destroy();
64  return *this;
65  }
66 }
67 
68 //=============================================================================
71 inline dgbmatrix& dgbmatrix::operator-=(const _dgbmatrix& mat)
72 {
73 #ifdef CPPL_VERBOSE
74  std::cerr << "# [MARK] dgbmatrix::operator-=(const _dgbmatrix&)"
75  << std::endl;
76 #endif//CPPL_VERBOSE
77 
78 #ifdef CPPL_DEBUG
79  if(N!=mat.N || M!=mat.M){
80  std::cerr << "[ERROR] dgbmatrix::operator-=(_dgbmatrix&)" << std::endl
81  << "These two matrises can not make a subtraction." << std::endl
82  << "Your input was"
83  << "(" << M <<"x"<< N <<","<< KL <<":"<< KU << ") "<< "-="
84  << "("<< mat.M <<"x"<< mat.N <<","<< mat.KL <<":"<< mat.KU <<") "
85  << std::endl;
86  exit(1);
87  }
88 #endif//CPPL_DEBUG
89 
90  if(KL>=mat.KL && KU>=mat.KU){
91  for(long i=0; i<M; i++){
92  for(long j=max(0,i-mat.KL); j<min(N,i+mat.KU+1); j++){
93  operator()(i,j)-=mat(i,j);
94  }
95  }
96 
97  mat.destroy();
98  return *this;
99  }
100  else{
101  dgbmatrix newmat(M,N,max(KL,mat.KL),max(KU,mat.KU));
102  newmat.zero();
103  for(long i=0; i<M; i++){
104  for(long j=max(0,i-KL); j<min(N,i+KU+1); j++){
105  newmat(i,j)+=operator()(i,j);
106  }
107  for(long j=max(0,i-mat.KL); j<min(mat.N,i+mat.KU+1); j++){
108  newmat(i,j)-=mat(i,j);
109  }
110  }
111 
112  swap(*this,newmat);
113  mat.destroy();
114  return *this;
115  }
116 }
117 
118 //=============================================================================
120 inline dgbmatrix& dgbmatrix::operator*=(const _dgbmatrix& mat)
121 {
122 #ifdef CPPL_VERBOSE
123  std::cerr << "# [MARK] dgbmatrix::operator*=(const _dgbmatrix&)"
124  << std::endl;
125 #endif//CPPL_VERBOSE
126 
127 #ifdef CPPL_DEBUG
128  if(N!=mat.M){
129  std::cerr << "[ERROR] dgbmatrix::operator*=(_dgbmatrix&)" << std::endl
130  << "These two matrises can not make a product." << std::endl
131  << "Your input was (" << M << "x" << N << ") * ("
132  << mat.M << "x" << mat.N << ")." << std::endl;
133  exit(1);
134  }
135 #endif//CPPL_DEBUG
136 
137  dgbmatrix newmat( M, mat.N, KU+mat.KL-1, KL+mat.KU+1 );
138  newmat.zero();
139 
140  for(long i=0; i<newmat.M; i++){
141  for(long j=max(0,i-newmat.KL); j<min(newmat.N,i+newmat.KU+1); j++){
142  for(long k=max( max(0,i-KL), max(0,j-mat.KU) );
143  k< min( min(N,i+KU+1), min(mat.M,j+mat.KL+1) ); k++){
144  newmat(i,j)+= operator()(i,k)*mat(k,j);
145  }
146  }
147  }
148 
149  swap(*this,newmat);
150  mat.destroy();
151  return *this;
152 }
153 
157 
158 //=============================================================================
160 inline _dgbmatrix operator+(const dgbmatrix& matA, const _dgbmatrix& matB)
161 {
162 #ifdef CPPL_VERBOSE
163  std::cerr << "# [MARK] operator+(const dgbmatrix&, const _dgbmatrix&)"
164  << std::endl;
165 #endif//CPPL_VERBOSE
166 
167 #ifdef CPPL_DEBUG
168  if(matA.N!=matB.N || matA.M!=matB.M){
169  std::cerr << "[ERROR] operator+(dgbmatrix&, _dgbmatrix&)" << std::endl
170  << "These two matrises can not make a summation." << std::endl
171  << "Your input was (" << matA.M << "x" << matA.N << ") + ("
172  << matB.M << "x" << matB.N << ")." << std::endl;
173  exit(1);
174  }
175 #endif//CPPL_DEBUG
176 
177  if(matB.KL>matA.KL && matB.KU>matA.KU){
178  for(long i=0; i<matA.M; i++){
179  for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
180  matB(i,j)+=matA(i,j);
181  }
182  }
183 
184  return matB;
185  }
186  else{
187  dgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
188  newmat.zero();
189 
190  for(long i=0; i<matA.M; i++){
191  for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
192  newmat(i,j)+=matA(i,j);
193  }
194  for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
195  newmat(i,j)+=matB(i,j);
196  }
197  }
198 
199  matB.destroy();
200  return _(newmat);
201  }
202 }
203 
204 //=============================================================================
206 inline _dgbmatrix operator-(const dgbmatrix& matA, const _dgbmatrix& matB)
207 {
208 #ifdef CPPL_VERBOSE
209  std::cerr << "# [MARK] operator-(const dgbmatrix&, const _dgbmatrix&)"
210  << std::endl;
211 #endif//CPPL_VERBOSE
212 
213 #ifdef CPPL_DEBUG
214  if(matA.N!=matB.N || matA.M!=matB.M){
215  std::cerr << "[ERROR] operator-(dgbmatrix&, _dgbmatrix&)" << std::endl
216  << "These two matrises can not make a subtraction." << std::endl
217  << "Your input was (" << matA.M << "x" << matA.N << ") - ("
218  << matB.M << "x" << matB.N << ")." << std::endl;
219  exit(1);
220  }
221 #endif//CPPL_DEBUG
222 
223  dgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
224  newmat.zero();
225 
226  for(long i=0; i<matA.M; i++){
227  for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
228  newmat(i,j)+=matA(i,j);
229  }
230  for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
231  newmat(i,j)-=matB(i,j);
232  }
233  }
234 
235  matB.destroy();
236  return _(newmat);
237 }
238 
239 //=============================================================================
241 inline _dgbmatrix operator*(const dgbmatrix& matA, const _dgbmatrix& matB)
242 {
243 #ifdef CPPL_VERBOSE
244  std::cerr << "# [MARK] operator*(const dgbmatrix&, const _dgbmatrix&)"
245  << std::endl;
246 #endif//CPPL_VERBOSE
247 
248 #ifdef CPPL_DEBUG
249  if(matA.N!=matB.M){
250  std::cerr << "[ERROR] operator*(dgbmatrix&, _dgbmatrix&)" << std::endl
251  << "These two matrises can not make a product." << std::endl
252  << "Your input was (" << matA.M << "x" << matA.N << ") * ("
253  << matB.M << "x" << matB.N << ")." << std::endl;
254  exit(1);
255  }
256 #endif//CPPL_DEBUG
257 
258  dgbmatrix newmat( matA.M, matB.N, matA.KU+matB.KL-1, matA.KL+matB.KU+1 );
259  newmat.zero();
260 
261  long i, j, k;
262 #pragma omp parallel for private(j,k)
263  for(i=0; i<newmat.M; i++){
264  for(j=max(0,i-newmat.KL); j<min(newmat.N,i+newmat.KU+1); j++){
265  for(k=max( max(0,i-matA.KL), max(0,j-matB.KU) );
266  k< min( min(matA.N,i+matA.KU+1), min(matB.M,j+matB.KL+1) ); k++){
267  newmat(i,j)+= matA(i,k)*matB(k,j);
268  }
269  }
270  }
271 
272  matB.destroy();
273  return _(newmat);
274 }
dgbmatrix & operator=(const dgbmatrix &)
Definition: dgbmatrix-dgbmatrix.hpp:4
friend _dgematrix i(const dgbmatrix &)
Definition: dgbmatrix-calc.hpp:22
long KL
lower band width
Definition: _dgbmatrix.hpp:9
void shallow_copy(const _dgbmatrix &)
Definition: dgbmatrix-misc.hpp:107
long M
matrix row size
Definition: _dgbmatrix.hpp:7
void destroy() const
Definition: _dgbmatrix-misc.hpp:3
double & operator()(const long &, const long &)
Definition: dgbmatrix-io.hpp:3
long N
matrix column size
Definition: _dgbmatrix.hpp:8
dgbmatrix & operator+=(const dgbmatrix &)
Definition: dgbmatrix-dgbmatrix.hpp:24
friend _drovector operator-(const _drovector &)
Definition: _drovector-unary.hpp:15
Real Double-precision General Band Matrix Class.
Definition: dgbmatrix.hpp:3
dgbmatrix & operator*=(const dgbmatrix &)
Definition: dgbmatrix-dgbmatrix.hpp:121
friend _drovector operator*(const drovector &, const dgematrix &)
Definition: drovector-dgematrix.hpp:3
(DO NOT USE) Smart-temporary Real Double-precision General Band Matrix Class
Definition: _dgbmatrix.hpp:3
dgbmatrix & operator-=(const dgbmatrix &)
Definition: dgbmatrix-dgbmatrix.hpp:73
long KU
upper band width
Definition: _dgbmatrix.hpp:10
friend const _drovector & operator+(const _drovector &)
Definition: _drovector-unary.hpp:3
long M
matrix row size
Definition: _dgematrix.hpp:7
friend void swap(dgbmatrix &, dgbmatrix &)
Definition: dgbmatrix-misc.hpp:167