<commit_before>// $Id: dense_matrix.C,v 1.26 2005-06-12 18:36:41 jwpeterson Exp $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2005  <NAME>, <NAME>
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


// C++ Includes
#include <cmath> // for sqrt

// Local Includes
#include "dense_matrix.h"
#include "dense_vector.h"
#include "libmesh.h"



// ------------------------------------------------------------
// Dense Matrix member functions
template<typename T>
void DenseMatrix<T>::left_multiply (const DenseMatrixBase<T>& M2)
{
  // (*this) <- M2 * M3
  // Where:
  // (*this) = (m x n),
  // M2      = (m x p),
  // M3      = (p x n)

  // M3 is a copy of *this before it gets resize()d
  DenseMatrix<T> M3(*this);

  // Resize *this so that the result can fit
  this->resize (M2.m(), M3.n());

  // Call the multiply function in the base class
  this->multiply(*this, M2, M3);
}




template<typename T>
void DenseMatrix<T>::left_multiply_transpose(const DenseMatrix<T>& A)
{

  DenseMatrix<T> B(*this);
  
  this->resize (A.n(), B.n());
      
  assert (A.m() == B.m());
  assert (this->m() == A.n());
  assert (this->n() == B.n());
      
  const unsigned int m_s = A.n();
  const unsigned int p_s = A.m(); 
  const unsigned int n_s = this->n();
  
  // Do it this way because there is a
  // decent chance (at least for constraint matrices)
  // that A.transpose(i,k) = 0.
  for (unsigned int i=0; i<m_s; i++)
    for (unsigned int k=0; k<p_s; k++)
      if (A.transpose(i,k) != 0.)
	for (unsigned int j=0; j<n_s; j++)
	  (*this)(i,j) += A.transpose(i,k)*B(k,j);  

}






template<typename T>
void DenseMatrix<T>::right_multiply (const DenseMatrixBase<T>& M3)
{
  // (*this) <- M2 * M3
  // Where:
  // (*this) = (m x n),
  // M2      = (m x p),
  // M3      = (p x n)

  // M2 is a copy of *this before it gets resize()d
  DenseMatrix<T> M2(*this);

  // Resize *this so that the result can fit
  this->resize (M2.m(), M3.n());
  
  this->multiply(*this, M2, M3);
}




template<typename T>
void DenseMatrix<T>::right_multiply_transpose (const DenseMatrix<T>& B)
{
  DenseMatrix<T> A(*this);
  
  this->resize (A.m(), B.m());
      
  assert (A.n() == B.n());
  assert (this->m() == A.m());
  assert (this->n() == B.m());
      
  const unsigned int m_s = A.m();
  const unsigned int p_s = A.n(); 
  const unsigned int n_s = this->n();

  // Do it this way because there is a
  // decent chance (at least for constraint matrices)
  // that B.transpose(k,j) = 0.
  for (unsigned int j=0; j<n_s; j++)
    for (unsigned int k=0; k<p_s; k++)
      if (B.transpose(k,j) != 0.)
	for (unsigned int i=0; i<m_s; i++)
	  (*this)(i,j) += A(i,k)*B.transpose(k,j);	            
}



