swh:1:snp:ffdd0a7d2c8ea15ad41d45b3b178f668bd942287
Tip revision: 57f09c35d67acf681aaca0a9b3869bacdcec0585 authored by Derek Young on 30 March 2009, 00:00:00 UTC
version 0.4.1
version 0.4.1
Tip revision: 57f09c3
new_svalues.c
#include <math.h> /* in order to make the sqrt function work */
/* First, create the function prototype (global declaration that there
exists a function called new_svalues returning void and taking these
arguments) */
void new_svalues (double *z, double *y, double *x, double *beta,
int *k, int *n, int *p, double *out, double *sz, double *runsum);
/* Next comes the function itself: */
void new_svalues (double *z, double *y, double *x, double *beta,
int *k, int *n, int *p, double *out, double *sz, double *runsum) {
int i, j, l;
double sum;
double xbeta;
double diff;
double diff2;
double zdiff;
double zdiff2;
double zdiff3;
/* Create the sz (column sum of z) vector */
for(j=0; j < *k ; j=j+1) {
sum=0.0;
for(i=0; i < *n; i=i+1) {
sum += z[j * (*n) + i];
}
sz[j]=sum;
}
for(j=0; j < *k ; j=j+1) {
zdiff=0.0;
for(i=0; i < *n; i=i+1) {
xbeta=0.0;
/* Calculate i,j component of x %*% beta */
for(l=0; l < *p; l++){
xbeta += x[l * (*n) + i] * beta[j * (*p) + l];
}
/* Subtract it from the i component of y */
diff = y[i] - xbeta;
/* square the difference */
diff2 = pow(diff,2);
/* multiply by the i,j component of z */
zdiff += z[j * (*n) + i] * diff2;
}
/* keep track of running sum */
runsum[j] = zdiff;
/* divide dot product by sz[j] */
zdiff2 = runsum[j] / sz[j];
/* take the square root */
zdiff3 = sqrt(zdiff2);
/* put result in out[j] */
out[j] = zdiff3;
}
}