void
nb_kernel120_power6
     (int *         nri,        int *         iinr,     
      int *         jindex,     int *         jjnr,   
      int *         shift,      real *        shiftvec,
      real *        fshift,     int *         gid, 
      real *        pos,        real *        faction,
      real *        charge,     real *        facel,
      real *        krf,        real *        crf,  
      real *        Vc,         int *         type,   
      int *         ntype,      real *        vdwparam,
      real *        Vvdw,       real *        tabscale,
      real *        VFtab,      real *        invsqrta, 
      real *        dvda,       real *        gbtabscale,
      real *        GBtab,      int *         nthreads, 
      int *         count,      void *        mtx,
      int *         outeriter,  int *         inneriter,
      real *        work)
{
  F77_FUNC(pwr6kernel120,PWR6KERNEL120)
    (nri,iinr,jindex,jjnr,shift,shiftvec,fshift,gid,pos,faction,
     charge,facel,krf,crf,Vc,type,ntype,vdwparam,Vvdw,tabscale,
     VFtab,invsqrta,dvda,gbtabscale,GBtab,nthreads,count,mtx,
     outeriter,inneriter,work);
}

void
nb_kernel233nf_f77_double
     (int *         nri,        int           iinr[],     
      int           jindex[],   int           jjnr[],   
      int           shift[],    double        shiftvec[],
      double        fshift[],   int           gid[], 
      double        pos[],      double        faction[],
      double        charge[],   double *      facel,
      double *      krf,        double *      crf,  
      double        Vc[],       int           type[],   
      int *         ntype,      double        vdwparam[],
      double        Vvdw[],     double *      tabscale,
      double        VFtab[],    double        invsqrta[], 
      double        dvda[],     double *      gbtabscale,
      double        GBtab[],    int *         nthreads, 
      int *         count,      void *        mtx,
      int *         outeriter,  int *         inneriter,
      double *      work)
{
  F77_FUNC(f77dkernel233nf,F77DKERNEL233NF)
    (nri,iinr,jindex,jjnr,shift,shiftvec,fshift,gid,pos,faction,
     charge,facel,krf,crf,Vc,type,ntype,vdwparam,Vvdw,tabscale,
     VFtab,invsqrta,dvda,gbtabscale,GBtab,nthreads,count,mtx,
     outeriter,inneriter,work);
}

void init_mopac(t_QMrec *qm)
{
    /* initializes the mopac routines ans sets up the semiempirical
     * computation by calling moldat(). The inline mopac routines can
     * only perform gradient operations. If one would like to optimize a
     * structure or find a transition state at PM3 level, gaussian is
     * used instead.
     */
    char
    *keywords;

    snew(keywords, 240);

    if (!qm->bSH)  /* if rerun then grad should not be done! */
    {
        sprintf(keywords, "PRECISE GEO-OK CHARGE=%d GRAD MMOK ANALYT %s\n",
                qm->QMcharge,
                eQMmethod_names[qm->QMmethod]);
    }
    else
    {
        sprintf(keywords, "PRECISE GEO-OK CHARGE=%d SINGLET GRAD %s C.I.=(%d,%d) root=2 MECI \n",
                qm->QMcharge,
                eQMmethod_names[qm->QMmethod],
                qm->CASorbitals, qm->CASelectrons/2);
    }
    F77_FUNC(domldt, DOMLDT) (&qm->nrQMatoms, qm->atomicnumberQM, keywords);
    fprintf(stderr, "keywords are: %s\n", keywords);
    free(keywords);

} /* init_mopac */

/* Normally, SSTEVR is the LAPACK wrapper which calls one
 * of the eigenvalue methods. However, our code includes a
 * version of SSTEGR which is never than LAPACK 3.0 and can
 * handle requests for a subset of eigenvalues/vectors too,
 * and it should not need to call SSTEIN.
 * Just in case somebody has a faster version in their lapack
 * library we still call the driver routine, but in our own
 * case this is just a wrapper to sstegr.
 */
void
F77_FUNC(sstevr,SSTEVR)(const char *jobz, 
                        const char *range,
                        int *n,
                        float *d,
                        float *e,
                        float *vl, 
                        float *vu,
                        int *il, 
                        int *iu, 
                        float *abstol,
                        int *m,
                        float *w, 
                        float *z,
                        int *ldz,
                        int *isuppz, 
                        float *work, 
                        int *lwork, 
                        int *iwork,
                        int *liwork, 
                        int *info)
{
  F77_FUNC(sstegr,SSTEGR)(jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w,
	  z, ldz, isuppz, work, lwork, iwork, liwork, info);
  

    return;

}

void
nb_kernel400nf_f77_single
     (int *         nri,        int           iinr[],     
      int           jindex[],   int           jjnr[],   
      int           shift[],    float         shiftvec[],
      float         fshift[],   int           gid[], 
      float         pos[],      float         faction[],
      float         charge[],   float *       facel,
      float *       krf,        float *       crf,  
      float         Vc[],       int           type[],   
      int *         ntype,      float         vdwparam[],
      float         Vvdw[],     float *       tabscale,
      float         VFtab[],    float         invsqrta[], 
      float         dvda[],     float *       gbtabscale,
      float         GBtab[],    int *         nthreads, 
      int *         count,      void *        mtx,
      int *         outeriter,  int *         inneriter,
      float *       work)
{
  F77_FUNC(f77skernel400nf,F77SKERNEL400NF)
    (nri,iinr,jindex,jjnr,shift,shiftvec,fshift,gid,pos,faction,
     charge,facel,krf,crf,Vc,type,ntype,vdwparam,Vvdw,tabscale,
     VFtab,invsqrta,dvda,gbtabscale,GBtab,nthreads,count,mtx,
     outeriter,inneriter,work);
}

  ESymSolverStatus Ma27TSolverInterface::Backsolve(Index nrhs,
      double *rhs_vals)
  {
    DBG_START_METH("Ma27TSolverInterface::Backsolve",dbg_verbosity);
    IpData().TimingStats().LinearSystemBackSolve().Start();

    ipfint N=dim_;
    double* W = new double[maxfrt_];
    ipfint* IW1 = new ipfint[nsteps_];

    // For each right hand side, call MA27CD
    for(Index irhs=0; irhs<nrhs; irhs++) {
      if (DBG_VERBOSITY()>=2) {
        for (Index i=0; i<dim_; i++) {
          DBG_PRINT((2, "rhs[%5d] = %23.15e\n", i, rhs_vals[irhs*dim_+i]));
        }
      }

      F77_FUNC(ma27cd,MA27CD)(&N, a_, &la_, iw_, &liw_, W, &maxfrt_,
                              &rhs_vals[irhs*dim_], IW1, &nsteps_,
                              icntl_, cntl_);

      if (DBG_VERBOSITY()>=2) {
        for (Index i=0; i<dim_; i++) {
          DBG_PRINT((2, "sol[%5d] = %23.15e\n", i, rhs_vals[irhs*dim_+i]));
        }
      }
    }
    delete [] W;
    delete [] IW1;

    IpData().TimingStats().LinearSystemBackSolve().End();
    return SYMSOLVER_SUCCESS;
  }

static int f77_fopt(integer ndim, const doublereal *u, const integer *icp,
		    const doublereal *par, integer ijac,
		    doublereal *fs, doublereal *dfdu, doublereal *dfdp)
{
  F77_FUNC(fopt,FOPT)(&ndim, u, icp, par, &ijac, fs, dfdu, dfdp);
  return 0;
}

static int f77_bcnd(integer ndim, const doublereal *par, const integer *icp,
		    integer nbc, const doublereal *u0, const doublereal *u1,
		    integer ijac, doublereal *fb, doublereal *dbc)
{
  F77_FUNC(bcnd,BCND)(&ndim, par, icp, &nbc, u0, u1, fb, &ijac, dbc);
  return 0;
}

static int f77_func(integer ndim, const doublereal *u, const integer *icp,
		    const doublereal *par, integer ijac, doublereal *f,
		    doublereal *dfdu, doublereal *dfdp)
{
  F77_FUNC(func,FUNC)(&ndim, u, icp, par, &ijac, f, dfdu, dfdp);
  return 0;
}

