newmat1.cpp
Go to the documentation of this file.
1 
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 
38 {
39  REPORT
40  int a = attribute & mt.attribute & ~(Symmetric | Skew);
41  a |= (a & Diagonal) * 63; // recognise diagonal
42  return MatrixType(a);
43 }
44 
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 
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 
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 
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
109 {
110  REPORT
111  return
112  ((a.attribute | b.attribute | c.attribute)
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 
GeneralMatrix * New() const
new matrix of given type
Miscellaneous exception (details in character string).
Definition: newmat.h:1947
MatrixType i() const
type of inverse
Definition: newmat1.cpp:82
#define REPORT
Definition: newmat1.cpp:22
MatrixType KP(const MatrixType &) const
Definition: newmat1.cpp:67
virtual MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:671
int attribute
Definition: newmat.h:139
MatrixType t() const
type of transpose
Definition: newmat1.cpp:90
Upper triangular band matrix.
Definition: newmat.h:1167
A matrix is not square exception.
Definition: newmat.h:1981
Upper triangular matrix.
Definition: newmat.h:799
const char * value() const
type as char string
Definition: newmat1.cpp:116
Square matrix.
Definition: newmat.h:679
MatrixType SP(const MatrixType &) const
Definition: newmat1.cpp:45
Band matrix.
Definition: newmat.h:1096
Base of the matrix classes.
Definition: newmat.h:292
#define Throw(E)
Definition: myexcept.h:191
The usual rectangular matrix.
Definition: newmat.h:625
FloatVector FloatVector * a
MatrixType MultRHS() const
type for rhs of multiply
Definition: newmat1.cpp:100
bool Rectangular(MatrixType a, MatrixType b, MatrixType c)
Definition: newmat1.cpp:108
void MatrixErrorNoSpace(const void *)
test for allocation fails
Definition: newmatex.cpp:301
Diagonal matrix.
Definition: newmat.h:896
Lower triangular matrix.
Definition: newmat.h:848
Symmetric band matrix.
Definition: newmat.h:1245
Lower triangular band matrix.
Definition: newmat.h:1206
Row vector.
Definition: newmat.h:953
void Protect()
Definition: newmat.h:509
Column vector.
Definition: newmat.h:1008
The classes for matrices that can contain data are derived from this.
Definition: newmat.h:447
MatrixType operator*(const MatrixType &) const
Definition: newmat1.cpp:37
Identity matrix.
Definition: newmat.h:1350
Symmetric matrix.
Definition: newmat.h:753


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:16