1 /*=============================================================================
2 File: matrix.cpp
3 Purpose:
4 Revision: $Id: matrix.cpp,v 1.2 2002/05/13 21:07:45 philosophil Exp $
5 Created by: Philippe Lavoie (3 Oct, 1996)
6 Modified by:
7
8 Copyright notice:
9 Copyright (C) 1996-1998 Philippe Lavoie
10
11 This library is free software; you can redistribute it and/or
12 modify it under the terms of the GNU Library General Public
13 License as published by the Free Software Foundation; either
14 version 2 of the License, or (at your option) any later version.
15
16 This library is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 Library General Public License for more details.
20
21 You should have received a copy of the GNU Library General Public
22 License along with this library; if not, write to the Free
23 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 =============================================================================*/
25
26 #ifndef MATRIX_SOURCES_
27 #define MATRIX_SOURCES_
28
29 #include "matrix_global.h"
30 #include <fstream>
31 #include <string.h>
32
33 #include "matrix.h"
34
35
36 /*!
37 */
38 namespace PLib {
39
40 /*!
41 \brief assignment operator
42
43 \param a the matrix to copy
44
45 \warning the matrix a must have compatible dimensions
46
47 \author Philippe Lavoie
48 \date 24 January 1997
49 */
50 template <class T>
operator =(const Matrix<T> & a)51 Matrix<T>& Matrix<T>::operator=(const Matrix<T> &a){
52 int i;
53
54 if ( this == &a )
55 return *this;
56
57 if ( a.rows() != this->rows() || a.cols() != this->cols() ){
58 this->resize(a.rows(),a.cols()) ;
59 }
60
61 int sze = this->rows()*this->cols() ;
62 T *ptr, *aptr ;
63 ptr = this->m-1 ;
64 aptr = a.m-1 ;
65
66 for (i = sze; i > 0; --i)
67 *(++ptr) = *(++aptr) ;
68
69 this->by_columns = a.by_columns;
70
71 return *this;
72 }
73
74 /*!
75 \brief sets the submatrix \a (s_r,s_c) to \a a .
76
77 \latexonly
78 The matrix can be viewed as
79 \[ \left[\begin{array}{ccc}
80 sub_{0,0} & \vdots & sub_{0,m} \\
81 \cdots & \cdots & \cdots \\
82 sub_{n,0} & \vdots & sub_{n,m} \end{array} \right] \]
83 where each sub matrix is of the size of $a$. This function
84 replaces $sub_{s_r,s_c}$ with $a$.
85 \endlatexonly
86
87 \param sr the row of the submatrix
88 \param sc the column of the submatrix
89 \param a the submatrix to copy from
90
91 \warning Since the size of \a a defines the size of the submatrices,
92 this size must be such that a submatrix located at \a (s_r,s_c)
93 exists.
94
95 \author Philippe Lavoie
96 \date 24 January 1997
97 */
98 template <class T>
submatrix(int sr,int sc,Matrix<T> & a)99 void Matrix<T>::submatrix(int sr, int sc, Matrix<T> &a)
100 {
101 int rwz,coz,i,j;
102
103 if ( this->rows() % a.rows() != 0 || this->cols() % a.cols() != 0 || this->rows() < a.rows() || this->cols() < a.cols() )
104 {
105 #ifdef USE_EXCEPTION
106 throw WrongSize2D(this->rows(),this->cols(),a.rows(),a.cols()) ;
107 #else
108 Error error("Matrix<T>::submatrix");
109 error << "Matrix and submatrix incommensurate" ;
110 error.fatal() ;
111 #endif
112 }
113
114 if ( sr >= this->rows()/a.rows() || sr < 0 || sc >= this->cols()/a.cols() || sc < 0 )
115 {
116 #ifdef USE_EXCEPTION
117 throw OutOfBound2D(sr,sc,0,this->rows()/a.rows()-1,0,this->cols()/a.cols()-1) ;
118 #else
119 Error error("Matrix<T>::submatrix");
120 error << "Submatrix location out of bounds.\nrowblock " << sr << ", " << rows()/a.rows() << " colblock " << sc << ", " << a.cols() << endl ;
121 error.fatal() ;
122 #endif
123 }
124 rwz = sr*a.rows();
125 coz = sc*a.cols();
126
127 #ifdef COLUMN_ORDER
128 for ( i = a.rows()-1; i >= 0; --i )
129 for(j=a.cols()-1;j>=0;--j)
130 elem(i+rwz,j+coz) = a(i,j) ;
131 #else
132 T *ptr, *aptr ;
133 aptr = a.m - 1;
134 for ( i = a.rows()-1; i >= 0; --i )
135 {
136 ptr = &this->m[(i+rwz)*this->cols()+coz]-1 ;
137 for ( j = a.cols(); j > 0; --j)
138 *(++ptr) = *(++aptr) ;
139 }
140 #endif
141 }
142
143 /*!
144 \brief copies a matrix into this matrix starting at index \a (rw,cl)
145
146 \param rw the row to insert the matrix at
147 \param cl the column to insert the matrix at
148 \param a the matrix to insert
149
150 \warning The matrix \a a must fit inside the matrix starting from
151 \a (rw,cl).
152
153 \author Philippe Lavoie
154 \date 24 January 1997
155 */
156 template <class T>
as(int rw,int cl,Matrix<T> & a)157 void Matrix<T>::as(int rw, int cl, Matrix<T>& a)
158 {
159 // Assign matrix a to this matrix at (i,j)
160 int i, j;
161
162 if ( (rw + a.rows()) > this->rows() || ( cl + a.cols()) > this->cols()) {
163 #ifdef USE_EXCEPTION
164 throw MatrixErr();
165 #else
166 Error error("Matrix<T>::as") ;
167 error << "Matrix A will not fit in this Matrix at " << rw << ", " << cl << endl ;
168 error.fatal() ;
169 #endif
170 }
171
172 #ifdef COLUMN_ORDER
173 for(i=0;i<a.rows();++i)
174 for(j=0;j<a.cols();++j)
175 elem(i+rw,j+cl) = a(i,j) ;
176 #else
177 T *pptr,*aptr ;
178 aptr = a.m-1 ;
179 for ( i = 0; i<a.rows(); ++i) {
180 pptr = &this->m[(i+rw)*this->cols()+cl]-1 ;
181 for ( j = 0; j < a.cols(); ++j)
182 *(++pptr) = *(++aptr);
183 }
184 #endif
185 }
186
187 /*!
188 \brief returns the matrix of size \a (nr,nc) starting at \a (rw,cl).
189
190 \latexonly
191 This generates a matrix of size $(n_r,n_c)$ with its first
192 point located at $rw,cl$.
193 \endlatexonly
194
195 \param rw the index of the row
196 \param cl the index of the column
197 \param nr the number of rows() of the matrix to generate
198 \param nc the number of coluns of the matrix to generate
199
200 \return the matrix of size \a (nr,nc) starting at index \a (rw,cl).
201
202 \warning The matrix to return must fit inside the original matrix.
203
204 \author Philippe Lavoie
205 \date 24 January 1997
206 */
207 template <class T>
get(int rw,int cl,int nr,int nc) const208 Matrix<T> Matrix<T>::get(int rw, int cl, int nr, int nc) const
209 {
210 Matrix<T> getmat(nr,nc) ;
211 if ( (rw+nr) > this->rows() || (cl+nc) > this->cols()) {
212 #ifdef USE_EXCEPTION
213 throw MatrixErr();
214 #else
215 Error error("Matrix<T>::get") ;
216 error << "Matrix of size "<< nr << ", " << nc << " not available at " << rw << ", " << cl << endl ;
217 error.fatal() ;
218 #endif
219 }
220
221 int i, j;
222
223 #ifdef COLUMN_ORDER
224 for(i=0;i<nr;++i)
225 for(j=0;j<nc;++j)
226 getmat(i,j) = elem(i+rw,j+cl) ;
227 #else
228 T *pptr,*aptr ;
229 aptr = getmat.m-1;
230 for (i = 0; i < nr; ++i) {
231 pptr = &this->m[(i+rw)*this->cols()+cl]-1 ;
232 for ( j = 0; j < nc; ++j)
233 *(++aptr) = *(++pptr) ;
234 }
235 #endif
236 return getmat ;
237 }
238
239
240
241 /*!
242 \brief Finds the first norm of the matrix
243
244 \return the norm of the matrix
245
246 \author Philippe Lavoie
247 \date 24 January 1997
248 */
249 template <class T>
norm(void)250 double Matrix<T>::norm(void){
251 int i,j ;
252 double sum, maxsum;
253 int init=0 ;
254 T *pptr ;
255 pptr = this->m-1 ;
256 maxsum = 0 ; // Silence the warning message
257 for(i=0;i<this->rows();++i){
258 sum = 0 ;
259 for ( j = 0; j < this->cols(); ++j)
260 sum += *(++pptr) ;
261 if(init)
262 maxsum = (maxsum>sum) ? maxsum : sum;
263 else{
264 maxsum = sum ;
265 init = 1;
266 }
267 }
268 return maxsum;
269 }
270
271
272 /*!
273 \brief sets the diagonal of the matrix to \a a
274
275 Sets the diagonal points of the matrix to \a a. The diagonal
276 points are \a (0,0),(1,1),(2,2),etc.
277
278 \param a the value to set the diagonal to
279
280 \author Philippe Lavoie
281 \date 24 January 1997
282 */
283 template <class T>
diag(const T a)284 void Matrix<T>::diag(const T a)
285 {
286 int i, iend;
287
288 iend = this->rows();
289 if ( iend > this->cols() )
290 iend = this->cols();
291
292 for (i = iend-1; i >=0; --i)
293 this->elem(i,i) = a;
294
295 }
296
297 /*!
298 \brief returns the diagonal of the matrix
299
300 Returns a vector with the component \a [i] being set to the
301 component \a (i,i) of the matrix.
302
303 \return the vector representing the diagonal of the matrix.
304
305 \author Philippe Lavoie
306 \date 24 January 1997
307 */
308 template <class T>
getDiag()309 Vector<T> Matrix<T>::getDiag(){
310 int i, iend;
311 Vector<T> vec(minimum(this->rows(),this->cols())) ;
312 iend = minimum(this->rows(),this->cols());
313 for (i = iend-1; i >=0; --i)
314 vec[i] = this->elem(i,i);
315 return vec ;
316 }
317
318 /*!
319 \brief increase every elements by a double
320
321 \param a the value to increase the elements by
322 \return a reference to itself
323
324 \author Philippe Lavoie
325 \date 1 June 1998
326 */
327 template <class T>
operator +=(double a)328 Matrix<T>& Matrix<T>::operator+=(double a)
329 {
330 T *p1 ;
331 p1 = this->m-1 ;
332 const int size = this->rows()*this->cols() ;
333 for(int i=size; i>0; --i)
334 *(++p1) += a ;
335 return *this ;
336 }
337
338 /*!
339 \brief decrease every elements by a double
340
341 \param a the value to decrease the elements by
342 \return a reference to itself
343
344 \author Philippe Lavoie
345 \date 1 June 1998
346 */
347 template <class T>
operator -=(double a)348 Matrix<T>& Matrix<T>::operator-=(double a)
349 {
350 T *p1 ;
351 p1 = this->m-1 ;
352 const int size = this->rows()*this->cols() ;
353 for(int i=size; i>0; --i)
354 *(++p1) -= a ;
355 return *this ;
356 }
357
358 /*!
359 \brief multiply every elements by a double
360
361 \param a the value to mutiply the elements with
362 \return a reference to itself
363
364 \author Philippe Lavoie
365 \date 1 June 1998
366 */
367 template <class T>
operator *=(double a)368 Matrix<T>& Matrix<T>::operator*=(double a)
369 {
370 T *p1 ;
371 p1 = this->m-1 ;
372 const int size = this->rows()*this->cols() ;
373 for(int i=size; i>0; --i)
374 *(++p1) *= a ;
375 return *this ;
376 }
377
378 /*!
379 \brief divide every elements by a double
380
381 \param a the value to divide the elements with
382 \return a reference to itself
383
384 \author Philippe Lavoie
385 \date 1 June 1998
386 */
387 template <class T>
operator /=(double a)388 Matrix<T>& Matrix<T>::operator/=(double a)
389 {
390 T *p1 ;
391 p1 = this->m-1 ;
392 const int size = this->rows()*this->cols() ;
393 for(int i=size; i>0; --i)
394 *(++p1) /= a ;
395 return *this ;
396 }
397
398 /*!
399 \brief adds a matrix to itself
400
401 \param a the matrix to increment itself with
402 \return a reference to itself
403 \warning the matrix a must have a size compatible with the matrix
404
405 \author Philippe Lavoie
406 \date 24 January 1997
407 */
408 template <class T>
operator +=(const Matrix<T> & a)409 Matrix<T>& Matrix<T>::operator+=(const Matrix<T> &a)
410 {
411 if ( a.rows() != this->rows() || a.cols() != this->cols() )
412 {
413 #ifdef USE_EXCEPTION
414 throw WrongSize2D(this->rows(),this->cols(),a.rows(),a.cols());
415 #else
416 Error error("Matrix<T>::operator+=") ;
417 if ( this->rows() != a.rows() )
418 error << "Matrices are of diferent size, a.rows() = " << rows() << " and b.rows() = " << a.rows() << endl ;
419 if ( this->cols() != a.cols())
420 error << "Matrices are of diferent size, a.cols() = " << cols() << " and b.cols() = " << a.cols() << endl ;
421 error.fatal() ;
422 #endif
423 }
424
425 int i, sze ;
426 T *aptr,*sptr ;
427 aptr = a.m - 1 ;
428 sptr = this->m - 1 ;
429 sze = this->rows()*this->cols() ;
430 for (i = sze; i > 0; --i){
431 *(++sptr) += *(++aptr) ;
432 }
433 return *this ;
434 }
435
436 /*!
437 \brief the addition operator
438
439 \param a the first matrix to add
440 \param b the second matrix to add
441 \return a matrix equal to $a+b$
442
443 \author Philippe Lavoie
444 \date 24 January 1997
445 */
446 template <class T>
operator +(const Matrix<T> & a,const Matrix<T> & b)447 Matrix<T> operator+(const Matrix<T> &a,const Matrix<T> &b)
448 {
449 Matrix<T> sum(a) ;
450 sum += b ;
451 return sum;
452 }
453
454 /*!
455 \brief self substraction
456
457 This will substract the matrix a from the matrix. The result is
458 thus matrix = matrix - a.
459
460 \param a the matrix to substract
461 \return a reference to itself
462
463 \warning The matrix $a$ must be compatible with this matrix
464
465 \author Philippe Lavoie
466 \date 24 January 1997
467 */
468 template <class T>
operator -=(const Matrix<T> & a)469 Matrix<T>& Matrix<T>::operator-=(const Matrix<T> &a)
470 {
471 if ( a.rows() != this->rows() || a.cols() != this->cols() )
472 {
473 #ifdef USE_EXCEPTION
474 throw WrongSize2D(this->rows(),this->cols(),a.rows(),a.cols());
475 #else
476 Error error("Matrix<T>::operator-=") ;
477 if ( this->rows() != a.rows() )
478 error << "Matrices are of diferent size, a.rows() = " << this->rows() << " and b.rows() = " << a.rows() << endl ;
479 if ( this->cols() != a.cols())
480 error << "Matrices are of diferent size, a.cols() = " << this->cols() << " and b.cols() = " << a.cols() << endl ;
481 error.fatal() ;
482 #endif
483 }
484
485 int i, size;
486 T *aptr,*sptr ;
487 aptr = a.m - 1 ;
488 sptr = this->m - 1 ;
489 size = this->rows()*this->cols() ;
490 for (i = size; i > 0; --i){
491 *(++sptr) -= *(++aptr) ;
492 }
493 return *this ;
494 }
495
496
497
498 /*!
499 \brief the substraction operator
500
501 \param a the matrix to substract from
502 \param b the matrix to substract
503 \return a matrix equal to \a a-b
504
505 \warning the matrix must be compatible
506
507 \author Philippe Lavoie
508 \date 24 January 1997
509 */
510 template <class T>
operator -(const Matrix<T> & a,const Matrix<T> & b)511 Matrix<T> operator-(const Matrix<T> &a,const Matrix<T> &b)
512 {
513 Matrix<T> diff(a) ;
514 diff -= b ;
515 return diff;
516 }
517
518 /*!
519 \brief the multiplication operator
520
521 \param a a matrix
522 \param b the matrix to multiply with
523 \return A matrix equal to $a b$
524
525 \warning The matrix must be compatible for the multiplication:
526 a.cols() == b.rows()
527
528 \author Philippe Lavoie
529 \date 24 January 1997
530 */
531 template <class T>
operator *(const Matrix<T> & a,const Matrix<T> & b)532 Matrix<T> operator*(const Matrix<T> &a,const Matrix<T> &b)
533 {
534 if ( a.cols() != b.rows() )
535 {
536 #ifdef USE_EXCEPTION
537 throw WrongSize2D(a.rows(),a.cols(),b.rows(),b.cols()) ;
538 #else
539 Error error("Matrix<T> operator*(Matrix<T>&,Matrix<T>&)");
540 error << "Matrix<T> a * Matrix<T> b incommensurate, a.cols() = " << a.cols() << ", b.rows() = " << b.rows() << endl ;
541 error.fatal() ;
542 #endif
543 }
544
545 int i, j, k, row=a.rows(), col=b.cols(),size = a.cols();
546 Matrix<T> prod(row,col);
547 T zero = (T)0;
548
549 #ifdef COLUMN_ORDER
550 for(i=row-1;i>=0;--i)
551 for(j=size-1;j>=0;--j)
552 if(a(i,j) != zero){
553 for(k=col-1;k>=0; --k)
554 prod(i,k) += a(i,j)* b(j,k) ;
555 }
556 #else
557 T *pptr,*aptr,*bptr ;
558 aptr = a.m ;
559 for (i = 0; i < row; ++i)
560 for (j = 0; j < size; ++j){
561 if ( *aptr != zero )
562 {
563 pptr = prod[i]-1;
564 bptr = b[j]-1;
565 for (k = col; k > 0; --k){
566 *(++pptr) += *aptr * *(++bptr);
567 }
568 }
569 ++aptr;
570 }
571 #endif
572 return prod;
573
574 }
575
576
577 /*!
578 \brief multiplication of a double and a matrix
579
580 \param d the double value to scale the matrix with
581 \param a the matrix to scale with the double value $d$
582 \return A matrix equal to \a d.A
583
584 \author Philippe Lavoie
585 \date 24 January 1997
586 */
587 template <class T>
operator *(const double d,const Matrix<T> & a)588 Matrix<T> operator*(const double d, const Matrix<T> &a)
589 {
590 int i, size=a.rows()*a.cols() ;
591 Matrix<T> b(a.rows(),a.cols());
592
593 T *bptr,*aptr ;
594 bptr = b.m - 1 ;
595 aptr = a.m - 1 ;
596 for (i = size; i > 0; --i)
597 *(++bptr) = (T)(d * (*(++aptr))) ;
598 return b;
599
600 }
601
602 /*!
603 \brief multiplies a matrix with a complex value
604
605 \param d a complex value
606 \param a a matrix
607
608 \return returns a matrix equal to \a d.A
609
610 \author Philippe Lavoie
611 \date 24 January 1997
612 */
613 template <class T>
operator *(const Complex & d,const Matrix<T> & a)614 Matrix<T> operator*(const Complex &d, const Matrix<T> &a)
615 {
616 int i, size=a.rows()*a.cols() ;
617 Matrix<T> b(a.rows(),a.cols());
618
619 T *bptr,*aptr ;
620 bptr = b.m - 1 ;
621 aptr = a.m - 1 ;
622 for (i = size; i > 0; --i)
623 *(++bptr) = (T)real(d) * *(++aptr) ;
624 return b;
625 }
626
627 /*!
628 \brief multiplies a matrix with a vector
629
630 \param a the matrix
631 \param x the vector
632 \return returns a vector representing \a a \a x
633
634 \warning The matrix and the vector must be of compatible sizes
635
636 \author Philippe Lavoie
637 \date 24 January 1997
638 */
639 template <class T>
operator *(const Matrix<T> & a,const Vector<T> & x)640 Vector<T> operator*(const Matrix<T> &a, const Vector<T> &x)
641 {
642 if ( a.cols() != x.size() )
643 {
644 #ifdef USE_EXCEPTION
645 throw WrongSize2D(a.rows(),a.cols(),x.size(),1);
646 #else
647 Error error("Matrix<T> operator*(Matrix<T>& a,Vector<T>& b)");
648 error << "a * b incommensurate, a.cols() = " << a.cols() << ", b.rows() = " << x.size() << endl ;
649 error.fatal() ;
650 #endif
651 }
652
653 int i, k, row = a.rows(), size = a.cols();
654 Vector<T> prod(row);
655 T zero = (T)0;
656
657 T *pptr,*aptr,*xptr ;
658 aptr = a.m - 1 ;
659 pptr = &prod[0] ;
660 for (i = row; i > 0; --i){
661 xptr = x.memory()-1 ;
662 for (k = size, *pptr = zero; k > 0 ; --k)
663 *pptr += *(++aptr) * *(++xptr) ;
664 ++pptr ;
665 }
666
667 return prod;
668 }
669
670
671 /*!
672 \brief the equality operator
673
674 Every elements are compared with each others. If one of them
675 in matrix \a a is not equal to the one in matrix \a b, then
676 the result is negative.
677
678 \param a the first matrix to compare
679 \param b the second matrix to compare
680 \return 1 if equal, 0 otherwise
681
682 \warning The matrices must be of compatible sizes
683
684 \author Philippe Lavoie
685 \date 24 January 1997
686 */
687 template <class T>
operator ==(const Matrix<T> & a,const Matrix<T> & b)688 int operator==(const Matrix<T> &a,const Matrix<T> &b)
689 {
690 if ( a.rows() != b.rows() || a.cols() != b.cols() )
691 {
692 #ifdef USE_EXCEPTION
693 throw WrongSize2D(a.rows(),a.cols(),b.rows(),b.cols()) ;
694 #else
695 Error error("operator==(const Matrix<T>&,const Matrix<T>&)");
696 if ( b.rows() != a.rows() )
697 error << "Matrices are of diferent size, a.rows() = " << a.rows() << " and b.rows() = " << b.rows() << endl ;
698 if ( b.cols() != a.cols())
699 error << "Matrices are of diferent size, a.cols() = " << a.cols() << " and b.cols() = " << b.cols() << endl ;
700 error.fatal() ;
701 #endif
702 }
703
704 int i, j, row = a.rows(), col = a.cols();
705 int l = 1;
706
707 for (i = 0; i < row; ++i)
708 for (j = 0; j < col; ++j)
709 l = l && ( a.elem(i,j) == b.elem(i,j) );
710
711 return l;
712 }
713
714 /*!
715 \brief computes a b - b a
716
717 \param a the a matrix
718 \param b the b matrix
719
720 \warning The a and b matrix must be compatible for the comm operation.
721
722 \author Philippe Lavoie
723 \date 24 January 1997
724 */
725 template <class T>
comm(const Matrix<T> & a,const Matrix<T> & b)726 Matrix<T> comm(const Matrix<T> &a,const Matrix<T> &b)
727 {
728 Matrix<T> r = a * b - b * a;
729
730 return r;
731 }
732
733 /*!
734 \brief The sum of all diagonal elements
735
736 \param a the matrix to trace
737 \return a value representing the sum of all diagonal elements
738
739 \author Philippe Lavoie
740 \date 24 January 1997
741 */
742 template <class T>
trace() const743 T Matrix<T>::trace() const
744 {
745 int size = this->rows();
746 T sum = (T)0;
747
748 if ( size > this->cols() )
749 size = this->cols();
750
751 for (int d = 0; d < size; ++d)
752 sum += this->elem(d,d) ;
753
754 return sum;
755 }
756
757
758 /*!
759 \brief computes the hermitian of the matrix
760
761 This functions returns a matrix for which every elements \a (i,j)
762 correspond to the element \a (j,i) of the original matrix.
763
764 \param a the a matrix
765 \return A matrix corresponding to the hermitian of $a$
766
767 \author Philippe Lavoie
768 \date 24 January 1997
769 */
770 template <class T>
herm() const771 Matrix<T> Matrix<T>::herm() const
772 {
773 int i, j, r = this->cols(), c = this->rows();
774 Matrix<T> adj(r,c);
775
776 for (i = 0; i < r; ++i)
777 for (j = 0; j < c; ++j)
778 adj.elem(i,j) = this->elem(j,i) ;
779
780 return adj;
781
782 }
783
784 /*!
785 \brief returns the matrix flopped
786
787 The flop pixel (i,j) = (i,cols-j-1)
788
789 \return the flop of the matrix
790
791 \author Philippe Lavoie
792 \date 2 May 1999
793 */
794 template <class T>
flop() const795 Matrix<T> Matrix<T>::flop() const
796 {
797 Matrix<T> f(this->rows(),this->cols()) ;
798 for(int i=this->rows()-1;i>=0;--i)
799 for(int j=this->cols()-1;j>=0;--j)
800 {
801 f(i,j) = this->elem(i,this->cols()-j-1);
802 }
803 return f;
804 }
805
806 /*!
807 \brief returns the transpose of the matrix
808
809 \param a the matrix to transpose
810 \return the transpose of the matrix
811
812 \author Philippe Lavoie
813 \date 24 January 1997
814 */
815 template <class T>
transpose() const816 Matrix<T> Matrix<T>::transpose() const
817 {
818 // same as hermitian for real Matrix<T>
819 int i, j;
820 const int& r = this->cols();
821 const int& c = this->rows();
822 Matrix<T> adj(r,c);
823
824 for (i = r-1; i >=0; --i)
825 for (j = c-1; j >=0; --j)
826 adj.elem(i,j) = this->elem(j,i) ;
827
828
829 return adj;
830 }
831
832
833 /*!
834 \brief read a matrix file
835 Reads a matrix file. The format of a file is
836 {\tt rows() cols() data...}, where rows() and cols() are int and data is a vector of the matrix type.
837 \param filename the name of the file to read
838 \return 1 if reading the file was successfull, 0 otherwise
839 \warning
840 \author Philippe Lavoie
841 \date 24 January 1997
842 */
843 template <class T>
read(char * filename)844 int Matrix<T>::read(char* filename) {
845 ifstream fin(filename) ;
846 if(!fin) {
847 this->resize(1,1) ;
848 return 0 ;
849 }
850 int r,c ;
851 char *type ;
852 type = new char[6] ;
853 if(!fin.read(type,sizeof(char)*6)) return 0 ;
854 r = strncmp(type,"matrix",6) ;
855 if(r) return 0 ;
856 if(!fin.read((char*)&r,sizeof(int))) return 0 ;
857 if(!fin.read((char*)&c,sizeof(int))) return 0 ;
858 this->resize(r,c) ;
859 if(!fin.read((char*)this->m,sizeof(T)*r*c)) return 0 ;
860
861 delete []type ;
862 return 1 ;
863 }
864
865 /*!
866 \brief read a raw file containing a matrix of size $(r,c)$
867 Reads a file containing raw data of a matrix of size $(r,c)$.
868 \param filename the name of the file to read
869 r the number of rows()
870 c the number of columns
871 \return 1 if reading the file was successfull, 0 otherwise
872 \warning
873 \author Philippe Lavoie
874 \date 24 January 1997
875 */
876 template <class T>
read(char * filename,int r,int c)877 int Matrix<T>::read(char* filename,int r, int c) {
878 ifstream fin(filename) ;
879 if(!fin) {
880 this->resize(1,1) ;
881 return 0 ;
882 }
883 this->resize(r,c) ;
884 if(!fin.read((char*)this->m,sizeof(T)*r*c)) return 0 ;
885
886 return 1 ;
887 }
888
889
890 /*!
891 \brief write a matrix into a Matrix file
892 Writes a matrix file. The format of the file is
893 {\tt rows() cols() data...}, where rows() and cols() are int and data
894 is a vector of the matrix type.
895 \param filename the name of the file to write to
896 \return
897 \warning 1 if reading the file was successfull, 0 otherwise
898 \author Philippe Lavoie
899 \date 24 January 1997
900 */
901 template <class T>
write(char * filename)902 int Matrix<T>::write(char* filename) {
903 ofstream fout(filename) ;
904 if(!fout)
905 return 0 ;
906 int r,c ;
907 r = this->rows() ; c = this->cols() ;
908 if(!fout.write((char*)&"matrix",sizeof(char)*6)) return 0 ;
909 if(!fout.write((char*)&r,sizeof(int))) return 0 ;
910 if(!fout.write((char*)&c,sizeof(int))) return 0 ;
911 if(!fout.write((char*)this->m,sizeof(T)*r*c)) return 0 ;
912 return 1;
913 }
914
915 /*!
916 \brief write the raw data to a file
917 Writes the raw data to a file. The size information is
918 {\em not} written to the file.
919 \param filename the name of the file to write to
920 \return 0 if an error occurs, 1 otherwise
921 \warning
922 \author Philippe Lavoie
923 \date 24 January 1997
924 */
925 template <class T>
writeRaw(char * filename)926 int Matrix<T>::writeRaw(char* filename) {
927 ofstream fout(filename) ;
928 if(!fout)
929 return 0 ;
930 if(!fout.write((char*)this->m,sizeof(T)*this->rows()*this->cols())) return 0 ;
931 return 1;
932 }
933
934
935 template <class T>
qSort()936 void Matrix<T>::qSort(){
937 #ifdef USE_EXCEPTION
938 throw MatrixErr();
939 #else
940 Error error("Matrix<T>::qSort()");
941 error << "qSort is not defined for that type.\nPlease defined it in your .cpp file!";
942 error.fatal() ;
943 #endif
944 }
945
946
947 }
948
949
950 #endif
951