template<typename T>
void DenseMatrix<T>::lu_solve (DenseVector<T>& b,
			       DenseVector<T>& x,
			       const bool partial_pivot)
{
  // Check for a previous decomposition
  switch(this->_decomposition_type)
    {
    case NONE:
      {
	this->_lu_decompose (partial_pivot);
	break;
      }
<commit_msg>Fixed left/right_multiply_transpose(A) to work in the case where (this==&A).<commit_after>// $Id: dense_matrix.C,v 1.27 2007-01-23 00:50:11 jwpeterson Exp $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2005  <NAME>, <NAME>
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


// C++ Includes
#include <cmath> // for sqrt

// Local Includes
#include "dense_matrix.h"
#include "dense_vector.h"
#include "libmesh.h"



// ------------------------------------------------------------
// Dense Matrix member functions
template<typename T>
void DenseMatrix<T>::left_multiply (const DenseMatrixBase<T>& M2)
{
  // (*this) <- M2 * M3
  // Where:
  // (*this) = (m x n),
  // M2      = (m x p),
  // M3      = (p x n)

  // M3 is a copy of *this before it gets resize()d
  DenseMatrix<T> M3(*this);

  // Resize *this so that the result can fit
  this->resize (M2.m(), M3.n());

  // Call the multiply function in the base class
  this->multiply(*this, M2, M3);
}




template<typename T>
void DenseMatrix<T>::left_multiply_transpose(const DenseMatrix<T>& A)
{
  //Check to see if we are doing (A^T)*A
  if (this == &A)
    {
      //here();
      DenseMatrix<T> B(*this);
      
      // Simple but inefficient way
      // return this->left_multiply_transpose(B);

      // More efficient, but more code way
      // If A is mxn, the result will be a square matrix of Size n x n.
      const unsigned int m = A.m();
      const unsigned int n = A.n();

      // resize() *this and also zero out all entries.
      this->resize(n,n);

      // Compute the lower-triangular part
      for (unsigned int i=0; i<n; ++i)
	for (unsigned int j=0; j<=i; ++j)
	  for (unsigned int k=0; k<m; ++k) // inner products are over m
	    (*this)(i,j) += B(k,i)*B(k,j);

      // Copy lower-triangular part into upper-triangular part
      for (unsigned int i=0; i<n; ++i)
	for (unsigned int j=i+1; j<n; ++j)
	  (*this)(i,j) = (*this)(j,i);
    }

  else
    {
      DenseMatrix<T> B(*this);
  
      this->resize (A.n(), B.n());
      
      assert (A.m() == B.m());
      assert (this->m() == A.n());
      assert (this->n() == B.n());
      
      const unsigned int m_s = A.n();
      const unsigned int p_s = A.m(); 
      const unsigned int n_s = this->n();
  
      // Do it this way because there is a
      // decent chance (at least for constraint matrices)
      // that A.transpose(i,k) = 0.
      for (unsigned int i=0; i<m_s; i++)
	for (unsigned int k=0; k<p_s; k++)
	  if (A.transpose(i,k) != 0.)
	    for (unsigned int j=0; j<n_s; j++)
	      (*this)(i,j) += A.transpose(i,k)*B(k,j);
    }

}






template<typename T>
void DenseMatrix<T>::right_multiply (const DenseMatrixBase<T>& M3)
{
  // (*this) <- M2 * M3
  // Where:
  // (*this) = (m x n),
  // M2      = (m x p),
  // M3      = (p x n)

  // M2 is a copy of *this before it gets resize()d
  DenseMatrix<T> M2(*this);

  // Resize *this so that the result can fit
  this->resize (M2.m(), M3.n());
  
  this->multiply(*this, M2, M3);
}




template<typename T>
void DenseMatrix<T>::right_multiply_transpose (const DenseMatrix<T>& B)
{
  //Check to see if we are doing B*(B^T)
  if (this == &B)
    {
      //here();
      DenseMatrix<T> A(*this);
      
      // Simple but inefficient way
      // return this->right_multiply_transpose(A);

      // More efficient, more code
      // If B is mxn, the result will be a square matrix of Size m x m.
      const unsigned int m = B.m();
      const unsigned int n = B.n();

      // resize() *this and also zero out all entries.
      this->resize(m,m);

      // Compute the lower-triangular part
      for (unsigned int i=0; i<m; ++i)
	for (unsigned int j=0; j<=i; ++j)
	  for (unsigned int k=0; k<n; ++k) // inner products are over n
	    (*this)(i,j) += A(i,k)*A(j,k);

      // Copy lower-triangular part into upper-triangular part
      for (unsigned int i=0; i<m; ++i)
	for (unsigned int j=i+1; j<m; ++j)
	  (*this)(i,j) = (*this)(j,i);
    }

  else
    {
      DenseMatrix<T> A(*this);
  
      this->resize (A.m(), B.m());
      
      assert (A.n() == B.n());
      assert (this->m() == A.m());
      assert (this->n() == B.m());
      
      const unsigned int m_s = A.m();
      const unsigned int p_s = A.n(); 
      const unsigned int n_s = this->n();

      // Do it this way because there is a
      // decent chance (at least for constraint matrices)
      // that B.transpose(k,j) = 0.
      for (unsigned int j=0; j<n_s; j++)
	for (unsigned int k=0; k<p_s; k++)
	  if (B.transpose(k,j) != 0.)
	    for (unsigned int i=0; i<m_s; i++)
	      (*this)(i,j) += A(i,k)*B.transpose(k,j);
    }
}



template<typename T>
void DenseMatrix<T>::lu_solve (DenseVector<T>& b,
			       DenseVector<T>& x,
			       const bool partial_pivot)
{
  // Check for a previous decomposition
  switch(this->_decomposition_type)
    {
    case NONE:
      {
	this->_lu_decompose (partial_pivot);
	break;
      }

{'fraction_non_alphanumeric': 0.16996344872070523, 'fraction_numerical': 0.01677058697054397, 'mean_word_length': 2.4725643896976486, 'pattern_counts': {'":': 0, '<': 77, '<?xml version=': 0, '>': 83, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4909976', 'n_tokens_mistral': 3441, 'n_tokens_neox': 3267, 'n_words': 1094}