#include <R.h>
#include <Rinternals.h>
#include <Rmath.h>

SEXP upTri(SEXP RinMatrix){

  int i,j,I,J;
  SEXP Rdim;
  SEXP Rval;

  Rdim = getAttrib(RinMatrix, R_DimSymbol);
  I = INTEGER(Rdim)[0];
  J = INTEGER(Rdim)[1];

  if (I != J)
    error("Input must be a square matrix");
    RinMatrix = coerceVector(RinMatrix, REALSXP);

    PROTECT(Rval = allocMatrix(REALSXP, I, J));

    for (i = 0; i < I; i++)
      for (j = 0; j < I; j++)
        if (i <= j)
          REAL(Rval)[i + I * j] =
          REAL(RinMatrix)[i + I * j];
        else
        REAL(Rval)[i + I * j] = 0;
    
    UNPROTECT(1);
    return Rval;
}