Line data Source code
1 : //# SquareMatrix.cc: Fast Square Matrix class with fixed (templated) size
2 : //# Copyright (C) 1996,1998,1999,2001,2002
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 :
26 : #ifndef SCIMATH_SQUAREMATRIX_TCC
27 : #define SCIMATH_SQUAREMATRIX_TCC
28 :
29 : //# Includes
30 : #include <casacore/scimath/Mathematics/SquareMatrix.h>
31 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
32 : #include <casacore/casa/BasicSL/Complex.h>
33 : #include <casacore/casa/IO/ArrayIO.h>
34 : #include <casacore/casa/Exceptions/Error.h>
35 : #include <casacore/casa/iostream.h>
36 :
37 : namespace casacore { //# NAMESPACE CASACORE - BEGIN
38 :
39 : template <class T, Int n>
40 : SquareMatrix<T,n>&
41 17672 : SquareMatrix<T,n>::operator=(const SquareMatrix<T,n>& m) {
42 17672 : type_p=m.type_p;
43 17672 : switch (type_p) {
44 3380 : case ScalarId: a_p[0][0]=m.a_p[0][0]; break;
45 0 : case Diagonal: {
46 0 : for (Int i=0; i<n; i++) a_p[i][i]=m.a_p[i][i];
47 0 : break;
48 : }
49 14292 : case General: {
50 14292 : const T* pm=&m.a_p[0][0];
51 14292 : T* pa_p=&a_p[0][0];
52 71460 : for (Int i=0; i<n*n; i++) *pa_p++=*pm++;
53 : }
54 : }
55 17672 : return *this;
56 : }
57 : //# not accepted out of line by native compiler- moved inline
58 : template <class T, Int n>
59 : SquareMatrix<T,n>&
60 : SquareMatrix<T,n>::operator=(const Vector<T>& v) {
61 : for (Int i=0; i<n; i++) a_p[i][i]=v(i);
62 : type_p=Diagonal;
63 : return *this;
64 : }
65 : template <class T, Int n>
66 : SquareMatrix<T,n>&
67 0 : SquareMatrix<T,n>::operator=(const Matrix<T>& m) {
68 0 : for (Int i=0; i<n; i++)
69 0 : for (Int j=0; j<n; j++) a_p[i][j]=m(i,j);
70 0 : type_p=General;
71 0 : return *this;
72 : }
73 : template <class T, Int n>
74 : SquareMatrix<T,n>&
75 : SquareMatrix<T,n>::operator+=(const SquareMatrix<T,n>& other) {
76 : switch (type_p) {
77 : case ScalarId:
78 : switch (other.type_p) {
79 : case ScalarId: {
80 : a_p[0][0]+=other.a_p[0][0];
81 : return *this;
82 : }
83 : case Diagonal: {
84 : T tmp=a_p[0][0];
85 : for (Int i=0; i<n; i++) a_p[i][i]=tmp+other.a_p[i][i];
86 : type_p=Diagonal;
87 : return *this;
88 : }
89 : case General: {
90 : T tmp=a_p[0][0];
91 : for (Int i=0; i<n; i++) {
92 : a_p[i][i]=tmp+other.a_p[i][i];
93 : for (Int j=0; j<n; j++)
94 : if (i!=j) a_p[i][j]=other.a_p[i][j];
95 : }
96 : type_p=General;
97 : return *this;
98 : }
99 : }
100 : case Diagonal:
101 : switch (other.type_p) {
102 : case ScalarId: {
103 : for (Int i=0; i<n; i++) a_p[i][i]+=other.a_p[0][0];
104 : return *this;
105 : }
106 : case Diagonal: {
107 : for (Int i=0; i<n; i++) a_p[i][i]+=other.a_p[i][i];
108 : return *this;
109 : }
110 : case General: {
111 : for (Int i=0; i<n; i++) {
112 : a_p[i][i]+=other.a_p[i][i];
113 : for (Int j=0; j<n; j++)
114 : if (i!=j) a_p[i][j]=other.a_p[i][j];
115 : }
116 : type_p=General;
117 : return *this;
118 : }
119 : }
120 : // case General:
121 : default:
122 : switch (other.type_p) {
123 : case ScalarId: {
124 : for (Int i=0; i<n; i++) a_p[i][i]+=other.a_p[0][0];
125 : return *this;
126 : }
127 : case Diagonal: {
128 : for (Int i=0; i<n; i++) a_p[i][i]+=other.a_p[i][i];
129 : return *this;
130 : }
131 : case General: {
132 : const T* po=&other.a_p[0][0];
133 : T* pa_p=&a_p[0][0];
134 : for (Int i=0; i<n*n; i++) *(pa_p++)+=*po++;
135 : return *this;
136 : }
137 : }
138 : }
139 : return *this;
140 : }
141 :
142 : template <class T, Int n>
143 : SquareMatrix<T,n>& SquareMatrix<T,n>::operator*=(const SquareMatrix<T,n>& other) {
144 : switch (type_p) {
145 : case ScalarId:
146 : switch (other.type_p) {
147 : case ScalarId: {
148 : a_p[0][0]*=other.a_p[0][0];
149 : return *this;
150 : }
151 : case Diagonal: {
152 : T tmp=a_p[0][0];
153 : for (Int i=0; i<n; i++) {
154 : a_p[i][i]=tmp; a_p[i][i]*=other.a_p[i][i];
155 : }
156 : type_p=Diagonal;
157 : return *this;
158 : }
159 : case General: {
160 : T tmp=a_p[0][0];
161 : for (Int i=0; i<n; i++)
162 : for (Int j=0; j<n; j++) {
163 : a_p[i][j]=tmp; a_p[i][j]*=other.a_p[i][j];
164 : }
165 : type_p=General;
166 : return *this;
167 : }
168 : }
169 : CASACORE_FALLTHROUGH;
170 : case Diagonal:
171 : switch (other.type_p) {
172 : case ScalarId: {
173 : for (Int i=0; i<n; i++) a_p[i][i]*=other.a_p[0][0];
174 : return *this;
175 : }
176 : case Diagonal: {
177 : for (Int i=0; i<n; i++) a_p[i][i]*=other.a_p[i][i];
178 : return *this;
179 : }
180 : case General: {
181 : T a[n];
182 : Int i;
183 : for (i=0; i<n; i++) a[i]=a_p[i][i];
184 : for (i=0; i<n; i++) {
185 : for (Int j=0; j<n; j++) {
186 : a_p[i][j]=a[i]; a_p[i][j]*=other.a_p[i][j];
187 : }
188 : }
189 : type_p=General;
190 : return *this;
191 : }
192 : }
193 : CASACORE_FALLTHROUGH;
194 : case General:
195 : switch (other.type_p) {
196 : case ScalarId: {
197 : for (Int i=0; i<n; i++)
198 : for (Int j=0; j<n; j++) a_p[i][j]*=other.a_p[0][0];
199 : return *this;
200 : }
201 : case Diagonal: {
202 : for (Int i=0; i<n; i++)
203 : for (Int j=0; j<n; j++) a_p[i][j]*=other.a_p[j][j];
204 : return *this;
205 : }
206 : // case General:
207 : default: {
208 : T a[n], tmp;
209 : for (Int i=0; i<n; i++) {
210 : Int j;
211 : for (j=0; j<n; j++) a[j]=a_p[i][j];
212 : for (j=0; j<n; j++) {
213 : a_p[i][j]=a[0]; a_p[i][j]*=other.a_p[0][j];
214 : for (Int k=1; k<n; k++) {
215 : //#a_p[i][j]+=a[k]*other.a_p[k][j]; inlining fails
216 : tmp=a[k]; tmp*=other.a_p[k][j]; a_p[i][j]+=tmp;
217 : }
218 : }
219 : }
220 : return *this;
221 : }
222 : }
223 : }
224 : return *this;
225 : }
226 : template <class T, Int n>
227 : SquareMatrix<T,n>& SquareMatrix<T,n>::operator*=(Float f)
228 : {
229 : switch (type_p) {
230 : case ScalarId: a_p[0][0]*=f; break;
231 : case Diagonal: {
232 : for (Int i=0; i<n; i++) a_p[i][i]*=f;
233 : break;
234 : }
235 : case General: {
236 : T* pa_p=&a_p[0][0];
237 : for (Int i=0; i<n*n; i++) *pa_p++*=f;
238 : }
239 : }
240 : return *this;
241 : }
242 :
243 : /* fails to compile - use explicitly instantiated global function instead
244 : template <class T, Int n>
245 : SquareMatrix<T,n*n>& SquareMatrix<T,n>::directProduct(SquareMatrix<T,n*n>& dp,
246 : const SquareMatrix<T,n>& other) const
247 : {
248 : switch (type_p) {
249 : case ScalarId:
250 : switch (other.type_p) {
251 : case ScalarId: {
252 : dp.a_p[0][0]=a_p[0][0]*other.a_p[0][0];
253 : dp.type_p=ScalarId;
254 : return dp;
255 : }
256 : case Diagonal: {
257 : T tmp=a_p[0][0];
258 : for (Int i=0; i<n*n; i++) dp.a_p[i][i]=tmp*other.a_p[i%n][i%n];
259 : dp.type_p=Diagonal;
260 : return dp;
261 : }
262 : case General: {
263 : T tmp=a_p[0][0];
264 : for (Int i=0; i<n*n; i++)
265 : for (Int j=0; j<n*n; j++) {
266 : if (i/n == j/n) dp.a_p[i][j]=tmp*other.a_p[i%n][j%n];
267 : else dp.a_p[i][j]=T();
268 : }
269 : dp.type_p=General
270 : return dp;
271 : }
272 : }
273 : case Diagonal:
274 : switch (other.type_p) {
275 : case ScalarId: {
276 : T tmp=other.a_p[0][0];
277 : for (Int i=0; i<n*n; i++) dp.a_p[i][i]=a_p[i/n][i/n]*tmp;
278 : dp.type_p=Diagonal;
279 : return dp;
280 : }
281 : case Diagonal: {
282 : for (Int i=0; i<n*n; i++)
283 : dp.a_p[i][i]=a_p[i/n][i/n]*other.a_p[i%n][i%n];
284 : dp.type_p=Diagonal;
285 : return dp;
286 : }
287 : case General: {
288 : for (Int i=0; i<n*n; i++) {
289 : for (Int j=0; j<n*n; j++) {
290 : if (i/n == j/n)
291 : dp.a_p[i][j]=a_p[i/n][i/n]*other.a_p[i%n][j%n];
292 : else dp.a_p[i][j]=T();
293 : }
294 : }
295 : dp.type_p=General;
296 : return dp;
297 : }
298 : }
299 : case General:
300 : switch (other.type_p) {
301 : case ScalarId: {
302 : T tmp=other.a_p[0][0];
303 : for (Int i=0; i<n*n; i++)
304 : for (Int j=0; j<n*n; j++) {
305 : if (i%n == j%n) dp.a_p[i][j]=a_p[i/n][j/n]*tmp;
306 : else dp.a_p[i][j]=T();
307 : }
308 : dp.type_p=General;
309 : return dp;
310 : }
311 : case Diagonal: {
312 : for (Int i=0; i<n*n; i++)
313 : for (Int j=0; j<n*n; j++) {
314 : if (i%n == j%n)
315 : dp.a_p[i][j]=a_p[i/n][j/n]*other.a_p[i%n][j%n];
316 : else dp.a_p[i][j]=T();
317 : }
318 : dp.type_p=General;
319 : return dp;
320 : }
321 : case General: {
322 : for (Int i=0; i<n*n; i++)
323 : for (Int j=0; j<n*n; j++)
324 : dp.a_p[i][j]=a_p[i/n][j/n]*other.a_p[i%n][j%n];
325 : dp.type_p=General;
326 : return dp;
327 : }
328 : }
329 : }
330 : }
331 : */
332 : //# above instantiated for T=Complex, n=2 in SquareMatrix2.cc
333 :
334 : template <class T, Int n>
335 : SquareMatrix<T,n>& SquareMatrix<T,n>::conj() {
336 : switch (type_p) {
337 : case ScalarId: {
338 : a_p[0][0]=std::conj(a_p[0][0]);
339 : return *this;
340 : }
341 : case Diagonal: {
342 : for (Int i=0; i<n; i++) a_p[i][i]=std::conj(a_p[i][i]);
343 : return *this;
344 : }
345 : // case General:
346 : default: {
347 : for (Int i=0; i<n; i++)
348 : for (Int j=0; j<n; j++) a_p[i][j]=std::conj(a_p[i][j]);
349 : return *this;
350 : }
351 : }
352 : }
353 :
354 : template <class T, Int n>
355 : SquareMatrix<T,n>& SquareMatrix<T,n>::adjoint() {
356 : switch (type_p) {
357 : case ScalarId: {
358 : a_p[0][0]=std::conj(a_p[0][0]);
359 : return *this;
360 : }
361 : case Diagonal: {
362 : for (Int i=0; i<n; i++) a_p[i][i]=std::conj(a_p[i][i]);
363 : return *this;
364 : }
365 : case General: {
366 : for (Int i=0; i<n; i++) {
367 : a_p[i][i]=std::conj(a_p[i][i]);
368 : for (Int j=i+1; j<n; j++) {
369 : T tmp=std::conj(a_p[i][j]);
370 : a_p[i][j]=std::conj(a_p[j][i]);
371 : a_p[j][i]=tmp;
372 : }
373 : }
374 : return *this;
375 : }
376 : }
377 : return *this;
378 : }
379 :
380 : template <class T, Int n>
381 : SquareMatrix<T,n>& SquareMatrix<T,n>::conj(SquareMatrix<T,n>& result) {
382 : result = *this;
383 : result.conj();
384 : return result;
385 : }
386 :
387 : template <class T, Int n>
388 : SquareMatrix<T,n>& SquareMatrix<T,n>::adjoint(SquareMatrix<T,n>& result) {
389 : result = *this;
390 : result.adjoint();
391 : return result;
392 : }
393 :
394 : template <class T, Int n>
395 : SquareMatrix<T,n>& SquareMatrix<T,n>::inverse(SquareMatrix<T,n>& result) const {
396 : switch (type_p) {
397 : case ScalarId: {
398 : result.a_p[0][0]=T(1)/a_p[0][0];
399 : result.type_p=ScalarId;
400 : return result;
401 : }
402 : case Diagonal: {
403 : for (Int i=0; i<n; i++) result.a_p[i][i]=T(1)/a_p[i][i];
404 : result.type_p=Diagonal;
405 : return result;
406 : }
407 : // case General:
408 : default: {
409 : switch (n) {
410 : case 2: {
411 : T det=a_p[0][0]*a_p[1][1]-a_p[1][0]*a_p[0][1];
412 : result.a_p[0][0]=a_p[1][1]/det;
413 : result.a_p[0][1]=-a_p[0][1]/det;
414 : result.a_p[1][0]=-a_p[1][0]/det;
415 : result.a_p[1][1]=a_p[0][0]/det;
416 : result.type_p=General;
417 : return result;
418 : }
419 : default: {
420 : // return result=invert(matrix());
421 : Matrix<T> mat=invert(matrix());
422 : if (mat.nelements()==0) {
423 : cerr<< "invert of singular matrix attempted:"<<
424 : matrix()
425 : << endl;
426 : result=T(1);
427 : }
428 : else result=mat;
429 : return result;
430 : }
431 : }
432 : }
433 : }
434 : }
435 :
436 : template <class T, Int n>
437 0 : Matrix<T>& SquareMatrix<T,n>::matrix(Matrix<T>& result) const {
438 0 : result.resize(n,n);
439 0 : switch (type_p) {
440 0 : case ScalarId: {
441 0 : result=T();
442 0 : for (Int i=0; i<n; i++) result(i,i)=a_p[0][0];
443 0 : return result;
444 : }
445 0 : case Diagonal: {
446 0 : result=T();
447 0 : for (Int i=0; i<n; i++) result(i,i)=a_p[i][i];
448 0 : return result;
449 : }
450 : // case General:
451 0 : default: {
452 0 : for (Int i=0; i<n; i++)
453 0 : for (Int j=0; j<n; j++) result(i,j)=a_p[i][j];
454 0 : return result;
455 : }
456 : }
457 : }
458 :
459 : template <class T, Int n>
460 0 : T& SquareMatrix<T,n>::throwInvAccess() {
461 : throw(AipsError("SquareMatrix - attempt to change element that is "
462 0 : "not available for this type of matrix"));
463 : // following just to make signature ok.
464 : return a_p[0][0];
465 : }
466 :
467 : } //# NAMESPACE CASACORE - END
468 :
469 :
470 : #endif
|
