//Ｒ用 Ｃ
//#include <stdio.h>
//#include <math.h>

#include <R.h>
#include <Rdefines.h>
#include <R_ext/Parse.h>

void convolve(double *a, int *na, double *b, int *nb, double *ab)
{
    int i, j, nab = *na + *nb - 1;
    for(i = 0; i < nab; i++)
        ab[i] = 0.0;
    for(i = 0; i < *na; i++)
        for(j = 0; j < *nb; j++)
            ab[i + j] += a[i] * b[j];
}

SEXP convolve2(SEXP a, SEXP b)
{
  int i, j, na, nb, nab;
  double *xa, *xb, *xab;
  SEXP ab;
  PROTECT(a = AS_NUMERIC(a));
  PROTECT(b = AS_NUMERIC(b));
  na = LENGTH(a); nb = LENGTH(b); nab = na + nb - 1;
  PROTECT(ab = NEW_NUMERIC(nab));
  xa = NUMERIC_POINTER(a); xb = NUMERIC_POINTER(b);
  xab = NUMERIC_POINTER(ab);
  for(i = 0; i < nab; i++) xab[i] = 0.0;
  for(i = 0; i < na; i++)
    for(j = 0; j < nb; j++) xab[i + j] += xa[i] * xb[j];
      UNPROTECT(3);
  return(ab);
}

SEXP test001(SEXP a){
  int n, m;
  n = INTEGER(GET_DIM(a))[0];
  m = INTEGER(GET_DIM(a))[1];
  printf("nrow of a : %d\n", n);
  printf("ncol of a : %d\n", m);
  return a;
}

SEXP Rdgesv01(SEXP A, SEXP b, SEXP info)
{
  int i;
  long int n = INTEGER(GET_DIM(A))[0];
  long int nrhs=1;
  long int incx=1;
  long int incy=1;
  long int piv[n];
  
  SEXP tmpA, tmpb;
  
  PROTECT(tmpA = allocMatrix(REALSXP, n, n));
  PROTECT(tmpb = allocVector(REALSXP, n));
  long int nn = n*n;
  dcopy_(&nn, REAL(A), &incx, REAL(tmpA), &incy);
  dcopy_(&n, REAL(b), &incx, REAL(tmpb), &incy);
  //A[0]=1.; A[1]=3.; A[2]=1.;
  //A[3]=1.; A[4]=1.; A[5]=-2.;
  //A[6]=1.; A[7]=-3.;A[8]=-5.;
  //x[0]=3.; x[1]=1.; x[2]=-6.;
  //y[0]=0.; y[1]=0.; y[2]=0.;

  printf("N = %ld\n", n);
  dgesv_(&n, &nrhs, REAL(tmpA), &n, piv, REAL(tmpb), &n, INTEGER(info));
  for(i=0; i<n; ++i) printf("%lf\n", REAL(b)[i]);
  UNPROTECT(2);
  return(tmpb);
}

SEXP eval_fun01(SEXP a, SEXP rho)
{
  SEXP s, t;
  SEXP ret;
  PROTECT(t = s = allocList(3));
  // Setting the type to LANGSXP makes this a call which can be evaluated
  SET_TYPEOF(s, LANGSXP);
  // set a symbol (pointing to the function to be called)
  SETCAR(t, install("print")); t = CDR(t);
  // set the first unnamed
  SETCAR(t, a); t = CDR(t);
  // set the second named
  SETCAR(t, allocVector(INTSXP, 1));
  INTEGER(CAR(t))[0] = 3;
  SET_TAG(t, install("digits"));
  PROTECT(ret = eval(s, rho));
  //PrintValue(ret);
  UNPROTECT(2);
  return R_NilValue;
}

SEXP eval_fun02(SEXP a, SEXP rho)
{
  SEXP call;
  SEXP ret;
  
  PROTECT(call = lang2(install("print"), a));
  eval(call, rho);
  UNPROTECT(1);
  return R_NilValue;
}

SEXP eval_fun03(SEXP fr, SEXP env)
{
  SEXP ans;
  
  if(!isEnvironment(env)) error(" 'env' should be an environment");
  PROTECT(ans = eval(fr, env));
  UNPROTECT(1);
  return(ans);
}

SEXP eval_fun04(SEXP fr, SEXP x, SEXP env)
{
  SEXP ans;
  
  if(!isEnvironment(env)) error(" 'env' should be an environment");
  defineVar(install("x"), x, env);
  PROTECT(ans = eval(fr, env));
  UNPROTECT(1);
  return(ans);
}

SEXP eval_fun05(SEXP a, SEXP env)
{
  SEXP cv, call, call2;
  SEXP ans;
  int i;
  ParseStatus status;
  char *cmd1 = "print(a)";  //なんでもOK．文字列を準備
  PROTECT(cv = allocVector(STRSXP, 1));  //文字列用のオブジェクトを用意
  SET_STRING_ELT(cv, 0, mkChar(cmd1));  //文字列変数に入れる
  PROTECT(call = R_ParseVector(cv, 1, &status, R_NilValue));  //お作法です．
  PROTECT(call2 = lang2(install("eval"), call));  //eval関数にセットする
  for(i=0; i<5; i++){
    INTEGER(a)[0] = INTEGER(a)[0] + i;
    defineVar(install("a"), a, env);
    PROTECT(ans = eval(call2, env));
    UNPROTECT(1);
  }
  UNPROTECT(3);
  return ans;
}




//////////////////////////////////////////////////
// ＢＬＡＳを使うテスト(dgemv)
// y := alpha*A*x + beta*y
long int gl_incx = 1;
long int gl_incy = 1;
double alphaOne = 1.0;
double betaOne = 1.0;
double betaZero = 0.0;
char TRANSN = 'N';
char TRANST = 'T';

void _dgemv(long int m, long int n, double *A, long int lda, double *x, double *ans){
  dgemv_(&TRANSN, &m, &n, &alphaOne, A, &lda, x, &gl_incx, &betaOne, ans, &gl_incy);
}

SEXP Rdgemv01(SEXP ex_A, SEXP ex_x, SEXP ex_y)
{
  SEXP ans;
  int m = INTEGER(GET_DIM(ex_A))[0];
  int n = INTEGER(GET_DIM(ex_A))[1];
  double *A;
  double *x;
  double *y;
  int lda = m;
  long int i;
  
  A = REAL(ex_A);
  x = REAL(ex_x);
  y = REAL(ex_y);
  PROTECT(ans = allocVector(REALSXP, m));
  for( i=0; i<m; i++){
    REAL(ans)[i] = y[i];
  }
  //dgemv_(&TRANSN, &m, &n, &alphaOne, A, &lda, x, &gl_incx, &betaOne, REAL(ans), &gl_incy);
  _dgemv(m, n, A, lda, x, REAL(ans));
  UNPROTECT(1);
  return(ans);
}