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  return matA;
28  }
29 
30  else{
31  dgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
32  newmat.zero();
33 
34  for(long i=0; i<matA.M; i++){
35  for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
36  newmat(i,j)+=matA(i,j);
37  }
38  for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
39  newmat(i,j)+=matB(i,j);
40  }
41  }
42 
43  matA.destroy();
44  return _(newmat);
45  }
46 }
47 
48 //=============================================================================
50 inline _dgbmatrix operator-(const _dgbmatrix& matA, const dgbmatrix& matB)
51 {
52 #ifdef CPPL_VERBOSE
53  std::cerr << "# [MARK] operator-(const _dgbmatrix&, const dgbmatrix&)"
54  << std::endl;
55 #endif//CPPL_VERBOSE
56 
57 #ifdef CPPL_DEBUG
58  if(matA.N!=matB.N || matA.M!=matB.M){
59  std::cerr << "[ERROR] operator-(_dgbmatrix&, dgbmatrix&)" << std::endl
60  << "These two matrises can not make a subtraction." << std::endl
61  << "Your input was (" << matA.M << "x" << matA.N << ") - ("
62  << matB.M << "x" << matB.N << ")." << std::endl;
63  exit(1);
64  }
65 #endif//CPPL_DEBUG
66 
67  if(matA.KL>=matB.KL && matA.KU>=matB.KU){
68  for(long i=0; i<matB.M; i++){
69  for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
70  matA(i,j)-=matB(i,j);
71  }
72  }
73 
74  return matA;
75  }
76 
77  else{
78  dgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
79  newmat.zero();
80 
81  for(long i=0; i<matA.M; i++){
82  for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
83  newmat(i,j)+=matA(i,j);
84  }
85  for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
86  newmat(i,j)-=matB(i,j);
87  }
88  }
89 
90  return _(newmat);
91  }
92 }
93 
94 //=============================================================================
96 inline _dgbmatrix operator*(const _dgbmatrix& matA, const dgbmatrix& matB)
97 {
98 #ifdef CPPL_VERBOSE
99  std::cerr << "# [MARK] operator*(const _dgbmatrix&, const dgbmatrix&)"
100  << std::endl;
101 #endif//CPPL_VERBOSE
102 
103 #ifdef CPPL_DEBUG
104  if(matA.N!=matB.M){
105  std::cerr << "[ERROR] operator*(_dgbmatrix&, dgbmatrix&)" << std::endl
106  << "These two matrises can not make a product." << std::endl
107  << "Your input was (" << matA.M << "x" << matA.N << ") * ("
108  << matB.M << "x" << matB.N << ")." << std::endl;
109  exit(1);
110  }
111 #endif//CPPL_DEBUG
112 
113  dgbmatrix newmat( matA.M, matB.N, matA.KU+matB.KL-1, matA.KL+matB.KU+1 );
114  newmat.zero();
115 
116  long i, j, k;
117 #pragma omp parallel for private(j,k)
118  for(i=0; i<newmat.m; i++){
119  for(j=max(0,i-newmat.kl); j<min(newmat.n,i+newmat.ku+1); j++){
120  for(k=max( max(0,i-matA.KL), max(0,j-matB.KU) );
121  k< min( min(matA.N,i+matA.KU+1), min(matB.M,j+matB.KL+1) ); k++){
122  newmat(i,j)+= matA(i,k)*matB(k,j);
123  }
124  }
125  }
126 
127  matA.destroy();
128  return _(newmat);
129 }
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