| //===================================================== |
| // File : STL_interface.hh |
| // Author : L. Plagne <laurent.plagne@edf.fr)> |
| // Copyright (C) EDF R&D, lun sep 30 14:23:24 CEST 2002 |
| //===================================================== |
| // |
| // This program is free software; you can redistribute it and/or |
| // modify it under the terms of the GNU General Public License |
| // as published by the Free Software Foundation; either version 2 |
| // of the License, or (at your option) any later version. |
| // |
| // This program is distributed in the hope that it will be useful, |
| // but WITHOUT ANY WARRANTY; without even the implied warranty of |
| // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| // GNU General Public License for more details. |
| // You should have received a copy of the GNU General Public License |
| // along with this program; if not, write to the Free Software |
| // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
| // |
| #ifndef STL_INTERFACE_HH |
| #define STL_INTERFACE_HH |
| #include <string> |
| #include <vector> |
| #include "utilities.h" |
| |
| using namespace std; |
| |
| template<class real> |
| class STL_interface{ |
| |
| public : |
| |
| typedef real real_type ; |
| |
| typedef std::vector<real> stl_vector; |
| typedef std::vector<stl_vector > stl_matrix; |
| |
| typedef stl_matrix gene_matrix; |
| |
| typedef stl_vector gene_vector; |
| |
| static inline std::string name( void ) |
| { |
| return "STL"; |
| } |
| |
| static void free_matrix(gene_matrix & /*A*/, int /*N*/){} |
| |
| static void free_vector(gene_vector & /*B*/){} |
| |
| static inline void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){ |
| A = A_stl; |
| } |
| |
| static inline void vector_from_stl(gene_vector & B, stl_vector & B_stl){ |
| B = B_stl; |
| } |
| |
| static inline void vector_to_stl(gene_vector & B, stl_vector & B_stl){ |
| B_stl = B ; |
| } |
| |
| |
| static inline void matrix_to_stl(gene_matrix & A, stl_matrix & A_stl){ |
| A_stl = A ; |
| } |
| |
| static inline void copy_vector(const gene_vector & source, gene_vector & cible, int N){ |
| for (int i=0;i<N;i++){ |
| cible[i]=source[i]; |
| } |
| } |
| |
| |
| static inline void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){ |
| for (int i=0;i<N;i++) |
| for (int j=0;j<N;j++) |
| cible[i][j]=source[i][j]; |
| } |
| |
| // static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N) |
| // { |
| // real somme; |
| // for (int j=0;j<N;j++){ |
| // for (int i=0;i<N;i++){ |
| // somme=0.0; |
| // for (int k=0;k<N;k++) |
| // somme += A[i][k]*A[j][k]; |
| // X[j][i]=somme; |
| // } |
| // } |
| // } |
| |
| static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N) |
| { |
| real somme; |
| for (int j=0;j<N;j++){ |
| for (int i=0;i<N;i++){ |
| somme=0.0; |
| if(i>=j) |
| { |
| for (int k=0;k<N;k++){ |
| somme+=A[k][i]*A[k][j]; |
| } |
| X[j][i]=somme; |
| } |
| } |
| } |
| } |
| |
| |
| static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N) |
| { |
| real somme; |
| for (int j=0;j<N;j++){ |
| for (int i=0;i<N;i++){ |
| somme=0.0; |
| for (int k=0;k<N;k++) |
| somme+=A[k][i]*B[j][k]; |
| X[j][i]=somme; |
| } |
| } |
| } |
| |
| static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N) |
| { |
| real somme; |
| for (int i=0;i<N;i++){ |
| somme=0.0; |
| for (int j=0;j<N;j++) |
| somme+=A[j][i]*B[j]; |
| X[i]=somme; |
| } |
| } |
| |
| static inline void symv(gene_matrix & A, gene_vector & B, gene_vector & X, int N) |
| { |
| for (int j=0; j<N; ++j) |
| X[j] = 0; |
| for (int j=0; j<N; ++j) |
| { |
| real t1 = B[j]; |
| real t2 = 0; |
| X[j] += t1 * A[j][j]; |
| for (int i=j+1; i<N; ++i) { |
| X[i] += t1 * A[j][i]; |
| t2 += A[j][i] * B[i]; |
| } |
| X[j] += t2; |
| } |
| } |
| |
| static inline void syr2(gene_matrix & A, gene_vector & B, gene_vector & X, int N) |
| { |
| for (int j=0; j<N; ++j) |
| { |
| for (int i=j; i<N; ++i) |
| A[j][i] += B[i]*X[j] + B[j]*X[i]; |
| } |
| } |
| |
| static inline void ger(gene_matrix & A, gene_vector & X, gene_vector & Y, int N) |
| { |
| for (int j=0; j<N; ++j) |
| { |
| for (int i=j; i<N; ++i) |
| A[j][i] += X[i]*Y[j]; |
| } |
| } |
| |
| static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N) |
| { |
| real somme; |
| for (int i=0;i<N;i++){ |
| somme = 0.0; |
| for (int j=0;j<N;j++) |
| somme += A[i][j]*B[j]; |
| X[i] = somme; |
| } |
| } |
| |
| static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N){ |
| for (int i=0;i<N;i++) |
| Y[i]+=coef*X[i]; |
| } |
| |
| static inline void axpby(real a, const gene_vector & X, real b, gene_vector & Y, int N){ |
| for (int i=0;i<N;i++) |
| Y[i] = a*X[i] + b*Y[i]; |
| } |
| |
| static inline void trisolve_lower(const gene_matrix & L, const gene_vector & B, gene_vector & X, int N){ |
| copy_vector(B,X,N); |
| for(int i=0; i<N; ++i) |
| { |
| X[i] /= L[i][i]; |
| real tmp = X[i]; |
| for (int j=i+1; j<N; ++j) |
| X[j] -= tmp * L[i][j]; |
| } |
| } |
| |
| static inline real norm_diff(const stl_vector & A, const stl_vector & B) |
| { |
| int N=A.size(); |
| real somme=0.0; |
| real somme2=0.0; |
| |
| for (int i=0;i<N;i++){ |
| real diff=A[i]-B[i]; |
| somme+=diff*diff; |
| somme2+=A[i]*A[i]; |
| } |
| return somme/somme2; |
| } |
| |
| static inline real norm_diff(const stl_matrix & A, const stl_matrix & B) |
| { |
| int N=A[0].size(); |
| real somme=0.0; |
| real somme2=0.0; |
| |
| for (int i=0;i<N;i++){ |
| for (int j=0;j<N;j++){ |
| real diff=A[i][j] - B[i][j]; |
| somme += diff*diff; |
| somme2 += A[i][j]*A[i][j]; |
| } |
| } |
| |
| return somme/somme2; |
| } |
| |
| static inline void display_vector(const stl_vector & A) |
| { |
| int N=A.size(); |
| for (int i=0;i<N;i++){ |
| INFOS("A["<<i<<"]="<<A[i]<<endl); |
| } |
| } |
| |
| }; |
| |
| #endif |