//output += trans(input)*input
		void add_symmetric_prod (const TSystemMatrixType& input,
			TSystemMatrixType& output) 
		{
			KRATOS_TRY
				typedef  unsigned int size_type;
			typedef  double value_type;


			for (size_type k = 0; k < input.size1 (); ++ k) 
			{
				size_type begin = input.index1_data () [k];
				size_type end = input.index1_data () [k + 1];


				for (size_type i = begin; i < end; ++ i)
				{
					unsigned int index_i = input.index2_data () [i];
					value_type data_i = input.value_data()[i];
					for (size_type j = begin; j < end; ++ j)
					{
						unsigned int index_j = input.index2_data () [j];
						//std::cout << index_i << " " << index_j << std::endl;
						output(index_i,index_j) += data_i * input.value_data()[j];
					}
				}
				//KRATOS_WATCH(output.non_zeros());
				//std::cout << k << " nnz = " << end << " " << begin << std::endl;
			}
			KRATOS_CATCH("");
		}