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
operator *(const MatrixType & mt) const37 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
SP(const MatrixType & mt) const45 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
KP(const MatrixType & mt) const67 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
i() const82 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
t() const90 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
MultRHS() const100 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
Rectangular(MatrixType a,MatrixType b,MatrixType c)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
value() const116 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
New(int nr,int nc,BaseMatrix * bm) const146 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 ///@}
210