My Project
_zgbmatrix-_zgbmatrix.hpp
1 //=============================================================================
3 inline _zgbmatrix operator+(const _zgbmatrix& matA, const _zgbmatrix& matB)
4 {
5 #ifdef CPPL_VERBOSE
6  std::cerr << "# [MARK] operator+(const _zgbmatrix&, const _zgbmatrix&)"
7  << std::endl;
8 #endif//CPPL_VERBOSE
9 
10 #ifdef CPPL_DEBUG
11  if(matA.N!=matB.N || matA.M!=matB.M){
12  std::cerr << "[ERROR] operator+(_zgbmatrix&, _zgbmatrix&)" << std::endl
13  << "These two matrises can not make a summation." << std::endl
14  << "Your input was (" << matA.M << "x" << matA.N << ") + ("
15  << matB.M << "x" << matB.N << ")." << std::endl;
16  exit(1);
17  }
18 #endif//CPPL_DEBUG
19 
20  if(matA.KL>matB.KL && matA.KU>matB.KU){
21  for(long i=0; i<matB.M; i++){
22  for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
23  matA(i,j)+=matB(i,j);
24  }
25  }
26 
27  matB.destroy();
28  return matA;
29  }
30 
31  else if(matB.KL>matA.KL && matB.KU>matA.KU){
32  for(long i=0; i<matA.M; i++){
33  for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
34  matB(i,j)+=matA(i,j);
35  }
36  }
37 
38  matA.destroy();
39  return matB;
40  }
41 
42  else{
43  zgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
44  newmat.zero();
45 
46  for(long i=0; i<matA.M; i++){
47  for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
48  newmat(i,j)+=matA(i,j);
49  }
50  for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
51  newmat(i,j)+=matB(i,j);
52  }
53  }
54 
55  matA.destroy();
56  matB.destroy();
57  return _(newmat);
58  }
59 }
60 
61 //=============================================================================
63 inline _zgbmatrix operator-(const _zgbmatrix& matA, const _zgbmatrix& matB)
64 {
65 #ifdef CPPL_VERBOSE
66  std::cerr << "# [MARK] operator-(const _zgbmatrix&, const _zgbmatrix&)"
67  << std::endl;
68 #endif//CPPL_VERBOSE
69 
70 #ifdef CPPL_DEBUG
71  if(matA.N!=matB.N || matA.M!=matB.M){
72  std::cerr << "[ERROR] operator-(_zgbmatrix&, _zgbmatrix&)" << std::endl
73  << "These two matrises can not make a subtraction." << std::endl
74  << "Your input was (" << matA.M << "x" << matA.N << ") - ("
75  << matB.M << "x" << matB.N << ")." << std::endl;
76  exit(1);
77  }
78 #endif//CPPL_DEBUG
79 
80  if(matA.KL>matB.KL && matA.KU>matB.KU){
81  for(long i=0; i<matB.M; i++){
82  for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
83  matA(i,j)-=matB(i,j);
84  }
85  }
86 
87  matB.destroy();
88  return matA;
89  }
90 
91  else{
92  zgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
93  newmat.zero();
94 
95  for(long i=0; i<matA.M; i++){
96  for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
97  newmat(i,j)+=matA(i,j);
98  }
99  for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
100  newmat(i,j)-=matB(i,j);
101  }
102  }
103 
104  matA.destroy();
105  matB.destroy();
106  return _(newmat);
107  }
108 }
109 
110 //=============================================================================
112 inline _zgbmatrix operator*(const _zgbmatrix& matA, const _zgbmatrix& matB)
113 {
114 #ifdef CPPL_VERBOSE
115  std::cerr << "# [MARK] operator*(const _zgbmatrix&, const _zgbmatrix&)"
116  << std::endl;
117 #endif//CPPL_VERBOSE
118 
119 #ifdef CPPL_DEBUG
120  if(matA.N!=matB.M){
121  std::cerr << "[ERROR] operator*(_zgbmatrix&, _zgbmatrix&)" << std::endl
122  << "These two matrises can not make a product." << std::endl
123  << "Your input was (" << matA.M << "x" << matA.N << ") * ("
124  << matB.M << "x" << matB.N << ")." << std::endl;
125  exit(1);
126  }
127 #endif//CPPL_DEBUG
128 
129  zgbmatrix newmat( matA.M, matB.N, matA.KU+matB.KL-1, matA.KL+matB.KU+1 );
130  newmat.zero();
131 
132  long i, j, k;
133 #pragma omp parallel for private(j,k)
134  for(i=0; i<newmat.m; i++){
135  for(j=max(0,i-newmat.kl); j<min(newmat.n,i+newmat.ku+1); j++){
136  for(k=max( max(0,i-matA.KL), max(0,j-matB.KU) );
137  k< min( min(matA.N,i+matA.KU+1), min(matB.M,j+matB.KL+1) ); k++){
138  newmat(i,j)+= matA(i,k)*matB(k,j);
139  }
140  }
141  }
142 
143  matA.destroy();
144  matB.destroy();
145  return _(newmat);
146 }
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
long KU
upper band width
Definition: _zgbmatrix.hpp:10
Complex Double-precision General Band Matrix Class.
Definition: zgbmatrix.hpp:3
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
friend const _zrovector & operator+(const _zrovector &)
Definition: _zrovector-unary.hpp:3
long KL
lower band width
Definition: _zgbmatrix.hpp:9