#include "bk_grid.h" 

void bk_grid__(double *xsw,
	     double *ysw,
	     double *xne,
	     double *yne,
	     double *angle,
	     int *nx,
	     int *ny,
	     double *dx,
	     double *dy,
	     double *xg,
	     double *yg,
	     double *zg,
	     double *varg,
	     int *dog,
	     double *lon,
	     double *lat,
	     double *z,
	     int *extrap,
	     int *n,
	     int *covtype,
	     double *covpar,
	     double *covmat,
	     int *ldcov,
	     int *extcov,
	     int *trend,
	     int *ntrend,
	     double *mupr,
	     int *ldmpr,
	     double *phipr,
	     int *ldphpr,
	     double *lonpr,
	     double *latpr,
	     int *npr,
	     int *typpr,
	     double *rsearch,
	     int *nsearch,
	     int *nsmin,
	     int *nsmax,
	     int *lwork,
	     int *mode,
	       double *mu,
	     double *lambda,
	     double *lambd0,
	     int *searchnb,
	     int *ierr,
	       int *retlm,
	     int *glsmth){
/* simple Fortran wrapper */
bk_grid(xsw,
	     ysw,
	     xne,
	     yne,
	     angle,
	     nx,
	     ny,
	     dx,
	     dy,
	     xg,
	     yg,
	     zg,
	     varg,
	     dog,
	     lon,
	     lat,
	     z,
	     extrap,
	     n,
	     covtype,
	     covpar,
	     covmat,
	     ldcov,
	     extcov,
	     trend,
	     ntrend,
	     mupr,
	     ldmpr,
	     phipr,
	     ldphpr,
	     lonpr,
	     latpr,
	     npr,
	     typpr,
	     rsearch,
	     nsearch,
	     nsmin,
	     nsmax,
	     lwork,
	     mode,
	mu,
	     lambda,
	     lambd0,
	     searchnb,
	     ierr,
	retlm,
	     glsmth);
}

