自学教程：C++ F77_FUNC函数代码示例

voidnb_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);}

voidnb_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. */voidF77_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;}

voidnb_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 */

voidnb_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);     }}

voidF77_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;}

/* Uses factorization to solve. */voidClpCholeskyWssmpKKT::solveKKT (double * region1, double * region2, const double * diagonal,                               double diagonalScaleFactor){    int numberRowsModel = model_->numberRows();    int numberColumns = model_->numberColumns();    int numberTotal = numberColumns + numberRowsModel;    double * array = new double [numberRows_];    CoinMemcpyN(region1, numberTotal, array);    CoinMemcpyN(region2, numberRowsModel, array + numberTotal);    int i1 = 1;    int i0 = 0;    integerParameters_[1] = 4;    integerParameters_[2] = 4;#if 0    integerParameters_[5] = 3;    doubleParameters_[5] = 1.0e-10;    integerParameters_[6] = 6;#endif    F77_FUNC(wssmp,WSSMP)(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,                          NULL, permute_, permuteInverse_, array, &numberRows_, &i1,                          NULL, &i0, NULL, integerParameters_, doubleParameters_);#if 0    int iRow;    for (iRow = 0; iRow < numberTotal; iRow++) {        if (rowsDropped_[iRow] && fabs(array[iRow]) > 1.0e-8) {            printf("row region1 %d dropped %g/n", iRow, array[iRow]);        }    }    for (; iRow < numberRows_; iRow++) {        if (rowsDropped_[iRow] && fabs(array[iRow]) > 1.0e-8) {            printf("row region2 %d dropped %g/n", iRow, array[iRow]);        }    }#endif    CoinMemcpyN(array + numberTotal, numberRowsModel, region2);#if 1    CoinMemcpyN(array, numberTotal, region1);#else    multiplyAdd(region2, numberRowsModel, -1.0, array + numberColumns, 0.0);    CoinZeroN(array, numberColumns);    model_->clpMatrix()->transposeTimes(1.0, region2, array);    for (int iColumn = 0; iColumn < numberTotal; iColumn++)        region1[iColumn] = diagonal[iColumn] * (array[iColumn] - region1[iColumn]);#endif    delete [] array;#if 0    if (integerParameters_[5]) {        std::cout << integerParameters_[5] << " refinements ";    }    std::cout << doubleParameters_[6] << std::endl;#endif}

voidF77_FUNC(slarnv,SLARNV)(int *idist, 	int *iseed, 	int *n, 	float *x){    int i__1, i__2, i__3;    int i__;    float u[128];    int il, iv, il2;    --x;    --iseed;    i__1 = *n;    for (iv = 1; iv <= i__1; iv += 64) {	i__2 = 64, i__3 = *n - iv + 1;	il = (i__2<i__3) ? i__2 : i__3;	if (*idist == 3) {	    il2 = il << 1;	} else {	    il2 = il;	}	F77_FUNC(slaruv,SLARUV)(&iseed[1], &il2, u);	if (*idist == 1) {	    i__2 = il;	    for (i__ = 1; i__ <= i__2; ++i__) {		x[iv + i__ - 1] = u[i__ - 1];	    }	} else if (*idist == 2) {	    i__2 = il;	    for (i__ = 1; i__ <= i__2; ++i__) {		x[iv + i__ - 1] = u[i__ - 1] * 2. - 1.;	    }	} else if (*idist == 3) {	    i__2 = il;	    for (i__ = 1; i__ <= i__2; ++i__) {		x[iv + i__ - 1] = sqrt(log(u[(i__ << 1) - 2]) * -2.) * 		  cos(u[(i__ << 1) - 1] * 		      (float)6.2831853071795864769252867663);	    }	}    }    return;}

  bool Ma27TSolverInterface::InitializeImpl(const OptionsList& options,      const std::string& prefix)  {    options.GetNumericValue("ma27_pivtol", pivtol_, prefix);    if (options.GetNumericValue("ma27_pivtolmax", pivtolmax_, prefix)) {      ASSERT_EXCEPTION(pivtolmax_>=pivtol_, OPTION_INVALID,                       "Option /"ma27_pivtolmax/": This value must be between "                       "ma27_pivtol and 1.");    }    else {      pivtolmax_ = Max(pivtolmax_, pivtol_);    }    options.GetNumericValue("ma27_liw_init_factor", liw_init_factor_, prefix);    options.GetNumericValue("ma27_la_init_factor", la_init_factor_, prefix);    options.GetNumericValue("ma27_meminc_factor", meminc_factor_, prefix);    options.GetBoolValue("ma27_skip_inertia_check",                         skip_inertia_check_, prefix);    options.GetBoolValue("ma27_ignore_singularity",                         ignore_singularity_, prefix);    // The following option is registered by OrigIpoptNLP    options.GetBoolValue("warm_start_same_structure",                         warm_start_same_structure_, prefix);    /* Set the default options for MA27 */    F77_FUNC(ma27id,MA27ID)(icntl_, cntl_);#if COIN_IPOPT_VERBOSITY == 0    icntl_[0] = 0;       // Suppress error messages    icntl_[1] = 0;       // Suppress diagnostic messages#endif    // Reset all private data    initialized_=false;    pivtol_changed_ = false;    refactorize_ = false;    la_increase_=false;    liw_increase_=false;    if (!warm_start_same_structure_) {      dim_=0;      nonzeros_=0;    }    else {      ASSERT_EXCEPTION(dim_>0 && nonzeros_>0, INVALID_WARMSTART,                       "Ma27TSolverInterface called with warm_start_same_structure, but the problem is solved for the first time.");    }    return true;  }

void IpLapackDgetrf(Index ndim, Number *a, Index *ipiv, Index lda, Index& info){#ifdef COIN_HAS_LAPACK    ipfint M=ndim, N=ndim, LDA=lda, INFO;    F77_FUNC(dgetrf,DGETRF)(&M, &N, a, &LDA, ipiv, &INFO);    info = INFO;#else    std::string msg = "Ipopt has been compiled without LAPACK routine DPOTRF, but options are chosen that require this dependency.  Abort.";    THROW_EXCEPTION(LAPACK_NOT_INCLUDED, msg);#endif}

CDM_FEA::CDM_FEA(void){	Indi = NULL;	Indj = NULL;	BondDir = NULL;	IndextoDOF = NULL;	FixedList = NULL;	F = NULL;	e = NULL;	MaxForces = NULL;	MaxDisps = NULL;	MaxReactions = NULL;	MaxStrains = NULL;	MaxSE = NULL;	//Pardiso Params!	//...that change	a = NULL;	ja = NULL;	ia = NULL;	b = NULL;	x = NULL;	DOF = -1;	//...that don't change	mtype = 2; // Real symmetric matrix 	nrhs = 1; // Number of right hand sides.	maxfct = 1; //Maximum number of numerical factorizations.	mnum = 1; //Which factorization to use.	msglvl = 1; //Print statistical information	error = 0; //Initialize error flag	int solver = 0; //use default (non-iterative) Pardiso solver	//	PARDISOINIT(pt, &mtype, &solver, iparm, dparm, &error); //initialize pardiso	//	F77_FUNC(PARDISOINIT)(pt, &mtype, &solver, iparm, dparm, &error); //initialize pardiso	F77_FUNC(pardisoinit)(pt, &mtype, &solver, iparm, dparm, &error); //initialize pardiso	pObj = NULL;		Element_type = FRAME; //the type of element! (default)	DOFperBlock = 0; //the dimension of each metablock	ELperDBlock = 0; //the number of elements per metablock 	ELperOBlock = 0; //the number of elements per metablock 	ResetFEA();}

  ESymSolverStatus WsmpSolverInterface::Solve(    const Index* ia,    const Index* ja,    Index nrhs,    double *rhs_vals)  {    DBG_START_METH("WsmpSolverInterface::Solve",dbg_verbosity);    IpData().TimingStats().LinearSystemBackSolve().Start();    // Call WSMP to solve for some right hand sides (including    // iterative refinement)    // ToDo: Make iterative refinement an option?    ipfint N = dim_;    ipfint LDB = dim_;    ipfint NRHS = nrhs;    ipfint NAUX = 0;    IPARM_[1] = 4; // Forward and Backward Elimintation    IPARM_[2] = 5; // Iterative refinement    IPARM_[5] = 1;    DPARM_[5] = 1e-12;    ipfint idmy;    double ddmy;    F77_FUNC(wssmp,WSSMP)(&N, ia, ja, a_, &ddmy, PERM_, INVP_,                          rhs_vals, &LDB, &NRHS, &ddmy, &NAUX,                          &idmy, IPARM_, DPARM_);    IpData().TimingStats().LinearSystemBackSolve().End();    Index ierror = IPARM_[63];    if (ierror!=0) {      if (ierror==-102) {        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,                       "Error: WSMP is not able to allocate sufficient amount of memory during ordering/symbolic factorization./n");      }      else {        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,                       "Error in WSMP during ordering/symbolic factorization phase./n     Error code is %d./n", ierror);      }      return SYMSOLVER_FATAL_ERROR;    }    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,                   "Number of iterative refinement steps in WSSMP: %d/n",                   IPARM_[5]);    return SYMSOLVER_SUCCESS;  }

/* Interface to FORTRAN routine DPOTRS. */void IpLapackDpotrs(Index ndim, Index nrhs, const Number *a, Index lda,                    Number *b, Index ldb){#ifdef COIN_HAS_LAPACK    ipfint N=ndim, NRHS=nrhs, LDA=lda, LDB=ldb, INFO;    char uplo = 'L';    F77_FUNC(dpotrs,DPOTRS)(&uplo, &N, &NRHS, a, &LDA, b, &LDB, &INFO, 1);    DBG_ASSERT(INFO==0);#else    std::string msg = "Ipopt has been compiled without LAPACK routine DPOTRS, but options are chosen that require this dependency.  Abort.";    THROW_EXCEPTION(LAPACK_NOT_INCLUDED, msg);#endif}

// solves system of symmertic positive definite Ax=b and store the solution// in x.void ColaModel::lapack_solve(double ** A, double * b, double * x, int dim) {  char uplo = 'L';  int num_rhs = 1;  int info;  // copy b to x  std::copy(b, b+dim, x);  int size = (dim*dim+dim)/2;  double * A_lower = new double[size];  for (int i=0; i<dim; ++i) {    std::copy(A[i], A[i]+i+1, A_lower+(i*i+i)/2);  }  F77_FUNC (dposv, DPOSV) (&uplo, &dim, &num_rhs, A_lower, &dim, x, &dim, &info);  if (info!=0) {    std::cerr << "Lapack dposv function failed." << std::endl;    throw std::exception();  }  delete[] A_lower;}