real call_mopac_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
                   rvec f[], rvec fshift[])
{
    /* do the actual SH QMMM calculation using directly linked mopac
       subroutines */

    double /* always double as the MOPAC routines are always compiled in
              double precission! */
    *qmcrd = NULL, *qmchrg = NULL, *mmcrd = NULL, *mmchrg = NULL,
    *qmgrad, *mmgrad = NULL, energy;
    int
        i, j;
    real
        QMener = 0.0;

    snew(qmcrd, 3*(qm->nrQMatoms));
    snew(qmgrad, 3*(qm->nrQMatoms));
    /* copy the data from qr into the arrays that are going to be used
     * in the fortran routines of MOPAC
     */
    for (i = 0; i < qm->nrQMatoms; i++)
    {
        for (j = 0; j < DIM; j++)
        {
            qmcrd[3*i+j] = (double)qm->xQM[i][j]*10;
        }
    }
    if (mm->nrMMatoms)
    {
        /* later we will add the point charges here. There are some
         * conceptual problems with semi-empirical QM in combination with
         * point charges that we need to solve first....
         */
        gmx_fatal(FARGS, "At present only ONIOM is allowed in combination with MOPAC\n");
    }
    else
    {
        /* now compute the energy and the gradients.
         */
        snew(qmchrg, qm->nrQMatoms);

        F77_FUNC(domop, DOMOP) (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms,
                                mmchrg, mmcrd, qmgrad, mmgrad, &energy, qmchrg);
        /* add the gradients to the f[] array, and also to the fshift[].
         * the mopac gradients are in kCal/angstrom.
         */
        for (i = 0; i < qm->nrQMatoms; i++)
        {
            for (j = 0; j < DIM; j++)
            {
                f[i][j]      = (real)10*CAL2JOULE*qmgrad[3*i+j];
                fshift[i][j] = (real)10*CAL2JOULE*qmgrad[3*i+j];
            }
        }
        QMener = (real)CAL2JOULE*energy;
    }
    free(qmgrad);
    free(qmcrd);
    return (QMener);
} /* call_mopac_SH */

void
nb_kernel203_f77_double
     (int *         nri,        int *         iinr,     
      int *         jindex,     int *         jjnr,   
      int *         shift,      double *      shiftvec,
      double *      fshift,     int *         gid, 
      double *      pos,        double*       faction,
      double *      charge,     double*       facel,
      double *      krf,        double*       crf,  
      double *      Vc,         int *         type,   
      int *         ntype,      double *      vdwparam,
      double *      Vvdw,       double*       tabscale,
      double *      VFtab,      double*       invsqrta, 
      double *      dvda,       double*       gbtabscale,
      double *      GBtab,      int *         nthreads, 
      int *         count,      void *        mtx,
      int *         outeriter,  int *         inneriter,
      double *      work)
{
  F77_FUNC(f77dkernel203,F77DKERNEL203)
    (nri,iinr,jindex,jjnr,shift,shiftvec,fshift,gid,pos,faction,
     charge,facel,krf,crf,Vc,type,ntype,vdwparam,Vvdw,tabscale,
     VFtab,invsqrta,dvda,gbtabscale,GBtab,nthreads,count,mtx,
     outeriter,inneriter,work);
}

void CDM_FEA::Solve() //formulates and solves system!
{
	RemoveDisconnected(); 
	CalcDOF();
	CalcBonds();
	CalcStiffness(); //jmc: think it crashes here
	ApplyForces();

	if (DOF != 0){
	    iparm[2]  = -1; //sets to defualt system value...

		double ddum = 0; //Double dummy var
		int idum = 0; //Integer dummy var

		//msglvl = 0; //don't output info!
		phase = 13;
		//		PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &DOF, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error, dparm);
		//		F77_FUNC(PARDISO)(pt, &maxfct, &mnum, &mtype, &phase, &DOF, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error, dparm);
		F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase, &DOF, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error, dparm);

		//if (error != 0) std::cout << "Pardiso error! (" << error << ") - Phase 1\n";
		if (error == -1) std::cout << "Pardiso error: Input inconsistent\n";
		else if (error == -2) std::cout << "Pardiso error: Not enough memory\n";	
		else if (error == -3) std::cout << "Pardiso error: Reodering Problem\n";	
		else if (error == -4) std::cout << "Pardiso error: Zero pivot, numerical factorization or iterative refinement problem\n";
		else if (error == -10) std::cout << "Pardiso error: No License file Pardiso.lic found\n";
		else if (error == -11) std::cout << "Pardiso error: License is expired\n";
		else if (error == -12) std::cout << "Pardiso error: Wrong username or hostname\n";

		phase = -1; /* Release internal memory. */
		//		PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &DOF, &ddum, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error, dparm);
		//		F77_FUNC(PARDISO)(pt, &maxfct, &mnum, &mtype, &phase, &DOF, &ddum, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error, dparm);
		F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase, &DOF, &ddum, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error, dparm);
	}

	//CalcMaxDisps();
	FindMaxOverall(&Disp, x, MaxDisps);

	if (WantForces)
		CalcForces();

