Logo Search packages:      
Sourcecode: yorick-soy version File versions  Download package

soy.c

/*
  SOY 1.2 (2005-Apr-11) Sparse Operations with Yorick
  Author: Ralf Flicker (rflicker@mac.com)
  Web site: http://homepage.mac.com/rflicker/soy.htm
  
  This work free software; you may redistribute and modify
  it under the terms of the GNU General Public License.
*/

// gcc -O2 -shared -o soy.so soy.c

#include <math.h>
#include <stdlib.h>
//#include <stdio.h>

int sprco_float(argc,argv)
int argc;
void *argv[];
{
  int i,j,k;
  int *r,*c,*ix,*jx;
  float *t,*x,*xn;
  x = (float *)argv[0];
  r = (int *)argv[1];
  c = (int *)argv[2];
  t = (float *)argv[3];
  ix = (int *)argv[4];
  jx = (int *)argv[5];
  xn = (float *)argv[6];
  k = 0;
  for (i=0; i < *r; i++) {
    for (j=0; j < *c; j++) {      
      if (fabs(x[*c*i+j]) > *t) {
      xn[k] = x[*c*i+j];
      jx[k] = j;
      //printf("%u  %u  %u  %u  %u  \n",i,j,k,ix[i],jx[k]);
      k++;
      }
    }
    ix[i+1] = k;
  }
  return k;
}
//------------------------------------------------
int sprco_double(argc,argv)
int argc;
void *argv[];
{
  int i,j,k;
  int *r,*c,*ix,*jx;
  double *t,*x,*xn;
  x = (double *)argv[0];
  r = (int *)argv[1];
  c = (int *)argv[2];
  t = (double *)argv[3];
  ix = (int *)argv[4];
  jx = (int *)argv[5];
  xn = (double *)argv[6];
  k = 0;
  for (i=0; i < *r; i++) {
    for (j=0; j < *c; j++) {
      if (fabs(x[*c*i+j]) > *t) {
      xn[k] = x[*c*i+j];
      jx[k] = j;
      k++;
      }
    }
    ix[i+1] = k;
  }
  return k;
}
//------------------------------------------------
int spruo_float(argc,argv)
int argc;
void *argv[];
{
  int i,j,k;
  int *r,*ix,*jx;
  float *x,*t,*xn,*xd;
  x = (float *)argv[0];
  r = (int *)argv[1];
  t = (float *)argv[2];
  ix = (int *)argv[3];
  jx = (int *)argv[4];
  xn = (float *)argv[5];
  xd = (float *)argv[6];
  k = 0;
  for (i=0; i < *r; i++) xd[i] = x[*r*i+i];
  for (i=0; i < *r-1; i++) {
    for (j=i+1; j < *r; j++) {
      if (fabs(x[*r*i+j]) > *t) {
      xn[k] = x[*r*i+j];
      jx[k] = j;
      k++;
      }
    }
    ix[i+1] = k;}
  return k;
}
//------------------------------------------------
int spruo_double(argc,argv)
int argc;
void *argv[];
{
  int i,j,k;
  int *r,*ix,*jx;
  double *x,*t,*xn,*xd;
  x = (double *)argv[0];
  r = (int *)argv[1];
  t = (double *)argv[2];
  ix = (int *)argv[3];
  jx = (int *)argv[4];
  xn = (double *)argv[5];
  xd = (double *)argv[6];
  k = 0;
  for (i=0; i < *r; i++) xd[i] = x[*r*i+i];
  for (i=0; i < *r-1; i++) {
    for (j=i+1; j < *r; j++) {
      if (fabs(x[*r*i+j]) > *t) {
      xn[k] = x[*r*i+j];
      jx[k] = j;
      k++;
      }
    }
    ix[i+1] = k;}
  return k;
}
//------------------------------------------------
int rcoxv_float(argc,argv)
int argc;
void *argv[];
{
  int i,j;
  float *v, *u, *xn;
  int *ix,*jx, *r;
  v = (float *)argv[0];
  u = (float *)argv[1];
  r = (int *)argv[2];
  ix = (int *)argv[3];
  jx = (int *)argv[4];
  xn = (float *)argv[5];
  /* struct timeval tod0,tod1; 
     gettimeofday(&tod0,NULL); */ 
  j = 0;
  for (i=0; i < *r; i++) {
    if (ix[i+1]-ix[i] > 0) {
      for (j=ix[i]; j < ix[i+1]; j++) {
      u[i] += xn[j]*v[jx[j]];
      }
    }
  }
  /* gettimeofday(&tod1,NULL); 
     printf("Time X: %u\n",(unsigned long)tod1.tv_usec-(unsigned long)tod0.tv_usec); */
  return j;
}
//------------------------------------------------
int rcoxv_double(argc,argv)
int argc;
void *argv[];
{
  int i,j;
  double *v, *u, *xn;
  int *ix,*jx, *r;
  v = (double *)argv[0];
  u = (double *)argv[1];
  r = (int *)argv[2];
  ix = (int *)argv[3];
  jx = (int *)argv[4];
  xn = (double *)argv[5];
  /* struct timeval tod0,tod1; 
     gettimeofday(&tod0,NULL); */ 
  j = 0;
  for (i=0; i < *r; i++) {
    if (ix[i+1]-ix[i] > 0) {
      for (j=ix[i]; j < ix[i+1]; j++) {
      u[i] += xn[j]*v[jx[j]];
      }
    }
  }
  /* gettimeofday(&tod1,NULL); 
     printf("Time X: %u\n",(unsigned long)tod1.tv_usec-(unsigned long)tod0.tv_usec); */
  return j;
}
//------------------------------------------------
int ruoxv_float(argc,argv)
int argc;
void *argv[];
{
  int i,j;
  int *ix,*jx, *r;
  float *v,*u,*w,*xn,*xd;
  v = (float *)argv[0];
  u = (float *)argv[1];
  w = (float *)argv[2];
  r = (int *)argv[3];
  ix = (int *)argv[4];
  jx = (int *)argv[5];
  xn = (float *)argv[6];
  xd = (float *)argv[7];
  j = 0;
  for (i=0; i < *r; i++) u[i] = xd[i]*v[i];
  for (i=0; i < *r-1; i++) {
    if (ix[i+1] > ix[i]) {
      for (j=ix[i]; j < ix[i+1]; j++) {
      u[i] += xn[j]*v[jx[j]];
      w[jx[j]] += xn[j]*v[i];
      }
    }
  }
  for (i=0; i < *r; i++) u[i] += w[i];
  return j;
}
//------------------------------------------------
int ruoxv_double(argc,argv)
int argc;
void *argv[];
{
  int i,j;
  int *ix,*jx, *r;
  double *v,*u,*w,*xn,*xd;
  v = (double *)argv[0];
  u = (double *)argv[1];
  w = (double *)argv[2];
  r = (int *)argv[3];
  ix = (int *)argv[4];
  jx = (int *)argv[5];
  xn = (double *)argv[6];
  xd = (double *)argv[7];
  j = 0;
  for (i=0; i < *r; i++) u[i] = xd[i]*v[i];
  for (i=0; i < *r-1; i++) {
    if (ix[i+1] > ix[i]) {
      for (j=ix[i]; j < ix[i+1]; j++) {
      u[i] += xn[j]*v[jx[j]];
      w[jx[j]] += xn[j]*v[i];
      }
    }
  }
  for (i=0; i < *r; i++) u[i] += w[i];
  return j;
}
//------------------------------------------------
int rcoadd_float(argc,argv)
int argc;
void *argv[];
{
  int i,j,p,*s;
  float *t,*axn,*bxn,*cxn;
  int *aix,*ajx,*ar,*ac;
  int *bix,*bjx,*br,*bc;
  int *cix,*cjx,*cr,*cc;
  ar = (int *)argv[0];
  ac = (int *)argv[1];
  aix = (int *)argv[2];
  ajx = (int *)argv[3];
  axn = (float *)argv[4];
  br = (int *)argv[5];
  bc = (int *)argv[6];
  bix = (int *)argv[7];
  bjx = (int *)argv[8];
  bxn = (float *)argv[9];
  cr = (int *)argv[10];
  cc = (int *)argv[11];
  cix = (int *)argv[12];
  cjx = (int *)argv[13];
  cxn = (float *)argv[14];
  t = (float *)argv[15];
  s = (int *)argv[16];
  p = 0;


  for (i=0; i < *ac; i++) s[i] = -1;
  for (i=0; i < *ar; i++) {
    cix[i] = p;
    if (aix[i+1] > aix[i]) {
      for (j=aix[i]; j < aix[i+1]; j++) {
      cjx[p++] = ajx[j];
      s[ajx[j]] = i;
      }
    }
    if (bix[i+1] > bix[i]) {
      for (j=bix[i]; j < bix[i+1]; j++) {
      if (s[bjx[j]] != i) cjx[p++] = bjx[j];
      }
    }
  }
  cix[*ar] = p;  
  for (i=0; i < *ar; i++) {
    if (cix[i+1] > cix[i]) {
      for (j=cix[i]; j < cix[i+1]; j++) {
      //printf("%u  %u  %u  %u  \n",j,cix[i],cix[i+1],cjx[j]);
      t[cjx[j]] = 0.0;
      }
      if (aix[i+1] > aix[i]) {
      for (j=aix[i]; j < aix[i+1]; j++) t[ajx[j]] = axn[j];
      }
      if (bix[i+1] > bix[i]) {
      for (j=bix[i]; j < bix[i+1]; j++) t[bjx[j]] = t[bjx[j]]+bxn[j];
      }
      for (j=cix[i]; j < cix[i+1]; j++) cxn[j] = t[cjx[j]];
    }
  }
  return p;
}
//------------------------------------------------
int rcoadd_double(argc,argv)
int argc;
void *argv[];
{
  int i,j,p,*s;
  double *t,*axn,*bxn,*cxn;
  int *aix,*ajx,*ar,*ac;
  int *bix,*bjx,*br,*bc;
  int *cix,*cjx,*cr,*cc;
  ar = (int *)argv[0];
  ac = (int *)argv[1];
  aix = (int *)argv[2];
  ajx = (int *)argv[3];
  axn = (double *)argv[4];
  br = (int *)argv[5];
  bc = (int *)argv[6];
  bix = (int *)argv[7];
  bjx = (int *)argv[8];
  bxn = (double *)argv[9];
  cr = (int *)argv[10];
  cc = (int *)argv[11];
  cix = (int *)argv[12];
  cjx = (int *)argv[13];
  cxn = (double *)argv[14];
  t = (double *)argv[15];
  s = (int *)argv[16];
  p = 0;

  for (i=0; i < *ac; i++) s[i] = -1;
  for (i=0; i < *ar; i++) {
    cix[i] = p;
    if (aix[i+1] > aix[i]) {
      for (j=aix[i]; j < aix[i+1]; j++) {
      cjx[p++] = ajx[j];
      s[ajx[j]] = i;
      }
    }
    if (bix[i+1] > bix[i]) {
      for (j=bix[i]; j < bix[i+1]; j++) {
      if (s[bjx[j]] != i) cjx[p++] = bjx[j];
      }
    }
  }
  cix[*ar] = p;  
  for (i=0; i < *ar; i++) {
    if (cix[i+1] > cix[i]) {
      for (j=cix[i]; j < cix[i+1]; j++) t[cjx[j]] = 0.0;
      if (aix[i+1] > aix[i]) {
      for (j=aix[i]; j < aix[i+1]; j++) t[ajx[j]] = axn[j];
      }
      if (bix[i+1] > bix[i]) {
      for (j=bix[i]; j < bix[i+1]; j++) t[bjx[j]] = t[bjx[j]]+bxn[j];
      }
      for (j=cix[i]; j < cix[i+1]; j++) cxn[j] = t[cjx[j]];
    }
  }
  return p;
}
//------------------------------------------------
int ruoadd_float(argc,argv)
int argc;
void *argv[];
{
  int i,j,p,*s;
  float *tt,*axn,*axd,*bxn,*bxd,*cxn,*cxd;
  int *aix,*ajx,*ar;
  int *bix,*bjx,*br;
  int *cix,*cjx,*cr;
  ar = (int *)argv[0];
  axd = (float *)argv[1];
  aix = (int *)argv[2];
  ajx = (int *)argv[3];
  axn = (float *)argv[4];
  br = (int *)argv[5];
  bxd = (float *)argv[6];
  bix = (int *)argv[7];
  bjx = (int *)argv[8];
  bxn = (float *)argv[9];
  cr = (int *)argv[10];
  cxd = (float *)argv[11];
  cix = (int *)argv[12];
  cjx = (int *)argv[13];
  cxn = (float *)argv[14];
  tt = (float *)argv[15];
  s = (int *)argv[16];

  for (i=0; i < *ar; i++) cxd[i] = axd[i]+bxd[i];
  for (i=0; i < *ar; i++) s[i] = -1;
  p = 0;
  for (i=0; i < *ar-1; i++) {
    cix[i] = p;
    if (aix[i+1] > aix[i]) {
      for (j=aix[i]; j < aix[i+1]; j++) {
      cjx[p++] = ajx[j];
      s[ajx[j]] = i;
      }
    }
    if (bix[i+1] > bix[i]) {
      for (j=bix[i]; j < bix[i+1]; j++) {
      if (s[bjx[j]] != i) cjx[p++] = bjx[j];
      }
    }
  }
  cix[*ar-1] = p;
  for (i=0; i < *ar-1; i++) {
    if (cix[i+1] > cix[i]) {
      for (j=cix[i]; j < cix[i+1]; j++) tt[cjx[j]] = 0.0;
      if (aix[i+1] > aix[i]) {
      for (j=aix[i]; j < aix[i+1]; j++) tt[ajx[j]] = axn[j];
      }
      if (bix[i+1] > bix[i]) {
      for (j=bix[i]; j < bix[i+1]; j++) tt[bjx[j]] = tt[bjx[j]]+bxn[j];
      }
      for (j=cix[i]; j < cix[i+1]; j++) cxn[j] = tt[cjx[j]];
    }
  }
  return p;
}
//------------------------------------------------
int ruoadd_double(argc,argv)
int argc;
void *argv[];
{
  int i,j,p,*s;
  double *tt,*axn,*axd,*bxn,*bxd,*cxn,*cxd;
  int *aix,*ajx,*ar;
  int *bix,*bjx,*br;
  int *cix,*cjx,*cr;
  ar = (int *)argv[0];
  axd = (double *)argv[1];
  aix = (int *)argv[2];
  ajx = (int *)argv[3];
  axn = (double *)argv[4];
  br = (int *)argv[5];
  bxd = (double *)argv[6];
  bix = (int *)argv[7];
  bjx = (int *)argv[8];
  bxn = (double *)argv[9];
  cr = (int *)argv[10];
  cxd = (double *)argv[11];
  cix = (int *)argv[12];
  cjx = (int *)argv[13];
  cxn = (double *)argv[14];
  tt = (double *)argv[15];
  s = (int *)argv[16];

  for (i=0; i < *ar; i++) cxd[i] = axd[i]+bxd[i];
  for (i=0; i < *ar; i++) s[i] = -1;
  p = 0;
  for (i=0; i < *ar-1; i++) {
    cix[i] = p;
    if (aix[i+1] > aix[i]) {
      for (j=aix[i]; j < aix[i+1]; j++) {
      cjx[p++] = ajx[j];
      s[ajx[j]] = i;
      }
    }
    if (bix[i+1] > bix[i]) {
      for (j=bix[i]; j < bix[i+1]; j++) {
      if (s[bjx[j]] != i) cjx[p++] = bjx[j];
      }
    }
  }
  cix[*ar-1] = p;
  for (i=0; i < *ar-1; i++) {
    if (cix[i+1] > cix[i]) {
      for (j=cix[i]; j < cix[i+1]; j++) tt[cjx[j]] = 0.0;
      if (aix[i+1] > aix[i]) {
      for (j=aix[i]; j < aix[i+1]; j++) tt[ajx[j]] = axn[j];
      }
      if (bix[i+1] > bix[i]) {
      for (j=bix[i]; j < bix[i+1]; j++) tt[bjx[j]] = tt[bjx[j]]+bxn[j];
      }
      for (j=cix[i]; j < cix[i+1]; j++) cxn[j] = tt[cjx[j]];
    }
  }
  return p;
}
//------------------------------------------------
int rcoata_float(argc,argv)
int argc;
void *argv[];
{
  int i,j,k,l,p,ni,nj;
  float mb=0.0;
  int *ar,*aix,*ajx,*bix,*bjx;
  float *bt,*axn,*bxn,*bxd;
  ar = (int *)argv[0];
  aix = (int *)argv[1];
  ajx = (int *)argv[2];
  axn = (float *)argv[3];
  bt = (float *)argv[4];
  bix = (int *)argv[5];
  bjx = (int *)argv[6];
  bxn = (float *)argv[7];
  bxd = (float *)argv[8];

  //printf("row ");
  for (i=0; i < *ar; i++) {
    //printf("%u ",i);;
    if (aix[i+1]-aix[i] > 0) {  
      for (k=aix[i]; k < aix[i+1]; k++) {
      bxd[i] += axn[k]*axn[k];
      }
    }
  }
  p = 0;
  //printf("row ");
  for (i=0; i < *ar-1; i++) {
    //printf("i%u ",i);
    for (j=i+1; j < *ar; j++) {
      //printf("j%u ",j);
      ni = aix[i+1]-aix[i];
      nj = aix[j+1]-aix[j];
      if (ni > 0 && nj > 0) {
      mb = 0.0;
      for (k=0; k < ni; k++) {
        //printf("k%u ",k);
        for (l=0; l < nj; l++) {
          //printf("l%u ",l);
          if (ajx[aix[i]+k] == ajx[aix[j]+l]) {
            mb += axn[aix[i]+k]*(axn[aix[j]+l]);
          }
        }
      }
      if (fabs(mb) > *bt) {
        bxn[p] = mb;
        bjx[p] = j;
        p++;
      }
      }
    }
    bix[i+1] = p;
  }
  return p;
}
//------------------------------------------------
int rcoata_double(argc,argv)
int argc;
void *argv[];
{
  int i,j,k,l,p,ni,nj;
  double mb=0.0;
  int *ar,*aix,*ajx,*bix,*bjx;
  double *bt,*axn,*bxn,*bxd;
  ar = (int *)argv[0];
  aix = (int *)argv[1];
  ajx = (int *)argv[2];
  axn = (double *)argv[3];
  bt = (double *)argv[4];
  bix = (int *)argv[5];
  bjx = (int *)argv[6];
  bxn = (double *)argv[7];
  bxd = (double *)argv[8];

  for (i=0; i < *ar; i++) {
    if (aix[i+1]-aix[i] > 0) {  
      for (k=aix[i]; k < aix[i+1]; k++) {
      bxd[i] += axn[k]*axn[k];
      }
    }
  }
  p = 0;
  for (i=0; i < *ar-1; i++) {
    for (j=i+1; j < *ar; j++) {
      ni = aix[i+1]-aix[i];
      nj = aix[j+1]-aix[j];
      if (ni > 0 && nj > 0) {
      mb = 0.0;
      for (k=0; k < ni; k++) {
        for (l=0; l < nj; l++) {
          if (ajx[aix[i]+k] == ajx[aix[j]+l])
            mb += axn[aix[i]+k]*(axn[aix[j]+l]);
        }
      }
      if (fabs(mb) > *bt) {
        bxn[p] = mb;
        bjx[p] = j;
        p++;
      }
      }
    }
    bix[i+1] = p;
  }
  return p;
}
//------------------------------------------------
int rcoatb_float(argc,argv)
int argc;
void *argv[];
{
  int i,j,k,l,p,ni,nj;
  float mb=0.0;
  int *ar,*aix,*ajx,*br,*bix,*bjx,*cix,*cjx;
  float *ct,*axn,*bxn,*cxn;
  ar = (int *)argv[0];
  aix = (int *)argv[1];
  ajx = (int *)argv[2];
  axn = (float *)argv[3];
  br = (int *)argv[4];
  bix = (int *)argv[5];
  bjx = (int *)argv[6];
  bxn = (float *)argv[7];
  ct = (float *)argv[8];
  cix = (int *)argv[9];
  cjx = (int *)argv[10];
  cxn = (float *)argv[11];

  p = 0;
  for (i=0; i < *ar; i++) {
    for (j=0; j < *br; j++) {
      ni = aix[i+1]-aix[i];
      nj = bix[j+1]-bix[j];
      if (ni > 0 && nj > 0) {
      mb = 0.0;
      for (k=0; k < ni; k++) {
        for (l=0; l < nj; l++) {
          if (ajx[aix[i]+k] == bjx[bix[j]+l])
            mb += axn[aix[i]+k]*(bxn[bix[j]+l]);
        }
      }
      if (fabs(mb) > *ct) {
        cxn[p] = mb;
        cjx[p] = j;
        p++;
      }
      }
    }
    cix[i+1] = p;
  }
  return p;
}
//------------------------------------------------
int rcoatb_double(argc,argv)
int argc;
void *argv[];
{
  int i,j,k,l,p,ni,nj;
  double mb=0.0;
  int *ar,*aix,*ajx,*br,*bix,*bjx,*cix,*cjx;
  double *ct,*axn,*bxn,*cxn;
  ar = (int *)argv[0];
  aix = (int *)argv[1];
  ajx = (int *)argv[2];
  axn = (double *)argv[3];
  br = (int *)argv[4];
  bix = (int *)argv[5];
  bjx = (int *)argv[6];
  bxn = (double *)argv[7];
  ct = (double *)argv[8];
  cix = (int *)argv[9];
  cjx = (int *)argv[10];
  cxn = (double *)argv[11];

  p = 0;
  for (i=0; i < *ar; i++) {
    for (j=0; j < *br; j++) {
      ni = aix[i+1]-aix[i];
      nj = bix[j+1]-bix[j];
      if (ni > 0 && nj > 0) {
      mb = 0.0;
      for (k=0; k < ni; k++) {
        for (l=0; l < nj; l++) {
          if (ajx[aix[i]+k] == bjx[bix[j]+l])
            mb += axn[aix[i]+k]*(bxn[bix[j]+l]);
        }
      }
      if (fabs(mb) > *ct) {
        cxn[p] = mb;
        cjx[p] = j;
        p++;
      }
      }
    }
    cix[i+1] = p;
  }
  return p;
}
//------------------------------------------------
int rcoatb2_float(argc,argv)
int argc;
void *argv[];
{
  int i,j,k,l,p,ni,nj;
  float mb=0.0;
  int *ar,*aix,*ajx,*br,*bix,*bjx,*cix,*cjx;
  float *ct,*axn,*bxn,*cxn,*cxd;
  ar = (int *)argv[0];
  aix = (int *)argv[1];
  ajx = (int *)argv[2];
  axn = (float *)argv[3];
  br = (int *)argv[4];
  bix = (int *)argv[5];
  bjx = (int *)argv[6];
  bxn = (float *)argv[7];
  ct = (float *)argv[8];
  cix = (int *)argv[9];
  cjx = (int *)argv[10];
  cxn = (float *)argv[11];
  cxd = (float *)argv[12];

  p = 0;
  for (i=0; i < *ar; i++) {
    for (j=i; j < *br; j++) {
      ni = aix[i+1]-aix[i];
      nj = bix[j+1]-bix[j];
      if (ni > 0 && nj > 0) {
      mb = 0.0;
      for (k=0; k < ni; k++) {
        for (l=0; l < nj; l++) {
          if (ajx[aix[i]+k] == bjx[bix[j]+l])
            mb += axn[aix[i]+k]*(bxn[bix[j]+l]);
        }
      }
      if (fabs(mb) > *ct && i != j) {
        cxn[p] = mb;
        cjx[p] = j;
        p++;
      } else if (i == j) {
        cxd[i] = mb;
      }
      }
    }
    if (i != j) {cix[i+1] = p;}
  }
  return p;
}
//------------------------------------------------
int rcotr_fix(argc,argv)
int argc;
void *argv[];
{
  int i,j,nz;
  int *ax,*acx,*hjx,*rind,*ar,*ac,*an,*aix;
  ax = (int *)argv[0];
  acx = (int *)argv[1];
  hjx = (int *)argv[2];
  rind = (int *)argv[3];
  ar = (int *)argv[4];
  ac = (int *)argv[5];
  an = (int *)argv[6];
  aix = (int *)argv[7];
  nz = 0;
  for (i=0; i < *an; i++) {ax[hjx[i]] += 1;}
  for (i=1; i < *ac+1; i++) {acx[i] = ax[i-1]+acx[i-1];}
  for (i=0; i < *ar; i++) {
    nz = aix[i+1]-aix[i];
    if (nz > 0) {
      for (j=0; j < nz; j++) {rind[aix[i]+j] = i;
      }
    }
  }
  return i;
}
//------------------------------------------------
int laplace1_float(argc,argv)
int argc;
void *argv[];
{
  int i, xc, yc;
  int *x, *y, *ind, *inds, *validx, *nx, *dim;
  x = (int *)argv[0];
  y = (int *)argv[1];
  ind = (int *)argv[2];
  inds = (int *)argv[3];
  validx = (int *)argv[4];
  nx = (int *)argv[5];
  dim = (int *)argv[6];
  for (i=0; i < *nx; i++) {
    //    printf("nx = %u",i);
    xc = x[ind[i]-1]-1;
    yc = y[ind[i]-1]-1;
    if (xc < *dim-1) validx[i*4] = inds[*dim*yc+xc+1] > 0;
    if (yc < *dim-1) validx[i*4+1] = inds[*dim*(yc+1)+xc] > 0;
    if (xc > 0) validx[i*4+2] = inds[*dim*yc+xc-1] > 0;
    if (yc > 0) validx[i*4+3] = inds[*dim*(yc-1)+xc] > 0;
  }
  return i;
}
//------------------------------------------------
int laplace2_float(argc,argv)
int argc;
void *argv[];
{
  int i,j,cn,cp,xc,yc,*nx,*dim;
  float *xn,*v5,*v4,*v3;
  int *ix,*jx,*validx,*inds,*aind;
  ix = (int *)argv[0];
  jx = (int *)argv[1];
  xn = (float *)argv[2];
  validx = (int *)argv[3];
  nx = (int *)argv[4];
  dim = (int *)argv[5];
  inds = (int *)argv[6];
  aind = (int *)argv[7];
  v5 = (float *)argv[8];
  v4 = (float *)argv[9];
  v3 = (float *)argv[10];
  cp = 0;
  for (i=0; i < *nx; i++) {
    cn = 0;
    jx[cp] = i;
    cn++;
    xc = aind[i]-1;
    yc = aind[*nx+i]-1;
    if (validx[i*4] == 1) {
      jx[cp+cn] = i+1;
      cn++;
    }
    if (validx[i*4+1] == 1) {
      jx[cp+cn] = inds[*dim*(yc+1)+xc]-1;
      cn++;
    }
    if (validx[i*4+2] == 1) {
      jx[cp+cn] = i-1;
      cn++;
    }
    if (validx[i*4+3] == 1) {
      jx[cp+cn] = inds[*dim*(yc-1)+xc]-1;
      cn++;
    }
    if (cn == 5) {
      for (j=0; j < cn; j++) {
      xn[cp+j] = v5[j];
      }
    } else if (cn == 4) {
      for (j=0; j < cn; j++) {
      xn[cp+j] = v4[j];
      }
    } else if (cn == 3){
      for (j=0; j < cn; j++) {
      xn[cp+j] = v3[j];
      }
    }
    cp += cn;
    ix[i+1] = cp;
  }
  return cp;
}

//------------------------------------------------

Generated by  Doxygen 1.6.0   Back to index