1 | //=====================================================
|
---|
2 | // File : STL_interface.hh
|
---|
3 | // Author : L. Plagne <laurent.plagne@edf.fr)>
|
---|
4 | // Copyright (C) EDF R&D, lun sep 30 14:23:24 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 STL_INTERFACE_HH
|
---|
21 | #define STL_INTERFACE_HH
|
---|
22 | #include <string>
|
---|
23 | #include <vector>
|
---|
24 | #include "utilities.h"
|
---|
25 |
|
---|
26 | using namespace std;
|
---|
27 |
|
---|
28 | template<class real>
|
---|
29 | class STL_interface{
|
---|
30 |
|
---|
31 | public :
|
---|
32 |
|
---|
33 | typedef real real_type ;
|
---|
34 |
|
---|
35 | typedef std::vector<real> stl_vector;
|
---|
36 | typedef std::vector<stl_vector > stl_matrix;
|
---|
37 |
|
---|
38 | typedef stl_matrix gene_matrix;
|
---|
39 |
|
---|
40 | typedef stl_vector gene_vector;
|
---|
41 |
|
---|
42 | static inline std::string name( void )
|
---|
43 | {
|
---|
44 | return "STL";
|
---|
45 | }
|
---|
46 |
|
---|
47 | static void free_matrix(gene_matrix & A, int N){}
|
---|
48 |
|
---|
49 | static void free_vector(gene_vector & B){}
|
---|
50 |
|
---|
51 | static inline void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){
|
---|
52 | A = A_stl;
|
---|
53 | }
|
---|
54 |
|
---|
55 | static inline void vector_from_stl(gene_vector & B, stl_vector & B_stl){
|
---|
56 | B = B_stl;
|
---|
57 | }
|
---|
58 |
|
---|
59 | static inline void vector_to_stl(gene_vector & B, stl_vector & B_stl){
|
---|
60 | B_stl = B ;
|
---|
61 | }
|
---|
62 |
|
---|
63 |
|
---|
64 | static inline void matrix_to_stl(gene_matrix & A, stl_matrix & A_stl){
|
---|
65 | A_stl = A ;
|
---|
66 | }
|
---|
67 |
|
---|
68 | static inline void copy_vector(const gene_vector & source, gene_vector & cible, int N){
|
---|
69 | for (int i=0;i<N;i++){
|
---|
70 | cible[i]=source[i];
|
---|
71 | }
|
---|
72 | }
|
---|
73 |
|
---|
74 |
|
---|
75 | static inline void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){
|
---|
76 | for (int i=0;i<N;i++)
|
---|
77 | for (int j=0;j<N;j++)
|
---|
78 | cible[i][j]=source[i][j];
|
---|
79 | }
|
---|
80 |
|
---|
81 | // static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N)
|
---|
82 | // {
|
---|
83 | // real somme;
|
---|
84 | // for (int j=0;j<N;j++){
|
---|
85 | // for (int i=0;i<N;i++){
|
---|
86 | // somme=0.0;
|
---|
87 | // for (int k=0;k<N;k++)
|
---|
88 | // somme += A[i][k]*A[j][k];
|
---|
89 | // X[j][i]=somme;
|
---|
90 | // }
|
---|
91 | // }
|
---|
92 | // }
|
---|
93 |
|
---|
94 | static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N)
|
---|
95 | {
|
---|
96 | real somme;
|
---|
97 | for (int j=0;j<N;j++){
|
---|
98 | for (int i=0;i<N;i++){
|
---|
99 | somme=0.0;
|
---|
100 | if(i>=j)
|
---|
101 | {
|
---|
102 | for (int k=0;k<N;k++){
|
---|
103 | somme+=A[k][i]*A[k][j];
|
---|
104 | }
|
---|
105 | X[j][i]=somme;
|
---|
106 | }
|
---|
107 | }
|
---|
108 | }
|
---|
109 | }
|
---|
110 |
|
---|
111 |
|
---|
112 | static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N)
|
---|
113 | {
|
---|
114 | real somme;
|
---|
115 | for (int j=0;j<N;j++){
|
---|
116 | for (int i=0;i<N;i++){
|
---|
117 | somme=0.0;
|
---|
118 | for (int k=0;k<N;k++)
|
---|
119 | somme+=A[k][i]*B[j][k];
|
---|
120 | X[j][i]=somme;
|
---|
121 | }
|
---|
122 | }
|
---|
123 | }
|
---|
124 |
|
---|
125 | static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
|
---|
126 | {
|
---|
127 | real somme;
|
---|
128 | for (int i=0;i<N;i++){
|
---|
129 | somme=0.0;
|
---|
130 | for (int j=0;j<N;j++)
|
---|
131 | somme+=A[j][i]*B[j];
|
---|
132 | X[i]=somme;
|
---|
133 | }
|
---|
134 | }
|
---|
135 |
|
---|
136 | static inline void symv(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
|
---|
137 | {
|
---|
138 | for (int j=0; j<N; ++j)
|
---|
139 | X[j] = 0;
|
---|
140 | for (int j=0; j<N; ++j)
|
---|
141 | {
|
---|
142 | real t1 = B[j];
|
---|
143 | real t2 = 0;
|
---|
144 | X[j] += t1 * A[j][j];
|
---|
145 | for (int i=j+1; i<N; ++i) {
|
---|
146 | X[i] += t1 * A[j][i];
|
---|
147 | t2 += A[j][i] * B[i];
|
---|
148 | }
|
---|
149 | X[j] += t2;
|
---|
150 | }
|
---|
151 | }
|
---|
152 |
|
---|
153 | static inline void syr2(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
|
---|
154 | {
|
---|
155 | for (int j=0; j<N; ++j)
|
---|
156 | {
|
---|
157 | for (int i=j; i<N; ++i)
|
---|
158 | A[j][i] += B[i]*X[j] + B[j]*X[i];
|
---|
159 | }
|
---|
160 | }
|
---|
161 |
|
---|
162 | static inline void ger(gene_matrix & A, gene_vector & X, gene_vector & Y, int N)
|
---|
163 | {
|
---|
164 | for (int j=0; j<N; ++j)
|
---|
165 | {
|
---|
166 | for (int i=j; i<N; ++i)
|
---|
167 | A[j][i] += X[i]*Y[j];
|
---|
168 | }
|
---|
169 | }
|
---|
170 |
|
---|
171 | static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
|
---|
172 | {
|
---|
173 | real somme;
|
---|
174 | for (int i=0;i<N;i++){
|
---|
175 | somme = 0.0;
|
---|
176 | for (int j=0;j<N;j++)
|
---|
177 | somme += A[i][j]*B[j];
|
---|
178 | X[i] = somme;
|
---|
179 | }
|
---|
180 | }
|
---|
181 |
|
---|
182 | static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N){
|
---|
183 | for (int i=0;i<N;i++)
|
---|
184 | Y[i]+=coef*X[i];
|
---|
185 | }
|
---|
186 |
|
---|
187 | static inline void axpby(real a, const gene_vector & X, real b, gene_vector & Y, int N){
|
---|
188 | for (int i=0;i<N;i++)
|
---|
189 | Y[i] = a*X[i] + b*Y[i];
|
---|
190 | }
|
---|
191 |
|
---|
192 | static inline void trisolve_lower(const gene_matrix & L, const gene_vector & B, gene_vector & X, int N){
|
---|
193 | copy_vector(B,X,N);
|
---|
194 | for(int i=0; i<N; ++i)
|
---|
195 | {
|
---|
196 | X[i] /= L[i][i];
|
---|
197 | real tmp = X[i];
|
---|
198 | for (int j=i+1; j<N; ++j)
|
---|
199 | X[j] -= tmp * L[i][j];
|
---|
200 | }
|
---|
201 | }
|
---|
202 |
|
---|
203 | static inline real norm_diff(const stl_vector & A, const stl_vector & B)
|
---|
204 | {
|
---|
205 | int N=A.size();
|
---|
206 | real somme=0.0;
|
---|
207 | real somme2=0.0;
|
---|
208 |
|
---|
209 | for (int i=0;i<N;i++){
|
---|
210 | real diff=A[i]-B[i];
|
---|
211 | somme+=diff*diff;
|
---|
212 | somme2+=A[i]*A[i];
|
---|
213 | }
|
---|
214 | return somme/somme2;
|
---|
215 | }
|
---|
216 |
|
---|
217 | static inline real norm_diff(const stl_matrix & A, const stl_matrix & B)
|
---|
218 | {
|
---|
219 | int N=A[0].size();
|
---|
220 | real somme=0.0;
|
---|
221 | real somme2=0.0;
|
---|
222 |
|
---|
223 | for (int i=0;i<N;i++){
|
---|
224 | for (int j=0;j<N;j++){
|
---|
225 | real diff=A[i][j] - B[i][j];
|
---|
226 | somme += diff*diff;
|
---|
227 | somme2 += A[i][j]*A[i][j];
|
---|
228 | }
|
---|
229 | }
|
---|
230 |
|
---|
231 | return somme/somme2;
|
---|
232 | }
|
---|
233 |
|
---|
234 | static inline void display_vector(const stl_vector & A)
|
---|
235 | {
|
---|
236 | int N=A.size();
|
---|
237 | for (int i=0;i<N;i++){
|
---|
238 | INFOS("A["<<i<<"]="<<A[i]<<endl);
|
---|
239 | }
|
---|
240 | }
|
---|
241 |
|
---|
242 | };
|
---|
243 |
|
---|
244 | #endif
|
---|