ATL / WTL / STL

#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
using namespace std;

template <class T>
void OUT(T c)
{
        cout << "!!" << c ;
}


int main()
{
        vector<string> st(10, " ");

        for_each(st.begin(), st.end(), OUT);
}

template <class T> class vector;  // vector space
template <class T> class matrix;  // derived from vector
template <typename T> T norm(vector<T> v);
template <typename T> T norm(matrix<T> m);
template <typename T> void transpose(matrix<T> &m);
template <typename T> void conjugate_transpose(matrix<T> &m);
template <typename T> bool hermetian(matrix<T> &a, matrix<T> &b);
template <typename T> matrix<T> tensor_product(vector<T> row, vector<T> col);
template <typename T> matrix<T> minor(matrix<T> &source, int row, int col);
template <typename T> bool diagonal_dominant(matrix<T> a);
template <typename T> matrix<T> gauss_elimination(matrix<T> &m);
template <typename T> matrix<T> gauss_jordan(matrix<T> &m);
template <typename T> void gauss_siedel(matrix<T> a, vector<T> b, vector<T> x, double tolerance);

template <typename T> matrix<T> gauss_elimination(matrix<T> &mat) {
	int i = 0; int j = 0; int maxi;
	matrix<T> a = mat;
	int m = (int)a.size();
	int n = (int)a[0].size();
	while ((i < m) && (j < n)) {
		maxi = i;
		// sweep down the column to find the pivot element
		for (int k = i + 1; k < m; k++) // partial pivoting
			if (std::abs(a[k][j]) > std::abs(a[maxi][j])) maxi = k;
		// eliminate
		if (a[maxi][j] != 0) {
			swap(a[i], a[maxi]); // row op 1, swap
			a[i] = a[i] / a[i][j];  // row op 2, scale the pivot to 1
			for (int u = i + 1; u < m; u++) 
				a[u] = a[u] - (a[i] * a[u][j]);  // row op 3, zero out the column below
			i++;
		}
		j++;
	}
	return a;
}