CPPLapack
 All Classes Files Functions Variables Friends Pages
dgbmatrix-dgbmatrix.hpp
Go to the documentation of this file.
1 //=============================================================================
2 /*! dgbmatrix=dgbmatrix operator\n
3  The left side matrix is overwritten thoroughly including band width. */
5 {CPPL_VERBOSE_REPORT;
6  if(array!=mat.array){ // if it is NOT self substitution
7  copy(mat);
8  }
9  return *this;
10 }
11 
12 ///////////////////////////////////////////////////////////////////////////////
13 ///////////////////////////////////////////////////////////////////////////////
14 ///////////////////////////////////////////////////////////////////////////////
15 
16 //=============================================================================
17 /*! dgbmatrix+=dgbmatrix operator\n
18  If the band width of the left side matrix is narrower than the right side matrix, the band width of the left side matrix become thicker as same as the right side matrix. */
20 {CPPL_VERBOSE_REPORT;
21 #ifdef CPPL_DEBUG
22  if(n!=mat.n || m!=mat.m){
23  ERROR_REPORT;
24  std::cerr << "These two matrises can not make a summation." << std::endl
25  << "Your input was" << "(" << m <<"x"<< n <<","<< kl <<":"<< ku << ") "<< "+=" << "("<< mat.m <<"x"<< mat.n <<","<< mat.kl <<":"<< mat.ku <<") " << std::endl;
26  exit(1);
27  }
28 #endif//CPPL_DEBUG
29 
30  if(kl>=mat.kl && ku>=mat.ku){
31  for(CPPL_INT i=0; i<m; i++){
32  const CPPL_INT jmax =std::min(n,i+mat.ku+1);
33  for(CPPL_INT j=std::max(CPPL_INT(0),i-mat.kl); j<jmax; j++){
34  operator()(i,j) += mat(i,j);
35  }
36  }
37 
38  return *this;
39  }
40 
41  else{
42  dgbmatrix newmat(m,n,std::max(kl,mat.kl),std::max(ku,mat.ku));
43  newmat.zero();
44  for(CPPL_INT i=0; i<m; i++){
45  const CPPL_INT jmax1 =std::min(n,i+ku+1);
46  for(CPPL_INT j=std::max(CPPL_INT(0),i-kl); j<jmax1; j++){
47  newmat(i,j) += operator()(i,j);
48  }
49  const CPPL_INT jmax2 =std::min(mat.n,i+mat.ku+1);
50  for(CPPL_INT j=std::max(CPPL_INT(0),i-mat.kl); j<jmax2; j++){
51  newmat(i,j) += mat(i,j);
52  }
53  }
54 
55  swap(*this,newmat);
56  return *this;
57  }
58 }
59 
60 //=============================================================================
61 /*! dgbmatrix-=dgbmatrix operator\n
62  If the band width of the left side matrix is narrower than the right side matrix, the band width of the left side matrix become thicker as same as the right side matrix. */
64 {CPPL_VERBOSE_REPORT;
65 #ifdef CPPL_DEBUG
66  if(n!=mat.n || m!=mat.m){
67  ERROR_REPORT;
68  std::cerr << "These two matrises can not make a subtraction." << std::endl
69  << "Your input was" << "(" << m <<"x"<< n <<","<< kl <<":"<< ku << ") "<< "-=" << "("<< mat.m <<"x"<< mat.n <<","<< mat.kl <<":"<< mat.ku <<") " << std::endl;
70  exit(1);
71  }
72 #endif//CPPL_DEBUG
73 
74  if(kl>=mat.kl && ku>=mat.ku){
75  for(CPPL_INT i=0; i<m; i++){
76  const CPPL_INT jmax =std::min(n,i+mat.ku+1);
77  for(CPPL_INT j=std::max(CPPL_INT(0),i-mat.kl); j<jmax; j++){
78  operator()(i,j) -= mat(i,j);
79  }
80  }
81 
82  return *this;
83  }
84 
85  else{
86  dgbmatrix newmat(m,n,std::max(kl,mat.kl),std::max(ku,mat.ku));
87  newmat.zero();
88  for(CPPL_INT i=0; i<m; i++){
89  const CPPL_INT jmax1 =std::min(n,i+ku+1);
90  for(CPPL_INT j=std::max(CPPL_INT(0),i-kl); j<jmax1; j++){
91  newmat(i,j) += operator()(i,j);
92  }
93  const CPPL_INT jmax2 =std::min(mat.n,i+mat.ku+1);
94  for(CPPL_INT j=std::max(CPPL_INT(0),i-mat.kl); j<jmax2; j++){
95  newmat(i,j) -= mat(i,j);
96  }
97  }
98 
99  swap(*this,newmat);
100  return *this;
101  }
102 }
103 
104 //=============================================================================
105 /*! dgbmatrix*=dgbmatrix operator */
107 {CPPL_VERBOSE_REPORT;
108 #ifdef CPPL_DEBUG
109  if(n!=mat.m){
110  ERROR_REPORT;
111  std::cerr << "These two matrises can not make a product." << std::endl
112  << "Your input was (" << m << "x" << n << ") * (" << mat.m << "x" << mat.n << ")." << std::endl;
113  exit(1);
114  }
115 #endif//CPPL_DEBUG
116 
117  dgbmatrix newmat( m, mat.n, std::min(kl+mat.kl, m-1), std::min(ku+mat.ku, mat.n-1) );
118  newmat.zero();
119 
120  for(CPPL_INT i=0; i<newmat.m; i++){
121  const CPPL_INT jmax =std::min(newmat.n,i+newmat.ku+1);
122  for(CPPL_INT j=std::max(CPPL_INT(0),i-newmat.kl); j<jmax; j++){
123  const CPPL_INT kmax =std::min( std::min(n,i+ku+1), std::min(mat.m,j+mat.kl+1) );
124  for(CPPL_INT k=std::max( std::max(CPPL_INT(0),i-kl), std::max(CPPL_INT(0),j-mat.ku) ); k<kmax; k++){
125  newmat(i,j) += operator()(i,k)*mat(k,j);
126  }
127  }
128  }
129 
130  swap(*this,newmat);
131  return *this;
132 }
133 
134 ///////////////////////////////////////////////////////////////////////////////
135 ///////////////////////////////////////////////////////////////////////////////
136 ///////////////////////////////////////////////////////////////////////////////
137 
138 //=============================================================================
139 /*! dgbmatrix+dgbmatrix operator */
140 inline _dgbmatrix operator+(const dgbmatrix& matA, const dgbmatrix& matB)
141 {CPPL_VERBOSE_REPORT;
142 #ifdef CPPL_DEBUG
143  if(matA.n!=matB.n || matA.m!=matB.m){
144  ERROR_REPORT;
145  std::cerr << "These two matrises can not make a summation." << std::endl
146  << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.m << "x" << matB.n << ")." << std::endl;
147  exit(1);
148  }
149 #endif//CPPL_DEBUG
150 
151  dgbmatrix newmat(matA.m,matA.n,std::max(matA.kl,matB.kl),std::max(matA.ku,matB.ku));
152  newmat.zero();
153 
154  for(CPPL_INT i=0; i<matA.m; i++){
155  const CPPL_INT jmax1 =std::min(matA.n,i+matA.ku+1);
156  for(CPPL_INT j=std::max(CPPL_INT(0),i-matA.kl); j<jmax1; j++){
157  newmat(i,j) += matA(i,j);
158  }
159  const CPPL_INT jmax2 =std::min(matB.n,i+matB.ku+1);
160  for(CPPL_INT j=std::max(CPPL_INT(0),i-matB.kl); j<jmax2; j++){
161  newmat(i,j) += matB(i,j);
162  }
163  }
164 
165  return _(newmat);
166 }
167 
168 //=============================================================================
169 /*! dgbmatrix-dgbmatrix operator */
170 inline _dgbmatrix operator-(const dgbmatrix& matA, const dgbmatrix& matB)
171 {CPPL_VERBOSE_REPORT;
172 #ifdef CPPL_DEBUG
173  if(matA.n!=matB.n || matA.m!=matB.m){
174  ERROR_REPORT;
175  std::cerr << "These two matrises can not make a subtraction." << std::endl
176  << "Your input was (" << matA.m << "x" << matA.n << ") - (" << matB.m << "x" << matB.n << ")." << std::endl;
177  exit(1);
178  }
179 #endif//CPPL_DEBUG
180 
181  dgbmatrix newmat(matA.m,matA.n,std::max(matA.kl,matB.kl),std::max(matA.ku,matB.ku));
182  newmat.zero();
183 
184  for(CPPL_INT i=0; i<matA.m; i++){
185  const CPPL_INT jmax1 =std::min(matA.n,i+matA.ku+1);
186  for(CPPL_INT j=std::max(CPPL_INT(0),i-matA.kl); j<jmax1; j++){
187  newmat(i,j) += matA(i,j);
188  }
189  const CPPL_INT jmax2 =std::min(matB.n,i+matB.ku+1);
190  for(CPPL_INT j=std::max(CPPL_INT(0),i-matB.kl); j<jmax2; j++){
191  newmat(i,j) -= matB(i,j);
192  }
193  }
194 
195  return _(newmat);
196 }
197 
198 //=============================================================================
199 /*! dgbmatrix*dgbmatrix operator */
200 inline _dgbmatrix operator*(const dgbmatrix& matA, const dgbmatrix& matB)
201 {CPPL_VERBOSE_REPORT;
202 #ifdef CPPL_DEBUG
203  if(matA.n!=matB.m){
204  ERROR_REPORT;
205  std::cerr << "These two matrises can not make a product." << std::endl
206  << "Your input was (" << matA.m << "x" << matA.n << ") * (" << matB.m << "x" << matB.n << ")." << std::endl;
207  exit(1);
208  }
209 #endif//CPPL_DEBUG
210 
211  dgbmatrix newmat( matA.m, matB.n, std::min(matA.kl+matB.kl,matA.m-1), std::min(matA.ku+matB.ku,matB.n-1) );
212  newmat.zero();
213 
214  for(CPPL_INT i=0; i<newmat.m; i++){
215  const CPPL_INT jmax =std::min(newmat.n,i+newmat.ku+1);
216  for(CPPL_INT j=std::max(CPPL_INT(0),i-newmat.kl); j<jmax; j++){
217  const CPPL_INT kmax =std::min( std::min(matA.n,i+matA.ku+1), std::min(matB.m,j+matB.kl+1) );
218  for(CPPL_INT k=std::max( std::max(CPPL_INT(0),i-matA.kl), std::max(CPPL_INT(0),j-matB.ku) ); k<kmax; k++){
219  newmat(i,j) += matA(i,k)*matB(k,j);
220  }
221  }
222  }
223 
224  return _(newmat);
225 }
dgbmatrix & operator=(const dgbmatrix &)
friend _dgematrix i(const dgbmatrix &)
CPPL_INT m
matrix row size
Definition: dgbmatrix.hpp:9
_dgematrix i(const _dgbmatrix &mat)
CPPL_INT kl
lower band width
Definition: dgbmatrix.hpp:11
CPPL_INT ku
upper band width
Definition: dgbmatrix.hpp:12
double * array
1D array to store matrix data
Definition: dgbmatrix.hpp:13
dgbmatrix & operator+=(const dgbmatrix &)
Real Double-precision General Band Matrix Class.
Definition: dgbmatrix.hpp:3
dgbmatrix & operator*=(const dgbmatrix &)
_dgbmatrix operator+(const dgbmatrix &matA, const dgbmatrix &matB)
_dgbmatrix operator-(const dgbmatrix &matA, const dgbmatrix &matB)
void copy(const dgbmatrix &)
(DO NOT USE) Smart-temporary Real Double-precision General Band Matrix Class
Definition: _dgbmatrix.hpp:3
dgbmatrix & zero()
double & operator()(const CPPL_INT &, const CPPL_INT &)
Definition: dgbmatrix-io.hpp:3
CPPL_INT n
matrix column size
Definition: dgbmatrix.hpp:10
dgbmatrix & operator-=(const dgbmatrix &)
_dcovector _(dcovector &vec)
_dgbmatrix operator*(const dgbmatrix &matA, const dgbmatrix &matB)
friend void swap(dgbmatrix &, dgbmatrix &)