cmplxmat.h 8.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361
  1. #ifndef _cmplxmat__h_
  2. #define _cmplxmat__h_
  3. #include <stdlib.h>
  4. #include <math.h>
  5. #include "cmplxvec.h"
  6. #include "vector.h"
  7. #include "matrix.h"
  8. class CmplxMatrix
  9. {
  10. friend class Matrix;
  11. friend ostream& operator<<(ostream&, const CmplxMatrix&);
  12. friend istream& operator>>(istream&, CmplxMatrix&);
  13. public:
  14. CmplxMatrix(int=0, int=0, double=0.0);
  15. CmplxMatrix(int, int, Complex);
  16. CmplxMatrix(const CmplxMatrix&);
  17. CmplxMatrix(const Matrix&);
  18. CmplxMatrix(const CmplxVector&);
  19. CmplxMatrix(const Vector&);
  20. CmplxMatrix& operator=(const CmplxMatrix&);
  21. CmplxMatrix& operator=(const Matrix&);
  22. ~CmplxMatrix();
  23. void free();
  24. int dim_i() const { return d_i; }
  25. int dim_j() const { return d_j; }
  26. CmplxVector& row(int) const;
  27. CmplxVector col(int i) const;
  28. CmplxMatrix trans() const;
  29. Matrix real() const;
  30. Matrix imag() const;
  31. Matrix abs() const;
  32. CmplxMatrix conjg() const;
  33. CmplxMatrix& resize(int,int);
  34. operator CmplxVector() const;
  35. int operator==(const CmplxMatrix&) const;
  36. int operator!=(const CmplxMatrix& x) const { return !(*this == x); }
  37. CmplxVector& operator[](int i) const { return row(i); }
  38. Complex& operator()(int, int);
  39. Complex& operator()(int, int) const;
  40. CmplxMatrix operator+(const CmplxMatrix&);
  41. CmplxMatrix operator-(const CmplxMatrix&);
  42. CmplxMatrix& operator+=(const CmplxMatrix&);
  43. CmplxMatrix& operator-=(const CmplxMatrix&);
  44. CmplxMatrix operator+(const Matrix&);
  45. CmplxMatrix operator-(const Matrix&);
  46. CmplxMatrix operator-();
  47. CmplxMatrix operator*(double);
  48. CmplxMatrix operator*(const Complex&);
  49. CmplxMatrix operator/(double);
  50. CmplxMatrix operator/(const Complex&);
  51. CmplxMatrix operator*(const CmplxMatrix&);
  52. CmplxMatrix operator*(const Matrix&);
  53. CmplxVector operator*(CmplxVector& w) { return CmplxVector(*this * CmplxMatrix(w)); }
  54. //CmplxVector operator*(CmplxVector& w) const { return CmplxVector(*this * CmplxMatrix(w)); }
  55. //CmplxVector operator*(const CmplxVector& w) { return CmplxVector(*this * CmplxMatrix(w)); }
  56. //CmplxVector operator*(const CmplxVector& w) const { return CmplxVector(*this * CmplxMatrix(w)); }
  57. CmplxVector operator*(Vector& w) { return CmplxVector(*this * CmplxMatrix(w)); }
  58. //CmplxVector operator*(Vector& w) const { return CmplxVector(*this * CmplxMatrix(w)); }
  59. CmplxVector operator*(const Vector& w) { return CmplxVector(*this * CmplxMatrix(w)); }
  60. //CmplxVector operator*(const Vector& w) const { return CmplxVector(*this * CmplxMatrix(w)); }
  61. private:
  62. CmplxVector** v;
  63. int d_i;
  64. int d_j;
  65. void flip_rows(int,int);
  66. void check_dimensions(const CmplxMatrix&) const;
  67. Complex& elem(int i, int j) const { return v[i]->v[j]; }
  68. };
  69. inline CmplxMatrix& CmplxMatrix::operator=(const CmplxMatrix& mat) {
  70. int i, j;
  71. if (d_i != mat.d_i || d_j != mat.d_j) {
  72. for(i = 0; i < d_i; i++) delete v[i];
  73. delete v;
  74. d_i = mat.d_i;
  75. d_j = mat.d_j;
  76. v = new CmplxVector*[d_i];
  77. for(i = 0; i < d_i; i++) v[i] = new CmplxVector(d_j);
  78. }
  79. for(i = 0; i < d_i; i++)
  80. for(j = 0; j < d_j; j++) elem(i,j) = mat.elem(i,j);
  81. return (*this);
  82. }
  83. inline CmplxMatrix& CmplxMatrix::operator=(const Matrix& mat) {
  84. int i, j;
  85. CmplxMatrix tmp(mat);
  86. if (d_i != mat.dim_i() || d_j != mat.dim_j()) {
  87. for(i = 0; i < d_i; i++) delete v[i];
  88. delete v;
  89. d_i = mat.dim_i();
  90. d_j = mat.dim_j();
  91. v = new CmplxVector*[d_i];
  92. for(i = 0; i < d_i; i++) v[i] = new CmplxVector(d_j);
  93. }
  94. for(i = 0; i < d_i; i++)
  95. for(j = 0; j < d_j; j++) elem(i,j) = tmp.elem(i,j);
  96. return (*this);
  97. }
  98. inline int CmplxMatrix::operator==(const CmplxMatrix& x) const {
  99. int i, j;
  100. if (d_i != x.d_i || d_j != x.d_j) return (0);
  101. for(i = 0; i < d_i; i++)
  102. for(j = 0; j < d_j; j++)
  103. if (elem(i,j) != x.elem(i,j)) return (0);
  104. return (1);
  105. }
  106. inline CmplxVector& CmplxMatrix::row(int i) const {
  107. if (i < 0 || i >= d_i) {
  108. cerr << "complex matrix: row index out of range" << endl;
  109. exit(1);
  110. }
  111. return (*v[i]);
  112. }
  113. inline Complex& CmplxMatrix::operator()(int i, int j) {
  114. if (i < 0 || i >= d_i) {
  115. cerr << "complex matrix: row index out of range" << endl;
  116. exit(1);
  117. }
  118. if (j < 0 || j >= d_j) {
  119. cerr << "complex matrix: col index out of range" << endl;
  120. exit(1);
  121. }
  122. return (elem(i,j));
  123. }
  124. inline Complex& CmplxMatrix::operator()(int i, int j) const {
  125. if (i < 0 || i >= d_i) {
  126. cerr << "complex matrix: row index out of range" << endl;
  127. exit(1);
  128. }
  129. if (j < 0 || j >= d_j) {
  130. cerr << "complex matrix: col index out of range" << endl;
  131. exit(1);
  132. }
  133. return (elem(i,j));
  134. }
  135. inline CmplxVector CmplxMatrix::col(int i) const {
  136. if (i < 0 || i >= d_j) {
  137. cerr << "complex matrix: col index out of range" << endl;
  138. exit(1);
  139. }
  140. CmplxVector result(d_i);
  141. int j = d_i;
  142. while (j--) result.v[j] = elem(j,i);
  143. return (result);
  144. }
  145. inline CmplxMatrix::operator CmplxVector() const {
  146. if (d_j != 1) {
  147. cerr << "error: cannot make complex vector from complex matrix" << endl;
  148. exit(1);
  149. }
  150. return (col(0));
  151. }
  152. inline CmplxMatrix CmplxMatrix::operator+(const CmplxMatrix& mat) {
  153. int i, j;
  154. check_dimensions(mat);
  155. CmplxMatrix result(d_i,d_j);
  156. for(i = 0; i < d_i; i++)
  157. for(j = 0; j < d_j; j++)
  158. result.elem(i,j) = elem(i,j) + mat.elem(i,j);
  159. return (result);
  160. }
  161. inline CmplxMatrix CmplxMatrix::operator+(const Matrix& mat) {
  162. int i, j;
  163. if (d_i != mat.dim_i() || d_j != mat.dim_j()) {
  164. cerr << "incompatible complex matrix types. +" << endl;
  165. exit(1);
  166. }
  167. CmplxMatrix result(d_i,d_j);
  168. CmplxMatrix tmp(mat);
  169. for(i = 0; i < d_i; i++)
  170. for(j = 0; j < d_j; j++)
  171. result.elem(i,j) = elem(i,j) + tmp.elem(i,j);
  172. return (result);
  173. }
  174. inline CmplxMatrix CmplxMatrix::operator-(const CmplxMatrix& mat) {
  175. int i, j;
  176. check_dimensions(mat);
  177. CmplxMatrix result(d_i,d_j);
  178. for(i = 0; i < d_i; i++)
  179. for(j = 0; j < d_j; j++)
  180. result.elem(i,j) = elem(i,j) - mat.elem(i,j);
  181. return (result);
  182. }
  183. inline CmplxMatrix CmplxMatrix::operator-(const Matrix& mat) {
  184. int i, j;
  185. if (d_i != mat.dim_i() || d_j != mat.dim_j()) {
  186. cerr << "incompatible complex matrix types in -." << endl;
  187. exit(1);
  188. }
  189. CmplxMatrix result(d_i,d_j);
  190. CmplxMatrix tmp(mat);
  191. for(i = 0; i < d_i; i++)
  192. for(j = 0; j < d_j; j++)
  193. result.elem(i,j) = elem(i,j) - tmp.elem(i,j);
  194. return (result);
  195. }
  196. inline CmplxMatrix& CmplxMatrix::operator+=(const CmplxMatrix& mat) {
  197. int i, j;
  198. check_dimensions(mat);
  199. for(i = 0; i < d_i; i++)
  200. for(j = 0; j < d_j; j++)
  201. elem(i,j) += mat.elem(i,j);
  202. return (*this);
  203. }
  204. inline CmplxMatrix& CmplxMatrix::operator-=(const CmplxMatrix& mat) {
  205. int i, j;
  206. check_dimensions(mat);
  207. for(i = 0; i < d_i; i++)
  208. for(j = 0; j < d_j; j++)
  209. elem(i,j) -= mat.elem(i,j);
  210. return (*this);
  211. }
  212. inline CmplxMatrix CmplxMatrix::operator-() {
  213. int i, j;
  214. CmplxMatrix result(d_i,d_j);
  215. for(i = 0; i < d_i; i++)
  216. for(j = 0; j < d_j; j++)
  217. result.elem(i,j) = -elem(i,j);
  218. return (result);
  219. }
  220. inline CmplxMatrix CmplxMatrix::operator*(double f) {
  221. int i, j;
  222. CmplxMatrix result(d_i,d_j);
  223. for(i = 0; i < d_i; i++)
  224. for(j = 0; j < d_j; j++)
  225. result.elem(i,j) = elem(i,j) * f;
  226. return (result);
  227. }
  228. inline CmplxMatrix CmplxMatrix::operator*(const Complex& f) {
  229. int i, j;
  230. CmplxMatrix result(d_i,d_j);
  231. for(i = 0; i < d_i; i++)
  232. for(j = 0; j < d_j; j++)
  233. result.elem(i,j) = elem(i,j) * f;
  234. return (result);
  235. }
  236. inline CmplxMatrix CmplxMatrix::operator*(const Matrix& mat) {
  237. if (d_j != mat.dim_i()) {
  238. cerr << "complex matrix: incompatible matrix types" << endl;
  239. exit(1);
  240. }
  241. CmplxMatrix result(d_i, mat.dim_j());
  242. int i, j;
  243. for(i = 0;i < mat.dim_j(); i++)
  244. for(j = 0; j < d_i; j++) result.elem(j,i) = *v[j] * mat.col(i);
  245. return (result);
  246. }
  247. inline CmplxMatrix CmplxMatrix::operator*(const CmplxMatrix& mat) {
  248. if (d_j != mat.d_i) {
  249. cerr << "complex matrix: incompatible matrix types" << endl;
  250. exit(1);
  251. }
  252. CmplxMatrix result(d_i, mat.d_j);
  253. int i,j;
  254. for(i = 0; i < mat.d_j; i++)
  255. for(j = 0; j < d_i; j++) result.elem(j,i) = *v[j] * mat.col(i);
  256. return (result);
  257. }
  258. inline CmplxMatrix CmplxMatrix::operator/(double a) {
  259. if (a==0) {
  260. cerr << "complex matrix: divided by zero" << endl;
  261. exit(1);
  262. }
  263. a = 1. / a;
  264. return (CmplxMatrix(*this * a));
  265. }
  266. inline CmplxMatrix CmplxMatrix::operator/(const Complex& a) {
  267. if (a==0) {
  268. cerr << "complex matrix: divided by zero" << endl;
  269. exit(1);
  270. }
  271. Complex tmp = 1. / a;
  272. return (CmplxMatrix(*this * a));
  273. }
  274. inline void CmplxMatrix::free()
  275. {
  276. if (v) {
  277. while(d_i--) delete v[d_i];
  278. delete v;
  279. }
  280. }
  281. inline void CmplxMatrix::check_dimensions(const CmplxMatrix& mat) const {
  282. if (d_i != mat.d_i || d_j != mat.d_j) {
  283. cerr << "incompatible complex matrix types." << endl;
  284. exit(1);
  285. }
  286. }
  287. inline void CmplxMatrix::flip_rows(int i,int j) {
  288. CmplxVector* p = v[i];
  289. v[i] = v[j];
  290. v[j] = p;
  291. }
  292. inline void Print(const CmplxMatrix& m, ostream& out=cout) { out << m; }
  293. inline void Read(CmplxMatrix& m, istream& in=cin) { in >> m; }
  294. #endif