source: ntrip/trunk/BNC/newmat/newmat1.cpp@ 6381

Last change on this file since 6381 was 2013, checked in by mervart, 15 years ago

* empty log message *

File size: 6.0 KB
RevLine 
[2013]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
16namespace 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
37MatrixType 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
45MatrixType 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
67MatrixType 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
82MatrixType 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
90MatrixType 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
100MatrixType 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
108bool 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
116const 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
146GeneralMatrix* 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///@}
Note: See TracBrowser for help on using the repository browser.