My Project
zgbmatrix-_zgbmatrix.hpp
1 //=============================================================================
4 {
5 #ifdef CPPL_VERBOSE
6  std::cerr << "# [MARK] zgbmatrix::operator=(const _zgbmatrix&)"
7  << std::endl;
8 #endif//CPPL_VERBOSE
9 
10  shallow_copy(mat);
11  return *this;
12 }
13 
17 
18 //=============================================================================
22 {
23 #ifdef CPPL_VERBOSE
24  std::cerr << "# [MARK] zgbmatrix::operator+=(const _zgbmatrix&)"
25  << std::endl;
26 #endif//CPPL_VERBOSE
27 
28 #ifdef CPPL_DEBUG
29  if(N!=mat.N || M!=mat.M){
30  std::cerr << "[ERROR] zgbmatrix::operator+=(_zgbmatrix&)" << 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  zgbmatrix 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 //=============================================================================
72 {
73 #ifdef CPPL_VERBOSE
74  std::cerr << "# [MARK] zgbmatrix::operator-=(const _zgbmatrix&)"
75  << std::endl;
76 #endif//CPPL_VERBOSE
77 
78 #ifdef CPPL_DEBUG
79  if(N!=mat.N || M!=mat.M){
80  std::cerr << "[ERROR] zgbmatrix::operator-=(_zgbmatrix&)" << 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  zgbmatrix 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 //=============================================================================
121 {
122 #ifdef CPPL_VERBOSE
123  std::cerr << "# [MARK] zgbmatrix::operator*=(const _zgbmatrix&)"
124  << std::endl;
125 #endif//CPPL_VERBOSE
126 
127 #ifdef CPPL_DEBUG
128  if(N!=mat.M){
129  std::cerr << "[ERROR] zgbmatrix::operator*=(_zgbmatrix&)" << 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  zgbmatrix 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 _zgbmatrix operator+(const zgbmatrix& matA, const _zgbmatrix& matB)
161 {
162 #ifdef CPPL_VERBOSE
163  std::cerr << "# [MARK] operator+(const zgbmatrix&, const _zgbmatrix&)"
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+(zgbmatrix&, _zgbmatrix&)" << 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  zgbmatrix 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 _zgbmatrix operator-(const zgbmatrix& matA, const _zgbmatrix& matB)
207 {
208 #ifdef CPPL_VERBOSE
209  std::cerr << "# [MARK] operator-(const zgbmatrix&, const _zgbmatrix&)"
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-(zgbmatrix&, _zgbmatrix&)" << 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  zgbmatrix 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 _zgbmatrix operator*(const zgbmatrix& matA, const _zgbmatrix& matB)
242 {
243 #ifdef CPPL_VERBOSE
244  std::cerr << "# [MARK] operator*(const zgbmatrix&, const _zgbmatrix&)"
245  << std::endl;
246 #endif//CPPL_VERBOSE
247 
248 #ifdef CPPL_DEBUG
249  if(matA.N!=matB.M){
250  std::cerr << "[ERROR] operator*(zgbmatrix&, _zgbmatrix&)" << 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  zgbmatrix 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 }
std::complex< double > & operator()(const long &, const long &)
Definition: zgbmatrix-io.hpp:3
friend _zgematrix i(const zgbmatrix &)
Definition: zgbmatrix-calc.hpp:22
zgbmatrix & operator*=(const zgbmatrix &)
Definition: zgbmatrix-zgbmatrix.hpp:121
friend _zrovector operator-(const _zrovector &)
Definition: _zrovector-unary.hpp:15
long N
matrix column size
Definition: _zgbmatrix.hpp:8
void destroy() const
Definition: _zgbmatrix-misc.hpp:3
zgbmatrix & operator+=(const zgbmatrix &)
Definition: zgbmatrix-zgbmatrix.hpp:24
long KU
upper band width
Definition: _zgbmatrix.hpp:10
zgbmatrix & operator=(const zgbmatrix &)
Definition: zgbmatrix-zgbmatrix.hpp:4
Complex Double-precision General Band Matrix Class.
Definition: zgbmatrix.hpp:3
zgbmatrix & operator-=(const zgbmatrix &)
Definition: zgbmatrix-zgbmatrix.hpp:73
friend _zrovector operator*(const zrovector &, const zgematrix &)
Definition: zrovector-zgematrix.hpp:3
long M
matrix row size
Definition: _zgbmatrix.hpp:7
(DO NOT USE) Smart-temporary Complex Double-precision General Band Matrix Class
Definition: _zgbmatrix.hpp:3
long const & kl
lower band width (readable)
Definition: zgbmatrix.hpp:18
void zero()
Definition: zgbmatrix-misc.hpp:29
long const & n
matrix column size (readable)
Definition: zgbmatrix.hpp:17
long const & m
matrix row size (readable)
Definition: zgbmatrix.hpp:16
friend const _zrovector & operator+(const _zrovector &)
Definition: _zrovector-unary.hpp:3
friend void swap(zgbmatrix &, zgbmatrix &)
Definition: zgbmatrix-misc.hpp:167
long KL
lower band width
Definition: _zgbmatrix.hpp:9
long const & ku
upper band width (readable)
Definition: zgbmatrix.hpp:19
void shallow_copy(const _zgbmatrix &)
Definition: zgbmatrix-misc.hpp:107