1 //===================================================== 2 // File : blitz_LU_solve_interface.hh 3 // Author : L. Plagne <laurent.plagne@edf.fr)> 4 // Copyright (C) EDF R&D, lun sep 30 14:23:31 CEST 2002 5 //===================================================== 6 // 7 // This program is free software; you can redistribute it and/or 8 // modify it under the terms of the GNU General Public License 9 // as published by the Free Software Foundation; either version 2 10 // of the License, or (at your option) any later version. 11 // 12 // This program is distributed in the hope that it will be useful, 13 // but WITHOUT ANY WARRANTY; without even the implied warranty of 14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 // GNU General Public License for more details. 16 // You should have received a copy of the GNU General Public License 17 // along with this program; if not, write to the Free Software 18 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 19 // 20 #ifndef BLITZ_LU_SOLVE_INTERFACE_HH 21 #define BLITZ_LU_SOLVE_INTERFACE_HH 22 23 #include "blitz/array.h" 24 #include <vector> 25 26 BZ_USING_NAMESPACE(blitz) 27 28 template<class real> 29 class blitz_LU_solve_interface : public blitz_interface<real> 30 { 31 32 public : 33 34 typedef typename blitz_interface<real>::gene_matrix gene_matrix; 35 typedef typename blitz_interface<real>::gene_vector gene_vector; 36 37 typedef blitz::Array<int,1> Pivot_Vector; 38 new_Pivot_Vector(Pivot_Vector & pivot,int N)39 inline static void new_Pivot_Vector(Pivot_Vector & pivot,int N) 40 { 41 42 pivot.resize(N); 43 44 } 45 free_Pivot_Vector(Pivot_Vector & pivot)46 inline static void free_Pivot_Vector(Pivot_Vector & pivot) 47 { 48 49 return; 50 51 } 52 53 matrix_vector_product_sliced(const gene_matrix & A,gene_vector B,int row,int col_start,int col_end)54 static inline real matrix_vector_product_sliced(const gene_matrix & A, gene_vector B, int row, int col_start, int col_end) 55 { 56 57 real somme=0.; 58 59 for (int j=col_start ; j<col_end+1 ; j++){ 60 61 somme+=A(row,j)*B(j); 62 63 } 64 65 return somme; 66 67 } 68 69 70 71 matrix_matrix_product_sliced(gene_matrix & A,int row,int col_start,int col_end,gene_matrix & B,int row_shift,int col)72 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 ) 73 { 74 75 real somme=0.; 76 77 for (int j=col_start ; j<col_end+1 ; j++){ 78 79 somme+=A(row,j)*B(j+row_shift,col); 80 81 } 82 83 return somme; 84 85 } 86 LU_factor(gene_matrix & LU,Pivot_Vector & pivot,int N)87 inline static void LU_factor(gene_matrix & LU, Pivot_Vector & pivot, int N) 88 { 89 90 ASSERT( LU.rows()==LU.cols() ) ; 91 int index_max = 0 ; 92 real big = 0. ; 93 real theSum = 0. ; 94 real dum = 0. ; 95 // Get the implicit scaling information : 96 gene_vector ImplicitScaling( N ) ; 97 for( int i=0; i<N; i++ ) { 98 big = 0. ; 99 for( int j=0; j<N; j++ ) { 100 if( abs( LU( i, j ) )>=big ) big = abs( LU( i, j ) ) ; 101 } 102 if( big==0. ) { 103 INFOS( "blitz_LU_factor::Singular matrix" ) ; 104 exit( 0 ) ; 105 } 106 ImplicitScaling( i ) = 1./big ; 107 } 108 // Loop over columns of Crout's method : 109 for( int j=0; j<N; j++ ) { 110 for( int i=0; i<j; i++ ) { 111 theSum = LU( i, j ) ; 112 theSum -= matrix_matrix_product_sliced(LU, i, 0, i-1, LU, 0, j) ; 113 // theSum -= sum( LU( i, Range( fromStart, i-1 ) )*LU( Range( fromStart, i-1 ), j ) ) ; 114 LU( i, j ) = theSum ; 115 } 116 117 // Search for the largest pivot element : 118 big = 0. ; 119 for( int i=j; i<N; i++ ) { 120 theSum = LU( i, j ) ; 121 theSum -= matrix_matrix_product_sliced(LU, i, 0, j-1, LU, 0, j) ; 122 // theSum -= sum( LU( i, Range( fromStart, j-1 ) )*LU( Range( fromStart, j-1 ), j ) ) ; 123 LU( i, j ) = theSum ; 124 if( (ImplicitScaling( i )*abs( theSum ))>=big ) { 125 dum = ImplicitScaling( i )*abs( theSum ) ; 126 big = dum ; 127 index_max = i ; 128 } 129 } 130 // Interchanging rows and the scale factor : 131 if( j!=index_max ) { 132 for( int k=0; k<N; k++ ) { 133 dum = LU( index_max, k ) ; 134 LU( index_max, k ) = LU( j, k ) ; 135 LU( j, k ) = dum ; 136 } 137 ImplicitScaling( index_max ) = ImplicitScaling( j ) ; 138 } 139 pivot( j ) = index_max ; 140 if ( LU( j, j )==0. ) LU( j, j ) = 1.e-20 ; 141 // Divide by the pivot element : 142 if( j<N ) { 143 dum = 1./LU( j, j ) ; 144 for( int i=j+1; i<N; i++ ) LU( i, j ) *= dum ; 145 } 146 } 147 148 } 149 LU_solve(const gene_matrix & LU,const Pivot_Vector pivot,gene_vector & B,gene_vector X,int N)150 inline static void LU_solve(const gene_matrix & LU, const Pivot_Vector pivot, gene_vector &B, gene_vector X, int N) 151 { 152 153 // Pour conserver le meme header, on travaille sur X, copie du second-membre B 154 X = B.copy() ; 155 ASSERT( LU.rows()==LU.cols() ) ; 156 firstIndex indI ; 157 // Forward substitution : 158 int ii = 0 ; 159 real theSum = 0. ; 160 for( int i=0; i<N; i++ ) { 161 int ip = pivot( i ) ; 162 theSum = X( ip ) ; 163 // theSum = B( ip ) ; 164 X( ip ) = X( i ) ; 165 // B( ip ) = B( i ) ; 166 if( ii ) { 167 theSum -= matrix_vector_product_sliced(LU, X, i, ii-1, i-1) ; 168 // theSum -= sum( LU( i, Range( ii-1, i-1 ) )*X( Range( ii-1, i-1 ) ) ) ; 169 // theSum -= sum( LU( i, Range( ii-1, i-1 ) )*B( Range( ii-1, i-1 ) ) ) ; 170 } else if( theSum ) { 171 ii = i+1 ; 172 } 173 X( i ) = theSum ; 174 // B( i ) = theSum ; 175 } 176 // Backsubstitution : 177 for( int i=N-1; i>=0; i-- ) { 178 theSum = X( i ) ; 179 // theSum = B( i ) ; 180 theSum -= matrix_vector_product_sliced(LU, X, i, i+1, N) ; 181 // theSum -= sum( LU( i, Range( i+1, toEnd ) )*X( Range( i+1, toEnd ) ) ) ; 182 // theSum -= sum( LU( i, Range( i+1, toEnd ) )*B( Range( i+1, toEnd ) ) ) ; 183 // Store a component of the solution vector : 184 X( i ) = theSum/LU( i, i ) ; 185 // B( i ) = theSum/LU( i, i ) ; 186 } 187 188 } 189 190 }; 191 192 #endif 193