VERB_code_2.3
_dgbmatrix-_dgbmatrix.hpp
1 //=============================================================================
3 inline _dgbmatrix operator+(const _dgbmatrix& matA, const _dgbmatrix& matB)
4 {
5 #ifdef CPPL_VERBOSE
6  std::cerr << "# [MARK] operator+(const _dgbmatrix&, const _dgbmatrix&)"
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+(_dgbmatrix&, _dgbmatrix&)" << 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  dgbmatrix 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 _dgbmatrix operator-(const _dgbmatrix& matA, const _dgbmatrix& matB)
64 {
65 #ifdef CPPL_VERBOSE
66  std::cerr << "# [MARK] operator-(const _dgbmatrix&, const _dgbmatrix&)"
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-(_dgbmatrix&, _dgbmatrix&)" << 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  dgbmatrix 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 _dgbmatrix operator*(const _dgbmatrix& matA, const _dgbmatrix& matB)
113 {
114 #ifdef CPPL_VERBOSE
115  std::cerr << "# [MARK] operator*(const _dgbmatrix&, const _dgbmatrix&)"
116  << std::endl;
117 #endif//CPPL_VERBOSE
118 
119 #ifdef CPPL_DEBUG
120  if(matA.N!=matB.M){
121  std::cerr << "[ERROR] operator*(_dgbmatrix&, _dgbmatrix&)" << 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  dgbmatrix 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 }
long const & kl
lower band width (readable)
Definition: dgbmatrix.hpp:18
double max(double v1, double v2)
Return maximum.
Definition: variousFunctions.cpp:355
long KL
lower band width
Definition: _dgbmatrix.hpp:9
long M
matrix row size
Definition: _dgbmatrix.hpp:7
void destroy() const
Definition: _dgbmatrix-misc.hpp:3
long N
matrix column size
Definition: _dgbmatrix.hpp:8
friend _drovector operator-(const _drovector &)
Definition: _drovector-unary.hpp:15
Real Double-precision General Band Matrix Class.
Definition: dgbmatrix.hpp:3
void zero()
Definition: dgbmatrix-misc.hpp:29
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
long const & m
matrix row size (readable)
Definition: dgbmatrix.hpp:16
long const & ku
upper band width (readable)
Definition: dgbmatrix.hpp:19
long KU
upper band width
Definition: _dgbmatrix.hpp:10
friend const _drovector & operator+(const _drovector &)
Definition: _drovector-unary.hpp:3
long const & n
matrix column size (readable)
Definition: dgbmatrix.hpp:17