/[escript]/trunk/paso/src/FCTSolver_util.c
ViewVC logotype

Contents of /trunk/paso/src/FCTSolver_util.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3642 - (show annotations)
Thu Oct 27 03:41:51 2011 UTC (7 years, 9 months ago) by caltinay
File MIME type: text/plain
File size: 18923 byte(s)
Assorted spelling/comment fixes in paso.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Paso: FluxControl */
18
19 /**************************************************************/
20
21 /* Author: l.gross@uq.edu.au */
22
23 /**************************************************************/
24
25
26 #include "Paso.h"
27 #include "FCTSolver.h"
28 #include "PasoUtil.h"
29
30 /**************************************************************/
31
32 /* Creates the low order transport matrix and stores its negative values
33 * into the iteration_matrix except for the main diagonal which is stored
34 * separately.
35 * If fc->iteration_matrix==NULL, fc->iteration_matrix is allocated
36 *
37 * a=transport_matrix
38 * b= low_order_transport_matrix = - iteration_matrix
39 * c=main diagonal low_order_transport_matrix
40 * initialize c[i] mit a[i,i]
41 *
42 * d_ij=max(0,-a[i,j],-a[j,i])
43 * b[i,j]=-(a[i,j]+d_ij)
44 * c[i]-=d_ij
45 */
46
47 void Paso_FCTSolver_setLowOrderOperator(Paso_TransportProblem * fc) {
48 dim_t n=Paso_SystemMatrix_getTotalNumRows(fc->transport_matrix),i;
49 const index_t* main_iptr=NULL;
50 index_t iptr_ij,j,iptr_ji;
51 Paso_SystemMatrixPattern *pattern;
52 register double d_ij, sum, rtmp1, rtmp2;
53
54 if (fc==NULL) return;
55 main_iptr=Paso_TransportProblem_borrowMainDiagonalPointer(fc);
56 if (fc->iteration_matrix==NULL) {
57 fc->iteration_matrix=Paso_SystemMatrix_alloc(fc->transport_matrix->type,
58 fc->transport_matrix->pattern,
59 fc->transport_matrix->row_block_size,
60 fc->transport_matrix->col_block_size, TRUE);
61 }
62
63 if (Esys_noError()) {
64 pattern=fc->iteration_matrix->pattern;
65 n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
66 #pragma omp parallel for private(i,sum,iptr_ij,j,iptr_ji,rtmp1, rtmp2,d_ij) schedule(static)
67 for (i = 0; i < n; ++i) {
68 sum=fc->transport_matrix->mainBlock->val[main_iptr[i]];
69
70 /* printf("sum[%d] = %e -> ", i, sum); */
71 /* look at a[i,j] */
72 for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
73 j=pattern->mainPattern->index[iptr_ij];
74 rtmp1=fc->transport_matrix->mainBlock->val[iptr_ij];
75 if (j!=i) {
76 /* find entry a[j,i] */
77 #pragma ivdep
78 for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_ji<pattern->mainPattern->ptr[j+1]; ++iptr_ji) {
79
80 if ( pattern->mainPattern->index[iptr_ji] == i) {
81 rtmp2=fc->transport_matrix->mainBlock->val[iptr_ji];
82 /*
83 printf("a[%d,%d]=%e\n",i,j,rtmp1);
84 printf("a[%d,%d]=%e\n",j,i,rtmp2);
85 */
86
87 d_ij=-MIN3(0.,rtmp1,rtmp2);
88 fc->iteration_matrix->mainBlock->val[iptr_ij]=-(rtmp1+d_ij);
89 /* printf("l[%d,%d]=%e\n",i,j,fc->iteration_matrix->mainBlock->val[iptr_ij]); */
90 sum-=d_ij;
91 break;
92 }
93 }
94 }
95 }
96 for (iptr_ij=pattern->col_couplePattern->ptr[i];iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
97 j=pattern->col_couplePattern->index[iptr_ij];
98 rtmp1=fc->transport_matrix->col_coupleBlock->val[iptr_ij];
99 /* find entry a[j,i] */
100 #pragma ivdep
101 for (iptr_ji=pattern->row_couplePattern->ptr[j]; iptr_ji<pattern->row_couplePattern->ptr[j+1]; ++iptr_ji) {
102 if (pattern->row_couplePattern->index[iptr_ji]==i) {
103 rtmp2=fc->transport_matrix->row_coupleBlock->val[iptr_ji];
104 d_ij=-MIN3(0.,rtmp1,rtmp2);
105 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]=-(rtmp1+d_ij);
106 fc->iteration_matrix->row_coupleBlock->val[iptr_ji]=-(rtmp2+d_ij);
107 sum-=d_ij;
108 break;
109 }
110 }
111 }
112 /* set main diagonal entry */
113 fc->main_diagonal_low_order_transport_matrix[i]=sum;
114 /* printf("%e \n", sum); */
115 }
116
117 }
118 }
119 /*
120 * out_i=m_i u_i + a * \sum_{j <> i} l_{ij} (u_j-u_i)
121 *
122 */
123 void Paso_FCTSolver_setMuPaLu(double* out,
124 const double* M,
125 const Paso_Coupler* u_coupler,
126 const double a,
127 const Paso_SystemMatrix *L)
128 {
129 dim_t n, i, j;
130 Paso_SystemMatrixPattern *pattern;
131 const double *u=Paso_Coupler_borrowLocalData(u_coupler);
132 const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
133 register double sum, u_i, l_ij;
134 register index_t iptr_ij;
135 n=Paso_SystemMatrix_getTotalNumRows(L);
136
137 #pragma omp parallel for private(i) schedule(static)
138 for (i = 0; i < n; ++i) {
139 out[i]=M[i]*u[i];
140 }
141 if (ABS(a)>0) {
142 pattern=L->pattern;
143 #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)
144 for (i = 0; i < n; ++i) {
145 sum=0;
146 u_i=u[i];
147 #pragma ivdep
148 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
149 j=pattern->mainPattern->index[iptr_ij];
150 l_ij=L->mainBlock->val[iptr_ij];
151 sum+=l_ij*(u[j]-u_i);
152 }
153 #pragma ivdep
154 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
155 j=pattern->col_couplePattern->index[iptr_ij];
156 l_ij=L->col_coupleBlock->val[iptr_ij];
157 sum+=l_ij*(remote_u[j]-u_i);
158 }
159 out[i]+=a*sum;
160 }
161 }
162 }
163 /*
164 * calculate QP[i] max_{i in L->pattern[i]} (u[j]-u[i] )
165 * QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )
166 *
167 */
168 void Paso_FCTSolver_setQs(const Paso_Coupler* u_coupler,double* QN, double* QP, const Paso_SystemMatrix *L)
169 {
170 dim_t n, i, j;
171 Paso_SystemMatrixPattern *pattern;
172 const double *u=Paso_Coupler_borrowLocalData(u_coupler);
173 const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
174 register double u_i, u_min_i, u_max_i, u_j;
175 register index_t iptr_ij;
176 n=Paso_SystemMatrix_getTotalNumRows(L);
177 pattern=L->pattern;
178 #pragma omp parallel for schedule(static) private(i, u_i, u_min_i, u_max_i, j, u_j, iptr_ij)
179 for (i = 0; i < n; ++i) {
180 u_i=u[i];
181 u_min_i=u_i;
182 u_max_i=u_i;
183 #pragma ivdep
184 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
185 j=pattern->mainPattern->index[iptr_ij];
186 u_j=u[j];
187 u_min_i=MIN(u_min_i,u_j);
188 u_max_i=MAX(u_max_i,u_j);
189 }
190 #pragma ivdep
191 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
192 j=pattern->col_couplePattern->index[iptr_ij];
193 u_j=remote_u[j];
194 u_min_i=MIN(u_min_i,u_j);
195 u_max_i=MAX(u_max_i,u_j);
196 }
197 QN[i]=u_min_i-u_i;
198 QP[i]=u_max_i-u_i;
199 }
200 }
201
202 /*
203 *
204 * f_{ij} = (m_{ij} - dt (1-theta) d_{ij}) (u_last[j]-u_last[i]) - (m_{ij} + dt theta d_{ij}) (u[j]-u[i])
205 *
206 * m=fc->mass matrix
207 * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
208 */
209 void Paso_FCTSolver_setAntiDiffusionFlux(const double dt, const Paso_TransportProblem * fc, Paso_SystemMatrix *flux_matrix, const Paso_Coupler* u_coupler)
210 {
211 dim_t n, j, i;
212 index_t iptr_ij;
213 const double *u=Paso_Coupler_borrowLocalData(u_coupler);
214 const double *u_last= Paso_Coupler_borrowLocalData(fc->u_coupler);
215 const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
216 const double *remote_u_last=Paso_Coupler_borrowRemoteData(fc->u_coupler);
217 const double theta= (fc->useBackwardEuler) ? 1. : 0.5;
218 const double f1= - dt * ( 1.- theta );
219 const double f2= dt * theta;
220 register double m_ij, d_ij, u_i, u_last_i, d_u_last, d_u;
221 const Paso_SystemMatrixPattern *pattern=fc->iteration_matrix->pattern;
222 n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
223
224 if ( (ABS(f1) >0 ) ) {
225 if ( (ABS(f2) >0 ) ) {
226 #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
227 for (i = 0; i < n; ++i) {
228 u_i=u[i];
229 u_last_i=u_last[i];
230 #pragma ivdep
231 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
232 j=pattern->mainPattern->index[iptr_ij];
233 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
234 d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
235 fc->iteration_matrix->mainBlock->val[iptr_ij]);
236 d_u=u[j]-u_i;
237 d_u_last=u_last[j]-u_last_i;
238
239 /* (m_{ij} - dt (1-theta) d_{ij}) (u_last[j]-u_last[i]) - (m_{ij} + dt theta d_{ij}) (u[j]-u[i]) */
240
241 flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
242 }
243 #pragma ivdep
244 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
245 j=pattern->col_couplePattern->index[iptr_ij];
246 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
247 d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
248 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
249 d_u=remote_u[j]-u_i;
250 d_u_last=remote_u_last[j]-u_last_i;
251 flux_matrix->col_coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
252 }
253 }
254 } else {
255 #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
256 for (i = 0; i < n; ++i) {
257 u_i=u[i];
258 u_last_i=u_last[i];
259 #pragma ivdep
260 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
261 j=pattern->mainPattern->index[iptr_ij];
262 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
263 d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
264 fc->iteration_matrix->mainBlock->val[iptr_ij]);
265 d_u=u[j]-u_i;
266 d_u_last=u_last[j]-u_last_i;
267 flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
268 }
269 #pragma ivdep
270 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
271 j=pattern->col_couplePattern->index[iptr_ij];
272 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
273 d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
274 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
275 d_u=remote_u[j]-u_i;
276 d_u_last=remote_u_last[j]-u_last_i;
277 flux_matrix->col_coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
278 }
279 }
280 }
281 } else {
282 if ( (ABS(f2) >0 ) ) {
283 #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
284 for (i = 0; i < n; ++i) {
285 u_i=u[i];
286 u_last_i=u_last[i];
287 #pragma ivdep
288 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
289 j=pattern->mainPattern->index[iptr_ij];
290 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
291 d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
292 fc->iteration_matrix->mainBlock->val[iptr_ij]);
293 d_u=u[j]-u_i;
294 d_u_last=u_last[j]-u_last_i;
295 flux_matrix->mainBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
296 }
297 #pragma ivdep
298 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
299 j=pattern->col_couplePattern->index[iptr_ij];
300 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
301 d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
302 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
303 d_u=remote_u[j]-u_i;
304 d_u_last=remote_u_last[j]-u_last_i;
305 flux_matrix->col_coupleBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
306 }
307 }
308 } else {
309 #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
310 for (i = 0; i < n; ++i) {
311 u_i=u[i];
312 u_last_i=u_last[i];
313 #pragma ivdep
314 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
315 j=pattern->mainPattern->index[iptr_ij];
316 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
317 d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
318 fc->iteration_matrix->mainBlock->val[iptr_ij]);
319 d_u=u[j]-u_i;
320 d_u_last=u_last[j]-u_last_i;
321 flux_matrix->mainBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
322 }
323 #pragma ivdep
324 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
325 j=pattern->col_couplePattern->index[iptr_ij];
326 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
327 d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
328 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
329 d_u=remote_u[j]-u_i;
330 d_u_last=remote_u_last[j]-u_last_i;
331 flux_matrix->col_coupleBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
332 }
333 }
334 }
335 }
336 }
337 /*
338 * apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0
339 *
340 */
341 void Paso_FCTSolver_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const Paso_Coupler* u_coupler)
342 {
343 dim_t n, i, j;
344 Paso_SystemMatrixPattern *pattern;
345 const double *u=Paso_Coupler_borrowLocalData(u_coupler);
346 const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
347 register double u_i, f_ij;
348 register index_t iptr_ij;
349
350 n=Paso_SystemMatrix_getTotalNumRows(f);
351 pattern=f->pattern;
352 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)
353 for (i = 0; i < n; ++i) {
354 u_i=u[i];
355 #pragma ivdep
356 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
357 j=pattern->mainPattern->index[iptr_ij];
358 f_ij=f->mainBlock->val[iptr_ij];
359 if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0;
360 }
361 #pragma ivdep
362 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
363 j=pattern->col_couplePattern->index[iptr_ij];
364 f_ij=f->col_coupleBlock->val[iptr_ij];
365 if (f_ij * (u_i-remote_u[j]) <= 0) f->col_coupleBlock->val[iptr_ij]=0;
366 }
367 }
368 }
369
370
371 void Paso_FCTSolver_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
372 const Paso_Coupler* QN_coupler,const Paso_Coupler* QP_coupler,double* RN,double* RP)
373 {
374 dim_t n, i, j;
375 Paso_SystemMatrixPattern *pattern;
376 const double *QN=Paso_Coupler_borrowLocalData(QN_coupler);
377 const double *QP=Paso_Coupler_borrowLocalData(QP_coupler);
378 register double f_ij, PP_i, PN_i;
379 register index_t iptr_ij;
380
381 n=Paso_SystemMatrix_getTotalNumRows(f);
382 pattern=f->pattern;
383 #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
384 for (i = 0; i < n; ++i) {
385 PP_i=0.;
386 PN_i=0.;
387 #pragma ivdep
388 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
389 j=pattern->mainPattern->index[iptr_ij];
390 if (i != j ) {
391 f_ij=f->mainBlock->val[iptr_ij];
392 if (f_ij <=0) {
393 PN_i+=f_ij;
394 } else {
395 PP_i+=f_ij;
396 }
397 }
398 }
399 #pragma ivdep
400 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
401 j=pattern->col_couplePattern->index[iptr_ij];
402 f_ij=f->col_coupleBlock->val[iptr_ij];
403 if (f_ij <=0 ) {
404 PN_i+=f_ij;
405 } else {
406 PP_i+=f_ij;
407 }
408 }
409 if (PN_i<0) {
410 RN[i]=MIN(1,QN[i]*lumped_mass_matrix[i]/PN_i);
411 } else {
412 RN[i]=0.;
413 }
414 if (PP_i>0) {
415 RP[i]=MIN(1,QP[i]*lumped_mass_matrix[i]/PP_i);
416 } else {
417 RP[i]=0.;
418 }
419 }
420
421 }
422
423 void Paso_FCTSolver_addCorrectedFluxes(double* f,const Paso_SystemMatrix *flux_matrix, const Paso_Coupler* RN_coupler, const Paso_Coupler* RP_coupler)
424 {
425 dim_t i, j;
426 Paso_SystemMatrixPattern *pattern;
427 register double RN_i, RP_i, f_i, f_ij, rtmp;
428 register index_t iptr_ij;
429 const double *RN=Paso_Coupler_borrowLocalData(RN_coupler);
430 const double *remote_RN=Paso_Coupler_borrowRemoteData(RN_coupler);
431 const double *RP=Paso_Coupler_borrowLocalData(RP_coupler);
432 const double *remote_RP=Paso_Coupler_borrowRemoteData(RP_coupler);
433 const dim_t n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
434
435 pattern=flux_matrix->pattern;
436 #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i, rtmp)
437 for (i = 0; i < n; ++i) {
438
439 RN_i=RN[i];
440 RP_i=RP[i];
441 f_i=0;
442 #pragma ivdep
443 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
444 j=pattern->mainPattern->index[iptr_ij];
445 f_ij=flux_matrix->mainBlock->val[iptr_ij];
446 rtmp=(f_ij >=0 ) ? MIN(RP_i,RN[j]) : MIN(RN_i,RP[j]) ;
447 f_i+=f_ij*rtmp;
448 }
449 #pragma ivdep
450 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
451 j=pattern->col_couplePattern->index[iptr_ij];
452 f_ij=flux_matrix->col_coupleBlock->val[iptr_ij];
453 rtmp=(f_ij >=0 ) ? MIN(RP_i,remote_RN[j]) : MIN(RN_i,remote_RP[j]) ;
454 f_i+=f_ij*rtmp;
455 }
456 f[i]+=f_i;
457
458 }
459 }
460

  ViewVC Help
Powered by ViewVC 1.1.26