grace_dofit.c

Go to the documentation of this file.

/* 
Copyright 2021 Arizona Board of Regents on behalf of the University of Arizona.

For commercial uses, to obtain a license to sell and or sublicense copies of
the Software please contact the University of Arizona at Tech Launch Arizona:
lewish@tla.arizona.edu.

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction for non-commercial use educational and
research purposes, including without limitation the rights to use, copy,
modify, merge, publish, distribute, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:

The entire above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

1. Additional Required Provisions

1.1.  Arbitration. The parties agree that if a dispute arises between them
concerning this Agreement, the parties may be required to submit the matter to
arbitration pursuant to Arizona law.

1.2.  Applicable Law and Venue. This Agreement shall be interpreted pursuant
to the laws of the State of Arizona. Any arbitration or litigation between the
Parties shall be conducted in Pima County, ARIZONA, and LICENSEE hereby
submits to venue and jurisdiction in Pima County, ARIZONA.

1.3.  Non-Discrimination. The Parties agree to be bound by state and federal
laws and regulations governing equal opportunity and non-discrimination and
immigration.

1.4.   Appropriation of Funds. The Parties recognize that performance by
ARIZONA may depend upon appropriation of funds by the State Legislature of
ARIZONA. If the Legislature fails to appropriate the necessary funds, or if
ARIZONA’S appropriation is reduced during the fiscal year, ARIZONA may cancel
this Agreement without further duty or obligation. ARIZONA will notify
LICENSEE as soon as reasonably possible after it knows of the loss of funds.

1.5.   Conflict of Interest. This Agreement is subject to the provisions of
A.R.S. 38-511 and other conflict of interest regulations. Within three years
of the EFFECTIVE DATE, ARIZONA may cancel this Agreement if any person
significantly involved in initiating, negotiating, drafting, securing, or
creating this Agreement for or on behalf of ARIZONA becomes an employee or
consultant in any capacity of LICENSEE with respect to the subject matter of
this Agreement.

 */

/*
 * Grace - Graphics for Exploratory Data Analysis
 * 
 * Home page: http://plasma-gate.weizmann.ac.il/Grace/
 * 
 * Copyright (c) 1991-95 Paul J Turner, Portland, OR
 * Copyright (c) 1996-98 GRACE Development Team
 * 
 * Maintained by Evgeny Stambulchik <fnevgeny@plasma-gate.weizmann.ac.il>
 * 
 * 
 *                           All Rights Reserved
 * 
 *    This program is free software; you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation; either version 2 of the License, or
 *    (at your option) any later version.
 * 
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 * 
 *    You should have received a copy of the GNU General Public License
 *    along with this program; if not, write to the Free Software
 *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include "extern.h"

/* ALGORITHM AS274.1  APPL. STATIST. (1992) VOL 41, NO. 2

   Calling this routine updates d, rbar, thetab and sserr by the
   inclusion of xrow, yelem with the specified weight.  This version has
   been modified to make it slightly faster when the early elements of
   XROW are not zeroes.  *** WARNING *** The elements of XROW are
   over-written. */
int gf_includ(np, nrbar, w, xrow, y, d, rbar, thetab, sserr)
int np, nrbar;
double w, *xrow, y, *d, *rbar, *thetab, *sserr;
{
  int ier=0, i, k, nextr;
  double cbar, sbar, di, xi, xk, dpi, wxi;


  --thetab; 
  --rbar; 
  --d; 
  --xrow;

    /* the function includ() is called for each obs, and is where all
       the time is spent.  To make it faster, we skip the checks about
       whether nrbar is big enough for np.  If you'd like to have the
       checks in anyway, then compile with -DSTRINGENT */
#ifdef STRINGENT /* Some checks. */
  if(np < 1)
  {
    ier = 1;
  }

  if(nrbar < np * (np - 1) / 2)
  {
     ier += 2;
  }

  if(ier)
  {
    return(ier);
  }

#endif /* STRINGENT */

  nextr = 1;

/* Skip unnecessary transformations.  Test on exact zeroes must be used or stability can be destroyed. */

  for(i=1;i<=np;++i)
  {
    if(w == 0)
    {
      return(ier);
    }

    xi = xrow[i];

    if(xi == 0)
    {
      nextr += np - i;
    }
    else
    {
      di   = d[i];
      wxi  = w * xi;
      dpi  = di + wxi * xi;
      cbar = di / dpi;
      sbar = wxi / dpi;
      w    = cbar * w;
      d[i] = dpi;

      if(i != np)
      {
        for(k=i+1;k<=np;++k)
        {
          xk          = xrow[k];
          xrow[k]     = xk - xi * rbar[nextr];
          rbar[nextr] = cbar * rbar[nextr] + sbar * xk;
          ++nextr;
        }
      }

      xk        = y;
      y         = xk - xi * thetab[i];
      thetab[i] = cbar * thetab[i] + sbar * xk;
    }
  }

/* Y * SQRT(W) is now equal to the Brown, Durbin & Evans recursive residual. */
  *sserr += w * y * y;

  return(ier);
}



/* ALGORITHM AS274.2  APPL. STATIST. (1992) VOL.41, NO. 2
   Sets arrays to zero prior to calling AS75.1 */
int gf_clear(np, nrbar, d, rbar, thetab, sserr)
int np, nrbar;
double *d, *rbar, *thetab, *sserr;
{
  int ier=0;

  if(np < 1)
  {
    ier = 1;
  }

  if(nrbar < (np * (np - 1) / 2))
  {
    ier += 2;
  }

  if(ier)
  {
    return(ier);
  }

  memset(d,      0, np    * sizeof(double));
  memset(thetab, 0, np    * sizeof(double));
  memset(rbar,   0, nrbar * sizeof(double));

  *sserr = 0;

  return(0);
}


int gf_regcf(int np, int nrbar, double *d, double *rbar,
          double *thetab, double *tol, double *beta,
          int nreq)
{
    int ier=0, i, j, nextr;

/* ALGORITHM AS274.3  APPL. STATIST. (1992) VOL 41, NO.2
   Modified version of AS75.4 to calculate regression coefficients
   for the first NREQ variables, given an orthogonal reduction from
   AS75.1. */

    --beta; --tol; --thetab; --rbar; --d;

    /*     Some checks. */
    if (np < 1) ier = 1;
    if (nrbar < np * (np - 1) / 2) ier += 2;
    if (nreq < 1 || nreq > np) ier += 4;
    if (ier) return ier;

    for (i = nreq; i >= 1; --i)
      if (sqrt(d[i]) < tol[i]) {
          beta[i] = d[i] = 0;
      } else {
          beta[i] = thetab[i];
          nextr = (i - 1) * (np + np - i) / 2 + 1;
          for (j=i+1; j<=nreq; ++j) beta[i] -= rbar[nextr++] * beta[j];
      }

    return ier;
} /* regcf */


int gf_tolset(int np, int nrbar, double *d, double *rbar, double *tol)
/* ALGORITHM AS274.4  APPL. STATIST. (1992) VOL 41, NO.2
   Sets up array TOL for testing for zeroes in an orthogonal
   reduction formed using AS75.1.
   EPS is a machine-dependent constant.   For compilers which use
   the IEEE format for floating-point numbers, recommended values
   are 1.E-06 for single precision and 1.E-15 for double precision. */
{
    const double eps = 1e-15;
    double d__1, sum, *work;
    int ier=0, col, pos, row;

    --tol; --rbar; --d;

    /* Some checks. */
    if (np < 1) ier = 1;
    if (nrbar < np * (np - 1) / 2) ier += 2;
    if (ier) return ier;

/* Set TOL(I) = sum of absolute values in column I of RBAR after
   scaling each element by the square root of its row multiplier. */

    work = -1 + (double *) malloc(np*sizeof(double));
    for (col = 1; col <= np; ++col) work[col] = sqrt(d[col]);
    for (col = 1; col <= np; ++col) {
        pos = col - 1;
        sum = work[col];
        for (row = 1; row <= col-1; ++row) {
            sum += (d__1 = rbar[pos], fabs(d__1)) * work[row];
            pos += np - row - 1;
        }
        tol[col] = eps * sum;
    }
    free(work+1);
    return 0;
} /* tolset */


