1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath//===================================================== 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// File : action_lu_solve.hh 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Author : L. Plagne <laurent.plagne@edf.fr)> 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) EDF R&D, lun sep 30 14:23:19 CEST 2002 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath//===================================================== 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This program is free software; you can redistribute it and/or 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// modify it under the terms of the GNU General Public License 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// as published by the Free Software Foundation; either version 2 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// of the License, or (at your option) any later version. 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This program is distributed in the hope that it will be useful, 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// but WITHOUT ANY WARRANTY; without even the implied warranty of 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// GNU General Public License for more details. 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// You should have received a copy of the GNU General Public License 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// along with this program; if not, write to the Free Software 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef ACTION_LU_SOLVE 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define ACTION_LU_SOLVE 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "utilities.h" 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "STL_interface.hh" 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <string> 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "init/init_function.hh" 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "init/init_vector.hh" 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "init/init_matrix.hh" 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathusing namespace std; 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<class Interface> 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass Action_lu_solve 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic : 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline std::string name( void ) 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return "lu_solve_"+Interface::name(); 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static double nb_op_base(int size){ 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 2.0*size*size*size/3.0; // questionable but not really important 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static double calculate( int nb_calc, int size ) { 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // STL matrix and vector initialization 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Interface::stl_matrix A_stl; 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Interface::stl_vector B_stl; 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Interface::stl_vector X_stl; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init_matrix<pseudo_random>(A_stl,size); 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init_vector<pseudo_random>(B_stl,size); 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init_vector<null_function>(X_stl,size); 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // generic matrix and vector initialization 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Interface::gene_matrix A; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Interface::gene_vector B; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Interface::gene_vector X; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Interface::gene_matrix LU; 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::matrix_from_stl(A,A_stl); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::vector_from_stl(B,B_stl); 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::vector_from_stl(X,X_stl); 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::matrix_from_stl(LU,A_stl); 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // local variable : 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Interface::Pivot_Vector pivot; // pivot vector 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::new_Pivot_Vector(pivot,size); 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // timer utilities 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Portable_Timer chronos; 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // time measurement 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath chronos.start(); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int ii=0;ii<nb_calc;ii++){ 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // LU factorization 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::copy_matrix(A,LU,size); 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::LU_factor(LU,pivot,size); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // LU solve 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::LU_solve(LU,pivot,B,X,size); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Time stop 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath chronos.stop(); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double time=chronos.user_time(); 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check result : 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Interface::stl_vector B_new_stl(size); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::vector_to_stl(X,X_stl); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath STL_interface<typename Interface::real_type>::matrix_vector_product(A_stl,X_stl,B_new_stl,size); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Interface::real_type error= 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath STL_interface<typename Interface::real_type>::norm_diff(B_stl,B_new_stl); 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (error>1.e-5){ 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath INFOS("WRONG CALCULATION...residual=" << error); 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath STL_interface<typename Interface::real_type>::display_vector(B_stl); 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath STL_interface<typename Interface::real_type>::display_vector(B_new_stl); 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath exit(0); 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // deallocation and return time 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::free_matrix(A,size); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::free_vector(B); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::free_vector(X); 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Interface::free_Pivot_Vector(pivot); 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return time; 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 137