Context Navigation

source: nanovis/trunk/newmat11/newmat.h @ 5590

Last change on this file since 5590 was 2096, checked in by ldelgass, 13 years ago
Normalize line endings, set eol-style to native on .cpp, .h files
Property svn:eol-style set to `native`
File size: 84.1 KB

Line
1	/// \defgroup newmat Newmat matrix manipulation library
2	///@{
3
4	/// \file newmat.h
5	/// Definition file for matrix library.
6
7	// Copyright (C) 2004: R B Davies
8
9	#ifndef NEWMAT_LIB
10	#define NEWMAT_LIB 0
11
12	#include "include.h"
13
14	#include "myexcept.h"
15
16
17	#ifdef use_namespace
18	namespace NEWMAT { using namespace RBD_COMMON; }
19	namespace RBD_LIBRARIES { using namespace NEWMAT; }
20	namespace NEWMAT {
21	#endif
22
23	//#define DO_REPORT // to activate REPORT
24
25	#ifdef NO_LONG_NAMES
26	#define UpperTriangularMatrix UTMatrix
27	#define LowerTriangularMatrix LTMatrix
28	#define SymmetricMatrix SMatrix
29	#define DiagonalMatrix DMatrix
30	#define BandMatrix BMatrix
31	#define UpperBandMatrix UBMatrix
32	#define LowerBandMatrix LBMatrix
33	#define SymmetricBandMatrix SBMatrix
34	#define BandLUMatrix BLUMatrix
35	#endif
36
37	// ************************ general utilities **************************/
38
39	class GeneralMatrix; // defined later
40	class BaseMatrix; // defined later
41	class MatrixInput; // defined later
42
43	void MatrixErrorNoSpace(const void*); ///< test for allocation fails
44
45	/// Return from LogDeterminant function.
46	/// Members are the log of the absolute value and the sign (+1, -1 or 0)
47	class LogAndSign
48	{
49	Real log_val;
50	int sign_val;
51	public:
52	LogAndSign() { log_val=0.0; sign_val=1; }
53	LogAndSign(Real);
54	void operator*=(Real); ///< multiply by a real
55	void pow_eq(int k); ///< raise to power of k
56	void PowEq(int k) { pow_eq(k); }
57	void ChangeSign() { sign_val = -sign_val; }
58	void change_sign() { sign_val = -sign_val; } ///< change sign
59	Real LogValue() const { return log_val; }
60	Real log_value() const { return log_val; } ///< log of the absolute value
61	int Sign() const { return sign_val; }
62	int sign() const { return sign_val; } ///< sign of the value
63	Real value() const; ///< the value
64	Real Value() const { return value(); }
65	FREE_CHECK(LogAndSign)
66	};
67
68	// the following class is for counting the number of times a piece of code
69	// is executed. It is used for locating any code not executed by test
70	// routines. Use turbo GREP locate all places this code is called and
71	// check which ones are not accessed.
72	// Somewhat implementation dependent as it relies on "cout" still being
73	// present when ExeCounter objects are destructed.
74
75	#ifdef DO_REPORT
76
77	class ExeCounter
78	{
79	int line; // code line number
80	int fileid; // file identifier
81	long nexe; // number of executions
82	static int nreports; // number of reports
83	public:
84	ExeCounter(int,int);
85	void operator++() { nexe++; }
86	~ExeCounter(); // prints out reports
87	};
88
89	#endif
90
91
92	// ************************ class MatrixType ***************************
93
94	/// Find the type of a matrix resulting from matrix operations.
95	/// Also identify what conversions are permissible.
96	/// This class must be updated when new matrix types are added.
97
98	class MatrixType
99	{
100	public:
101	enum Attribute { Valid = 1,
102	Diagonal = 2, // order of these is important
103	Symmetric = 4,
104	Band = 8,
105	Lower = 16,
106	Upper = 32,
107	Square = 64,
108	Skew = 128,
109	LUDeco = 256,
110	Ones = 512 };
111
112	enum { US = 0,
113	UT = Valid + Upper + Square,
114	LT = Valid + Lower + Square,
115	Rt = Valid,
116	Sq = Valid + Square,
117	Sm = Valid + Symmetric + Square,
118	Sk = Valid + Skew + Square,
119	Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric
120	+ Square,
121	Id = Valid + Diagonal + Band + Lower + Upper + Symmetric
122	+ Square + Ones,
123	RV = Valid, // do not separate out
124	CV = Valid, // vectors
125	BM = Valid + Band + Square,
126	UB = Valid + Band + Upper + Square,
127	LB = Valid + Band + Lower + Square,
128	SB = Valid + Band + Symmetric + Square,
129	KB = Valid + Band + Skew + Square,
130	Ct = Valid + LUDeco + Square,
131	BC = Valid + Band + LUDeco + Square,
132	Mask = ~Square
133	};
134
135
136	static int nTypes() { return 13; } // number of different types
137	// exclude Ct, US, BC
138	public:
139	int attribute;
140	bool DataLossOK; // true if data loss is OK when
141	// this represents a destination
142	public:
143	MatrixType () : DataLossOK(false) {}
144	MatrixType (int i) : attribute(i), DataLossOK(false) {}
145	MatrixType (int i, bool dlok) : attribute(i), DataLossOK(dlok) {}
146	MatrixType (const MatrixType& mt)
147	: attribute(mt.attribute), DataLossOK(mt.DataLossOK) {}
148	void operator=(const MatrixType& mt)
149	{ attribute = mt.attribute; DataLossOK = mt.DataLossOK; }
150	void SetDataLossOK() { DataLossOK = true; }
151	int operator+() const { return attribute; }
152	MatrixType operator+(MatrixType mt) const
153	{ return MatrixType(attribute & mt.attribute); }
154	MatrixType operator*(const MatrixType&) const;
155	MatrixType SP(const MatrixType&) const;
156	MatrixType KP(const MatrixType&) const;
157	MatrixType operator\|(const MatrixType& mt) const
158	{ return MatrixType(attribute & mt.attribute & Valid); }
159	MatrixType operator&(const MatrixType& mt) const
160	{ return MatrixType(attribute & mt.attribute & Valid); }
161	bool operator>=(MatrixType mt) const
162	{ return ( attribute & ~mt.attribute & Mask ) == 0; }
163	bool operator<(MatrixType mt) const // for MS Visual C++ 4
164	{ return ( attribute & ~mt.attribute & Mask ) != 0; }
165	bool operator==(MatrixType t) const
166	{ return (attribute == t.attribute); }
167	bool operator!=(MatrixType t) const
168	{ return (attribute != t.attribute); }
169	bool operator!() const { return (attribute & Valid) == 0; }
170	MatrixType i() const; ///< type of inverse
171	MatrixType t() const; ///< type of transpose
172	MatrixType AddEqualEl() const ///< add constant to matrix
173	{ return MatrixType(attribute & (Valid + Symmetric + Square)); }
174	MatrixType MultRHS() const; ///< type for rhs of multiply
175	MatrixType sub() const ///< type of submatrix
176	{ return MatrixType(attribute & Valid); }
177	MatrixType ssub() const ///< type of sym submatrix
178	{ return MatrixType(attribute); } // not for selection matrix
179	GeneralMatrix* New() const; ///< new matrix of given type
180	GeneralMatrix* New(int,int,BaseMatrix*) const;
181	///< new matrix of given type
182	const char* value() const; ///< type as char string
183	const char* Value() const { return value(); }
184	friend bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
185	friend bool Compare(const MatrixType&, MatrixType&);
186	///< compare and check conversion
187	bool is_band() const { return (attribute & Band) != 0; }
188	bool is_diagonal() const { return (attribute & Diagonal) != 0; }
189	bool is_symmetric() const { return (attribute & Symmetric) != 0; }
190	bool CannotConvert() const { return (attribute & LUDeco) != 0; }
191	// used by operator==
192	FREE_CHECK(MatrixType)
193	};
194
195
196	// ********************* class MatrixBandWidth *********************/
197
198	///Upper and lower bandwidths of a matrix.
199	///That is number of diagonals strictly above or below main diagonal,
200	///e.g. diagonal matrix has 0 upper and lower bandwiths.
201	///-1 means the matrix may have the maximum bandwidth.
202	class MatrixBandWidth
203	{
204	public:
205	int lower_val;
206	int upper_val;
207	MatrixBandWidth(const int l, const int u) : lower_val(l), upper_val(u) {}
208	MatrixBandWidth(const int i) : lower_val(i), upper_val(i) {}
209	MatrixBandWidth operator+(const MatrixBandWidth&) const;
210	MatrixBandWidth operator*(const MatrixBandWidth&) const;
211	MatrixBandWidth minimum(const MatrixBandWidth&) const;
212	MatrixBandWidth t() const { return MatrixBandWidth(upper_val,lower_val); }
213	bool operator==(const MatrixBandWidth& bw) const
214	{ return (lower_val == bw.lower_val) && (upper_val == bw.upper_val); }
215	bool operator!=(const MatrixBandWidth& bw) const { return !operator==(bw); }
216	int Upper() const { return upper_val; }
217	int upper() const { return upper_val; }
218	int Lower() const { return lower_val; }
219	int lower() const { return lower_val; }
220	FREE_CHECK(MatrixBandWidth)
221	};
222
223
224	// ******************* Array length specifier **********************/
225
226	/// This class is used to avoid constructors such as
227	/// ColumnVector(int) being used for conversions.
228	/// Eventually this should be replaced by the use of the keyword "explicit".
229
230	class ArrayLengthSpecifier
231	{
232	int v;
233	public:
234	int Value() const { return v; }
235	int value() const { return v; }
236	ArrayLengthSpecifier(int l) : v(l) {}
237	};
238
239	// *********************** Matrix routines *************************/
240
241
242	class MatrixRowCol; // defined later
243	class MatrixRow;
244	class MatrixCol;
245	class MatrixColX;
246
247	class GeneralMatrix; // defined later
248	class AddedMatrix;
249	class MultipliedMatrix;
250	class SubtractedMatrix;
251	class SPMatrix;
252	class KPMatrix;
253	class ConcatenatedMatrix;
254	class StackedMatrix;
255	class SolvedMatrix;
256	class ShiftedMatrix;
257	class NegShiftedMatrix;
258	class ScaledMatrix;
259	class TransposedMatrix;
260	class ReversedMatrix;
261	class NegatedMatrix;
262	class InvertedMatrix;
263	class RowedMatrix;
264	class ColedMatrix;
265	class DiagedMatrix;
266	class MatedMatrix;
267	class GetSubMatrix;
268	class ReturnMatrix;
269	class Matrix;
270	class SquareMatrix;
271	class nricMatrix;
272	class RowVector;
273	class ColumnVector;
274	class SymmetricMatrix;
275	class UpperTriangularMatrix;
276	class LowerTriangularMatrix;
277	class DiagonalMatrix;
278	class CroutMatrix;
279	class BandMatrix;
280	class LowerBandMatrix;
281	class UpperBandMatrix;
282	class SymmetricBandMatrix;
283	class LinearEquationSolver;
284	class GenericMatrix;
285
286
287	#define MatrixTypeUnSp 0
288	//static MatrixType MatrixTypeUnSp(MatrixType::US);
289	// // AT&T needs this
290
291	/// Base of the matrix classes.
292	class BaseMatrix : public Janitor
293	{
294	protected:
295	virtual int search(const BaseMatrix*) const = 0;
296	// count number of times matrix is referred to
297	public:
298	virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
299	// evaluate temporary
300	// for old version of G++
301	// virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
302	// GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
303	AddedMatrix operator+(const BaseMatrix&) const; // results of operations
304	MultipliedMatrix operator*(const BaseMatrix&) const;
305	SubtractedMatrix operator-(const BaseMatrix&) const;
306	ConcatenatedMatrix operator\|(const BaseMatrix&) const;
307	StackedMatrix operator&(const BaseMatrix&) const;
308	ShiftedMatrix operator+(Real) const;
309	ScaledMatrix operator*(Real) const;
310	ScaledMatrix operator/(Real) const;
311	ShiftedMatrix operator-(Real) const;
312	TransposedMatrix t() const;
313	// TransposedMatrix t;
314	NegatedMatrix operator-() const; // change sign of elements
315	ReversedMatrix reverse() const;
316	ReversedMatrix Reverse() const;
317	InvertedMatrix i() const;
318	// InvertedMatrix i;
319	RowedMatrix as_row() const;
320	RowedMatrix AsRow() const;
321	ColedMatrix as_column() const;
322	ColedMatrix AsColumn() const;
323	DiagedMatrix as_diagonal() const;
324	DiagedMatrix AsDiagonal() const;
325	MatedMatrix as_matrix(int,int) const;
326	MatedMatrix AsMatrix(int m, int n) const;
327	GetSubMatrix submatrix(int,int,int,int) const;
328	GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const;
329	GetSubMatrix sym_submatrix(int,int) const;
330	GetSubMatrix SymSubMatrix(int f, int l) const;
331	GetSubMatrix row(int) const;
332	GetSubMatrix rows(int,int) const;
333	GetSubMatrix column(int) const;
334	GetSubMatrix columns(int,int) const;
335	GetSubMatrix Row(int f) const;
336	GetSubMatrix Rows(int f, int l) const;
337	GetSubMatrix Column(int f) const;
338	GetSubMatrix Columns(int f, int l) const;
339	Real as_scalar() const; // conversion of 1 x 1 matrix
340	Real AsScalar() const;
341	virtual LogAndSign log_determinant() const;
342	LogAndSign LogDeterminant() const { return log_determinant(); }
343	Real determinant() const;
344	Real Determinant() const { return determinant(); }
345	virtual Real sum_square() const;
346	Real SumSquare() const { return sum_square(); }
347	Real norm_Frobenius() const;
348	Real norm_frobenius() const { return norm_Frobenius(); }
349	Real NormFrobenius() const { return norm_Frobenius(); }
350	virtual Real sum_absolute_value() const;
351	Real SumAbsoluteValue() const { return sum_absolute_value(); }
352	virtual Real sum() const;
353	virtual Real Sum() const { return sum(); }
354	virtual Real maximum_absolute_value() const;
355	Real MaximumAbsoluteValue() const { return maximum_absolute_value(); }
356	virtual Real maximum_absolute_value1(int& i) const;
357	Real MaximumAbsoluteValue1(int& i) const
358	{ return maximum_absolute_value1(i); }
359	virtual Real maximum_absolute_value2(int& i, int& j) const;
360	Real MaximumAbsoluteValue2(int& i, int& j) const
361	{ return maximum_absolute_value2(i,j); }
362	virtual Real minimum_absolute_value() const;
363	Real MinimumAbsoluteValue() const { return minimum_absolute_value(); }
364	virtual Real minimum_absolute_value1(int& i) const;
365	Real MinimumAbsoluteValue1(int& i) const
366	{ return minimum_absolute_value1(i); }
367	virtual Real minimum_absolute_value2(int& i, int& j) const;
368	Real MinimumAbsoluteValue2(int& i, int& j) const
369	{ return minimum_absolute_value2(i,j); }
370	virtual Real maximum() const;
371	Real Maximum() const { return maximum(); }
372	virtual Real maximum1(int& i) const;
373	Real Maximum1(int& i) const { return maximum1(i); }
374	virtual Real maximum2(int& i, int& j) const;
375	Real Maximum2(int& i, int& j) const { return maximum2(i,j); }
376	virtual Real minimum() const;
377	Real Minimum() const { return minimum(); }
378	virtual Real minimum1(int& i) const;
379	Real Minimum1(int& i) const { return minimum1(i); }
380	virtual Real minimum2(int& i, int& j) const;
381	Real Minimum2(int& i, int& j) const { return minimum2(i,j); }
382	virtual Real trace() const;
383	Real Trace() const { return trace(); }
384	Real norm1() const;
385	Real Norm1() const { return norm1(); }
386	Real norm_infinity() const;
387	Real NormInfinity() const { return norm_infinity(); }
388	virtual MatrixBandWidth bandwidth() const; // bandwidths of band matrix
389	virtual MatrixBandWidth BandWidth() const { return bandwidth(); }
390	void IEQND() const; // called by ineq. ops
391	ReturnMatrix sum_square_columns() const;
392	ReturnMatrix sum_square_rows() const;
393	ReturnMatrix sum_columns() const;
394	ReturnMatrix sum_rows() const;
395	virtual void cleanup() {}
396	void CleanUp() { cleanup(); }
397
398	// virtual ReturnMatrix Reverse() const; // reverse order of elements
399	//protected:
400	// BaseMatrix() : t(this), i(this) {}
401
402	friend class GeneralMatrix;
403	friend class Matrix;
404	friend class SquareMatrix;
405	friend class nricMatrix;
406	friend class RowVector;
407	friend class ColumnVector;
408	friend class SymmetricMatrix;
409	friend class UpperTriangularMatrix;
410	friend class LowerTriangularMatrix;
411	friend class DiagonalMatrix;
412	friend class CroutMatrix;
413	friend class BandMatrix;
414	friend class LowerBandMatrix;
415	friend class UpperBandMatrix;
416	friend class SymmetricBandMatrix;
417	friend class AddedMatrix;
418	friend class MultipliedMatrix;
419	friend class SubtractedMatrix;
420	friend class SPMatrix;
421	friend class KPMatrix;
422	friend class ConcatenatedMatrix;
423	friend class StackedMatrix;
424	friend class SolvedMatrix;
425	friend class ShiftedMatrix;
426	friend class NegShiftedMatrix;
427	friend class ScaledMatrix;
428	friend class TransposedMatrix;
429	friend class ReversedMatrix;
430	friend class NegatedMatrix;
431	friend class InvertedMatrix;
432	friend class RowedMatrix;
433	friend class ColedMatrix;
434	friend class DiagedMatrix;
435	friend class MatedMatrix;
436	friend class GetSubMatrix;
437	friend class ReturnMatrix;
438	friend class LinearEquationSolver;
439	friend class GenericMatrix;
440	NEW_DELETE(BaseMatrix)
441	};
442
443
444	// *************************** working classes ************************/
445
446	/// The classes for matrices that can contain data are derived from this.
447	class GeneralMatrix : public BaseMatrix // declarable matrix types
448	{
449	virtual GeneralMatrix* Image() const; // copy of matrix
450	protected:
451	int tag_val; // shows whether can reuse
452	int nrows_val, ncols_val; // dimensions
453	int storage; // total store required
454	Real* store; // point to store (0=not set)
455	GeneralMatrix(); // initialise with no store
456	GeneralMatrix(ArrayLengthSpecifier); // constructor getting store
457	void Add(GeneralMatrix*, Real); // sum of GM and Real
458	void Add(Real); // add Real to this
459	void NegAdd(GeneralMatrix*, Real); // Real - GM
460	void NegAdd(Real); // this = this - Real
461	void Multiply(GeneralMatrix*, Real); // product of GM and Real
462	void Multiply(Real); // multiply this by Real
463	void Negate(GeneralMatrix*); // change sign
464	void Negate(); // change sign
465	void ReverseElements(); // internal reverse of elements
466	void ReverseElements(GeneralMatrix*); // reverse order of elements
467	void operator=(Real); // set matrix to constant
468	Real* GetStore(); // get store or copy
469	GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
470	// temporarily access store
471	void GetMatrix(const GeneralMatrix*); // used by = and initialise
472	void Eq(const BaseMatrix&, MatrixType); // used by =
473	void Eq(const GeneralMatrix&); // version with no conversion
474	void Eq(const BaseMatrix&, MatrixType, bool);// used by <<
475	void Eq2(const BaseMatrix&, MatrixType); // cut down version of Eq
476	int search(const BaseMatrix*) const;
477	virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
478	void CheckConversion(const BaseMatrix&); // check conversion OK
479	void resize(int, int, int); // change dimensions
480	virtual short SimpleAddOK(const GeneralMatrix*) { return 0; }
481	// see bandmat.cpp for explanation
482	virtual void MiniCleanUp()
483	{ store = 0; storage = 0; nrows_val = 0; ncols_val = 0; tag_val = -1;}
484	// CleanUp when the data array has already been deleted
485	void PlusEqual(const GeneralMatrix& gm);
486	void MinusEqual(const GeneralMatrix& gm);
487	void PlusEqual(Real f);
488	void MinusEqual(Real f);
489	void swap(GeneralMatrix& gm); // swap values
490	public:
491	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
492	virtual MatrixType type() const = 0; // type of a matrix
493	MatrixType Type() const { return type(); }
494	int Nrows() const { return nrows_val; } // get dimensions
495	int Ncols() const { return ncols_val; }
496	int Storage() const { return storage; }
497	Real* Store() const { return store; }
498	// updated names
499	int nrows() const { return nrows_val; } // get dimensions
500	int ncols() const { return ncols_val; }
501	int size() const { return storage; }
502	Real* data() { return store; }
503	const Real* data() const { return store; }
504	const Real* const_data() const { return store; }
505	virtual ~GeneralMatrix(); // delete store if set
506	void tDelete(); // delete if tag_val permits
507	bool reuse(); // true if tag_val allows reuse
508	void protect() { tag_val=-1; } // cannot delete or reuse
509	void Protect() { tag_val=-1; } // cannot delete or reuse
510	int tag() const { return tag_val; }
511	int Tag() const { return tag_val; }
512	bool is_zero() const; // test matrix has all zeros
513	bool IsZero() const { return is_zero(); } // test matrix has all zeros
514	void Release() { tag_val=1; } // del store after next use
515	void Release(int t) { tag_val=t; } // del store after t accesses
516	void ReleaseAndDelete() { tag_val=0; } // delete matrix after use
517	void release() { tag_val=1; } // del store after next use
518	void release(int t) { tag_val=t; } // del store after t accesses
519	void release_and_delete() { tag_val=0; } // delete matrix after use
520	void operator<<(const double*); // assignment from an array
521	void operator<<(const float*); // assignment from an array
522	void operator<<(const int*); // assignment from an array
523	void operator<<(const BaseMatrix& X)
524	{ Eq(X,this->type(),true); } // = without checking type
525	void inject(const GeneralMatrix&); // copy stored els only
526	void Inject(const GeneralMatrix& GM) { inject(GM); }
527	void operator+=(const BaseMatrix&);
528	void operator-=(const BaseMatrix&);
529	void operator*=(const BaseMatrix&);
530	void operator\|=(const BaseMatrix&);
531	void operator&=(const BaseMatrix&);
532	void operator+=(Real);
533	void operator-=(Real r) { operator+=(-r); }
534	void operator*=(Real);
535	void operator/=(Real r) { operator*=(1.0/r); }
536	virtual GeneralMatrix* MakeSolver(); // for solving
537	virtual void Solver(MatrixColX&, const MatrixColX&) {}
538	virtual void GetRow(MatrixRowCol&) = 0; // Get matrix row
539	virtual void RestoreRow(MatrixRowCol&) {} // Restore matrix row
540	virtual void NextRow(MatrixRowCol&); // Go to next row
541	virtual void GetCol(MatrixRowCol&) = 0; // Get matrix col
542	virtual void GetCol(MatrixColX&) = 0; // Get matrix col
543	virtual void RestoreCol(MatrixRowCol&) {} // Restore matrix col
544	virtual void RestoreCol(MatrixColX&) {} // Restore matrix col
545	virtual void NextCol(MatrixRowCol&); // Go to next col
546	virtual void NextCol(MatrixColX&); // Go to next col
547	Real sum_square() const;
548	Real sum_absolute_value() const;
549	Real sum() const;
550	Real maximum_absolute_value1(int& i) const;
551	Real minimum_absolute_value1(int& i) const;
552	Real maximum1(int& i) const;
553	Real minimum1(int& i) const;
554	Real maximum_absolute_value() const;
555	Real maximum_absolute_value2(int& i, int& j) const;
556	Real minimum_absolute_value() const;
557	Real minimum_absolute_value2(int& i, int& j) const;
558	Real maximum() const;
559	Real maximum2(int& i, int& j) const;
560	Real minimum() const;
561	Real minimum2(int& i, int& j) const;
562	LogAndSign log_determinant() const;
563	virtual bool IsEqual(const GeneralMatrix&) const;
564	// same type, same values
565	void CheckStore() const; // check store is non-zero
566	virtual void SetParameters(const GeneralMatrix*) {}
567	// set parameters in GetMatrix
568	operator ReturnMatrix() const; // for building a ReturnMatrix
569	ReturnMatrix for_return() const;
570	ReturnMatrix ForReturn() const;
571	//virtual bool SameStorageType(const GeneralMatrix& A) const;
572	//virtual void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
573	//virtual void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
574	virtual void resize(const GeneralMatrix& A);
575	virtual void ReSize(const GeneralMatrix& A) { resize(A); }
576	MatrixInput operator<<(double); // for loading a list
577	MatrixInput operator<<(float); // for loading a list
578	MatrixInput operator<<(int f);
579	// ReturnMatrix Reverse() const; // reverse order of elements
580	void cleanup(); // to clear store
581
582	friend class Matrix;
583	friend class SquareMatrix;
584	friend class nricMatrix;
585	friend class SymmetricMatrix;
586	friend class UpperTriangularMatrix;
587	friend class LowerTriangularMatrix;
588	friend class DiagonalMatrix;
589	friend class CroutMatrix;
590	friend class RowVector;
591	friend class ColumnVector;
592	friend class BandMatrix;
593	friend class LowerBandMatrix;
594	friend class UpperBandMatrix;
595	friend class SymmetricBandMatrix;
596	friend class BaseMatrix;
597	friend class AddedMatrix;
598	friend class MultipliedMatrix;
599	friend class SubtractedMatrix;
600	friend class SPMatrix;
601	friend class KPMatrix;
602	friend class ConcatenatedMatrix;
603	friend class StackedMatrix;
604	friend class SolvedMatrix;
605	friend class ShiftedMatrix;
606	friend class NegShiftedMatrix;
607	friend class ScaledMatrix;
608	friend class TransposedMatrix;
609	friend class ReversedMatrix;
610	friend class NegatedMatrix;
611	friend class InvertedMatrix;
612	friend class RowedMatrix;
613	friend class ColedMatrix;
614	friend class DiagedMatrix;
615	friend class MatedMatrix;
616	friend class GetSubMatrix;
617	friend class ReturnMatrix;
618	friend class LinearEquationSolver;
619	friend class GenericMatrix;
620	NEW_DELETE(GeneralMatrix)
621	};
622
623
624	/// The usual rectangular matrix.
625	class Matrix : public GeneralMatrix
626	{
627	GeneralMatrix* Image() const; // copy of matrix
628	public:
629	Matrix() {}
630	~Matrix() {}
631	Matrix(int, int); // standard declaration
632	Matrix(const BaseMatrix&); // evaluate BaseMatrix
633	void operator=(const BaseMatrix&);
634	void operator=(Real f) { GeneralMatrix::operator=(f); }
635	void operator=(const Matrix& m) { Eq(m); }
636	MatrixType type() const;
637	Real& operator()(int, int); // access element
638	Real& element(int, int); // access element
639	Real operator()(int, int) const; // access element
640	Real element(int, int) const; // access element
641	#ifdef SETUP_C_SUBSCRIPTS
642	Real* operator[](int m) { return store+m*ncols_val; }
643	const Real* operator[](int m) const { return store+m*ncols_val; }
644	// following for Numerical Recipes in C++
645	Matrix(Real, int, int);
646	Matrix(const Real*, int, int);
647	#endif
648	Matrix(const Matrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
649	GeneralMatrix* MakeSolver();
650	Real trace() const;
651	void GetRow(MatrixRowCol&);
652	void GetCol(MatrixRowCol&);
653	void GetCol(MatrixColX&);
654	void RestoreCol(MatrixRowCol&);
655	void RestoreCol(MatrixColX&);
656	void NextRow(MatrixRowCol&);
657	void NextCol(MatrixRowCol&);
658	void NextCol(MatrixColX&);
659	virtual void resize(int,int); // change dimensions
660	// virtual so we will catch it being used in a vector called as a matrix
661	virtual void resize_keep(int,int);
662	virtual void ReSize(int m, int n) { resize(m, n); }
663	void resize(const GeneralMatrix& A);
664	void ReSize(const GeneralMatrix& A) { resize(A); }
665	Real maximum_absolute_value2(int& i, int& j) const;
666	Real minimum_absolute_value2(int& i, int& j) const;
667	Real maximum2(int& i, int& j) const;
668	Real minimum2(int& i, int& j) const;
669	void operator+=(const Matrix& M) { PlusEqual(M); }
670	void operator-=(const Matrix& M) { MinusEqual(M); }
671	void operator+=(Real f) { GeneralMatrix::Add(f); }
672	void operator-=(Real f) { GeneralMatrix::Add(-f); }
673	void swap(Matrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
674	friend Real dotproduct(const Matrix& A, const Matrix& B);
675	NEW_DELETE(Matrix)
676	};
677
678	/// Square matrix.
679	class SquareMatrix : public Matrix
680	{
681	GeneralMatrix* Image() const; // copy of matrix
682	public:
683	SquareMatrix() {}
684	~SquareMatrix() {}
685	SquareMatrix(ArrayLengthSpecifier); // standard declaration
686	SquareMatrix(const BaseMatrix&); // evaluate BaseMatrix
687	void operator=(const BaseMatrix&);
688	void operator=(Real f) { GeneralMatrix::operator=(f); }
689	void operator=(const SquareMatrix& m) { Eq(m); }
690	void operator=(const Matrix& m);
691	MatrixType type() const;
692	SquareMatrix(const SquareMatrix& gm) : Matrix() { GetMatrix(&gm); }
693	SquareMatrix(const Matrix& gm);
694	void resize(int); // change dimensions
695	void ReSize(int m) { resize(m); }
696	void resize_keep(int);
697	void resize_keep(int,int);
698	void resize(int,int); // change dimensions
699	void ReSize(int m, int n) { resize(m, n); }
700	void resize(const GeneralMatrix& A);
701	void ReSize(const GeneralMatrix& A) { resize(A); }
702	void operator+=(const Matrix& M) { PlusEqual(M); }
703	void operator-=(const Matrix& M) { MinusEqual(M); }
704	void operator+=(Real f) { GeneralMatrix::Add(f); }
705	void operator-=(Real f) { GeneralMatrix::Add(-f); }
706	void swap(SquareMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
707	NEW_DELETE(SquareMatrix)
708	};
709
710	/// Rectangular matrix for use with Numerical Recipes in C.
711	class nricMatrix : public Matrix
712	{
713	GeneralMatrix* Image() const; // copy of matrix
714	Real** row_pointer; // points to rows
715	void MakeRowPointer(); // build rowpointer
716	void DeleteRowPointer();
717	public:
718	nricMatrix() {}
719	nricMatrix(int m, int n) // standard declaration
720	: Matrix(m,n) { MakeRowPointer(); }
721	nricMatrix(const BaseMatrix& bm) // evaluate BaseMatrix
722	: Matrix(bm) { MakeRowPointer(); }
723	void operator=(const BaseMatrix& bm)
724	{ DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
725	void operator=(Real f) { GeneralMatrix::operator=(f); }
726	void operator=(const nricMatrix& m)
727	{ DeleteRowPointer(); Eq(m); MakeRowPointer(); }
728	void operator<<(const BaseMatrix& X)
729	{ DeleteRowPointer(); Eq(X,this->type(),true); MakeRowPointer(); }
730	nricMatrix(const nricMatrix& gm) : Matrix()
731	{ GetMatrix(&gm); MakeRowPointer(); }
732	void resize(int m, int n) // change dimensions
733	{ DeleteRowPointer(); Matrix::resize(m,n); MakeRowPointer(); }
734	void resize_keep(int m, int n) // change dimensions
735	{ DeleteRowPointer(); Matrix::resize_keep(m,n); MakeRowPointer(); }
736	void ReSize(int m, int n) // change dimensions
737	{ DeleteRowPointer(); Matrix::resize(m,n); MakeRowPointer(); }
738	void resize(const GeneralMatrix& A);
739	void ReSize(const GeneralMatrix& A) { resize(A); }
740	~nricMatrix() { DeleteRowPointer(); }
741	Real** nric() const { CheckStore(); return row_pointer-1; }
742	void cleanup(); // to clear store
743	void MiniCleanUp();
744	void operator+=(const Matrix& M) { PlusEqual(M); }
745	void operator-=(const Matrix& M) { MinusEqual(M); }
746	void operator+=(Real f) { GeneralMatrix::Add(f); }
747	void operator-=(Real f) { GeneralMatrix::Add(-f); }
748	void swap(nricMatrix& gm);
749	NEW_DELETE(nricMatrix)
750	};
751
752	/// Symmetric matrix.
753	class SymmetricMatrix : public GeneralMatrix
754	{
755	GeneralMatrix* Image() const; // copy of matrix
756	public:
757	SymmetricMatrix() {}
758	~SymmetricMatrix() {}
759	SymmetricMatrix(ArrayLengthSpecifier);
760	SymmetricMatrix(const BaseMatrix&);
761	void operator=(const BaseMatrix&);
762	void operator=(Real f) { GeneralMatrix::operator=(f); }
763	void operator=(const SymmetricMatrix& m) { Eq(m); }
764	Real& operator()(int, int); // access element
765	Real& element(int, int); // access element
766	Real operator()(int, int) const; // access element
767	Real element(int, int) const; // access element
768	#ifdef SETUP_C_SUBSCRIPTS
769	Real* operator[](int m) { return store+(m*(m+1))/2; }
770	const Real* operator[](int m) const { return store+(m*(m+1))/2; }
771	#endif
772	MatrixType type() const;
773	SymmetricMatrix(const SymmetricMatrix& gm)
774	: GeneralMatrix() { GetMatrix(&gm); }
775	Real sum_square() const;
776	Real sum_absolute_value() const;
777	Real sum() const;
778	Real trace() const;
779	void GetRow(MatrixRowCol&);
780	void GetCol(MatrixRowCol&);
781	void GetCol(MatrixColX&);
782	void RestoreCol(MatrixRowCol&) {}
783	void RestoreCol(MatrixColX&);
784	GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
785	void resize(int); // change dimensions
786	void ReSize(int m) { resize(m); }
787	void resize_keep(int);
788	void resize(const GeneralMatrix& A);
789	void ReSize(const GeneralMatrix& A) { resize(A); }
790	void operator+=(const SymmetricMatrix& M) { PlusEqual(M); }
791	void operator-=(const SymmetricMatrix& M) { MinusEqual(M); }
792	void operator+=(Real f) { GeneralMatrix::Add(f); }
793	void operator-=(Real f) { GeneralMatrix::Add(-f); }
794	void swap(SymmetricMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
795	NEW_DELETE(SymmetricMatrix)
796	};
797
798	/// Upper triangular matrix.
799	class UpperTriangularMatrix : public GeneralMatrix
800	{
801	GeneralMatrix* Image() const; // copy of matrix
802	public:
803	UpperTriangularMatrix() {}
804	~UpperTriangularMatrix() {}
805	UpperTriangularMatrix(ArrayLengthSpecifier);
806	void operator=(const BaseMatrix&);
807	void operator=(const UpperTriangularMatrix& m) { Eq(m); }
808	UpperTriangularMatrix(const BaseMatrix&);
809	UpperTriangularMatrix(const UpperTriangularMatrix& gm)
810	: GeneralMatrix() { GetMatrix(&gm); }
811	void operator=(Real f) { GeneralMatrix::operator=(f); }
812	Real& operator()(int, int); // access element
813	Real& element(int, int); // access element
814	Real operator()(int, int) const; // access element
815	Real element(int, int) const; // access element
816	#ifdef SETUP_C_SUBSCRIPTS
817	Real* operator[](int m) { return store+mncols_val-(m(m+1))/2; }
818	const Real* operator[](int m) const
819	{ return store+mncols_val-(m(m+1))/2; }
820	#endif
821	MatrixType type() const;
822	GeneralMatrix* MakeSolver() { return this; } // for solving
823	void Solver(MatrixColX&, const MatrixColX&);
824	LogAndSign log_determinant() const;
825	Real trace() const;
826	void GetRow(MatrixRowCol&);
827	void GetCol(MatrixRowCol&);
828	void GetCol(MatrixColX&);
829	void RestoreCol(MatrixRowCol&);
830	void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
831	void NextRow(MatrixRowCol&);
832	void resize(int); // change dimensions
833	void ReSize(int m) { resize(m); }
834	void resize(const GeneralMatrix& A);
835	void ReSize(const GeneralMatrix& A) { resize(A); }
836	void resize_keep(int);
837	MatrixBandWidth bandwidth() const;
838	void operator+=(const UpperTriangularMatrix& M) { PlusEqual(M); }
839	void operator-=(const UpperTriangularMatrix& M) { MinusEqual(M); }
840	void operator+=(Real f) { GeneralMatrix::operator+=(f); }
841	void operator-=(Real f) { GeneralMatrix::operator-=(f); }
842	void swap(UpperTriangularMatrix& gm)
843	{ GeneralMatrix::swap((GeneralMatrix&)gm); }
844	NEW_DELETE(UpperTriangularMatrix)
845	};
846
847	/// Lower triangular matrix.
848	class LowerTriangularMatrix : public GeneralMatrix
849	{
850	GeneralMatrix* Image() const; // copy of matrix
851	public:
852	LowerTriangularMatrix() {}
853	~LowerTriangularMatrix() {}
854	LowerTriangularMatrix(ArrayLengthSpecifier);
855	LowerTriangularMatrix(const LowerTriangularMatrix& gm)
856	: GeneralMatrix() { GetMatrix(&gm); }
857	LowerTriangularMatrix(const BaseMatrix& M);
858	void operator=(const BaseMatrix&);
859	void operator=(Real f) { GeneralMatrix::operator=(f); }
860	void operator=(const LowerTriangularMatrix& m) { Eq(m); }
861	Real& operator()(int, int); // access element
862	Real& element(int, int); // access element
863	Real operator()(int, int) const; // access element
864	Real element(int, int) const; // access element
865	#ifdef SETUP_C_SUBSCRIPTS
866	Real* operator[](int m) { return store+(m*(m+1))/2; }
867	const Real* operator[](int m) const { return store+(m*(m+1))/2; }
868	#endif
869	MatrixType type() const;
870	GeneralMatrix* MakeSolver() { return this; } // for solving
871	void Solver(MatrixColX&, const MatrixColX&);
872	LogAndSign log_determinant() const;
873	Real trace() const;
874	void GetRow(MatrixRowCol&);
875	void GetCol(MatrixRowCol&);
876	void GetCol(MatrixColX&);
877	void RestoreCol(MatrixRowCol&);
878	void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
879	void NextRow(MatrixRowCol&);
880	void resize(int); // change dimensions
881	void ReSize(int m) { resize(m); }
882	void resize_keep(int);
883	void resize(const GeneralMatrix& A);
884	void ReSize(const GeneralMatrix& A) { resize(A); }
885	MatrixBandWidth bandwidth() const;
886	void operator+=(const LowerTriangularMatrix& M) { PlusEqual(M); }
887	void operator-=(const LowerTriangularMatrix& M) { MinusEqual(M); }
888	void operator+=(Real f) { GeneralMatrix::operator+=(f); }
889	void operator-=(Real f) { GeneralMatrix::operator-=(f); }
890	void swap(LowerTriangularMatrix& gm)
891	{ GeneralMatrix::swap((GeneralMatrix&)gm); }
892	NEW_DELETE(LowerTriangularMatrix)
893	};
894
895	/// Diagonal matrix.
896	class DiagonalMatrix : public GeneralMatrix
897	{
898	GeneralMatrix* Image() const; // copy of matrix
899	public:
900	DiagonalMatrix() {}
901	~DiagonalMatrix() {}
902	DiagonalMatrix(ArrayLengthSpecifier);
903	DiagonalMatrix(const BaseMatrix&);
904	DiagonalMatrix(const DiagonalMatrix& gm)
905	: GeneralMatrix() { GetMatrix(&gm); }
906	void operator=(const BaseMatrix&);
907	void operator=(Real f) { GeneralMatrix::operator=(f); }
908	void operator=(const DiagonalMatrix& m) { Eq(m); }
909	Real& operator()(int, int); // access element
910	Real& operator()(int); // access element
911	Real operator()(int, int) const; // access element
912	Real operator()(int) const;
913	Real& element(int, int); // access element
914	Real& element(int); // access element
915	Real element(int, int) const; // access element
916	Real element(int) const; // access element
917	#ifdef SETUP_C_SUBSCRIPTS
918	Real& operator[](int m) { return store[m]; }
919	const Real& operator[](int m) const { return store[m]; }
920	#endif
921	MatrixType type() const;
922
923	LogAndSign log_determinant() const;
924	Real trace() const;
925	void GetRow(MatrixRowCol&);
926	void GetCol(MatrixRowCol&);
927	void GetCol(MatrixColX&);
928	void NextRow(MatrixRowCol&);
929	void NextCol(MatrixRowCol&);
930	void NextCol(MatrixColX&);
931	GeneralMatrix* MakeSolver() { return this; } // for solving
932	void Solver(MatrixColX&, const MatrixColX&);
933	GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
934	void resize(int); // change dimensions
935	void ReSize(int m) { resize(m); }
936	void resize_keep(int);
937	void resize(const GeneralMatrix& A);
938	void ReSize(const GeneralMatrix& A) { resize(A); }
939	Real* nric() const
940	{ CheckStore(); return store-1; } // for use by NRIC
941	MatrixBandWidth bandwidth() const;
942	// ReturnMatrix Reverse() const; // reverse order of elements
943	void operator+=(const DiagonalMatrix& M) { PlusEqual(M); }
944	void operator-=(const DiagonalMatrix& M) { MinusEqual(M); }
945	void operator+=(Real f) { GeneralMatrix::operator+=(f); }
946	void operator-=(Real f) { GeneralMatrix::operator-=(f); }
947	void swap(DiagonalMatrix& gm)
948	{ GeneralMatrix::swap((GeneralMatrix&)gm); }
949	NEW_DELETE(DiagonalMatrix)
950	};
951
952	/// Row vector.
953	class RowVector : public Matrix
954	{
955	GeneralMatrix* Image() const; // copy of matrix
956	public:
957	RowVector() { nrows_val = 1; }
958	~RowVector() {}
959	RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
960	RowVector(const BaseMatrix&);
961	RowVector(const RowVector& gm) : Matrix() { GetMatrix(&gm); }
962	void operator=(const BaseMatrix&);
963	void operator=(Real f) { GeneralMatrix::operator=(f); }
964	void operator=(const RowVector& m) { Eq(m); }
965	Real& operator()(int); // access element
966	Real& element(int); // access element
967	Real operator()(int) const; // access element
968	Real element(int) const; // access element
969	#ifdef SETUP_C_SUBSCRIPTS
970	Real& operator[](int m) { return store[m]; }
971	const Real& operator[](int m) const { return store[m]; }
972	// following for Numerical Recipes in C++
973	RowVector(Real a, int n) : Matrix(a, 1, n) {}
974	RowVector(const Real* a, int n) : Matrix(a, 1, n) {}
975	#endif
976	MatrixType type() const;
977	void GetCol(MatrixRowCol&);
978	void GetCol(MatrixColX&);
979	void NextCol(MatrixRowCol&);
980	void NextCol(MatrixColX&);
981	void RestoreCol(MatrixRowCol&) {}
982	void RestoreCol(MatrixColX& c);
983	GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
984	void resize(int); // change dimensions
985	void ReSize(int m) { resize(m); }
986	void resize_keep(int);
987	void resize_keep(int,int);
988	void resize(int,int); // in case access is matrix
989	void ReSize(int m,int n) { resize(m, n); }
990	void resize(const GeneralMatrix& A);
991	void ReSize(const GeneralMatrix& A) { resize(A); }
992	Real* nric() const
993	{ CheckStore(); return store-1; } // for use by NRIC
994	void cleanup(); // to clear store
995	void MiniCleanUp()
996	{ store = 0; storage = 0; nrows_val = 1; ncols_val = 0; tag_val = -1; }
997	// friend ReturnMatrix GetMatrixRow(Matrix& A, int row);
998	void operator+=(const Matrix& M) { PlusEqual(M); }
999	void operator-=(const Matrix& M) { MinusEqual(M); }
1000	void operator+=(Real f) { GeneralMatrix::Add(f); }
1001	void operator-=(Real f) { GeneralMatrix::Add(-f); }
1002	void swap(RowVector& gm)
1003	{ GeneralMatrix::swap((GeneralMatrix&)gm); }
1004	NEW_DELETE(RowVector)
1005	};
1006
1007	/// Column vector.
1008	class ColumnVector : public Matrix
1009	{
1010	GeneralMatrix* Image() const; // copy of matrix
1011	public:
1012	ColumnVector() { ncols_val = 1; }
1013	~ColumnVector() {}
1014	ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
1015	ColumnVector(const BaseMatrix&);
1016	ColumnVector(const ColumnVector& gm) : Matrix() { GetMatrix(&gm); }
1017	void operator=(const BaseMatrix&);
1018	void operator=(Real f) { GeneralMatrix::operator=(f); }
1019	void operator=(const ColumnVector& m) { Eq(m); }
1020	Real& operator()(int); // access element
1021	Real& element(int); // access element
1022	Real operator()(int) const; // access element
1023	Real element(int) const; // access element
1024	#ifdef SETUP_C_SUBSCRIPTS
1025	Real& operator[](int m) { return store[m]; }
1026	const Real& operator[](int m) const { return store[m]; }
1027	// following for Numerical Recipes in C++
1028	ColumnVector(Real a, int m) : Matrix(a, m, 1) {}
1029	ColumnVector(const Real* a, int m) : Matrix(a, m, 1) {}
1030	#endif
1031	MatrixType type() const;
1032	GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1033	void resize(int); // change dimensions
1034	void ReSize(int m) { resize(m); }
1035	void resize_keep(int);
1036	void resize_keep(int,int);
1037	void resize(int,int); // in case access is matrix
1038	void ReSize(int m,int n) { resize(m, n); }
1039	void resize(const GeneralMatrix& A);
1040	void ReSize(const GeneralMatrix& A) { resize(A); }
1041	Real* nric() const
1042	{ CheckStore(); return store-1; } // for use by NRIC
1043	void cleanup(); // to clear store
1044	void MiniCleanUp()
1045	{ store = 0; storage = 0; nrows_val = 0; ncols_val = 1; tag_val = -1; }
1046	// ReturnMatrix Reverse() const; // reverse order of elements
1047	void operator+=(const Matrix& M) { PlusEqual(M); }
1048	void operator-=(const Matrix& M) { MinusEqual(M); }
1049	void operator+=(Real f) { GeneralMatrix::Add(f); }
1050	void operator-=(Real f) { GeneralMatrix::Add(-f); }
1051	void swap(ColumnVector& gm)
1052	{ GeneralMatrix::swap((GeneralMatrix&)gm); }
1053	NEW_DELETE(ColumnVector)
1054	};
1055
1056	/// LU matrix.
1057	/// A square matrix decomposed into upper and lower triangular
1058	/// in preparation for inverting or solving equations.
1059	class CroutMatrix : public GeneralMatrix
1060	{
1061	int* indx;
1062	bool d; // number of row swaps = even or odd
1063	bool sing;
1064	void ludcmp();
1065	void get_aux(CroutMatrix&); // for copying indx[] etc
1066	GeneralMatrix* Image() const; // copy of matrix
1067	public:
1068	CroutMatrix(const BaseMatrix&);
1069	CroutMatrix() : indx(0), d(true), sing(true) {}
1070	CroutMatrix(const CroutMatrix&);
1071	void operator=(const CroutMatrix&);
1072	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1073	MatrixType type() const;
1074	void lubksb(Real*, int=0);
1075	~CroutMatrix();
1076	GeneralMatrix* MakeSolver() { return this; } // for solving
1077	LogAndSign log_determinant() const;
1078	void Solver(MatrixColX&, const MatrixColX&);
1079	void GetRow(MatrixRowCol&);
1080	void GetCol(MatrixRowCol&);
1081	void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
1082	void cleanup(); // to clear store
1083	void MiniCleanUp();
1084	bool IsEqual(const GeneralMatrix&) const;
1085	bool is_singular() const { return sing; }
1086	bool IsSingular() const { return sing; }
1087	const int* const_data_indx() const { return indx; }
1088	bool even_exchanges() const { return d; }
1089	void swap(CroutMatrix& gm);
1090	NEW_DELETE(CroutMatrix)
1091	};
1092
1093	// *************************** band matrices *************************/
1094
1095	/// Band matrix.
1096	class BandMatrix : public GeneralMatrix
1097	{
1098	GeneralMatrix* Image() const; // copy of matrix
1099	protected:
1100	void CornerClear() const; // set unused elements to zero
1101	short SimpleAddOK(const GeneralMatrix* gm);
1102	public:
1103	int lower_val, upper_val; // band widths
1104	BandMatrix() { lower_val=0; upper_val=0; CornerClear(); }
1105	~BandMatrix() {}
1106	BandMatrix(int n,int lb,int ub) { resize(n,lb,ub); CornerClear(); }
1107	// standard declaration
1108	BandMatrix(const BaseMatrix&); // evaluate BaseMatrix
1109	void operator=(const BaseMatrix&);
1110	void operator=(Real f) { GeneralMatrix::operator=(f); }
1111	void operator=(const BandMatrix& m) { Eq(m); }
1112	MatrixType type() const;
1113	Real& operator()(int, int); // access element
1114	Real& element(int, int); // access element
1115	Real operator()(int, int) const; // access element
1116	Real element(int, int) const; // access element
1117	#ifdef SETUP_C_SUBSCRIPTS
1118	Real* operator[](int m) { return store+(upper_val+lower_val)*m+lower_val; }
1119	const Real* operator[](int m) const
1120	{ return store+(upper_val+lower_val)*m+lower_val; }
1121	#endif
1122	BandMatrix(const BandMatrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
1123	LogAndSign log_determinant() const;
1124	GeneralMatrix* MakeSolver();
1125	Real trace() const;
1126	Real sum_square() const
1127	{ CornerClear(); return GeneralMatrix::sum_square(); }
1128	Real sum_absolute_value() const
1129	{ CornerClear(); return GeneralMatrix::sum_absolute_value(); }
1130	Real sum() const
1131	{ CornerClear(); return GeneralMatrix::sum(); }
1132	Real maximum_absolute_value() const
1133	{ CornerClear(); return GeneralMatrix::maximum_absolute_value(); }
1134	Real minimum_absolute_value() const
1135	{ int i, j; return GeneralMatrix::minimum_absolute_value2(i, j); }
1136	Real maximum() const { int i, j; return GeneralMatrix::maximum2(i, j); }
1137	Real minimum() const { int i, j; return GeneralMatrix::minimum2(i, j); }
1138	void GetRow(MatrixRowCol&);
1139	void GetCol(MatrixRowCol&);
1140	void GetCol(MatrixColX&);
1141	void RestoreCol(MatrixRowCol&);
1142	void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
1143	void NextRow(MatrixRowCol&);
1144	virtual void resize(int, int, int); // change dimensions
1145	virtual void ReSize(int m, int n, int b) { resize(m, n, b); }
1146	void resize(const GeneralMatrix& A);
1147	void ReSize(const GeneralMatrix& A) { resize(A); }
1148	//bool SameStorageType(const GeneralMatrix& A) const;
1149	//void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
1150	//void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
1151	MatrixBandWidth bandwidth() const;
1152	void SetParameters(const GeneralMatrix*);
1153	MatrixInput operator<<(double); // will give error
1154	MatrixInput operator<<(float); // will give error
1155	MatrixInput operator<<(int f);
1156	void operator<<(const double* r); // will give error
1157	void operator<<(const float* r); // will give error
1158	void operator<<(const int* r); // will give error
1159	// the next is included because Zortech and Borland
1160	// cannot find the copy in GeneralMatrix
1161	void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
1162	void swap(BandMatrix& gm);
1163	NEW_DELETE(BandMatrix)
1164	};
1165
1166	/// Upper triangular band matrix.
1167	class UpperBandMatrix : public BandMatrix
1168	{
1169	GeneralMatrix* Image() const; // copy of matrix
1170	public:
1171	UpperBandMatrix() {}
1172	~UpperBandMatrix() {}
1173	UpperBandMatrix(int n, int ubw) // standard declaration
1174	: BandMatrix(n, 0, ubw) {}
1175	UpperBandMatrix(const BaseMatrix&); // evaluate BaseMatrix
1176	void operator=(const BaseMatrix&);
1177	void operator=(Real f) { GeneralMatrix::operator=(f); }
1178	void operator=(const UpperBandMatrix& m) { Eq(m); }
1179	MatrixType type() const;
1180	UpperBandMatrix(const UpperBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
1181	GeneralMatrix* MakeSolver() { return this; }
1182	void Solver(MatrixColX&, const MatrixColX&);
1183	LogAndSign log_determinant() const;
1184	void resize(int, int, int); // change dimensions
1185	void ReSize(int m, int n, int b) { resize(m, n, b); }
1186	void resize(int n,int ubw) // change dimensions
1187	{ BandMatrix::resize(n,0,ubw); }
1188	void ReSize(int n,int ubw) // change dimensions
1189	{ BandMatrix::resize(n,0,ubw); }
1190	void resize(const GeneralMatrix& A) { BandMatrix::resize(A); }
1191	void ReSize(const GeneralMatrix& A) { BandMatrix::resize(A); }
1192	Real& operator()(int, int);
1193	Real operator()(int, int) const;
1194	Real& element(int, int);
1195	Real element(int, int) const;
1196	#ifdef SETUP_C_SUBSCRIPTS
1197	Real* operator[](int m) { return store+upper_val*m; }
1198	const Real* operator[](int m) const { return store+upper_val*m; }
1199	#endif
1200	void swap(UpperBandMatrix& gm)
1201	{ BandMatrix::swap((BandMatrix&)gm); }
1202	NEW_DELETE(UpperBandMatrix)
1203	};
1204
1205	/// Lower triangular band matrix.
1206	class LowerBandMatrix : public BandMatrix
1207	{
1208	GeneralMatrix* Image() const; // copy of matrix
1209	public:
1210	LowerBandMatrix() {}
1211	~LowerBandMatrix() {}
1212	LowerBandMatrix(int n, int lbw) // standard declaration
1213	: BandMatrix(n, lbw, 0) {}
1214	LowerBandMatrix(const BaseMatrix&); // evaluate BaseMatrix
1215	void operator=(const BaseMatrix&);
1216	void operator=(Real f) { GeneralMatrix::operator=(f); }
1217	void operator=(const LowerBandMatrix& m) { Eq(m); }
1218	MatrixType type() const;
1219	LowerBandMatrix(const LowerBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
1220	GeneralMatrix* MakeSolver() { return this; }
1221	void Solver(MatrixColX&, const MatrixColX&);
1222	LogAndSign log_determinant() const;
1223	void resize(int, int, int); // change dimensions
1224	void ReSize(int m, int n, int b) { resize(m, n, b); }
1225	void resize(int n,int lbw) // change dimensions
1226	{ BandMatrix::resize(n,lbw,0); }
1227	void ReSize(int n,int lbw) // change dimensions
1228	{ BandMatrix::resize(n,lbw,0); }
1229	void resize(const GeneralMatrix& A) { BandMatrix::resize(A); }
1230	void ReSize(const GeneralMatrix& A) { BandMatrix::resize(A); }
1231	Real& operator()(int, int);
1232	Real operator()(int, int) const;
1233	Real& element(int, int);
1234	Real element(int, int) const;
1235	#ifdef SETUP_C_SUBSCRIPTS
1236	Real* operator[](int m) { return store+lower_val*(m+1); }
1237	const Real* operator[](int m) const { return store+lower_val*(m+1); }
1238	#endif
1239	void swap(LowerBandMatrix& gm)
1240	{ BandMatrix::swap((BandMatrix&)gm); }
1241	NEW_DELETE(LowerBandMatrix)
1242	};
1243
1244	/// Symmetric band matrix.
1245	class SymmetricBandMatrix : public GeneralMatrix
1246	{
1247	GeneralMatrix* Image() const; // copy of matrix
1248	void CornerClear() const; // set unused elements to zero
1249	short SimpleAddOK(const GeneralMatrix* gm);
1250	public:
1251	int lower_val; // lower band width
1252	SymmetricBandMatrix() { lower_val=0; CornerClear(); }
1253	~SymmetricBandMatrix() {}
1254	SymmetricBandMatrix(int n, int lb) { resize(n,lb); CornerClear(); }
1255	SymmetricBandMatrix(const BaseMatrix&);
1256	void operator=(const BaseMatrix&);
1257	void operator=(Real f) { GeneralMatrix::operator=(f); }
1258	void operator=(const SymmetricBandMatrix& m) { Eq(m); }
1259	Real& operator()(int, int); // access element
1260	Real& element(int, int); // access element
1261	Real operator()(int, int) const; // access element
1262	Real element(int, int) const; // access element
1263	#ifdef SETUP_C_SUBSCRIPTS
1264	Real* operator[](int m) { return store+lower_val*(m+1); }
1265	const Real* operator[](int m) const { return store+lower_val*(m+1); }
1266	#endif
1267	MatrixType type() const;
1268	SymmetricBandMatrix(const SymmetricBandMatrix& gm)
1269	: GeneralMatrix() { GetMatrix(&gm); }
1270	GeneralMatrix* MakeSolver();
1271	Real sum_square() const;
1272	Real sum_absolute_value() const;
1273	Real sum() const;
1274	Real maximum_absolute_value() const
1275	{ CornerClear(); return GeneralMatrix::maximum_absolute_value(); }
1276	Real minimum_absolute_value() const
1277	{ int i, j; return GeneralMatrix::minimum_absolute_value2(i, j); }
1278	Real maximum() const { int i, j; return GeneralMatrix::maximum2(i, j); }
1279	Real minimum() const { int i, j; return GeneralMatrix::minimum2(i, j); }
1280	Real trace() const;
1281	LogAndSign log_determinant() const;
1282	void GetRow(MatrixRowCol&);
1283	void GetCol(MatrixRowCol&);
1284	void GetCol(MatrixColX&);
1285	void RestoreCol(MatrixRowCol&) {}
1286	void RestoreCol(MatrixColX&);
1287	GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1288	void resize(int,int); // change dimensions
1289	void ReSize(int m,int b) { resize(m, b); }
1290	void resize(const GeneralMatrix& A);
1291	void ReSize(const GeneralMatrix& A) { resize(A); }
1292	//bool SameStorageType(const GeneralMatrix& A) const;
1293	//void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
1294	//void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
1295	MatrixBandWidth bandwidth() const;
1296	void SetParameters(const GeneralMatrix*);
1297	void operator<<(const double* r); // will give error
1298	void operator<<(const float* r); // will give error
1299	void operator<<(const int* r); // will give error
1300	void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
1301	void swap(SymmetricBandMatrix& gm);
1302	NEW_DELETE(SymmetricBandMatrix)
1303	};
1304
1305	/// LU decomposition of a band matrix.
1306	class BandLUMatrix : public GeneralMatrix
1307	{
1308	int* indx;
1309	bool d;
1310	bool sing; // true if singular
1311	Real* store2;
1312	int storage2;
1313	int m1,m2; // lower and upper
1314	void ludcmp();
1315	void get_aux(BandLUMatrix&); // for copying indx[] etc
1316	GeneralMatrix* Image() const; // copy of matrix
1317	public:
1318	BandLUMatrix(const BaseMatrix&);
1319	BandLUMatrix()
1320	: indx(0), d(true), sing(true), store2(0), storage2(0), m1(0), m2(0) {}
1321	BandLUMatrix(const BandLUMatrix&);
1322	void operator=(const BandLUMatrix&);
1323	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1324	MatrixType type() const;
1325	void lubksb(Real*, int=0);
1326	~BandLUMatrix();
1327	GeneralMatrix* MakeSolver() { return this; } // for solving
1328	LogAndSign log_determinant() const;
1329	void Solver(MatrixColX&, const MatrixColX&);
1330	void GetRow(MatrixRowCol&);
1331	void GetCol(MatrixRowCol&);
1332	void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
1333	void cleanup(); // to clear store
1334	void MiniCleanUp();
1335	bool IsEqual(const GeneralMatrix&) const;
1336	bool is_singular() const { return sing; }
1337	bool IsSingular() const { return sing; }
1338	const Real* const_data2() const { return store2; }
1339	int size2() const { return storage2; }
1340	const int* const_data_indx() const { return indx; }
1341	bool even_exchanges() const { return d; }
1342	MatrixBandWidth bandwidth() const;
1343	void swap(BandLUMatrix& gm);
1344	NEW_DELETE(BandLUMatrix)
1345	};
1346
1347	// ************************ special matrices **************************
1348
1349	/// Identity matrix.
1350	class IdentityMatrix : public GeneralMatrix
1351	{
1352	GeneralMatrix* Image() const; // copy of matrix
1353	public:
1354	IdentityMatrix() {}
1355	~IdentityMatrix() {}
1356	IdentityMatrix(ArrayLengthSpecifier n) : GeneralMatrix(1)
1357	{ nrows_val = ncols_val = n.Value(); *store = 1; }
1358	IdentityMatrix(const IdentityMatrix& gm)
1359	: GeneralMatrix() { GetMatrix(&gm); }
1360	IdentityMatrix(const BaseMatrix&);
1361	void operator=(const BaseMatrix&);
1362	void operator=(const IdentityMatrix& m) { Eq(m); }
1363	void operator=(Real f) { GeneralMatrix::operator=(f); }
1364	MatrixType type() const;
1365
1366	LogAndSign log_determinant() const;
1367	Real trace() const;
1368	Real sum_square() const;
1369	Real sum_absolute_value() const;
1370	Real sum() const { return trace(); }
1371	void GetRow(MatrixRowCol&);
1372	void GetCol(MatrixRowCol&);
1373	void GetCol(MatrixColX&);
1374	void NextRow(MatrixRowCol&);
1375	void NextCol(MatrixRowCol&);
1376	void NextCol(MatrixColX&);
1377	GeneralMatrix* MakeSolver() { return this; } // for solving
1378	void Solver(MatrixColX&, const MatrixColX&);
1379	GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1380	void resize(int n);
1381	void ReSize(int n) { resize(n); }
1382	void resize(const GeneralMatrix& A);
1383	void ReSize(const GeneralMatrix& A) { resize(A); }
1384	MatrixBandWidth bandwidth() const;
1385	// ReturnMatrix Reverse() const; // reverse order of elements
1386	void swap(IdentityMatrix& gm)
1387	{ GeneralMatrix::swap((GeneralMatrix&)gm); }
1388	NEW_DELETE(IdentityMatrix)
1389	};
1390
1391
1392
1393
1394	// ************************ GenericMatrix class **********************/
1395
1396	/// A matrix which can be of any GeneralMatrix type.
1397	class GenericMatrix : public BaseMatrix
1398	{
1399	GeneralMatrix* gm;
1400	int search(const BaseMatrix* bm) const;
1401	friend class BaseMatrix;
1402	public:
1403	GenericMatrix() : gm(0) {}
1404	GenericMatrix(const BaseMatrix& bm)
1405	{ gm = ((BaseMatrix&)bm).Evaluate(); gm = gm->Image(); }
1406	GenericMatrix(const GenericMatrix& bm) : BaseMatrix()
1407	{ gm = bm.gm->Image(); }
1408	void operator=(const GenericMatrix&);
1409	void operator=(const BaseMatrix&);
1410	void operator+=(const BaseMatrix&);
1411	void operator-=(const BaseMatrix&);
1412	void operator*=(const BaseMatrix&);
1413	void operator\|=(const BaseMatrix&);
1414	void operator&=(const BaseMatrix&);
1415	void operator+=(Real);
1416	void operator-=(Real r) { operator+=(-r); }
1417	void operator*=(Real);
1418	void operator/=(Real r) { operator*=(1.0/r); }
1419	~GenericMatrix() { delete gm; }
1420	void cleanup() { delete gm; gm = 0; }
1421	void Release() { gm->Release(); }
1422	void release() { gm->release(); }
1423	GeneralMatrix* Evaluate(MatrixType = MatrixTypeUnSp);
1424	MatrixBandWidth bandwidth() const;
1425	void swap(GenericMatrix& x);
1426	NEW_DELETE(GenericMatrix)
1427	};
1428
1429	// ************************* temporary classes ***********************/
1430
1431	/// Product of two matrices.
1432	/// \internal
1433	class MultipliedMatrix : public BaseMatrix
1434	{
1435	protected:
1436	// if these union statements cause problems, simply remove them
1437	// and declare the items individually
1438	union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
1439	// pointers to summands
1440	union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
1441	MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1442	: bm1(bm1x),bm2(bm2x) {}
1443	int search(const BaseMatrix*) const;
1444	friend class BaseMatrix;
1445	friend class GeneralMatrix;
1446	friend class GenericMatrix;
1447	public:
1448	~MultipliedMatrix() {}
1449	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1450	MatrixBandWidth bandwidth() const;
1451	NEW_DELETE(MultipliedMatrix)
1452	};
1453
1454	/// Sum of two matrices.
1455	/// \internal
1456	class AddedMatrix : public MultipliedMatrix
1457	{
1458	protected:
1459	AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1460	: MultipliedMatrix(bm1x,bm2x) {}
1461
1462	friend class BaseMatrix;
1463	friend class GeneralMatrix;
1464	friend class GenericMatrix;
1465	public:
1466	~AddedMatrix() {}
1467	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1468	MatrixBandWidth bandwidth() const;
1469	NEW_DELETE(AddedMatrix)
1470	};
1471
1472	/// Schur (elementwise) product of two matrices.
1473	/// \internal
1474	class SPMatrix : public AddedMatrix
1475	{
1476	protected:
1477	SPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1478	: AddedMatrix(bm1x,bm2x) {}
1479
1480	friend class BaseMatrix;
1481	friend class GeneralMatrix;
1482	friend class GenericMatrix;
1483	public:
1484	~SPMatrix() {}
1485	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1486	MatrixBandWidth bandwidth() const;
1487
1488	friend SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
1489
1490	NEW_DELETE(SPMatrix)
1491	};
1492
1493	/// Kronecker product of two matrices.
1494	/// \internal
1495	class KPMatrix : public MultipliedMatrix
1496	{
1497	protected:
1498	KPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1499	: MultipliedMatrix(bm1x,bm2x) {}
1500
1501	friend class BaseMatrix;
1502	friend class GeneralMatrix;
1503	friend class GenericMatrix;
1504	public:
1505	~KPMatrix() {}
1506	MatrixBandWidth bandwidth() const;
1507	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1508	friend KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
1509	NEW_DELETE(KPMatrix)
1510	};
1511
1512	/// Two matrices horizontally concatenated.
1513	/// \internal
1514	class ConcatenatedMatrix : public MultipliedMatrix
1515	{
1516	protected:
1517	ConcatenatedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1518	: MultipliedMatrix(bm1x,bm2x) {}
1519
1520	friend class BaseMatrix;
1521	friend class GeneralMatrix;
1522	friend class GenericMatrix;
1523	public:
1524	~ConcatenatedMatrix() {}
1525	MatrixBandWidth bandwidth() const;
1526	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1527	NEW_DELETE(ConcatenatedMatrix)
1528	};
1529
1530	/// Two matrices vertically concatenated.
1531	/// \internal
1532	class StackedMatrix : public ConcatenatedMatrix
1533	{
1534	protected:
1535	StackedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1536	: ConcatenatedMatrix(bm1x,bm2x) {}
1537
1538	friend class BaseMatrix;
1539	friend class GeneralMatrix;
1540	friend class GenericMatrix;
1541	public:
1542	~StackedMatrix() {}
1543	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1544	NEW_DELETE(StackedMatrix)
1545	};
1546
1547	/// Inverted matrix times matrix.
1548	/// \internal
1549	class SolvedMatrix : public MultipliedMatrix
1550	{
1551	SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1552	: MultipliedMatrix(bm1x,bm2x) {}
1553	friend class BaseMatrix;
1554	friend class InvertedMatrix; // for operator*
1555	public:
1556	~SolvedMatrix() {}
1557	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1558	MatrixBandWidth bandwidth() const;
1559	NEW_DELETE(SolvedMatrix)
1560	};
1561
1562	/// Difference between two matrices.
1563	/// \internal
1564	class SubtractedMatrix : public AddedMatrix
1565	{
1566	SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1567	: AddedMatrix(bm1x,bm2x) {}
1568	friend class BaseMatrix;
1569	friend class GeneralMatrix;
1570	friend class GenericMatrix;
1571	public:
1572	~SubtractedMatrix() {}
1573	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1574	NEW_DELETE(SubtractedMatrix)
1575	};
1576
1577	/// Any type of matrix plus Real.
1578	/// \internal
1579	class ShiftedMatrix : public BaseMatrix
1580	{
1581	protected:
1582	union { const BaseMatrix* bm; GeneralMatrix* gm; };
1583	Real f;
1584	ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx),f(fx) {}
1585	int search(const BaseMatrix*) const;
1586	friend class BaseMatrix;
1587	friend class GeneralMatrix;
1588	friend class GenericMatrix;
1589	public:
1590	~ShiftedMatrix() {}
1591	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1592	friend ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
1593	NEW_DELETE(ShiftedMatrix)
1594	};
1595
1596	/// Real minus matrix.
1597	/// \internal
1598	class NegShiftedMatrix : public ShiftedMatrix
1599	{
1600	protected:
1601	NegShiftedMatrix(Real fx, const BaseMatrix* bmx) : ShiftedMatrix(bmx,fx) {}
1602	friend class BaseMatrix;
1603	friend class GeneralMatrix;
1604	friend class GenericMatrix;
1605	public:
1606	~NegShiftedMatrix() {}
1607	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1608	friend NegShiftedMatrix operator-(Real, const BaseMatrix&);
1609	NEW_DELETE(NegShiftedMatrix)
1610	};
1611
1612	/// Any type of matrix times Real.
1613	/// \internal
1614	class ScaledMatrix : public ShiftedMatrix
1615	{
1616	ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
1617	friend class BaseMatrix;
1618	friend class GeneralMatrix;
1619	friend class GenericMatrix;
1620	public:
1621	~ScaledMatrix() {}
1622	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1623	MatrixBandWidth bandwidth() const;
1624	friend ScaledMatrix operator*(Real f, const BaseMatrix& BM);
1625	NEW_DELETE(ScaledMatrix)
1626	};
1627
1628	/// Any type of matrix times -1.
1629	/// \internal
1630	class NegatedMatrix : public BaseMatrix
1631	{
1632	protected:
1633	union { const BaseMatrix* bm; GeneralMatrix* gm; };
1634	NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
1635	int search(const BaseMatrix*) const;
1636	private:
1637	friend class BaseMatrix;
1638	public:
1639	~NegatedMatrix() {}
1640	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1641	MatrixBandWidth bandwidth() const;
1642	NEW_DELETE(NegatedMatrix)
1643	};
1644
1645	/// Transposed matrix.
1646	/// \internal
1647	class TransposedMatrix : public NegatedMatrix
1648	{
1649	TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1650	friend class BaseMatrix;
1651	public:
1652	~TransposedMatrix() {}
1653	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1654	MatrixBandWidth bandwidth() const;
1655	NEW_DELETE(TransposedMatrix)
1656	};
1657
1658	/// Any type of matrix with order of elements reversed.
1659	/// \internal
1660	class ReversedMatrix : public NegatedMatrix
1661	{
1662	ReversedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1663	friend class BaseMatrix;
1664	public:
1665	~ReversedMatrix() {}
1666	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1667	NEW_DELETE(ReversedMatrix)
1668	};
1669
1670	/// Inverse of matrix.
1671	/// \internal
1672	class InvertedMatrix : public NegatedMatrix
1673	{
1674	InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1675	public:
1676	~InvertedMatrix() {}
1677	SolvedMatrix operator(const BaseMatrix&) const; // inverse(A) B
1678	ScaledMatrix operator(Real t) const { return BaseMatrix::operator(t); }
1679	friend class BaseMatrix;
1680	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1681	MatrixBandWidth bandwidth() const;
1682	NEW_DELETE(InvertedMatrix)
1683	};
1684
1685	/// Any type of matrix interpreted as a RowVector.
1686	/// \internal
1687	class RowedMatrix : public NegatedMatrix
1688	{
1689	RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1690	friend class BaseMatrix;
1691	public:
1692	~RowedMatrix() {}
1693	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1694	MatrixBandWidth bandwidth() const;
1695	NEW_DELETE(RowedMatrix)
1696	};
1697
1698	/// Any type of matrix interpreted as a ColumnVector.
1699	/// \internal
1700	class ColedMatrix : public NegatedMatrix
1701	{
1702	ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1703	friend class BaseMatrix;
1704	public:
1705	~ColedMatrix() {}
1706	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1707	MatrixBandWidth bandwidth() const;
1708	NEW_DELETE(ColedMatrix)
1709	};
1710
1711	/// Any type of matrix interpreted as a DiagonalMatrix.
1712	/// \internal
1713	class DiagedMatrix : public NegatedMatrix
1714	{
1715	DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1716	friend class BaseMatrix;
1717	public:
1718	~DiagedMatrix() {}
1719	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1720	MatrixBandWidth bandwidth() const;
1721	NEW_DELETE(DiagedMatrix)
1722	};
1723
1724	/// Any type of matrix interpreted as a (rectangular) Matrix.
1725	/// \internal
1726	class MatedMatrix : public NegatedMatrix
1727	{
1728	int nr, nc;
1729	MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
1730	: NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
1731	friend class BaseMatrix;
1732	public:
1733	~MatedMatrix() {}
1734	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1735	MatrixBandWidth bandwidth() const;
1736	NEW_DELETE(MatedMatrix)
1737	};
1738
1739	/// A matrix in an "envelope' for return from a function.
1740	/// \internal
1741	class ReturnMatrix : public BaseMatrix
1742	{
1743	GeneralMatrix* gm;
1744	int search(const BaseMatrix*) const;
1745	public:
1746	~ReturnMatrix() {}
1747	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1748	friend class BaseMatrix;
1749	ReturnMatrix(const ReturnMatrix& tm) : BaseMatrix(), gm(tm.gm) {}
1750	ReturnMatrix(const GeneralMatrix* gmx) : gm((GeneralMatrix*&)gmx) {}
1751	// ReturnMatrix(GeneralMatrix&);
1752	MatrixBandWidth bandwidth() const;
1753	NEW_DELETE(ReturnMatrix)
1754	};
1755
1756
1757	// ************************ submatrices ****************************/
1758
1759	/// A submatrix of a matrix.
1760	/// \internal
1761	class GetSubMatrix : public NegatedMatrix
1762	{
1763	int row_skip;
1764	int row_number;
1765	int col_skip;
1766	int col_number;
1767	bool IsSym;
1768
1769	GetSubMatrix
1770	(const BaseMatrix* bmx, int rs, int rn, int cs, int cn, bool is)
1771	: NegatedMatrix(bmx),
1772	row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), IsSym(is) {}
1773	void SetUpLHS();
1774	friend class BaseMatrix;
1775	public:
1776	GetSubMatrix(const GetSubMatrix& g)
1777	: NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
1778	col_skip(g.col_skip), col_number(g.col_number), IsSym(g.IsSym) {}
1779	~GetSubMatrix() {}
1780	GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1781	void operator=(const BaseMatrix&);
1782	void operator+=(const BaseMatrix&);
1783	void operator-=(const BaseMatrix&);
1784	void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }
1785	void operator<<(const BaseMatrix&);
1786	void operator<<(const double*); // copy from array
1787	void operator<<(const float*); // copy from array
1788	void operator<<(const int*); // copy from array
1789	MatrixInput operator<<(double); // for loading a list
1790	MatrixInput operator<<(float); // for loading a list
1791	MatrixInput operator<<(int f);
1792	void operator=(Real); // copy from constant
1793	void operator+=(Real); // add constant
1794	void operator-=(Real r) { operator+=(-r); } // subtract constant
1795	void operator*=(Real); // multiply by constant
1796	void operator/=(Real r) { operator*=(1.0/r); } // divide by constant
1797	void inject(const GeneralMatrix&); // copy stored els only
1798	void Inject(const GeneralMatrix& GM) { inject(GM); }
1799	MatrixBandWidth bandwidth() const;
1800	NEW_DELETE(GetSubMatrix)
1801	};
1802
1803	// ****************** linear equation solving **************************/
1804
1805	/// A class for finding A.i() * B.
1806	/// This is supposed to choose the appropriate method depending on the
1807	/// type A. Not very satisfactory as it doesn't know about Cholesky for
1808	/// for positive definite matrices.
1809	class LinearEquationSolver : public BaseMatrix
1810	{
1811	GeneralMatrix* gm;
1812	int search(const BaseMatrix*) const { return 0; }
1813	friend class BaseMatrix;
1814	public:
1815	LinearEquationSolver(const BaseMatrix& bm);
1816	~LinearEquationSolver() { delete gm; }
1817	void cleanup() { delete gm; }
1818	GeneralMatrix* Evaluate(MatrixType) { return gm; }
1819	// probably should have an error message if MatrixType != UnSp
1820	NEW_DELETE(LinearEquationSolver)
1821	};
1822
1823	// ************************ matrix input *****************************/
1824
1825	/// Class for reading values into a (small) matrix within a program.
1826	/// \internal
1827	/// Is able to detect a mismatch in the number of elements.
1828
1829	class MatrixInput
1830	{
1831	int n; // number values still to be read
1832	Real* r; // pointer to next location to be read to
1833	public:
1834	MatrixInput(const MatrixInput& mi) : n(mi.n), r(mi.r) {}
1835	MatrixInput(int nx, Real* rx) : n(nx), r(rx) {}
1836	~MatrixInput();
1837	MatrixInput operator<<(double);
1838	MatrixInput operator<<(float);
1839	MatrixInput operator<<(int f);
1840	friend class GeneralMatrix;
1841	};
1842
1843
1844
1845	// ************** a very simple integer array class ******************/
1846
1847	/// A very simple integer array class.
1848	/// A minimal array class to imitate a C style array but giving dynamic storage
1849	/// mostly intended for internal use by newmat.
1850	/// Probably to be replaced by a templated class when I start using templates.
1851
1852	class SimpleIntArray : public Janitor
1853	{
1854	protected:
1855	int* a; ///< pointer to the array
1856	int n; ///< length of the array
1857	public:
1858	SimpleIntArray(int xn); ///< build an array length xn
1859	SimpleIntArray() : a(0), n(0) {} ///< build an array length 0
1860	~SimpleIntArray(); ///< return the space to memory
1861	int& operator[](int i); ///< access element of the array - start at 0
1862	int operator[](int i) const;
1863	///< access element of constant array
1864	void operator=(int ai); ///< set the array equal to a constant
1865	void operator=(const SimpleIntArray& b);
1866	///< copy the elements of an array
1867	SimpleIntArray(const SimpleIntArray& b);
1868	///< make a new array equal to an existing one
1869	int Size() const { return n; }
1870	///< return the size of the array
1871	int size() const { return n; }
1872	///< return the size of the array
1873	int* Data() { return a; } ///< pointer to the data
1874	const int* Data() const { return a; } ///< pointer to the data
1875	int* data() { return a; } ///< pointer to the data
1876	const int* data() const { return a; } ///< pointer to the data
1877	const int* const_data() const { return a; } ///< pointer to the data
1878	void resize(int i, bool keep = false);
1879	///< change length, keep data if keep = true
1880	void ReSize(int i, bool keep = false) { resize(i, keep); }
1881	///< change length, keep data if keep = true
1882	void resize_keep(int i) { resize(i, true); }
1883	///< change length, keep data
1884	void cleanup() { resize(0); } ///< set length to zero
1885	void CleanUp() { resize(0); } ///< set length to zero
1886	NEW_DELETE(SimpleIntArray)
1887	};
1888
1889	// ******************** C subscript classes **************************
1890
1891	/// Let matrix simulate a C type two dimensional array
1892	class RealStarStar
1893	{
1894	Real** a;
1895	public:
1896	RealStarStar(Matrix& A);
1897	~RealStarStar() { delete [] a; }
1898	operator Real**() { return a; }
1899	};
1900
1901	/// Let matrix simulate a C type const two dimensional array
1902	class ConstRealStarStar
1903	{
1904	const Real** a;
1905	public:
1906	ConstRealStarStar(const Matrix& A);
1907	~ConstRealStarStar() { delete [] a; }
1908	operator const Real**() { return a; }
1909	};
1910
1911	// ************************* exceptions ******************************/
1912
1913	/// Not positive definite exception.
1914	class NPDException : public Runtime_error
1915	{
1916	public:
1917	static unsigned long Select;
1918	NPDException(const GeneralMatrix&);
1919	};
1920
1921	/// Covergence failure exception.
1922	class ConvergenceException : public Runtime_error
1923	{
1924	public:
1925	static unsigned long Select;
1926	ConvergenceException(const GeneralMatrix& A);
1927	ConvergenceException(const char* c);
1928	};
1929
1930	/// Singular matrix exception.
1931	class SingularException : public Runtime_error
1932	{
1933	public:
1934	static unsigned long Select;
1935	SingularException(const GeneralMatrix& A);
1936	};
1937
1938	/// Real overflow exception.
1939	class OverflowException : public Runtime_error
1940	{
1941	public:
1942	static unsigned long Select;
1943	OverflowException(const char* c);
1944	};
1945
1946	/// Miscellaneous exception (details in character string).
1947	class ProgramException : public Logic_error
1948	{
1949	protected:
1950	ProgramException();
1951	public:
1952	static unsigned long Select;
1953	ProgramException(const char* c);
1954	ProgramException(const char* c, const GeneralMatrix&);
1955	ProgramException(const char* c, const GeneralMatrix&, const GeneralMatrix&);
1956	ProgramException(const char* c, MatrixType, MatrixType);
1957	};
1958
1959	/// Index exception.
1960	class IndexException : public Logic_error
1961	{
1962	public:
1963	static unsigned long Select;
1964	IndexException(int i, const GeneralMatrix& A);
1965	IndexException(int i, int j, const GeneralMatrix& A);
1966	// next two are for access via element function
1967	IndexException(int i, const GeneralMatrix& A, bool);
1968	IndexException(int i, int j, const GeneralMatrix& A, bool);
1969	};
1970
1971	/// Cannot convert to vector exception.
1972	class VectorException : public Logic_error
1973	{
1974	public:
1975	static unsigned long Select;
1976	VectorException();
1977	VectorException(const GeneralMatrix& A);
1978	};
1979
1980	/// A matrix is not square exception.
1981	class NotSquareException : public Logic_error
1982	{
1983	public:
1984	static unsigned long Select;
1985	NotSquareException(const GeneralMatrix& A);
1986	NotSquareException();
1987	};
1988
1989	/// Submatrix dimension exception.
1990	class SubMatrixDimensionException : public Logic_error
1991	{
1992	public:
1993	static unsigned long Select;
1994	SubMatrixDimensionException();
1995	};
1996
1997	/// Incompatible dimensions exception.
1998	class IncompatibleDimensionsException : public Logic_error
1999	{
2000	public:
2001	static unsigned long Select;
2002	IncompatibleDimensionsException();
2003	IncompatibleDimensionsException(const GeneralMatrix&);
2004	IncompatibleDimensionsException(const GeneralMatrix&, const GeneralMatrix&);
2005	};
2006
2007	/// Not defined exception.
2008	class NotDefinedException : public Logic_error
2009	{
2010	public:
2011	static unsigned long Select;
2012	NotDefinedException(const char* op, const char* matrix);
2013	};
2014
2015	/// Cannot build matrix with these properties exception.
2016	class CannotBuildException : public Logic_error
2017	{
2018	public:
2019	static unsigned long Select;
2020	CannotBuildException(const char* matrix);
2021	};
2022
2023
2024	/// Internal newmat exception - shouldn't happen.
2025	class InternalException : public Logic_error
2026	{
2027	public:
2028	static unsigned long Select; // for identifying exception
2029	InternalException(const char* c);
2030	};
2031
2032	// ********************** functions ********************************** //
2033
2034	bool operator==(const GeneralMatrix& A, const GeneralMatrix& B);
2035	bool operator==(const BaseMatrix& A, const BaseMatrix& B);
2036	inline bool operator!=(const GeneralMatrix& A, const GeneralMatrix& B)
2037	{ return ! (A==B); }
2038	inline bool operator!=(const BaseMatrix& A, const BaseMatrix& B)
2039	{ return ! (A==B); }
2040
2041	// inequality operators are dummies included for compatibility
2042	// with STL. They throw an exception if actually called.
2043	inline bool operator<=(const BaseMatrix& A, const BaseMatrix&)
2044	{ A.IEQND(); return true; }
2045	inline bool operator>=(const BaseMatrix& A, const BaseMatrix&)
2046	{ A.IEQND(); return true; }
2047	inline bool operator<(const BaseMatrix& A, const BaseMatrix&)
2048	{ A.IEQND(); return true; }
2049	inline bool operator>(const BaseMatrix& A, const BaseMatrix&)
2050	{ A.IEQND(); return true; }
2051
2052	bool is_zero(const BaseMatrix& A);
2053	inline bool IsZero(const BaseMatrix& A) { return is_zero(A); }
2054
2055	Real dotproduct(const Matrix& A, const Matrix& B);
2056	Matrix crossproduct(const Matrix& A, const Matrix& B);
2057	ReturnMatrix crossproduct_rows(const Matrix& A, const Matrix& B);
2058	ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B);
2059
2060	inline Real DotProduct(const Matrix& A, const Matrix& B)
2061	{ return dotproduct(A, B); }
2062	inline Matrix CrossProduct(const Matrix& A, const Matrix& B)
2063	{ return crossproduct(A, B); }
2064	inline ReturnMatrix CrossProductRows(const Matrix& A, const Matrix& B)
2065	{ return crossproduct_rows(A, B); }
2066	inline ReturnMatrix CrossProductColumns(const Matrix& A, const Matrix& B)
2067	{ return crossproduct_columns(A, B); }
2068
2069	void newmat_block_copy(int n, Real* from, Real* to);
2070
2071	// ******************* friend functions ****************************** //
2072
2073	// Functions declared as friends - G++ wants them declared externally as well
2074
2075	bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
2076	bool Compare(const MatrixType&, MatrixType&);
2077	Real dotproduct(const Matrix& A, const Matrix& B);
2078	SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
2079	KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
2080	ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
2081	NegShiftedMatrix operator-(Real, const BaseMatrix&);
2082	ScaledMatrix operator*(Real f, const BaseMatrix& BM);
2083
2084	// ******************* inline functions ****************************** //
2085
2086	inline LogAndSign log_determinant(const BaseMatrix& B)
2087	{ return B.log_determinant(); }
2088	inline LogAndSign LogDeterminant(const BaseMatrix& B)
2089	{ return B.log_determinant(); }
2090	inline Real determinant(const BaseMatrix& B)
2091	{ return B.determinant(); }
2092	inline Real Determinant(const BaseMatrix& B)
2093	{ return B.determinant(); }
2094	inline Real sum_square(const BaseMatrix& B) { return B.sum_square(); }
2095	inline Real SumSquare(const BaseMatrix& B) { return B.sum_square(); }
2096	inline Real norm_Frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
2097	inline Real norm_frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
2098	inline Real NormFrobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
2099	inline Real trace(const BaseMatrix& B) { return B.trace(); }
2100	inline Real Trace(const BaseMatrix& B) { return B.trace(); }
2101	inline Real sum_absolute_value(const BaseMatrix& B)
2102	{ return B.sum_absolute_value(); }
2103	inline Real SumAbsoluteValue(const BaseMatrix& B)
2104	{ return B.sum_absolute_value(); }
2105	inline Real sum(const BaseMatrix& B)
2106	{ return B.sum(); }
2107	inline Real Sum(const BaseMatrix& B)
2108	{ return B.sum(); }
2109	inline Real maximum_absolute_value(const BaseMatrix& B)
2110	{ return B.maximum_absolute_value(); }
2111	inline Real MaximumAbsoluteValue(const BaseMatrix& B)
2112	{ return B.maximum_absolute_value(); }
2113	inline Real minimum_absolute_value(const BaseMatrix& B)
2114	{ return B.minimum_absolute_value(); }
2115	inline Real MinimumAbsoluteValue(const BaseMatrix& B)
2116	{ return B.minimum_absolute_value(); }
2117	inline Real maximum(const BaseMatrix& B) { return B.maximum(); }
2118	inline Real Maximum(const BaseMatrix& B) { return B.maximum(); }
2119	inline Real minimum(const BaseMatrix& B) { return B.minimum(); }
2120	inline Real Minimum(const BaseMatrix& B) { return B.minimum(); }
2121	inline Real norm1(const BaseMatrix& B) { return B.norm1(); }
2122	inline Real Norm1(const BaseMatrix& B) { return B.norm1(); }
2123	inline Real norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
2124	inline Real Norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
2125	inline Real norm_infinity(const BaseMatrix& B) { return B.norm_infinity(); }
2126	inline Real NormInfinity(const BaseMatrix& B) { return B.norm_infinity(); }
2127	inline Real norm_infinity(ColumnVector& CV)
2128	{ return CV.maximum_absolute_value(); }
2129	inline Real NormInfinity(ColumnVector& CV)
2130	{ return CV.maximum_absolute_value(); }
2131	inline bool IsZero(const GeneralMatrix& A) { return A.IsZero(); }
2132	inline bool is_zero(const GeneralMatrix& A) { return A.is_zero(); }
2133
2134
2135	inline MatrixInput MatrixInput::operator<<(int f) { return *this << (Real)f; }
2136	inline MatrixInput GeneralMatrix::operator<<(int f) { return *this << (Real)f; }
2137	inline MatrixInput BandMatrix::operator<<(int f) { return *this << (Real)f; }
2138	inline MatrixInput GetSubMatrix::operator<<(int f) { return *this << (Real)f; }
2139
2140	inline ReversedMatrix BaseMatrix::Reverse() const { return reverse(); }
2141	inline RowedMatrix BaseMatrix::AsRow() const { return as_row(); }
2142	inline ColedMatrix BaseMatrix::AsColumn() const { return as_column(); }
2143	inline DiagedMatrix BaseMatrix::AsDiagonal() const { return as_diagonal(); }
2144	inline MatedMatrix BaseMatrix::AsMatrix(int m, int n) const
2145	{ return as_matrix(m, n); }
2146	inline GetSubMatrix BaseMatrix::SubMatrix(int fr, int lr, int fc, int lc) const
2147	{ return submatrix(fr, lr, fc, lc); }
2148	inline GetSubMatrix BaseMatrix::SymSubMatrix(int f, int l) const
2149	{ return sym_submatrix(f, l); }
2150	inline GetSubMatrix BaseMatrix::Row(int f) const { return row(f); }
2151	inline GetSubMatrix BaseMatrix::Rows(int f, int l) const { return rows(f, l); }
2152	inline GetSubMatrix BaseMatrix::Column(int f) const { return column(f); }
2153	inline GetSubMatrix BaseMatrix::Columns(int f, int l) const
2154	{ return columns(f, l); }
2155	inline Real BaseMatrix::AsScalar() const { return as_scalar(); }
2156
2157	inline ReturnMatrix GeneralMatrix::ForReturn() const { return for_return(); }
2158
2159	inline void swap(Matrix& A, Matrix& B) { A.swap(B); }
2160	inline void swap(SquareMatrix& A, SquareMatrix& B) { A.swap(B); }
2161	inline void swap(nricMatrix& A, nricMatrix& B) { A.swap(B); }
2162	inline void swap(UpperTriangularMatrix& A, UpperTriangularMatrix& B)
2163	{ A.swap(B); }
2164	inline void swap(LowerTriangularMatrix& A, LowerTriangularMatrix& B)
2165	{ A.swap(B); }
2166	inline void swap(SymmetricMatrix& A, SymmetricMatrix& B) { A.swap(B); }
2167	inline void swap(DiagonalMatrix& A, DiagonalMatrix& B) { A.swap(B); }
2168	inline void swap(RowVector& A, RowVector& B) { A.swap(B); }
2169	inline void swap(ColumnVector& A, ColumnVector& B) { A.swap(B); }
2170	inline void swap(CroutMatrix& A, CroutMatrix& B) { A.swap(B); }
2171	inline void swap(BandMatrix& A, BandMatrix& B) { A.swap(B); }
2172	inline void swap(UpperBandMatrix& A, UpperBandMatrix& B) { A.swap(B); }
2173	inline void swap(LowerBandMatrix& A, LowerBandMatrix& B) { A.swap(B); }
2174	inline void swap(SymmetricBandMatrix& A, SymmetricBandMatrix& B) { A.swap(B); }
2175	inline void swap(BandLUMatrix& A, BandLUMatrix& B) { A.swap(B); }
2176	inline void swap(IdentityMatrix& A, IdentityMatrix& B) { A.swap(B); }
2177	inline void swap(GenericMatrix& A, GenericMatrix& B) { A.swap(B); }
2178
2179	#ifdef OPT_COMPATIBLE // for compatibility with opt++
2180
2181	inline Real Norm2(const ColumnVector& CV) { return CV.norm_Frobenius(); }
2182	inline Real Dot(ColumnVector& CV1, ColumnVector& CV2)
2183	{ return dotproduct(CV1, CV2); }
2184
2185	#endif
2186
2187
2188	#ifdef use_namespace
2189	}
2190	#endif
2191
2192
2193	#endif
2194
2195	// body file: newmat1.cpp
2196	// body file: newmat2.cpp
2197	// body file: newmat3.cpp
2198	// body file: newmat4.cpp
2199	// body file: newmat5.cpp
2200	// body file: newmat6.cpp
2201	// body file: newmat7.cpp
2202	// body file: newmat8.cpp
2203	// body file: newmatex.cpp
2204	// body file: bandmat.cpp
2205	// body file: submat.cpp
2206
2207
2208
2209	///@}
2210
2211
2212
2213

Note: See TracBrowser for help on using the repository browser.

Download in other formats: