TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
linearalgebra.c
Go to the documentation of this file.
1 #include "turtle.h"
2 #include "linearalgebra.h"
3 
4 /*-----------------------------------------------------*/
5 
6 void convlv(FLOATVECTOR *data,unsigned long n,FLOATVECTOR *respns,unsigned long m,int isign,
7  double delta_time,FLOATVECTOR *ans)
8 
9 {
10  unsigned long i,no2;
11  float dum,mag2;
12  FLOATVECTOR *fft1;
13  FLOATVECTOR *fft;
14 
15  fft1=new_floatvector(n<<1);
16 
17  fft=fft1;
18  for (i=1;i<=(m-1)/2;i++){
19  respns->co[n+1-i]=respns->co[m+1-i];
20  }
21  for (i=(m+3)/2;i<=n-(m-1)/2;i++){
22  respns->co[i]=0.0;
23  }
24  twofft(data,respns,fft,ans,n);
25  no2=n>>1;
26  for (i=2;i<=n+2;i+=2) {
27  if (isign == 1) {
28  ans->co[i-1]=(fft->co[i-1]*(dum=ans->co[i-1])-fft->co[i]*ans->co[i])/no2*delta_time;
29  ans->co[i]=(fft->co[i]*dum+fft->co[i-1]*ans->co[i])/no2*delta_time;
30  } else if (isign == -1) {
31  if ((mag2=SQR(ans->co[i-1])+SQR(ans->co[i])) == 0.0)
32  t_error("Deconvolving at response zero in convlv");
33  ans->co[i-1]=(fft->co[i-1]*(dum=ans->co[i-1])+fft->co[i]*ans->co[i])/mag2/no2;
34  ans->co[i]=(fft->co[i]*dum-fft->co[i-1]*ans->co[i])/mag2/no2;
35  } else t_error("No meaning for isign in convlv");
36  }
37  ans->co[2]=ans->co[n+1];
38  realft(ans,n,-1);
39 
40  for(i=1;i<=n;i++){
41  ans->co[i]=ans->co[i];
42  }
43  free_floatvector(fft1);
44 }
45 
46 /*-----------------------------------------------------*/
47 
48 void four1(FLOATVECTOR *data,unsigned long nn,int isign)
49 
50 {
51  unsigned long n,mmax,m,j,istep,i;
52  float wtemp,wr,wpr,wpi,wi,theta;
53  float tempr,tempi;
54 
55  n=nn << 1;
56  j=1;
57  for (i=1;i<n;i+=2) {
58  if (j > i) {
59  SWAP(data->co[j],data->co[i]);
60  SWAP(data->co[j+1],data->co[i+1]);
61  }
62  m=n >> 1;
63  while (m >= 2 && j > m) {
64  j -= m;
65  m >>= 1;
66  }
67  j += m;
68  }
69  mmax=2;
70  while (n > mmax) {
71  istep=mmax << 1;
72  theta=isign*(6.28318530717959/mmax);
73  wtemp=sin(0.5*theta);
74  wpr = -2.0*wtemp*wtemp;
75  wpi=sin(theta);
76  wr=1.0;
77  wi=0.0;
78  for (m=1;m<mmax;m+=2) {
79  for (i=m;i<=n;i+=istep) {
80  j=i+mmax;
81  tempr=wr*data->co[j]-wi*data->co[j+1];
82  tempi=wr*data->co[j+1]+wi*data->co[j];
83  data->co[j]=data->co[i]-tempr;
84  data->co[j+1]=data->co[i+1]-tempi;
85  data->co[i] += tempr;
86  data->co[i+1] += tempi;
87  }
88  wr=(wtemp=wr)*wpr-wi*wpi+wr;
89  wi=wi*wpr+wtemp*wpi+wi;
90  }
91  mmax=istep;
92  }
93 }
94 
95 /*-----------------------------------------------------*/
96 
97 void ludcmp(SHORTVECTOR *indx, DOUBLEMATRIX *var)
98 
99 {
100 short mode;
101 long i,j,imax,k,n;
102 double aamax,dum,sum;
103 DOUBLEVECTOR *vv;
104 mode=1;
105 n=var->nrh;
106 vv=new_doublevector(n);
107 /*calcolo l'elemento di modulo massimo nella riga i e divido tutti gli elementi
108  per quel numero*/
109 for(i=1;i<=n;i++){
110  aamax=0;
111  for(j=1;j<=n;j++){
112  if(fabs(var->co[i][j])>aamax)aamax=fabs(var->co[i][j]);
113  }
114  if(aamax==0){
115  printf("aamax e' nullo nella riga (%ld)\n",i);
116  aamax=1.0;
117  }
118  vv->co[i]=1.0/aamax;
119 }
120 /*Questo ciclo va eseguito per tutte le colonne della matrice*/
121 for(j=1;j<=n;j++){
122 /*Questo ciclo cambia il valore della matrice modificando secondo*/
123  for(i=1;i<=j-1;i++){
124  sum=var->co[i][j];
125  for(k=1;k<=i-1;k++){
126  sum-=var->co[i][k]*var->co[k][j];
127  }
128  var->co[i][j]=sum;
129  }
130  aamax=0.0;
131  for(i=j;i<=n;i++){
132  sum=var->co[i][j];
133  for(k=1;k<=j-1;k++){
134  sum-=var->co[i][k]*var->co[k][j];
135  }
136  var->co[i][j]=sum;
137  dum=vv->co[i]*fabs(sum);
138  if(dum>aamax){
139  imax=i;
140  aamax=dum;
141  }
142  }
143  if(j!=imax){
144  for(i=1;i<=n;i++){
145  dum=var->co[imax][i];
146  var->co[imax][i]=var->co[j][i];
147  var->co[j][i]=dum;
148  }
149  vv->co[imax]=vv->co[j];
150  }
151  indx->co[j]=imax;
152  if(var->co[j][j]==0){
153  var->co[j][j]=TINY;
154  printf("var e' nullo\n");
155  }
156 
157  if(j!=n){
158  dum=1.0/var->co[j][j];
159  for(i=j+1;i<=n;i++){
160  var->co[i][j]*=dum;
161  }
162  }
163 }
164 free_doublevector(vv);
165 }
166 
167 /*-----------------------------------------------------*/
168 
169 void lubksb(DOUBLEMATRIX *var, SHORTVECTOR *indx,DOUBLEVECTOR *gam)
170 
171 {
172 long i,j,k,ii,n;
173 double sum;
174 ii=0;
175 n=var->nrh;
176 for(i=1;i<=n;i++){
177  k=indx->co[i];
178  sum=gam->co[k];
179  gam->co[k]=gam->co[i];
180  if(ii!=0){
181 /*calcolo il valore del termine modificato*/
182  for(j=ii;j<=i-1;j++){
183  sum-=var->co[i][j]*gam->co[j];
184  }
185  }else if(sum!=0){
186  ii=i;
187  }
188  gam->co[i]=sum;
189 }
190 /*Valuto il vettore incognito eseguendo una back substitution*/
191 for(i=n;i>=1;i--){
192  sum=gam->co[i];
193  for(j=i+1;j<=n;j++){
194  sum-=var->co[i][j]*gam->co[j];
195  }
196  gam->co[i]=sum/var->co[i][i];
197 }
198 }
199 
200 /*------------------------------------------*/
201 
202 void realft(FLOATVECTOR *data,unsigned long n,int isign)
203 {
204  unsigned long i,i1,i2,i3,i4,np3;
205  float c1=0.5,c2,h1r,h1i,h2r,h2i;
206  float wr,wi,wpr,wpi,wtemp,theta;
207 
208  theta=3.141592653589793/(float) (n>>1);
209  if (isign == 1) {
210  c2 = -0.5;
211  four1(data,n>>1,1);
212  } else {
213  c2=0.5;
214  theta = -theta;
215  }
216  wtemp=sin(0.5*theta);
217  wpr = -2.0*wtemp*wtemp;
218  wpi=sin(theta);
219  wr=1.0+wpr;
220  wi=wpi;
221  np3=n+3;
222  for (i=2;i<=(n>>2);i++) {
223  i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
224  h1r=c1*(data->co[i1]+data->co[i3]);
225  h1i=c1*(data->co[i2]-data->co[i4]);
226  h2r = -c2*(data->co[i2]+data->co[i4]);
227  h2i=c2*(data->co[i1]-data->co[i3]);
228  data->co[i1]=h1r+wr*h2r-wi*h2i;
229  data->co[i2]=h1i+wr*h2i+wi*h2r;
230  data->co[i3]=h1r-wr*h2r+wi*h2i;
231  data->co[i4] = -h1i+wr*h2i+wi*h2r;
232  wr=(wtemp=wr)*wpr-wi*wpi+wr;
233  wi=wi*wpr+wtemp*wpi+wi;
234  }
235  if (isign == 1) {
236  data->co[1] = (h1r=data->co[1])+data->co[2];
237  data->co[2] = h1r-data->co[2];
238  } else {
239  data->co[1]=c1*((h1r=data->co[1])+data->co[2]);
240  data->co[2]=c1*(h1r-data->co[2]);
241  four1(data,n>>1,-1);
242  }
243 }
244 
245 /*-----------------------------------------------------*/
246 
247 void twofft(FLOATVECTOR *data1,FLOATVECTOR *data2,FLOATVECTOR *fft1,FLOATVECTOR *fft2,unsigned long n)
248 
249 {
250 
251  unsigned long nn3,nn2,jj,j;
252 
253  float rep,rem,aip,aim;
254 
255 
256 
257  nn3=1+(nn2=2+n+n);
258 
259  for (j=1,jj=2;j<=n;j++,jj+=2) {
260 
261  fft1->co[jj-1]=data1->co[j];
262 
263  fft1->co[jj]=data2->co[j];
264 
265  }
266 
267  four1(fft1,n,1);
268 
269  fft2->co[1]=fft1->co[2];
270 
271  fft1->co[2]=fft2->co[2]=0.0;
272 
273  for (j=3;j<=n+1;j+=2) {
274 
275  rep=0.5*(fft1->co[j]+fft1->co[nn2-j]);
276 
277  rem=0.5*(fft1->co[j]-fft1->co[nn2-j]);
278 
279  aip=0.5*(fft1->co[j+1]+fft1->co[nn3-j]);
280 
281  aim=0.5*(fft1->co[j+1]-fft1->co[nn3-j]);
282 
283  fft1->co[j]=rep;
284 
285  fft1->co[j+1]=aim;
286 
287  fft1->co[nn2-j]=rep;
288 
289  fft1->co[nn3-j] = -aim;
290 
291  fft2->co[j]=aip;
292 
293  fft2->co[j+1] = -rem;
294 
295  fft2->co[nn2-j]=aip;
296 
297  fft2->co[nn3-j]=rem;
298 
299  }
300 
301 }
302 
303 
304 /*-----------------------------------------------------*/
305 /* commentata per l'errore (dovuto ad una errata definizione di asolve):
306 error: conflicting types for `asolve'
307 void ris_sistema (double d[], double ds[], double di[], double b[], double x[], int n)
308 {
309 
310 int iter;
311 double err;
312 
313 LONGVECTOR *ija;
314 
315 DOUBLEVECTOR *sa;
316 
317 DOUBLEMATRIX *A;*/
318 
319 /* Alloco una matrice normale a partire da tre vettori, uno per diagonale */
320 /* commentata per l'errore (dovuto ad una errata definizione di asolve):
321 error: conflicting types for `asolve'
322 A=vett_mat(d,ds,di,n);*/
323 
324 /* Alloco i vettori che memorizzano la matrice tridiagonale
325  come una matrice sparsa alla maniera di N.R.*/
326 /* commentata per l'errore (dovuto ad una errata definizione di asolve):
327 error: conflicting types for `asolve'
328 sa=new_doublevector((3*n)-1);
329 ija=new_longvector((3*n)-1);
330 */
331 /* Definisco gli elementi dei vettori sa[] e ija[] */
332 /* commentata per l'errore (dovuto ad una errata definizione di asolve):
333 error: conflicting types for `asolve'
334 sprsin(A->co,A->nrh,THRESH,(3*A->nrh)-1,sa->co,ija->co); */
335 
336 /* Calcolo la soluzione del sistema lineare */
337 /* commentata per l'errore (dovuto ad una errata definizione di asolve):
338 error: conflicting types for `asolve'
339 linbcg(n,b,x,ITOL,TOL,ITMAX,&iter,&err,sa->co,ija->co);*/
340 
341 /* Deallocazioni */
342 /* commentata per l'errore (dovuto ad una errata definizione di asolve):
343 error: conflicting types for `asolve'
344 free_longvector(ija);
345 free_doublevector(sa);
346 free_doublematrix(A);
347 
348 }
349 
350 */
351 
352 /*
353 ___________________________________________________________________________
354 ===========================================================================
355  FUNZIONE vett_mat
356 ___________________________________________________________________________
357 */
358 /* Questa funzione converte tre vettori in una matrice quadrata tridiagonale.
359  Bisogna passare alla funzione i puntatori agli elementi dei vettori
360  - d elementi delle diagonale principale;
361  - ds elementi della diagonale superiore;
362  - di elementi della diagonale inferiore;
363  inoltre bisogna passare n, che e' la dimensione della matrice che
364  si vuole generare.
365 
366  La funzione restituisce il puntatore alla matrice.
367  */
368 
369 DOUBLEMATRIX *vett_mat (double *d,double *ds,double *di,int n)
370 {
371  int i,j;
372 
373  DOUBLEMATRIX *A;
374 
375  A=new_doublematrix(n,n);
376 
377  for (i=1; i<=n; i++)
378  { for (j=1; j<=n; j++) A->co[i][j]=0; }
379 
380  A->co[1][2]=ds[1];
381  A->co[n][n-1]=di[n-1];
382 
383  for (i=1; i<=n; i++) A->co[i][i]=d[i];
384  for (i=2; i<=n-1; i++)
385  {
386  A->co[i][i+1]=ds[i];
387  A->co[i][i-1]=di[i-1];
388  }
389 
390  return A;
391 }
392 
393 /*DA NUMERICAL RECEPIS IN C. (Second Edition - Cambridge Univ. Press). pag 79
394 ______________________________________________________________________________
395 ==============================================================================
396  FUNZIONE sprsin
397 ______________________________________________________________________________
398 */
399 /*
400  Questa funzione converte una matrice memorizzata nel modo convenzionale
401  in un vettore sa[] che contiene solo i valori non nulli della matrice
402  e in un vettore ija[] che permette di individuare la posizione originale
403  degli elementi di sa[].
404 
405  La funzione richiede:
406  - **a, un puntatore agli elementi della matrice originale;
407  - n, dimensione della matrice;
408  - thresh, gli elementi della matrice minori di thresh non vengono
409  letti;
410  - nmax, la lunghezza dei vettori sa[] e ija[].
411 */
412 
413 
414 void sprsin(double **a, int n, float thresh, long nmax, double sa[],
415  long ija[])
416 {
417  int i,j;
418  long k;
419 
420  for (j=1;j<=n;j++) sa[j]=a[j][j];
421  ija[1]=n+2;
422  k=n+1;
423  for (i=1;i<=n;i++) {
424  for (j=1;j<=n;j++) {
425  if (fabs(a[i][j]) > thresh && i != j) {
426  if (++k > nmax) t_error("sprsin: nmax too small");
427  sa[k]=a[i][j];
428  ija[k]=j;
429  }
430  }
431  ija[i+1]=k+1;
432  }
433 }
434 
435 /* DA NUMERICAL RECEPIS IN C. (Second Edition - Cambridge Univ. Press). pag 86-88
436 _________________________________________________________________________________
437 =================================================================================
438  FUNZIONE linbcg
439 _________________________________________________________________________________
440 */
441 /* Questa funzione consente di risolvere di risolvere un sistema lineare del tipo
442  A x = b con il metodo iterativo del gradiente coniugato.
443 
444  Alla funzione devono essere passati i seguenti argomenti:
445  - n: dimensione del sistema;
446  - sa[] e ija[]: vettori generati dalla funzione sprssin() che memorizzano la
447  matrice;
448  - b[]: elementi del vettore dei termini noti;
449  - x[]: elementi del vettore soluzione (in ingresso questo vettore deve contenere
450  una soluzione di primo tentativo);
451  - itol, tol, itmax: parametri gli definiti sopra.
452 
453  Oltre alla soluzione la funzione calcola anche il numero di iterazione effetuate
454  ( iter ) e l'errore commesso ( err ).
455 */
456 
457 
458 /* commentata per l'errore (dovuto ad una errata definizione di asolve):
459 error: conflicting types for `asolve'
460 
461 void linbcg(long n, double b[], double x[], int itol, double tol,
462  int itmax, int *iter, double *err, double sa[], long ija[])
463 {
464  void asolve(long n, double b[], double x[], int itrnsp,double sa[], long ija[]);
465  void atimes(long n, double x[], double r[], int itrnsp,double sa[], long ija[]);
466  double snrm(long n, double sx[], int itol);
467  long j;
468  double ak,akden,bk,bkden,bknum,bnrm,dxnrm,xnrm,zm1nrm,znrm;
469 
470  DOUBLEVECTOR *P,*PP,*R,*RR,*Z,*ZZ;
471  double *p,*pp,*r,*rr,*z,*zz;
472 
473 
474  P=new_doublevector(n);
475  PP=new_doublevector(n);
476  R=new_doublevector(n);
477  RR=new_doublevector(n);
478  Z=new_doublevector(n);
479  ZZ=new_doublevector(n);
480 
481  p=P->co;
482  pp=PP->co;
483  r=R->co;
484  rr=RR->co;
485  z=Z->co;
486  zz=ZZ->co;
487 
488  *iter=0;
489  atimes(n,x,r,0,sa,ija);
490  for (j=1;j<=n;j++) {
491  r[j]=b[j]-r[j];
492  rr[j]=r[j];
493  }
494  if (itol == 1) {
495  bnrm=snrm(n,b,itol);
496  asolve(n,r,z,0,sa,ija);
497  }
498  else if (itol == 2) {
499  asolve(n,b,z,0,sa,ija);
500  bnrm=snrm(n,z,itol);
501  asolve(n,r,z,0,sa,ija);
502  }
503  else if (itol == 3 || itol == 4) {
504  asolve(n,b,z,0,sa,ija);
505  bnrm=snrm(n,z,itol);
506  asolve(n,r,z,0,sa,ija);
507  znrm=snrm(n,z,itol);
508  } else t_error("illegal itol in linbcg");
509  while (*iter <= itmax) {
510  ++(*iter);
511  asolve(n,rr,zz,1,sa,ija);
512  for (bknum=0.0,j=1;j<=n;j++) bknum += z[j]*rr[j];
513  if (*iter == 1) {
514  for (j=1;j<=n;j++) {
515  p[j]=z[j];
516  pp[j]=zz[j];
517  }
518  }
519  else {
520  bk=bknum/bkden;
521  for (j=1;j<=n;j++) {
522  p[j]=bk*p[j]+z[j];
523  pp[j]=bk*pp[j]+zz[j];
524  }
525  }
526  bkden=bknum;
527  atimes(n,p,z,0,sa,ija);
528  for (akden=0.0,j=1;j<=n;j++) akden += z[j]*pp[j];
529  ak=bknum/akden;
530  atimes(n,pp,zz,1,sa,ija);
531  for (j=1;j<=n;j++) {
532  x[j] += ak*p[j];
533  r[j] -= ak*z[j];
534  rr[j] -= ak*zz[j];
535  }
536  asolve(n,r,z,0,sa,ija);
537  if (itol == 1)
538  *err=snrm(n,r,itol)/bnrm;
539  else if (itol == 2)
540  *err=snrm(n,z,itol)/bnrm;
541  else if (itol == 3 || itol == 4) {
542  zm1nrm=znrm;
543  znrm=snrm(n,z,itol);
544  if (fabs(zm1nrm-znrm) > EPS*znrm) {
545  dxnrm=fabs(ak)*snrm(n,p,itol);
546  *err=znrm/fabs(zm1nrm-znrm)*dxnrm;
547  } else {
548  *err=znrm/bnrm;
549  continue;
550  }
551  xnrm=snrm(n,x,itol);
552  if (*err <= 0.5*xnrm) *err /= xnrm;
553  else {
554  *err=znrm/bnrm;
555  continue;
556  }
557  }
558  if (*err <= tol) break;
559  }
560 
561  free_doublevector(P);
562  free_doublevector(PP);
563  free_doublevector(R);
564  free_doublevector(RR);
565  free_doublevector(Z);
566  free_doublevector(ZZ);
567 }
568 */
569 /* DA NUMERICAL RECEPIS IN C. (Second Edition - Cambridge Univ. Press). pag 88
570 _________________________________________________________________________________
571 =================================================================================
572  FUNZIONE snrm
573 _________________________________________________________________________________
574 */
575 /* Questa funzione calcola la norma di un vettore con la modalita' specificata
576  dal parametro itol */
577 
578 
579 double snrm(long n, double sx[], int itol)
580 {
581  long i,isamax;
582  double ans;
583 
584  if (itol <= 3) {
585  ans = 0.0;
586  for (i=1;i<=n;i++) ans += sx[i]*sx[i];
587  return sqrt(ans);
588  } else {
589  isamax=1;
590  for (i=1;i<=n;i++) {
591  if (fabs(sx[i]) > fabs(sx[isamax])) isamax=i;
592  }
593  return fabs(sx[isamax]);
594  }
595 }
596 
597 
598 
599 
600 
601 
602 
603 
604 
605 
606 /* DA NUMERICAL RECEPIS IN C. (Second Edition - Cambridge Univ. Press). pag 88
607 _________________________________________________________________________________
608 =================================================================================
609  FUNZIONE atimes
610 _________________________________________________________________________________
611 */
612 
613 
614 
615 void atimes(long n, double x[], double r[], int itrnsp,double sa[], long ija[])
616 {
617  void dsprsax(double sa[], long ija[], double x[], double b[],
618  long n);
619  void dsprstx(double sa[], long ija[], double x[], double b[],
620  long n);
621 
622  if (itrnsp) dsprstx(sa,ija,x,r,n);
623  else dsprsax(sa,ija,x,r,n);
624 }
625 
626 
627 /*DA NUMERICAL RECEPIS IN C. (Second Edition - Cambridge Univ. Press). pag 89
628 _________________________________________________________________________________
629 =================================================================================
630  FUNZIONE asolve
631 _________________________________________________________________________________
632 */
633 
634 
635 void asolve(long n, double b[], double x[],double sa[])
636 {
637  long i;
638 
639  for(i=1;i<=n;i++) x[i]=(sa[i] != 0.0 ? b[i]/sa[i] : b[i]);
640 }
641 
642 
643 /* DA NUMERICAL RECEPIS IN C. (Second Edition - Cambridge Univ. Press). pag 79
644 _________________________________________________________________________________
645 =================================================================================
646  FUNZIONE dsprsax
647 _________________________________________________________________________________
648 */
649 /* Questa funzione moltiplica una matrice, memoririzzata alla maniera di N.R.,
650  per un vettore x[]. Il risultato e' un vettore b[].
651 */
652 
653 
654 void dsprsax(double sa[], long ija[], double x[], double b[], long n)
655 {
656  long i,k;
657 
658  if (ija[1] != n+2) t_error("dsprsax: mismatched vector and matrix");
659  for (i=1;i<=n;i++) {
660  b[i]=sa[i]*x[i];
661  for (k=ija[i];k<=ija[i+1]-1;k++) b[i] += sa[k]*x[ija[k]];
662  }
663 }
664 
665 
666 /* DA NUMERICAL RECEPIS IN C. (Second Edition - Cambridge Univ. Press). pag 80
667 _________________________________________________________________________________
668 =================================================================================
669  FUNZIONE dsprstx
670 _________________________________________________________________________________
671 */
672 /* Questa funzione moltiplica la trasposta di una matrice, memoririzzata alla
673  maniera di N.R., per un vettore x[]. Il risultato e' un vettore b[].
674 */
675 
676 
677 void dsprstx(double sa[], long ija[], double x[], double b[], long n)
678 {
679  long i,j,k;
680  if (ija[1] != n+2) t_error("mismatched vector and matrix in dsprstx");
681  for (i=1;i<=n;i++) b[i]=sa[i]*x[i];
682  for (i=1;i<=n;i++) {
683  for (k=ija[i];k<=ija[i+1]-1;k++) {
684  j=ija[k];
685  b[j] += sa[k]*x[i];
686  }
687  }
688 }
689 /*
690 _________________________________________________________________________________
691 =================================================================================
692  FUNZIONE integral
693 _________________________________________________________________________________
694  Questa funzione calcola l'integrale come sommatoria della funz func
695 
696  */
697  float integration(float (*func)(float), float a, float b, int n)
698 {
699  float x=0.0,tnm=0.0,sum=0.0,del=0.0;
700  static float s=0.0;
701  int it=0,j=0;
702  if (n == 1) {
703  return (s=0.5*(b-a)*(func(a)+func(b)));
704  } else {
705  for (it=1,j=1;j<n-1;j++) {
706  it <<= 1;
707 
708  }
709 
710  tnm=it;
711  del=(b-a)/tnm;
712  x=a+0.5*del;
713  for (sum=0.0,s=0.0,j=1;j<=it;j++,x+=del) {
714  sum += func(x);
715 
716  }
717  s=(s+(b-a)*sum/tnm);
718  return s;
719  }
720 }