Bart Janssens / SVL

Simple Vector Library 1.5 http://www.cs.cmu.edu/~ajw/doc/svl.html

Mat.cpp@1:e25ff4b06ed2, 2016-01-05 (annotated)

Committer:
BartJanssens
Date:
Tue Jan 05 13:37:50 2016 +0000
Revision:
1:e25ff4b06ed2
Parent:
0:785cff1e5a7c 
fixed va_list bug

Who changed what in which revision?

UserRevisionLine numberNew contents of line
BartJanssens 0:785cff1e5a7c 1 /*
BartJanssens 0:785cff1e5a7c 2 File: Mat.cpp
BartJanssens 0:785cff1e5a7c 3
BartJanssens 0:785cff1e5a7c 4 Function: Implements Mat.h
BartJanssens 0:785cff1e5a7c 5
BartJanssens 0:785cff1e5a7c 6 Author(s): Andrew Willmott
BartJanssens 0:785cff1e5a7c 7
BartJanssens 0:785cff1e5a7c 8 Copyright: (c) 1995-2001, Andrew Willmott
BartJanssens 0:785cff1e5a7c 9
BartJanssens 0:785cff1e5a7c 10 */
BartJanssens 0:785cff1e5a7c 11
BartJanssens 0:785cff1e5a7c 12 #include "Mat.h"
BartJanssens 0:785cff1e5a7c 13
BartJanssens 0:785cff1e5a7c 14 //#include <cctype>
BartJanssens 0:785cff1e5a7c 15 //#include <cstring>
BartJanssens 1:e25ff4b06ed2 16 #include <cstdarg>
BartJanssens 0:785cff1e5a7c 17 //#include <iomanip>
BartJanssens 0:785cff1e5a7c 18
BartJanssens 0:785cff1e5a7c 19
BartJanssens 0:785cff1e5a7c 20
BartJanssens 0:785cff1e5a7c 21 // --- Mat Constructors & Destructors -----------------------------------------
BartJanssens 0:785cff1e5a7c 22
BartJanssens 0:785cff1e5a7c 23
BartJanssens 0:785cff1e5a7c 24 Mat::Mat(Int rows, Int cols, ZeroOrOne k) : rows(rows), cols(cols)
BartJanssens 0:785cff1e5a7c 25 {
BartJanssens 0:785cff1e5a7c 26 Assert(rows > 0 && cols > 0, "(Mat) illegal matrix size");
BartJanssens 0:785cff1e5a7c 27
BartJanssens 0:785cff1e5a7c 28 data = new Real[rows * cols];
BartJanssens 0:785cff1e5a7c 29
BartJanssens 0:785cff1e5a7c 30 MakeDiag(k);
BartJanssens 0:785cff1e5a7c 31 }
BartJanssens 0:785cff1e5a7c 32
BartJanssens 0:785cff1e5a7c 33 Mat::Mat(Int rows, Int cols, Block k) : rows(rows), cols(cols)
BartJanssens 0:785cff1e5a7c 34 {
BartJanssens 0:785cff1e5a7c 35 Assert(rows > 0 && cols > 0, "(Mat) illegal matrix size");
BartJanssens 0:785cff1e5a7c 36
BartJanssens 0:785cff1e5a7c 37 data = new Real[rows * cols];
BartJanssens 0:785cff1e5a7c 38
BartJanssens 0:785cff1e5a7c 39 MakeBlock(k);
BartJanssens 0:785cff1e5a7c 40 }
BartJanssens 0:785cff1e5a7c 41
BartJanssens 0:785cff1e5a7c 42 Mat::Mat(Int rows, Int cols, double elt0, ...) : rows(rows), cols(cols)
BartJanssens 0:785cff1e5a7c 43 // The double is hardwired here because it is the only type that will work
BartJanssens 0:785cff1e5a7c 44 // with var args and C++ real numbers.
BartJanssens 0:785cff1e5a7c 45 {
BartJanssens 0:785cff1e5a7c 46 Assert(rows > 0 && cols > 0, "(Mat) illegal matrix size");
BartJanssens 0:785cff1e5a7c 47
BartJanssens 0:785cff1e5a7c 48 va_list ap;
BartJanssens 0:785cff1e5a7c 49 Int i, j;
BartJanssens 0:785cff1e5a7c 50
BartJanssens 0:785cff1e5a7c 51 data = new Real[rows * cols];
BartJanssens 0:785cff1e5a7c 52 va_start(ap, elt0);
BartJanssens 0:785cff1e5a7c 53
BartJanssens 0:785cff1e5a7c 54 SetReal(data[0], elt0);
BartJanssens 0:785cff1e5a7c 55
BartJanssens 0:785cff1e5a7c 56 for (i = 1; i < cols; i++)
BartJanssens 0:785cff1e5a7c 57 SetReal(Elt(0, i), va_arg(ap, double));
BartJanssens 0:785cff1e5a7c 58
BartJanssens 0:785cff1e5a7c 59 for (i = 1; i < rows; i++)
BartJanssens 0:785cff1e5a7c 60 for (j = 0; j < cols; j++)
BartJanssens 0:785cff1e5a7c 61 SetReal(Elt(i, j), va_arg(ap, double));
BartJanssens 0:785cff1e5a7c 62
BartJanssens 0:785cff1e5a7c 63 va_end(ap);
BartJanssens 0:785cff1e5a7c 64 }
BartJanssens 0:785cff1e5a7c 65
BartJanssens 0:785cff1e5a7c 66 Mat::Mat(const Mat &m) : cols(m.cols)
BartJanssens 0:785cff1e5a7c 67 {
BartJanssens 0:785cff1e5a7c 68 Assert(m.data != 0, "(Mat) Can't construct from null matrix");
BartJanssens 0:785cff1e5a7c 69 rows = m.Rows();
BartJanssens 0:785cff1e5a7c 70
BartJanssens 0:785cff1e5a7c 71 UInt elts = rows * cols;
BartJanssens 0:785cff1e5a7c 72
BartJanssens 0:785cff1e5a7c 73 data = new Real[elts];
BartJanssens 0:785cff1e5a7c 74 #ifdef VL_USE_MEMCPY
BartJanssens 0:785cff1e5a7c 75 memcpy(data, m.data, elts * sizeof(Real));
BartJanssens 0:785cff1e5a7c 76 #else
BartJanssens 0:785cff1e5a7c 77 for (UInt i = 0; i < elts; i++)
BartJanssens 0:785cff1e5a7c 78 data[i] = m.data[i];
BartJanssens 0:785cff1e5a7c 79 #endif
BartJanssens 0:785cff1e5a7c 80 }
BartJanssens 0:785cff1e5a7c 81
BartJanssens 0:785cff1e5a7c 82 Mat::Mat(const Mat2 &m) : data(m.Ref()), rows(2 | VL_REF_FLAG), cols(2)
BartJanssens 0:785cff1e5a7c 83 {
BartJanssens 0:785cff1e5a7c 84 }
BartJanssens 0:785cff1e5a7c 85
BartJanssens 0:785cff1e5a7c 86 Mat::Mat(const Mat3 &m) : data(m.Ref()), rows(3 | VL_REF_FLAG), cols(3)
BartJanssens 0:785cff1e5a7c 87 {
BartJanssens 0:785cff1e5a7c 88 }
BartJanssens 0:785cff1e5a7c 89
BartJanssens 0:785cff1e5a7c 90 Mat::Mat(const Mat4 &m) : data(m.Ref()), rows(4 | VL_REF_FLAG), cols(4)
BartJanssens 0:785cff1e5a7c 91 {
BartJanssens 0:785cff1e5a7c 92 }
BartJanssens 0:785cff1e5a7c 93
BartJanssens 0:785cff1e5a7c 94
BartJanssens 0:785cff1e5a7c 95 // --- Mat Assignment Operators -----------------------------------------------
BartJanssens 0:785cff1e5a7c 96
BartJanssens 0:785cff1e5a7c 97
BartJanssens 0:785cff1e5a7c 98 Mat &Mat::operator = (const Mat &m)
BartJanssens 0:785cff1e5a7c 99 {
BartJanssens 0:785cff1e5a7c 100 if (!IsRef())
BartJanssens 0:785cff1e5a7c 101 SetSize(m.Rows(), m.Cols());
BartJanssens 0:785cff1e5a7c 102 else
BartJanssens 0:785cff1e5a7c 103 Assert(Rows() == m.Rows(), "(Mat::=) Matrix rows don't match");
BartJanssens 0:785cff1e5a7c 104 for (Int i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 105 SELF[i] = m[i];
BartJanssens 0:785cff1e5a7c 106
BartJanssens 0:785cff1e5a7c 107 return(SELF);
BartJanssens 0:785cff1e5a7c 108 }
BartJanssens 0:785cff1e5a7c 109
BartJanssens 0:785cff1e5a7c 110 Mat &Mat::operator = (const Mat2 &m)
BartJanssens 0:785cff1e5a7c 111 {
BartJanssens 0:785cff1e5a7c 112 if (!IsRef())
BartJanssens 0:785cff1e5a7c 113 SetSize(m.Rows(), m.Cols());
BartJanssens 0:785cff1e5a7c 114 else
BartJanssens 0:785cff1e5a7c 115 Assert(Rows() == m.Rows(), "(Mat::=) Matrix rows don't match");
BartJanssens 0:785cff1e5a7c 116 for (Int i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 117 SELF[i] = m[i];
BartJanssens 0:785cff1e5a7c 118
BartJanssens 0:785cff1e5a7c 119 return(SELF);
BartJanssens 0:785cff1e5a7c 120 }
BartJanssens 0:785cff1e5a7c 121
BartJanssens 0:785cff1e5a7c 122 Mat &Mat::operator = (const Mat3 &m)
BartJanssens 0:785cff1e5a7c 123 {
BartJanssens 0:785cff1e5a7c 124 if (!IsRef())
BartJanssens 0:785cff1e5a7c 125 SetSize(m.Rows(), m.Cols());
BartJanssens 0:785cff1e5a7c 126 else
BartJanssens 0:785cff1e5a7c 127 Assert(Rows() == m.Rows(), "(Mat::=) Matrix rows don't match");
BartJanssens 0:785cff1e5a7c 128 for (Int i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 129 SELF[i] = m[i];
BartJanssens 0:785cff1e5a7c 130
BartJanssens 0:785cff1e5a7c 131 return(SELF);
BartJanssens 0:785cff1e5a7c 132 }
BartJanssens 0:785cff1e5a7c 133
BartJanssens 0:785cff1e5a7c 134 Mat &Mat::operator = (const Mat4 &m)
BartJanssens 0:785cff1e5a7c 135 {
BartJanssens 0:785cff1e5a7c 136 if (!IsRef())
BartJanssens 0:785cff1e5a7c 137 SetSize(m.Rows(), m.Cols());
BartJanssens 0:785cff1e5a7c 138 else
BartJanssens 0:785cff1e5a7c 139 Assert(Rows() == m.Rows(), "(Mat::=) Matrix rows don't match");
BartJanssens 0:785cff1e5a7c 140
BartJanssens 0:785cff1e5a7c 141 for (Int i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 142 SELF[i] = m[i];
BartJanssens 0:785cff1e5a7c 143
BartJanssens 0:785cff1e5a7c 144 return(SELF);
BartJanssens 0:785cff1e5a7c 145 }
BartJanssens 0:785cff1e5a7c 146
BartJanssens 0:785cff1e5a7c 147 Void Mat::SetSize(Int nrows, Int ncols)
BartJanssens 0:785cff1e5a7c 148 {
BartJanssens 0:785cff1e5a7c 149 UInt elts = nrows * ncols;
BartJanssens 0:785cff1e5a7c 150 Assert(nrows > 0 && ncols > 0, "(Mat::SetSize) Illegal matrix size.");
BartJanssens 0:785cff1e5a7c 151 UInt oldElts = Rows() * Cols();
BartJanssens 0:785cff1e5a7c 152
BartJanssens 0:785cff1e5a7c 153 if (IsRef())
BartJanssens 0:785cff1e5a7c 154 {
BartJanssens 0:785cff1e5a7c 155 // Abort! We don't allow this operation on references.
BartJanssens 0:785cff1e5a7c 156 _Error("(Mat::SetSize) Trying to resize a matrix reference");
BartJanssens 0:785cff1e5a7c 157 }
BartJanssens 0:785cff1e5a7c 158
BartJanssens 0:785cff1e5a7c 159 rows = nrows;
BartJanssens 0:785cff1e5a7c 160 cols = ncols;
BartJanssens 0:785cff1e5a7c 161
BartJanssens 0:785cff1e5a7c 162 // Don't reallocate if we already have enough storage
BartJanssens 0:785cff1e5a7c 163 if (elts <= oldElts)
BartJanssens 0:785cff1e5a7c 164 return;
BartJanssens 0:785cff1e5a7c 165
BartJanssens 0:785cff1e5a7c 166 // Otherwise, delete old storage and reallocate
BartJanssens 0:785cff1e5a7c 167 delete[] data;
BartJanssens 0:785cff1e5a7c 168 data = 0;
BartJanssens 0:785cff1e5a7c 169 data = new Real[elts]; // may throw an exception
BartJanssens 0:785cff1e5a7c 170 }
BartJanssens 0:785cff1e5a7c 171
BartJanssens 0:785cff1e5a7c 172 Void Mat::SetSize(const Mat &m)
BartJanssens 0:785cff1e5a7c 173 {
BartJanssens 0:785cff1e5a7c 174 SetSize(m.Rows(), m.Cols());
BartJanssens 0:785cff1e5a7c 175 }
BartJanssens 0:785cff1e5a7c 176
BartJanssens 0:785cff1e5a7c 177 Void Mat::MakeZero()
BartJanssens 0:785cff1e5a7c 178 {
BartJanssens 0:785cff1e5a7c 179 #ifdef VL_USE_MEMCPY
BartJanssens 0:785cff1e5a7c 180 memset(data, 0, sizeof(Real) * Rows() * Cols());
BartJanssens 0:785cff1e5a7c 181 #else
BartJanssens 0:785cff1e5a7c 182 Int i;
BartJanssens 0:785cff1e5a7c 183
BartJanssens 0:785cff1e5a7c 184 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 185 SELF[i] = vl_zero;
BartJanssens 0:785cff1e5a7c 186 #endif
BartJanssens 0:785cff1e5a7c 187 }
BartJanssens 0:785cff1e5a7c 188
BartJanssens 0:785cff1e5a7c 189 Void Mat::MakeDiag(Real k)
BartJanssens 0:785cff1e5a7c 190 {
BartJanssens 0:785cff1e5a7c 191 Int i, j;
BartJanssens 0:785cff1e5a7c 192
BartJanssens 0:785cff1e5a7c 193 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 194 for (j = 0; j < Cols(); j++)
BartJanssens 0:785cff1e5a7c 195 if (i == j)
BartJanssens 0:785cff1e5a7c 196 Elt(i,j) = k;
BartJanssens 0:785cff1e5a7c 197 else
BartJanssens 0:785cff1e5a7c 198 Elt(i,j) = vl_zero;
BartJanssens 0:785cff1e5a7c 199 }
BartJanssens 0:785cff1e5a7c 200
BartJanssens 0:785cff1e5a7c 201 Void Mat::MakeDiag()
BartJanssens 0:785cff1e5a7c 202 {
BartJanssens 0:785cff1e5a7c 203 Int i, j;
BartJanssens 0:785cff1e5a7c 204
BartJanssens 0:785cff1e5a7c 205 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 206 for (j = 0; j < Cols(); j++)
BartJanssens 0:785cff1e5a7c 207 Elt(i,j) = (i == j) ? vl_one : vl_zero;
BartJanssens 0:785cff1e5a7c 208 }
BartJanssens 0:785cff1e5a7c 209
BartJanssens 0:785cff1e5a7c 210 Void Mat::MakeBlock(Real k)
BartJanssens 0:785cff1e5a7c 211 {
BartJanssens 0:785cff1e5a7c 212 Int i;
BartJanssens 0:785cff1e5a7c 213
BartJanssens 0:785cff1e5a7c 214 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 215 SELF[i].MakeBlock(k);
BartJanssens 0:785cff1e5a7c 216 }
BartJanssens 0:785cff1e5a7c 217
BartJanssens 0:785cff1e5a7c 218 Void Mat::MakeBlock()
BartJanssens 0:785cff1e5a7c 219 {
BartJanssens 0:785cff1e5a7c 220 Int i, j;
BartJanssens 0:785cff1e5a7c 221
BartJanssens 0:785cff1e5a7c 222 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 223 for (j = 0; j < Cols(); j++)
BartJanssens 0:785cff1e5a7c 224 Elt(i,j) = vl_one;
BartJanssens 0:785cff1e5a7c 225 }
BartJanssens 0:785cff1e5a7c 226
BartJanssens 0:785cff1e5a7c 227
BartJanssens 0:785cff1e5a7c 228 // --- Mat Assignment Operators -----------------------------------------------
BartJanssens 0:785cff1e5a7c 229
BartJanssens 0:785cff1e5a7c 230
BartJanssens 0:785cff1e5a7c 231 Mat &Mat::operator += (const Mat &m)
BartJanssens 0:785cff1e5a7c 232 {
BartJanssens 0:785cff1e5a7c 233 Assert(Rows() == m.Rows(), "(Mat::+=) matrix rows don't match");
BartJanssens 0:785cff1e5a7c 234
BartJanssens 0:785cff1e5a7c 235 Int i;
BartJanssens 0:785cff1e5a7c 236
BartJanssens 0:785cff1e5a7c 237 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 238 SELF[i] += m[i];
BartJanssens 0:785cff1e5a7c 239
BartJanssens 0:785cff1e5a7c 240 return(SELF);
BartJanssens 0:785cff1e5a7c 241 }
BartJanssens 0:785cff1e5a7c 242
BartJanssens 0:785cff1e5a7c 243 Mat &Mat::operator -= (const Mat &m)
BartJanssens 0:785cff1e5a7c 244 {
BartJanssens 0:785cff1e5a7c 245 Assert(Rows() == m.Rows(), "(Mat::-=) matrix rows don't match");
BartJanssens 0:785cff1e5a7c 246
BartJanssens 0:785cff1e5a7c 247 Int i;
BartJanssens 0:785cff1e5a7c 248
BartJanssens 0:785cff1e5a7c 249 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 250 SELF[i] -= m[i];
BartJanssens 0:785cff1e5a7c 251
BartJanssens 0:785cff1e5a7c 252 return(SELF);
BartJanssens 0:785cff1e5a7c 253 }
BartJanssens 0:785cff1e5a7c 254
BartJanssens 0:785cff1e5a7c 255 Mat &Mat::operator *= (const Mat &m)
BartJanssens 0:785cff1e5a7c 256 {
BartJanssens 0:785cff1e5a7c 257 Assert(Cols() == m.Cols(), "(Mat::*=) matrix columns don't match");
BartJanssens 0:785cff1e5a7c 258
BartJanssens 0:785cff1e5a7c 259 Int i;
BartJanssens 0:785cff1e5a7c 260
BartJanssens 0:785cff1e5a7c 261 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 262 SELF[i] = SELF[i] * m;
BartJanssens 0:785cff1e5a7c 263
BartJanssens 0:785cff1e5a7c 264 return(SELF);
BartJanssens 0:785cff1e5a7c 265 }
BartJanssens 0:785cff1e5a7c 266
BartJanssens 0:785cff1e5a7c 267 Mat &Mat::operator *= (Real s)
BartJanssens 0:785cff1e5a7c 268 {
BartJanssens 0:785cff1e5a7c 269 Int i;
BartJanssens 0:785cff1e5a7c 270
BartJanssens 0:785cff1e5a7c 271 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 272 SELF[i] *= s;
BartJanssens 0:785cff1e5a7c 273
BartJanssens 0:785cff1e5a7c 274 return(SELF);
BartJanssens 0:785cff1e5a7c 275 }
BartJanssens 0:785cff1e5a7c 276
BartJanssens 0:785cff1e5a7c 277 Mat &Mat::operator /= (Real s)
BartJanssens 0:785cff1e5a7c 278 {
BartJanssens 0:785cff1e5a7c 279 Int i;
BartJanssens 0:785cff1e5a7c 280
BartJanssens 0:785cff1e5a7c 281 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 282 SELF[i] /= s;
BartJanssens 0:785cff1e5a7c 283
BartJanssens 0:785cff1e5a7c 284 return(SELF);
BartJanssens 0:785cff1e5a7c 285 }
BartJanssens 0:785cff1e5a7c 286
BartJanssens 0:785cff1e5a7c 287
BartJanssens 0:785cff1e5a7c 288 // --- Mat Comparison Operators -----------------------------------------------
BartJanssens 0:785cff1e5a7c 289
BartJanssens 0:785cff1e5a7c 290
BartJanssens 0:785cff1e5a7c 291 Bool operator == (const Mat &m, const Mat &n)
BartJanssens 0:785cff1e5a7c 292 {
BartJanssens 0:785cff1e5a7c 293 Assert(n.Rows() == m.Rows(), "(Mat::==) matrix rows don't match");
BartJanssens 0:785cff1e5a7c 294
BartJanssens 0:785cff1e5a7c 295 Int i;
BartJanssens 0:785cff1e5a7c 296
BartJanssens 0:785cff1e5a7c 297 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 298 if (m[i] != n[i])
BartJanssens 0:785cff1e5a7c 299 return(0);
BartJanssens 0:785cff1e5a7c 300
BartJanssens 0:785cff1e5a7c 301 return(1);
BartJanssens 0:785cff1e5a7c 302 }
BartJanssens 0:785cff1e5a7c 303
BartJanssens 0:785cff1e5a7c 304 Bool operator != (const Mat &m, const Mat &n)
BartJanssens 0:785cff1e5a7c 305 {
BartJanssens 0:785cff1e5a7c 306 Assert(n.Rows() == m.Rows(), "(Mat::!=) matrix rows don't match");
BartJanssens 0:785cff1e5a7c 307
BartJanssens 0:785cff1e5a7c 308 Int i;
BartJanssens 0:785cff1e5a7c 309
BartJanssens 0:785cff1e5a7c 310 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 311 if (m[i] != n[i])
BartJanssens 0:785cff1e5a7c 312 return(1);
BartJanssens 0:785cff1e5a7c 313
BartJanssens 0:785cff1e5a7c 314 return(0);
BartJanssens 0:785cff1e5a7c 315 }
BartJanssens 0:785cff1e5a7c 316
BartJanssens 0:785cff1e5a7c 317
BartJanssens 0:785cff1e5a7c 318 // --- Mat Arithmetic Operators -----------------------------------------------
BartJanssens 0:785cff1e5a7c 319
BartJanssens 0:785cff1e5a7c 320
BartJanssens 0:785cff1e5a7c 321 Mat operator + (const Mat &m, const Mat &n)
BartJanssens 0:785cff1e5a7c 322 {
BartJanssens 0:785cff1e5a7c 323 Assert(n.Rows() == m.Rows(), "(Mat::+) matrix rows don't match");
BartJanssens 0:785cff1e5a7c 324
BartJanssens 0:785cff1e5a7c 325 Mat result(m.Rows(), m.Cols());
BartJanssens 0:785cff1e5a7c 326 Int i;
BartJanssens 0:785cff1e5a7c 327
BartJanssens 0:785cff1e5a7c 328 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 329 result[i] = m[i] + n[i];
BartJanssens 0:785cff1e5a7c 330
BartJanssens 0:785cff1e5a7c 331 return(result);
BartJanssens 0:785cff1e5a7c 332 }
BartJanssens 0:785cff1e5a7c 333
BartJanssens 0:785cff1e5a7c 334 Mat operator - (const Mat &m, const Mat &n)
BartJanssens 0:785cff1e5a7c 335 {
BartJanssens 0:785cff1e5a7c 336 Assert(n.Rows() == m.Rows(), "(Mat::-) matrix rows don't match");
BartJanssens 0:785cff1e5a7c 337
BartJanssens 0:785cff1e5a7c 338 Mat result(m.Rows(), m.Cols());
BartJanssens 0:785cff1e5a7c 339 Int i;
BartJanssens 0:785cff1e5a7c 340
BartJanssens 0:785cff1e5a7c 341 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 342 result[i] = m[i] - n[i];
BartJanssens 0:785cff1e5a7c 343
BartJanssens 0:785cff1e5a7c 344 return(result);
BartJanssens 0:785cff1e5a7c 345 }
BartJanssens 0:785cff1e5a7c 346
BartJanssens 0:785cff1e5a7c 347 Mat operator - (const Mat &m)
BartJanssens 0:785cff1e5a7c 348 {
BartJanssens 0:785cff1e5a7c 349 Mat result(m.Rows(), m.Cols());
BartJanssens 0:785cff1e5a7c 350 Int i;
BartJanssens 0:785cff1e5a7c 351
BartJanssens 0:785cff1e5a7c 352 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 353 result[i] = -m[i];
BartJanssens 0:785cff1e5a7c 354
BartJanssens 0:785cff1e5a7c 355 return(result);
BartJanssens 0:785cff1e5a7c 356 }
BartJanssens 0:785cff1e5a7c 357
BartJanssens 0:785cff1e5a7c 358 Mat operator * (const Mat &m, const Mat &n)
BartJanssens 0:785cff1e5a7c 359 {
BartJanssens 0:785cff1e5a7c 360 Assert(m.Cols() == n.Rows(), "(Mat::*m) matrix cols don't match");
BartJanssens 0:785cff1e5a7c 361
BartJanssens 0:785cff1e5a7c 362 Mat result(m.Rows(), n.Cols());
BartJanssens 0:785cff1e5a7c 363 Int i;
BartJanssens 0:785cff1e5a7c 364
BartJanssens 0:785cff1e5a7c 365 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 366 result[i] = m[i] * n;
BartJanssens 0:785cff1e5a7c 367
BartJanssens 0:785cff1e5a7c 368 return(result);
BartJanssens 0:785cff1e5a7c 369 }
BartJanssens 0:785cff1e5a7c 370
BartJanssens 0:785cff1e5a7c 371 Vec operator * (const Mat &m, const Vec &v)
BartJanssens 0:785cff1e5a7c 372 {
BartJanssens 0:785cff1e5a7c 373 Assert(m.Cols() == v.Elts(), "(Mat::*v) matrix and vector sizes don't match");
BartJanssens 0:785cff1e5a7c 374
BartJanssens 0:785cff1e5a7c 375 Int i;
BartJanssens 0:785cff1e5a7c 376 Vec result(m.Rows());
BartJanssens 0:785cff1e5a7c 377
BartJanssens 0:785cff1e5a7c 378 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 379 result[i] = dot(v, m[i]);
BartJanssens 0:785cff1e5a7c 380
BartJanssens 0:785cff1e5a7c 381 return(result);
BartJanssens 0:785cff1e5a7c 382 }
BartJanssens 0:785cff1e5a7c 383
BartJanssens 0:785cff1e5a7c 384 Mat operator * (const Mat &m, Real s)
BartJanssens 0:785cff1e5a7c 385 {
BartJanssens 0:785cff1e5a7c 386 Int i;
BartJanssens 0:785cff1e5a7c 387 Mat result(m.Rows(), m.Cols());
BartJanssens 0:785cff1e5a7c 388
BartJanssens 0:785cff1e5a7c 389 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 390 result[i] = m[i] * s;
BartJanssens 0:785cff1e5a7c 391
BartJanssens 0:785cff1e5a7c 392 return(result);
BartJanssens 0:785cff1e5a7c 393 }
BartJanssens 0:785cff1e5a7c 394
BartJanssens 0:785cff1e5a7c 395 Mat operator / (const Mat &m, Real s)
BartJanssens 0:785cff1e5a7c 396 {
BartJanssens 0:785cff1e5a7c 397 Int i;
BartJanssens 0:785cff1e5a7c 398 Mat result(m.Rows(), m.Cols());
BartJanssens 0:785cff1e5a7c 399
BartJanssens 0:785cff1e5a7c 400 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 401 result[i] = m[i] / s;
BartJanssens 0:785cff1e5a7c 402
BartJanssens 0:785cff1e5a7c 403 return(result);
BartJanssens 0:785cff1e5a7c 404 }
BartJanssens 0:785cff1e5a7c 405
BartJanssens 0:785cff1e5a7c 406
BartJanssens 0:785cff1e5a7c 407 // --- Mat Mat-Vec Functions --------------------------------------------------
BartJanssens 0:785cff1e5a7c 408
BartJanssens 0:785cff1e5a7c 409
BartJanssens 0:785cff1e5a7c 410 Vec operator * (const Vec &v, const Mat &m) // v * m
BartJanssens 0:785cff1e5a7c 411 {
BartJanssens 0:785cff1e5a7c 412 Assert(v.Elts() == m.Rows(), "(Mat::v*m) vector/matrix sizes don't match");
BartJanssens 0:785cff1e5a7c 413
BartJanssens 0:785cff1e5a7c 414 Vec result(m.Cols(), vl_zero);
BartJanssens 0:785cff1e5a7c 415 Int i;
BartJanssens 0:785cff1e5a7c 416
BartJanssens 0:785cff1e5a7c 417 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 418 result += m[i] * v[i];
BartJanssens 0:785cff1e5a7c 419
BartJanssens 0:785cff1e5a7c 420 return(result);
BartJanssens 0:785cff1e5a7c 421 }
BartJanssens 0:785cff1e5a7c 422
BartJanssens 0:785cff1e5a7c 423
BartJanssens 0:785cff1e5a7c 424 // --- Mat Special Functions --------------------------------------------------
BartJanssens 0:785cff1e5a7c 425
BartJanssens 0:785cff1e5a7c 426
BartJanssens 0:785cff1e5a7c 427 Mat trans(const Mat &m)
BartJanssens 0:785cff1e5a7c 428 {
BartJanssens 0:785cff1e5a7c 429 Int i,j;
BartJanssens 0:785cff1e5a7c 430 Mat result(m.Cols(), m.Rows());
BartJanssens 0:785cff1e5a7c 431
BartJanssens 0:785cff1e5a7c 432 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 433 for (j = 0; j < m.Cols(); j++)
BartJanssens 0:785cff1e5a7c 434 result.Elt(j,i) = m.Elt(i,j);
BartJanssens 0:785cff1e5a7c 435
BartJanssens 0:785cff1e5a7c 436 return(result);
BartJanssens 0:785cff1e5a7c 437 }
BartJanssens 0:785cff1e5a7c 438
BartJanssens 0:785cff1e5a7c 439 Real trace(const Mat &m)
BartJanssens 0:785cff1e5a7c 440 {
BartJanssens 0:785cff1e5a7c 441 Int i;
BartJanssens 0:785cff1e5a7c 442 Real result = vl_0;
BartJanssens 0:785cff1e5a7c 443
BartJanssens 0:785cff1e5a7c 444 for (i = 0; i < m.Rows(); i++)
BartJanssens 0:785cff1e5a7c 445 result += m.Elt(i,i);
BartJanssens 0:785cff1e5a7c 446
BartJanssens 0:785cff1e5a7c 447 return(result);
BartJanssens 0:785cff1e5a7c 448 }
BartJanssens 0:785cff1e5a7c 449
BartJanssens 0:785cff1e5a7c 450 Mat &Mat::Clamp(Real fuzz)
BartJanssens 0:785cff1e5a7c 451 // clamps all values of the matrix with a magnitude
BartJanssens 0:785cff1e5a7c 452 // smaller than fuzz to zero.
BartJanssens 0:785cff1e5a7c 453 {
BartJanssens 0:785cff1e5a7c 454 Int i;
BartJanssens 0:785cff1e5a7c 455
BartJanssens 0:785cff1e5a7c 456 for (i = 0; i < Rows(); i++)
BartJanssens 0:785cff1e5a7c 457 SELF[i].Clamp(fuzz);
BartJanssens 0:785cff1e5a7c 458
BartJanssens 0:785cff1e5a7c 459 return(SELF);
BartJanssens 0:785cff1e5a7c 460 }
BartJanssens 0:785cff1e5a7c 461
BartJanssens 0:785cff1e5a7c 462 Mat &Mat::Clamp()
BartJanssens 0:785cff1e5a7c 463 {
BartJanssens 0:785cff1e5a7c 464 return(Clamp(1e-7));
BartJanssens 0:785cff1e5a7c 465 }
BartJanssens 0:785cff1e5a7c 466
BartJanssens 0:785cff1e5a7c 467 Mat clamped(const Mat &m, Real fuzz)
BartJanssens 0:785cff1e5a7c 468 // clamps all values of the matrix with a magnitude
BartJanssens 0:785cff1e5a7c 469 // smaller than fuzz to zero.
BartJanssens 0:785cff1e5a7c 470 {
BartJanssens 0:785cff1e5a7c 471 Mat result(m);
BartJanssens 0:785cff1e5a7c 472
BartJanssens 0:785cff1e5a7c 473 return(result.Clamp(fuzz));
BartJanssens 0:785cff1e5a7c 474 }
BartJanssens 0:785cff1e5a7c 475
BartJanssens 0:785cff1e5a7c 476 Mat clamped(const Mat &m)
BartJanssens 0:785cff1e5a7c 477 {
BartJanssens 0:785cff1e5a7c 478 return(clamped(m, 1e-7));
BartJanssens 0:785cff1e5a7c 479 }
BartJanssens 0:785cff1e5a7c 480
BartJanssens 0:785cff1e5a7c 481 Mat oprod(const Vec &a, const Vec &b)
BartJanssens 0:785cff1e5a7c 482 // returns outerproduct of a and b: a * trans(b)
BartJanssens 0:785cff1e5a7c 483 {
BartJanssens 0:785cff1e5a7c 484 Mat result;
BartJanssens 0:785cff1e5a7c 485 Int i;
BartJanssens 0:785cff1e5a7c 486
BartJanssens 0:785cff1e5a7c 487 result.SetSize(a.Elts(), b.Elts());
BartJanssens 0:785cff1e5a7c 488 for (i = 0; i < a.Elts(); i++)
BartJanssens 0:785cff1e5a7c 489 result[i] = a[i] * b;
BartJanssens 0:785cff1e5a7c 490
BartJanssens 0:785cff1e5a7c 491 return(result);
BartJanssens 0:785cff1e5a7c 492 }
BartJanssens 0:785cff1e5a7c 493
BartJanssens 0:785cff1e5a7c 494
BartJanssens 0:785cff1e5a7c 495 // --- Mat Input & Output -----------------------------------------------------
BartJanssens 0:785cff1e5a7c 496
BartJanssens 0:785cff1e5a7c 497 void printMat(const Mat &m)
BartJanssens 0:785cff1e5a7c 498 {
BartJanssens 0:785cff1e5a7c 499 Int i;
BartJanssens 0:785cff1e5a7c 500
BartJanssens 0:785cff1e5a7c 501 printf("[");
BartJanssens 0:785cff1e5a7c 502 for (i = 0; i < m.Rows() - 1; i++){
BartJanssens 0:785cff1e5a7c 503 printVec(m[i]);
BartJanssens 0:785cff1e5a7c 504 printf("\r\n");
BartJanssens 0:785cff1e5a7c 505 }
BartJanssens 0:785cff1e5a7c 506 printVec(m[i]);
BartJanssens 0:785cff1e5a7c 507 printf("]\r\n");
BartJanssens 0:785cff1e5a7c 508 }
BartJanssens 0:785cff1e5a7c 509
BartJanssens 0:785cff1e5a7c 510 /*
BartJanssens 0:785cff1e5a7c 511 ostream &operator << (ostream &s, const Mat &m)
BartJanssens 0:785cff1e5a7c 512 {
BartJanssens 0:785cff1e5a7c 513 Int i, w = s.width();
BartJanssens 0:785cff1e5a7c 514
BartJanssens 0:785cff1e5a7c 515 s << '[';
BartJanssens 0:785cff1e5a7c 516 for (i = 0; i < m.Rows() - 1; i++)
BartJanssens 0:785cff1e5a7c 517 s << setw(w) << m[i] << "\r\n";
BartJanssens 0:785cff1e5a7c 518 s << setw(w) << m[i] << ']' << "\r\n";
BartJanssens 0:785cff1e5a7c 519 return(s);
BartJanssens 0:785cff1e5a7c 520 }
BartJanssens 0:785cff1e5a7c 521 */
BartJanssens 0:785cff1e5a7c 522
BartJanssens 0:785cff1e5a7c 523 inline Void CopyPartialMat(const Mat &m, Mat &n, Int numRows)
BartJanssens 0:785cff1e5a7c 524 {
BartJanssens 0:785cff1e5a7c 525 for (Int i = 0; i < numRows; i++)
BartJanssens 0:785cff1e5a7c 526 n[i] = m[i];
BartJanssens 0:785cff1e5a7c 527 }
BartJanssens 0:785cff1e5a7c 528
BartJanssens 0:785cff1e5a7c 529 /*
BartJanssens 0:785cff1e5a7c 530 istream &operator >> (istream &s, Mat &m)
BartJanssens 0:785cff1e5a7c 531 {
BartJanssens 0:785cff1e5a7c 532 Int size = 1;
BartJanssens 0:785cff1e5a7c 533 Char c;
BartJanssens 0:785cff1e5a7c 534 Mat inMat;
BartJanssens 0:785cff1e5a7c 535
BartJanssens 0:785cff1e5a7c 536 // Expected format: [row0 row1 row2 ...]
BartJanssens 0:785cff1e5a7c 537
BartJanssens 0:785cff1e5a7c 538 while (isspace(s.peek())) // chomp white space
BartJanssens 0:785cff1e5a7c 539 s.get(c);
BartJanssens 0:785cff1e5a7c 540
BartJanssens 0:785cff1e5a7c 541 if (s.peek() == '[')
BartJanssens 0:785cff1e5a7c 542 {
BartJanssens 0:785cff1e5a7c 543 Vec row;
BartJanssens 0:785cff1e5a7c 544
BartJanssens 0:785cff1e5a7c 545 s.get(c);
BartJanssens 0:785cff1e5a7c 546
BartJanssens 0:785cff1e5a7c 547 s >> row;
BartJanssens 0:785cff1e5a7c 548 inMat.SetSize(2 * row.Elts(), row.Elts());
BartJanssens 0:785cff1e5a7c 549 inMat[0] = row;
BartJanssens 0:785cff1e5a7c 550
BartJanssens 0:785cff1e5a7c 551 while (isspace(s.peek())) // chomp white space
BartJanssens 0:785cff1e5a7c 552 s.get(c);
BartJanssens 0:785cff1e5a7c 553
BartJanssens 0:785cff1e5a7c 554 while (s.peek() != ']') // resize if needed
BartJanssens 0:785cff1e5a7c 555 {
BartJanssens 0:785cff1e5a7c 556 if (size == inMat.Rows())
BartJanssens 0:785cff1e5a7c 557 {
BartJanssens 0:785cff1e5a7c 558 Mat holdMat(inMat);
BartJanssens 0:785cff1e5a7c 559
BartJanssens 0:785cff1e5a7c 560 inMat.SetSize(size * 2, inMat.Cols());
BartJanssens 0:785cff1e5a7c 561 CopyPartialMat(holdMat, inMat, size);
BartJanssens 0:785cff1e5a7c 562 }
BartJanssens 0:785cff1e5a7c 563
BartJanssens 0:785cff1e5a7c 564 s >> row; // read a row
BartJanssens 0:785cff1e5a7c 565 inMat[size++] = row;
BartJanssens 0:785cff1e5a7c 566
BartJanssens 0:785cff1e5a7c 567 if (!s)
BartJanssens 0:785cff1e5a7c 568 {
BartJanssens 0:785cff1e5a7c 569 _Warning("Couldn't read matrix row");
BartJanssens 0:785cff1e5a7c 570 return(s);
BartJanssens 0:785cff1e5a7c 571 }
BartJanssens 0:785cff1e5a7c 572
BartJanssens 0:785cff1e5a7c 573 while (isspace(s.peek())) // chomp white space
BartJanssens 0:785cff1e5a7c 574 s.get(c);
BartJanssens 0:785cff1e5a7c 575 }
BartJanssens 0:785cff1e5a7c 576 s.get(c);
BartJanssens 0:785cff1e5a7c 577 }
BartJanssens 0:785cff1e5a7c 578 else
BartJanssens 0:785cff1e5a7c 579 {
BartJanssens 0:785cff1e5a7c 580 s.clear(ios::failbit);
BartJanssens 0:785cff1e5a7c 581 _Warning("Error: Expected '[' while reading matrix");
BartJanssens 0:785cff1e5a7c 582 return(s);
BartJanssens 0:785cff1e5a7c 583 }
BartJanssens 0:785cff1e5a7c 584
BartJanssens 0:785cff1e5a7c 585 m.SetSize(size, inMat.Cols());
BartJanssens 0:785cff1e5a7c 586 CopyPartialMat(inMat, m, size);
BartJanssens 0:785cff1e5a7c 587
BartJanssens 0:785cff1e5a7c 588 return(s);
BartJanssens 0:785cff1e5a7c 589 }
BartJanssens 0:785cff1e5a7c 590 */
BartJanssens 0:785cff1e5a7c 591
BartJanssens 0:785cff1e5a7c 592 // --- Matrix Inversion -------------------------------------------------------
BartJanssens 0:785cff1e5a7c 593
BartJanssens 0:785cff1e5a7c 594 #if !defined(CL_CHECKING) && !defined(VL_CHECKING)
BartJanssens 0:785cff1e5a7c 595 // we #define away pAssertEps if it is not used, to avoid
BartJanssens 0:785cff1e5a7c 596 // compiler warnings.
BartJanssens 0:785cff1e5a7c 597 #define pAssertEps
BartJanssens 0:785cff1e5a7c 598 #endif
BartJanssens 0:785cff1e5a7c 599
BartJanssens 0:785cff1e5a7c 600 Mat inv(const Mat &m, Real *determinant, Real pAssertEps)
BartJanssens 0:785cff1e5a7c 601 // matrix inversion using Gaussian pivoting
BartJanssens 0:785cff1e5a7c 602 {
BartJanssens 0:785cff1e5a7c 603 Assert(m.IsSquare(), "(inv) Matrix not square");
BartJanssens 0:785cff1e5a7c 604
BartJanssens 0:785cff1e5a7c 605 Int i, j, k;
BartJanssens 0:785cff1e5a7c 606 Int n = m.Rows();
BartJanssens 0:785cff1e5a7c 607 Real t, pivot, det;
BartJanssens 0:785cff1e5a7c 608 Real max;
BartJanssens 0:785cff1e5a7c 609 Mat A(m);
BartJanssens 0:785cff1e5a7c 610 Mat B(n, n, vl_I);
BartJanssens 0:785cff1e5a7c 611
BartJanssens 0:785cff1e5a7c 612 // ---------- Forward elimination ---------- ------------------------------
BartJanssens 0:785cff1e5a7c 613
BartJanssens 0:785cff1e5a7c 614 det = vl_1;
BartJanssens 0:785cff1e5a7c 615
BartJanssens 0:785cff1e5a7c 616 for (i = 0; i < n; i++) // Eliminate in column i, below diag
BartJanssens 0:785cff1e5a7c 617 {
BartJanssens 0:785cff1e5a7c 618 max = -1.0;
BartJanssens 0:785cff1e5a7c 619
BartJanssens 0:785cff1e5a7c 620 for (k = i; k < n; k++) // Find a pivot for column i
BartJanssens 0:785cff1e5a7c 621 if (len(A[k][i]) > max)
BartJanssens 0:785cff1e5a7c 622 {
BartJanssens 0:785cff1e5a7c 623 max = len(A[k][i]);
BartJanssens 0:785cff1e5a7c 624 j = k;
BartJanssens 0:785cff1e5a7c 625 }
BartJanssens 0:785cff1e5a7c 626
BartJanssens 0:785cff1e5a7c 627 Assert(max > pAssertEps, "(inv) Matrix not invertible");
BartJanssens 0:785cff1e5a7c 628
BartJanssens 0:785cff1e5a7c 629 if (j != i) // Swap rows i and j
BartJanssens 0:785cff1e5a7c 630 {
BartJanssens 0:785cff1e5a7c 631 for (k = i; k < n; k++)
BartJanssens 0:785cff1e5a7c 632 Swap(A.Elt(i, k), A.Elt(j, k));
BartJanssens 0:785cff1e5a7c 633 for (k = 0; k < n; k++)
BartJanssens 0:785cff1e5a7c 634 Swap(B.Elt(i, k), B.Elt(j, k));
BartJanssens 0:785cff1e5a7c 635
BartJanssens 0:785cff1e5a7c 636 det = -det;
BartJanssens 0:785cff1e5a7c 637 }
BartJanssens 0:785cff1e5a7c 638
BartJanssens 0:785cff1e5a7c 639 pivot = A.Elt(i, i);
BartJanssens 0:785cff1e5a7c 640 Assert(abs(pivot) > pAssertEps, "(inv) Matrix not invertible");
BartJanssens 0:785cff1e5a7c 641 det *= pivot;
BartJanssens 0:785cff1e5a7c 642
BartJanssens 0:785cff1e5a7c 643 for (k = i + 1; k < n; k++) // Only do elements to the right of the pivot
BartJanssens 0:785cff1e5a7c 644 A.Elt(i, k) /= pivot;
BartJanssens 0:785cff1e5a7c 645
BartJanssens 0:785cff1e5a7c 646 for (k = 0; k < n; k++)
BartJanssens 0:785cff1e5a7c 647 B.Elt(i, k) /= pivot;
BartJanssens 0:785cff1e5a7c 648
BartJanssens 0:785cff1e5a7c 649 // We know that A(i, i) will be set to 1, so don't bother to do it
BartJanssens 0:785cff1e5a7c 650
BartJanssens 0:785cff1e5a7c 651 for (j = i + 1; j < n; j++)
BartJanssens 0:785cff1e5a7c 652 { // Eliminate in rows below i
BartJanssens 0:785cff1e5a7c 653 t = A.Elt(j, i); // We're gonna zero this guy
BartJanssens 0:785cff1e5a7c 654 for (k = i + 1; k < n; k++) // Subtract scaled row i from row j
BartJanssens 0:785cff1e5a7c 655 A.Elt(j, k) -= A.Elt(i, k) * t; // (Ignore k <= i, we know they're 0)
BartJanssens 0:785cff1e5a7c 656 for (k = 0; k < n; k++)
BartJanssens 0:785cff1e5a7c 657 B.Elt(j, k) -= B.Elt(i, k) * t;
BartJanssens 0:785cff1e5a7c 658 }
BartJanssens 0:785cff1e5a7c 659 }
BartJanssens 0:785cff1e5a7c 660
BartJanssens 0:785cff1e5a7c 661 // ---------- Backward elimination ---------- -----------------------------
BartJanssens 0:785cff1e5a7c 662
BartJanssens 0:785cff1e5a7c 663 for (i = n - 1; i > 0; i--) // Eliminate in column i, above diag
BartJanssens 0:785cff1e5a7c 664 {
BartJanssens 0:785cff1e5a7c 665 for (j = 0; j < i; j++) // Eliminate in rows above i
BartJanssens 0:785cff1e5a7c 666 {
BartJanssens 0:785cff1e5a7c 667 t = A.Elt(j, i); // We're gonna zero this guy
BartJanssens 0:785cff1e5a7c 668 for (k = 0; k < n; k++) // Subtract scaled row i from row j
BartJanssens 0:785cff1e5a7c 669 B.Elt(j, k) -= B.Elt(i, k) * t;
BartJanssens 0:785cff1e5a7c 670 }
BartJanssens 0:785cff1e5a7c 671 }
BartJanssens 0:785cff1e5a7c 672 if (determinant)
BartJanssens 0:785cff1e5a7c 673 *determinant = det;
BartJanssens 0:785cff1e5a7c 674 return(B);
BartJanssens 0:785cff1e5a7c 675 }
BartJanssens 0:785cff1e5a7c 676