//	OutputMatrices();
	if (a != NULL) {delete [] a; a = NULL;}
	if (ia != NULL) {delete [] ia; ia = NULL;}
	if (ja != NULL) {delete [] ja; ja = NULL;}

}

int 
lapack_dgelqf (const int M, const int N, double *A, const int ldA,
               double *tau, double *work, const int lwork)
{
    int info = 0;
    F77_FUNC(dgelqf) (&M, &N, A, &ldA, tau, work, &lwork, &info);
    return info;
}

static int f77_icnd(integer ndim, const doublereal *par, const integer *icp,
		    integer nint, const doublereal *u, const doublereal *uold,
		    const doublereal *udot, const doublereal *upold, integer ijac,
		    doublereal *fi, doublereal *dint)
{
  F77_FUNC(icnd,ICND)(&ndim, par, icp, &nint, u, uold, udot, upold, fi, &ijac, dint);
  return 0;
}

void doit(int iter, struct problem *p)
{
     int i;

     for (i = 0; i < iter; ++i) {
	  F77_FUNC(fft4, FFT4)(&p->n[0], &m);
     }
}

void
F77_FUNC(dorml2,DORML2)(const char *side,
	const char *trans,
	int *m,
	int *n,
	int *k,
	double *a,
	int *lda,
	double *tau,
	double *c,
	int *ldc,
	double *work,
    int gmx_unused *info)
{
  const char xside=std::toupper(*side);
  const char xtrans=std::toupper(*trans);
  int i,i1,i2,i3,ni,mi,ic,jc;
  double aii;

  if(*m<=0 || *n<=0 || *k<=0)
    return;

  ic = jc = 0;

  if((xside=='L' && xtrans=='N') || (xside!='L' && xtrans!='N')) {
    i1 = 0;
    i2 = *k;
    i3 = 1;
  } else {
    i1 = *k-1;
    i2 = -1;
    i3 = -1;
  }
  
  if(xside=='L') {
    ni = *n;
    jc = 0;
  } else {
    mi = *m;
    ic = 0;
  }

  for(i=i1;i!=i2;i+=i3) {
    if(xside=='L') {
      mi = *m - i;
      ic = i;
    } else {
      ni = *n - i;
      jc = i;
    }
    aii = a[i*(*lda)+i];
    a[i*(*lda)+i] = 1.0;
    F77_FUNC(dlarf,DLARF)(side,&mi,&ni,&(a[i*(*lda)+i]),lda,tau+i,
	   &(c[jc*(*ldc)+ic]),ldc,work);
    a[i*(*lda)+i] = aii;
  }
  return;
}

