source: branches/blt4/packages/vizservers/nanovis/newmat11/newmat1.cpp @ 3892

Last change on this file since 3892 was 3892, checked in by gah, 11 years ago
File size: 6.0 KB
Line 
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.