Revision 8ad4e6ec8732718f8759ef4472695e9bfb5ead73 authored by Douglas Nychka on 23 April 2019, 11:50:03 UTC, committed by cran-robot on 23 April 2019, 11:50:03 UTC
1 parent 9d93234
multebC.c
/*
c**** # fields is a package for analysis of spatial data written for
c**** # the R software environment .
c**** # Copyright (C) 2018
c**** # University Corporation for Atmospheric Research (UCAR)
c**** # Contact: Douglas Nychka, nychka@ucar.edu,
c**** # National Center for Atmospheric Research, PO Box 3000, Boulder, CO 80307-3000
c**** #
c**** # This program is free software; you can redistribute it and/or modify
c**** # it under the terms of the GNU General Public License as published by
c**** # the Free Software Foundation; either version 2 of the License, or
c**** # (at your option) any later version.
c**** # This program is distributed in the hope that it will be useful,
c**** # but WITHOUT ANY WARRANTY; without even the implied warranty of
c**** # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c**** # GNU General Public License for more details.
*/
#include <R.h>
#include <Rinternals.h>
#include <R_ext/Arith.h>
#include <Rmath.h>
#include <float.h>
#include <math.h>
#include <R_ext/BLAS.h>
#include <R_ext/Print.h>
#include <unistd.h>
SEXP multebC(SEXP nd, SEXP x1, SEXP n1, SEXP x2, SEXP n2, SEXP par, SEXP c, SEXP work) {
// NOTE: in original multeb function in fields, "h" was passed but not used and returned as an answer.
// h is not an argument to this function, but "ans" is equivalent to h and is allocated in C and returned.
int Ind, In1, In2, r1, r2, d;
double sum;
double *Px1, *Px2, *Pc, *Pwork, *cans;
SEXP ans;
void expfnC( SEXP, SEXP, SEXP);
// cast R variables to C variables
Ind = INTEGER(nd)[0];
In1 = INTEGER(n1)[0];
In2 = INTEGER(n2)[0];
Px1 = REAL(x1);
Px2 = REAL(x2);
Pc = REAL(c);
Pwork = REAL(work);
// const int tmp = 1;
const int cn2 = In2;
// allocate answer vector (corresponds to h in fields' multeb)
ans = PROTECT(allocVector(REALSXP, In1));
cans = REAL(ans);
// work aray must be dimensioned to size n2
// outer most loop over columns of x1 and x2 should reduce paging
for(r1 = 0; r1 < In1; r1++) {
// evaluate all basis functions at x1[r2,.]
for(r2 = 0; r2 < In2; r2++) {
// zero out sum accumulator
sum=0.0;
for(d = 0; d < Ind; d++) {
sum += pow(fabs(Px1[In1*d+r1] - Px2[In2*d+r2]), 2.0);
}
Pwork[r2]=sum;
}
// evaluate squared distances with basis functions.
expfnC(n2, work, par);
// now the dot product you have all been waiting for!
sum=0.0;
for(d = 0; d < cn2; d++) {
sum += Pwork[d]*Pc[d] ;
}
// cans[r1] = ddot_(&cn2, Pwork, &tmp, Pc, &tmp);
cans[r1] = sum;
}
UNPROTECT(1);
return(ans);
}
Computing file changes ...