00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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
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
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
00114 LU( i, j ) = theSum ;
00115 }
00116
00117
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
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
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
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
00154 X = B.copy() ;
00155 ASSERT( LU.rows()==LU.cols() ) ;
00156 firstIndex indI ;
00157
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
00164 X( ip ) = X( i ) ;
00165
00166 if( ii ) {
00167 theSum -= matrix_vector_product_sliced(LU, X, i, ii-1, i-1) ;
00168
00169
00170 } else if( theSum ) {
00171 ii = i+1 ;
00172 }
00173 X( i ) = theSum ;
00174
00175 }
00176
00177 for( int i=N-1; i>=0; i-- ) {
00178 theSum = X( i ) ;
00179
00180 theSum -= matrix_vector_product_sliced(LU, X, i, i+1, N) ;
00181
00182
00183
00184 X( i ) = theSum/LU( i, i ) ;
00185
00186 }
00187
00188 }
00189
00190 };
00191
00192 #endif