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