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

  ViewVC Help
Powered by ViewVC 1.1.26