src/gromacs/gmxana/levenmar.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <math.h>
  42 #include <stdio.h>
  43 #include <stdlib.h>
  44
  45 #include "types/simple.h"
  46
  47 static void nrerror(const char error_text[], gmx_bool bExit)
  48 {
  49     fprintf(stderr, "Numerical Recipes run-time error...\n");
  50     fprintf(stderr, "%s\n", error_text);
  51     if (bExit)
  52     {
  53         fprintf(stderr, "...now exiting to system...\n");
  54         exit(1);
  55     }
  56 }
  57
  58 /* dont use the keyword vector - it will clash with the
  59  * altivec extensions used for powerpc processors.
  60  */
  61
  62 static real *rvector(int nl, int nh)
  63 {
  64     real *v;
  65
  66     v = (real *)malloc((unsigned) (nh-nl+1)*sizeof(real));
  67     if (!v)
  68     {
  69         nrerror("allocation failure in rvector()", TRUE);
  70     }
  71     /* cppcheck-suppress memleak
  72      * free_vector does the same vector arithmetic */
  73     return v-nl;
  74 }
  75
  76 static int *ivector(int nl, int nh)
  77 {
  78     int *v;
  79
  80     v = (int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
  81     if (!v)
  82     {
  83         nrerror("allocation failure in ivector()", TRUE);
  84     }
  85     /* cppcheck-suppress memleak
  86      * free_vector does the same vector arithmetic */
  87     return v-nl;
  88 }
  89
  90
  91 static real **matrix1(int nrl, int nrh, int ncl, int nch)
  92 {
  93     int    i;
  94     real **m;
  95
  96     m = (real **) malloc((unsigned) (nrh-nrl+1)*sizeof(real*));
  97     if (!m)
  98     {
  99         nrerror("allocation failure 1 in matrix1()", TRUE);
 100     }
 101     m -= nrl;
 102
 103     for (i = nrl; i <= nrh; i++)
 104     {
 105         m[i] = (real *) malloc((unsigned) (nch-ncl+1)*sizeof(real));
 106         if (!m[i])
 107         {
 108             nrerror("allocation failure 2 in matrix1()", TRUE);
 109         }
 110         m[i] -= ncl;
 111     }
 112     return m;
 113 }
 114
 115 static double **dmatrix(int nrl, int nrh, int ncl, int nch)
 116 {
 117     int      i;
 118     double **m;
 119
 120     m = (double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
 121     if (!m)
 122     {
 123         nrerror("allocation failure 1 in dmatrix()", TRUE);
 124     }
 125     m -= nrl;
 126
 127     for (i = nrl; i <= nrh; i++)
 128     {
 129         m[i] = (double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
 130         if (!m[i])
 131         {
 132             nrerror("allocation failure 2 in dmatrix()", TRUE);
 133         }
 134         m[i] -= ncl;
 135     }
 136     return m;
 137 }
 138
 139 static int **imatrix1(int nrl, int nrh, int ncl, int nch)
 140 {
 141     int i, **m;
 142
 143     m = (int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
 144     if (!m)
 145     {
 146         nrerror("allocation failure 1 in imatrix1()", TRUE);
 147     }
 148     m -= nrl;
 149
 150     for (i = nrl; i <= nrh; i++)
 151     {
 152         m[i] = (int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
 153         if (!m[i])
 154         {
 155             nrerror("allocation failure 2 in imatrix1()", TRUE);
 156         }
 157         m[i] -= ncl;
 158     }
 159     return m;
 160 }
 161
 162
 163
 164 static real **submatrix(real **a, int oldrl, int oldrh, int oldcl,
 165                         int newrl, int newcl)
 166 {
 167     int    i, j;
 168     real **m;
 169
 170     m = (real **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(real*));
 171     if (!m)
 172     {
 173         nrerror("allocation failure in submatrix()", TRUE);
 174     }
 175     m -= newrl;
 176
 177     for (i = oldrl, j = newrl; i <= oldrh; i++, j++)
 178     {
 179         m[j] = a[i]+oldcl-newcl;
 180     }
 181
 182     return m;
 183 }
 184
 185
 186
 187 static void free_vector(real *v, int nl)
 188 {
 189     free((char*) (v+nl));
 190 }
 191
 192 static void free_ivector(int *v, int nl)
 193 {
 194     free((char*) (v+nl));
 195 }
 196
 197 static void free_dvector(int *v, int nl)
 198 {
 199     free((char*) (v+nl));
 200 }
 201
 202
 203
 204 static void free_matrix(real **m, int nrl, int nrh, int ncl)
 205 {
 206     int i;
 207
 208     for (i = nrh; i >= nrl; i--)
 209     {
 210         free((char*) (m[i]+ncl));
 211     }
 212     free((char*) (m+nrl));
 213 }
 214
 215 static real **convert_matrix(real *a, int nrl, int nrh, int ncl, int nch)
 216 {
 217     int    i, j, nrow, ncol;
 218     real **m;
 219
 220     nrow = nrh-nrl+1;
 221     ncol = nch-ncl+1;
 222     m    = (real **) malloc((unsigned) (nrow)*sizeof(real*));
 223     if (!m)
 224     {
 225         nrerror("allocation failure in convert_matrix()", TRUE);
 226     }
 227     m -= nrl;
 228     for (i = 0, j = nrl; i <= nrow-1; i++, j++)
 229     {
 230         m[j] = a+ncol*i-ncl;
 231     }
 232     return m;
 233 }
 234
 235
 236
 237 static void free_convert_matrix(real **b, int nrl)
 238 {
 239     free((char*) (b+nrl));
 240 }
 241
 242 #define SWAP(a, b) {real temp = (a); (a) = (b); (b) = temp; }
 243
 244 static void dump_mat(int n, real **a)
 245 {
 246     int i, j;
 247
 248     for (i = 1; (i <= n); i++)
 249     {
 250         for (j = 1; (j <= n); j++)
 251         {
 252             fprintf(stderr, "  %10.3f", a[i][j]);
 253         }
 254         fprintf(stderr, "\n");
 255     }
 256 }
 257
 258 gmx_bool gaussj(real **a, int n, real **b, int m)
 259 {
 260     int *indxc, *indxr, *ipiv;
 261     int  i, icol = 0, irow = 0, j, k, l, ll;
 262     real big, dum, pivinv;
 263
 264     indxc = ivector(1, n);
 265     indxr = ivector(1, n);
 266     ipiv  = ivector(1, n);
 267     for (j = 1; j <= n; j++)
 268     {
 269         ipiv[j] = 0;
 270     }
 271     for (i = 1; i <= n; i++)
 272     {
 273         big = 0.0;
 274         for (j = 1; j <= n; j++)
 275         {
 276             if (ipiv[j] != 1)
 277             {
 278                 for (k = 1; k <= n; k++)
 279                 {
 280                     if (ipiv[k] == 0)
 281                     {
 282                         if (fabs(a[j][k]) >= big)
 283                         {
 284                             big  = fabs(a[j][k]);
 285                             irow = j;
 286                             icol = k;
 287                         }
 288                     }
 289                     else if (ipiv[k] > 1)
 290                     {
 291                         nrerror("GAUSSJ: Singular Matrix-1", FALSE);
 292                         return FALSE;
 293                     }
 294                 }
 295             }
 296         }
 297         ++(ipiv[icol]);
 298         if (irow != icol)
 299         {
 300             for (l = 1; l <= n; l++)
 301             {
 302                 SWAP(a[irow][l], a[icol][l]);
 303             }
 304             for (l = 1; l <= m; l++)
 305             {
 306                 SWAP(b[irow][l], b[icol][l]);
 307             }
 308         }
 309         indxr[i] = irow;
 310         indxc[i] = icol;
 311         if (a[icol][icol] == 0.0)
 312         {
 313             fprintf(stderr, "irow = %d, icol = %d\n", irow, icol);
 314             dump_mat(n, a);
 315             nrerror("GAUSSJ: Singular Matrix-2", FALSE);
 316             return FALSE;
 317         }
 318         pivinv        = 1.0/a[icol][icol];
 319         a[icol][icol] = 1.0;
 320         for (l = 1; l <= n; l++)
 321         {
 322             a[icol][l] *= pivinv;
 323         }
 324         for (l = 1; l <= m; l++)
 325         {
 326             b[icol][l] *= pivinv;
 327         }
 328         for (ll = 1; ll <= n; ll++)
 329         {
 330             if (ll != icol)
 331             {
 332                 dum         = a[ll][icol];
 333                 a[ll][icol] = 0.0;
 334                 for (l = 1; l <= n; l++)
 335                 {
 336                     a[ll][l] -= a[icol][l]*dum;
 337                 }
 338                 for (l = 1; l <= m; l++)
 339                 {
 340                     b[ll][l] -= b[icol][l]*dum;
 341                 }
 342             }
 343         }
 344     }
 345     for (l = n; l >= 1; l--)
 346     {
 347         if (indxr[l] != indxc[l])
 348         {
 349             for (k = 1; k <= n; k++)
 350             {
 351                 SWAP(a[k][indxr[l]], a[k][indxc[l]]);
 352             }
 353         }
 354     }
 355     free_ivector(ipiv, 1);
 356     free_ivector(indxr, 1);
 357     free_ivector(indxc, 1);
 358
 359     return TRUE;
 360 }
 361
 362 #undef SWAP
 363
 364
 365 static void covsrt(real **covar, int ma, int lista[], int mfit)
 366 {
 367     int  i, j;
 368     real swap;
 369
 370     for (j = 1; j < ma; j++)
 371     {
 372         for (i = j+1; i <= ma; i++)
 373         {
 374             covar[i][j] = 0.0;
 375         }
 376     }
 377     for (i = 1; i < mfit; i++)
 378     {
 379         for (j = i+1; j <= mfit; j++)
 380         {
 381             if (lista[j] > lista[i])
 382             {
 383                 covar[lista[j]][lista[i]] = covar[i][j];
 384             }
 385             else
 386             {
 387                 covar[lista[i]][lista[j]] = covar[i][j];
 388             }
 389         }
 390     }
 391     swap = covar[1][1];
 392     for (j = 1; j <= ma; j++)
 393     {
 394         covar[1][j] = covar[j][j];
 395         covar[j][j] = 0.0;
 396     }
 397     covar[lista[1]][lista[1]] = swap;
 398     for (j = 2; j <= mfit; j++)
 399     {
 400         covar[lista[j]][lista[j]] = covar[1][j];
 401     }
 402     for (j = 2; j <= ma; j++)
 403     {
 404         for (i = 1; i <= j-1; i++)
 405         {
 406             covar[i][j] = covar[j][i];
 407         }
 408     }
 409 }
 410
 411 #define SWAP(a, b) {swap = (a); (a) = (b); (b) = swap; }
 412
 413 static void covsrt_new(real **covar, int ma, int ia[], int mfit)
 414 /* Expand in storage the covariance matrix covar, so as to take
 415  * into account parameters that are being held fixed. (For the
 416  * latter, return zero covariances.)
 417  */
 418 {
 419     int  i, j, k;
 420     real swap;
 421     for (i = mfit+1; i <= ma; i++)
 422     {
 423         for (j = 1; j <= i; j++)
 424         {
 425             covar[i][j] = covar[j][i] = 0.0;
 426         }
 427     }
 428     k = mfit;
 429     for (j = ma; j >= 1; j--)
 430     {
 431         if (ia[j])
 432         {
 433             for (i = 1; i <= ma; i++)
 434             {
 435                 SWAP(covar[i][k], covar[i][j]);
 436             }
 437             for (i = 1; i <= ma; i++)
 438             {
 439                 SWAP(covar[k][i], covar[j][i]);
 440             }
 441             k--;
 442         }
 443     }
 444 }
 445 #undef SWAP
 446
 447 static void mrqcof(real x[], real y[], real sig[], int ndata, real a[],
 448                    int ma, int lista[], int mfit,
 449                    real **alpha, real beta[], real *chisq,
 450                    void (*funcs)(real, real *, real *, real *))
 451 {
 452     int  k, j, i;
 453     real ymod, wt, sig2i, dy, *dyda;
 454
 455     dyda = rvector(1, ma);
 456     for (j = 1; j <= mfit; j++)
 457     {
 458         for (k = 1; k <= j; k++)
 459         {
 460             alpha[j][k] = 0.0;
 461         }
 462         beta[j] = 0.0;
 463     }
 464     *chisq = 0.0;
 465     for (i = 1; i <= ndata; i++)
 466     {
 467         (*funcs)(x[i], a, &ymod, dyda);
 468         sig2i = 1.0/(sig[i]*sig[i]);
 469         dy    = y[i]-ymod;
 470         for (j = 1; j <= mfit; j++)
 471         {
 472             wt = dyda[lista[j]]*sig2i;
 473             for (k = 1; k <= j; k++)
 474             {
 475                 alpha[j][k] += wt*dyda[lista[k]];
 476             }
 477             beta[j] += dy*wt;
 478         }
 479         (*chisq) += dy*dy*sig2i;
 480     }
 481     for (j = 2; j <= mfit; j++)
 482     {
 483         for (k = 1; k <= j-1; k++)
 484         {
 485             alpha[k][j] = alpha[j][k];
 486         }
 487     }
 488     free_vector(dyda, 1);
 489 }
 490
 491
 492 gmx_bool mrqmin(real x[], real y[], real sig[], int ndata, real a[],
 493                 int ma, int lista[], int mfit,
 494                 real **covar, real **alpha, real *chisq,
 495                 void (*funcs)(real, real *, real *, real *),
 496                 real *alamda)
 497 {
 498     int          k, kk, j, ihit;
 499     static real *da, *atry, **oneda, *beta, ochisq;
 500
 501     if (*alamda < 0.0)
 502     {
 503         oneda = matrix1(1, mfit, 1, 1);
 504         atry  = rvector(1, ma);
 505         da    = rvector(1, ma);
 506         beta  = rvector(1, ma);
 507         kk    = mfit+1;
 508         for (j = 1; j <= ma; j++)
 509         {
 510             ihit = 0;
 511             for (k = 1; k <= mfit; k++)
 512             {
 513                 if (lista[k] == j)
 514                 {
 515                     ihit++;
 516                 }
 517             }
 518             if (ihit == 0)
 519             {
 520                 lista[kk++] = j;
 521             }
 522             else if (ihit > 1)
 523             {
 524                 nrerror("Bad LISTA permutation in MRQMIN-1", FALSE);
 525                 return FALSE;
 526             }
 527         }
 528         if (kk != ma+1)
 529         {
 530             nrerror("Bad LISTA permutation in MRQMIN-2", FALSE);
 531             return FALSE;
 532         }
 533         *alamda = 0.001;
 534         mrqcof(x, y, sig, ndata, a, ma, lista, mfit, alpha, beta, chisq, funcs);
 535         ochisq = (*chisq);
 536     }
 537     for (j = 1; j <= mfit; j++)
 538     {
 539         for (k = 1; k <= mfit; k++)
 540         {
 541             covar[j][k] = alpha[j][k];
 542         }
 543         covar[j][j] = alpha[j][j]*(1.0+(*alamda));
 544         oneda[j][1] = beta[j];
 545     }
 546     if (!gaussj(covar, mfit, oneda, 1))
 547     {
 548         return FALSE;
 549     }
 550     for (j = 1; j <= mfit; j++)
 551     {
 552         da[j] = oneda[j][1];
 553     }
 554     if (*alamda == 0.0)
 555     {
 556         covsrt(covar, ma, lista, mfit);
 557         free_vector(beta, 1);
 558         free_vector(da, 1);
 559         free_vector(atry, 1);
 560         free_matrix(oneda, 1, mfit, 1);
 561         return TRUE;
 562     }
 563     for (j = 1; j <= ma; j++)
 564     {
 565         atry[j] = a[j];
 566     }
 567     for (j = 1; j <= mfit; j++)
 568     {
 569         atry[lista[j]] = a[lista[j]]+da[j];
 570     }
 571     mrqcof(x, y, sig, ndata, atry, ma, lista, mfit, covar, da, chisq, funcs);
 572     if (*chisq < ochisq)
 573     {
 574         *alamda *= 0.1;
 575         ochisq   = (*chisq);
 576         for (j = 1; j <= mfit; j++)
 577         {
 578             for (k = 1; k <= mfit; k++)
 579             {
 580                 alpha[j][k] = covar[j][k];
 581             }
 582             beta[j]     = da[j];
 583             a[lista[j]] = atry[lista[j]];
 584         }
 585     }
 586     else
 587     {
 588         *alamda *= 10.0;
 589         *chisq   = ochisq;
 590     }
 591     return TRUE;
 592 }
 593
 594
 595 gmx_bool mrqmin_new(real x[], real y[], real sig[], int ndata, real a[],
 596                     int ia[], int ma, real **covar, real **alpha, real *chisq,
 597                     void (*funcs)(real, real [], real *, real []),
 598                     real *alamda)
 599 /* Levenberg-Marquardt method, attempting to reduce the value Chi^2
 600  * of a fit between a set of data points x[1..ndata], y[1..ndata]
 601  * with individual standard deviations sig[1..ndata], and a nonlinear
 602  * function dependent on ma coefficients a[1..ma]. The input array
 603  * ia[1..ma] indicates by nonzero entries those components of a that
 604  * should be fitted for, and by zero entries those components that
 605  * should be held fixed at their input values. The program returns
 606  * current best-fit values for the parameters a[1..ma], and
 607  * Chi^2 = chisq. The arrays covar[1..ma][1..ma], alpha[1..ma][1..ma]
 608  * are used as working space during most iterations. Supply a routine
 609  * funcs(x,a,yfit,dyda,ma) that evaluates the fitting function yfit,
 610  * and its derivatives dyda[1..ma] with respect to the fitting
 611  * parameters a at x. On the first call provide an initial guess for
 612  * the parameters a, and set alamda < 0 for initialization (which then
 613  * sets alamda=.001). If a step succeeds chisq becomes smaller and
 614  * alamda de-creases by a factor of 10. If a step fails alamda grows by
 615  * a factor of 10. You must call this routine repeatedly until
 616  * convergence is achieved. Then, make one final call with alamda=0,
 617  * so that covar[1..ma][1..ma] returns the covariance matrix, and alpha
 618  * the curvature matrix.
 619  * (Parameters held fixed will return zero covariances.)
 620  */
 621 {
 622     void covsrt(real **covar, int ma, int ia[], int mfit);
 623     gmx_bool gaussj(real **a, int n, real **b, int m);
 624     void mrqcof_new(real x[], real y[], real sig[], int ndata, real a[],
 625                     int ia[], int ma, real **alpha, real beta[], real *chisq,
 626                     void (*funcs)(real, real [], real *, real []));
 627     int         j, k, l;
 628     static int  mfit;
 629     static real ochisq, *atry, *beta, *da, **oneda;
 630
 631     if (*alamda < 0.0)                    /* Initialization. */
 632     {
 633         atry = rvector(1, ma);
 634         beta = rvector(1, ma);
 635         da   = rvector(1, ma);
 636         for (mfit = 0, j = 1; j <= ma; j++)
 637         {
 638             if (ia[j])
 639             {
 640                 mfit++;
 641             }
 642         }
 643         oneda   = matrix1(1, mfit, 1, 1);
 644         *alamda = 0.001;
 645         mrqcof_new(x, y, sig, ndata, a, ia, ma, alpha, beta, chisq, funcs);
 646         ochisq = (*chisq);
 647         for (j = 1; j <= ma; j++)
 648         {
 649             atry[j] = a[j];
 650         }
 651     }
 652     for (j = 1; j <= mfit; j++) /* Alter linearized fitting matrix, by augmenting. */
 653     {
 654         for (k = 1; k <= mfit; k++)
 655         {
 656             covar[j][k] = alpha[j][k]; /* diagonal elements. */
 657         }
 658         covar[j][j] = alpha[j][j]*(1.0+(*alamda));
 659         oneda[j][1] = beta[j];
 660     }
 661     if (!gaussj(covar, mfit, oneda, 1)) /* Matrix solution. */
 662     {
 663         return FALSE;
 664     }
 665     for (j = 1; j <= mfit; j++)
 666     {
 667         da[j] = oneda[j][1];
 668     }
 669     if (*alamda == 0.0) /* Once converged, evaluate covariance matrix. */
 670     {
 671         covsrt_new(covar, ma, ia, mfit);
 672         free_matrix(oneda, 1, mfit, 1);
 673         free_vector(da, 1);
 674         free_vector(beta, 1);
 675         free_vector(atry, 1);
 676         return TRUE;
 677     }
 678     for (j = 0, l = 1; l <= ma; l++) /* Did the trial succeed? */
 679     {
 680         if (ia[l])
 681         {
 682             atry[l] = a[l]+da[++j];
 683         }
 684     }
 685     mrqcof_new(x, y, sig, ndata, atry, ia, ma, covar, da, chisq, funcs);
 686     if (*chisq < ochisq)
 687     {
 688         /* Success, accept the new solution. */
 689         *alamda *= 0.1;
 690         ochisq   = (*chisq);
 691         for (j = 1; j <= mfit; j++)
 692         {
 693             for (k = 1; k <= mfit; k++)
 694             {
 695                 alpha[j][k] = covar[j][k];
 696             }
 697             beta[j] = da[j];
 698         }
 699         for (l = 1; l <= ma; l++)
 700         {
 701             a[l] = atry[l];
 702         }
 703     }
 704     else                 /* Failure, increase alamda and return. */
 705     {
 706         *alamda *= 10.0;
 707         *chisq   = ochisq;
 708     }
 709     return TRUE;
 710 }
 711
 712 void mrqcof_new(real x[], real y[], real sig[], int ndata, real a[],
 713                 int ia[], int ma, real **alpha, real beta[], real *chisq,
 714                 void (*funcs)(real, real [], real *, real[]))
 715 /* Used by mrqmin to evaluate the linearized fitting matrix alpha, and
 716  * vector beta as in (15.5.8), and calculate Chi^2.
 717  */
 718 {
 719     int  i, j, k, l, m, mfit = 0;
 720     real ymod, wt, sig2i, dy, *dyda;
 721
 722     dyda = rvector(1, ma);
 723     for (j = 1; j <= ma; j++)
 724     {
 725         if (ia[j])
 726         {
 727             mfit++;
 728         }
 729     }
 730     for (j = 1; j <= mfit; j++) /* Initialize (symmetric) alpha), beta. */
 731     {
 732         for (k = 1; k <= j; k++)
 733         {
 734             alpha[j][k] = 0.0;
 735         }
 736         beta[j] = 0.0;
 737     }
 738     *chisq = 0.0;
 739     for (i = 1; i <= ndata; i++) /* Summation loop over all data. */
 740     {
 741         (*funcs)(x[i], a, &ymod, dyda);
 742         sig2i = 1.0/(sig[i]*sig[i]);
 743         dy    = y[i]-ymod;
 744         for (j = 0, l = 1; l <= ma; l++)
 745         {
 746             if (ia[l])
 747             {
 748                 wt = dyda[l]*sig2i;
 749                 for (j++, k = 0, m = 1; m <= l; m++)
 750                 {
 751                     if (ia[m])
 752                     {
 753                         alpha[j][++k] += wt*dyda[m];
 754                     }
 755                 }
 756                 beta[j] += dy*wt;
 757             }
 758         }
 759         *chisq += dy*dy*sig2i;  /* And find Chi^2. */
 760     }
 761     for (j = 2; j <= mfit; j++) /* Fill in the symmetric side. */
 762     {
 763         for (k = 1; k < j; k++)
 764         {
 765             alpha[k][j] = alpha[j][k];
 766         }
 767     }
 768     free_vector(dyda, 1);
 769 }