int gf_sing(int np, int nrbar, double *d,
         double *rbar, double *thetab, double *sserr,
         double *tol, int *lindep)
/* ALGORITHM AS274.5  APPL. STATIST. (1992) VOL 41, NO. 2
   Checks for singularities, reports, and adjusts orthogonal
   reductions produced by AS75.1. */
{
    double d__1, temp, *work;
    int ier=0, np2, col, pos, row, pos2;

    --lindep; --tol; --thetab; --rbar; --d;

    if (np <= 0) ier = 1;
    if (nrbar < np * (np - 1) / 2) ier += 2;
    if (ier) return ier;
    work = -1 + (double *) malloc(np*sizeof(double));

/* Set elements within RBAR to zero if they are less than TOL(COL) in
   absolute value after being scaled by the square root of their row
   multiplier. */
    for (col = 1; col <= np; ++col) work[col] = sqrt(d[col]);
    for (col = 1; col <= np; ++col) {
        temp = tol[col];
        pos = col - 1;
        for (row = 1; row < col; ++row) {
            if ((d__1 = rbar[pos], fabs(d__1)) * work[row] < temp)
              rbar[pos] = 0;
            pos += np - row - 1;
        }

/* If diagonal element is near zero, set it to zero, set appropriate
   element of LINDEP, and use INCLUD to augment the projections in
   the lower rows of the orthogonalization. */

        lindep[col] = 0;
        if (work[col] < temp) {
            lindep[col] = 1;
            --ier;
            if (col < np) {
                np2 = np - col;
                pos2 = pos + np2 + 1;
                ier = gf_includ(np2, np2*(np2-1)/2, d[col], rbar+pos+1,
                             thetab[col], d+col+1, rbar+pos2,
                             thetab+col+1, sserr);
            } else {
                d__1 = thetab[col]; /* for squaring */
                *sserr += d[col] * (d__1 * d__1);
            }
            d[col] = work[col] = thetab[col] = 0;
        }
    }
    free(work+1);
    return ier;
} /* sing */


int gf_ss(int np, double *d, double *thetab, double *sserr, double *rss)
/* ALGORITHM AS274.6  APPL. STATIST. (1992) VOL 41, NO. 2
   Calculates partial residual sums of squares from an orthogonal
   reduction from AS75.1. */
{
    double sum;
    int i;

    --rss; --thetab; --d;

    if (np < 1) return 1;
    sum = rss[np] = *sserr;
    for (i = np; i >= 2; --i) {
        sum += d[i] * thetab[i] * thetab[i];
        rss[i - 1] = sum;
    }
    return 0;
} /* ss */


int gf_cov(int np, int nrbar, double *d, 
         double *rbar, int nreq, double *rinv, double *var, 
         double *covmat, int dimcov, double *sterr)
/* ALGORITHM AS274.7  APPL. STATIST. (1992) VOL 41, NO.2 
   Calculate covariance matrix for regression coefficients for the 
   first NREQ variables, from an orthogonal reduction produced from 
   AS75.1. */
{
    int ier=0, k, start, col, pos, row, pos1, pos2;
    double sum;

    --sterr; --covmat; --rinv; --rbar; --d;

    if (np < 1) ier = 1;
    if (nrbar < np * (np - 1) / 2) ier += 2;
    if (dimcov < nreq * (nreq + 1) / 2) ier += 4;

    for (row = 1; row <= nreq; ++row) if (d[row] == 0) ier = -row;
    if (ier) return ier;

    gf_inv(np, nrbar, rbar+1, nreq, rinv+1);
    start = pos = 1;
    for (row = 1; row <= nreq; ++row) {
        pos2 = start;
        for (col = row; col <= nreq; ++col) {
            pos1 = start + col - row;
            sum = (row == col ? 1.0 : rinv[pos1-1])/d[col];
            for (k = col + 1; k <= nreq; ++k)
              sum += rinv[pos1++] * rinv[pos2++] / d[k];
            covmat[pos] = sum * *var;
            if (row == col) sterr[row] = sqrt(covmat[pos]);
            ++pos;
        }
        start += nreq - row;
    }
    return 0;
} /* cov */


void gf_inv(int np, int nrbar, double *rbar, int nreq, double *rinv)
/* ALGORITHM AS274.8  APPL. STATIST. (1992) VOL 41, NO. 2
   Invert first NREQ rows and columns of Cholesky factorization
   produced by AS75.1. */
{
    int pos, row;

    --rinv; --rbar;

    /* Invert RBAR ignoring row multipliers, from the bottom up. */
    pos = nreq * (nreq - 1) / 2;
    for (row = nreq - 1; row >= 1; --row) {
        int col, start = (row - 1) * (np + np - row) / 2 + 1;
        for (col = nreq; col >= row+1; --col) {
            double sum=0;
            int pos1=start, pos2=pos, k;
            for (k = row + 1; k <= col-1; ++k) {
                pos2 += nreq - k;
                sum -= rbar[pos1++] * rinv[pos2];
            }
            rinv[pos] = sum - rbar[pos1];
            --pos;
        }
    }
} /* inv */


int gf_pcorr(int np, int nrbar, double *d, 
          double *rbar, double *thetab, double *sserr, int in, 
          double *cormat, int dimc, double *ycorr)
/* ALGORITHM AS274.9  APPL. STATIST. (1992) VOL 41, NO. 2
   Calculate partial correlations after the first IN variables
   have been forced into the regression. */
{
    int ier=0, i, start, in1;
    double *work;

    --cormat; --thetab; --rbar; --d;

    if (np < 1) ier = 1;
    if (nrbar < np * (np - 1) / 2) ier += 2;
    if (in < 0 || in > np - 1) ier += 4;
    if (dimc < (np - in) * (np - in - 1) / 2) ier += 8;
    if (ier) return ier;

    start = in * (np + np - in - 1) / 2 + 1;
    in1 = in + 1;
    work = (double *) malloc(np*sizeof(double));
    gf_cor(np-in, d+in1, rbar+start, thetab+in1, sserr, work, cormat+1, ycorr);

    /* Check for zeroes. */
    for (i=0; i<(np-in); ++i) if (work[i] <= 0) ier = -1-i;
    free(work);
    return ier;
} /* pcorr */


void gf_cor(int np, double *d, double *rbar, 
         double *thetab, double *sserr, double *work,
         double *cormat, double *ycorr)
/* ALGORITHM AS274.10  APPL. STATIST. (1992) VOL 41, NO.2
   Calculate correlations from an orthogonal reduction.   This
   routine will usually be called from PCORR, which will have
   removed the appropriate number of rows at the start. */
{
    int i__1, i__2, diff, pos, row, col1, col2, pos1, pos2;
    double sumy, sum;

    --ycorr; --cormat; --work; --thetab; --rbar; --d;

/* process by columns, including the projections of the dependent
   variable (THETAB). */
    sumy = *sserr;
    i__1 = np;
    for (row = 1; row <= i__1; ++row)
      sumy += d[row] * thetab[row] * thetab[row];
    sumy = sqrt(sumy);

    pos = np * (np - 1) / 2;
    for (col1 = np; col1 >= 1; --col1) { /* find length of column COL1. */
        sum = d[col1];
        i__1 = pos1 = col1 - 1;
        for (row = 1; row <= i__1; ++row) {
            sum += d[row] * rbar[pos1] * rbar[pos1];
            pos1 += np - row - 1;
        }
        work[col1] = sqrt(sum);

/* If SUM = 0, set all correlations with this variable to zero. */
        if (sum == 0) {
            ycorr[col1] = 0;
            i__1 = col1 + 1;
            for (col2 = np; col2 >= i__1; --col2) {
                cormat[pos] = 0; --pos;
            }
            goto L70;
        }

/* Form cross-products, then divide by product of column lengths. */
        sum = d[col1] * thetab[col1];
        i__1 = pos1 = col1 - 1;
        for (row = 1; row <= i__1; ++row) {
            sum += d[row] * rbar[pos1] * thetab[row];
            pos1 += np - row - 1;
        }
        ycorr[col1] = sum / (sumy * work[col1]);

        i__1 = col1 + 1;
        for (col2 = np; col2 >= i__1; --col2) {
            if (work[col2] > 0) {
                pos1 = col1 - 1;
                pos2 = col2 - 1;
                diff = col2 - col1;
                i__2 = col1 - 1;
                for (sum=0, row = 1; row <= i__2; ++row) {
                    sum += d[row] * rbar[pos1] * rbar[pos2];
                    pos1 += np - row - 1;
                    pos2 = pos1 + diff;
                }
                sum += d[col1] * rbar[pos2];
                cormat[pos] = sum / (work[col1] * work[col2]);
            } else cormat[pos] = 0;
            --pos;
        }
L70:
        ;
    }
} /* cor */




/* ALGORITHM AS274.11 APPL. STATIST. (1992) VOL 41, NO. 2
 * Move variable from position FROM to position TO in an
 * orthogonal reduction produced by AS75.1.
 */
int gf_vmove(np, nrbar, vorder, d, rbar, thetab, rss, from, to, tol)
int np, nrbar, *vorder;
double *d, *rbar, *thetab, *rss;
int from, to;
double *tol;
{
  int ier=0, last, m;
  double d__1, cbar, sbar, d1new, d2new, x, y, d1, d2;
  int first, m1, m2, mp1, inc, col, pos, row;

  --tol;
  --rss;
  --thetab;
  --rbar;
  --d;
  --vorder;

  if(np <= 0)
  {
    ier = 1;
  }

  if(nrbar < np * (np - 1) / 2)
  {
    ier += 2;
  }

  if(from < 1 || from > np)
  {
    ier += 4;
  }

  if(to < 1 || to > np)
  {
    ier += 8;
  }

  if(ier)
  {
    return ier;
  }

  if(from == to)
  {
    return ier;
  }

  if(from < to)
  {
    first = from;
    last  = to - 1;
    inc   = 1;
  }
  else
  {
    first = from - 1;
    last  = to;
    inc   = -1;
  }

/* Find addresses of first elements of RBAR in rows M and (M+1). */    
  for (m = first; inc < 0 ? m >= last : m <= last; m += inc)
  {
    m1  = (m - 1) * (np + np - m) / 2 + 1;
    m2  = m1 + np - m;
    mp1 = m + 1;
    d1  = d[m];
    d2  = d[mp1];

        /* Special cases. */
    if(d1 == 0 && d2 == 0)
    {
      goto L40;
    }

    x = rbar[m1];
    if(fabs(x) * sqrt(d1) < tol[mp1])
    {
      x = 0;
    }

    if(d1 == 0 || x == 0)
    {
      d[m] = d2;
      d[mp1] = d1;
      rbar[m1] = 0;

      for (col = m + 2; col <= np; ++col)
      {
        x          = rbar[++m1];
        rbar[m1]   = rbar[m2];
        rbar[m2++] = x;
      }

      x           = thetab[m];
      thetab[m]   = thetab[mp1];
      thetab[mp1] = x;

      goto L40;
    }
    else
    if(d2 == 0)
    {
      d__1 = x;
      d[m] = d1 * d__1 * d__1;
      rbar[m1] = 1.0/x;

      for(col = m + 2; col <= np; ++col)
      {
        rbar[++m1] /= x;
        thetab[m] /= x;
      }

      goto L40;
    }

/*  Planar rotation in regular case. */
    d__1     = x;
    d1new    = d2 + d1 * (d__1 * d__1);
    cbar     = d2 / d1new;
    sbar     = x * d1 / d1new;
    d2new    = d1 * cbar;
    d[m]     = d1new;
    d[mp1]   = d2new;
    rbar[m1] = sbar;

    for(col = m + 2; col <= np; ++col)
    {
      y        = rbar[++m1];
      rbar[m1] = cbar * rbar[m2] + sbar * y;
      rbar[m2] = y - x * rbar[m2];
      ++m2;
    }

    y           = thetab[m];
    thetab[m]   = cbar * thetab[mp1] + sbar * y;
    thetab[mp1] = y - x * thetab[mp1];

/*     Swap columns M and (M+1) down to row (M-1). */
L40:
    if(m == 1)
    {
      goto L60;
    }

    pos = m;

    for(row = 1; row < m; ++row)
    {
      x             = rbar[pos];
      rbar[pos]     = rbar[pos - 1];
      rbar[pos - 1] = x;
      pos          += np - row - 1;
    }

/* Adjust variable order (VORDER), the tolerances (TOL) and the vector
   of residual sums of squares (RSS). */
L60:
    m1          = vorder[m];
    vorder[m]   = vorder[mp1];
    vorder[mp1] = m1;
    x           = tol[m];
    tol[m]      = tol[mp1];
    tol[mp1]    = x;
    d__1        = thetab[mp1];
    rss[m]      = rss[mp1] + d[mp1] * (d__1 * d__1);
  }

  return(0);
} /* gf_vmove */


/* ALGORITHM AS274.12 APPL. STATIST. (1992) VOL 41, NO.2 Re-order the
 * variables in an orthogonal reduction produced by AS75.1 so that the N
 * variables in LIST start at position POS1, though will not necessarily
 * be in the same order as in LIST.  Any variables in VORDER before
 * position POS1 are not moved.  Auxiliary routine called: VMOVE
 */
int gf_reordr(np, nrbar, vorder, d, rbar, thetab, rss, tol, list, n, pos1)
int np, nrbar, *vorder;
double *d, *rbar, *thetab, *rss, *tol;
int *list, n, pos1;
{
  int ier=0, next, i, j;

  --list;
  --tol;
  --rss;
  --thetab;
  --rbar;
  --d;
  --vorder;

  if(np < 1)
  {
    ier = 1;
  }

  if(nrbar < np * (np - 1) / 2)
  {
    ier += 2;
  }

  if(n < 1 || n >= np + 1 - pos1)
  {
    ier += 4;
  }

  if(ier)
  {
    return(ier);
  }

/* Work through VORDER finding variables which are in LIST. */
  i = next = pos1;
  do
  {
    int l=vorder[i];

    for (j=1; j<=n; j++)
    {
      if (l==list[j])
      {
        if (i > next) 
        {
          gf_vmove(np, nrbar, vorder+1, d+1, rbar+1, thetab+1, rss+1, i, next, tol+1);
        }

        if(++next >= n + pos1)
        {
          return 0;
        }

        break;
      }
    }
  }while(++i <= np);

  return(8);
} /* gf_reordr */



/* ALGORITHM AS274.13  APPL. STATIST. (1992) VOL.41, NO.2 */     
int gf_hdiag(double *xrow, int np, int nrbar, double *d, double *rbar, double *tol, int nreq, double *hii)
{
  int i__1, i__2, ier=0, col, pos, row;
  double d__1, sum, *wk;

  --tol;
  --rbar;
  --d;
  --xrow;

  if(np < 1)
  {
    ier = 1;
  }

  if(nrbar < np * (np - 1) / 2)
  {
    ier += 2;
  }

  if(nreq > np)
  {
    ier += 4;
  }

  if(ier)
  {
    return ier;
  }


/* The elements of XROW.inv(RBAR).sqrt(D) are calculated and stored in WK. */
  wk = -1 + (double *) malloc(np*sizeof(double));
  *hii = 0;
  i__1 = nreq;

  for(col=1;col<= i__1;++col)
  {
    if(sqrt(d[col]) <= tol[col])
    {
      wk[col] = 0;
    }
    else
    {
      i__2 = pos = col - 1;
      sum = xrow[col];

      for(row = 1; row <= i__2; ++row)
      {
        sum -= wk[row] * rbar[pos];
        pos += np - row - 1;
      }

      d__1 = wk[col] = sum;
      *hii += d__1 * d__1 / d[col];
    }
  }

  free(wk+1);

  return(0);
} /* gf_hdiag */




/* "print a upper triangular double matrix stored as a vector"
 *  The matrix is N x N, the vector has N*(N+1)/2 elements.
 *  Each element is formatted using width and precision.
 *  There are no sanity checks at all.
 */
void gf_pr_utdm_v(double *x, int N, int width, int precision)
{
  int pos=0, i, j, leavespace;
  char s[100], fmt[100];      /* will be used in making printf() formats */
  char buf[256];

  sprintf(fmt, "%%%d.%dg", width, precision);

  for(i=0;i<N;i++)
  {
    leavespace = i*width;
    sprintf(s, "%%%ds", leavespace); 
    sprintf(buf, s, "");
    st_print(buf);

    for(j=i;j<N;j++)
    {
      sprintf(buf, fmt, x[pos++]);
      st_print(buf);
    }

    st_print("\n");
  }
}



double *gf_dvector(int l, int h)
{
  double *block;
  int size;

  size = h - l + 1;
  block = (double *) malloc(sizeof(double)*size);

  if(block == NULL)
  {
    fprintf(stderr, "malloc failure in dvector()\n");

    return NULL;
  }

  return(block-l);
}



int *gf_ivector(int l, int h)
{
  int *block, size;

  size  = h - l + 1;
  block = (int *) malloc(sizeof(int)*size);

  if(block == NULL)
  {
    fprintf(stderr, "malloc failure in ivector()\n");

    return NULL;
  }

  return(block-l);
}



double **gf_dmatrix(int rl, int rh, int cl, int ch)
{
  double *block;
  double **m;
  int size, i, rowsize, numrows;

  rowsize = ch - cl + 1;        /* #locations consumed by 1 row */
  numrows = rh -rl + 1;
  size    = numrows*rowsize;
  block   = (double *) malloc((unsigned) sizeof(double)*size);

  if(block == NULL)
  {
    fprintf(stderr, "malloc failure in matrix allocation\n");

    return(NULL);
  }
    /* so we have the matrix. */

  /* Now for the row pointers */
  m = (double **) malloc((unsigned) sizeof(double *)*numrows);

  if(m == NULL)
  {
    fprintf(stderr, "malloc failure in matrix allocation\n");

    return(NULL);
  }
  m -= rl;      /* fixup m pointer so m[rl] == old m[0] */

    /* Finally, setup pointers to rows */
  block -= cl;

  for(i=rl;i<=rh;i++)
  {
    m[i] = block; block += rowsize;
  }

  return(m);
}



void gf_putdvec(const char *s, double *x, int l, int h)
{
  int i;
  char buf[512];
  sprintf(buf, "Vector %-10s: \n", s);

  st_print(buf);

  for(i=l;i<=h;i++)
  {
    sprintf(buf, " %d: %.4g \n", i, x[i]);
    st_print(buf);
  }
}


int dofitcurve(int cnt, double *xd, double *yd, int nd, double *c)
{
  double **x, *y, *d, *rbar, *thetab, *xrow, sserr, *beta;
  double *tol, *rss;
  int i, j, nvars = nd + 1, nrbar = 200, *vorder, *lindep, error;
  double xval;
  char buf[256];

  x      = gf_dmatrix(0, cnt - 1, 0, nvars - 1);
  y      = gf_dvector(0, cnt - 1);
  d      = gf_dvector(0, nvars - 1);
  rbar   = gf_dvector(0, nrbar);
  thetab = gf_dvector(0, nvars - 1);
  xrow   = gf_dvector(0, nvars - 1);
  beta   = gf_dvector(0, nvars - 1);
  tol    = gf_dvector(0, nvars - 1);
  rss    = gf_dvector(0, nvars - 1);
  vorder = gf_ivector(0, nvars - 1);
  lindep = gf_ivector(0, nvars - 1);

  for(i=0;i<nvars;i++)
  {
    vorder[i] = i;
  }

  for(i=0;i<cnt;i++)
  {
    xval = xd[i];

    for (j=1;j<nvars;j++)
    {
      x[i][j - 1] = pow(xval, (double) j);
    }

    y[i] = yd[i];
  }

  error = gf_clear(nvars, nrbar, d, rbar, thetab, &sserr);

  if(error)
  {
    sprintf(buf, "as274c: clear() returned %d", error);
    st_fail(buf);
    goto bustout;
  }

  for (i=0; i<cnt; i++)
  {
    xrow[0] = 1; /* include constant */

    for(j=0;j<nvars-1;j++)
    {
      xrow[j+1] = x[i][j];
    }

    error = gf_includ(nvars, nrbar, 1.0, xrow, y[i], d, rbar, thetab, &sserr);

    if(error)
    {
      sprintf(buf, "as274c: includ() returned %d at row %d", error, i);
      st_fail(buf);
      goto bustout;
    }
  }

  error = gf_tolset(nvars, nrbar, d, rbar, tol);

  if(error)
  {
    sprintf(buf, "as274c: tolset() returned %d\n", error);
    st_fail(buf);
    goto bustout;
  }

  error = gf_sing(nvars, nrbar, d, rbar, thetab, &sserr, tol, lindep);

  if(error)
  {
    sprintf(buf, "as274c: sing() returned %d", error);
    st_fail(buf);
    goto bustout;
  }

  if(isSet(PF_V_FLAG))
  {
    sprintf(buf, "\nSSerr = %17g\n", sserr);
    st_print(buf);
  }

  error = gf_regcf(nvars, nrbar, d, rbar, thetab, tol, beta, nvars);

  if(error)
  {
    sprintf(buf, "as274c: regcf() returned %d", error);
    st_fail(buf);
    goto bustout;
  }

  if(isSet(PF_V_FLAG))
  {
    st_print("\nVariable order:\n ");

    for(j=0;j<nvars;j++)
    {
      sprintf(buf, "   %d ", vorder[j]);
      st_print(buf);
    }

    st_print("\n");

    gf_putdvec("Beta", beta, 0, nvars - 1);
    gf_putdvec("d", d, 0, nvars - 1);

    sprintf(buf, "rbar matrix:\n"); 
    st_print(buf);

    gf_pr_utdm_v(rbar, nvars - 1, 13, 5);
    gf_putdvec("thetab", thetab, 0, nvars - 1);
  }

  for(j=0;j<nvars;j++)
  {
    c[j] = beta[j];
  }

bustout:;
  free(x[0]);
  free(x);
  free(y);
  free(d);
  free(rbar);
  free(thetab);
  free(xrow);
  free(beta);
  free(tol);
  free(rss);
  free(vorder);
  free(lindep);

  return(0);
}