void bk_grid(double *xsw,
	     double *ysw,
	     double *xne,
	     double *yne,
	     double *angle,
	     int *nx,
	     int *ny,
	     double *dx,
	     double *dy,
	     double *xg,
	     double *yg,
	     double *zg,
	     double *varg,
	     int *dog,
	     double *lon,
	     double *lat,
	     double *z,
	     int *extrap,
	     int *n,
	     int *covtype,
	     double *covpar,
	     double *covmat,
	     int *ldcov,
	     int *extcov,
	     int *trend,
	     int *ntrend,
	     double *mupr,
	     int *ldmpr,
	     double *phipr,
	     int *ldphpr,
	     double *lonpr,
	     double *latpr,
	     int *npr,
	     int *typpr,
	     double *rsearch,
	     int *nsearch,
	     int *nsmin,
	     int *nsmax,
	     int *lwork,
	     int *mode,
	     double *mu,
	     double *lambda,
	     double *lambd0,
	     int *searchnb,
	     int *ierr,
	     int *retlm,
	     int *glsmth){
    
    int nz=(*nx)*(*ny), ldc0=(*n),ldphwk=(*ntrend),
	ldfwrk=(*n), ldlmbd=(*n), ldkwrk=(*n), 
	ldcvbt=(*ntrend), 
	ldclup=(*n),ldcinv=(*n),ldzg=(*nx);
    double *c0vec, *muwrk, *phiwrk, *beta, errbta,
	*dev, errdev, *covbta, *cvsrnb, *zsrnb,
	*fwork, *fwrk2, *f0work, *dist, *kwork,
	*rhswork, *fpwork, *fpfwork, *fpf0wrk, *chlup, 
	*cminv, *work, *ferr, *berr,
	cov0;
    int *indsnb, *indsnw, *indsrt, *ipiv, *ipvt, *iwork;


    /* #if 0 */
#ifndef TRANSIENT
    c0vec  =Calloc((size_t)(*n),double);  
    muwrk  =Calloc((size_t)(*ntrend),double); 
    phiwrk =Calloc((size_t)((*ntrend)*(*ntrend)),double); 
    covbta =Calloc((size_t)((*ntrend)*(*ntrend)),double); 
    beta   =Calloc((size_t)(*ntrend),double); 
    dev    =Calloc((size_t)(*n),double); 
    cvsrnb =Calloc((size_t)(*n)*(*n),double);
    zsrnb  =Calloc((size_t)(*n),double);
    fwork  =Calloc((size_t)(*n)*(*ntrend),double);    
    fwrk2 =Calloc((size_t)(*n)*(*ntrend),double);    
    f0work =Calloc((size_t)(*ntrend),double);    
    dist   =Calloc((size_t)(*n),double);    
    indsnb =Calloc((size_t)(*n),int);    
    indsnw =Calloc((size_t)(*n),int);    
    indsrt =Calloc((size_t)(*n),int); 
    kwork  =Calloc((size_t)(*n)*(*n),double); 
    rhswork=Calloc((size_t)(*n),double);
    fpwork =Calloc((size_t)(*n)*(*ntrend),double);    
    fpfwork=Calloc((size_t)(*n)*(*n),double);    
    fpf0wrk=Calloc((size_t)(*n),double);    
    chlup  =Calloc((size_t)(*n)*(*n),double);    
    cminv  =Calloc((size_t)(*n)*(*n),double);    
    work   =Calloc((size_t)(*lwork),double);    
    ipvt   =Calloc((size_t)(*n),int);
    ipiv   =Calloc((size_t)(*n+*ntrend),int);
    ferr   =Calloc((size_t)(*n),double);
    berr   =Calloc((size_t)(*n),double);
    iwork  =Calloc((size_t)(3*(*n)),int);
#else
    c0vec  =(double *) R_alloc((*n),sizeof(double));
    muwrk  =(double *) R_alloc((*ntrend),sizeof(double));
    phiwrk =(double *) R_alloc((*ntrend)*(*ntrend),sizeof(double)); 
    covbta =(double *) R_alloc((*ntrend)*(*ntrend),sizeof(double)); 
    beta   =(double *) R_alloc((*ntrend),sizeof(double));
    dev    =(double *) R_alloc((*n),sizeof(double));
    cvsrnb =(double *) R_alloc((*n)*(*n),sizeof(double));
    zsrnb  =(double *) R_alloc((*n),sizeof(double));
    fwork  =(double *) R_alloc((*n)*(*ntrend),sizeof(double)); 
    fwrk2 =(double *) R_alloc((*n)*(*ntrend),sizeof(double)); 
    f0work =(double *) R_alloc((*ntrend),sizeof(double));    
    dist   =(double *) R_alloc((*n),sizeof(double));
    indsnb =(int *) R_alloc((*n),sizeof(int));
    indsnw =(int *) R_alloc((*n),sizeof(int));
    indsrt =(int *) R_alloc((*n),sizeof(int));
    kwork  =(double *) R_alloc((*n+*ntrend)*(*n+*ntrend),sizeof(double)); 
    rhswork=(double *) R_alloc((*n+*ntrend),sizeof(double)); 
    fpwork =(double *) R_alloc((*n)*(*ntrend),sizeof(double)); 
    fpfwork=(double *) R_alloc((*n)*(*n),sizeof(double)); 
    fpf0wrk=(double *) R_alloc((*n),sizeof(double)); 
    chlup  =(double *) R_alloc((*n)*(*n),sizeof(double)); 
    cminv  =(double *) R_alloc((*n)*(*n),sizeof(double)); 
    work   =(double *) R_alloc((*lwork),sizeof(double)); 
    ipvt   =(int *) R_alloc((*n),sizeof(int)); 
    ipiv   =(int *) R_alloc((*n+*ntrend),sizeof(int)); 
    ferr   =(double *) R_alloc((*n),sizeof(double)); 
    berr   =(double *) R_alloc((*n),sizeof(double)); 
    iwork  =(int *) R_alloc(3*(*n),sizeof(int)); 
#endif

    F77_CALL(bkgrid)(xsw,
		     ysw,
		     xne,
		     yne,
		     angle,
		     nx,
		     ny,
		     &nz,
		     dx,
		     dy,
		     xg,
		     yg,
		     zg,
		     &ldzg,
		     varg,
		     dog,
		     lon,
		     lat,
		     z,
		     extrap,
		     n,
		     covtype,
		     covpar,
		     c0vec,
		     &cov0,
		     covmat,
		     ldcov,
		     extcov,
		     trend,
		     ntrend,
		     mupr,
		     ldmpr,
		     phipr,
		     ldphpr,
		     muwrk,
		     phiwrk,
		     &ldphwk,
		     lonpr,
		     latpr,
		     beta,
		     &errbta,
		     covbta,
		     &ldcvbt,
		     dev,
		     &errdev,
		     cvsrnb,
		     zsrnb,
		     npr,
		     typpr,
		     rsearch,
		     nsearch,
		     nsmin,
		     nsmax,
		     fwork,
		     fwrk2,
		     &ldfwrk,
		     f0work,
		     dist,
		     indsnb,
		     indsnw,
		     indsrt,
		     kwork,
		     &ldkwrk,
		     rhswork,
		     fpwork,
		     fpfwork,
		     fpf0wrk,
		     chlup,
		     &ldclup,
		     cminv,
		     &ldcinv,
		     work,
		     lwork,
		     ipvt,
		     ferr,
		     berr,
		     ipiv,
		     iwork,
		     mode,
		     mu,
		     lambda,
		     &ldlmbd,
		     lambd0,
		     searchnb,
		     ierr,
		     retlm,
		     glsmth);

#ifndef TRANSIENT
    Free(iwork);
    Free(berr);
    Free(ferr);
    Free(ipiv);
    Free(ipvt);
    Free(work);
    Free(cminv);
    Free(chlup);
    Free(fpf0wrk);
    Free(fpfwork);
    Free(fpwork);
    Free(rhswork);
    Free(kwork);
    Free(indsrt);
    Free(indsnw);
    Free(indsnb);
    Free(dist);
    Free(f0work);

    Free(fwrk2);  
    Free(fwork);
    Free(zsrnb);
    Free(cvsrnb);
    Free(dev);
    Free(beta);
    Free(covbta); 
    Free(phiwrk);
    Free(muwrk);
    Free(c0vec);


#else
    /*
    free(iwork);
    free(berr);
    free(ferr);
    free(ipiv);
    free(ipvt);
    free(work);
    free(cminv);
    free(chlup);
    free(fpf0wrk);
    free(fpfwork);
    free(fpwork);
    free(rhswork);
    free(kwork);
    free(indsrt);
    free(indsnw);
    free(indsnb);
    free(dist);
    free(f0work);
    free(fwrk2);
    free(fwork);
    free(zsrnb);
    free(cvsrnb);
    free(dev);
    free(beta);
    free(covbta);
    free(phiwrk);
    free(muwrk);
    free(c0vec);
    */
#endif
    /* #endif */

}