CPPLapack
 All Classes Files Functions Variables Friends Pages
dgematrix_small-functions.hpp
Go to the documentation of this file.
1 //=============================================================================
2 /*! convert dgematrix_small to dgematrix */
3 template<CPPL_INT m, CPPL_INT n>
5 {CPPL_VERBOSE_REPORT;
6  dgematrix mat(m,n);
7  for(CPPL_INT i=0; i<m; i++){
8  for(CPPL_INT j=0; j<n; j++){
9  mat(i,j) =(*this)(i,j);
10  }
11  }
12  return _(mat);
13 }
14 
15 ///////////////////////////////////////////////////////////////////////////////
16 ///////////////////////////////////////////////////////////////////////////////
17 ///////////////////////////////////////////////////////////////////////////////
18 
19 //=============================================================================
20 /*! operator() */
21 template<CPPL_INT m, CPPL_INT n>
22 inline double& dgematrix_small<m,n>::operator()(const CPPL_INT& i, const CPPL_INT& j)
23 {CPPL_VERBOSE_REPORT;
24 #ifdef CPPL_DEBUG
25  if( i<0 || j<0 || m<=i || n<=j ){
26  ERROR_REPORT;
27  std::cerr << "The required component is out of the matrix size." << std::endl
28  << "Your input is (" << i << "," << j << "), whereas the matrix size is " << m << "x" << n << "." << std::endl;
29  exit(1);
30  }
31 #endif//CPPL_DEBUG
32 
33  return array[i+m*j];
34 }
35 
36 //=============================================================================
37 /*! operator() for const */
38 template<CPPL_INT m, CPPL_INT n>
39 inline double dgematrix_small<m,n>::operator()(const CPPL_INT& i, const CPPL_INT& j) const
40 {CPPL_VERBOSE_REPORT;
41 #ifdef CPPL_DEBUG
42  if( i<0 || j<0 || m<=i || n<=j ){
43  ERROR_REPORT;
44  std::cerr << "The required component is out of the matrix size." << std::endl
45  << "Your input is (" << i << "," << j << "), whereas the matrix size is " << m << "x" << n << "." << std::endl;
46  exit(1);
47  }
48 #endif//CPPL_DEBUG
49 
50  return array[i+m*j];
51 }
52 
53 //=============================================================================
54 /*! set function */
55 template<CPPL_INT m, CPPL_INT n>
56 inline dgematrix_small<m,n>& dgematrix_small<m,n>::set(const CPPL_INT& i, const CPPL_INT& j, const double& v)
57 {CPPL_VERBOSE_REPORT;
58  (*this)(i,j) =v;
59  return *this;
60 }
61 
62 //=============================================================================
63 /*! operator<< */
64 template<CPPL_INT m, CPPL_INT n>
65 inline std::ostream& operator<<(std::ostream& s, const dgematrix_small<m,n>& A)
66 {CPPL_VERBOSE_REPORT;
67  s << std::setiosflags(std::ios::showpos);
68  for(CPPL_INT i=0; i<m; i++){
69  for(CPPL_INT j=0; j<n; j++){
70  s << " " << A(i,j);
71  }
72  s << std::endl;
73  }
74  return s;
75 }
76 
77 //=============================================================================
78 /*! write to file */
79 template<CPPL_INT m, CPPL_INT n>
80 inline void dgematrix_small<m,n>::write(const char* filename) const
81 {CPPL_VERBOSE_REPORT;
82  std::ofstream ofs(filename, std::ios::trunc);
83  ofs.setf(std::cout.flags());
84  ofs.precision(std::cout.precision());
85  ofs.width(std::cout.width());
86  ofs.fill(std::cout.fill());
87  ofs << "#dgematrix" << " " << m << " " << n << std::endl;
88  for(CPPL_INT i=0; i<m; i++){
89  for(CPPL_INT j=0; j<n; j++){
90  ofs << (*this)(i,j) << " ";
91  }
92  ofs << std::endl;
93  }
94  ofs.close();
95 }
96 
97 //=============================================================================
98 /*! read from file */
99 template<CPPL_INT m, CPPL_INT n>
100 inline void dgematrix_small<m,n>::read(const char* filename)
101 {CPPL_VERBOSE_REPORT;
102  std::ifstream s( filename );
103  if(!s){
104  ERROR_REPORT;
105  std::cerr << "The file \"" << filename << "\" can not be opened." << std::endl;
106  exit(1);
107  }
108 
109  std::string id;
110  s >> id;
111  if( id != "dgematrix" && id != "#dgematrix" ){
112  ERROR_REPORT;
113  std::cerr << "The type name of the file \"" << filename << "\" is not dgematrix." << std::endl
114  << "Its type name was " << id << " ." << std::endl;
115  exit(1);
116  }
117 
118  CPPL_INT _m, _n;
119  s >> _m >> _n;
120  if(m!=_m || n!=_n){
121  ERROR_REPORT;
122  std::cerr << "Matrix size is invalid." << std::endl;
123  exit(1);
124  }
125  for(CPPL_INT i=0; i<m; i++){
126  for(CPPL_INT j=0; j<n; j++ ){
127  s >> operator()(i,j);
128  }
129  }
130  if(s.eof()){
131  ERROR_REPORT;
132  std::cerr << "There is something is wrong with the file \"" << filename << "\"." << std::endl
133  << "Most likely, there is not enough data components, or a linefeed code or space code is missing at the end of the last line." << std::endl;
134  exit(1);
135  }
136 
137  s >> id;
138  if(!s.eof()){
139  ERROR_REPORT;
140  std::cerr << "There is something is wrong with the file \"" << filename << "\"." << std::endl
141  << "Most likely, there are extra data components." << std::endl;
142  exit(1);
143  }
144 
145  s.close();
146 }
147 
148 ///////////////////////////////////////////////////////////////////////////////
149 ///////////////////////////////////////////////////////////////////////////////
150 ///////////////////////////////////////////////////////////////////////////////
151 ///////////////////////////////////////////////////////////////////////////////
152 ///////////////////////////////////////////////////////////////////////////////
153 ///////////////////////////////////////////////////////////////////////////////
154 
155 //=============================================================================
156 /*! return transposed matrix */
157 template<CPPL_INT m, CPPL_INT n>
159 {CPPL_VERBOSE_REPORT;
161  for(CPPL_INT i=0; i<m; i++){
162  for(CPPL_INT j=0; j<n; j++){
163  X(j,i) =A(i,j);
164  }
165  }
166  return X;
167 }
168 
169 //=============================================================================
170 /*! return its trace */
171 template<CPPL_INT m, CPPL_INT n>
172 inline double trace(const dgematrix_small<m,n>& A)
173 {CPPL_VERBOSE_REPORT;
174  double trace =0.;
175 
176  const CPPL_INT imax =std::min(m,n);
177  for(CPPL_INT i=0; i<imax; i++){
178  trace +=A(i,i);
179  }
180 
181  return trace;
182 }
183 
184 //=============================================================================
185 /*! find index of the maximum component */
186 template<CPPL_INT m, CPPL_INT n>
187 inline void idamax(CPPL_INT& I, CPPL_INT& J, const dgematrix_small<m,n>& A)
188 {CPPL_VERBOSE_REPORT;
189  double max(-1.);
190  for(CPPL_INT i=0; i<m; i++){
191  for(CPPL_INT j=0; j<n; j++){
192  if( max<fabs(A(i,j)) ){
193  I=i;
194  J=j;
195  max =fabs(A(i,j));
196  }
197  }
198  }
199  return;
200 }
201 
202 //=============================================================================
203 /*! return the maximum component */
204 template<CPPL_INT m, CPPL_INT n>
205 inline double damax(const dgematrix_small<m,n>& A)
206 {CPPL_VERBOSE_REPORT;
207  CPPL_INT i(0), j(0);
208  idamax(i,j,A);
209  return A(i,j);
210 }
211 
212 ///////////////////////////////////////////////////////////////////////////////
213 ///////////////////////////////////////////////////////////////////////////////
214 ///////////////////////////////////////////////////////////////////////////////
215 
216 //=============================================================================
217 /*! zero */
218 template<CPPL_INT m, CPPL_INT n>
220 {CPPL_VERBOSE_REPORT;
221  for(CPPL_INT k=0; k<m*n; k++){ array[k]=0.; }
222  return *this;
223 }
224 
225 //=============================================================================
226 /*! identity */
227 template<CPPL_INT m, CPPL_INT n>
229 {CPPL_VERBOSE_REPORT;
230  zero();
231 
232  const CPPL_INT kmax =std::min(m,n);
233  for(CPPL_INT k=0; k<kmax; k++){
234  (*this)(k,k)=1.;
235  }
236 
237  return *this;
238 }
239 
240 //=============================================================================
241 /*! return the j-th column vector */
242 template<CPPL_INT m, CPPL_INT n>
243 inline dcovector_small<m> dgematrix_small<m,n>::col(const CPPL_INT& j) const
244 {CPPL_VERBOSE_REPORT;
245  dcovector_small<m> vec;
246  for(CPPL_INT i=0; i<m; i++){ vec(i)=(*this)(i,j); }
247  return vec;
248 }
249 
250 //=============================================================================
251 /*! return the i-th row vector */
252 template<CPPL_INT m, CPPL_INT n>
253 inline drovector_small<n> dgematrix_small<m,n>::row(const CPPL_INT& i) const
254 {CPPL_VERBOSE_REPORT;
255  drovector_small<n> vec;
256  for(CPPL_INT j=0; j<n; j++){ vec(j)=(*this)(i,j); }
257  return vec;
258 }
259 
260 ///////////////////////////////////////////////////////////////////////////////
261 ///////////////////////////////////////////////////////////////////////////////
262 ///////////////////////////////////////////////////////////////////////////////
263 
264 //=============================================================================
265 /*! dgematrix_small+=dgematrix_small operator */
266 template<CPPL_INT m, CPPL_INT n>
268 {CPPL_VERBOSE_REPORT;
269  for(CPPL_INT k=0; k<m*n; k++){
270  A.array[k] +=B.array[k];
271  }
272  return A;
273 }
274 
275 //=============================================================================
276 /*! dgematrix_small-=dgematrix_small operator */
277 template<CPPL_INT m, CPPL_INT n>
279 {CPPL_VERBOSE_REPORT;
280  for(CPPL_INT k=0; k<m*n; k++){
281  A.array[k] -=B.array[k];
282  }
283  return A;
284 }
285 
286 //=============================================================================
287 /*! dgematrix_small*=dgematrix_small operator */
288 template<CPPL_INT m, CPPL_INT l, CPPL_INT n>
290 {CPPL_VERBOSE_REPORT;
292  X.zero();
293  for(CPPL_INT i=0; i<m; i++){
294  for(CPPL_INT j=0; j<n; j++){
295  for(CPPL_INT k=0; k<l; k++){
296  X(i,j) += A(i,k)*B(k,j);
297  }
298  }
299  }
300  A =X;
301  return A;
302 }
303 
304 //=============================================================================
305 /*! dgematrix_small*=double operator */
306 template<CPPL_INT m, CPPL_INT n>
308 {CPPL_VERBOSE_REPORT;
309  for(CPPL_INT k=0; k<m*n; k++){
310  A.array[k] *=d;
311  }
312  return A;
313 }
314 
315 //=============================================================================
316 /*! dgematrix_small/double operator */
317 template<CPPL_INT m, CPPL_INT n>
319 {CPPL_VERBOSE_REPORT;
320  for(CPPL_INT k=0; k<m*n; k++){
321  A.array[k] /=d;
322  }
323  return A;
324 }
325 
326 ///////////////////////////////////////////////////////////////////////////////
327 ///////////////////////////////////////////////////////////////////////////////
328 ///////////////////////////////////////////////////////////////////////////////
329 
330 //=============================================================================
331 /*! unary + operator */
332 template<CPPL_INT m, CPPL_INT n>
334 {CPPL_VERBOSE_REPORT;
335  return A;
336 }
337 
338 //=============================================================================
339 /*! unary - operator */
340 template<CPPL_INT m, CPPL_INT n>
342 {CPPL_VERBOSE_REPORT;
344  for(CPPL_INT i=0; i<m; i++){
345  for(CPPL_INT j=0; j<n; j++){
346  X(i,j) =-A(i,j);
347  }
348  }
349  return X;
350 }
351 
352 ///////////////////////////////////////////////////////////////////////////////
353 ///////////////////////////////////////////////////////////////////////////////
354 ///////////////////////////////////////////////////////////////////////////////
355 
356 //=============================================================================
357 /*! dgematrix_small+dgematrix_small operator */
358 template<CPPL_INT m, CPPL_INT n>
360 {CPPL_VERBOSE_REPORT;
362  for(int i=0; i<m; i++){
363  for(int j=0; j<n; j++){
364  C(i,j) =A(i,j)+B(i,j);
365  }
366  }
367  return C;
368 }
369 
370 //=============================================================================
371 /*! dgematrix_small+dsymatrix_small operator */
372 template<CPPL_INT n>
374 {CPPL_VERBOSE_REPORT;
376  for(CPPL_INT i=0; i<n; i++){
377  for(CPPL_INT j=0; j<i; j++){
378  X(i,j) =A(i,j)+B(i,j);
379  }
380  for(CPPL_INT j=i; j<n; j++){
381  X(i,j) =A(i,j)+B(j,i);
382  }
383  }
384  return X;
385 }
386 
387 //=============================================================================
388 /*! dgematrix_small-dgematrix_small operator */
389 template<CPPL_INT m, CPPL_INT n>
391 {CPPL_VERBOSE_REPORT;
393  for(int i=0; i<m; i++){
394  for(int j=0; j<n; j++){
395  C(i,j)=A(i,j)-B(i,j);
396  }
397  }
398  return C;
399 }
400 
401 //=============================================================================
402 /*! dgematrix_small-dsymatrix_small operator */
403 template<CPPL_INT n>
405 {CPPL_VERBOSE_REPORT;
407  for(CPPL_INT i=0; i<n; i++){
408  for(CPPL_INT j=0; j<=i; j++){ X(i,j)=A(i,j)-B(i,j); }
409  for(CPPL_INT j=i+1; j<n; j++){ X(i,j)=A(i,j)-B(j,i); }
410  }
411  return X;
412 }
413 
414 //=============================================================================
415 /*! dgematrix_small*dcovector_small operator */
416 template<CPPL_INT m, CPPL_INT n>
418 {CPPL_VERBOSE_REPORT;
420  C.zero();
421  for(CPPL_INT i=0; i<m; i++){
422  for(CPPL_INT j=0; j<n; j++){
423  C(i) +=A(i,j)*B(j);
424  }
425  }
426  return C;
427 }
428 
429 //=============================================================================
430 /*! dgematrix_small*dgematrix_small operator */
431 template<CPPL_INT m, CPPL_INT l, CPPL_INT n>
433 {CPPL_VERBOSE_REPORT;
435  C.zero();
436  for(int i=0; i<m; i++){
437  for(int j=0; j<n; j++){
438  for(int k=0; k<l; k++){
439  C(i,j) +=A(i,k)*B(k,j);
440  }
441  }
442  }
443  return C;
444 }
445 
446 //=============================================================================
447 /*! dgematrix_small*dsymatrix_small operator */
448 template<CPPL_INT m, CPPL_INT n>
450 {CPPL_VERBOSE_REPORT;
452  X.zero();
453  for(CPPL_INT i=0; i<m; i++){
454  for(CPPL_INT j=0; j<n; j++){
455  for(CPPL_INT k=0; k<j; k++){ X(i,j)+=A(i,k)*B(j,k); }
456  for(CPPL_INT k=j; k<n; k++){ X(i,j)+=A(i,k)*B(k,j); }
457  }
458  }
459  return X;
460 }
461 
462 //=============================================================================
463 /*! dgematrix_small*double operator */
464 template<CPPL_INT m, CPPL_INT n>
465 inline dgematrix_small<m,n> operator*(const dgematrix_small<m,n>& A, const double& v)
466 {CPPL_VERBOSE_REPORT;
468  for(CPPL_INT i=0; i<m; i++){
469  for(CPPL_INT j=0; j<n; j++){
470  C(i,j) =A(i,j)*v;
471  }
472  }
473  return C;
474 }
475 
476 //=============================================================================
477 /*! dgematrix_small/double operator */
478 template<CPPL_INT m, CPPL_INT n>
479 inline dgematrix_small<m,n> operator/(const dgematrix_small<m,n>& A, const double& v)
480 {CPPL_VERBOSE_REPORT;
482  for(CPPL_INT i=0; i<m; i++){
483  for(CPPL_INT j=0; j<n; j++){
484  C(i,j) =A(i,j)/v;
485  }
486  }
487  return C;
488 }
489 
490 ///////////////////////////////////////////////////////////////////////////////
491 ///////////////////////////////////////////////////////////////////////////////
492 ///////////////////////////////////////////////////////////////////////////////
493 
494 //=============================================================================
495 /*! Hadamard product */
496 template<CPPL_INT m, CPPL_INT n>
498 {CPPL_VERBOSE_REPORT;
500  for(CPPL_INT i=0; i<m; i++){
501  for(CPPL_INT j=0; j<n; j++){
502  C(i,j) =A(i,j)*B(i,j);
503  }
504  }
505  return C;
506 }
507 
508 //=============================================================================
509 /*! Hadamard product */
510 template<CPPL_INT n>
512 {CPPL_VERBOSE_REPORT;
514  for(CPPL_INT i=0; i<n; i++){
515  for(CPPL_INT j=0; j<=i; j++){
516  C(i,j) =A(i,j)*B(i,j);
517  }
518  for(CPPL_INT j=i+1; j<n; j++){
519  C(i,j) =A(i,j)*B(j,i);
520  }
521  }
522  return C;
523 }
dcovector_small< m > operator*(const dgematrix_small< m, n > &A, const dcovector_small< n > &B)
void write(const char *filename) const
dgematrix_small< m, n > & identity()
dgematrix_small< m, n > & set(const CPPL_INT &i, const CPPL_INT &j, const double &v)
double & operator()(const CPPL_INT &i, const CPPL_INT &j)
dcovector_small< l > & zero()
Samll Real Double-precision Symmetric Matrix Class.
void read(const char *filename)
Samll Real Double-precision General Dence Matrix Class.
_dgematrix i(const _dgbmatrix &mat)
dgematrix_small< m, n > & zero()
dgematrix_small< m, n > & operator*=(dgematrix_small< m, l > &A, const dgematrix_small< l, n > &B)
Real Double-precision General Dence Matrix Class.
Definition: dgematrix.hpp:3
dgematrix_small< m, n > operator-(const dgematrix_small< m, n > &A)
Samll Real Double-precision Row Vector Class.
dgematrix_small< m, n > & operator-=(dgematrix_small< m, n > &A, const dgematrix_small< m, n > &B)
const dgematrix_small< m, n > & operator+(const dgematrix_small< m, n > &A)
dcovector_small< m > col(const CPPL_INT &j) const
dgematrix_small< m, n > & operator+=(dgematrix_small< m, n > &A, const dgematrix_small< m, n > &B)
(DO NOT USE) Smart-temporary Real Double-precision General Dence Matrix Class
Definition: _dgematrix.hpp:3
dgematrix_small< m, n > & operator/=(dgematrix_small< m, n > &A, const double &d)
double damax(const dgematrix_small< m, n > &A)
dgematrix_small< m, n > operator/(const dgematrix_small< m, n > &A, const double &v)
double trace(const dgematrix_small< m, n > &A)
void idamax(CPPL_INT &I, CPPL_INT &J, const dgematrix_small< m, n > &A)
Samll Real Double-precision Column Vector Class.
dgematrix_small< n, m > t(const dgematrix_small< m, n > &A)
double array[m *n]
1D array to store vector data
drovector_small< n > row(const CPPL_INT &i) const
_dcovector _(dcovector &vec)
_dgematrix to_dgematrix() const
dgematrix_small< m, n > hadamard(const dgematrix_small< m, n > &A, const dgematrix_small< m, n > &B)