Context Navigation

source: trunk/linbox/linbox/matrix/matrix-domain.h @ 4086

Revision 4086, 39.8 KB checked in by bboyer, 19 months ago (diff)
clang++ emits no error (almost) and compiles linbox \!
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`

Line
1	/* -- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -- */
2	// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
3	/* linbox/matrix/matrix-domain.h
4	* Copyright (C) 2002 Zhendong Wan, Bradford Hovinen
5	*
6	* Written by Zhendong Wan <wan@mail.eecis.udel.edu>,
7	* Bradford Hovinen <bghovinen@math.uwaterloo.ca>
8	*
9	* ------------------------------------------------------------
10	* 2002-11-26 Bradford Hovinen <bghovinen@math.uwaterloo.ca>
11	*
12	* Added detailed documentation, cleaned up the interface slightly, and added
13	* support for matrix traits. Added read, write, neg, negin, axpy, and
14	* matrix-vector and matrix-black box operations.
15	* ------------------------------------------------------------
16	*
17	* See COPYING for license information.
18	*/
19
20	#ifndef __LINBOX_matrix_domain_H
21	#define __LINBOX_matrix_domain_H
22
23	#include <iostream>
24
25	#include "linbox/blackbox/archetype.h"
26	#include "linbox/vector/vector-domain.h"
27
28	namespace LinBox
29	{
30
31	/** \brief For specializing matrix arithmetic
32	*
33	* This class defines matrix categories that allow us to specialize the matrix
34	* arithmetic in \ref MatrixDomain for different matrix representations. For
35	* example, a sparse matrix may have an efficient iterator over row vectors but
36	* not over column vectors. Therefore, an algorithm that tries to iterate over
37	* column vectors will run very slowly. Hence a specialization that avoids using
38	* column vectors is used instead.
39	*/
40
41	struct MatrixCategories {
42	struct BlackboxTag { };
43	struct RowMatrixTag : public virtual BlackboxTag { };
44	struct ColMatrixTag : public virtual BlackboxTag { };
45	struct RowColMatrixTag : public RowMatrixTag, public ColMatrixTag { };
46	};
47
48	template <class Matrix>
49	struct MatrixTraits {
50	typedef Matrix MatrixType;
51	typedef typename Matrix::MatrixCategory MatrixCategory;
52	};
53
54	/** \brief Helper class to allow specializations of certain matrix-vector products
55	*
56	* This class implements a method mulColSPD that multiplies a
57	* column-represented matrix by a dense vector
58	*/
59	template <class Field>
60	class MVProductDomain {
61	public:
62	typedef typename Field::Element Element;
63
64	MVProductDomain () {}
65
66	protected:
67	template <class Vector1, class Matrix, class Vector2>
68	inline Vector1 &mulColDense (const VectorDomain<Field> &VD, Vector1 &w, const Matrix &A, const Vector2 &v) const;
69	};
70
71	/** Class of matrix arithmetic functions.
72	*
73	* This class encapuslated matrix-matrix and matrix-vector operations, roughly
74	* equivalent to BLAS levels 2 and 3. The arithmetic methods are parameterized
75	* by matrix type so that they may be used the same way with sparse matrices,
76	* dense matrices, and dense submatrices. Except where otherwise noted, they
77	* require the matrix inputs to meet the \ref DenseMatrix archetype.
78	*
79	* These methods are specialized so that they can run efficiently with different
80	* matrix representations. If a matrix has an efficient row iterator, but not an
81	* efficient column iterator, a specialization that makes use of the former will
82	* be selected. This allows a great deal of flexibility when dealing with sparse
83	* matrix arithmetic.
84	*
85	* For all of the arithmetic operations that output matrices, it is assumed that
86	* the output matrix has an efficient row iterator. In typical use, the output
87	* matrix will be a \ref BlasMatrix or a \ref BlasSubmatrix, which has
88	* efficient row and column iterators. In particular, one should not perform
89	* these arithmetic operations outputting to a \ref SparseMatrixBase.
90	*
91	* There are other restrictions. See the method-specific documentation for more
92	* details.
93	*/
94	template <class Field>
95	class MatrixDomain : public MVProductDomain<Field> {
96	public:
97
98	/// Constructor.
99	//! @param F field for MatrixDomain operations.
100	MatrixDomain (const Field &F) :
101	_field (F), _VD (F)
102	{}
103
104	/// Copy operator.
105	MatrixDomain& operator= (const MatrixDomain& MD)
106	{
107	_field = MD._field;
108	_VD = MD._VD;
109	return *this;
110	}
111
112	/** Retrieve the underlying field.
113	* Return a reference to the field that this matrix domain
114	* object uses
115	* @returns reference to field
116	*/
117	//@{
118	const Field &field () const
119	{
120	return _field;
121	}
122	Field &field ()
123	{
124	return _field;
125	}
126	//@}
127
128	/** Print matrix.
129	* @param os Output stream to which matrix is written.
130	* @param A Matrix.
131	* @returns reference to os.
132	*/
133	template <class Matrix>
134	inline std::ostream &write (std::ostream &os, const Matrix &A) const
135	{
136	return A.write (os);
137	}
138
139	/** Read matrix.
140	* @param is Input stream from which matrix is read.
141	* @param A Matrix.
142	* @returns reference to is.
143	*/
144	template <class Matrix>
145	inline std::istream &read (std::istream &is, Matrix &A) const
146	{
147	return A.read (is, _field);
148	}
149
150	/** Matrix copy
151	* B <- A.
152	* Copy the contents of the matrix B to the matrix A
153	*
154	* Both matrices must support the same iterators, row or column.
155	*
156	* @param B Matrix B
157	* @param A Matrix A
158	* @returns Reference to B
159	*/
160	template <class Matrix1, class Matrix2>
161	inline Matrix1 &copy (Matrix1 &B, const Matrix2 &A) const
162	{
163	return copySpecialized (B, A,
164	typename MatrixTraits<Matrix1>::MatrixCategory (),
165	typename MatrixTraits<Matrix2>::MatrixCategory ());
166	}
167
168	/** Matrix equality.
169	* Test whether the matrices A and B are equal
170	* @param A Input vector
171	* @param B Input vector
172	* @returns true if and only if the matrices A and B are equal
173	*/
174	template <class Matrix1, class Matrix2>
175	bool areEqual (const Matrix1 &A, const Matrix2 &B) const
176	{
177	return areEqualSpecialized (B, A,
178	typename MatrixTraits<Matrix1>::MatrixCategory (),
179	typename MatrixTraits<Matrix2>::MatrixCategory ());
180	}
181
182	/** Matrix equality with zero.
183	* @param A Input matrix
184	* @returns true if and only if the matrix A is zero
185	*/
186	template <class Matrix>
187	inline bool isZero (const Matrix &A) const
188	{
189	return isZeroSpecialized (A, typename MatrixTraits<Matrix>::MatrixCategory ());
190	}
191
192	/** Matrix-matrix addition
193	* C <- A + B.
194	*
195	* Each of A, B, and C must support the same iterator, either row or
196	* column
197	*
198	* @param C Output matrix C
199	* @param A Input matrix A
200	* @param B Input matrix B
201	* @returns Reference to C
202	*/
203	template <class Matrix1, class Matrix2, class Matrix3>
204	inline Matrix1& add (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const
205	{
206	return addSpecialized (C, A, B,
207	typename MatrixTraits<Matrix1>::MatrixCategory (),
208	typename MatrixTraits<Matrix2>::MatrixCategory (),
209	typename MatrixTraits<Matrix3>::MatrixCategory ());
210	}
211
212	/** Matrix-matrix in-place addition
213	* A <- A + B.
214	*
215	* Each of A and B must support the same iterator, either row or column
216	*
217	* @param A Input matrix A
218	* @param B Input matrix B
219	* @returns Reference to A
220	*/
221	template <class Matrix1, class Matrix2>
222	inline Matrix1& addin (Matrix1 &A, const Matrix2 &B) const
223	{
224	return addinSpecialized (A, B,
225	typename MatrixTraits<Matrix1>::MatrixCategory (),
226	typename MatrixTraits<Matrix2>::MatrixCategory ());
227	}
228
229	/** Matrix-matrix subtraction
230	* C <- A - B.
231	*
232	* Each of A, B, and C must support the same iterator, either row or
233	* column
234	*
235	* @param C Output matrix C
236	* @param A Input matrix A
237	* @param B Input matrix B
238	* @returns Reference to C
239	*/
240	template <class Matrix1, class Matrix2, class Matrix3>
241	inline Matrix1 &sub (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const
242	{
243	return subSpecialized (C, A, B,
244	typename MatrixTraits<Matrix1>::MatrixCategory (),
245	typename MatrixTraits<Matrix2>::MatrixCategory (),
246	typename MatrixTraits<Matrix3>::MatrixCategory ());
247	}
248
249	/** Matrix-matrix in-place subtraction
250	* A <- A - B.
251	*
252	* Each of A and B must support the same iterator, either row or column
253	*
254	* @param A Input matrix A
255	* @param B Input matrix B
256	* @returns Reference to A
257	*/
258	template <class Matrix1, class Matrix2>
259	inline Matrix1 &subin (Matrix1 &A, const Matrix2 &B) const
260	{
261	return subinSpecialized (A, B,
262	typename MatrixTraits<Matrix1>::MatrixCategory (),
263	typename MatrixTraits<Matrix2>::MatrixCategory ());
264	}
265
266	/** Matrix negate
267	* B <- -A.
268	*
269	* Each of A and B must support the same iterator, either row or column
270	*
271	* @param B Output matrix B
272	* @param A Input matrix A
273	* @returns reference to B
274	*/
275	template <class Matrix1, class Matrix2>
276	inline Matrix1 &neg (Matrix1 &B, const Matrix2 &A) const
277	{
278	return negSpecialized (B, A,
279	typename MatrixTraits<Matrix1>::MatrixCategory (),
280	typename MatrixTraits<Matrix2>::MatrixCategory ());
281	}
282
283	/** Matrix in-place negate
284	* A <- -A.
285	* @param A Input matrix A; result is stored here
286	*/
287	template <class Matrix>
288	inline Matrix &negin (Matrix &A) const
289	{
290	return neginSpecialized (A, typename MatrixTraits<Matrix>::MatrixCategory ());
291	}
292
293	/** Matrix-matrix multiply
294	* C <- A * B.
295	*
296	* C must support both row and column iterators, and the vector
297	* representations must be dense. Examples of supported matrices are
298	* \ref BlasMatrix and \ref BlasSubmatrix.
299	*
300	* Either A or B, or both, may have limited iterators. However, either A
301	* must support row iterators or B must support column iterators. If
302	* both A and B lack support for an iterator (either row or column),
303	* then C must support the same type of iterator as A and B.
304	*
305	* @param C Output matrix C
306	* @param A Input matrix A
307	* @param B Input matrix B
308	* @returns Reference to C
309	*/
310	template <class Matrix1, class Matrix2, class Matrix3>
311	inline Matrix1 &mul (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const
312	{
313	return mulSpecialized (C, A, B,
314	typename MatrixTraits<Matrix1>::MatrixCategory (),
315	typename MatrixTraits<Matrix2>::MatrixCategory (),
316	typename MatrixTraits<Matrix3>::MatrixCategory ());
317	}
318
319	/** Matrix-matrix in-place multiply on the left
320	* B <- A * B.
321	*
322	* B should support both row and column iterators, and must be dense. A
323	* must support row iterators.
324	*
325	* @param A Input matrix A
326	* @param B Input matrix B
327	* @returns Reference to B
328	*/
329	template <class Matrix1, class Matrix2>
330	inline Matrix2 &leftMulin (const Matrix1 &A, Matrix2 &B) const;
331
332	/** Matrix-matrix in-place multiply on the right
333	* A <- A * B.
334	*
335	* A should support both row and column iterators, and must be dense. B
336	* must support column iterators.
337	*
338	* @param A Input matrix A
339	* @param B Input matrix B
340	* @returns Reference to A
341	*/
342	template <class Matrix1, class Matrix2>
343	inline Matrix1 &rightMulin (Matrix1 &A, const Matrix2 &B) const;
344
345	/** Matrix-matrix in-place multiply
346	* A <- A * B.
347	*
348	* This is an alias for \ref rightMulin
349	*
350	* @param A Input matrix A
351	* @param B Input matrix B
352	* @returns Reference to A
353	*/
354	template <class Matrix1, class Matrix2>
355	inline Matrix1 &mulin (Matrix1 &A, const Matrix2 &B) const
356	{
357	return rightMulin (A, B);
358	}
359
360	/** Matrix-scalar multiply
361	* C <- B * a.
362	*
363	* Multiply B by the scalar element a and store the result in C. B and C
364	* must support the same iterators.
365	*
366	* @param C Output matrix C
367	* @param B Input matrix B
368	* @param a Input scalar a
369	* @returns Reference to C
370	*/
371	template <class Matrix1, class Matrix2>
372	inline Matrix1 &mul (Matrix1 &C, const Matrix2 &B, const typename Field::Element &a) const
373	{
374	return mulSpecialized (C, B, a,
375	typename MatrixTraits<Matrix1>::MatrixCategory (),
376	typename MatrixTraits<Matrix2>::MatrixCategory ());
377	}
378
379	/** Matrix-scalar in-place multiply
380	* B <- B * a.
381	*
382	* Multiply B by the scalar element a in-place.
383	*
384	* @param B Input matrix B
385	* @param a Input scalar a
386	* @returns Reference to B
387	*/
388	template <class Matrix>
389	inline Matrix &mulin (Matrix &B, const typename Field::Element &a) const
390	{
391	return mulinSpecialized (B, a, typename MatrixTraits<Matrix>::MatrixCategory ());
392	}
393
394	/** Matrix-matrix in-place axpy
395	* Y <- Y + A*X.
396	*
397	* This function combines \ref mul and \ref add, eliminating the need
398	* for an additional temporary in expressions of the form $Y = Y +
399	* AX$. Only one row of additional storage is required. Y may have
400	* either efficient row iterators or efficient column iterators, and the
401	* same restrictions on A and X apply as in \ref mul.
402	*
403	* Note that no out-of-place axpy is provided, since it gives no
404	* benefit. One may just as easily multiply into the result and call
405	* \ref addin.
406	*
407	* @param Y Input matrix Y; result is stored here
408	* @param A Input matrix A
409	* @param X Input matrix X
410	*/
411	template <class Matrix1, class Matrix2, class Matrix3>
412	inline Matrix1 &axpyin (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X) const
413	{
414	return axpyinSpecialized (Y, A, X,
415	typename MatrixTraits<Matrix1>::MatrixCategory (),
416	typename MatrixTraits<Matrix2>::MatrixCategory (),
417	typename MatrixTraits<Matrix3>::MatrixCategory ());
418	}
419
420	//! Y <- AX-Y
421	template <class Matrix1, class Matrix2, class Matrix3>
422	inline Matrix1 &axmyin (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X) const
423	{
424	negin(Y);
425	axpyin(Y,A,X);
426	return Y;
427	}
428
429
430	/*! General matrix multiply
431	* \f$ D \gets \alpha A B + \beta C\f$.
432	* @todo not efficient...
433	*/
434	template <class Matrix1, class Matrix2, class Matrix3>
435	inline Matrix1 &muladd (Matrix1 & D,
436	const typename Field::Element & beta,
437	const Matrix1 & C,
438	const typename Field::Element & alpha,
439	const Matrix2 & A,
440	const Matrix3 & B) const
441	{
442	mul(D,A,B); // D = AB
443	mulin(D,alpha); // D = alpha D
444	Matrix1 CC(C);
445	mulin(CC,beta); // C = beta C
446	addin(D,CC); // D = D+C
447	return D;
448	}
449
450	/! @todo Need documentation of these methods /
451	//@{
452	template<class Matrix1, class Matrix2>
453	Matrix1 &pow_apply (Matrix1 &M1, const Matrix2 &M2, unsigned long int k) const;
454
455	template<class Matrix1, class Matrix2>
456	Matrix1 &pow_horn (Matrix1 &M1, const Matrix2 &M2, unsigned long int k) const;
457	//@}
458
459
460	/*! @name Matrix-vector arithmetic operations
461	* These operations take a matrix satisfying the \ref DenseMatrix
462	* archetype and LinBox vectors as inputs. They involve matrix-vector
463	* product and matrix-vector AXPY
464	*/
465	//@{
466	/** Matrix-vector multiply
467	* w <- A * v.
468	*
469	* The vectors v and w must be of the same representation (dense, sparse
470	* sequence, sparse associative, or sparse parallel), but they may be of
471	* different types. The matrix A may have any representation.
472	*
473	* @param w Output vector w
474	* @param A Input matrix A
475	* @param v Input vector v
476	* @returns Reference to w
477	*/
478	template <class Vector1, class Matrix, class Vector2>
479	inline Vector1 &vectorMul (Vector1 &w, const Matrix &A, const Vector2 &v) const
480	{
481	return mulSpecialized (w, A, v, typename MatrixTraits<Matrix>::MatrixCategory ());
482	}
483
484	/** Matrix-vector in-place axpy
485	* \f$y \gets y + A x\f$.
486	*
487	* This function eliminates the requirement for temporary storage when
488	* one is computing an expression of the form given above.
489	*
490	* The vectors y and x must be of the same representation (dense, sparse
491	* sequence, sparse associative, or sparse parallel), but they may be of
492	* different types. The matrix A may have any representation.
493	*
494	* Note that out-of-place axpy is not provided since it provides no
495	* benefit -- one can use mul and then addin to exactly the same effect,
496	* with no additional storage or performance cost.
497	*
498	* @param y Input vector y; result is stored here
499	* @param A Input matrix A
500	* @param x Input vector x
501	*/
502	template <class Vector1, class Matrix, class Vector2>
503	inline Vector1 &vectorAxpyin (Vector1 &y, const Matrix &A, const Vector2 &x) const
504	{
505	return axpyinSpecialized (y, A, x, typename MatrixTraits<Matrix>::MatrixCategory ());
506	}
507	//@}
508
509	/*! @name Matrix-black box arithmetic operations
510	* These operations mimic the matrix-matrix arithmetic operations above,
511	* but one of the parameters is a \ref BlackboxArchetype.
512	*/
513	//@{
514	/** Matrix-black box left-multiply
515	* C <- A * B.
516	*
517	* Both C and B must support column iterators
518	*
519	* @param C Output matrix
520	* @param A Black box for A
521	* @param B Matrix B
522	*/
523	template <class Matrix1, class Blackbox, class Matrix2>
524	inline Matrix1 &blackboxMulLeft (Matrix1 &C, const Blackbox &A, const Matrix2 &B) const;
525
526	/** Matrix-black box right-multiply
527	* C <- A * B.
528	*
529	* Both C and A must support row iterators
530	*
531	* @param C Output matrix
532	* @param A Matrix A
533	* @param B Black box for B
534	*/
535	template <class Matrix1, class Matrix2, class Blackbox>
536	inline Matrix1 &blackboxMulRight (Matrix1 &C, const Matrix2 &A, const Blackbox &B) const;
537	//@}
538
539	/*! @name Matrix permutations
540	* @brief
541	* These operations permute the rows or columns of a matrix based on
542	* the given permutation. They are intended for use with Gauss-Jordan
543	* elimination
544	*/
545	//@{
546	/// Transposition.
547	typedef std::pair<unsigned int, unsigned int> Transposition;
548	/** Permutation.
549	*
550	* A permutation is represented as a vector of pairs, each
551	* pair representing a transposition.
552	*/
553	typedef std::vector<Transposition> Permutation;
554
555
556	/** Permute the rows of the given matrix.
557	*
558	* @param A Output matrix
559	* @param P_start Start of permutation
560	* @param P_end End of permutation
561	* @returns Reference to A
562	*/
563	template <class Matrix, class Iterator>
564	inline Matrix &permuteRows (Matrix &A,
565	Iterator P_start,
566	Iterator P_end) const
567	{
568	return permuteRowsSpecialized (A, P_start, P_end,
569	typename MatrixTraits<Matrix>::MatrixCategory ());
570	}
571
572	/** Permute the columns of the given matrix.
573	*
574	* @param A Output matrix
575	* @param P_start Start of permutation
576	* @param P_end End of permutation
577	* @returns Reference to A
578	*/
579	template <class Matrix, class Iterator>
580	inline Matrix &permuteColumns (Matrix &A,
581	Iterator P_start,
582	Iterator P_end) const
583	{
584	return permuteColsSpecialized (A, P_start, P_end,
585	typename MatrixTraits<Matrix>::MatrixCategory ());
586	}
587	//@}
588
589	private:
590
591	// Specialized function implementations
592	template <class Matrix1, class Matrix2> Matrix1 &copyRow (Matrix1 &B, const Matrix2 &A) const;
593	template <class Matrix1, class Matrix2> Matrix1 &copyCol (Matrix1 &B, const Matrix2 &A) const;
594
595	template <class Matrix1, class Matrix2>
596	inline Matrix1 &copySpecialized (Matrix1 &B, const Matrix2 &A,
597	MatrixCategories::RowMatrixTag,
598	MatrixCategories::RowMatrixTag) const
599	{
600	return copyRow (B, A);
601	}
602	template <class Matrix1, class Matrix2>
603	inline Matrix1 &copySpecialized (Matrix1 &B, const Matrix2 &A,
604	MatrixCategories::ColMatrixTag,
605	MatrixCategories::ColMatrixTag) const
606	{
607	return copyCol (B, A);
608	}
609	template <class Matrix1, class Matrix2>
610	inline Matrix1 &copySpecialized (Matrix1 &B, const Matrix2 &A,
611	MatrixCategories::RowColMatrixTag,
612	MatrixCategories::RowColMatrixTag) const
613	{
614	return copyRow (B, A);
615	}
616
617	template <class Matrix1, class Matrix2> bool areEqualRow (const Matrix1 &A, const Matrix2 &B) const;
618	template <class Matrix1, class Matrix2> bool areEqualCol (const Matrix1 &A, const Matrix2 &B) const;
619
620	template <class Matrix1, class Matrix2>
621	inline bool areEqualSpecialized (const Matrix1 &A, const Matrix2 &B,
622	MatrixCategories::RowMatrixTag,
623	MatrixCategories::RowMatrixTag) const
624	{
625	return areEqualRow (A, B);
626	}
627	template <class Matrix1, class Matrix2>
628	inline bool areEqualSpecialized (const Matrix1 &A, const Matrix2 &B,
629	MatrixCategories::ColMatrixTag,
630	MatrixCategories::ColMatrixTag) const
631	{
632	return areEqualCol (A, B);
633	}
634	template <class Matrix1, class Matrix2>
635	inline bool areEqualSpecialized (const Matrix1 &A, const Matrix2 &B,
636	MatrixCategories::RowColMatrixTag,
637	MatrixCategories::RowColMatrixTag) const
638	{
639	return areEqualRow (A, B);
640	}
641
642	template <class Matrix> bool isZeroRow (const Matrix &v) const;
643	template <class Matrix> bool isZeroCol (const Matrix &v) const;
644
645	template <class Matrix>
646	bool isZeroSpecialized (const Matrix &A, MatrixCategories::RowMatrixTag) const
647	{
648	return isZeroRow (A);
649	}
650	template <class Matrix>
651	bool isZeroSpecialized (const Matrix &A, MatrixCategories::ColMatrixTag) const
652	{
653	return isZeroCol (A);
654	}
655	template <class Matrix>
656	bool isZeroSpecialized (const Matrix &A, MatrixCategories::RowColMatrixTag) const
657	{
658	return isZeroRow (A);
659	}
660
661	template <class Matrix1, class Matrix2, class Matrix3>
662	Matrix1& addRow (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const;
663	template <class Matrix1, class Matrix2, class Matrix3>
664	Matrix1& addCol (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const;
665
666	template <class Matrix1, class Matrix2, class Matrix3>
667	Matrix1& addSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
668	MatrixCategories::RowMatrixTag,
669	MatrixCategories::RowMatrixTag,
670	MatrixCategories::RowMatrixTag) const
671	{
672	return addRow (C, A, B);
673	}
674	template <class Matrix1, class Matrix2, class Matrix3>
675	Matrix1& addSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
676	MatrixCategories::ColMatrixTag,
677	MatrixCategories::ColMatrixTag,
678	MatrixCategories::ColMatrixTag) const
679	{
680	return addCol (C, A, B);
681	}
682	template <class Matrix1, class Matrix2, class Matrix3>
683	Matrix1& addSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
684	MatrixCategories::RowColMatrixTag,
685	MatrixCategories::RowColMatrixTag,
686	MatrixCategories::RowColMatrixTag) const
687	{
688	return addRow (C, A, B);
689	}
690
691	template <class Matrix1, class Matrix2> Matrix1& addinRow (Matrix1 &A, const Matrix2 &B) const;
692	template <class Matrix1, class Matrix2> Matrix1& addinCol (Matrix1 &A, const Matrix2 &B) const;
693
694	template <class Matrix1, class Matrix2>
695	inline Matrix1& addinSpecialized (Matrix1 &A, const Matrix2 &B,
696	MatrixCategories::RowMatrixTag,
697	MatrixCategories::RowMatrixTag) const
698	{
699	return addinRow (A, B);
700	}
701	template <class Matrix1, class Matrix2>
702	inline Matrix1& addinSpecialized (Matrix1 &A, const Matrix2 &B,
703	MatrixCategories::ColMatrixTag,
704	MatrixCategories::ColMatrixTag) const
705	{
706	return addinCol (A, B);
707	}
708	template <class Matrix1, class Matrix2>
709	inline Matrix1& addinSpecialized (Matrix1 &A, const Matrix2 &B,
710	MatrixCategories::RowColMatrixTag,
711	MatrixCategories::RowColMatrixTag) const
712	{
713	return addinRow (A, B);
714	}
715
716	template <class Matrix1, class Matrix2, class Matrix3>
717	Matrix1& subRow (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const;
718	template <class Matrix1, class Matrix2, class Matrix3>
719	Matrix1& subCol (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const;
720
721	template <class Matrix1, class Matrix2, class Matrix3>
722	Matrix1& subSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
723	MatrixCategories::RowMatrixTag,
724	MatrixCategories::RowMatrixTag,
725	MatrixCategories::RowMatrixTag) const
726	{
727	return subRow (C, A, B);
728	}
729	template <class Matrix1, class Matrix2, class Matrix3>
730	Matrix1& subSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
731	MatrixCategories::ColMatrixTag,
732	MatrixCategories::ColMatrixTag,
733	MatrixCategories::ColMatrixTag) const
734	{
735	return subCol (C, A, B);
736	}
737	template <class Matrix1, class Matrix2, class Matrix3>
738	Matrix1& subSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
739	MatrixCategories::RowColMatrixTag,
740	MatrixCategories::RowColMatrixTag,
741	MatrixCategories::RowColMatrixTag) const
742	{
743	return subRow (C, A, B);
744	}
745
746	template <class Matrix1, class Matrix2> Matrix1& subinRow (Matrix1 &A, const Matrix2 &B) const;
747	template <class Matrix1, class Matrix2> Matrix1& subinCol (Matrix1 &A, const Matrix2 &B) const;
748
749	template <class Matrix1, class Matrix2>
750	Matrix1& subinSpecialized (Matrix1 &A, const Matrix2 &B,
751	MatrixCategories::RowMatrixTag,
752	MatrixCategories::RowMatrixTag) const
753	{
754	return subinRow (A, B);
755	}
756	template <class Matrix1, class Matrix2>
757	Matrix1& subinSpecialized (Matrix1 &A, const Matrix2 &B,
758	MatrixCategories::ColMatrixTag,
759	MatrixCategories::ColMatrixTag) const
760	{
761	return subinCol (A, B);
762	}
763	template <class Matrix1, class Matrix2>
764	Matrix1& subinSpecialized (Matrix1 &A, const Matrix2 &B,
765	MatrixCategories::RowColMatrixTag,
766	MatrixCategories::RowColMatrixTag) const
767	{
768	return subinRow (A, B);
769	}
770
771	template <class Matrix1, class Matrix2> Matrix1& negRow (Matrix1 &A, const Matrix2 &B) const;
772	template <class Matrix1, class Matrix2> Matrix1& negCol (Matrix1 &A, const Matrix2 &B) const;
773
774	template <class Matrix1, class Matrix2>
775	inline Matrix1& negSpecialized (Matrix1 &A, const Matrix2 &B,
776	MatrixCategories::RowMatrixTag,
777	MatrixCategories::RowMatrixTag) const
778	{
779	return negRow (A, B);
780	}
781	template <class Matrix1, class Matrix2>
782	inline Matrix1& negSpecialized (Matrix1 &A, const Matrix2 &B,
783	MatrixCategories::ColMatrixTag,
784	MatrixCategories::ColMatrixTag) const
785	{
786	return negCol (A, B);
787	}
788	template <class Matrix1, class Matrix2>
789	inline Matrix1& negSpecialized (Matrix1 &A, const Matrix2 &B,
790	MatrixCategories::RowColMatrixTag,
791	MatrixCategories::RowColMatrixTag) const
792	{
793	return negRow (A, B);
794	}
795
796	template <class Matrix> Matrix &neginRow (Matrix &A) const;
797	template <class Matrix> Matrix &neginCol (Matrix &A) const;
798
799	template <class Matrix>
800	Matrix &neginSpecialized (Matrix &A, MatrixCategories::RowMatrixTag) const
801	{
802	return neginRow (A);
803	}
804	template <class Matrix>
805	Matrix &neginSpecialized (Matrix &A, MatrixCategories::ColMatrixTag) const
806	{
807	return neginCol (A);
808	}
809	template <class Matrix>
810	Matrix &neginSpecialized (Matrix &A, MatrixCategories::RowColMatrixTag) const
811	{
812	return neginRow (A);
813	}
814
815	template <class Matrix1, class Matrix2, class Matrix3>
816	Matrix1 &mulRowRowCol (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const;
817	template <class Matrix1, class Matrix2, class Matrix3>
818	Matrix1 &mulColRowCol (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const;
819	template <class Matrix1, class Matrix2, class Matrix3>
820	Matrix1 &mulRowRowRow (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const;
821	template <class Matrix1, class Matrix2, class Matrix3>
822	Matrix1 &mulColColCol (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const;
823
824	template <class Matrix1, class Matrix2, class Matrix3>
825	Matrix1 &mulSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
826	MatrixCategories::RowMatrixTag,
827	MatrixCategories::RowMatrixTag,
828	MatrixCategories::ColMatrixTag) const
829	{
830	return mulRowRowCol (C, A, B);
831	}
832	template <class Matrix1, class Matrix2, class Matrix3>
833	Matrix1 &mulSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
834	MatrixCategories::ColMatrixTag,
835	MatrixCategories::RowMatrixTag,
836	MatrixCategories::ColMatrixTag) const
837	{
838	return mulColRowCol (C, A, B);
839	}
840	template <class Matrix1, class Matrix2, class Matrix3>
841	Matrix1 &mulSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
842	MatrixCategories::RowColMatrixTag,
843	MatrixCategories::RowMatrixTag,
844	MatrixCategories::ColMatrixTag) const
845	{
846	return mulRowRowCol (C, A, B);
847	}
848	template <class Matrix1, class Matrix2, class Matrix3>
849	Matrix1 &mulSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
850	MatrixCategories::RowMatrixTag,
851	MatrixCategories::RowMatrixTag,
852	MatrixCategories::RowMatrixTag) const
853	{
854	return mulRowRowRow (C, A, B);
855	}
856	template <class Matrix1, class Matrix2, class Matrix3>
857	Matrix1 &mulSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
858	MatrixCategories::ColMatrixTag,
859	MatrixCategories::ColMatrixTag,
860	MatrixCategories::ColMatrixTag) const
861	{
862	return mulColColCol (C, A, B);
863	}
864	template <class Matrix1, class Matrix2, class Matrix3>
865	Matrix1 &mulSpecialized (Matrix1 &C, const Matrix2 &A, const Matrix3 &B,
866	MatrixCategories::RowColMatrixTag,
867	MatrixCategories::RowColMatrixTag,
868	MatrixCategories::RowColMatrixTag) const
869	{
870	return mulRowRowCol (C, A, B);
871	}
872
873	template <class Matrix1, class Matrix2>
874	Matrix1 &mulRow (Matrix1 &C, const Matrix2 &B, const typename Field::Element &a) const;
875	template <class Matrix1, class Matrix2>
876	Matrix1 &mulCol (Matrix1 &C, const Matrix2 &B, const typename Field::Element &a) const;
877
878	template <class Matrix1, class Matrix2>
879	Matrix1 &mulSpecialized (Matrix1 &C, const Matrix2 &B, const typename Field::Element &a,
880	MatrixCategories::RowMatrixTag,
881	MatrixCategories::RowMatrixTag) const
882	{
883	return mulRow (C, B, a);
884	}
885	template <class Matrix1, class Matrix2>
886	Matrix1 &mulSpecialized (Matrix1 &C, const Matrix2 &B, const typename Field::Element &a,
887	MatrixCategories::ColMatrixTag,
888	MatrixCategories::ColMatrixTag) const
889	{
890	return mulCol (C, B, a);
891	}
892	template <class Matrix1, class Matrix2>
893	Matrix1 &mulSpecialized (Matrix1 &C, const Matrix2 &B, const typename Field::Element &a,
894	MatrixCategories::RowColMatrixTag,
895	MatrixCategories::RowColMatrixTag) const
896	{
897	return mulRow (C, B, a);
898	}
899
900	template <class Matrix> Matrix &mulinRow (Matrix &B, const typename Field::Element &a) const;
901	template <class Matrix> Matrix &mulinCol (Matrix &B, const typename Field::Element &a) const;
902
903	template <class Matrix>
904	Matrix &mulinSpecialized (Matrix &B, const typename Field::Element &a,
905	MatrixCategories::RowMatrixTag) const
906	{
907	return mulinRow (B, a);
908	}
909	template <class Matrix>
910	Matrix &mulinSpecialized (Matrix &B, const typename Field::Element &a,
911	MatrixCategories::ColMatrixTag) const
912	{
913	return mulinCol (B, a);
914	}
915	template <class Matrix>
916	Matrix &mulinSpecialized (Matrix &B, const typename Field::Element &a,
917	MatrixCategories::RowColMatrixTag) const
918	{
919	return mulinRow (B, a);
920	}
921
922	template <class Matrix1, class Matrix2, class Matrix3>
923	Matrix1 &axpyinRowRowCol (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X) const;
924	template <class Matrix1, class Matrix2, class Matrix3>
925	Matrix1 &axpyinColRowCol (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X) const;
926	template <class Matrix1, class Matrix2, class Matrix3>
927	Matrix1 &axpyinRowRowRow (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X) const;
928	template <class Matrix1, class Matrix2, class Matrix3>
929	Matrix1 &axpyinColColCol (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X) const;
930
931	template <class Matrix1, class Matrix2, class Matrix3>
932	Matrix1 &axpyinSpecialized (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X,
933	MatrixCategories::RowMatrixTag,
934	MatrixCategories::RowMatrixTag,
935	MatrixCategories::ColMatrixTag) const
936	{
937	return axpyinRowRowCol (Y, A, X);
938	}
939	template <class Matrix1, class Matrix2, class Matrix3>
940	Matrix1 &axpyinSpecialized (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X,
941	MatrixCategories::ColMatrixTag,
942	MatrixCategories::RowMatrixTag,
943	MatrixCategories::ColMatrixTag) const
944	{
945	return axpyinColRowCol (Y, A, X);
946	}
947	template <class Matrix1, class Matrix2, class Matrix3>
948	Matrix1 &axpyinSpecialized (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X,
949	MatrixCategories::RowColMatrixTag,
950	MatrixCategories::RowMatrixTag,
951	MatrixCategories::ColMatrixTag) const
952	{
953	return axpyinRowRowCol (Y, A, X);
954	}
955	template <class Matrix1, class Matrix2, class Matrix3>
956	Matrix1 &axpyinSpecialized (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X,
957	MatrixCategories::RowMatrixTag,
958	MatrixCategories::RowMatrixTag,
959	MatrixCategories::RowMatrixTag) const
960	{
961	return axpyinRowRowRow (Y, A, X);
962	}
963
964	template <class Matrix1, class Matrix2, class Matrix3>
965	Matrix1 &axpyinSpecialized (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X,
966	MatrixCategories::ColMatrixTag,
967	MatrixCategories::ColMatrixTag,
968	MatrixCategories::ColMatrixTag) const
969	{
970	return axpyinColColCol (Y, A, X);
971	}
972
973	template <class Matrix1, class Matrix2, class Matrix3>
974	Matrix1 &axpyinSpecialized (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X,
975	MatrixCategories::RowColMatrixTag,
976	MatrixCategories::RowColMatrixTag,
977	MatrixCategories::RowColMatrixTag) const
978	{
979	return axpyinRowRowCol (Y, A, X);
980	}
981
982	template <class Vector1, class Matrix, class Vector2>
983	Vector1 &mulRowSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
984	VectorCategories::DenseVectorTag) const;
985	template <class Vector1, class Matrix, class Vector2>
986	Vector1 &mulRowSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
987	VectorCategories::SparseSequenceVectorTag) const;
988	template <class Vector1, class Matrix, class Vector2>
989	Vector1 &mulRowSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
990	VectorCategories::SparseAssociativeVectorTag) const;
991	template <class Vector1, class Matrix, class Vector2>
992	Vector1 &mulRowSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
993	VectorCategories::SparseParallelVectorTag) const;
994
995	template <class Vector1, class Matrix, class Vector2>
996	Vector1 &mulColSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
997	VectorCategories::DenseVectorTag,
998	VectorCategories::DenseVectorTag) const
999	{
1000	return this->mulColDense (_VD, w, A, v);
1001	}
1002	template <class Vector1, class Matrix, class Vector2>
1003	Vector1 &mulColSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
1004	VectorCategories::DenseVectorTag,
1005	VectorCategories::SparseSequenceVectorTag) const;
1006	template <class Vector1, class Matrix, class Vector2>
1007	Vector1 &mulColSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
1008	VectorCategories::DenseVectorTag,
1009	VectorCategories::SparseAssociativeVectorTag) const;
1010	template <class Vector1, class Matrix, class Vector2>
1011	Vector1 &mulColSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
1012	VectorCategories::DenseVectorTag,
1013	VectorCategories::SparseParallelVectorTag) const;
1014
1015	template <class Vector1, class Matrix, class Vector2>
1016	inline Vector1 &mulColSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
1017	VectorCategories::GenericVectorTag,
1018	VectorCategories::GenericVectorTag) const
1019	{
1020	typename LinBox::Vector<Field>::Dense y;
1021
1022	VectorWrapper::ensureDim (y, w.size ());
1023
1024	VectorWrapper::ensureDim (y, w.size ());
1025
1026	vectorMul (y, A, v);
1027	_VD.copy (w, y);
1028
1029	return w;
1030	}
1031
1032	template <class Vector1, class Matrix, class Vector2>
1033	Vector1 &mulSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
1034	MatrixCategories::RowMatrixTag) const
1035	{
1036	return mulRowSpecialized (w, A, v, typename VectorTraits<Vector1>::VectorCategory ());
1037	}
1038	template <class Vector1, class Matrix, class Vector2>
1039	Vector1 &mulSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
1040	MatrixCategories::ColMatrixTag) const
1041	{
1042	return mulColSpecialized (w, A, v,
1043	typename VectorTraits<Vector1>::VectorCategory (),
1044	typename VectorTraits<Vector2>::VectorCategory ());
1045	}
1046	template <class Vector1, class Matrix, class Vector2>
1047	Vector1 &mulSpecialized (Vector1 &w, const Matrix &A, const Vector2 &v,
1048	MatrixCategories::RowColMatrixTag) const
1049	{
1050	return mulRowSpecialized (w, A, v, typename VectorTraits<Vector1>::VectorCategory ());
1051	}
1052
1053	template <class Vector1, class Matrix, class Vector2>
1054	Vector1 &axpyinRowSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1055	VectorCategories::DenseVectorTag) const;
1056	template <class Vector1, class Matrix, class Vector2>
1057	Vector1 &axpyinRowSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1058	VectorCategories::SparseSequenceVectorTag) const;
1059	template <class Vector1, class Matrix, class Vector2>
1060	Vector1 &axpyinRowSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1061	VectorCategories::SparseAssociativeVectorTag) const;
1062	template <class Vector1, class Matrix, class Vector2>
1063	Vector1 &axpyinRowSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1064	VectorCategories::SparseParallelVectorTag) const;
1065
1066	template <class Vector1, class Matrix, class Vector2>
1067	Vector1 &axpyinColSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1068	VectorCategories::DenseVectorTag) const;
1069	template <class Vector1, class Matrix, class Vector2>
1070	Vector1 &axpyinColSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1071	VectorCategories::SparseSequenceVectorTag) const;
1072	template <class Vector1, class Matrix, class Vector2>
1073	Vector1 &axpyinColSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1074	VectorCategories::SparseAssociativeVectorTag) const;
1075	template <class Vector1, class Matrix, class Vector2>
1076	Vector1 &axpyinColSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1077	VectorCategories::SparseParallelVectorTag) const;
1078
1079	template <class Vector1, class Matrix, class Vector2>
1080	Vector1 &axpyinSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1081	MatrixCategories::RowMatrixTag) const
1082	{
1083	return axpyinRowSpecialized (y, A, x, typename VectorTraits<Vector1>::VectorCategory ());
1084	}
1085	template <class Vector1, class Matrix, class Vector2>
1086	Vector1 &axpyinSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1087	MatrixCategories::ColMatrixTag) const
1088	{
1089	return axpyinColSpecialized (y, A, x, typename VectorTraits<Vector1>::VectorCategory ());
1090	}
1091	template <class Vector1, class Matrix, class Vector2>
1092	Vector1 &axpyinSpecialized (Vector1 &y, const Matrix &A, const Vector2 &x,
1093	MatrixCategories::RowColMatrixTag) const
1094	{
1095	return axpyinRowSpecialized (y, A, x, typename VectorTraits<Vector1>::VectorCategory ());
1096	}
1097
1098	template <class Matrix, class Iterator>
1099	inline Matrix &permuteRowsByRow (Matrix &A,
1100	Iterator P_start,
1101	Iterator P_end) const;
1102
1103	template <class Matrix, class Iterator>
1104	inline Matrix &permuteRowsByCol (Matrix &A,
1105	Iterator P_start,
1106	Iterator P_end) const;
1107
1108	template <class Matrix, class Iterator>
1109	inline Matrix &permuteRowsSpecialized (Matrix &A,
1110	Iterator P_start,
1111	Iterator P_end,
1112	MatrixCategories::RowColMatrixTag) const
1113	{
1114	return permuteRowsByCol (A, P_start, P_end);
1115	}
1116
1117	template <class Matrix, class Iterator>
1118	inline Matrix &permuteRowsSpecialized (Matrix &A,
1119	Iterator P_start,
1120	Iterator P_end,
1121	MatrixCategories::RowMatrixTag) const
1122	{
1123	return permuteRowsByRow (A, P_start, P_end);
1124	}
1125
1126	template <class Matrix, class Iterator>
1127	inline Matrix &permuteRowsSpecialized (Matrix &A,
1128	Iterator P_start,
1129	Iterator P_end,
1130	MatrixCategories::ColMatrixTag) const
1131	{
1132	return permuteRowsByCol (A, P_start, P_end);
1133	}
1134
1135	template <class Matrix, class Iterator>
1136	inline Matrix &permuteColsByRow (Matrix &A,
1137	Iterator P_start,
1138	Iterator P_end) const;
1139
1140	template <class Matrix, class Iterator>
1141	inline Matrix &permuteColsByCol (Matrix &A,
1142	Iterator P_start,
1143	Iterator P_end) const;
1144
1145	template <class Matrix, class Iterator>
1146	inline Matrix &permuteColsSpecialized (Matrix &A,
1147	Iterator P_start,
1148	Iterator P_end,
1149	MatrixCategories::RowColMatrixTag) const
1150	{
1151	return permuteColsByRow (A, P_start, P_end);
1152	}
1153
1154	template <class Matrix, class Iterator>
1155	inline Matrix &permuteColsSpecialized (Matrix &A,
1156	Iterator P_start,
1157	Iterator P_end,
1158	MatrixCategories::RowMatrixTag) const
1159	{
1160	return permuteColsByRow (A, P_start, P_end);
1161	}
1162
1163	template <class Matrix, class Iterator>
1164	inline Matrix &permuteColsSpecialized (Matrix &A,
1165	Iterator P_start,
1166	Iterator P_end,
1167	MatrixCategories::ColMatrixTag) const
1168	{
1169	return permuteColsByCol (A, P_start, P_end);
1170	}
1171
1172	Field _field;
1173	VectorDomain<Field> _VD;
1174	};
1175
1176	}
1177
1178	#include "linbox/matrix/matrix-domain.inl"
1179
1180	#endif // __LINBOX_matrix_domain_H
1181

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

Download in other formats: