    // initialize the "hessian" or whatever
    void hcalcs0(int n,dfsdmat & h,double& dmin)
    {
      int i,j;
      h.elem(1,1) = 1;
      for (i=2; i<=n; i++)
      {
        for ( j=1; j<i; j++)
        {
           h.elem(i,j)=0;
        }
        h.elem(i,i)=1;
      }
      dmin=h.elem(1,1);
      for ( i=2; i<=n; i++)
      {
         if(h.elem(i,i)<dmin)
            dmin=h.elem(i,i);
      }
    }

    void hcalcs1(int n,int n1,int np,int is,
      dfsdmat & h,dvector & g,dvector & w)
    {
      int i,j,i1;
      double z; 
      for (i=2; i<=n; i++)
      {
	 i1=i-1;
         z=-g.elem(i);
         double * pd=&(h.elem(i,1));
         double * pw=&(w.elem(1));
         for (j=1; j<=i1; j++)
         {
            z-=*pd++ * *pw++;
         }
         w.elem(i)=z;
      }
      w.elem(is+n)=w.elem(n)/h.elem(n,n);
      {
        dvector tmp(1,n);
        tmp.initialize();
        for (i=1; i<=n1; i++)
        {
          j=i;
          double * pd=&(h.elem(n-j+1,n-1));
          double qd=w.elem(is+np-j);
          double * pt=&(tmp(1));
          for (int ii=1; ii<=n1; ii++)
          {
            *pt++ +=*pd-- * qd;
          }
          w.elem(is+n-i)=w.elem(n-i)/h.elem(n-i,n-i)-tmp(i);
        }
      }
    }

    void hcalcs2(int n,int n1,int np,int iv,int iu,int ib, 
      dfsdmat & h,dvector & g,dvector & w,double sig,double dmin)
    {
      int i,j,i1;
      double z; 
      for (i=2;i<=n;i++)
      {
         i1=i-1;
         z=w.elem(iu+i);
         double * pd=&(h.elem(i,1));
         double * pw=&(w.elem(iv+1));
         for (j=1;j<=i1;j++)
         {
           z-=*pd++ * *pw++;
         }
         w.elem(iv+i)=z;
      }
      for (i=1;i<=n;i++)
      {
         z=h.elem(i,i)+sig*w.elem(iv+i)*w.elem(iv+i);
         if(z <= 0.0)
            z=dmin;
         if(z<dmin)
            dmin=z;
         h.elem(i,i)=z;
         w.elem(ib+i)=w.elem(iv+i)*sig/z;
         sig-=w.elem(ib+i)*w.elem(ib+i)*z;
       }
      for (j=2;j<=n;j++)
      {
         double * pd=&(h.elem(j,1));
         double * qd=&(w.elem(iu+j));
         double * rd=&(w.elem(iv+1));
         double * wd=&(w.elem(ib+1));
         for (i=1;i<j;i++)
         {
            *qd-=*pd * *rd++;
            *pd++ +=*wd++ * *qd;
         }
      }
    }