ESymSolverStatus Ma57TSolverInterface::Backsolve(
    Index     nrhs,
    double    *rhs_vals)
{
    DBG_START_METH("Ma27TSolverInterface::Backsolve",dbg_verbosity);
    if (HaveIpData()) {
        IpData().TimingStats().LinearSystemBackSolve().Start();
    }

    ipfint  n      = dim_;
    ipfint  job    = 1;

    ipfint  nrhs_X = nrhs;
    ipfint  lrhs   = n;

    ipfint  lwork;
    double* work;

    lwork = n * nrhs;
    work = new double[lwork];

    // For each right hand side, call MA57CD
    // XXX MH: MA57 can do several RHSs; just do one solve...
    // AW: Ok is the following correct?
    if (DBG_VERBOSITY()>=2) {
        for (Index irhs=0; irhs<nrhs; irhs++) {
            for (Index i=0; i<dim_; i++) {
                DBG_PRINT((2, "rhs[%2d,%5d] = %23.15e\n", irhs, i, rhs_vals[irhs*dim_+i]));
            }
        }
    }

    F77_FUNC (ma57cd, MA57CD)
    (&job, &n, wd_fact_, &wd_lfact_, wd_ifact_, &wd_lifact_,
     &nrhs_X, rhs_vals, &lrhs,
     work, &lwork, wd_iwork_,
     wd_icntl_, wd_info_);

    if (wd_info_[0] != 0)
        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
                       "Error in MA57CD:  %d.\n", wd_info_[0]);

    if (DBG_VERBOSITY()>=2) {
        for (Index irhs=0; irhs<nrhs; irhs++) {
            for (Index i=0; i<dim_; i++) {
                DBG_PRINT((2, "sol[%2d,%5d] = %23.15e\n", irhs, i, rhs_vals[irhs*dim_+i]));
            }
        }
    }

    delete [] work;

    if (HaveIpData()) {
        IpData().TimingStats().LinearSystemBackSolve().End();
    }
    return SYMSOLVER_SUCCESS;
}

  ESymSolverStatus
  IterativeWsmpSolverInterface::InternalSymFact(
    const Index* ia,
    const Index* ja)
  {
    if (HaveIpData()) {
      IpData().TimingStats().LinearSystemSymbolicFactorization().Start();
    }

    // Call WISMP for ordering and symbolic factorization
    ipfint N = dim_;
    IPARM_[1] = 1; // ordering
    IPARM_[2] = 1; // symbolic factorization
    ipfint idmy;
    double ddmy;
    Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
                   "Calling WISMP-1-1 for symbolic analysis at cpu time %10.3f (wall %10.3f).\n", CpuTime(), WallclockTime());
    F77_FUNC(wismp,WISMP)(&N, ia, ja, a_, &ddmy, &idmy, &ddmy, &idmy, &idmy,
                          &ddmy, &ddmy, IPARM_, DPARM_);
    Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
                   "Done with WISMP-1-1 for symbolic analysis at cpu time %10.3f (wall %10.3f).\n", CpuTime(), WallclockTime());

    Index ierror = IPARM_[63];
    if (ierror!=0) {
      if (ierror==-102) {
        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
                       "Error: WISMP is not able to allocate sufficient amount of memory during ordering/symbolic factorization.\n");
      }
      else if (ierror>0) {
        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                       "Matrix appears to be singular (with ierror = %d).\n",
                       ierror);
        if (HaveIpData()) {
          IpData().TimingStats().LinearSystemSymbolicFactorization().End();
        }
        return SYMSOLVER_SINGULAR;
      }
      else {
        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
                       "Error in WISMP during ordering/symbolic factorization phase.\n     Error code is %d.\n", ierror);
      }
      if (HaveIpData()) {
        IpData().TimingStats().LinearSystemSymbolicFactorization().End();
      }
      return SYMSOLVER_FATAL_ERROR;
    }
    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                   "Predicted memory usage for WISMP after symbolic factorization IPARM(23)= %d.\n",
                   IPARM_[22]);

    if (HaveIpData()) {
      IpData().TimingStats().LinearSystemSymbolicFactorization().End();
    }

    return SYMSOLVER_SUCCESS;
  }

void IpLapackDsyev(bool compute_eigenvectors, Index ndim, Number *a,
                   Index lda, Number *w, Index& info)
{
#ifdef COIN_HAS_LAPACK
    ipfint N=ndim, LDA=lda, INFO;

    char JOBZ;
    if (compute_eigenvectors) {
        JOBZ = 'V';
    }
    else {
        JOBZ = 'N';
    }
    char UPLO = 'L';

    // First we find out how large LWORK should be
    ipfint LWORK = -1;
    double WORK_PROBE;
    F77_FUNC(dsyev,DSYEV)(&JOBZ, &UPLO, &N, a, &LDA, w,
                          &WORK_PROBE, &LWORK, &INFO, 1, 1);
    DBG_ASSERT(INFO==0);

    LWORK = (ipfint) WORK_PROBE;
    DBG_ASSERT(LWORK>0);

    double* WORK = new double[LWORK];
    for (Index i=0; i<LWORK; i++) {
        WORK[i] = i;
    }
    F77_FUNC(dsyev,DSYEV)(&JOBZ, &UPLO, &N, a, &LDA, w,
                          WORK, &LWORK, &INFO, 1, 1);

    DBG_ASSERT(INFO>=0);
    info = INFO;

    delete [] WORK;
#else

    std::string msg = "Ipopt has been compiled without LAPACK routine DSYEV, but options are chosen that require this dependency.  Abort.";
    THROW_EXCEPTION(LAPACK_NOT_INCLUDED, msg);
#endif

}

int 
lapack_dormlq (const enum BLAS_SIDE side, const enum BLAS_TRANSPOSE trans,
               const int M, const int N, const int K, const double *A,
               const int ldA, double *tau, double *C, const int ldC,
               double *Work, const int ldWork)
{
    int info = 0;
    F77_FUNC(dormlq) (SIDE(side), TRANS(trans), &M, &N, &K, A, &ldA, tau, C,
                      &ldC, Work, &ldWork, &info);
    return info;
}