gmm_LU_solve_interface.hh
Go to the documentation of this file.
00001 //=====================================================
00002 // File   :  blitz_LU_solve_interface.hh
00003 // Author :  L. Plagne <laurent.plagne@edf.fr)>        
00004 // Copyright (C) EDF R&D,  lun sep 30 14:23:31 CEST 2002
00005 //=====================================================
00006 // 
00007 // This program is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU General Public License
00009 // as published by the Free Software Foundation; either version 2
00010 // of the License, or (at your option) any later version.
00011 // 
00012 // This program is distributed in the hope that it will be useful,
00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 // GNU General Public License for more details.
00016 // You should have received a copy of the GNU General Public License
00017 // along with this program; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00019 // 
00020 #ifndef BLITZ_LU_SOLVE_INTERFACE_HH
00021 #define BLITZ_LU_SOLVE_INTERFACE_HH
00022 
00023 #include "blitz/array.h"
00024 #include <vector>
00025 
00026 BZ_USING_NAMESPACE(blitz)
00027 
00028 template<class real>
00029 class blitz_LU_solve_interface : public blitz_interface<real>
00030 {
00031 
00032 public :
00033 
00034   typedef typename blitz_interface<real>::gene_matrix gene_matrix;
00035   typedef typename blitz_interface<real>::gene_vector gene_vector;
00036 
00037   typedef blitz::Array<int,1> Pivot_Vector;
00038 
00039   inline static void new_Pivot_Vector(Pivot_Vector & pivot,int N)
00040   {
00041 
00042     pivot.resize(N);
00043 
00044   }
00045 
00046   inline static void free_Pivot_Vector(Pivot_Vector & pivot)
00047   {
00048     
00049     return;
00050 
00051   }
00052 
00053 
00054   static inline real matrix_vector_product_sliced(const gene_matrix & A, gene_vector B, int row, int col_start, int col_end)
00055   {
00056     
00057     real somme=0.;
00058     
00059     for (int j=col_start ; j<col_end+1 ; j++){
00060         
00061         somme+=A(row,j)*B(j);
00062         
00063     }
00064 
00065     return somme;
00066 
00067   }
00068 
00069 
00070 
00071 
00072   static inline real matrix_matrix_product_sliced(gene_matrix & A, int row, int col_start, int col_end, gene_matrix & B, int row_shift, int col )
00073   {
00074     
00075     real somme=0.;
00076     
00077     for (int j=col_start ; j<col_end+1 ; j++){
00078         
00079         somme+=A(row,j)*B(j+row_shift,col);
00080         
00081     }
00082 
00083     return somme;
00084 
00085   }
00086 
00087   inline static void LU_factor(gene_matrix & LU, Pivot_Vector & pivot, int N)
00088   {
00089 
00090     ASSERT( LU.rows()==LU.cols() ) ;
00091     int index_max = 0 ;
00092     real big = 0. ;
00093     real theSum = 0. ;
00094     real dum = 0. ;
00095     // Get the implicit scaling information :
00096     gene_vector ImplicitScaling( N ) ;
00097     for( int i=0; i<N; i++ ) {
00098       big = 0. ;
00099       for( int j=0; j<N; j++ ) {
00100         if( abs( LU( i, j ) )>=big ) big = abs( LU( i, j ) ) ;
00101       }
00102       if( big==0. ) {
00103         INFOS( "blitz_LU_factor::Singular matrix" ) ;
00104         exit( 0 ) ;
00105       }
00106       ImplicitScaling( i ) = 1./big ;
00107     }
00108     // Loop over columns of Crout's method :
00109     for( int j=0; j<N; j++ ) {
00110       for( int i=0; i<j; i++ ) {
00111         theSum = LU( i, j ) ;
00112         theSum -= matrix_matrix_product_sliced(LU, i, 0, i-1, LU, 0, j) ;
00113         //      theSum -= sum( LU( i, Range( fromStart, i-1 ) )*LU( Range( fromStart, i-1 ), j ) ) ;
00114         LU( i, j ) = theSum ;
00115       }
00116       
00117       // Search for the largest pivot element :
00118       big = 0. ;
00119       for( int i=j; i<N; i++ ) {
00120         theSum = LU( i, j ) ;
00121         theSum -= matrix_matrix_product_sliced(LU, i, 0, j-1, LU, 0, j) ;
00122         //      theSum -= sum( LU( i, Range( fromStart, j-1 ) )*LU( Range( fromStart, j-1 ), j ) ) ;
00123         LU( i, j ) = theSum ;
00124         if( (ImplicitScaling( i )*abs( theSum ))>=big ) {
00125           dum = ImplicitScaling( i )*abs( theSum ) ;
00126           big = dum ;
00127           index_max = i ;
00128         }
00129       }
00130       // Interchanging rows and the scale factor :
00131       if( j!=index_max ) {
00132         for( int k=0; k<N; k++ ) {
00133           dum = LU( index_max, k ) ;
00134           LU( index_max, k ) = LU( j, k ) ;
00135           LU( j, k ) = dum ;
00136         }
00137         ImplicitScaling( index_max ) = ImplicitScaling( j ) ;
00138       }
00139       pivot( j ) = index_max ;
00140       if ( LU( j, j )==0. ) LU( j, j ) = 1.e-20 ;
00141       // Divide by the pivot element :
00142       if( j<N ) {
00143         dum = 1./LU( j, j ) ;
00144         for( int i=j+1; i<N; i++ ) LU( i, j ) *= dum ;
00145       }
00146     }
00147 
00148   }
00149 
00150   inline static void LU_solve(const gene_matrix & LU, const Pivot_Vector pivot, gene_vector &B, gene_vector X, int N)
00151   {
00152 
00153     // Pour conserver le meme header, on travaille sur X, copie du second-membre B
00154     X = B.copy() ;
00155     ASSERT( LU.rows()==LU.cols() ) ;
00156     firstIndex indI ;
00157     // Forward substitution :
00158     int ii = 0 ;
00159     real theSum = 0. ;
00160     for( int i=0; i<N; i++ ) {
00161       int ip = pivot( i ) ;
00162       theSum = X( ip ) ;
00163       //      theSum = B( ip ) ;
00164       X( ip ) = X( i ) ;
00165       //      B( ip ) = B( i ) ;
00166       if( ii ) {
00167         theSum -= matrix_vector_product_sliced(LU, X, i, ii-1, i-1) ;
00168         //      theSum -= sum( LU( i, Range( ii-1, i-1 ) )*X( Range( ii-1, i-1 ) ) ) ;
00169         //      theSum -= sum( LU( i, Range( ii-1, i-1 ) )*B( Range( ii-1, i-1 ) ) ) ;
00170       } else if( theSum ) {
00171         ii = i+1 ;
00172       }
00173       X( i ) = theSum ;
00174       //      B( i ) = theSum ;
00175     }
00176     // Backsubstitution :
00177     for( int i=N-1; i>=0; i-- ) {
00178       theSum = X( i ) ;
00179       //      theSum = B( i ) ;
00180       theSum -= matrix_vector_product_sliced(LU, X, i, i+1, N) ;
00181       //      theSum -= sum( LU( i, Range( i+1, toEnd ) )*X( Range( i+1, toEnd ) ) ) ;
00182       //      theSum -= sum( LU( i, Range( i+1, toEnd ) )*B( Range( i+1, toEnd ) ) ) ;
00183       // Store a component of the solution vector :
00184       X( i ) = theSum/LU( i, i ) ;
00185       //      B( i ) = theSum/LU( i, i ) ;
00186     }
00187 
00188   }
00189 
00190 };
00191 
00192 #endif


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:44