1 | /// \ingroup newmat |
---|
2 | ///@{ |
---|
3 | |
---|
4 | /// \file newmat1.cpp |
---|
5 | /// MatrixType functions. |
---|
6 | /// Find the type of a matrix resulting from a multiply, add etc |
---|
7 | /// Make a new matrix corresponding to a MatrixType |
---|
8 | |
---|
9 | // Copyright (C) 1991,2,3,4: R B Davies |
---|
10 | |
---|
11 | //#define WANT_STREAM |
---|
12 | |
---|
13 | #include "newmat.h" |
---|
14 | |
---|
15 | #ifdef use_namespace |
---|
16 | namespace NEWMAT { |
---|
17 | #endif |
---|
18 | |
---|
19 | #ifdef DO_REPORT |
---|
20 | #define REPORT { static ExeCounter ExeCount(__LINE__,1); ++ExeCount; } |
---|
21 | #else |
---|
22 | #define REPORT {} |
---|
23 | #endif |
---|
24 | |
---|
25 | |
---|
26 | /************************* MatrixType functions *****************************/ |
---|
27 | |
---|
28 | |
---|
29 | // Skew needs more work <<<<<<<<< |
---|
30 | |
---|
31 | // all operations to return MatrixTypes which correspond to valid types |
---|
32 | // of matrices. |
---|
33 | // Eg: if it has the Diagonal attribute, then it must also have |
---|
34 | // Upper, Lower, Band, Square and Symmetric. |
---|
35 | |
---|
36 | |
---|
37 | MatrixType MatrixType::operator*(const MatrixType& mt) const |
---|
38 | { |
---|
39 | REPORT |
---|
40 | int a = attribute & mt.attribute & ~(Symmetric | Skew); |
---|
41 | a |= (a & Diagonal) * 63; // recognise diagonal |
---|
42 | return MatrixType(a); |
---|
43 | } |
---|
44 | |
---|
45 | MatrixType MatrixType::SP(const MatrixType& mt) const |
---|
46 | // elementwise product |
---|
47 | // Lower, Upper, Diag, Band if only one is |
---|
48 | // Symmetric, Ones, Valid (and Real) if both are |
---|
49 | // Need to include Lower & Upper => Diagonal |
---|
50 | // Will need to include both Skew => Symmetric |
---|
51 | { |
---|
52 | REPORT |
---|
53 | int a = ((attribute | mt.attribute) & ~(Symmetric + Skew + Valid + Ones)) |
---|
54 | | (attribute & mt.attribute); |
---|
55 | if ((a & Lower) != 0 && (a & Upper) != 0) a |= Diagonal; |
---|
56 | if ((attribute & Skew) != 0) |
---|
57 | { |
---|
58 | if ((mt.attribute & Symmetric) != 0) a |= Skew; |
---|
59 | if ((mt.attribute & Skew) != 0) { a &= ~Skew; a |= Symmetric; } |
---|
60 | } |
---|
61 | else if ((mt.attribute & Skew) != 0 && (attribute & Symmetric) != 0) |
---|
62 | a |= Skew; |
---|
63 | a |= (a & Diagonal) * 63; // recognise diagonal |
---|
64 | return MatrixType(a); |
---|
65 | } |
---|
66 | |
---|
67 | MatrixType MatrixType::KP(const MatrixType& mt) const |
---|
68 | // Kronecker product |
---|
69 | // Lower, Upper, Diag, Symmetric, Band, Valid if both are |
---|
70 | // Band if LHS is band & other is square |
---|
71 | // Ones is complicated so leave this out |
---|
72 | { |
---|
73 | REPORT |
---|
74 | int a = (attribute & mt.attribute) & ~Ones; |
---|
75 | if ((attribute & Band) != 0 && (mt.attribute & Square) != 0) |
---|
76 | a |= Band; |
---|
77 | //int a = ((attribute & mt.attribute) | (attribute & Band)) & ~Ones; |
---|
78 | |
---|
79 | return MatrixType(a); |
---|
80 | } |
---|
81 | |
---|
82 | MatrixType MatrixType::i() const // type of inverse |
---|
83 | { |
---|
84 | REPORT |
---|
85 | int a = attribute & ~(Band+LUDeco); |
---|
86 | a |= (a & Diagonal) * 63; // recognise diagonal |
---|
87 | return MatrixType(a); |
---|
88 | } |
---|
89 | |
---|
90 | MatrixType MatrixType::t() const |
---|
91 | // swap lower and upper attributes |
---|
92 | // assume Upper is in bit above Lower |
---|
93 | { |
---|
94 | REPORT |
---|
95 | int a = attribute; |
---|
96 | a ^= (((a >> 1) ^ a) & Lower) * 3; |
---|
97 | return MatrixType(a); |
---|
98 | } |
---|
99 | |
---|
100 | MatrixType MatrixType::MultRHS() const |
---|
101 | { |
---|
102 | REPORT |
---|
103 | // remove symmetric attribute unless diagonal |
---|
104 | return (attribute >= Dg) ? attribute : (attribute & ~Symmetric); |
---|
105 | } |
---|
106 | |
---|
107 | // this is used for deciding type of multiplication |
---|
108 | bool Rectangular(MatrixType a, MatrixType b, MatrixType c) |
---|
109 | { |
---|
110 | REPORT |
---|
111 | return |
---|
112 | ((a.attribute | b.attribute | c.attribute) |
---|
113 | & ~(MatrixType::Valid | MatrixType::Square)) == 0; |
---|
114 | } |
---|
115 | |
---|
116 | const char* MatrixType::value() const |
---|
117 | { |
---|
118 | // make a string with the name of matrix with the given attributes |
---|
119 | switch (attribute) |
---|
120 | { |
---|
121 | case Valid: REPORT return "Rect "; |
---|
122 | case Valid+Square: REPORT return "Squ "; |
---|
123 | case Valid+Symmetric+Square: REPORT return "Sym "; |
---|
124 | case Valid+Skew+Square: REPORT return "Skew "; |
---|
125 | case Valid+Band+Square: REPORT return "Band "; |
---|
126 | case Valid+Symmetric+Band+Square: REPORT return "SmBnd"; |
---|
127 | case Valid+Skew+Band+Square: REPORT return "SkBnd"; |
---|
128 | case Valid+Upper+Square: REPORT return "UT "; |
---|
129 | case Valid+Diagonal+Symmetric+Band+Upper+Lower+Square: |
---|
130 | REPORT return "Diag "; |
---|
131 | case Valid+Diagonal+Symmetric+Band+Upper+Lower+Ones+Square: |
---|
132 | REPORT return "Ident"; |
---|
133 | case Valid+Band+Upper+Square: REPORT return "UpBnd"; |
---|
134 | case Valid+Lower+Square: REPORT return "LT "; |
---|
135 | case Valid+Band+Lower+Square: REPORT return "LwBnd"; |
---|
136 | default: |
---|
137 | REPORT |
---|
138 | if (!(attribute & Valid)) return "UnSp "; |
---|
139 | if (attribute & LUDeco) |
---|
140 | return (attribute & Band) ? "BndLU" : "Crout"; |
---|
141 | return "?????"; |
---|
142 | } |
---|
143 | } |
---|
144 | |
---|
145 | |
---|
146 | GeneralMatrix* MatrixType::New(int nr, int nc, BaseMatrix* bm) const |
---|
147 | { |
---|
148 | // make a new matrix with the given attributes |
---|
149 | |
---|
150 | Tracer tr("New"); GeneralMatrix* gm=0; // initialised to keep gnu happy |
---|
151 | switch (attribute) |
---|
152 | { |
---|
153 | case Valid: |
---|
154 | REPORT |
---|
155 | if (nc==1) { gm = new ColumnVector(nr); break; } |
---|
156 | if (nr==1) { gm = new RowVector(nc); break; } |
---|
157 | gm = new Matrix(nr, nc); break; |
---|
158 | |
---|
159 | case Valid+Square: |
---|
160 | REPORT |
---|
161 | if (nc!=nr) { Throw(NotSquareException()); } |
---|
162 | gm = new SquareMatrix(nr); break; |
---|
163 | |
---|
164 | case Valid+Symmetric+Square: |
---|
165 | REPORT gm = new SymmetricMatrix(nr); break; |
---|
166 | |
---|
167 | case Valid+Band+Square: |
---|
168 | { |
---|
169 | REPORT |
---|
170 | MatrixBandWidth bw = bm->bandwidth(); |
---|
171 | gm = new BandMatrix(nr,bw.lower_val,bw.upper_val); break; |
---|
172 | } |
---|
173 | |
---|
174 | case Valid+Symmetric+Band+Square: |
---|
175 | REPORT gm = new SymmetricBandMatrix(nr,bm->bandwidth().lower_val); break; |
---|
176 | |
---|
177 | case Valid+Upper+Square: |
---|
178 | REPORT gm = new UpperTriangularMatrix(nr); break; |
---|
179 | |
---|
180 | case Valid+Diagonal+Symmetric+Band+Upper+Lower+Square: |
---|
181 | REPORT gm = new DiagonalMatrix(nr); break; |
---|
182 | |
---|
183 | case Valid+Band+Upper+Square: |
---|
184 | REPORT gm = new UpperBandMatrix(nr,bm->bandwidth().upper_val); break; |
---|
185 | |
---|
186 | case Valid+Lower+Square: |
---|
187 | REPORT gm = new LowerTriangularMatrix(nr); break; |
---|
188 | |
---|
189 | case Valid+Band+Lower+Square: |
---|
190 | REPORT gm = new LowerBandMatrix(nr,bm->bandwidth().lower_val); break; |
---|
191 | |
---|
192 | case Valid+Diagonal+Symmetric+Band+Upper+Lower+Ones+Square: |
---|
193 | REPORT gm = new IdentityMatrix(nr); break; |
---|
194 | |
---|
195 | default: |
---|
196 | Throw(ProgramException("Invalid matrix type")); |
---|
197 | } |
---|
198 | |
---|
199 | MatrixErrorNoSpace(gm); gm->Protect(); return gm; |
---|
200 | } |
---|
201 | |
---|
202 | |
---|
203 | #ifdef use_namespace |
---|
204 | } |
---|
205 | #endif |
---|
206 | |
---|
207 | |
---|
208 | |
---|
209 | ///@} |
---|