matrix_multiply_test.cc File Reference

#include <iostream>
#include "generic.h"
#include "navier_stokes.h"
#include "meshes/simple_rectangular_quadmesh.h"

Classes
class	RectangularDrivenCavityProblem< ELEMENT >
	Driven cavity problem in rectangular domain. More...

Namespaces
	Global_Variables
	Namespace for physical parameters.

Functions
int	main ()

Function Documentation

main()

int main

(

)

///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// Driver for RectangularDrivenCavity test problem – test drive with two different types of element.

{
 // build the problems and
 typedef QCrouzeixRaviartElement<2> ELEMENT;
 
 Global_Variables::Nx=5;
 Global_Variables::Ny=5;
 RectangularDrivenCavityProblem<ELEMENT>* problem1_pt
  = new RectangularDrivenCavityProblem<ELEMENT>;
 problem1_pt->actions_before_newton_solve();
 
 Global_Variables::Nx=5;
 Global_Variables::Ny=5;
 RectangularDrivenCavityProblem<ELEMENT>* problem2_pt
     = new RectangularDrivenCavityProblem<ELEMENT>;
 problem2_pt->actions_before_newton_solve();
 
 // get Jacobians
 DoubleVector residual;
 CRDoubleMatrix* CR_matrix1_pt = new CRDoubleMatrix;
 CRDoubleMatrix* CR_matrix2_pt = new CRDoubleMatrix;
 CCDoubleMatrix* CC_matrix1_pt = new CCDoubleMatrix;
 CCDoubleMatrix* CC_matrix2_pt = new CCDoubleMatrix;
 problem1_pt->get_jacobian(residual,*CR_matrix1_pt);
 residual.clear();
 problem1_pt->get_jacobian(residual,*CR_matrix2_pt);
 residual.clear();
 problem1_pt->get_jacobian(residual,*CC_matrix1_pt);
 residual.clear();
 problem1_pt->get_jacobian(residual,*CC_matrix2_pt);
 residual.clear();
 
 DenseDoubleMatrix* D_matrix1_pt = new DenseDoubleMatrix;
 DenseDoubleMatrix* D_matrix2_pt = new DenseDoubleMatrix;
 problem2_pt->get_jacobian(residual,*D_matrix1_pt);
 problem2_pt->get_jacobian(residual,*D_matrix2_pt);
 
    // test matrix-matrix multiplication
    CRDoubleMatrix* CR_result1_pt = new CRDoubleMatrix;
 CRDoubleMatrix* CR_result2_pt = new CRDoubleMatrix;
 CRDoubleMatrix* CR_result3_pt = new CRDoubleMatrix;
 CCDoubleMatrix* CC_result1_pt = new CCDoubleMatrix;
 CCDoubleMatrix* CC_result2_pt = new CCDoubleMatrix;
 CCDoubleMatrix* CC_result3_pt = new CCDoubleMatrix;
 DenseDoubleMatrix* D_result_pt = new DenseDoubleMatrix;
 
 cout << "Running through matrix-matrix multiplication methods" << endl;
 cout << "----------------------------------------------------" << endl;
 
 cout << "Multiplying " << CR_matrix1_pt->ncol()
                    << " x " << CR_matrix1_pt->nrow()
      << " matrix by a " << CR_matrix2_pt->ncol()
      << " x " << CR_matrix1_pt->nrow()
      << " matrix" << endl;
 
 cout << "CRDoubleMatrix method 1 " << endl;
 CR_matrix1_pt->serial_matrix_matrix_multiply_method()=1;
 CR_matrix1_pt->multiply(*CR_matrix2_pt, *CR_result1_pt);
 
 cout << "CRDoubleMatrix method 2 " << endl;
 CR_matrix1_pt->serial_matrix_matrix_multiply_method()=2;
 CR_matrix1_pt->multiply(*CR_matrix2_pt, *CR_result2_pt);
 
 cout << "CRDoubleMatrix method 3 " << endl;
 CR_matrix1_pt->serial_matrix_matrix_multiply_method()=3;
 CR_matrix1_pt->multiply(*CR_matrix2_pt, *CR_result3_pt);
 
 cout << "CCDoubleMatrix method 1 " << endl;
 CC_matrix1_pt->matrix_matrix_multiply_method()=1;
 CC_matrix1_pt->multiply(*CC_matrix2_pt, *CC_result1_pt);
 
 cout << "CCDoubleMatrix method 2 " << endl;
 CC_matrix1_pt->matrix_matrix_multiply_method()=2;
 CC_matrix1_pt->multiply(*CC_matrix2_pt, *CC_result2_pt);
 
 cout << "CCDoubleMatrix method 3 " << endl;
 CC_matrix1_pt->matrix_matrix_multiply_method()=3;
 CC_matrix1_pt->multiply(*CC_matrix2_pt, *CC_result3_pt);
 
 cout << "Multiplying " << D_matrix1_pt->ncol()
                    << " x " << D_matrix1_pt->nrow()
      << " matrix by a " << D_matrix2_pt->ncol()
      << " x " << D_matrix1_pt->nrow()
      << " matrix" << endl;
 
 cout << "DenseDoubleMatrix " << endl;
 D_matrix1_pt->multiply(*D_matrix2_pt, *D_result_pt);
 
 // output matrices
 bool output = true;
 if (output)
 {
  ofstream file;
  file.open("CR_matrix1.dat", ios_base::out);
  CR_matrix1_pt->sparse_indexed_output(file);
  file.close();
  file.open("CR_matrix2.dat", ios_base::out);
  CR_matrix2_pt->sparse_indexed_output(file);
  file.close();
 
  file.open("CR_result1.dat", ios_base::out);
  CR_result1_pt->sparse_indexed_output(file);
  file.close();
  file.open("CR_result2.dat", ios_base::out);
  CR_result2_pt->sparse_indexed_output(file);
  file.close();
  file.open("CR_result3.dat", ios_base::out);
  CR_result3_pt->sparse_indexed_output(file);
  file.close();
 
  file.open("CC_matrix1.dat", ios_base::out);
  CC_matrix1_pt->sparse_indexed_output(file);
  file.close();
  file.open("CC_matrix2.dat", ios_base::out);
  CR_matrix2_pt->sparse_indexed_output(file);
  file.close();
 
  file.open("CC_result1.dat", ios_base::out);
  CC_result1_pt->sparse_indexed_output(file);
  file.close();
  file.open("CC_result2.dat", ios_base::out);
  CC_result2_pt->sparse_indexed_output(file);
  file.close();
  file.open("CC_result3.dat", ios_base::out);
  CC_result3_pt->sparse_indexed_output(file);
  file.close();
 
  file.open("D_matrix1.dat", ios_base::out);
  D_matrix1_pt->sparse_indexed_output(file);
  file.close();
  file.open("D_matrix2.dat", ios_base::out);
  D_matrix2_pt->sparse_indexed_output(file);
  file.close();
  file.open("D_result.dat", ios_base::out);
  D_result_pt->sparse_indexed_output(file);
  file.close();
 }
 
 // clear up memory
 delete problem1_pt;
 delete problem2_pt;
 
 delete CR_matrix1_pt;
 delete CR_matrix2_pt;
 delete CR_result1_pt;
 delete CR_result2_pt;
 delete CR_result3_pt;
 
 delete CC_matrix1_pt;
 delete CC_matrix2_pt;
 delete CC_result1_pt;
 delete CC_result2_pt;
 delete CC_result3_pt;
 
 delete D_matrix1_pt;
 delete D_matrix2_pt;
 delete D_result_pt;
 
} // end_of_main

References RectangularDrivenCavityProblem< ELEMENT >::actions_before_newton_solve(), oomph::DoubleVector::clear(), oomph::Problem::get_jacobian(), oomph::CCDoubleMatrix::matrix_matrix_multiply_method(), oomph::CRDoubleMatrix::multiply(), oomph::DenseDoubleMatrix::multiply(), oomph::CCDoubleMatrix::multiply(), oomph::CRDoubleMatrix::ncol(), oomph::DenseDoubleMatrix::ncol(), oomph::CRDoubleMatrix::nrow(), oomph::DenseDoubleMatrix::nrow(), Global_Variables::Nx, Global_Variables::Ny, out(), output(), oomph::CRDoubleMatrix::serial_matrix_matrix_multiply_method(), and oomph::Matrix< T, MATRIX_TYPE >::sparse_indexed_output().