einzelpunkt vorhersage

src/bk.f

@@ -20,8 +20,6 @@
      .              LDMPR,
      .              PHIPR,
      .              LDPHPR,
-c     .              PRINV,
-c     .              LDPRIV,
      .              MUWRK,
      .              PHIWRK,
      .              LDPHWK,

src/bk_point.c

+#include "bk_point.h" 
+
+void bk_point(double *xp,
+	     double *yp,
+	     double *zp,
+	     double *varp,
+	     int *dop,
+	      int* np,
+	     double *lon,
+	     double *lat,
+	     double *z,
+	     int *extrap,
+	     int *n,
+	     int *covtype,
+	     double *covpar,
+	     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 *lambda,
+	     double *lambd0,
+	     int *searchnb,
+	     int *ierr,
+	     int *glsmth){
+    
+    int ldcov=(*n), ldc0=(*n),ldphwk=(*ntrend),
+	ldfwrk=(*n), ldlmbd=(*n), ldkwrk=(*n), // ldcvbt=(*ntrend)*(*ntrend),
+	ldcvbt=(*ntrend), 
+	ldclup=(*n),ldcinv=(*n);
+    double *covmat, *c0vec, *muwrk, *phiwrk, *beta, errbta,
+	*dev, errdev, *covbta, *cvsrnb, *zsrnb,
+	*fwork, *fwrk2, *f0work, *dist, *kwork,
+	*rhswork, *fpwork, *fpfwork, *fpf0wrk, *chlup, 
+	*cminv, *work, *ferr, *berr,
+	*mu, // *lambda, 
+	cov0; // lambd0
+    int *indsnb, *indsnw, *indsrt, *ipiv, *ipvt, *iwork;
+
+
+    /* #if 0 */
+#ifndef TRANSIENT
+    covmat =Calloc((size_t)(*n)*(*n),double);
+    c0vec  =Calloc((size_t)(*n),double);  
+    muwrk  =Calloc((size_t)((*ntrend)*(*npr)),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); 
+    inddop =Calloc((size_t)(*np),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);
+    mu     =Calloc((size_t)(*ntrend),double);
+    //    lambda =Calloc((size_t)(*n),double); 
+#else
+    covmat =(double *) R_alloc((*n)*(*n),sizeof(double));
+    c0vec  =(double *) R_alloc((*n),sizeof(double));
+    muwrk  =(double *) R_alloc((*ntrend)*(*npr),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));
+    inddop =(int *) R_alloc((*np),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)); 
+    mu     =(double *) R_alloc((*ntrend),sizeof(double));
+    //    lambda =(double *) R_alloc((*n),sizeof(double));
+#endif
+
+    F77_CALL(bk)(xp,
+		 yp,
+		 dop,
+		 inddop, // 
+		 np,
+		     lon,
+		     lat,
+		     z,
+		     n,
+		     covtype,
+		     covpar,
+		     covmat,
+		     &ldcov,
+		     c0vec,
+		 n, //
+		     &cov0,
+		     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,
+		 ntrend, //
+		     dist,
+		     indsnb,
+		     indsnw,
+		     indsrt,
+		     kwork,
+		     &ldkwrk,
+		     rhswork,
+		     fpwork,
+		     fpfwork,
+		     fpf0wrk,
+		     chlup,
+		     &ldclup,
+		     cminv,
+		     &ldcinv,
+		     work,
+		     lwork,
+		     ipvt,
+		     ferr,
+		     berr,
+		     ipiv,
+		     iwork,
+		     mode,
+		     mu,
+		 ntrend, //
+		 zp,
+		 nap, //
+		     lambda,
+		     &ldlmbd,
+		     lambd0,
+		 varp,
+		 searchnb, //
+		 usesbbt, //
+		     ierr,
+		     glsmth);
+
+#ifndef TRANSIENT
+    //    Free(lambda);
+    Free(mu);
+    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(inddop);
+    Free(indsrt);
+    Free(indsnw);
+    Free(indsnb);
+    Free(dist);
+    Free(f0work);
+
+    Free(fwrk2);  // crash bei trend=0
+    Free(fwork);
+    Free(zsrnb);
+    Free(cvsrnb);
+    Free(dev);
+    Free(beta);
+    Free(covbta); // crash bei trend=1
+    Free(phiwrk);
+    Free(muwrk);
+    Free(c0vec);
+    Free(covmat);
+
+
+#else
+    /*
+    // free(lambda);
+    free(mu);
+    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(inddop);
+    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);
+    free(covmat);
+    */
+#endif
+    /* #endif */
+
+    printf("huhu\n");
+    
+}
+ 

src/bk_point.h

+#include <R.h>
+#include <Rinternals.h>
+
+// #define TRANSIENT
+
+void bk_point(double *xp,
+	     double *yp,
+	     double *zp,
+	     double *varp,
+	     int *dop,
+	      int *np,
+	     double *lon,
+	     double *lat,
+	     double *z,
+	     int *extrap,
+	     int *n,
+	     int *covtype,
+	     double *covpar,
+	     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 *lambda,
+	     double *lambd0,
+	     int *searchnb,
+	     int *ierr,
+	      int *glsmth);
+
+
+void F77_NAME(bk)(double *xp,
+		      double *yp,
+		  int *dop,
+		  int *inddop,
+		  int *np,
+		      double *lon,
+		      double *lat,
+		      double *z,
+		      int *n,
+		      int *covtype,
+		      double *covpar,
+		      double *covmat,
+		      int *ldcov,
+		      double *c0vec,
+		  int *ldc0,
+		      double *cov0,
+		      int *trend,
+		      int *ntrend,
+		      double *mupr,
+		      int *ldmpr,
+		      double *phipr,
+		      int *ldphpr,
+		      double *muwrk,
+		      double *phiwrk,
+		      int *ldphwk,
+		      double *lonpr,
+		      double *latpr,
+		      double *beta,
+		      double *errbta,
+		      double *covbta,
+		      int *ldcvbt,
+		      double *dev,
+		      double *errdev,
+		      double *cvsrnb,
+		      double *zsrnb,
+		      int *npr,
+		      int *typpr,
+		      double *rsearch,
+		      int *nsearch,
+		      int *nsmin,
+		      int *nsmax,
+		      double *fwork,
+		      double *fwrk2,
+		      int *ldfwrk,
+		      double *f0work,
+		  int *ldf0wk,
+		      double *dist,
+		      int *indsnb,
+		      int *indsnw,
+		      int *indsrt,
+		      double *kwork,
+		      int *ldkwrk,
+		      double *rhswork,
+		      double *fpwork,
+		      double *fpfwork,
+		      double *fpf0wrk,
+		      double *chlup,
+		      int *ldclup,
+		      double *cminv,
+		      int *ldcinv,
+		      double *work,
+		      int *lwork,
+		      int *ipvt,
+		      double *ferr,
+		      double *berr,
+		      int *ipiv,
+		      int *iwork,
+		      int *mode,
+		      double *mu,
+		  int *ldmu,
+		  double *zp,
+		  int *nap,
+		      double *lambda,
+		      int *ldlmbd,
+		      double *lambd0,
+		  double *varp,
+		      int *searchnb,
+		  int *usesbbt,
+		      int *ierr,
+		      int *glsmth);
+    
+